diff options
Diffstat (limited to 'src/modules/audio_processing/aecm/main/matlab/compsup.m')
| -rw-r--r-- | src/modules/audio_processing/aecm/main/matlab/compsup.m | 447 |
1 files changed, 447 insertions, 0 deletions
diff --git a/src/modules/audio_processing/aecm/main/matlab/compsup.m b/src/modules/audio_processing/aecm/main/matlab/compsup.m new file mode 100644 index 0000000..9575ec4 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/compsup.m @@ -0,0 +1,447 @@ +function [emicrophone,aaa]=compsup(microphone,TheFarEnd,avtime,samplingfreq); +% microphone = microphone signal +% aaa = nonlinearity input variable +% TheFarEnd = far end signal +% avtime = interval to compute suppression from (seconds) +% samplingfreq = sampling frequency + +%if(nargin==6) +% fprintf(1,'suppress has received a delay sequence\n'); +%end + + +Ap500=[ 1.00, -4.95, 9.801, -9.70299, 4.80298005, -0.9509900499]; +Bp500=[ 0.662743088639636, -2.5841655608125, 3.77668102146288, -2.45182477425154, 0.596566274575251, 0.0]; + + +Ap200=[ 1.00, -4.875, 9.50625, -9.26859375, 4.518439453125, -0.881095693359375]; +Bp200=[ 0.862545460994275, -3.2832804496114, 4.67892032308828, -2.95798023879133, 0.699796870041299, 0.0]; + +maxDelay=0.4; %[s] +histLen=1; %[s] + + +% CONSTANTS THAT YOU CAN EXPERIMENT WITH +A_GAIN=10.0; % for the suppress case +oversampling = 2; % must be power of 2; minimum is 2; 4 works +% fine for support=64, but for support=128, +% 8 gives better results. +support=64; %512 % fft support (frequency resolution; at low +% settings you can hear more distortion +% (e.g. pitch that is left-over from far-end)) +% 128 works well, 64 is ok) + +lowlevel = mean(abs(microphone))*0.0001; + +G_ol = 0; % Use overlapping sets of estimates + +% ECHO SUPPRESSION SPECIFIC PARAMETERS +suppress_overdrive=1.0; % overdrive factor for suppression 1.4 is good +gamma_echo=1.0; % same as suppress_overdrive but at different place +de_echo_bound=0.0; +mLim=10; % rank of matrix G +%limBW = 1; % use bandwidth-limited response for G +if mLim > (support/2+1) + error('mLim in suppress.m too large\n'); +end + + +dynrange=1.0000e-004; + +% other, constants +hsupport = support/2; +hsupport1 = hsupport+1; +factor = 2 / oversampling; +updatel = support/oversampling; +win=sqrt(designwindow(0,support)); +estLen = round(avtime * samplingfreq/updatel) + +runningfmean =0.0; + +mLim = floor(hsupport1/2); +V = sqrt(2/hsupport1)*cos(pi/hsupport1*(repmat((0:hsupport1-1) + 0.5, mLim, 1).* ... + repmat((0:mLim-1)' + 0.5, 1, hsupport1))); + +fprintf(1,'updatel is %5.3f s\n', updatel/samplingfreq); + + + +bandfirst=8; bandlast=25; +dosmooth=0; % to get rid of wavy bin counts (can be worse or better) + +% compute some constants +blockLen = support/oversampling; +maxDelayb = floor(samplingfreq*maxDelay/updatel); % in blocks +histLenb = floor(samplingfreq*histLen/updatel); % in blocks + +x0=TheFarEnd; +y0=microphone; + + +%input +tlength=min([length(microphone),length(TheFarEnd)]); +updateno=floor(tlength/updatel); +tlength=updatel*updateno; +updateno = updateno - oversampling + 1; + +TheFarEnd =TheFarEnd(1:tlength); +microphone =microphone(1:tlength); + +TheFarEnd =[zeros(hsupport,1);TheFarEnd(1:tlength)]; +microphone =[zeros(hsupport,1);microphone(1:tlength)]; + + +% signal length +n = min([floor(length(x0)/support)*support,floor(length(y0)/support)*support]); +nb = n/blockLen - oversampling + 1; % in blocks + +% initialize space +win = sqrt([0 ; hanning(support-1)]); +sxAll2 = zeros(hsupport1,nb); +syAll2 = zeros(hsupport1,nb); + +z500=zeros(5,maxDelayb+1); +z200=zeros(5,hsupport1); + +bxspectrum=uint32(zeros(nb,1)); +bxhist=uint32(zeros(maxDelayb+1,1)); +byspectrum=uint32(zeros(nb,1)); +bcount=zeros(1+maxDelayb,nb); +fcount=zeros(1+maxDelayb,nb); +fout=zeros(1+maxDelayb,nb); +delay=zeros(nb,1); +tdelay=zeros(nb,1); +nlgains=zeros(nb,1); + +% create space (mainly for debugging) +emicrophone=zeros(tlength,1); +femicrophone=complex(zeros(hsupport1,updateno)); +thefilter=zeros(hsupport1,updateno); +thelimiter=ones(hsupport1,updateno); +fTheFarEnd=complex(zeros(hsupport1,updateno)); +afTheFarEnd=zeros(hsupport1,updateno); +fmicrophone=complex(zeros(hsupport1,updateno)); +afmicrophone=zeros(hsupport1,updateno); + +G = zeros(hsupport1, hsupport1); +zerovec = zeros(hsupport1,1); +zeromat = zeros(hsupport1); + +% Reset sums +mmxs_a = zerovec; +mmys_a = zerovec; +s2xs_a = zerovec; +s2ys_a = zerovec; +Rxxs_a = zeromat; +Ryxs_a = zeromat; +count_a = 1; + +mmxs_b = zerovec; +mmys_b = zerovec; +s2xs_b = zerovec; +s2ys_b = zerovec; +Rxxs_b = zeromat; +Ryxs_b = zeromat; +count_b = 1; + +nog=0; + +aaa=zeros(size(TheFarEnd)); + +% loop over signal blocks +fprintf(1,'.. Suppression; averaging G over %5.1f seconds; file length %5.1f seconds ..\n',avtime, length(microphone)/samplingfreq); +fprintf(1,'.. SUPPRESSING ONLY AFTER %5.1f SECONDS! ..\n',avtime); +fprintf(1,'.. 20 seconds is good ..\n'); +hh = waitbar_j(0,'Please wait...'); + + +for i=1:updateno + + sb = (i-1)*updatel + 1; + se=sb+support-1; + + % analysis FFTs + temp=fft(win .* TheFarEnd(sb:se)); + fTheFarEnd(:,i)=temp(1:hsupport1); + xf=fTheFarEnd(:,i); + afTheFarEnd(:,i)= abs(fTheFarEnd(:,i)); + + temp=win .* microphone(sb:se); + + temp=fft(win .* microphone(sb:se)); + fmicrophone(:,i)=temp(1:hsupport1); + yf=fmicrophone(:,i); + afmicrophone(:,i)= abs(fmicrophone(:,i)); + + + ener_orig = afmicrophone(:,i)'*afmicrophone(:,i); + if( ener_orig == 0) + afmicrophone(:,i)=lowlevel*ones(size(afmicrophone(:,i))); + end + + + % use log domain (showed improved performance) +xxf= sqrt(real(xf.*conj(xf))+1e-20); +yyf= sqrt(real(yf.*conj(yf))+1e-20); + sxAll2(:,i) = 20*log10(xxf); + syAll2(:,i) = 20*log10(yyf); + + mD=min(i-1,maxDelayb); + xthreshold = sum(sxAll2(:,i-mD:i),2)/(maxDelayb+1); + + [yout, z200] = filter(Bp200,Ap200,syAll2(:,i),z200,2); + yout=yout/(maxDelayb+1); + ythreshold = mean(syAll2(:,i-mD:i),2); + + + bxspectrum(i)=getBspectrum(sxAll2(:,i),xthreshold,bandfirst,bandlast); + byspectrum(i)=getBspectrum(syAll2(:,i),yout,bandfirst,bandlast); + + bxhist(end-mD:end)=bxspectrum(i-mD:i); + + bcount(:,i)=hisser2( ... + byspectrum(i),flipud(bxhist),bandfirst,bandlast); + + + [fout(:,i), z500] = filter(Bp500,Ap500,bcount(:,i),z500,2); + fcount(:,i)=sum(bcount(:,max(1,i-histLenb+1):i),2); % using the history range + fout(:,i)=round(fout(:,i)); + [value,delay(i)]=min(fout(:,i),[],1); + tdelay(i)=(delay(i)-1)*support/(samplingfreq*oversampling); + + % compensate + + idel = max(i - delay(i) + 1,1); + + + % echo suppression + + noisyspec = afmicrophone(:,i); + + % Estimate G using covariance matrices + + % Cumulative estimates + xx = afTheFarEnd(:,idel); + yy = afmicrophone(:,i); + + % Means + mmxs_a = mmxs_a + xx; + mmys_a = mmys_a + yy; + if (G_ol) + mmxs_b = mmxs_b + xx; + mmys_b = mmys_b + yy; + mmy = mean([mmys_a/count_a mmys_b/count_b],2); + mmx = mean([mmxs_a/count_a mmxs_b/count_b],2); + else + mmx = mmxs_a/count_a; + mmy = mmys_a/count_a; + end + count_a = count_a + 1; + count_b = count_b + 1; + + % Mean removal + xxm = xx - mmx; + yym = yy - mmy; + + % Variances + s2xs_a = s2xs_a + xxm .* xxm; + s2ys_a = s2ys_a + yym .* yym; + s2xs_b = s2xs_b + xxm .* xxm; + s2ys_b = s2ys_b + yym .* yym; + + % Correlation matrices + Rxxs_a = Rxxs_a + xxm * xxm'; + Ryxs_a = Ryxs_a + yym * xxm'; + Rxxs_b = Rxxs_b + xxm * xxm'; + Ryxs_b = Ryxs_b + yym * xxm'; + + + % Gain matrix A + + if mod(i, estLen) == 0 + + + % Cumulative based estimates + Rxxf = Rxxs_a / (estLen - 1); + Ryxf = Ryxs_a / (estLen - 1); + + % Variance normalization + s2x2 = s2xs_a / (estLen - 1); + s2x2 = sqrt(s2x2); + % Sx = diag(max(s2x2,dynrange*max(s2x2))); + Sx = diag(s2x2); + if (sum(s2x2) > 0) + iSx = inv(Sx); + else + iSx= Sx + 0.01; + end + + s2y2 = s2ys_a / (estLen - 1); + s2y2 = sqrt(s2y2); + % Sy = diag(max(s2y2,dynrange*max(s2y2))); + Sy = diag(s2y2); + iSy = inv(Sy); + rx = iSx * Rxxf * iSx; + ryx = iSy * Ryxf * iSx; + + + + dbd= 7; % Us less than the full matrix + + % k x m + % Bandlimited structure on G + LSEon = 0; % Default is using MMSE + if (LSEon) + ryx = ryx*rx; + rx = rx*rx; + end + p = dbd-1; + gaj = min(min(hsupport1,2*p+1),min([p+(1:hsupport1); hsupport1+p+1-(1:hsupport1)])); + cgaj = [0 cumsum(gaj)]; + + G3 = zeros(hsupport1); + for kk=1:hsupport1 + ki = max(0,kk-p-1); + if (sum(sum(rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk))))>0) + G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk))/rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk)); + else + G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk)); + end + end + % End Bandlimited structure + + G = G3; + G(abs(G)<0.01)=0; + G = suppress_overdrive * Sy * G * iSx; + + if 1 + figure(32); mi=2; + surf(max(min(G,mi),-mi)); view(2) + title('Unscaled Masked Limited-bandwidth G'); + end + pause(0.05); + + % Reset sums + mmxs_a = zerovec; + mmys_a = zerovec; + s2xs_a = zerovec; + s2ys_a = zerovec; + Rxxs_a = zeromat; + Ryxs_a = zeromat; + count_a = 1; + + end + + if (G_ol) + % Gain matrix B + + if ((mod((i-estLen/2), estLen) == 0) & i>estLen) + + + % Cumulative based estimates + Rxxf = Rxxs_b / (estLen - 1); + Ryxf = Ryxs_b / (estLen - 1); + + % Variance normalization + s2x2 = s2xs_b / (estLen - 1); + s2x2 = sqrt(s2x2); + Sx = diag(max(s2x2,dynrange*max(s2x2))); + iSx = inv(Sx); + s2y2 = s2ys_b / (estLen - 1); + s2y2 = sqrt(s2y2); + Sy = diag(max(s2y2,dynrange*max(s2y2))); + iSy = inv(Sy); + rx = iSx * Rxxf * iSx; + ryx = iSy * Ryxf * iSx; + + + % Bandlimited structure on G + LSEon = 0; % Default is using MMSE + if (LSEon) + ryx = ryx*rx; + rx = rx*rx; + end + p = dbd-1; + gaj = min(min(hsupport1,2*p+1),min([p+(1:hsupport1); hsupport1+p+1-(1:hsupport1)])); + cgaj = [0 cumsum(gaj)]; + + G3 = zeros(hsupport1); + for kk=1:hsupport1 + ki = max(0,kk-p-1); + G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk))/rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk)); + end + % End Bandlimited structure + + G = G3; + G(abs(G)<0.01)=0; + G = suppress_overdrive * Sy * G * iSx; + + if 1 + figure(32); mi=2; + surf(max(min(G,mi),-mi)); view(2) + title('Unscaled Masked Limited-bandwidth G'); + end + pause(0.05); + + + % Reset sums + mmxs_b = zerovec; + mmys_b = zerovec; + s2xs_b = zerovec; + s2ys_b = zerovec; + Rxxs_b = zeromat; + Ryxs_b = zeromat; + count_b = 1; + + end + + end + + FECestimate2 = G*afTheFarEnd(:,idel); + + % compute Wiener filter and suppressor function + thefilter(:,i) = (noisyspec - gamma_echo*FECestimate2) ./ noisyspec; + ix0 = find(thefilter(:,i)<de_echo_bound); % bounding trick 1 + thefilter(ix0,i) = de_echo_bound; % bounding trick 2 + ix0 = find(thefilter(:,i)>1); % bounding in reasonable range + thefilter(ix0,i) = 1; + + % NONLINEARITY + nl_alpha=0.8; % memory; seems not very critical + nlSeverity=0.3; % nonlinearity severity: 0 does nothing; 1 suppresses all + thefmean=mean(thefilter(8:16,i)); + if (thefmean<1) + disp(''); + end + runningfmean = nl_alpha*runningfmean + (1-nl_alpha)*thefmean; + aaa(sb+20+1:sb+20+updatel)=10000*runningfmean* ones(updatel,1); % debug + slope0=1.0/(1.0-nlSeverity); % + thegain = max(0.0,min(1.0,slope0*(runningfmean-nlSeverity))); + % END NONLINEARITY + thefilter(:,i) = thegain*thefilter(:,i); + + + % Wiener filtering + femicrophone(:,i) = fmicrophone(:,i) .* thefilter(:,i); + thelimiter(:,i) = (noisyspec - A_GAIN*FECestimate2) ./ noisyspec; + index = find(thelimiter(:,i)>1.0); + thelimiter(index,i) = 1.0; + index = find(thelimiter(:,i)<0.0); + thelimiter(index,i) = 0.0; + + if (rem(i,floor(updateno/20))==0) + fprintf(1,'.'); + end + if mod(i,50)==0 + waitbar_j(i/updateno,hh); + end + + + % reconstruction; first make spectrum odd + temp=[femicrophone(:,i);flipud(conj(femicrophone(2:hsupport,i)))]; + emicrophone(sb:se) = emicrophone(sb:se) + factor * win .* real(ifft(temp)); + +end +fprintf(1,'\n'); + +close(hh);
\ No newline at end of file |