diff options
Diffstat (limited to 'src/modules/audio_processing/aecm/main/matlab')
15 files changed, 1780 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 0000000000..9575ec40fc --- /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 diff --git a/src/modules/audio_processing/aecm/main/matlab/getBspectrum.m b/src/modules/audio_processing/aecm/main/matlab/getBspectrum.m new file mode 100644 index 0000000000..a4a533d600 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/getBspectrum.m @@ -0,0 +1,22 @@ +function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast) +% function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast) +% compute binary spectrum using threshold spectrum as pivot +% bspectrum = binary spectrum (binary) +% ps=current power spectrum (float) +% threshold=threshold spectrum (float) +% bandfirst = first band considered +% bandlast = last band considered + +% initialization stuff + if( length(ps)<bandlast | bandlast>32 | length(ps)~=length(threshold)) + error('BinDelayEst:spectrum:invalid','Dimensionality error'); +end + +% get current binary spectrum +diff = ps - threshold; +bspectrum=uint32(0); +for(i=bandfirst:bandlast) + if( diff(i)>0 ) + bspectrum = bitset(bspectrum,i); + end +end diff --git a/src/modules/audio_processing/aecm/main/matlab/hisser2.m b/src/modules/audio_processing/aecm/main/matlab/hisser2.m new file mode 100644 index 0000000000..5a414f9da8 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/hisser2.m @@ -0,0 +1,21 @@ +function bcount=hisser2(bs,bsr,bandfirst,bandlast) +% function bcount=hisser(bspectrum,bandfirst,bandlast) +% histogram for the binary spectra +% bcount= array of bit counts +% bs=binary spectrum (one int32 number each) +% bsr=reference binary spectra (one int32 number each) +% blockSize = histogram over blocksize blocks +% bandfirst = first band considered +% bandlast = last band considered + +% weight all delays equally +maxDelay = length(bsr); + +% compute counts (two methods; the first works better and is operational) +bcount=zeros(maxDelay,1); +for(i=1:maxDelay) + % the delay should have low count for low-near&high-far and high-near&low-far + bcount(i)= sum(bitget(bitxor(bs,bsr(i)),bandfirst:bandlast)); + % the delay should have low count for low-near&high-far (works less well) +% bcount(i)= sum(bitget(bitand(bsr(i),bitxor(bs,bsr(i))),bandfirst:bandlast)); +end diff --git a/src/modules/audio_processing/aecm/main/matlab/main2.m b/src/modules/audio_processing/aecm/main/matlab/main2.m new file mode 100644 index 0000000000..7e24c69ccf --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/main2.m @@ -0,0 +1,19 @@ + +fid=fopen('aecfar.pcm'); far=fread(fid,'short'); fclose(fid); +fid=fopen('aecnear.pcm'); mic=fread(fid,'short'); fclose(fid); + +%fid=fopen('QA1far.pcm'); far=fread(fid,'short'); fclose(fid); +%fid=fopen('QA1near.pcm'); mic=fread(fid,'short'); fclose(fid); + +start=0 * 8000+1; +stop= 30 * 8000; +microphone=mic(start:stop); +TheFarEnd=far(start:stop); +avtime=1; + +% 16000 to make it compatible with the C-version +[emicrophone,tdel]=compsup(microphone,TheFarEnd,avtime,16000); + +spclab(8000,TheFarEnd,microphone,emicrophone); + + diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/AECMobile.m b/src/modules/audio_processing/aecm/main/matlab/matlab/AECMobile.m new file mode 100644 index 0000000000..2d3e6867df --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/AECMobile.m @@ -0,0 +1,269 @@ +function [femicrophone, aecmStructNew, enerNear, enerFar] = AECMobile(fmicrophone, afTheFarEnd, setupStruct, aecmStruct) +global NEARENDFFT; +global F; + +aecmStructNew = aecmStruct; + +% Magnitude spectrum of near end signal +afmicrophone = abs(fmicrophone); +%afmicrophone = NEARENDFFT(setupStruct.currentBlock,:)'/2^F(setupStruct.currentBlock,end); +% Near end energy level +ener_orig = afmicrophone'*afmicrophone; +if( ener_orig == 0) + lowlevel = 0.01; + afmicrophone = lowlevel*ones(size(afmicrophone)); +end +%adiff = max(abs(afmicrophone - afTheFarEnd)); +%if (adiff > 0) +% disp([setupStruct.currentBlock adiff]) +%end + +% Store the near end energy +%aecmStructNew.enerNear(setupStruct.currentBlock) = log(afmicrophone'*afmicrophone); +aecmStructNew.enerNear(setupStruct.currentBlock) = log(sum(afmicrophone)); +% Store the far end energy +%aecmStructNew.enerFar(setupStruct.currentBlock) = log(afTheFarEnd'*afTheFarEnd); +aecmStructNew.enerFar(setupStruct.currentBlock) = log(sum(afTheFarEnd)); + +% Update subbands (We currently use all frequency bins, hence .useSubBand is turned off) +if aecmStructNew.useSubBand + internalIndex = 1; + for kk=1:setupStruct.subBandLength+1 + ySubBand(kk) = mean(afmicrophone(internalIndex:internalIndex+setupStruct.numInBand(kk)-1).^aecmStructNew.bandFactor); + xSubBand(kk) = mean(afTheFarEnd(internalIndex:internalIndex+setupStruct.numInBand(kk)-1).^aecmStructNew.bandFactor); + internalIndex = internalIndex + setupStruct.numInBand(kk); + end +else + ySubBand = afmicrophone.^aecmStructNew.bandFactor; + xSubBand = afTheFarEnd.^aecmStructNew.bandFactor; +end + +% Estimated echo energy +if (aecmStructNew.bandFactor == 1) + %aecmStructNew.enerEcho(setupStruct.currentBlock) = log((aecmStructNew.H.*xSubBand)'*(aecmStructNew.H.*xSubBand)); + %aecmStructNew.enerEchoStored(setupStruct.currentBlock) = log((aecmStructNew.HStored.*xSubBand)'*(aecmStructNew.HStored.*xSubBand)); + aecmStructNew.enerEcho(setupStruct.currentBlock) = log(sum(aecmStructNew.H.*xSubBand)); + aecmStructNew.enerEchoStored(setupStruct.currentBlock) = log(sum(aecmStructNew.HStored.*xSubBand)); +elseif (aecmStructNew.bandFactor == 2) + aecmStructNew.enerEcho(setupStruct.currentBlock) = log(aecmStructNew.H'*xSubBand); + aecmStructNew.enerEchoStored(setupStruct.currentBlock) = log(aecmStructNew.HStored'*xSubBand); +end + +% Last 100 blocks of data, used for plotting +n100 = max(1,setupStruct.currentBlock-99):setupStruct.currentBlock; +enerError = aecmStructNew.enerNear(n100)-aecmStructNew.enerEcho(n100); +enerErrorStored = aecmStructNew.enerNear(n100)-aecmStructNew.enerEchoStored(n100); + +% Store the far end sub band. This is needed if we use LSE instead of NLMS +aecmStructNew.X = [xSubBand aecmStructNew.X(:,1:end-1)]; + +% Update energy levels, which control the VAD +if ((aecmStructNew.enerFar(setupStruct.currentBlock) < aecmStructNew.energyMin) & (aecmStructNew.enerFar(setupStruct.currentBlock) >= aecmStruct.FAR_ENERGY_MIN)) + aecmStructNew.energyMin = aecmStructNew.enerFar(setupStruct.currentBlock); + %aecmStructNew.energyMin = max(aecmStructNew.energyMin,12); + aecmStructNew.energyMin = max(aecmStructNew.energyMin,aecmStruct.FAR_ENERGY_MIN); + aecmStructNew.energyLevel = (aecmStructNew.energyMax-aecmStructNew.energyMin)*aecmStructNew.energyThres+aecmStructNew.energyMin; + aecmStructNew.energyLevelMSE = (aecmStructNew.energyMax-aecmStructNew.energyMin)*aecmStructNew.energyThresMSE+aecmStructNew.energyMin; +end +if (aecmStructNew.enerFar(setupStruct.currentBlock) > aecmStructNew.energyMax) + aecmStructNew.energyMax = aecmStructNew.enerFar(setupStruct.currentBlock); + aecmStructNew.energyLevel = (aecmStructNew.energyMax-aecmStructNew.energyMin)*aecmStructNew.energyThres+aecmStructNew.energyMin; + aecmStructNew.energyLevelMSE = (aecmStructNew.energyMax-aecmStructNew.energyMin)*aecmStructNew.energyThresMSE+aecmStructNew.energyMin; +end + +% Calculate current energy error in near end (estimated echo vs. near end) +dE = aecmStructNew.enerNear(setupStruct.currentBlock)-aecmStructNew.enerEcho(setupStruct.currentBlock); + +%%%%%%%% +% Calculate step size used in LMS algorithm, based on current far end energy and near end energy error (dE) +%%%%%%%% +if setupStruct.stepSize_flag + [mu, aecmStructNew] = calcStepSize(aecmStructNew.enerFar(setupStruct.currentBlock), dE, aecmStructNew, setupStruct.currentBlock, 1); +else + mu = 0.25; +end +aecmStructNew.muLog(setupStruct.currentBlock) = mu; % Store the step size + +% Estimate Echo Spectral Shape +[U, aecmStructNew.H] = fallerEstimator(ySubBand,aecmStructNew.X,aecmStructNew.H,mu); + +%%%%% +% Determine if we should store or restore the channel +%%%%% +if ((setupStruct.currentBlock <= aecmStructNew.convLength) | (~setupStruct.channelUpdate_flag)) + aecmStructNew.HStored = aecmStructNew.H; % Store what you have after startup +elseif ((setupStruct.currentBlock > aecmStructNew.convLength) & (setupStruct.channelUpdate_flag)) + if ((aecmStructNew.enerFar(setupStruct.currentBlock) < aecmStructNew.energyLevelMSE) & (aecmStructNew.enerFar(setupStruct.currentBlock-1) >= aecmStructNew.energyLevelMSE)) + xxx = aecmStructNew.countMseH; + if (xxx > 20) + mseStored = mean(abs(aecmStructNew.enerEchoStored(setupStruct.currentBlock-xxx:setupStruct.currentBlock-1)-aecmStructNew.enerNear(setupStruct.currentBlock-xxx:setupStruct.currentBlock-1))); + mseLatest = mean(abs(aecmStructNew.enerEcho(setupStruct.currentBlock-xxx:setupStruct.currentBlock-1)-aecmStructNew.enerNear(setupStruct.currentBlock-xxx:setupStruct.currentBlock-1))); + %fprintf('Stored: %4f Latest: %4f\n', mseStored, mseLatest) % Uncomment if you want to display the MSE values + if ((mseStored < 0.8*mseLatest) & (aecmStructNew.mseHStoredOld < 0.8*aecmStructNew.mseHLatestOld)) + aecmStructNew.H = aecmStructNew.HStored; + fprintf('Restored H at block %d\n',setupStruct.currentBlock) + elseif (((0.8*mseStored > mseLatest) & (mseLatest < aecmStructNew.mseHThreshold) & (aecmStructNew.mseHLatestOld < aecmStructNew.mseHThreshold)) | (mseStored == Inf)) + aecmStructNew.HStored = aecmStructNew.H; + fprintf('Stored new H at block %d\n',setupStruct.currentBlock) + end + aecmStructNew.mseHStoredOld = mseStored; + aecmStructNew.mseHLatestOld = mseLatest; + end + elseif ((aecmStructNew.enerFar(setupStruct.currentBlock) >= aecmStructNew.energyLevelMSE) & (aecmStructNew.enerFar(setupStruct.currentBlock-1) < aecmStructNew.energyLevelMSE)) + aecmStructNew.countMseH = 1; + elseif (aecmStructNew.enerFar(setupStruct.currentBlock) >= aecmStructNew.energyLevelMSE) + aecmStructNew.countMseH = aecmStructNew.countMseH + 1; + end +end + +%%%%% +% Check delay (calculate the delay offset (if we can)) +% The algorithm is not tuned and should be used with care. It runs separately from Bastiaan's algorithm. +%%%%% +yyy = 31; % Correlation buffer length (currently unfortunately hard coded) +dxxx = 25; % Maximum offset (currently unfortunately hard coded) +if (setupStruct.currentBlock > aecmStructNew.convLength) + if (aecmStructNew.enerFar(setupStruct.currentBlock-(yyy+2*dxxx-1):setupStruct.currentBlock) > aecmStructNew.energyLevelMSE) + for xxx = -dxxx:dxxx + aecmStructNew.delayLatestS(xxx+dxxx+1) = sum(sign(aecmStructNew.enerEcho(setupStruct.currentBlock-(yyy+dxxx-xxx)+1:setupStruct.currentBlock+xxx-dxxx)-mean(aecmStructNew.enerEcho(setupStruct.currentBlock-(yyy++dxxx-xxx)+1:setupStruct.currentBlock+xxx-dxxx))).*sign(aecmStructNew.enerNear(setupStruct.currentBlock-yyy-dxxx+1:setupStruct.currentBlock-dxxx)-mean(aecmStructNew.enerNear(setupStruct.currentBlock-yyy-dxxx+1:setupStruct.currentBlock-dxxx)))); + end + aecmStructNew.newDelayCurve = 1; + end +end +if ((setupStruct.currentBlock > 2*aecmStructNew.convLength) & ~rem(setupStruct.currentBlock,yyy*2) & aecmStructNew.newDelayCurve) + [maxV,maxP] = max(aecmStructNew.delayLatestS); + if ((maxP > 2) & (maxP < 2*dxxx)) + maxVLeft = aecmStructNew.delayLatestS(max(1,maxP-4)); + maxVRight = aecmStructNew.delayLatestS(min(2*dxxx+1,maxP+4)); + %fprintf('Max %d, Left %d, Right %d\n',maxV,maxVLeft,maxVRight) % Uncomment if you want to see max value + if ((maxV > 24) & (maxVLeft < maxV - 10) & (maxVRight < maxV - 10)) + aecmStructNew.feedbackDelay = maxP-dxxx-1; + aecmStructNew.newDelayCurve = 0; + aecmStructNew.feedbackDelayUpdate = 1; + fprintf('Feedback Update at block %d\n',setupStruct.currentBlock) + end + end +end +% End of "Check delay" +%%%%%%%% + +%%%%% +% Calculate suppression gain, based on far end energy and near end energy error (dE) +if (setupStruct.supGain_flag) + [gamma_echo, aecmStructNew.cntIn, aecmStructNew.cntOut] = calcFilterGain(aecmStructNew.enerFar(setupStruct.currentBlock), dE, aecmStructNew, setupStruct.currentBlock, aecmStructNew.convLength, aecmStructNew.cntIn, aecmStructNew.cntOut); +else + gamma_echo = 1; +end +aecmStructNew.gammaLog(setupStruct.currentBlock) = gamma_echo; % Store the gain +gamma_use = gamma_echo; + +% Use the stored channel +U = aecmStructNew.HStored.*xSubBand; + +% compute Wiener filter and suppressor function +Iy = find(ySubBand); +subBandFilter = zeros(size(ySubBand)); +if (aecmStructNew.bandFactor == 2) + subBandFilter(Iy) = (1 - gamma_use*sqrt(U(Iy)./ySubBand(Iy))); % For Faller +else + subBandFilter(Iy) = (1 - gamma_use*(U(Iy)./ySubBand(Iy))); % For COV +end +ix0 = find(subBandFilter < 0); % bounding trick 1 +subBandFilter(ix0) = 0; +ix0 = find(subBandFilter > 1); % bounding trick 1 +subBandFilter(ix0) = 1; + +% Interpolate back to normal frequency bins if we use sub bands +if aecmStructNew.useSubBand + thefilter = interp1(setupStruct.centerFreq,subBandFilter,linspace(0,setupStruct.samplingfreq/2,setupStruct.hsupport1)','nearest'); + testfilter = interp1(setupStruct.centerFreq,subBandFilter,linspace(0,setupStruct.samplingfreq/2,1000),'nearest'); + thefilter(end) = subBandFilter(end); + + internalIndex = 1; + for kk=1:setupStruct.subBandLength+1 + internalIndex:internalIndex+setupStruct.numInBand(kk)-1; + thefilter(internalIndex:internalIndex+setupStruct.numInBand(kk)-1) = subBandFilter(kk); + internalIndex = internalIndex + setupStruct.numInBand(kk); + end +else + thefilter = subBandFilter; + testfilter = subBandFilter; +end + +% Bound the filter +ix0 = find(thefilter < setupStruct.de_echo_bound); % bounding trick 1 +thefilter(ix0) = setupStruct.de_echo_bound; % bounding trick 2 +ix0 = find(thefilter > 1); % bounding in reasonable range +thefilter(ix0) = 1; + +%%%% +% NLP +%%%% +thefmean = mean(thefilter(8:16)); +if (thefmean < 1) + disp(''); +end +aecmStructNew.runningfmean = setupStruct.nl_alpha*aecmStructNew.runningfmean + (1-setupStruct.nl_alpha)*thefmean; +slope0 = 1.0/(1.0 - setupStruct.nlSeverity); % +thegain = max(0.0, min(1.0, slope0*(aecmStructNew.runningfmean - setupStruct.nlSeverity))); +if ~setupStruct.nlp_flag + thegain = 1; +end +% END NONLINEARITY +thefilter = thegain*thefilter; + +%%%% +% The suppression +%%%% +femicrophone = fmicrophone .* thefilter; +% Store the output energy (used for plotting) +%aecmStructNew.enerOut(setupStruct.currentBlock) = log(abs(femicrophone)'*abs(femicrophone)); +aecmStructNew.enerOut(setupStruct.currentBlock) = log(sum(abs(femicrophone))); + +if aecmStructNew.plotIt + figure(13) + subplot(311) + %plot(n100,enerFar(n100),'b-',n100,enerNear(n100),'k--',n100,enerEcho(n100),'r-',[n100(1) n100(end)],[1 1]*vadThNew,'b:',[n100(1) n100(end)],[1 1]*((energyMax-energyMin)/4+energyMin),'r-.',[n100(1) n100(end)],[1 1]*vadNearThNew,'g:',[n100(1) n100(end)],[1 1]*energyMax,'r-.',[n100(1) n100(end)],[1 1]*energyMin,'r-.','LineWidth',2) + plot(n100,aecmStructNew.enerFar(n100),'b-',n100,aecmStructNew.enerNear(n100),'k--',n100,aecmStructNew.enerOut(n100),'r-.',n100,aecmStructNew.enerEcho(n100),'r-',n100,aecmStructNew.enerEchoStored(n100),'c-',[n100(1) n100(end)],[1 1]*((aecmStructNew.energyMax-aecmStructNew.energyMin)/4+aecmStructNew.energyMin),'g-.',[n100(1) n100(end)],[1 1]*aecmStructNew.energyMax,'g-.',[n100(1) n100(end)],[1 1]*aecmStructNew.energyMin,'g-.','LineWidth',2) + %title(['Frame ',int2str(i),' av ',int2str(setupStruct.updateno),' State = ',int2str(speechState),' \mu = ',num2str(mu)]) + title(['\gamma = ',num2str(gamma_echo),' \mu = ',num2str(mu)]) + subplot(312) + %plot(n100,enerError,'b-',[n100(1) n100(end)],[1 1]*vadNearTh,'r:',[n100(1) n100(end)],[-1.5 -1.5]*vadNearTh,'r:','LineWidth',2) + %plot(n100,enerError,'b-',[n100(1) n100(end)],[1 1],'r:',[n100(1) n100(end)],[-2 -2],'r:','LineWidth',2) + plot(n100,enerError,'b-',n100,enerErrorStored,'c-',[n100(1) n100(end)],[1 1]*aecmStructNew.varMean,'k--',[n100(1) n100(end)],[1 1],'r:',[n100(1) n100(end)],[-2 -2],'r:','LineWidth',2) + % Plot mu + %plot(n100,log2(aecmStructNew.muLog(n100)),'b-','LineWidth',2) + %plot(n100,log2(aecmStructNew.HGain(n100)),'b-',[n100(1) n100(end)],[1 1]*log2(sum(aecmStructNew.HStored)),'r:','LineWidth',2) + title(['Block ',int2str(setupStruct.currentBlock),' av ',int2str(setupStruct.updateno)]) + subplot(313) + %plot(n100,enerVar(n100),'b-',[n100(1) n100(end)],[1 1],'r:',[n100(1) n100(end)],[-2 -2],'r:','LineWidth',2) + %plot(n100,enerVar(n100),'b-','LineWidth',2) + % Plot correlation curve + + %plot(-25:25,aecmStructNew.delayStored/max(aecmStructNew.delayStored),'c-',-25:25,aecmStructNew.delayLatest/max(aecmStructNew.delayLatest),'r-',-25:25,(max(aecmStructNew.delayStoredS)-aecmStructNew.delayStoredS)/(max(aecmStructNew.delayStoredS)-min(aecmStructNew.delayStoredS)),'c:',-25:25,(max(aecmStructNew.delayLatestS)-aecmStructNew.delayLatestS)/(max(aecmStructNew.delayLatestS)-min(aecmStructNew.delayLatestS)),'r:','LineWidth',2) + %plot(-25:25,aecmStructNew.delayStored,'c-',-25:25,aecmStructNew.delayLatest,'r-',-25:25,(max(aecmStructNew.delayStoredS)-aecmStructNew.delayStoredS)/(max(aecmStructNew.delayStoredS)-min(aecmStructNew.delayStoredS)),'c:',-25:25,(max(aecmStructNew.delayLatestS)-aecmStructNew.delayLatestS)/(max(aecmStructNew.delayLatestS)-min(aecmStructNew.delayLatestS)),'r:','LineWidth',2) + %plot(-25:25,aecmStructNew.delayLatest,'r-',-25:25,(50-aecmStructNew.delayLatestS)/100,'r:','LineWidth',2) + plot(-25:25,aecmStructNew.delayLatestS,'r:','LineWidth',2) + %plot(-25:25,aecmStructNew.delayStored,'c-',-25:25,aecmStructNew.delayLatest,'r-','LineWidth',2) + plot(0:32,aecmStruct.HStored,'bo-','LineWidth',2) + %title(['\gamma | In = ',int2str(aecmStructNew.muStruct.countInInterval),' | Out High = ',int2str(aecmStructNew.muStruct.countOutHighInterval),' | Out Low = ',int2str(aecmStructNew.muStruct.countOutLowInterval)]) + pause(1) + %if ((setupStruct.currentBlock == 860) | (setupStruct.currentBlock == 420) | (setupStruct.currentBlock == 960)) + if 0%(setupStruct.currentBlock == 960) + figure(60) + plot(n100,aecmStructNew.enerNear(n100),'k--',n100,aecmStructNew.enerEcho(n100),'k:','LineWidth',2) + legend('Near End','Estimated Echo') + title('Signal Energy witH offset compensation') + figure(61) + subplot(211) + stem(sign(aecmStructNew.enerNear(n100)-mean(aecmStructNew.enerNear(n100)))) + title('Near End Energy Pattern (around mean value)') + subplot(212) + stem(sign(aecmStructNew.enerEcho(n100)-mean(aecmStructNew.enerEcho(n100)))) + title('Estimated Echo Energy Pattern (around mean value)') + pause + end + drawnow%,pause +elseif ~rem(setupStruct.currentBlock,100) + fprintf('Block %d of %d\n',setupStruct.currentBlock,setupStruct.updateno) +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/align.m b/src/modules/audio_processing/aecm/main/matlab/matlab/align.m new file mode 100644 index 0000000000..9b9c0baf3b --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/align.m @@ -0,0 +1,98 @@ +function [delayStructNew] = align(xf, yf, delayStruct, i, trueDelay); + +%%%%%%% +% Bastiaan's algorithm copied +%%%%%%% +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]; + +oldMethod = 1; % Turn on or off the old method. The new one is Bastiaan's August 2008 updates +THReSHoLD = 2.0; % ADJUSTABLE threshold factor; 4.0 seems good +%%%%%%%%%%%%%%%%%%% +% use log domain (showed improved performance) +xxf = sqrt(real(xf.*conj(xf))+1e-20); +yyf = sqrt(real(yf.*conj(yf))+1e-20); +delayStruct.sxAll2(:,i) = 20*log10(xxf); +delayStruct.syAll2(:,i) = 20*log10(yyf); + +mD = min(i-1,delayStruct.maxDelayb); +if oldMethod + factor = 1.0; + histLenb = 250; + xthreshold = factor*median(delayStruct.sxAll2(:,i-mD:i),2); + ythreshold = factor*median(delayStruct.syAll2(:,i-mD:i),2); +else + xthreshold = sum(delayStruct.sxAll2(:,i-mD:i),2)/(delayStruct.maxDelayb+1); + + [yout, delayStruct.z200] = filter(Bp200, Ap200, delayStruct.syAll2(:,i), delayStruct.z200, 2); + yout = yout/(delayStruct.maxDelayb+1); + ythreshold = mean(delayStruct.syAll2(:,i-mD:i),2); + ythreshold = yout; +end + +delayStruct.bxspectrum(i) = getBspectrum(delayStruct.sxAll2(:,i), xthreshold, delayStruct.bandfirst, delayStruct.bandlast); +delayStruct.byspectrum(i) = getBspectrum(delayStruct.syAll2(:,i), ythreshold, delayStruct.bandfirst, delayStruct.bandlast); + +delayStruct.bxhist(end-mD:end) = delayStruct.bxspectrum(i-mD:i); + +delayStruct.bcount(:,i) = hisser2(delayStruct.byspectrum(i), flipud(delayStruct.bxhist), delayStruct.bandfirst, delayStruct.bandlast); +[delayStruct.fout(:,i), delayStruct.z500] = filter(Bp500, Ap500, delayStruct.bcount(:,i), delayStruct.z500, 2); +if oldMethod + %delayStruct.new(:,i) = sum(delayStruct.bcount(:,max(1,i-histLenb+1):i),2); % using the history range + tmpVec = [delayStruct.fout(1,i)*ones(2,1); delayStruct.fout(:,i); delayStruct.fout(end,i)*ones(2,1)]; % using the history range + tmpVec = filter(ones(1,5), 1, tmpVec); + delayStruct.new(:,i) = tmpVec(5:end); + %delayStruct.new(:,i) = delayStruct.fout(:,i); % using the history range +else + [delayStruct.fout(:,i), delayStruct.z500] = filter(Bp500, Ap500, delayStruct.bcount(:,i), delayStruct.z500, 2); + % NEW CODE + delayStruct.new(:,i) = filter([-1,-2,1,4,1,-2,-1], 1, delayStruct.fout(:,i)); %remv smth component + delayStruct.new(1:end-3,i) = delayStruct.new(1+3:end,i); + delayStruct.new(1:6,i) = 0.0; + delayStruct.new(end-6:end,i) = 0.0; % ends are no good +end +[valuen, tempdelay] = min(delayStruct.new(:,i)); % find minimum +if oldMethod + threshold = valuen + (max(delayStruct.new(:,i)) - valuen)/4; + thIndex = find(delayStruct.new(:,i) <= threshold); + if (i > 1) + delayDiff = abs(delayStruct.delay(i-1)-tempdelay+1); + if (delayStruct.oneGoodEstimate & (max(diff(thIndex)) > 1) & (delayDiff < 10)) + % We consider this minimum to be significant, hence update the delay + delayStruct.delay(i) = tempdelay; + elseif (~delayStruct.oneGoodEstimate & (max(diff(thIndex)) > 1)) + delayStruct.delay(i) = tempdelay; + if (i > histLenb) + delayStruct.oneGoodEstimate = 1; + end + else + delayStruct.delay(i) = delayStruct.delay(i-1); + end + else + delayStruct.delay(i) = tempdelay; + end +else + threshold = THReSHoLD*std(delayStruct.new(:,i)); % set updata threshold + if ((-valuen > threshold) | (i < delayStruct.smlength)) % see if you want to update delay + delayStruct.delay(i) = tempdelay; + else + delayStruct.delay(i) = delayStruct.delay(i-1); + end + % END NEW CODE +end +delayStructNew = delayStruct; + +% administrative and plotting stuff +if( 0) + figure(10); + plot([1:length(delayStructNew.new(:,i))],delayStructNew.new(:,i),trueDelay*[1 1],[min(delayStructNew.new(:,i)),max(delayStructNew.new(:,i))],'r',[1 length(delayStructNew.new(:,i))],threshold*[1 1],'r:', 'LineWidth',2); + %plot([1:length(delayStructNew.bcount(:,i))],delayStructNew.bcount(:,i),trueDelay*[1 1],[min(delayStructNew.bcount(:,i)),max(delayStructNew.bcount(:,i))],'r','LineWidth',2); + %plot([thedelay,thedelay],[min(fcount(:,i)),max(fcount(:,i))],'r'); + %title(sprintf('bin count and known delay at time %5.1f s\n',(i-1)*(support/(fs*oversampling)))); + title(delayStructNew.oneGoodEstimate) + xlabel('delay in frames'); + %hold off; + drawnow +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/calcFilterGain.m b/src/modules/audio_processing/aecm/main/matlab/matlab/calcFilterGain.m new file mode 100644 index 0000000000..a09a7f2225 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/calcFilterGain.m @@ -0,0 +1,88 @@ +function [gam, cntIn2, cntOut2] = calcFilterGain(energy, dE, aecmStruct, t, T, cntIn, cntOut) + +defaultLevel = 1.2; +cntIn2 = cntIn; +cntOut2 = cntOut; +if (t < T) + gam = 1; +else + dE1 = -5; + dE2 = 1; + gamMid = 0.2; + gam = max(0,min((energy - aecmStruct.energyMin)/(aecmStruct.energyLevel - aecmStruct.energyMin), 1-(1-gamMid)*(aecmStruct.energyMax-energy)/(aecmStruct.energyMax-aecmStruct.energyLevel))); + + dEOffset = -0.5; + dEWidth = 1.5; + %gam2 = max(1,2-((dE-dEOffset)/(dE2-dEOffset)).^2); + gam2 = 1+(abs(dE-dEOffset)<(dE2-dEOffset)); + + gam = gam*gam2; + + + if (energy < aecmStruct.energyLevel) + gam = 0; + else + gam = defaultLevel; + end + dEVec = aecmStruct.enerNear(t-63:t)-aecmStruct.enerEcho(t-63:t); + %dEVec = aecmStruct.enerNear(t-20:t)-aecmStruct.enerEcho(t-20:t); + numCross = 0; + currentState = 0; + for ii=1:64 + if (currentState == 0) + currentState = (dEVec(ii) > dE2) - (dEVec(ii) < -2); + elseif ((currentState == 1) & (dEVec(ii) < -2)) + numCross = numCross + 1; + currentState = -1; + elseif ((currentState == -1) & (dEVec(ii) > dE2)) + numCross = numCross + 1; + currentState = 1; + end + end + gam = max(0, gam - numCross/25); + gam = 1; + + ener_A = 1; + ener_B = 0.8; + ener_C = aecmStruct.energyLevel + (aecmStruct.energyMax-aecmStruct.energyLevel)/5; + dE_A = 4;%2; + dE_B = 3.6;%1.8; + dE_C = 0.9*dEWidth; + dE_D = 1; + timeFactorLength = 10; + ddE = abs(dE-dEOffset); + if (energy < aecmStruct.energyLevel) + gam = 0; + else + gam = 1; + gam2 = max(0, min(ener_B*(energy-aecmStruct.energyLevel)/(ener_C-aecmStruct.energyLevel), ener_B+(ener_A-ener_B)*(energy-ener_C)/(aecmStruct.energyMax-ener_C))); + if (ddE < dEWidth) + % Update counters + cntIn2 = cntIn2 + 1; + if (cntIn2 > 2) + cntOut2 = 0; + end + gam3 = max(dE_D, min(dE_A-(dE_A-dE_B)*(ddE/dE_C), dE_D+(dE_B-dE_D)*(dEWidth-ddE)/(dEWidth-dE_C))); + gam3 = dE_A; + else + % Update counters + cntOut2 = cntOut2 + 1; + if (cntOut2 > 2) + cntIn2 = 0; + end + %gam2 = 1; + gam3 = dE_D; + end + timeFactor = min(1, cntIn2/timeFactorLength); + gam = gam*(1-timeFactor) + timeFactor*gam2*gam3; + end + %gam = gam/floor(numCross/2+1); +end +if isempty(gam) + numCross + timeFactor + cntIn2 + cntOut2 + gam2 + gam3 +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/calcStepSize.m b/src/modules/audio_processing/aecm/main/matlab/matlab/calcStepSize.m new file mode 100644 index 0000000000..ae1365fa48 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/calcStepSize.m @@ -0,0 +1,105 @@ +function [mu, aecmStructNew] = calcStepSize(energy, dE, aecmStruct, t, logscale) + +if (nargin < 4) + t = 1; + logscale = 1; +elseif (nargin == 4) + logscale = 1; +end +T = aecmStruct.convLength; + +if logscale + currentMuMax = aecmStruct.MU_MIN + (aecmStruct.MU_MAX-aecmStruct.MU_MIN)*min(t,T)/T; + if (aecmStruct.energyMin >= aecmStruct.energyMax) + mu = aecmStruct.MU_MIN; + else + mu = (energy - aecmStruct.energyMin)/(aecmStruct.energyMax - aecmStruct.energyMin)*(currentMuMax-aecmStruct.MU_MIN) + aecmStruct.MU_MIN; + end + mu = 2^mu; + if (energy < aecmStruct.energyLevel) + mu = 0; + end +else + muMin = 0; + muMax = 0.5; + currentMuMax = muMin + (muMax-muMin)*min(t,T)/T; + if (aecmStruct.energyMin >= aecmStruct.energyMax) + mu = muMin; + else + mu = (energy - aecmStruct.energyMin)/(aecmStruct.energyMax - aecmStruct.energyMin)*(currentMuMax-muMin) + muMin; + end +end +dE2 = 1; +dEOffset = -0.5; +offBoost = 5; +if (mu > 0) + if (abs(dE-aecmStruct.ENERGY_DEV_OFFSET) > aecmStruct.ENERGY_DEV_TOL) + aecmStruct.muStruct.countInInterval = 0; + else + aecmStruct.muStruct.countInInterval = aecmStruct.muStruct.countInInterval + 1; + end + if (dE < aecmStruct.ENERGY_DEV_OFFSET - aecmStruct.ENERGY_DEV_TOL) + aecmStruct.muStruct.countOutLowInterval = aecmStruct.muStruct.countOutLowInterval + 1; + else + aecmStruct.muStruct.countOutLowInterval = 0; + end + if (dE > aecmStruct.ENERGY_DEV_OFFSET + aecmStruct.ENERGY_DEV_TOL) + aecmStruct.muStruct.countOutHighInterval = aecmStruct.muStruct.countOutHighInterval + 1; + else + aecmStruct.muStruct.countOutHighInterval = 0; + end +end +muVar = 2^min(-3,5/50*aecmStruct.muStruct.countInInterval-3); +muOff = 2^max(offBoost,min(0,offBoost*(aecmStruct.muStruct.countOutLowInterval-aecmStruct.muStruct.minOutLowInterval)/(aecmStruct.muStruct.maxOutLowInterval-aecmStruct.muStruct.minOutLowInterval))); + +muLow = 1/64; +muVar = 1; +if (t < 2*T) + muDT = 1; + muVar = 1; + mdEVec = 0; + numCross = 0; +else + muDT = min(1,max(muLow,1-(1-muLow)*(dE-aecmStruct.ENERGY_DEV_OFFSET)/aecmStruct.ENERGY_DEV_TOL)); + dEVec = aecmStruct.enerNear(t-63:t)-aecmStruct.enerEcho(t-63:t); + %dEVec = aecmStruct.enerNear(t-20:t)-aecmStruct.enerEcho(t-20:t); + numCross = 0; + currentState = 0; + for ii=1:64 + if (currentState == 0) + currentState = (dEVec(ii) > dE2) - (dEVec(ii) < -2); + elseif ((currentState == 1) & (dEVec(ii) < -2)) + numCross = numCross + 1; + currentState = -1; + elseif ((currentState == -1) & (dEVec(ii) > dE2)) + numCross = numCross + 1; + currentState = 1; + end + end + + %logicDEVec = (dEVec > dE2) - (dEVec < -2); + %numCross = sum(abs(diff(logicDEVec))); + %mdEVec = mean(abs(dEVec-dEOffset)); + %mdEVec = mean(abs(dEVec-mean(dEVec))); + %mdEVec = max(dEVec)-min(dEVec); + %if (mdEVec > 4)%1.5) + % muVar = 0; + %end + muVar = 2^(-floor(numCross/2)); + muVar = 2^(-numCross); +end +%muVar = 1; + + +% if (eStd > (dE2-dEOffset)) +% muVar = 1/8; +% else +% muVar = 1; +% end + +%mu = mu*muDT*muVar*muOff; +mu = mu*muDT*muVar; +mu = min(mu,0.25); +aecmStructNew = aecmStruct; +%aecmStructNew.varMean = mdEVec; +aecmStructNew.varMean = numCross; diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/fallerEstimator.m b/src/modules/audio_processing/aecm/main/matlab/matlab/fallerEstimator.m new file mode 100644 index 0000000000..d038b519c0 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/fallerEstimator.m @@ -0,0 +1,42 @@ +function [U, Hnew] = fallerEstimator(Y, X, H, mu) + +% Near end signal is stacked frame by frame columnwise in matrix Y and far end in X +% +% Possible estimation procedures are +% 1) LSE +% 2) NLMS +% 3) Separated numerator and denomerator filters +regParam = 1; +[numFreqs, numFrames] = size(Y); +[numFreqs, Q] = size(X); +U = zeros(numFreqs, 1); + +if ((nargin == 3) | (nargin == 5)) + dtd = 0; +end +if (nargin == 4) + dtd = H; +end +Emax = 7; +dEH = Emax-sum(sum(H)); +nu = 2*mu; +% if (nargin < 5) +% H = zeros(numFreqs, Q); +% for kk = 1:numFreqs +% Xmatrix = hankel(X(kk,1:Q),X(kk,Q:end)); +% y = Y(kk,1:end-Q+1)'; +% H(kk,:) = (y'*Xmatrix')*inv(Xmatrix*Xmatrix'+regParam); +% U(kk,1) = H(kk,:)*Xmatrix(:,1); +% end +% else + for kk = 1:numFreqs + x = X(kk,1:Q)'; + y = Y(kk,1); + Htmp = mu*(y-H(kk,:)*x)/(x'*x+regParam)*x; + %Htmp = (mu*(y-H(kk,:)*x)/(x'*x+regParam) - nu/dEH)*x; + H(kk,:) = H(kk,:) + Htmp'; + U(kk,1) = H(kk,:)*x; + end +% end + +Hnew = H; diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/getBspectrum.m b/src/modules/audio_processing/aecm/main/matlab/matlab/getBspectrum.m new file mode 100644 index 0000000000..a4a533d600 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/getBspectrum.m @@ -0,0 +1,22 @@ +function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast) +% function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast) +% compute binary spectrum using threshold spectrum as pivot +% bspectrum = binary spectrum (binary) +% ps=current power spectrum (float) +% threshold=threshold spectrum (float) +% bandfirst = first band considered +% bandlast = last band considered + +% initialization stuff + if( length(ps)<bandlast | bandlast>32 | length(ps)~=length(threshold)) + error('BinDelayEst:spectrum:invalid','Dimensionality error'); +end + +% get current binary spectrum +diff = ps - threshold; +bspectrum=uint32(0); +for(i=bandfirst:bandlast) + if( diff(i)>0 ) + bspectrum = bitset(bspectrum,i); + end +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/hisser2.m b/src/modules/audio_processing/aecm/main/matlab/matlab/hisser2.m new file mode 100644 index 0000000000..5a414f9da8 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/hisser2.m @@ -0,0 +1,21 @@ +function bcount=hisser2(bs,bsr,bandfirst,bandlast) +% function bcount=hisser(bspectrum,bandfirst,bandlast) +% histogram for the binary spectra +% bcount= array of bit counts +% bs=binary spectrum (one int32 number each) +% bsr=reference binary spectra (one int32 number each) +% blockSize = histogram over blocksize blocks +% bandfirst = first band considered +% bandlast = last band considered + +% weight all delays equally +maxDelay = length(bsr); + +% compute counts (two methods; the first works better and is operational) +bcount=zeros(maxDelay,1); +for(i=1:maxDelay) + % the delay should have low count for low-near&high-far and high-near&low-far + bcount(i)= sum(bitget(bitxor(bs,bsr(i)),bandfirst:bandlast)); + % the delay should have low count for low-near&high-far (works less well) +% bcount(i)= sum(bitget(bitand(bsr(i),bitxor(bs,bsr(i))),bandfirst:bandlast)); +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/mainProgram.m b/src/modules/audio_processing/aecm/main/matlab/matlab/mainProgram.m new file mode 100644 index 0000000000..eeb2aaa79c --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/mainProgram.m @@ -0,0 +1,283 @@ +useHTC = 1; % Set this if you want to run a single file and set file names below. Otherwise use simEnvironment to run from several scenarios in a row +delayCompensation_flag = 0; % Set this flag to one if you want to turn on the delay compensation/enhancement +global FARENDFFT; +global NEARENDFFT; +global F; + +if useHTC +% fid=fopen('./htcTouchHd/nb/aecFar.pcm'); xFar=fread(fid,'short'); fclose(fid); +% fid=fopen('./htcTouchHd/nb/aecNear.pcm'); yNear=fread(fid,'short'); fclose(fid); +% fid=fopen('./samsungBlackjack/nb/aecFar.pcm'); xFar=fread(fid,'short'); fclose(fid); +% fid=fopen('./samsungBlackjack/nb/aecNear.pcm'); yNear=fread(fid,'short'); fclose(fid); +% fid=fopen('aecFarPoor.pcm'); xFar=fread(fid,'short'); fclose(fid); +% fid=fopen('aecNearPoor.pcm'); yNear=fread(fid,'short'); fclose(fid); +% fid=fopen('out_aes.pcm'); outAES=fread(fid,'short'); fclose(fid); + fid=fopen('aecFar4.pcm'); xFar=fread(fid,'short'); fclose(fid); + fid=fopen('aecNear4.pcm'); yNear=fread(fid,'short'); fclose(fid); + yNearSpeech = zeros(size(xFar)); + fs = 8000; + frameSize = 64; +% frameSize = 128; + fs = 16000; +% frameSize = 256; +%F = load('fftValues.txt'); +%FARENDFFT = F(:,1:33); +%NEARENDFFT = F(:,34:66); + +else + loadFileFar = [speakerType, '_s_',scenario,'_far_b.wav']; + [xFar,fs,nbits] = wavread(loadFileFar); + xFar = xFar*2^(nbits-1); + loadFileNear = [speakerType, '_s_',scenario,'_near_b.wav']; + [yNear,fs,nbits] = wavread(loadFileNear); + yNear = yNear*2^(nbits-1); + loadFileNearSpeech = [speakerType, '_s_',scenario,'_nearSpeech_b.wav']; + [yNearSpeech,fs,nbits] = wavread(loadFileNearSpeech); + yNearSpeech = yNearSpeech*2^(nbits-1); + frameSize = 256; +end + +dtRegions = []; + +% General settings for the AECM +setupStruct = struct(... + 'stepSize_flag', 1,... % This flag turns on the step size calculation. If turned off, mu = 0.25. + 'supGain_flag', 0,... % This flag turns on the suppression gain calculation. If turned off, gam = 1. + 'channelUpdate_flag', 0,... % This flag turns on the channel update. If turned off, H is updated for convLength and then kept constant. + 'nlp_flag', 0,... % Turn on/off NLP + 'withVAD_flag', 0,... % Turn on/off NLP + 'useSubBand', 0,... % Set to 1 if to use subBands + 'useDelayEstimation', 1,... % Set to 1 if to use delay estimation + 'support', frameSize,... % # of samples per frame + 'samplingfreq',fs,... % Sampling frequency + 'oversampling', 2,... % Overlap between blocks/frames + 'updatel', 0,... % # of samples between blocks + 'hsupport1', 0,... % # of bins in frequency domain + 'factor', 0,... % synthesis window amplification + 'tlength', 0,... % # of samples of entire file + 'updateno', 0,... % # of updates + 'nb', 1,... % # of blocks + 'currentBlock', 0,... % + 'win', zeros(frameSize,1),...% Window to apply for fft and synthesis + 'avtime', 1,... % Time (in sec.) to perform averaging + 'estLen', 0,... % Averaging in # of blocks + 'A_GAIN', 10.0,... % + '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,... % + 'nl_alpha', 0.4,... % memory; seems not very critical + 'nlSeverity', 0.2,... % nonlinearity severity: 0 does nothing; 1 suppresses all + 'numInBand', [],... % # of frequency bins in resp. subBand + 'centerFreq', [],... % Center frequency of resp. subBand + 'dtRegions', dtRegions,... % Regions where we have DT + 'subBandLength', frameSize/2);%All bins + %'subBandLength', 11); %Something's wrong when subBandLength even + %'nl_alpha', 0.8,... % memory; seems not very critical + +delayStruct = struct(... + 'bandfirst', 8,... + 'bandlast', 25,... + 'smlength', 600,... + 'maxDelay', 0.4,... + 'oneGoodEstimate', 0,... + 'delayAdjust', 0,... + 'maxDelayb', 0); +% More parameters in delayStruct are constructed in "updateSettings" below + +% Make struct settings +[setupStruct, delayStruct] = updateSettings(yNear, xFar, setupStruct, delayStruct); +setupStruct.numInBand = ones(setupStruct.hsupport1,1); + +Q = 1; % Time diversity in channel +% General settings for the step size calculation +muStruct = struct(... + 'countInInterval', 0,... + 'countOutHighInterval', 0,... + 'countOutLowInterval', 0,... + 'minInInterval', 50,... + 'minOutHighInterval', 10,... + 'minOutLowInterval', 10,... + 'maxOutLowInterval', 50); +% General settings for the AECM +aecmStruct = struct(... + 'plotIt', 0,... % Set to 0 to turn off plotting + 'useSubBand', 0,... + 'bandFactor', 1,... + 'H', zeros(setupStruct.subBandLength+1,Q),... + 'HStored', zeros(setupStruct.subBandLength+1,Q),... + 'X', zeros(setupStruct.subBandLength+1,Q),... + 'energyThres', 0.28,... + 'energyThresMSE', 0.4,... + 'energyMin', inf,... + 'energyMax', -inf,... + 'energyLevel', 0,... + 'energyLevelMSE', 0,... + 'convLength', 100,... + 'gammaLog', ones(setupStruct.updateno,1),... + 'muLog', ones(setupStruct.updateno,1),... + 'enerFar', zeros(setupStruct.updateno,1),... + 'enerNear', zeros(setupStruct.updateno,1),... + 'enerEcho', zeros(setupStruct.updateno,1),... + 'enerEchoStored', zeros(setupStruct.updateno,1),... + 'enerOut', zeros(setupStruct.updateno,1),... + 'runningfmean', 0,... + 'muStruct', muStruct,... + 'varMean', 0,... + 'countMseH', 0,... + 'mseHThreshold', 1.1,... + 'mseHStoredOld', inf,... + 'mseHLatestOld', inf,... + 'delayLatestS', zeros(1,51),... + 'feedbackDelay', 0,... + 'feedbackDelayUpdate', 0,... + 'cntIn', 0,... + 'cntOut', 0,... + 'FAR_ENERGY_MIN', 1,... + 'ENERGY_DEV_OFFSET', 0.5,... + 'ENERGY_DEV_TOL', 1.5,... + 'MU_MIN', -16,... + 'MU_MAX', -2,... + 'newDelayCurve', 0); + +% Adjust speech signals +xFar = [zeros(setupStruct.hsupport1-1,1);xFar(1:setupStruct.tlength)]; +yNear = [zeros(setupStruct.hsupport1-1,1);yNear(1:setupStruct.tlength)]; +yNearSpeech = [zeros(setupStruct.hsupport1-1,1);yNearSpeech(1:setupStruct.tlength)]; +xFar = xFar(1:setupStruct.tlength); +yNear = yNear(1:setupStruct.tlength); + +% Set figure settings +if aecmStruct.plotIt + figure(13) + set(gcf,'doublebuffer','on') +end +%%%%%%%%%% +% Here starts the algorithm +% Dividing into frames and then estimating the near end speech +%%%%%%%%%% +fTheFarEnd = complex(zeros(setupStruct.hsupport1,1)); +afTheFarEnd = zeros(setupStruct.hsupport1,setupStruct.updateno+1); +fFar = zeros(setupStruct.hsupport1,setupStruct.updateno+1); +fmicrophone = complex(zeros(setupStruct.hsupport1,1)); +afmicrophone = zeros(setupStruct.hsupport1,setupStruct.updateno+1); +fNear = zeros(setupStruct.hsupport1,setupStruct.updateno+1); +femicrophone = complex(zeros(setupStruct.hsupport1,1)); +emicrophone = zeros(setupStruct.tlength,1); + +if (setupStruct.useDelayEstimation == 2) + delSamples = [1641 1895 2032 1895 2311 2000 2350 2222 NaN 2332 2330 2290 2401 2415 NaN 2393 2305 2381 2398]; + delBlocks = round(delSamples/setupStruct.updatel); + delStarts = floor([25138 46844 105991 169901 195739 218536 241803 333905 347703 362660 373753 745135 765887 788078 806257 823835 842443 860139 881869]/setupStruct.updatel); +else + delStarts = []; +end + +for i=1:setupStruct.updateno + setupStruct.currentBlock = i; + + sb = (i-1)*setupStruct.updatel + 1; + se = sb + setupStruct.support - 1; + + %%%%%%% + % Analysis FFTs + %%%%%%% + % Far end signal + temp = fft(setupStruct.win .* xFar(sb:se))/frameSize; + fTheFarEnd = temp(1:setupStruct.hsupport1); + afTheFarEnd(:,i) = abs(fTheFarEnd); + fFar(:,i) = fTheFarEnd; + % Near end signal + temp = fft(setupStruct.win .* yNear(sb:se))/frameSize;%,pause + fmicrophone = temp(1:setupStruct.hsupport1); + afmicrophone(:,i) = abs(fmicrophone); + fNear(:,i) = fmicrophone; + %abs(fmicrophone),pause + % The true near end speaker (if we have such info) + temp = fft(setupStruct.win .* yNearSpeech(sb:se)); + aftrueSpeech = abs(temp(1:setupStruct.hsupport1)); + + if(i == 1000) + %break; + end + + % Perform delay estimation + if (setupStruct.useDelayEstimation == 1) + % Delay Estimation + delayStruct = align(fTheFarEnd, fmicrophone, delayStruct, i); + %delayStruct.delay(i) = 39;%19; + idel = max(i - delayStruct.delay(i) + 1,1); + + if delayCompensation_flag + % If we have a new delay estimate from Bastiaan's alg. update the offset + if (delayStruct.delay(i) ~= delayStruct.delay(max(1,i-1))) + delayStruct.delayAdjust = delayStruct.delayAdjust + delayStruct.delay(i) - delayStruct.delay(i-1); + end + % Store the compensated delay + delayStruct.delayNew(i) = delayStruct.delay(i) - delayStruct.delayAdjust; + if (delayStruct.delayNew(i) < 1) + % Something's wrong + pause,break + end + % Compensate with the offset estimate + idel = idel + delayStruct.delayAdjust; + end + if 0%aecmStruct.plotIt + figure(1) + plot(1:i,delayStruct.delay(1:i),'k:',1:i,delayStruct.delayNew(1:i),'k--','LineWidth',2),drawnow + end + elseif (setupStruct.useDelayEstimation == 2) + % Use "manual delay" + delIndex = find(delStarts<i); + if isempty(delIndex) + idel = i; + else + idel = i - delBlocks(max(delIndex)); + if isnan(idel) + idel = i - delBlocks(max(delIndex)-1); + end + end + else + % No delay estimation + %idel = max(i - 18, 1); + idel = max(i - 50, 1); + end + + %%%%%%%% + % This is the AECM algorithm + % + % Output is the new frequency domain signal (hopefully) echo compensated + %%%%%%%% + [femicrophone, aecmStruct] = AECMobile(fmicrophone, afTheFarEnd(:,idel), setupStruct, aecmStruct); + %[femicrophone, aecmStruct] = AECMobile(fmicrophone, FARENDFFT(idel,:)'/2^F(idel,end-1), setupStruct, aecmStruct); + + if aecmStruct.feedbackDelayUpdate + % If the feedback tells us there is a new offset out there update the enhancement + delayStruct.delayAdjust = delayStruct.delayAdjust + aecmStruct.feedbackDelay; + aecmStruct.feedbackDelayUpdate = 0; + end + + % reconstruction; first make spectrum odd + temp = [femicrophone; flipud(conj(femicrophone(2:(setupStruct.hsupport1-1))))]; + emicrophone(sb:se) = emicrophone(sb:se) + setupStruct.factor * setupStruct.win .* real(ifft(temp))*frameSize; + if max(isnan(emicrophone(sb:se))) + % Something's wrong with the output at block i + i + break + end +end + + +if useHTC + fid=fopen('aecOutMatlabC.pcm','w');fwrite(fid,int16(emicrophone),'short');fclose(fid); + %fid=fopen('farendFFT.txt','w');fwrite(fid,int16(afTheFarEnd(:)),'short');fclose(fid); + %fid=fopen('farendFFTreal.txt','w');fwrite(fid,int16(imag(fFar(:))),'short');fclose(fid); + %fid=fopen('farendFFTimag.txt','w');fwrite(fid,int16(real(fFar(:))),'short');fclose(fid); + %fid=fopen('nearendFFT.txt','w');fwrite(fid,int16(afmicrophone(:)),'short');fclose(fid); + %fid=fopen('nearendFFTreal.txt','w');fwrite(fid,int16(real(fNear(:))),'short');fclose(fid); + %fid=fopen('nearendFFTimag.txt','w');fwrite(fid,int16(imag(fNear(:))),'short');fclose(fid); +end +if useHTC + %spclab(setupStruct.samplingfreq,xFar,yNear,emicrophone) +else + spclab(setupStruct.samplingfreq,xFar,yNear,emicrophone,yNearSpeech) +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/simEnvironment.m b/src/modules/audio_processing/aecm/main/matlab/matlab/simEnvironment.m new file mode 100644 index 0000000000..3ebe701dfd --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/simEnvironment.m @@ -0,0 +1,15 @@ +speakerType = 'fm'; +%for k=2:5 +%for k=[2 4 5] +for k=3 + scenario = int2str(k); + fprintf('Current scenario: %d\n',k) + mainProgram + %saveFile = [speakerType, '_s_',scenario,'_delayEst_v2_vad_man.wav']; + %wavwrite(emic,fs,nbits,saveFile); + %saveFile = ['P:\Engineering_share\BjornV\AECM\',speakerType, '_s_',scenario,'_delayEst_v2_vad_man.pcm']; + %saveFile = [speakerType, '_s_',scenario,'_adaptMu_adaptGamma_withVar_gammFilt_HSt.pcm']; + saveFile = ['scenario_',scenario,'_090417_backupH_nlp.pcm']; + fid=fopen(saveFile,'w');fwrite(fid,int16(emicrophone),'short');fclose(fid); + %pause +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/updateSettings.m b/src/modules/audio_processing/aecm/main/matlab/matlab/updateSettings.m new file mode 100644 index 0000000000..c805f1d09f --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/updateSettings.m @@ -0,0 +1,94 @@ +function [setupStructNew, delayStructNew] = updateSettings(microphone, TheFarEnd, setupStruct, delayStruct); + +% other, constants +setupStruct.hsupport1 = setupStruct.support/2 + 1; +setupStruct.factor = 2 / setupStruct.oversampling; +setupStruct.updatel = setupStruct.support/setupStruct.oversampling; +setupStruct.estLen = round(setupStruct.avtime * setupStruct.samplingfreq/setupStruct.updatel); + +% compute some constants +blockLen = setupStruct.support/setupStruct.oversampling; +delayStruct.maxDelayb = floor(setupStruct.samplingfreq*delayStruct.maxDelay/setupStruct.updatel); % in blocks + +%input +tlength = min([length(microphone),length(TheFarEnd)]); +updateno = floor(tlength/setupStruct.updatel); +setupStruct.tlength = setupStruct.updatel*updateno; +setupStruct.updateno = updateno - setupStruct.oversampling + 1; + +% signal length +n = floor(min([length(TheFarEnd), length(microphone)])/setupStruct.support)*setupStruct.support; +setupStruct.nb = n/blockLen - setupStruct.oversampling + 1; % in blocks + +setupStruct.win = sqrt([0 ; hanning(setupStruct.support-1)]); + +% Construct filterbank in Bark-scale + +K = setupStruct.subBandLength; %Something's wrong when K even +erbs = 21.4*log10(0.00437*setupStruct.samplingfreq/2+1); +fe = (10.^((0:K)'*erbs/K/21.4)-1)/0.00437; +setupStruct.centerFreq = fe; +H = diag(ones(1,K-1))+diag(ones(1,K-2),-1); +Hinv = inv(H); +aty = 2*Hinv(end,:)*fe(2:end-1); +boundary = aty - (setupStruct.samplingfreq/2 + fe(end-1))/2; +if rem(K,2) + x1 = min([fe(2)/2, -boundary]); +else + x1 = max([0, boundary]); +end +%x1 +g = fe(2:end-1); +g(1) = g(1) - x1/2; +x = 2*Hinv*g; +x = [x1;x]; +%figure(42), clf +xy = zeros((K+1)*4,1); +yy = zeros((K+1)*4,1); +xy(1:4) = [fe(1) fe(1) x(1) x(1)]'; +yy(1:4) = [0 1 1 0]'/x(1); +for kk=2:K + xy((kk-1)*4+(1:4)) = [x(kk-1) x(kk-1) x(kk) x(kk)]'; + yy((kk-1)*4+(1:4)) = [0 1 1 0]'/(x(kk)-x(kk-1)); +end +xy(end-3:end) = [x(K) x(K) fe(end) fe(end)]'; +yy(end-3:end) = [0 1 1 0]'/(fe(end)*2-2*x(K)); +%plot(xy,yy,'LineWidth',2) +%fill(xy,yy,'y') + +x = [0;x]; +xk = x*setupStruct.hsupport1/setupStruct.samplingfreq*2; +%setupStruct.erbBoundaries = xk; +numInBand = zeros(length(xk),1); +xh = (0:setupStruct.hsupport1-1); + +for kk=1:length(xk) + if (kk==length(xk)) + numInBand(kk) = length(find(xh>=xk(kk))); + else + numInBand(kk) = length(intersect(find(xh>=xk(kk)),find(xh<xk(kk+1)))); + end +end +setupStruct.numInBand = numInBand; + +setupStructNew = setupStruct; + +delayStructNew = struct(... + 'sxAll2',zeros(setupStructNew.hsupport1,setupStructNew.nb),... + 'syAll2',zeros(setupStructNew.hsupport1,setupStructNew.nb),... + 'z200',zeros(5,setupStructNew.hsupport1),... + 'z500',zeros(5,delayStruct.maxDelayb+1),... + 'bxspectrum',uint32(zeros(setupStructNew.nb,1)),... + 'byspectrum',uint32(zeros(setupStructNew.nb,1)),... + 'bandfirst',delayStruct.bandfirst,'bandlast',delayStruct.bandlast,... + 'bxhist',uint32(zeros(delayStruct.maxDelayb+1,1)),... + 'bcount',zeros(1+delayStruct.maxDelayb,setupStructNew.nb),... + 'fout',zeros(1+delayStruct.maxDelayb,setupStructNew.nb),... + 'new',zeros(1+delayStruct.maxDelayb,setupStructNew.nb),... + 'smlength',delayStruct.smlength,... + 'maxDelay', delayStruct.maxDelay,... + 'maxDelayb', delayStruct.maxDelayb,... + 'oneGoodEstimate', 0,... + 'delayAdjust', 0,... + 'delayNew',zeros(setupStructNew.nb,1),... + 'delay',zeros(setupStructNew.nb,1)); diff --git a/src/modules/audio_processing/aecm/main/matlab/waitbar_j.m b/src/modules/audio_processing/aecm/main/matlab/waitbar_j.m new file mode 100644 index 0000000000..50b9ccf309 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/waitbar_j.m @@ -0,0 +1,234 @@ +function fout = waitbar_j(x,whichbar, varargin) +%WAITBAR Display wait bar. +% H = WAITBAR(X,'title', property, value, property, value, ...) +% creates and displays a waitbar of fractional length X. The +% handle to the waitbar figure is returned in H. +% X should be between 0 and 1. Optional arguments property and +% value allow to set corresponding waitbar figure properties. +% Property can also be an action keyword 'CreateCancelBtn', in +% which case a cancel button will be added to the figure, and +% the passed value string will be executed upon clicking on the +% cancel button or the close figure button. +% +% WAITBAR(X) will set the length of the bar in the most recently +% created waitbar window to the fractional length X. +% +% WAITBAR(X,H) will set the length of the bar in waitbar H +% to the fractional length X. +% +% WAITBAR(X,H,'updated title') will update the title text in +% the waitbar figure, in addition to setting the fractional +% length to X. +% +% WAITBAR is typically used inside a FOR loop that performs a +% lengthy computation. A sample usage is shown below: +% +% h = waitbar(0,'Please wait...'); +% for i=1:100, +% % computation here % +% waitbar(i/100,h) +% end +% close(h) + +% Clay M. Thompson 11-9-92 +% Vlad Kolesnikov 06-7-99 +% Copyright 1984-2001 The MathWorks, Inc. +% $Revision: 1.22 $ $Date: 2001/04/15 12:03:29 $ + +if nargin>=2 + if ischar(whichbar) + type=2; %we are initializing + name=whichbar; + elseif isnumeric(whichbar) + type=1; %we are updating, given a handle + f=whichbar; + else + error(['Input arguments of type ' class(whichbar) ' not valid.']) + end +elseif nargin==1 + f = findobj(allchild(0),'flat','Tag','TMWWaitbar'); + + if isempty(f) + type=2; + name='Waitbar'; + else + type=1; + f=f(1); + end +else + error('Input arguments not valid.'); +end + +x = max(0,min(100*x,100)); + +switch type + case 1, % waitbar(x) update + p = findobj(f,'Type','patch'); + l = findobj(f,'Type','line'); + if isempty(f) | isempty(p) | isempty(l), + error('Couldn''t find waitbar handles.'); + end + xpatch = get(p,'XData'); + xpatch = [0 x x 0]; + set(p,'XData',xpatch) + xline = get(l,'XData'); + set(l,'XData',xline); + + if nargin>2, + % Update waitbar title: + hAxes = findobj(f,'type','axes'); + hTitle = get(hAxes,'title'); + set(hTitle,'string',varargin{1}); + end + + case 2, % waitbar(x,name) initialize + vertMargin = 0; + if nargin > 2, + % we have optional arguments: property-value pairs + if rem (nargin, 2 ) ~= 0 + error( 'Optional initialization arguments must be passed in pairs' ); + end + end + + oldRootUnits = get(0,'Units'); + + set(0, 'Units', 'points'); + screenSize = get(0,'ScreenSize'); + + axFontSize=get(0,'FactoryAxesFontSize'); + + pointsPerPixel = 72/get(0,'ScreenPixelsPerInch'); + + width = 360 * pointsPerPixel; + height = 75 * pointsPerPixel; + pos = [screenSize(3)/2-width/2 screenSize(4)/2-height/2 width height]; + +%pos= [501.75 589.5 393.75 52.5]; + f = figure(... + 'Units', 'points', ... + 'BusyAction', 'queue', ... + 'Position', pos, ... + 'Resize','on', ... + 'CreateFcn','', ... + 'NumberTitle','off', ... + 'IntegerHandle','off', ... + 'MenuBar', 'none', ... + 'Tag','TMWWaitbar',... + 'Interruptible', 'off', ... + 'Visible','on'); + + %%%%%%%%%%%%%%%%%%%%% + % set figure properties as passed to the fcn + % pay special attention to the 'cancel' request + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if nargin > 2, + propList = varargin(1:2:end); + valueList = varargin(2:2:end); + cancelBtnCreated = 0; + for ii = 1:length( propList ) + try + if strcmp(lower(propList{ii}), 'createcancelbtn' ) & ~cancelBtnCreated + cancelBtnHeight = 23 * pointsPerPixel; + cancelBtnWidth = 60 * pointsPerPixel; + newPos = pos; + vertMargin = vertMargin + cancelBtnHeight; + newPos(4) = newPos(4)+vertMargin; + callbackFcn = [valueList{ii}]; + set( f, 'Position', newPos, 'CloseRequestFcn', callbackFcn ); + cancelButt = uicontrol('Parent',f, ... + 'Units','points', ... + 'Callback',callbackFcn, ... + 'ButtonDownFcn', callbackFcn, ... + 'Enable','on', ... + 'Interruptible','off', ... + 'Position', [pos(3)-cancelBtnWidth*1.4, 7, ... + cancelBtnWidth, cancelBtnHeight], ... + 'String','Cancel', ... + 'Tag','TMWWaitbarCancelButton'); + cancelBtnCreated = 1; + else + % simply set the prop/value pair of the figure + set( f, propList{ii}, valueList{ii}); + end + catch + disp ( ['Warning: could not set property ''' propList{ii} ''' with value ''' num2str(valueList{ii}) '''' ] ); + end + end + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + colormap([]); + + axNorm=[.05 .3 .9 .2]; + % axNorm=[1 1 1 1]; + axPos=axNorm.*[pos(3:4),pos(3:4)] + [0 vertMargin 0 0]; + + h = axes('XLim',[0 100],... + 'YLim',[0 1],... + 'Box','on', ... + 'Units','Points',... + 'FontSize', axFontSize,... + 'Position',axPos,... + 'XTickMode','manual',... + 'YTickMode','manual',... + 'XTick',[],... + 'YTick',[],... + 'XTickLabelMode','manual',... + 'XTickLabel',[],... + 'YTickLabelMode','manual',... + 'YTickLabel',[]); + + tHandle=title(name); + tHandle=get(h,'title'); + oldTitleUnits=get(tHandle,'Units'); + set(tHandle,... + 'Units', 'points',... + 'String', name); + + tExtent=get(tHandle,'Extent'); + set(tHandle,'Units',oldTitleUnits); + + titleHeight=tExtent(4)+axPos(2)+axPos(4)+5; + if titleHeight>pos(4) + pos(4)=titleHeight; + pos(2)=screenSize(4)/2-pos(4)/2; + figPosDirty=logical(1); + else + figPosDirty=logical(0); + end + + if tExtent(3)>pos(3)*1.10; + pos(3)=min(tExtent(3)*1.10,screenSize(3)); + pos(1)=screenSize(3)/2-pos(3)/2; + + axPos([1,3])=axNorm([1,3])*pos(3); + set(h,'Position',axPos); + + figPosDirty=logical(1); + end + + if figPosDirty + set(f,'Position',pos); + end + + xpatch = [0 x x 0]; + ypatch = [0 0 1 1]; + xline = [100 0 0 100 100]; + yline = [0 0 1 1 0]; + + p = patch(xpatch,ypatch,'r','EdgeColor','r','EraseMode','none'); + l = line(xline,yline,'EraseMode','none'); + set(l,'Color',get(gca,'XColor')); + + + set(f,'HandleVisibility','callback','visible','on', 'resize','off'); + + set(0, 'Units', oldRootUnits); +end % case +drawnow; + +if nargout==1, + fout = f; +end |