diff options
| author | Arun Raghavan <arun.raghavan@collabora.co.uk> | 2011-10-20 13:10:08 +0530 |
|---|---|---|
| committer | Arun Raghavan <arun.raghavan@collabora.co.uk> | 2011-10-20 13:23:11 +0530 |
| commit | 7e71fffb59132109ef1eff43c5c52386df6fc759 (patch) | |
| tree | e5c59b423e8c5cba8d90dddc8080e984b77455c4 /src/modules/audio_processing/aecm | |
| parent | 139f0b6dc06618516a0658fa24bb2adbef8e0369 (diff) | |
Update code to upstream revision r766
Removes matlab tests, adds delay estimation logging, and some other
minor fixes/improvements.
Diffstat (limited to 'src/modules/audio_processing/aecm')
21 files changed, 35 insertions, 2529 deletions
diff --git a/src/modules/audio_processing/aecm/Makefile.am b/src/modules/audio_processing/aecm/Makefile.am index c36bb40..b25fd09 100644 --- a/src/modules/audio_processing/aecm/Makefile.am +++ b/src/modules/audio_processing/aecm/Makefile.am @@ -3,9 +3,7 @@ noinst_LTLIBRARIES = libaecm.la libaecm_la_SOURCES = main/interface/echo_control_mobile.h \ main/source/echo_control_mobile.c \ main/source/aecm_core.c \ - main/source/aecm_core.h \ - main/source/aecm_delay_estimator.c \ - main/source/aecm_delay_estimator.h + main/source/aecm_core.h libaecm_la_CFLAGS = $(AM_CFLAGS) $(COMMON_CFLAGS) \ -I$(top_srcdir)/src/common_audio/signal_processing_library/main/interface \ -I$(top_srcdir)/src/modules/audio_processing/utility diff --git a/src/modules/audio_processing/aecm/main/matlab/compsup.m b/src/modules/audio_processing/aecm/main/matlab/compsup.m deleted file mode 100644 index 9575ec4..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/compsup.m +++ /dev/null @@ -1,447 +0,0 @@ -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 deleted file mode 100644 index a4a533d..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/getBspectrum.m +++ /dev/null @@ -1,22 +0,0 @@ -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 deleted file mode 100644 index 5a414f9..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/hisser2.m +++ /dev/null @@ -1,21 +0,0 @@ -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 deleted file mode 100644 index 7e24c69..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/main2.m +++ /dev/null @@ -1,19 +0,0 @@ - -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 deleted file mode 100644 index 2d3e686..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/AECMobile.m +++ /dev/null @@ -1,269 +0,0 @@ -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 deleted file mode 100644 index 9b9c0ba..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/align.m +++ /dev/null @@ -1,98 +0,0 @@ -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 deleted file mode 100644 index a09a7f2..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/calcFilterGain.m +++ /dev/null @@ -1,88 +0,0 @@ -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 deleted file mode 100644 index ae1365f..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/calcStepSize.m +++ /dev/null @@ -1,105 +0,0 @@ -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 deleted file mode 100644 index d038b51..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/fallerEstimator.m +++ /dev/null @@ -1,42 +0,0 @@ -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 deleted file mode 100644 index a4a533d..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/getBspectrum.m +++ /dev/null @@ -1,22 +0,0 @@ -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 deleted file mode 100644 index 5a414f9..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/hisser2.m +++ /dev/null @@ -1,21 +0,0 @@ -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 deleted file mode 100644 index eeb2aaa..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/mainProgram.m +++ /dev/null @@ -1,283 +0,0 @@ -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 deleted file mode 100644 index 3ebe701..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/simEnvironment.m +++ /dev/null @@ -1,15 +0,0 @@ -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 deleted file mode 100644 index c805f1d..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/matlab/updateSettings.m +++ /dev/null @@ -1,94 +0,0 @@ -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 deleted file mode 100644 index 50b9ccf..0000000 --- a/src/modules/audio_processing/aecm/main/matlab/waitbar_j.m +++ /dev/null @@ -1,234 +0,0 @@ -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 diff --git a/src/modules/audio_processing/aecm/main/source/aecm.gypi b/src/modules/audio_processing/aecm/main/source/aecm.gypi index 3c63c52..bbfb1ca 100644 --- a/src/modules/audio_processing/aecm/main/source/aecm.gypi +++ b/src/modules/audio_processing/aecm/main/source/aecm.gypi @@ -28,8 +28,6 @@ 'echo_control_mobile.c', 'aecm_core.c', 'aecm_core.h', - 'aecm_delay_estimator.c', - 'aecm_delay_estimator.h', ], }, ], diff --git a/src/modules/audio_processing/aecm/main/source/aecm_core.c b/src/modules/audio_processing/aecm/main/source/aecm_core.c index b7dae90..13bffae 100644 --- a/src/modules/audio_processing/aecm/main/source/aecm_core.c +++ b/src/modules/audio_processing/aecm/main/source/aecm_core.c @@ -13,8 +13,8 @@ #include <assert.h> #include <stdlib.h> -#include "aecm_delay_estimator.h" #include "echo_control_mobile.h" +#include "delay_estimator.h" #include "ring_buffer.h" #include "typedefs.h" @@ -153,11 +153,13 @@ int WebRtcAecm_CreateCore(AecmCore_t **aecmInst) return -1; } - if (WebRtcAecm_CreateDelayEstimator(&aecm->delay_estimator, PART_LEN1, MAX_DELAY) == -1) - { - WebRtcAecm_FreeCore(aecm); - aecm = NULL; - return -1; + if (WebRtc_CreateDelayEstimator(&aecm->delay_estimator, + PART_LEN1, + MAX_DELAY, + 1) == -1) { + WebRtcAecm_FreeCore(aecm); + aecm = NULL; + return -1; } // Init some aecm pointers. 16 and 32 byte alignment is only necessary @@ -242,9 +244,8 @@ int WebRtcAecm_InitCore(AecmCore_t * const aecm, int samplingFreq) aecm->seed = 666; aecm->totCount = 0; - if (WebRtcAecm_InitDelayEstimator(aecm->delay_estimator) != 0) - { - return -1; + if (WebRtc_InitDelayEstimator(aecm->delay_estimator) != 0) { + return -1; } // Initialize to reasonable values @@ -339,7 +340,7 @@ int WebRtcAecm_FreeCore(AecmCore_t *aecm) WebRtcApm_FreeBuffer(aecm->nearCleanFrameBuf); WebRtcApm_FreeBuffer(aecm->outFrameBuf); - WebRtcAecm_FreeDelayEstimator(aecm->delay_estimator); + WebRtc_FreeDelayEstimator(aecm->delay_estimator); free(aecm); return 0; @@ -1161,6 +1162,7 @@ int WebRtcAecm_ProcessBlock(AecmCore_t * aecm, WebRtc_Word16 supGain; WebRtc_Word16 zeros32, zeros16; WebRtc_Word16 zerosDBufNoisy, zerosDBufClean, zerosXBuf; + int far_q; WebRtc_Word16 resolutionDiff, qDomainDiff; const int kMinPrefBand = 4; @@ -1200,10 +1202,10 @@ int WebRtcAecm_ProcessBlock(AecmCore_t * aecm, #endif // Transform far end signal from time domain to frequency domain. - zerosXBuf = TimeToFrequencyDomain(aecm->xBuf, - dfw, - xfa, - &xfaSum); + far_q = TimeToFrequencyDomain(aecm->xBuf, + dfw, + xfa, + &xfaSum); // Transform noisy near end signal from time domain to frequency domain. zerosDBufNoisy = TimeToFrequencyDomain(aecm->dBufNoisy, @@ -1211,7 +1213,7 @@ int WebRtcAecm_ProcessBlock(AecmCore_t * aecm, dfaNoisy, &dfaNoisySum); aecm->dfaNoisyQDomainOld = aecm->dfaNoisyQDomain; - aecm->dfaNoisyQDomain = zerosDBufNoisy; + aecm->dfaNoisyQDomain = (WebRtc_Word16)zerosDBufNoisy; if (nearendClean == NULL) @@ -1228,7 +1230,7 @@ int WebRtcAecm_ProcessBlock(AecmCore_t * aecm, dfaClean, &dfaCleanSum); aecm->dfaCleanQDomainOld = aecm->dfaCleanQDomain; - aecm->dfaCleanQDomain = zerosDBufClean; + aecm->dfaCleanQDomain = (WebRtc_Word16)zerosDBufClean; } #ifdef ARM_WINM_LOG_ @@ -1243,12 +1245,12 @@ int WebRtcAecm_ProcessBlock(AecmCore_t * aecm, // Get the delay // Save far-end history and estimate delay - delay = WebRtcAecm_DelayEstimatorProcess(aecm->delay_estimator, - xfa, - dfaNoisy, - PART_LEN1, - zerosXBuf, - aecm->currentVADValue); + delay = WebRtc_DelayEstimatorProcess(aecm->delay_estimator, + xfa, + dfaNoisy, + PART_LEN1, + far_q, + aecm->currentVADValue); if (delay < 0) { return -1; @@ -1272,16 +1274,21 @@ int WebRtcAecm_ProcessBlock(AecmCore_t * aecm, QueryPerformanceCounter((LARGE_INTEGER*)&start); #endif // Get aligned far end spectrum - far_spectrum_ptr = WebRtcAecm_GetAlignedFarend(aecm->delay_estimator, - PART_LEN1, - &zerosXBuf); + far_spectrum_ptr = WebRtc_AlignedFarend(aecm->delay_estimator, + PART_LEN1, + &far_q); + zerosXBuf = (WebRtc_Word16) far_q; if (far_spectrum_ptr == NULL) { return -1; } // Calculate log(energy) and update energy threshold levels - WebRtcAecm_CalcEnergies(aecm, far_spectrum_ptr, zerosXBuf, dfaNoisySum, echoEst32); + WebRtcAecm_CalcEnergies(aecm, + far_spectrum_ptr, + zerosXBuf, + dfaNoisySum, + echoEst32); // Calculate stepsize mu = WebRtcAecm_CalcStepSize(aecm); @@ -1923,4 +1930,3 @@ void WebRtcAecm_ResetAdaptiveChannel(AecmCore_t* aecm) } #endif // !(defined(WEBRTC_ANDROID) && defined(WEBRTC_ARCH_ARM_NEON)) - diff --git a/src/modules/audio_processing/aecm/main/source/aecm_core.h b/src/modules/audio_processing/aecm/main/source/aecm_core.h index e431c71..0dfdb04 100644 --- a/src/modules/audio_processing/aecm/main/source/aecm_core.h +++ b/src/modules/audio_processing/aecm/main/source/aecm_core.h @@ -178,7 +178,7 @@ typedef struct WebRtc_Word16 farEnergyMaxMin; WebRtc_Word16 farEnergyVAD; WebRtc_Word16 farEnergyMSE; - WebRtc_Word16 currentVADValue; + int currentVADValue; WebRtc_Word16 vadUpdateCount; WebRtc_Word16 startupState; diff --git a/src/modules/audio_processing/aecm/main/source/aecm_delay_estimator.c b/src/modules/audio_processing/aecm/main/source/aecm_delay_estimator.c deleted file mode 100644 index 2ed9e03..0000000 --- a/src/modules/audio_processing/aecm/main/source/aecm_delay_estimator.c +++ /dev/null @@ -1,604 +0,0 @@ -/* - * Copyright (c) 2011 The WebRTC project authors. All Rights Reserved. - * - * Use of this source code is governed by a BSD-style license - * that can be found in the LICENSE file in the root of the source - * tree. An additional intellectual property rights grant can be found - * in the file PATENTS. All contributing project authors may - * be found in the AUTHORS file in the root of the source tree. - */ - -#include "aecm_delay_estimator.h" - -#include <assert.h> -#include <stdlib.h> - -#include "signal_processing_library.h" -#include "typedefs.h" - -typedef struct -{ - // Pointers to mean values of spectrum and bit counts - WebRtc_Word32* mean_far_spectrum; - WebRtc_Word32* mean_near_spectrum; - WebRtc_Word32* mean_bit_counts; - - // Arrays only used locally in DelayEstimatorProcess() but whose size - // is determined at run-time. - WebRtc_Word32* bit_counts; - WebRtc_Word32* far_spectrum_32; - WebRtc_Word32* near_spectrum_32; - - // Binary history variables - WebRtc_UWord32* binary_far_history; - - // Far end history variables - WebRtc_UWord16* far_history; - int far_history_position; - WebRtc_Word16* far_q_domains; - - // Delay histogram variables - WebRtc_Word16* delay_histogram; - WebRtc_Word16 vad_counter; - - // Delay memory - int last_delay; - - // Buffer size parameters - int history_size; - int spectrum_size; - -} DelayEstimator_t; - -// Only bit |kBandFirst| through bit |kBandLast| are processed -// |kBandFirst| - |kBandLast| must be < 32 -static const int kBandFirst = 12; -static const int kBandLast = 43; - -static __inline WebRtc_UWord32 SetBit(WebRtc_UWord32 in, - WebRtc_Word32 pos) -{ - WebRtc_UWord32 mask = WEBRTC_SPL_LSHIFT_W32(1, pos); - WebRtc_UWord32 out = (in | mask); - - return out; -} - -// Compares the binary vector |binary_vector| with all rows of the binary -// matrix |binary_matrix| and counts per row the number of times they have the -// same value. -// Input: -// - binary_vector : binary "vector" stored in a long -// - binary_matrix : binary "matrix" stored as a vector of long -// - matrix_size : size of binary "matrix" -// Output: -// - bit_counts : "Vector" stored as a long, containing for each -// row the number of times the matrix row and the -// input vector have the same value -// -static void BitCountComparison(const WebRtc_UWord32 binary_vector, - const WebRtc_UWord32* binary_matrix, - int matrix_size, - WebRtc_Word32* bit_counts) -{ - int n = 0; - WebRtc_UWord32 a = binary_vector; - register WebRtc_UWord32 tmp; - - // compare binary vector |binary_vector| with all rows of the binary matrix - // |binary_matrix| - for (; n < matrix_size; n++) - { - a = (binary_vector ^ binary_matrix[n]); - // Returns bit counts in tmp - tmp = a - ((a >> 1) & 033333333333) - ((a >> 2) & 011111111111); - tmp = ((tmp + (tmp >> 3)) & 030707070707); - tmp = (tmp + (tmp >> 6)); - tmp = (tmp + (tmp >> 12) + (tmp >> 24)) & 077; - - bit_counts[n] = (WebRtc_Word32)tmp; - } -} - -// Computes the binary spectrum by comparing the input |spectrum| with a -// |threshold_spectrum|. -// -// Input: -// - spectrum : Spectrum of which the binary spectrum should -// be calculated. -// - threshold_spectrum : Threshold spectrum with which the input -// spectrum is compared. -// Return: -// - out : Binary spectrum -// -static WebRtc_UWord32 GetBinarySpectrum(WebRtc_Word32* spectrum, - WebRtc_Word32* threshold_spectrum) -{ - int k = kBandFirst; - WebRtc_UWord32 out = 0; - - for (; k <= kBandLast; k++) - { - if (spectrum[k] > threshold_spectrum[k]) - { - out = SetBit(out, k - kBandFirst); - } - } - - return out; -} - -// Calculates the mean recursively. -// -// Input: -// - new_value : new additional value -// - factor : factor for smoothing -// -// Input/Output: -// - mean_value : pointer to the mean value that should be updated -// -static void MeanEstimator(const WebRtc_Word32 new_value, - int factor, - WebRtc_Word32* mean_value) -{ - WebRtc_Word32 mean_new = *mean_value; - WebRtc_Word32 diff = new_value - mean_new; - - // mean_new = mean_value + ((new_value - mean_value) >> factor); - if (diff < 0) - { - diff = -WEBRTC_SPL_RSHIFT_W32(-diff, factor); - } - else - { - diff = WEBRTC_SPL_RSHIFT_W32(diff, factor); - } - mean_new += diff; - - *mean_value = mean_new; -} - -// Moves the pointer to the next entry and inserts new far end spectrum and -// corresponding Q-domain in its buffer. -// -// Input: -// - handle : Pointer to the delay estimation instance -// - far_spectrum : Pointer to the far end spectrum -// - far_q : Q-domain of far end spectrum -// -static void UpdateFarHistory(DelayEstimator_t* self, - WebRtc_UWord16* far_spectrum, - WebRtc_Word16 far_q) -{ - // Get new buffer position - self->far_history_position++; - if (self->far_history_position >= self->history_size) - { - self->far_history_position = 0; - } - // Update Q-domain buffer - self->far_q_domains[self->far_history_position] = far_q; - // Update far end spectrum buffer - memcpy(&(self->far_history[self->far_history_position * self->spectrum_size]), - far_spectrum, - sizeof(WebRtc_UWord16) * self->spectrum_size); -} - -int WebRtcAecm_FreeDelayEstimator(void* handle) -{ - DelayEstimator_t* self = (DelayEstimator_t*)handle; - - if (self == NULL) - { - return -1; - } - - if (self->mean_far_spectrum != NULL) - { - free(self->mean_far_spectrum); - self->mean_far_spectrum = NULL; - } - if (self->mean_near_spectrum != NULL) - { - free(self->mean_near_spectrum); - self->mean_near_spectrum = NULL; - } - if (self->mean_bit_counts != NULL) - { - free(self->mean_bit_counts); - self->mean_bit_counts = NULL; - } - if (self->bit_counts != NULL) - { - free(self->bit_counts); - self->bit_counts = NULL; - } - if (self->far_spectrum_32 != NULL) - { - free(self->far_spectrum_32); - self->far_spectrum_32 = NULL; - } - if (self->near_spectrum_32 != NULL) - { - free(self->near_spectrum_32); - self->near_spectrum_32 = NULL; - } - if (self->far_history != NULL) - { - free(self->far_history); - self->far_history = NULL; - } - if (self->binary_far_history != NULL) - { - free(self->binary_far_history); - self->binary_far_history = NULL; - } - if (self->far_q_domains != NULL) - { - free(self->far_q_domains); - self->far_q_domains = NULL; - } - if (self->delay_histogram != NULL) - { - free(self->delay_histogram); - self->delay_histogram = NULL; - } - - free(self); - - return 0; -} - -int WebRtcAecm_CreateDelayEstimator(void** handle, - int spectrum_size, - int history_size) -{ - DelayEstimator_t *self = NULL; - // Check if the sub band used in the delay estimation is small enough to - // fit in a Word32. - assert(kBandLast - kBandFirst < 32); - if (spectrum_size < kBandLast) - { - return -1; - } - if (history_size < 0) - { - return -1; - } - - self = malloc(sizeof(DelayEstimator_t)); - *handle = self; - if (self == NULL) - { - return -1; - } - - self->mean_far_spectrum = NULL; - self->mean_near_spectrum = NULL; - self->bit_counts = NULL; - self->far_spectrum_32 = NULL; - self->near_spectrum_32 = NULL; - self->far_history = NULL; - self->mean_bit_counts = NULL; - self->binary_far_history = NULL; - self->far_q_domains = NULL; - self->delay_histogram = NULL; - - // Allocate memory for spectrum buffers - self->mean_far_spectrum = malloc(spectrum_size * sizeof(WebRtc_Word32)); - if (self->mean_far_spectrum == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - self->mean_near_spectrum = malloc(spectrum_size * sizeof(WebRtc_Word32)); - if (self->mean_near_spectrum == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - self->mean_bit_counts = malloc(history_size * sizeof(WebRtc_Word32)); - if (self->mean_bit_counts == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - self->bit_counts = malloc(history_size * sizeof(WebRtc_Word32)); - if (self->bit_counts == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - self->far_spectrum_32 = malloc(spectrum_size * sizeof(WebRtc_Word32)); - if (self->far_spectrum_32 == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - self->near_spectrum_32 = malloc(spectrum_size * sizeof(WebRtc_Word32)); - if (self->near_spectrum_32 == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - // Allocate memory for history buffers - self->far_history = malloc(spectrum_size * history_size * - sizeof(WebRtc_UWord16)); - if (self->far_history == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - self->binary_far_history = malloc(history_size * sizeof(WebRtc_UWord32)); - if (self->binary_far_history == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - self->far_q_domains = malloc(history_size * sizeof(WebRtc_Word16)); - if (self->far_q_domains == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - self->delay_histogram = malloc(history_size * sizeof(WebRtc_Word16)); - if (self->delay_histogram == NULL) - { - WebRtcAecm_FreeDelayEstimator(self); - self = NULL; - return -1; - } - - self->spectrum_size = spectrum_size; - self->history_size = history_size; - - return 0; -} - -int WebRtcAecm_InitDelayEstimator(void* handle) -{ - DelayEstimator_t* self = (DelayEstimator_t*)handle; - - if (self == NULL) - { - return -1; - } - // Set averaged far and near end spectra to zero - memset(self->mean_far_spectrum, - 0, - sizeof(WebRtc_Word32) * self->spectrum_size); - memset(self->mean_near_spectrum, - 0, - sizeof(WebRtc_Word32) * self->spectrum_size); - // Set averaged bit counts to zero - memset(self->mean_bit_counts, - 0, - sizeof(WebRtc_Word32) * self->history_size); - memset(self->bit_counts, - 0, - sizeof(WebRtc_Word32) * self->history_size); - memset(self->far_spectrum_32, - 0, - sizeof(WebRtc_Word32) * self->spectrum_size); - memset(self->near_spectrum_32, - 0, - sizeof(WebRtc_Word32) * self->spectrum_size); - // Set far end histories to zero - memset(self->binary_far_history, - 0, - sizeof(WebRtc_UWord32) * self->history_size); - memset(self->far_history, - 0, - sizeof(WebRtc_UWord16) * self->spectrum_size * - self->history_size); - memset(self->far_q_domains, - 0, - sizeof(WebRtc_Word16) * self->history_size); - - self->far_history_position = self->history_size; - // Set delay histogram to zero - memset(self->delay_histogram, - 0, - sizeof(WebRtc_Word16) * self->history_size); - // Set VAD counter to zero - self->vad_counter = 0; - // Set delay memory to zero - self->last_delay = 0; - - return 0; -} - -int WebRtcAecm_DelayEstimatorProcess(void* handle, - WebRtc_UWord16* far_spectrum, - WebRtc_UWord16* near_spectrum, - int spectrum_size, - WebRtc_Word16 far_q, - WebRtc_Word16 vad_value) -{ - DelayEstimator_t* self = (DelayEstimator_t*)handle; - - WebRtc_UWord32 bxspectrum, byspectrum; - - int i; - - WebRtc_Word32 dtmp1; - - WebRtc_Word16 maxHistLvl = 0; - WebRtc_Word16 minpos = -1; - - const int kVadCountThreshold = 25; - const int kMaxHistogram = 600; - - if (self == NULL) - { - return -1; - } - - if (spectrum_size != self->spectrum_size) - { - // Data sizes don't match - return -1; - } - if (far_q > 15) - { - // If far_Q is larger than 15 we can not guarantee no wrap around - return -1; - } - - // Update far end history - UpdateFarHistory(self, far_spectrum, far_q); - // Update the far and near end means - for (i = 0; i < self->spectrum_size; i++) - { - self->far_spectrum_32[i] = (WebRtc_Word32)far_spectrum[i]; - MeanEstimator(self->far_spectrum_32[i], 6, &(self->mean_far_spectrum[i])); - - self->near_spectrum_32[i] = (WebRtc_Word32)near_spectrum[i]; - MeanEstimator(self->near_spectrum_32[i], 6, &(self->mean_near_spectrum[i])); - } - - // Shift binary spectrum history - memmove(&(self->binary_far_history[1]), - &(self->binary_far_history[0]), - (self->history_size - 1) * sizeof(WebRtc_UWord32)); - - // Get binary spectra - bxspectrum = GetBinarySpectrum(self->far_spectrum_32, self->mean_far_spectrum); - byspectrum = GetBinarySpectrum(self->near_spectrum_32, self->mean_near_spectrum); - // Insert new binary spectrum - self->binary_far_history[0] = bxspectrum; - - // Compare with delayed spectra - BitCountComparison(byspectrum, - self->binary_far_history, - self->history_size, - self->bit_counts); - - // Smooth bit count curve - for (i = 0; i < self->history_size; i++) - { - // Update sum - // |bit_counts| is constrained to [0, 32], meaning we can smooth with a - // factor up to 2^26. We use Q9. - dtmp1 = WEBRTC_SPL_LSHIFT_W32(self->bit_counts[i], 9); // Q9 - MeanEstimator(dtmp1, 9, &(self->mean_bit_counts[i])); - } - - // Find minimum position of bit count curve - minpos = WebRtcSpl_MinIndexW32(self->mean_bit_counts, self->history_size); - - // If the farend has been active sufficiently long, begin accumulating a - // histogram of the minimum positions. Search for the maximum bin to - // determine the delay. - if (vad_value == 1) - { - if (self->vad_counter >= kVadCountThreshold) - { - // Increment the histogram at the current minimum position. - if (self->delay_histogram[minpos] < kMaxHistogram) - { - self->delay_histogram[minpos] += 3; - } - -#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) - // Decrement the entire histogram. - // Select the histogram index corresponding to the maximum bin as - // the delay. - self->last_delay = 0; - for (i = 0; i < self->history_size; i++) - { - if (self->delay_histogram[i] > 0) - { - self->delay_histogram[i]--; - } - if (self->delay_histogram[i] > maxHistLvl) - { - maxHistLvl = self->delay_histogram[i]; - self->last_delay = i; - } - } -#else - self->last_delay = 0; - - for (i = 0; i < self->history_size; i++) - { - WebRtc_Word16 tempVar = self->delay_histogram[i]; - - // Decrement the entire histogram. - if (tempVar > 0) - { - tempVar--; - self->delay_histogram[i] = tempVar; - - // Select the histogram index corresponding to the maximum - // bin as the delay. - if (tempVar > maxHistLvl) - { - maxHistLvl = tempVar; - self->last_delay = i; - } - } - } -#endif - } else - { - self->vad_counter++; - } - } else - { - self->vad_counter = 0; - } - - return self->last_delay; -} - -const WebRtc_UWord16* WebRtcAecm_GetAlignedFarend(void* handle, - int far_spectrum_size, - WebRtc_Word16* far_q) -{ - DelayEstimator_t* self = (DelayEstimator_t*)handle; - int buffer_position = 0; - - if (self == NULL) - { - return NULL; - } - if (far_spectrum_size != self->spectrum_size) - { - return NULL; - } - - // Get buffer position - buffer_position = self->far_history_position - self->last_delay; - if (buffer_position < 0) - { - buffer_position += self->history_size; - } - // Get Q-domain - *far_q = self->far_q_domains[buffer_position]; - // Return far end spectrum - return (self->far_history + (buffer_position * self->spectrum_size)); - -} - -int WebRtcAecm_GetLastDelay(void* handle) -{ - DelayEstimator_t* self = (DelayEstimator_t*)handle; - - if (self == NULL) - { - return -1; - } - - // Return last calculated delay - return self->last_delay; -} diff --git a/src/modules/audio_processing/aecm/main/source/aecm_delay_estimator.h b/src/modules/audio_processing/aecm/main/source/aecm_delay_estimator.h deleted file mode 100644 index 5ce57fa..0000000 --- a/src/modules/audio_processing/aecm/main/source/aecm_delay_estimator.h +++ /dev/null @@ -1,112 +0,0 @@ -/* - * Copyright (c) 2011 The WebRTC project authors. All Rights Reserved. - * - * Use of this source code is governed by a BSD-style license - * that can be found in the LICENSE file in the root of the source - * tree. An additional intellectual property rights grant can be found - * in the file PATENTS. All contributing project authors may - * be found in the AUTHORS file in the root of the source tree. - */ - -// Performs delay estimation on a block by block basis -// The return value is 0 - OK and -1 - Error, unless otherwise stated. - -#ifndef WEBRTC_MODULES_AUDIO_PROCESSING_AECM_MAIN_SOURCE_AECM_DELAY_ESTIMATOR_H_ -#define WEBRTC_MODULES_AUDIO_PROCESSING_AECM_MAIN_SOURCE_AECM_DELAY_ESTIMATOR_H_ - -#include "typedefs.h" - -// Releases the memory allocated by WebRtcAecm_CreateDelayEstimator(...) -// Input: -// - handle : Pointer to the delay estimation instance -// -int WebRtcAecm_FreeDelayEstimator(void* handle); - -// Allocates the memory needed by the delay estimation. The memory needs to be -// initialized separately using the WebRtcAecm_InitDelayEstimator(...) function. -// -// Input: -// - handle : Instance that should be created -// - spectrum_size : Size of the spectrum used both in far end and near -// end. Used to allocate memory for spectrum specific -// buffers. -// - history_size : Size of the far end history used to estimate the -// delay from. Used to allocate memory for history -// specific buffers. -// -// Output: -// - handle : Created instance -// -int WebRtcAecm_CreateDelayEstimator(void** handle, - int spectrum_size, - int history_size); - -// Initializes the delay estimation instance created with -// WebRtcAecm_CreateDelayEstimator(...) -// Input: -// - handle : Pointer to the delay estimation instance -// -// Output: -// - handle : Initialized instance -// -int WebRtcAecm_InitDelayEstimator(void* handle); - -// Estimates and returns the delay between the far end and near end blocks. -// Input: -// - handle : Pointer to the delay estimation instance -// - far_spectrum : Pointer to the far end spectrum data -// - near_spectrum : Pointer to the near end spectrum data of the current -// block -// - spectrum_size : The size of the data arrays (same for both far and -// near end) -// - far_q : The Q-domain of the far end data -// - vad_value : The VAD decision of the current block -// -// Output: -// - handle : Updated instance -// -// Return value: -// - delay : >= 0 - Calculated delay value -// -1 - Error -// -int WebRtcAecm_DelayEstimatorProcess(void* handle, - WebRtc_UWord16* far_spectrum, - WebRtc_UWord16* near_spectrum, - int spectrum_size, - WebRtc_Word16 far_q, - WebRtc_Word16 vad_value); - -// Returns a pointer to the far end spectrum aligned to current near end -// spectrum. The function WebRtcAecm_DelayEstimatorProcess(...) should -// have been called before WebRtcAecm_GetAlignedFarend(...). Otherwise, you get -// the pointer to the previous frame. The memory is only valid until the next -// call of WebRtcAecm_DelayEstimatorProcess(...). -// -// Inputs: -// - handle : Pointer to the delay estimation instance -// -// Output: -// - far_spectrum_size : Size of far_spectrum allocated by the caller -// - far_q : The Q-domain of the aligned far end spectrum -// -// Return value: -// - far_spectrum : Pointer to the aligned far end spectrum -// NULL - Error -// -const WebRtc_UWord16* WebRtcAecm_GetAlignedFarend(void* handle, - int far_spectrum_size, - WebRtc_Word16* far_q); - -// Returns the last calculated delay updated by the function -// WebRtcAecm_DelayEstimatorProcess(...) -// -// Inputs: -// - handle : Pointer to the delay estimation instance -// -// Return value: -// - delay : >= 0 - Last calculated delay value -// -1 - Error -// -int WebRtcAecm_GetLastDelay(void* handle); - -#endif // WEBRTC_MODULES_AUDIO_PROCESSING_AECM_MAIN_SOURCE_AECM_DELAY_ESTIMATOR_H_ |