diff options
Diffstat (limited to 'src/modules/audio_processing/aecm/main')
21 files changed, 5689 insertions, 0 deletions
diff --git a/src/modules/audio_processing/aecm/main/interface/echo_control_mobile.h b/src/modules/audio_processing/aecm/main/interface/echo_control_mobile.h new file mode 100644 index 0000000000..26b1172726 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/interface/echo_control_mobile.h @@ -0,0 +1,206 @@ +/* + * 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. + */ + +#ifndef WEBRTC_MODULES_AUDIO_PROCESSING_AECM_MAIN_INTERFACE_ECHO_CONTROL_MOBILE_H_ +#define WEBRTC_MODULES_AUDIO_PROCESSING_AECM_MAIN_INTERFACE_ECHO_CONTROL_MOBILE_H_ + +#include "typedefs.h" + +enum { + AecmFalse = 0, + AecmTrue +}; + +// Errors +#define AECM_UNSPECIFIED_ERROR 12000 +#define AECM_UNSUPPORTED_FUNCTION_ERROR 12001 +#define AECM_UNINITIALIZED_ERROR 12002 +#define AECM_NULL_POINTER_ERROR 12003 +#define AECM_BAD_PARAMETER_ERROR 12004 + +// Warnings +#define AECM_BAD_PARAMETER_WARNING 12100 + +typedef struct { + WebRtc_Word16 cngMode; // AECM_FALSE, AECM_TRUE (default) + WebRtc_Word16 echoMode; // 0, 1, 2, 3 (default), 4 +} AecmConfig; + +#ifdef __cplusplus +extern "C" { +#endif + +/* + * Allocates the memory needed by the AECM. The memory needs to be + * initialized separately using the WebRtcAecm_Init() function. + * + * Inputs Description + * ------------------------------------------------------------------- + * void **aecmInst Pointer to the AECM instance to be + * created and initialized + * + * Outputs Description + * ------------------------------------------------------------------- + * WebRtc_Word32 return 0: OK + * -1: error + */ +WebRtc_Word32 WebRtcAecm_Create(void **aecmInst); + +/* + * This function releases the memory allocated by WebRtcAecm_Create() + * + * Inputs Description + * ------------------------------------------------------------------- + * void *aecmInst Pointer to the AECM instance + * + * Outputs Description + * ------------------------------------------------------------------- + * WebRtc_Word32 return 0: OK + * -1: error + */ +WebRtc_Word32 WebRtcAecm_Free(void *aecmInst); + +/* + * Initializes an AECM instance. + * + * Inputs Description + * ------------------------------------------------------------------- + * void *aecmInst Pointer to the AECM instance + * WebRtc_Word32 sampFreq Sampling frequency of data + * WebRtc_Word32 scSampFreq Soundcard sampling frequency + * + * Outputs Description + * ------------------------------------------------------------------- + * WebRtc_Word32 return 0: OK + * -1: error + */ +WebRtc_Word32 WebRtcAecm_Init(void* aecmInst, + WebRtc_Word32 sampFreq, + WebRtc_Word32 scSampFreq); + +/* + * Inserts an 80 or 160 sample block of data into the farend buffer. + * + * Inputs Description + * ------------------------------------------------------------------- + * void *aecmInst Pointer to the AECM instance + * WebRtc_Word16 *farend In buffer containing one frame of + * farend signal + * WebRtc_Word16 nrOfSamples Number of samples in farend buffer + * + * Outputs Description + * ------------------------------------------------------------------- + * WebRtc_Word32 return 0: OK + * -1: error + */ +WebRtc_Word32 WebRtcAecm_BufferFarend(void* aecmInst, + const WebRtc_Word16* farend, + WebRtc_Word16 nrOfSamples); + +/* + * Runs the AECM on an 80 or 160 sample blocks of data. + * + * Inputs Description + * ------------------------------------------------------------------- + * void *aecmInst Pointer to the AECM instance + * WebRtc_Word16 *nearendNoisy In buffer containing one frame of + * reference nearend+echo signal. If + * noise reduction is active, provide + * the noisy signal here. + * WebRtc_Word16 *nearendClean In buffer containing one frame of + * nearend+echo signal. If noise + * reduction is active, provide the + * clean signal here. Otherwise pass a + * NULL pointer. + * WebRtc_Word16 nrOfSamples Number of samples in nearend buffer + * WebRtc_Word16 msInSndCardBuf Delay estimate for sound card and + * system buffers + * + * Outputs Description + * ------------------------------------------------------------------- + * WebRtc_Word16 *out Out buffer, one frame of processed nearend + * WebRtc_Word32 return 0: OK + * -1: error + */ +WebRtc_Word32 WebRtcAecm_Process(void* aecmInst, + const WebRtc_Word16* nearendNoisy, + const WebRtc_Word16* nearendClean, + WebRtc_Word16* out, + WebRtc_Word16 nrOfSamples, + WebRtc_Word16 msInSndCardBuf); + +/* + * This function enables the user to set certain parameters on-the-fly + * + * Inputs Description + * ------------------------------------------------------------------- + * void *aecmInst Pointer to the AECM instance + * AecmConfig config Config instance that contains all + * properties to be set + * + * Outputs Description + * ------------------------------------------------------------------- + * WebRtc_Word32 return 0: OK + * -1: error + */ +WebRtc_Word32 WebRtcAecm_set_config(void* aecmInst, + AecmConfig config); + +/* + * This function enables the user to set certain parameters on-the-fly + * + * Inputs Description + * ------------------------------------------------------------------- + * void *aecmInst Pointer to the AECM instance + * + * Outputs Description + * ------------------------------------------------------------------- + * AecmConfig *config Pointer to the config instance that + * all properties will be written to + * WebRtc_Word32 return 0: OK + * -1: error + */ +WebRtc_Word32 WebRtcAecm_get_config(void *aecmInst, + AecmConfig *config); + +/* + * Gets the last error code. + * + * Inputs Description + * ------------------------------------------------------------------- + * void *aecmInst Pointer to the AECM instance + * + * Outputs Description + * ------------------------------------------------------------------- + * WebRtc_Word32 return 11000-11100: error code + */ +WebRtc_Word32 WebRtcAecm_get_error_code(void *aecmInst); + +/* + * Gets a version string + * + * Inputs Description + * ------------------------------------------------------------------- + * char *versionStr Pointer to a string array + * WebRtc_Word16 len The maximum length of the string + * + * Outputs Description + * ------------------------------------------------------------------- + * WebRtc_Word8 *versionStr Pointer to a string array + * WebRtc_Word32 return 0: OK + * -1: error + */ +WebRtc_Word32 WebRtcAecm_get_version(WebRtc_Word8 *versionStr, + WebRtc_Word16 len); + +#ifdef __cplusplus +} +#endif +#endif /* WEBRTC_MODULES_AUDIO_PROCESSING_AECM_MAIN_INTERFACE_ECHO_CONTROL_MOBILE_H_ */ diff --git a/src/modules/audio_processing/aecm/main/matlab/compsup.m b/src/modules/audio_processing/aecm/main/matlab/compsup.m new file mode 100644 index 0000000000..9575ec40fc --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/compsup.m @@ -0,0 +1,447 @@ +function [emicrophone,aaa]=compsup(microphone,TheFarEnd,avtime,samplingfreq); +% microphone = microphone signal +% aaa = nonlinearity input variable +% TheFarEnd = far end signal +% avtime = interval to compute suppression from (seconds) +% samplingfreq = sampling frequency + +%if(nargin==6) +% fprintf(1,'suppress has received a delay sequence\n'); +%end + + +Ap500=[ 1.00, -4.95, 9.801, -9.70299, 4.80298005, -0.9509900499]; +Bp500=[ 0.662743088639636, -2.5841655608125, 3.77668102146288, -2.45182477425154, 0.596566274575251, 0.0]; + + +Ap200=[ 1.00, -4.875, 9.50625, -9.26859375, 4.518439453125, -0.881095693359375]; +Bp200=[ 0.862545460994275, -3.2832804496114, 4.67892032308828, -2.95798023879133, 0.699796870041299, 0.0]; + +maxDelay=0.4; %[s] +histLen=1; %[s] + + +% CONSTANTS THAT YOU CAN EXPERIMENT WITH +A_GAIN=10.0; % for the suppress case +oversampling = 2; % must be power of 2; minimum is 2; 4 works +% fine for support=64, but for support=128, +% 8 gives better results. +support=64; %512 % fft support (frequency resolution; at low +% settings you can hear more distortion +% (e.g. pitch that is left-over from far-end)) +% 128 works well, 64 is ok) + +lowlevel = mean(abs(microphone))*0.0001; + +G_ol = 0; % Use overlapping sets of estimates + +% ECHO SUPPRESSION SPECIFIC PARAMETERS +suppress_overdrive=1.0; % overdrive factor for suppression 1.4 is good +gamma_echo=1.0; % same as suppress_overdrive but at different place +de_echo_bound=0.0; +mLim=10; % rank of matrix G +%limBW = 1; % use bandwidth-limited response for G +if mLim > (support/2+1) + error('mLim in suppress.m too large\n'); +end + + +dynrange=1.0000e-004; + +% other, constants +hsupport = support/2; +hsupport1 = hsupport+1; +factor = 2 / oversampling; +updatel = support/oversampling; +win=sqrt(designwindow(0,support)); +estLen = round(avtime * samplingfreq/updatel) + +runningfmean =0.0; + +mLim = floor(hsupport1/2); +V = sqrt(2/hsupport1)*cos(pi/hsupport1*(repmat((0:hsupport1-1) + 0.5, mLim, 1).* ... + repmat((0:mLim-1)' + 0.5, 1, hsupport1))); + +fprintf(1,'updatel is %5.3f s\n', updatel/samplingfreq); + + + +bandfirst=8; bandlast=25; +dosmooth=0; % to get rid of wavy bin counts (can be worse or better) + +% compute some constants +blockLen = support/oversampling; +maxDelayb = floor(samplingfreq*maxDelay/updatel); % in blocks +histLenb = floor(samplingfreq*histLen/updatel); % in blocks + +x0=TheFarEnd; +y0=microphone; + + +%input +tlength=min([length(microphone),length(TheFarEnd)]); +updateno=floor(tlength/updatel); +tlength=updatel*updateno; +updateno = updateno - oversampling + 1; + +TheFarEnd =TheFarEnd(1:tlength); +microphone =microphone(1:tlength); + +TheFarEnd =[zeros(hsupport,1);TheFarEnd(1:tlength)]; +microphone =[zeros(hsupport,1);microphone(1:tlength)]; + + +% signal length +n = min([floor(length(x0)/support)*support,floor(length(y0)/support)*support]); +nb = n/blockLen - oversampling + 1; % in blocks + +% initialize space +win = sqrt([0 ; hanning(support-1)]); +sxAll2 = zeros(hsupport1,nb); +syAll2 = zeros(hsupport1,nb); + +z500=zeros(5,maxDelayb+1); +z200=zeros(5,hsupport1); + +bxspectrum=uint32(zeros(nb,1)); +bxhist=uint32(zeros(maxDelayb+1,1)); +byspectrum=uint32(zeros(nb,1)); +bcount=zeros(1+maxDelayb,nb); +fcount=zeros(1+maxDelayb,nb); +fout=zeros(1+maxDelayb,nb); +delay=zeros(nb,1); +tdelay=zeros(nb,1); +nlgains=zeros(nb,1); + +% create space (mainly for debugging) +emicrophone=zeros(tlength,1); +femicrophone=complex(zeros(hsupport1,updateno)); +thefilter=zeros(hsupport1,updateno); +thelimiter=ones(hsupport1,updateno); +fTheFarEnd=complex(zeros(hsupport1,updateno)); +afTheFarEnd=zeros(hsupport1,updateno); +fmicrophone=complex(zeros(hsupport1,updateno)); +afmicrophone=zeros(hsupport1,updateno); + +G = zeros(hsupport1, hsupport1); +zerovec = zeros(hsupport1,1); +zeromat = zeros(hsupport1); + +% Reset sums +mmxs_a = zerovec; +mmys_a = zerovec; +s2xs_a = zerovec; +s2ys_a = zerovec; +Rxxs_a = zeromat; +Ryxs_a = zeromat; +count_a = 1; + +mmxs_b = zerovec; +mmys_b = zerovec; +s2xs_b = zerovec; +s2ys_b = zerovec; +Rxxs_b = zeromat; +Ryxs_b = zeromat; +count_b = 1; + +nog=0; + +aaa=zeros(size(TheFarEnd)); + +% loop over signal blocks +fprintf(1,'.. Suppression; averaging G over %5.1f seconds; file length %5.1f seconds ..\n',avtime, length(microphone)/samplingfreq); +fprintf(1,'.. SUPPRESSING ONLY AFTER %5.1f SECONDS! ..\n',avtime); +fprintf(1,'.. 20 seconds is good ..\n'); +hh = waitbar_j(0,'Please wait...'); + + +for i=1:updateno + + sb = (i-1)*updatel + 1; + se=sb+support-1; + + % analysis FFTs + temp=fft(win .* TheFarEnd(sb:se)); + fTheFarEnd(:,i)=temp(1:hsupport1); + xf=fTheFarEnd(:,i); + afTheFarEnd(:,i)= abs(fTheFarEnd(:,i)); + + temp=win .* microphone(sb:se); + + temp=fft(win .* microphone(sb:se)); + fmicrophone(:,i)=temp(1:hsupport1); + yf=fmicrophone(:,i); + afmicrophone(:,i)= abs(fmicrophone(:,i)); + + + ener_orig = afmicrophone(:,i)'*afmicrophone(:,i); + if( ener_orig == 0) + afmicrophone(:,i)=lowlevel*ones(size(afmicrophone(:,i))); + end + + + % use log domain (showed improved performance) +xxf= sqrt(real(xf.*conj(xf))+1e-20); +yyf= sqrt(real(yf.*conj(yf))+1e-20); + sxAll2(:,i) = 20*log10(xxf); + syAll2(:,i) = 20*log10(yyf); + + mD=min(i-1,maxDelayb); + xthreshold = sum(sxAll2(:,i-mD:i),2)/(maxDelayb+1); + + [yout, z200] = filter(Bp200,Ap200,syAll2(:,i),z200,2); + yout=yout/(maxDelayb+1); + ythreshold = mean(syAll2(:,i-mD:i),2); + + + bxspectrum(i)=getBspectrum(sxAll2(:,i),xthreshold,bandfirst,bandlast); + byspectrum(i)=getBspectrum(syAll2(:,i),yout,bandfirst,bandlast); + + bxhist(end-mD:end)=bxspectrum(i-mD:i); + + bcount(:,i)=hisser2( ... + byspectrum(i),flipud(bxhist),bandfirst,bandlast); + + + [fout(:,i), z500] = filter(Bp500,Ap500,bcount(:,i),z500,2); + fcount(:,i)=sum(bcount(:,max(1,i-histLenb+1):i),2); % using the history range + fout(:,i)=round(fout(:,i)); + [value,delay(i)]=min(fout(:,i),[],1); + tdelay(i)=(delay(i)-1)*support/(samplingfreq*oversampling); + + % compensate + + idel = max(i - delay(i) + 1,1); + + + % echo suppression + + noisyspec = afmicrophone(:,i); + + % Estimate G using covariance matrices + + % Cumulative estimates + xx = afTheFarEnd(:,idel); + yy = afmicrophone(:,i); + + % Means + mmxs_a = mmxs_a + xx; + mmys_a = mmys_a + yy; + if (G_ol) + mmxs_b = mmxs_b + xx; + mmys_b = mmys_b + yy; + mmy = mean([mmys_a/count_a mmys_b/count_b],2); + mmx = mean([mmxs_a/count_a mmxs_b/count_b],2); + else + mmx = mmxs_a/count_a; + mmy = mmys_a/count_a; + end + count_a = count_a + 1; + count_b = count_b + 1; + + % Mean removal + xxm = xx - mmx; + yym = yy - mmy; + + % Variances + s2xs_a = s2xs_a + xxm .* xxm; + s2ys_a = s2ys_a + yym .* yym; + s2xs_b = s2xs_b + xxm .* xxm; + s2ys_b = s2ys_b + yym .* yym; + + % Correlation matrices + Rxxs_a = Rxxs_a + xxm * xxm'; + Ryxs_a = Ryxs_a + yym * xxm'; + Rxxs_b = Rxxs_b + xxm * xxm'; + Ryxs_b = Ryxs_b + yym * xxm'; + + + % Gain matrix A + + if mod(i, estLen) == 0 + + + % Cumulative based estimates + Rxxf = Rxxs_a / (estLen - 1); + Ryxf = Ryxs_a / (estLen - 1); + + % Variance normalization + s2x2 = s2xs_a / (estLen - 1); + s2x2 = sqrt(s2x2); + % Sx = diag(max(s2x2,dynrange*max(s2x2))); + Sx = diag(s2x2); + if (sum(s2x2) > 0) + iSx = inv(Sx); + else + iSx= Sx + 0.01; + end + + s2y2 = s2ys_a / (estLen - 1); + s2y2 = sqrt(s2y2); + % Sy = diag(max(s2y2,dynrange*max(s2y2))); + Sy = diag(s2y2); + iSy = inv(Sy); + rx = iSx * Rxxf * iSx; + ryx = iSy * Ryxf * iSx; + + + + dbd= 7; % Us less than the full matrix + + % k x m + % Bandlimited structure on G + LSEon = 0; % Default is using MMSE + if (LSEon) + ryx = ryx*rx; + rx = rx*rx; + end + p = dbd-1; + gaj = min(min(hsupport1,2*p+1),min([p+(1:hsupport1); hsupport1+p+1-(1:hsupport1)])); + cgaj = [0 cumsum(gaj)]; + + G3 = zeros(hsupport1); + for kk=1:hsupport1 + ki = max(0,kk-p-1); + if (sum(sum(rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk))))>0) + G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk))/rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk)); + else + G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk)); + end + end + % End Bandlimited structure + + G = G3; + G(abs(G)<0.01)=0; + G = suppress_overdrive * Sy * G * iSx; + + if 1 + figure(32); mi=2; + surf(max(min(G,mi),-mi)); view(2) + title('Unscaled Masked Limited-bandwidth G'); + end + pause(0.05); + + % Reset sums + mmxs_a = zerovec; + mmys_a = zerovec; + s2xs_a = zerovec; + s2ys_a = zerovec; + Rxxs_a = zeromat; + Ryxs_a = zeromat; + count_a = 1; + + end + + if (G_ol) + % Gain matrix B + + if ((mod((i-estLen/2), estLen) == 0) & i>estLen) + + + % Cumulative based estimates + Rxxf = Rxxs_b / (estLen - 1); + Ryxf = Ryxs_b / (estLen - 1); + + % Variance normalization + s2x2 = s2xs_b / (estLen - 1); + s2x2 = sqrt(s2x2); + Sx = diag(max(s2x2,dynrange*max(s2x2))); + iSx = inv(Sx); + s2y2 = s2ys_b / (estLen - 1); + s2y2 = sqrt(s2y2); + Sy = diag(max(s2y2,dynrange*max(s2y2))); + iSy = inv(Sy); + rx = iSx * Rxxf * iSx; + ryx = iSy * Ryxf * iSx; + + + % Bandlimited structure on G + LSEon = 0; % Default is using MMSE + if (LSEon) + ryx = ryx*rx; + rx = rx*rx; + end + p = dbd-1; + gaj = min(min(hsupport1,2*p+1),min([p+(1:hsupport1); hsupport1+p+1-(1:hsupport1)])); + cgaj = [0 cumsum(gaj)]; + + G3 = zeros(hsupport1); + for kk=1:hsupport1 + ki = max(0,kk-p-1); + G3(kk,ki+1:ki+gaj(kk)) = ryx(kk,ki+1:ki+gaj(kk))/rx(ki+1:ki+gaj(kk),ki+1:ki+gaj(kk)); + end + % End Bandlimited structure + + G = G3; + G(abs(G)<0.01)=0; + G = suppress_overdrive * Sy * G * iSx; + + if 1 + figure(32); mi=2; + surf(max(min(G,mi),-mi)); view(2) + title('Unscaled Masked Limited-bandwidth G'); + end + pause(0.05); + + + % Reset sums + mmxs_b = zerovec; + mmys_b = zerovec; + s2xs_b = zerovec; + s2ys_b = zerovec; + Rxxs_b = zeromat; + Ryxs_b = zeromat; + count_b = 1; + + end + + end + + FECestimate2 = G*afTheFarEnd(:,idel); + + % compute Wiener filter and suppressor function + thefilter(:,i) = (noisyspec - gamma_echo*FECestimate2) ./ noisyspec; + ix0 = find(thefilter(:,i)<de_echo_bound); % bounding trick 1 + thefilter(ix0,i) = de_echo_bound; % bounding trick 2 + ix0 = find(thefilter(:,i)>1); % bounding in reasonable range + thefilter(ix0,i) = 1; + + % NONLINEARITY + nl_alpha=0.8; % memory; seems not very critical + nlSeverity=0.3; % nonlinearity severity: 0 does nothing; 1 suppresses all + thefmean=mean(thefilter(8:16,i)); + if (thefmean<1) + disp(''); + end + runningfmean = nl_alpha*runningfmean + (1-nl_alpha)*thefmean; + aaa(sb+20+1:sb+20+updatel)=10000*runningfmean* ones(updatel,1); % debug + slope0=1.0/(1.0-nlSeverity); % + thegain = max(0.0,min(1.0,slope0*(runningfmean-nlSeverity))); + % END NONLINEARITY + thefilter(:,i) = thegain*thefilter(:,i); + + + % Wiener filtering + femicrophone(:,i) = fmicrophone(:,i) .* thefilter(:,i); + thelimiter(:,i) = (noisyspec - A_GAIN*FECestimate2) ./ noisyspec; + index = find(thelimiter(:,i)>1.0); + thelimiter(index,i) = 1.0; + index = find(thelimiter(:,i)<0.0); + thelimiter(index,i) = 0.0; + + if (rem(i,floor(updateno/20))==0) + fprintf(1,'.'); + end + if mod(i,50)==0 + waitbar_j(i/updateno,hh); + end + + + % reconstruction; first make spectrum odd + temp=[femicrophone(:,i);flipud(conj(femicrophone(2:hsupport,i)))]; + emicrophone(sb:se) = emicrophone(sb:se) + factor * win .* real(ifft(temp)); + +end +fprintf(1,'\n'); + +close(hh);
\ No newline at end of file diff --git a/src/modules/audio_processing/aecm/main/matlab/getBspectrum.m b/src/modules/audio_processing/aecm/main/matlab/getBspectrum.m new file mode 100644 index 0000000000..a4a533d600 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/getBspectrum.m @@ -0,0 +1,22 @@ +function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast) +% function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast) +% compute binary spectrum using threshold spectrum as pivot +% bspectrum = binary spectrum (binary) +% ps=current power spectrum (float) +% threshold=threshold spectrum (float) +% bandfirst = first band considered +% bandlast = last band considered + +% initialization stuff + if( length(ps)<bandlast | bandlast>32 | length(ps)~=length(threshold)) + error('BinDelayEst:spectrum:invalid','Dimensionality error'); +end + +% get current binary spectrum +diff = ps - threshold; +bspectrum=uint32(0); +for(i=bandfirst:bandlast) + if( diff(i)>0 ) + bspectrum = bitset(bspectrum,i); + end +end diff --git a/src/modules/audio_processing/aecm/main/matlab/hisser2.m b/src/modules/audio_processing/aecm/main/matlab/hisser2.m new file mode 100644 index 0000000000..5a414f9da8 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/hisser2.m @@ -0,0 +1,21 @@ +function bcount=hisser2(bs,bsr,bandfirst,bandlast) +% function bcount=hisser(bspectrum,bandfirst,bandlast) +% histogram for the binary spectra +% bcount= array of bit counts +% bs=binary spectrum (one int32 number each) +% bsr=reference binary spectra (one int32 number each) +% blockSize = histogram over blocksize blocks +% bandfirst = first band considered +% bandlast = last band considered + +% weight all delays equally +maxDelay = length(bsr); + +% compute counts (two methods; the first works better and is operational) +bcount=zeros(maxDelay,1); +for(i=1:maxDelay) + % the delay should have low count for low-near&high-far and high-near&low-far + bcount(i)= sum(bitget(bitxor(bs,bsr(i)),bandfirst:bandlast)); + % the delay should have low count for low-near&high-far (works less well) +% bcount(i)= sum(bitget(bitand(bsr(i),bitxor(bs,bsr(i))),bandfirst:bandlast)); +end diff --git a/src/modules/audio_processing/aecm/main/matlab/main2.m b/src/modules/audio_processing/aecm/main/matlab/main2.m new file mode 100644 index 0000000000..7e24c69ccf --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/main2.m @@ -0,0 +1,19 @@ + +fid=fopen('aecfar.pcm'); far=fread(fid,'short'); fclose(fid); +fid=fopen('aecnear.pcm'); mic=fread(fid,'short'); fclose(fid); + +%fid=fopen('QA1far.pcm'); far=fread(fid,'short'); fclose(fid); +%fid=fopen('QA1near.pcm'); mic=fread(fid,'short'); fclose(fid); + +start=0 * 8000+1; +stop= 30 * 8000; +microphone=mic(start:stop); +TheFarEnd=far(start:stop); +avtime=1; + +% 16000 to make it compatible with the C-version +[emicrophone,tdel]=compsup(microphone,TheFarEnd,avtime,16000); + +spclab(8000,TheFarEnd,microphone,emicrophone); + + diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/AECMobile.m b/src/modules/audio_processing/aecm/main/matlab/matlab/AECMobile.m new file mode 100644 index 0000000000..2d3e6867df --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/AECMobile.m @@ -0,0 +1,269 @@ +function [femicrophone, aecmStructNew, enerNear, enerFar] = AECMobile(fmicrophone, afTheFarEnd, setupStruct, aecmStruct) +global NEARENDFFT; +global F; + +aecmStructNew = aecmStruct; + +% Magnitude spectrum of near end signal +afmicrophone = abs(fmicrophone); +%afmicrophone = NEARENDFFT(setupStruct.currentBlock,:)'/2^F(setupStruct.currentBlock,end); +% Near end energy level +ener_orig = afmicrophone'*afmicrophone; +if( ener_orig == 0) + lowlevel = 0.01; + afmicrophone = lowlevel*ones(size(afmicrophone)); +end +%adiff = max(abs(afmicrophone - afTheFarEnd)); +%if (adiff > 0) +% disp([setupStruct.currentBlock adiff]) +%end + +% Store the near end energy +%aecmStructNew.enerNear(setupStruct.currentBlock) = log(afmicrophone'*afmicrophone); +aecmStructNew.enerNear(setupStruct.currentBlock) = log(sum(afmicrophone)); +% Store the far end energy +%aecmStructNew.enerFar(setupStruct.currentBlock) = log(afTheFarEnd'*afTheFarEnd); +aecmStructNew.enerFar(setupStruct.currentBlock) = log(sum(afTheFarEnd)); + +% Update subbands (We currently use all frequency bins, hence .useSubBand is turned off) +if aecmStructNew.useSubBand + internalIndex = 1; + for kk=1:setupStruct.subBandLength+1 + ySubBand(kk) = mean(afmicrophone(internalIndex:internalIndex+setupStruct.numInBand(kk)-1).^aecmStructNew.bandFactor); + xSubBand(kk) = mean(afTheFarEnd(internalIndex:internalIndex+setupStruct.numInBand(kk)-1).^aecmStructNew.bandFactor); + internalIndex = internalIndex + setupStruct.numInBand(kk); + end +else + ySubBand = afmicrophone.^aecmStructNew.bandFactor; + xSubBand = afTheFarEnd.^aecmStructNew.bandFactor; +end + +% Estimated echo energy +if (aecmStructNew.bandFactor == 1) + %aecmStructNew.enerEcho(setupStruct.currentBlock) = log((aecmStructNew.H.*xSubBand)'*(aecmStructNew.H.*xSubBand)); + %aecmStructNew.enerEchoStored(setupStruct.currentBlock) = log((aecmStructNew.HStored.*xSubBand)'*(aecmStructNew.HStored.*xSubBand)); + aecmStructNew.enerEcho(setupStruct.currentBlock) = log(sum(aecmStructNew.H.*xSubBand)); + aecmStructNew.enerEchoStored(setupStruct.currentBlock) = log(sum(aecmStructNew.HStored.*xSubBand)); +elseif (aecmStructNew.bandFactor == 2) + aecmStructNew.enerEcho(setupStruct.currentBlock) = log(aecmStructNew.H'*xSubBand); + aecmStructNew.enerEchoStored(setupStruct.currentBlock) = log(aecmStructNew.HStored'*xSubBand); +end + +% Last 100 blocks of data, used for plotting +n100 = max(1,setupStruct.currentBlock-99):setupStruct.currentBlock; +enerError = aecmStructNew.enerNear(n100)-aecmStructNew.enerEcho(n100); +enerErrorStored = aecmStructNew.enerNear(n100)-aecmStructNew.enerEchoStored(n100); + +% Store the far end sub band. This is needed if we use LSE instead of NLMS +aecmStructNew.X = [xSubBand aecmStructNew.X(:,1:end-1)]; + +% Update energy levels, which control the VAD +if ((aecmStructNew.enerFar(setupStruct.currentBlock) < aecmStructNew.energyMin) & (aecmStructNew.enerFar(setupStruct.currentBlock) >= aecmStruct.FAR_ENERGY_MIN)) + aecmStructNew.energyMin = aecmStructNew.enerFar(setupStruct.currentBlock); + %aecmStructNew.energyMin = max(aecmStructNew.energyMin,12); + aecmStructNew.energyMin = max(aecmStructNew.energyMin,aecmStruct.FAR_ENERGY_MIN); + aecmStructNew.energyLevel = (aecmStructNew.energyMax-aecmStructNew.energyMin)*aecmStructNew.energyThres+aecmStructNew.energyMin; + aecmStructNew.energyLevelMSE = (aecmStructNew.energyMax-aecmStructNew.energyMin)*aecmStructNew.energyThresMSE+aecmStructNew.energyMin; +end +if (aecmStructNew.enerFar(setupStruct.currentBlock) > aecmStructNew.energyMax) + aecmStructNew.energyMax = aecmStructNew.enerFar(setupStruct.currentBlock); + aecmStructNew.energyLevel = (aecmStructNew.energyMax-aecmStructNew.energyMin)*aecmStructNew.energyThres+aecmStructNew.energyMin; + aecmStructNew.energyLevelMSE = (aecmStructNew.energyMax-aecmStructNew.energyMin)*aecmStructNew.energyThresMSE+aecmStructNew.energyMin; +end + +% Calculate current energy error in near end (estimated echo vs. near end) +dE = aecmStructNew.enerNear(setupStruct.currentBlock)-aecmStructNew.enerEcho(setupStruct.currentBlock); + +%%%%%%%% +% Calculate step size used in LMS algorithm, based on current far end energy and near end energy error (dE) +%%%%%%%% +if setupStruct.stepSize_flag + [mu, aecmStructNew] = calcStepSize(aecmStructNew.enerFar(setupStruct.currentBlock), dE, aecmStructNew, setupStruct.currentBlock, 1); +else + mu = 0.25; +end +aecmStructNew.muLog(setupStruct.currentBlock) = mu; % Store the step size + +% Estimate Echo Spectral Shape +[U, aecmStructNew.H] = fallerEstimator(ySubBand,aecmStructNew.X,aecmStructNew.H,mu); + +%%%%% +% Determine if we should store or restore the channel +%%%%% +if ((setupStruct.currentBlock <= aecmStructNew.convLength) | (~setupStruct.channelUpdate_flag)) + aecmStructNew.HStored = aecmStructNew.H; % Store what you have after startup +elseif ((setupStruct.currentBlock > aecmStructNew.convLength) & (setupStruct.channelUpdate_flag)) + if ((aecmStructNew.enerFar(setupStruct.currentBlock) < aecmStructNew.energyLevelMSE) & (aecmStructNew.enerFar(setupStruct.currentBlock-1) >= aecmStructNew.energyLevelMSE)) + xxx = aecmStructNew.countMseH; + if (xxx > 20) + mseStored = mean(abs(aecmStructNew.enerEchoStored(setupStruct.currentBlock-xxx:setupStruct.currentBlock-1)-aecmStructNew.enerNear(setupStruct.currentBlock-xxx:setupStruct.currentBlock-1))); + mseLatest = mean(abs(aecmStructNew.enerEcho(setupStruct.currentBlock-xxx:setupStruct.currentBlock-1)-aecmStructNew.enerNear(setupStruct.currentBlock-xxx:setupStruct.currentBlock-1))); + %fprintf('Stored: %4f Latest: %4f\n', mseStored, mseLatest) % Uncomment if you want to display the MSE values + if ((mseStored < 0.8*mseLatest) & (aecmStructNew.mseHStoredOld < 0.8*aecmStructNew.mseHLatestOld)) + aecmStructNew.H = aecmStructNew.HStored; + fprintf('Restored H at block %d\n',setupStruct.currentBlock) + elseif (((0.8*mseStored > mseLatest) & (mseLatest < aecmStructNew.mseHThreshold) & (aecmStructNew.mseHLatestOld < aecmStructNew.mseHThreshold)) | (mseStored == Inf)) + aecmStructNew.HStored = aecmStructNew.H; + fprintf('Stored new H at block %d\n',setupStruct.currentBlock) + end + aecmStructNew.mseHStoredOld = mseStored; + aecmStructNew.mseHLatestOld = mseLatest; + end + elseif ((aecmStructNew.enerFar(setupStruct.currentBlock) >= aecmStructNew.energyLevelMSE) & (aecmStructNew.enerFar(setupStruct.currentBlock-1) < aecmStructNew.energyLevelMSE)) + aecmStructNew.countMseH = 1; + elseif (aecmStructNew.enerFar(setupStruct.currentBlock) >= aecmStructNew.energyLevelMSE) + aecmStructNew.countMseH = aecmStructNew.countMseH + 1; + end +end + +%%%%% +% Check delay (calculate the delay offset (if we can)) +% The algorithm is not tuned and should be used with care. It runs separately from Bastiaan's algorithm. +%%%%% +yyy = 31; % Correlation buffer length (currently unfortunately hard coded) +dxxx = 25; % Maximum offset (currently unfortunately hard coded) +if (setupStruct.currentBlock > aecmStructNew.convLength) + if (aecmStructNew.enerFar(setupStruct.currentBlock-(yyy+2*dxxx-1):setupStruct.currentBlock) > aecmStructNew.energyLevelMSE) + for xxx = -dxxx:dxxx + aecmStructNew.delayLatestS(xxx+dxxx+1) = sum(sign(aecmStructNew.enerEcho(setupStruct.currentBlock-(yyy+dxxx-xxx)+1:setupStruct.currentBlock+xxx-dxxx)-mean(aecmStructNew.enerEcho(setupStruct.currentBlock-(yyy++dxxx-xxx)+1:setupStruct.currentBlock+xxx-dxxx))).*sign(aecmStructNew.enerNear(setupStruct.currentBlock-yyy-dxxx+1:setupStruct.currentBlock-dxxx)-mean(aecmStructNew.enerNear(setupStruct.currentBlock-yyy-dxxx+1:setupStruct.currentBlock-dxxx)))); + end + aecmStructNew.newDelayCurve = 1; + end +end +if ((setupStruct.currentBlock > 2*aecmStructNew.convLength) & ~rem(setupStruct.currentBlock,yyy*2) & aecmStructNew.newDelayCurve) + [maxV,maxP] = max(aecmStructNew.delayLatestS); + if ((maxP > 2) & (maxP < 2*dxxx)) + maxVLeft = aecmStructNew.delayLatestS(max(1,maxP-4)); + maxVRight = aecmStructNew.delayLatestS(min(2*dxxx+1,maxP+4)); + %fprintf('Max %d, Left %d, Right %d\n',maxV,maxVLeft,maxVRight) % Uncomment if you want to see max value + if ((maxV > 24) & (maxVLeft < maxV - 10) & (maxVRight < maxV - 10)) + aecmStructNew.feedbackDelay = maxP-dxxx-1; + aecmStructNew.newDelayCurve = 0; + aecmStructNew.feedbackDelayUpdate = 1; + fprintf('Feedback Update at block %d\n',setupStruct.currentBlock) + end + end +end +% End of "Check delay" +%%%%%%%% + +%%%%% +% Calculate suppression gain, based on far end energy and near end energy error (dE) +if (setupStruct.supGain_flag) + [gamma_echo, aecmStructNew.cntIn, aecmStructNew.cntOut] = calcFilterGain(aecmStructNew.enerFar(setupStruct.currentBlock), dE, aecmStructNew, setupStruct.currentBlock, aecmStructNew.convLength, aecmStructNew.cntIn, aecmStructNew.cntOut); +else + gamma_echo = 1; +end +aecmStructNew.gammaLog(setupStruct.currentBlock) = gamma_echo; % Store the gain +gamma_use = gamma_echo; + +% Use the stored channel +U = aecmStructNew.HStored.*xSubBand; + +% compute Wiener filter and suppressor function +Iy = find(ySubBand); +subBandFilter = zeros(size(ySubBand)); +if (aecmStructNew.bandFactor == 2) + subBandFilter(Iy) = (1 - gamma_use*sqrt(U(Iy)./ySubBand(Iy))); % For Faller +else + subBandFilter(Iy) = (1 - gamma_use*(U(Iy)./ySubBand(Iy))); % For COV +end +ix0 = find(subBandFilter < 0); % bounding trick 1 +subBandFilter(ix0) = 0; +ix0 = find(subBandFilter > 1); % bounding trick 1 +subBandFilter(ix0) = 1; + +% Interpolate back to normal frequency bins if we use sub bands +if aecmStructNew.useSubBand + thefilter = interp1(setupStruct.centerFreq,subBandFilter,linspace(0,setupStruct.samplingfreq/2,setupStruct.hsupport1)','nearest'); + testfilter = interp1(setupStruct.centerFreq,subBandFilter,linspace(0,setupStruct.samplingfreq/2,1000),'nearest'); + thefilter(end) = subBandFilter(end); + + internalIndex = 1; + for kk=1:setupStruct.subBandLength+1 + internalIndex:internalIndex+setupStruct.numInBand(kk)-1; + thefilter(internalIndex:internalIndex+setupStruct.numInBand(kk)-1) = subBandFilter(kk); + internalIndex = internalIndex + setupStruct.numInBand(kk); + end +else + thefilter = subBandFilter; + testfilter = subBandFilter; +end + +% Bound the filter +ix0 = find(thefilter < setupStruct.de_echo_bound); % bounding trick 1 +thefilter(ix0) = setupStruct.de_echo_bound; % bounding trick 2 +ix0 = find(thefilter > 1); % bounding in reasonable range +thefilter(ix0) = 1; + +%%%% +% NLP +%%%% +thefmean = mean(thefilter(8:16)); +if (thefmean < 1) + disp(''); +end +aecmStructNew.runningfmean = setupStruct.nl_alpha*aecmStructNew.runningfmean + (1-setupStruct.nl_alpha)*thefmean; +slope0 = 1.0/(1.0 - setupStruct.nlSeverity); % +thegain = max(0.0, min(1.0, slope0*(aecmStructNew.runningfmean - setupStruct.nlSeverity))); +if ~setupStruct.nlp_flag + thegain = 1; +end +% END NONLINEARITY +thefilter = thegain*thefilter; + +%%%% +% The suppression +%%%% +femicrophone = fmicrophone .* thefilter; +% Store the output energy (used for plotting) +%aecmStructNew.enerOut(setupStruct.currentBlock) = log(abs(femicrophone)'*abs(femicrophone)); +aecmStructNew.enerOut(setupStruct.currentBlock) = log(sum(abs(femicrophone))); + +if aecmStructNew.plotIt + figure(13) + subplot(311) + %plot(n100,enerFar(n100),'b-',n100,enerNear(n100),'k--',n100,enerEcho(n100),'r-',[n100(1) n100(end)],[1 1]*vadThNew,'b:',[n100(1) n100(end)],[1 1]*((energyMax-energyMin)/4+energyMin),'r-.',[n100(1) n100(end)],[1 1]*vadNearThNew,'g:',[n100(1) n100(end)],[1 1]*energyMax,'r-.',[n100(1) n100(end)],[1 1]*energyMin,'r-.','LineWidth',2) + plot(n100,aecmStructNew.enerFar(n100),'b-',n100,aecmStructNew.enerNear(n100),'k--',n100,aecmStructNew.enerOut(n100),'r-.',n100,aecmStructNew.enerEcho(n100),'r-',n100,aecmStructNew.enerEchoStored(n100),'c-',[n100(1) n100(end)],[1 1]*((aecmStructNew.energyMax-aecmStructNew.energyMin)/4+aecmStructNew.energyMin),'g-.',[n100(1) n100(end)],[1 1]*aecmStructNew.energyMax,'g-.',[n100(1) n100(end)],[1 1]*aecmStructNew.energyMin,'g-.','LineWidth',2) + %title(['Frame ',int2str(i),' av ',int2str(setupStruct.updateno),' State = ',int2str(speechState),' \mu = ',num2str(mu)]) + title(['\gamma = ',num2str(gamma_echo),' \mu = ',num2str(mu)]) + subplot(312) + %plot(n100,enerError,'b-',[n100(1) n100(end)],[1 1]*vadNearTh,'r:',[n100(1) n100(end)],[-1.5 -1.5]*vadNearTh,'r:','LineWidth',2) + %plot(n100,enerError,'b-',[n100(1) n100(end)],[1 1],'r:',[n100(1) n100(end)],[-2 -2],'r:','LineWidth',2) + plot(n100,enerError,'b-',n100,enerErrorStored,'c-',[n100(1) n100(end)],[1 1]*aecmStructNew.varMean,'k--',[n100(1) n100(end)],[1 1],'r:',[n100(1) n100(end)],[-2 -2],'r:','LineWidth',2) + % Plot mu + %plot(n100,log2(aecmStructNew.muLog(n100)),'b-','LineWidth',2) + %plot(n100,log2(aecmStructNew.HGain(n100)),'b-',[n100(1) n100(end)],[1 1]*log2(sum(aecmStructNew.HStored)),'r:','LineWidth',2) + title(['Block ',int2str(setupStruct.currentBlock),' av ',int2str(setupStruct.updateno)]) + subplot(313) + %plot(n100,enerVar(n100),'b-',[n100(1) n100(end)],[1 1],'r:',[n100(1) n100(end)],[-2 -2],'r:','LineWidth',2) + %plot(n100,enerVar(n100),'b-','LineWidth',2) + % Plot correlation curve + + %plot(-25:25,aecmStructNew.delayStored/max(aecmStructNew.delayStored),'c-',-25:25,aecmStructNew.delayLatest/max(aecmStructNew.delayLatest),'r-',-25:25,(max(aecmStructNew.delayStoredS)-aecmStructNew.delayStoredS)/(max(aecmStructNew.delayStoredS)-min(aecmStructNew.delayStoredS)),'c:',-25:25,(max(aecmStructNew.delayLatestS)-aecmStructNew.delayLatestS)/(max(aecmStructNew.delayLatestS)-min(aecmStructNew.delayLatestS)),'r:','LineWidth',2) + %plot(-25:25,aecmStructNew.delayStored,'c-',-25:25,aecmStructNew.delayLatest,'r-',-25:25,(max(aecmStructNew.delayStoredS)-aecmStructNew.delayStoredS)/(max(aecmStructNew.delayStoredS)-min(aecmStructNew.delayStoredS)),'c:',-25:25,(max(aecmStructNew.delayLatestS)-aecmStructNew.delayLatestS)/(max(aecmStructNew.delayLatestS)-min(aecmStructNew.delayLatestS)),'r:','LineWidth',2) + %plot(-25:25,aecmStructNew.delayLatest,'r-',-25:25,(50-aecmStructNew.delayLatestS)/100,'r:','LineWidth',2) + plot(-25:25,aecmStructNew.delayLatestS,'r:','LineWidth',2) + %plot(-25:25,aecmStructNew.delayStored,'c-',-25:25,aecmStructNew.delayLatest,'r-','LineWidth',2) + plot(0:32,aecmStruct.HStored,'bo-','LineWidth',2) + %title(['\gamma | In = ',int2str(aecmStructNew.muStruct.countInInterval),' | Out High = ',int2str(aecmStructNew.muStruct.countOutHighInterval),' | Out Low = ',int2str(aecmStructNew.muStruct.countOutLowInterval)]) + pause(1) + %if ((setupStruct.currentBlock == 860) | (setupStruct.currentBlock == 420) | (setupStruct.currentBlock == 960)) + if 0%(setupStruct.currentBlock == 960) + figure(60) + plot(n100,aecmStructNew.enerNear(n100),'k--',n100,aecmStructNew.enerEcho(n100),'k:','LineWidth',2) + legend('Near End','Estimated Echo') + title('Signal Energy witH offset compensation') + figure(61) + subplot(211) + stem(sign(aecmStructNew.enerNear(n100)-mean(aecmStructNew.enerNear(n100)))) + title('Near End Energy Pattern (around mean value)') + subplot(212) + stem(sign(aecmStructNew.enerEcho(n100)-mean(aecmStructNew.enerEcho(n100)))) + title('Estimated Echo Energy Pattern (around mean value)') + pause + end + drawnow%,pause +elseif ~rem(setupStruct.currentBlock,100) + fprintf('Block %d of %d\n',setupStruct.currentBlock,setupStruct.updateno) +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/align.m b/src/modules/audio_processing/aecm/main/matlab/matlab/align.m new file mode 100644 index 0000000000..9b9c0baf3b --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/align.m @@ -0,0 +1,98 @@ +function [delayStructNew] = align(xf, yf, delayStruct, i, trueDelay); + +%%%%%%% +% Bastiaan's algorithm copied +%%%%%%% +Ap500 = [1.00, -4.95, 9.801, -9.70299, 4.80298005, -0.9509900499]; +Bp500 = [0.662743088639636, -2.5841655608125, 3.77668102146288, -2.45182477425154, 0.596566274575251, 0.0]; +Ap200 = [1.00, -4.875, 9.50625, -9.26859375, 4.518439453125, -0.881095693359375]; +Bp200 = [0.862545460994275, -3.2832804496114, 4.67892032308828, -2.95798023879133, 0.699796870041299, 0.0]; + +oldMethod = 1; % Turn on or off the old method. The new one is Bastiaan's August 2008 updates +THReSHoLD = 2.0; % ADJUSTABLE threshold factor; 4.0 seems good +%%%%%%%%%%%%%%%%%%% +% use log domain (showed improved performance) +xxf = sqrt(real(xf.*conj(xf))+1e-20); +yyf = sqrt(real(yf.*conj(yf))+1e-20); +delayStruct.sxAll2(:,i) = 20*log10(xxf); +delayStruct.syAll2(:,i) = 20*log10(yyf); + +mD = min(i-1,delayStruct.maxDelayb); +if oldMethod + factor = 1.0; + histLenb = 250; + xthreshold = factor*median(delayStruct.sxAll2(:,i-mD:i),2); + ythreshold = factor*median(delayStruct.syAll2(:,i-mD:i),2); +else + xthreshold = sum(delayStruct.sxAll2(:,i-mD:i),2)/(delayStruct.maxDelayb+1); + + [yout, delayStruct.z200] = filter(Bp200, Ap200, delayStruct.syAll2(:,i), delayStruct.z200, 2); + yout = yout/(delayStruct.maxDelayb+1); + ythreshold = mean(delayStruct.syAll2(:,i-mD:i),2); + ythreshold = yout; +end + +delayStruct.bxspectrum(i) = getBspectrum(delayStruct.sxAll2(:,i), xthreshold, delayStruct.bandfirst, delayStruct.bandlast); +delayStruct.byspectrum(i) = getBspectrum(delayStruct.syAll2(:,i), ythreshold, delayStruct.bandfirst, delayStruct.bandlast); + +delayStruct.bxhist(end-mD:end) = delayStruct.bxspectrum(i-mD:i); + +delayStruct.bcount(:,i) = hisser2(delayStruct.byspectrum(i), flipud(delayStruct.bxhist), delayStruct.bandfirst, delayStruct.bandlast); +[delayStruct.fout(:,i), delayStruct.z500] = filter(Bp500, Ap500, delayStruct.bcount(:,i), delayStruct.z500, 2); +if oldMethod + %delayStruct.new(:,i) = sum(delayStruct.bcount(:,max(1,i-histLenb+1):i),2); % using the history range + tmpVec = [delayStruct.fout(1,i)*ones(2,1); delayStruct.fout(:,i); delayStruct.fout(end,i)*ones(2,1)]; % using the history range + tmpVec = filter(ones(1,5), 1, tmpVec); + delayStruct.new(:,i) = tmpVec(5:end); + %delayStruct.new(:,i) = delayStruct.fout(:,i); % using the history range +else + [delayStruct.fout(:,i), delayStruct.z500] = filter(Bp500, Ap500, delayStruct.bcount(:,i), delayStruct.z500, 2); + % NEW CODE + delayStruct.new(:,i) = filter([-1,-2,1,4,1,-2,-1], 1, delayStruct.fout(:,i)); %remv smth component + delayStruct.new(1:end-3,i) = delayStruct.new(1+3:end,i); + delayStruct.new(1:6,i) = 0.0; + delayStruct.new(end-6:end,i) = 0.0; % ends are no good +end +[valuen, tempdelay] = min(delayStruct.new(:,i)); % find minimum +if oldMethod + threshold = valuen + (max(delayStruct.new(:,i)) - valuen)/4; + thIndex = find(delayStruct.new(:,i) <= threshold); + if (i > 1) + delayDiff = abs(delayStruct.delay(i-1)-tempdelay+1); + if (delayStruct.oneGoodEstimate & (max(diff(thIndex)) > 1) & (delayDiff < 10)) + % We consider this minimum to be significant, hence update the delay + delayStruct.delay(i) = tempdelay; + elseif (~delayStruct.oneGoodEstimate & (max(diff(thIndex)) > 1)) + delayStruct.delay(i) = tempdelay; + if (i > histLenb) + delayStruct.oneGoodEstimate = 1; + end + else + delayStruct.delay(i) = delayStruct.delay(i-1); + end + else + delayStruct.delay(i) = tempdelay; + end +else + threshold = THReSHoLD*std(delayStruct.new(:,i)); % set updata threshold + if ((-valuen > threshold) | (i < delayStruct.smlength)) % see if you want to update delay + delayStruct.delay(i) = tempdelay; + else + delayStruct.delay(i) = delayStruct.delay(i-1); + end + % END NEW CODE +end +delayStructNew = delayStruct; + +% administrative and plotting stuff +if( 0) + figure(10); + plot([1:length(delayStructNew.new(:,i))],delayStructNew.new(:,i),trueDelay*[1 1],[min(delayStructNew.new(:,i)),max(delayStructNew.new(:,i))],'r',[1 length(delayStructNew.new(:,i))],threshold*[1 1],'r:', 'LineWidth',2); + %plot([1:length(delayStructNew.bcount(:,i))],delayStructNew.bcount(:,i),trueDelay*[1 1],[min(delayStructNew.bcount(:,i)),max(delayStructNew.bcount(:,i))],'r','LineWidth',2); + %plot([thedelay,thedelay],[min(fcount(:,i)),max(fcount(:,i))],'r'); + %title(sprintf('bin count and known delay at time %5.1f s\n',(i-1)*(support/(fs*oversampling)))); + title(delayStructNew.oneGoodEstimate) + xlabel('delay in frames'); + %hold off; + drawnow +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/calcFilterGain.m b/src/modules/audio_processing/aecm/main/matlab/matlab/calcFilterGain.m new file mode 100644 index 0000000000..a09a7f2225 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/calcFilterGain.m @@ -0,0 +1,88 @@ +function [gam, cntIn2, cntOut2] = calcFilterGain(energy, dE, aecmStruct, t, T, cntIn, cntOut) + +defaultLevel = 1.2; +cntIn2 = cntIn; +cntOut2 = cntOut; +if (t < T) + gam = 1; +else + dE1 = -5; + dE2 = 1; + gamMid = 0.2; + gam = max(0,min((energy - aecmStruct.energyMin)/(aecmStruct.energyLevel - aecmStruct.energyMin), 1-(1-gamMid)*(aecmStruct.energyMax-energy)/(aecmStruct.energyMax-aecmStruct.energyLevel))); + + dEOffset = -0.5; + dEWidth = 1.5; + %gam2 = max(1,2-((dE-dEOffset)/(dE2-dEOffset)).^2); + gam2 = 1+(abs(dE-dEOffset)<(dE2-dEOffset)); + + gam = gam*gam2; + + + if (energy < aecmStruct.energyLevel) + gam = 0; + else + gam = defaultLevel; + end + dEVec = aecmStruct.enerNear(t-63:t)-aecmStruct.enerEcho(t-63:t); + %dEVec = aecmStruct.enerNear(t-20:t)-aecmStruct.enerEcho(t-20:t); + numCross = 0; + currentState = 0; + for ii=1:64 + if (currentState == 0) + currentState = (dEVec(ii) > dE2) - (dEVec(ii) < -2); + elseif ((currentState == 1) & (dEVec(ii) < -2)) + numCross = numCross + 1; + currentState = -1; + elseif ((currentState == -1) & (dEVec(ii) > dE2)) + numCross = numCross + 1; + currentState = 1; + end + end + gam = max(0, gam - numCross/25); + gam = 1; + + ener_A = 1; + ener_B = 0.8; + ener_C = aecmStruct.energyLevel + (aecmStruct.energyMax-aecmStruct.energyLevel)/5; + dE_A = 4;%2; + dE_B = 3.6;%1.8; + dE_C = 0.9*dEWidth; + dE_D = 1; + timeFactorLength = 10; + ddE = abs(dE-dEOffset); + if (energy < aecmStruct.energyLevel) + gam = 0; + else + gam = 1; + gam2 = max(0, min(ener_B*(energy-aecmStruct.energyLevel)/(ener_C-aecmStruct.energyLevel), ener_B+(ener_A-ener_B)*(energy-ener_C)/(aecmStruct.energyMax-ener_C))); + if (ddE < dEWidth) + % Update counters + cntIn2 = cntIn2 + 1; + if (cntIn2 > 2) + cntOut2 = 0; + end + gam3 = max(dE_D, min(dE_A-(dE_A-dE_B)*(ddE/dE_C), dE_D+(dE_B-dE_D)*(dEWidth-ddE)/(dEWidth-dE_C))); + gam3 = dE_A; + else + % Update counters + cntOut2 = cntOut2 + 1; + if (cntOut2 > 2) + cntIn2 = 0; + end + %gam2 = 1; + gam3 = dE_D; + end + timeFactor = min(1, cntIn2/timeFactorLength); + gam = gam*(1-timeFactor) + timeFactor*gam2*gam3; + end + %gam = gam/floor(numCross/2+1); +end +if isempty(gam) + numCross + timeFactor + cntIn2 + cntOut2 + gam2 + gam3 +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/calcStepSize.m b/src/modules/audio_processing/aecm/main/matlab/matlab/calcStepSize.m new file mode 100644 index 0000000000..ae1365fa48 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/calcStepSize.m @@ -0,0 +1,105 @@ +function [mu, aecmStructNew] = calcStepSize(energy, dE, aecmStruct, t, logscale) + +if (nargin < 4) + t = 1; + logscale = 1; +elseif (nargin == 4) + logscale = 1; +end +T = aecmStruct.convLength; + +if logscale + currentMuMax = aecmStruct.MU_MIN + (aecmStruct.MU_MAX-aecmStruct.MU_MIN)*min(t,T)/T; + if (aecmStruct.energyMin >= aecmStruct.energyMax) + mu = aecmStruct.MU_MIN; + else + mu = (energy - aecmStruct.energyMin)/(aecmStruct.energyMax - aecmStruct.energyMin)*(currentMuMax-aecmStruct.MU_MIN) + aecmStruct.MU_MIN; + end + mu = 2^mu; + if (energy < aecmStruct.energyLevel) + mu = 0; + end +else + muMin = 0; + muMax = 0.5; + currentMuMax = muMin + (muMax-muMin)*min(t,T)/T; + if (aecmStruct.energyMin >= aecmStruct.energyMax) + mu = muMin; + else + mu = (energy - aecmStruct.energyMin)/(aecmStruct.energyMax - aecmStruct.energyMin)*(currentMuMax-muMin) + muMin; + end +end +dE2 = 1; +dEOffset = -0.5; +offBoost = 5; +if (mu > 0) + if (abs(dE-aecmStruct.ENERGY_DEV_OFFSET) > aecmStruct.ENERGY_DEV_TOL) + aecmStruct.muStruct.countInInterval = 0; + else + aecmStruct.muStruct.countInInterval = aecmStruct.muStruct.countInInterval + 1; + end + if (dE < aecmStruct.ENERGY_DEV_OFFSET - aecmStruct.ENERGY_DEV_TOL) + aecmStruct.muStruct.countOutLowInterval = aecmStruct.muStruct.countOutLowInterval + 1; + else + aecmStruct.muStruct.countOutLowInterval = 0; + end + if (dE > aecmStruct.ENERGY_DEV_OFFSET + aecmStruct.ENERGY_DEV_TOL) + aecmStruct.muStruct.countOutHighInterval = aecmStruct.muStruct.countOutHighInterval + 1; + else + aecmStruct.muStruct.countOutHighInterval = 0; + end +end +muVar = 2^min(-3,5/50*aecmStruct.muStruct.countInInterval-3); +muOff = 2^max(offBoost,min(0,offBoost*(aecmStruct.muStruct.countOutLowInterval-aecmStruct.muStruct.minOutLowInterval)/(aecmStruct.muStruct.maxOutLowInterval-aecmStruct.muStruct.minOutLowInterval))); + +muLow = 1/64; +muVar = 1; +if (t < 2*T) + muDT = 1; + muVar = 1; + mdEVec = 0; + numCross = 0; +else + muDT = min(1,max(muLow,1-(1-muLow)*(dE-aecmStruct.ENERGY_DEV_OFFSET)/aecmStruct.ENERGY_DEV_TOL)); + dEVec = aecmStruct.enerNear(t-63:t)-aecmStruct.enerEcho(t-63:t); + %dEVec = aecmStruct.enerNear(t-20:t)-aecmStruct.enerEcho(t-20:t); + numCross = 0; + currentState = 0; + for ii=1:64 + if (currentState == 0) + currentState = (dEVec(ii) > dE2) - (dEVec(ii) < -2); + elseif ((currentState == 1) & (dEVec(ii) < -2)) + numCross = numCross + 1; + currentState = -1; + elseif ((currentState == -1) & (dEVec(ii) > dE2)) + numCross = numCross + 1; + currentState = 1; + end + end + + %logicDEVec = (dEVec > dE2) - (dEVec < -2); + %numCross = sum(abs(diff(logicDEVec))); + %mdEVec = mean(abs(dEVec-dEOffset)); + %mdEVec = mean(abs(dEVec-mean(dEVec))); + %mdEVec = max(dEVec)-min(dEVec); + %if (mdEVec > 4)%1.5) + % muVar = 0; + %end + muVar = 2^(-floor(numCross/2)); + muVar = 2^(-numCross); +end +%muVar = 1; + + +% if (eStd > (dE2-dEOffset)) +% muVar = 1/8; +% else +% muVar = 1; +% end + +%mu = mu*muDT*muVar*muOff; +mu = mu*muDT*muVar; +mu = min(mu,0.25); +aecmStructNew = aecmStruct; +%aecmStructNew.varMean = mdEVec; +aecmStructNew.varMean = numCross; diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/fallerEstimator.m b/src/modules/audio_processing/aecm/main/matlab/matlab/fallerEstimator.m new file mode 100644 index 0000000000..d038b519c0 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/fallerEstimator.m @@ -0,0 +1,42 @@ +function [U, Hnew] = fallerEstimator(Y, X, H, mu) + +% Near end signal is stacked frame by frame columnwise in matrix Y and far end in X +% +% Possible estimation procedures are +% 1) LSE +% 2) NLMS +% 3) Separated numerator and denomerator filters +regParam = 1; +[numFreqs, numFrames] = size(Y); +[numFreqs, Q] = size(X); +U = zeros(numFreqs, 1); + +if ((nargin == 3) | (nargin == 5)) + dtd = 0; +end +if (nargin == 4) + dtd = H; +end +Emax = 7; +dEH = Emax-sum(sum(H)); +nu = 2*mu; +% if (nargin < 5) +% H = zeros(numFreqs, Q); +% for kk = 1:numFreqs +% Xmatrix = hankel(X(kk,1:Q),X(kk,Q:end)); +% y = Y(kk,1:end-Q+1)'; +% H(kk,:) = (y'*Xmatrix')*inv(Xmatrix*Xmatrix'+regParam); +% U(kk,1) = H(kk,:)*Xmatrix(:,1); +% end +% else + for kk = 1:numFreqs + x = X(kk,1:Q)'; + y = Y(kk,1); + Htmp = mu*(y-H(kk,:)*x)/(x'*x+regParam)*x; + %Htmp = (mu*(y-H(kk,:)*x)/(x'*x+regParam) - nu/dEH)*x; + H(kk,:) = H(kk,:) + Htmp'; + U(kk,1) = H(kk,:)*x; + end +% end + +Hnew = H; diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/getBspectrum.m b/src/modules/audio_processing/aecm/main/matlab/matlab/getBspectrum.m new file mode 100644 index 0000000000..a4a533d600 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/getBspectrum.m @@ -0,0 +1,22 @@ +function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast) +% function bspectrum=getBspectrum(ps,threshold,bandfirst,bandlast) +% compute binary spectrum using threshold spectrum as pivot +% bspectrum = binary spectrum (binary) +% ps=current power spectrum (float) +% threshold=threshold spectrum (float) +% bandfirst = first band considered +% bandlast = last band considered + +% initialization stuff + if( length(ps)<bandlast | bandlast>32 | length(ps)~=length(threshold)) + error('BinDelayEst:spectrum:invalid','Dimensionality error'); +end + +% get current binary spectrum +diff = ps - threshold; +bspectrum=uint32(0); +for(i=bandfirst:bandlast) + if( diff(i)>0 ) + bspectrum = bitset(bspectrum,i); + end +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/hisser2.m b/src/modules/audio_processing/aecm/main/matlab/matlab/hisser2.m new file mode 100644 index 0000000000..5a414f9da8 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/hisser2.m @@ -0,0 +1,21 @@ +function bcount=hisser2(bs,bsr,bandfirst,bandlast) +% function bcount=hisser(bspectrum,bandfirst,bandlast) +% histogram for the binary spectra +% bcount= array of bit counts +% bs=binary spectrum (one int32 number each) +% bsr=reference binary spectra (one int32 number each) +% blockSize = histogram over blocksize blocks +% bandfirst = first band considered +% bandlast = last band considered + +% weight all delays equally +maxDelay = length(bsr); + +% compute counts (two methods; the first works better and is operational) +bcount=zeros(maxDelay,1); +for(i=1:maxDelay) + % the delay should have low count for low-near&high-far and high-near&low-far + bcount(i)= sum(bitget(bitxor(bs,bsr(i)),bandfirst:bandlast)); + % the delay should have low count for low-near&high-far (works less well) +% bcount(i)= sum(bitget(bitand(bsr(i),bitxor(bs,bsr(i))),bandfirst:bandlast)); +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/mainProgram.m b/src/modules/audio_processing/aecm/main/matlab/matlab/mainProgram.m new file mode 100644 index 0000000000..eeb2aaa79c --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/mainProgram.m @@ -0,0 +1,283 @@ +useHTC = 1; % Set this if you want to run a single file and set file names below. Otherwise use simEnvironment to run from several scenarios in a row +delayCompensation_flag = 0; % Set this flag to one if you want to turn on the delay compensation/enhancement +global FARENDFFT; +global NEARENDFFT; +global F; + +if useHTC +% fid=fopen('./htcTouchHd/nb/aecFar.pcm'); xFar=fread(fid,'short'); fclose(fid); +% fid=fopen('./htcTouchHd/nb/aecNear.pcm'); yNear=fread(fid,'short'); fclose(fid); +% fid=fopen('./samsungBlackjack/nb/aecFar.pcm'); xFar=fread(fid,'short'); fclose(fid); +% fid=fopen('./samsungBlackjack/nb/aecNear.pcm'); yNear=fread(fid,'short'); fclose(fid); +% fid=fopen('aecFarPoor.pcm'); xFar=fread(fid,'short'); fclose(fid); +% fid=fopen('aecNearPoor.pcm'); yNear=fread(fid,'short'); fclose(fid); +% fid=fopen('out_aes.pcm'); outAES=fread(fid,'short'); fclose(fid); + fid=fopen('aecFar4.pcm'); xFar=fread(fid,'short'); fclose(fid); + fid=fopen('aecNear4.pcm'); yNear=fread(fid,'short'); fclose(fid); + yNearSpeech = zeros(size(xFar)); + fs = 8000; + frameSize = 64; +% frameSize = 128; + fs = 16000; +% frameSize = 256; +%F = load('fftValues.txt'); +%FARENDFFT = F(:,1:33); +%NEARENDFFT = F(:,34:66); + +else + loadFileFar = [speakerType, '_s_',scenario,'_far_b.wav']; + [xFar,fs,nbits] = wavread(loadFileFar); + xFar = xFar*2^(nbits-1); + loadFileNear = [speakerType, '_s_',scenario,'_near_b.wav']; + [yNear,fs,nbits] = wavread(loadFileNear); + yNear = yNear*2^(nbits-1); + loadFileNearSpeech = [speakerType, '_s_',scenario,'_nearSpeech_b.wav']; + [yNearSpeech,fs,nbits] = wavread(loadFileNearSpeech); + yNearSpeech = yNearSpeech*2^(nbits-1); + frameSize = 256; +end + +dtRegions = []; + +% General settings for the AECM +setupStruct = struct(... + 'stepSize_flag', 1,... % This flag turns on the step size calculation. If turned off, mu = 0.25. + 'supGain_flag', 0,... % This flag turns on the suppression gain calculation. If turned off, gam = 1. + 'channelUpdate_flag', 0,... % This flag turns on the channel update. If turned off, H is updated for convLength and then kept constant. + 'nlp_flag', 0,... % Turn on/off NLP + 'withVAD_flag', 0,... % Turn on/off NLP + 'useSubBand', 0,... % Set to 1 if to use subBands + 'useDelayEstimation', 1,... % Set to 1 if to use delay estimation + 'support', frameSize,... % # of samples per frame + 'samplingfreq',fs,... % Sampling frequency + 'oversampling', 2,... % Overlap between blocks/frames + 'updatel', 0,... % # of samples between blocks + 'hsupport1', 0,... % # of bins in frequency domain + 'factor', 0,... % synthesis window amplification + 'tlength', 0,... % # of samples of entire file + 'updateno', 0,... % # of updates + 'nb', 1,... % # of blocks + 'currentBlock', 0,... % + 'win', zeros(frameSize,1),...% Window to apply for fft and synthesis + 'avtime', 1,... % Time (in sec.) to perform averaging + 'estLen', 0,... % Averaging in # of blocks + 'A_GAIN', 10.0,... % + 'suppress_overdrive', 1.0,... % overdrive factor for suppression 1.4 is good + 'gamma_echo', 1.0,... % same as suppress_overdrive but at different place + 'de_echo_bound', 0.0,... % + 'nl_alpha', 0.4,... % memory; seems not very critical + 'nlSeverity', 0.2,... % nonlinearity severity: 0 does nothing; 1 suppresses all + 'numInBand', [],... % # of frequency bins in resp. subBand + 'centerFreq', [],... % Center frequency of resp. subBand + 'dtRegions', dtRegions,... % Regions where we have DT + 'subBandLength', frameSize/2);%All bins + %'subBandLength', 11); %Something's wrong when subBandLength even + %'nl_alpha', 0.8,... % memory; seems not very critical + +delayStruct = struct(... + 'bandfirst', 8,... + 'bandlast', 25,... + 'smlength', 600,... + 'maxDelay', 0.4,... + 'oneGoodEstimate', 0,... + 'delayAdjust', 0,... + 'maxDelayb', 0); +% More parameters in delayStruct are constructed in "updateSettings" below + +% Make struct settings +[setupStruct, delayStruct] = updateSettings(yNear, xFar, setupStruct, delayStruct); +setupStruct.numInBand = ones(setupStruct.hsupport1,1); + +Q = 1; % Time diversity in channel +% General settings for the step size calculation +muStruct = struct(... + 'countInInterval', 0,... + 'countOutHighInterval', 0,... + 'countOutLowInterval', 0,... + 'minInInterval', 50,... + 'minOutHighInterval', 10,... + 'minOutLowInterval', 10,... + 'maxOutLowInterval', 50); +% General settings for the AECM +aecmStruct = struct(... + 'plotIt', 0,... % Set to 0 to turn off plotting + 'useSubBand', 0,... + 'bandFactor', 1,... + 'H', zeros(setupStruct.subBandLength+1,Q),... + 'HStored', zeros(setupStruct.subBandLength+1,Q),... + 'X', zeros(setupStruct.subBandLength+1,Q),... + 'energyThres', 0.28,... + 'energyThresMSE', 0.4,... + 'energyMin', inf,... + 'energyMax', -inf,... + 'energyLevel', 0,... + 'energyLevelMSE', 0,... + 'convLength', 100,... + 'gammaLog', ones(setupStruct.updateno,1),... + 'muLog', ones(setupStruct.updateno,1),... + 'enerFar', zeros(setupStruct.updateno,1),... + 'enerNear', zeros(setupStruct.updateno,1),... + 'enerEcho', zeros(setupStruct.updateno,1),... + 'enerEchoStored', zeros(setupStruct.updateno,1),... + 'enerOut', zeros(setupStruct.updateno,1),... + 'runningfmean', 0,... + 'muStruct', muStruct,... + 'varMean', 0,... + 'countMseH', 0,... + 'mseHThreshold', 1.1,... + 'mseHStoredOld', inf,... + 'mseHLatestOld', inf,... + 'delayLatestS', zeros(1,51),... + 'feedbackDelay', 0,... + 'feedbackDelayUpdate', 0,... + 'cntIn', 0,... + 'cntOut', 0,... + 'FAR_ENERGY_MIN', 1,... + 'ENERGY_DEV_OFFSET', 0.5,... + 'ENERGY_DEV_TOL', 1.5,... + 'MU_MIN', -16,... + 'MU_MAX', -2,... + 'newDelayCurve', 0); + +% Adjust speech signals +xFar = [zeros(setupStruct.hsupport1-1,1);xFar(1:setupStruct.tlength)]; +yNear = [zeros(setupStruct.hsupport1-1,1);yNear(1:setupStruct.tlength)]; +yNearSpeech = [zeros(setupStruct.hsupport1-1,1);yNearSpeech(1:setupStruct.tlength)]; +xFar = xFar(1:setupStruct.tlength); +yNear = yNear(1:setupStruct.tlength); + +% Set figure settings +if aecmStruct.plotIt + figure(13) + set(gcf,'doublebuffer','on') +end +%%%%%%%%%% +% Here starts the algorithm +% Dividing into frames and then estimating the near end speech +%%%%%%%%%% +fTheFarEnd = complex(zeros(setupStruct.hsupport1,1)); +afTheFarEnd = zeros(setupStruct.hsupport1,setupStruct.updateno+1); +fFar = zeros(setupStruct.hsupport1,setupStruct.updateno+1); +fmicrophone = complex(zeros(setupStruct.hsupport1,1)); +afmicrophone = zeros(setupStruct.hsupport1,setupStruct.updateno+1); +fNear = zeros(setupStruct.hsupport1,setupStruct.updateno+1); +femicrophone = complex(zeros(setupStruct.hsupport1,1)); +emicrophone = zeros(setupStruct.tlength,1); + +if (setupStruct.useDelayEstimation == 2) + delSamples = [1641 1895 2032 1895 2311 2000 2350 2222 NaN 2332 2330 2290 2401 2415 NaN 2393 2305 2381 2398]; + delBlocks = round(delSamples/setupStruct.updatel); + delStarts = floor([25138 46844 105991 169901 195739 218536 241803 333905 347703 362660 373753 745135 765887 788078 806257 823835 842443 860139 881869]/setupStruct.updatel); +else + delStarts = []; +end + +for i=1:setupStruct.updateno + setupStruct.currentBlock = i; + + sb = (i-1)*setupStruct.updatel + 1; + se = sb + setupStruct.support - 1; + + %%%%%%% + % Analysis FFTs + %%%%%%% + % Far end signal + temp = fft(setupStruct.win .* xFar(sb:se))/frameSize; + fTheFarEnd = temp(1:setupStruct.hsupport1); + afTheFarEnd(:,i) = abs(fTheFarEnd); + fFar(:,i) = fTheFarEnd; + % Near end signal + temp = fft(setupStruct.win .* yNear(sb:se))/frameSize;%,pause + fmicrophone = temp(1:setupStruct.hsupport1); + afmicrophone(:,i) = abs(fmicrophone); + fNear(:,i) = fmicrophone; + %abs(fmicrophone),pause + % The true near end speaker (if we have such info) + temp = fft(setupStruct.win .* yNearSpeech(sb:se)); + aftrueSpeech = abs(temp(1:setupStruct.hsupport1)); + + if(i == 1000) + %break; + end + + % Perform delay estimation + if (setupStruct.useDelayEstimation == 1) + % Delay Estimation + delayStruct = align(fTheFarEnd, fmicrophone, delayStruct, i); + %delayStruct.delay(i) = 39;%19; + idel = max(i - delayStruct.delay(i) + 1,1); + + if delayCompensation_flag + % If we have a new delay estimate from Bastiaan's alg. update the offset + if (delayStruct.delay(i) ~= delayStruct.delay(max(1,i-1))) + delayStruct.delayAdjust = delayStruct.delayAdjust + delayStruct.delay(i) - delayStruct.delay(i-1); + end + % Store the compensated delay + delayStruct.delayNew(i) = delayStruct.delay(i) - delayStruct.delayAdjust; + if (delayStruct.delayNew(i) < 1) + % Something's wrong + pause,break + end + % Compensate with the offset estimate + idel = idel + delayStruct.delayAdjust; + end + if 0%aecmStruct.plotIt + figure(1) + plot(1:i,delayStruct.delay(1:i),'k:',1:i,delayStruct.delayNew(1:i),'k--','LineWidth',2),drawnow + end + elseif (setupStruct.useDelayEstimation == 2) + % Use "manual delay" + delIndex = find(delStarts<i); + if isempty(delIndex) + idel = i; + else + idel = i - delBlocks(max(delIndex)); + if isnan(idel) + idel = i - delBlocks(max(delIndex)-1); + end + end + else + % No delay estimation + %idel = max(i - 18, 1); + idel = max(i - 50, 1); + end + + %%%%%%%% + % This is the AECM algorithm + % + % Output is the new frequency domain signal (hopefully) echo compensated + %%%%%%%% + [femicrophone, aecmStruct] = AECMobile(fmicrophone, afTheFarEnd(:,idel), setupStruct, aecmStruct); + %[femicrophone, aecmStruct] = AECMobile(fmicrophone, FARENDFFT(idel,:)'/2^F(idel,end-1), setupStruct, aecmStruct); + + if aecmStruct.feedbackDelayUpdate + % If the feedback tells us there is a new offset out there update the enhancement + delayStruct.delayAdjust = delayStruct.delayAdjust + aecmStruct.feedbackDelay; + aecmStruct.feedbackDelayUpdate = 0; + end + + % reconstruction; first make spectrum odd + temp = [femicrophone; flipud(conj(femicrophone(2:(setupStruct.hsupport1-1))))]; + emicrophone(sb:se) = emicrophone(sb:se) + setupStruct.factor * setupStruct.win .* real(ifft(temp))*frameSize; + if max(isnan(emicrophone(sb:se))) + % Something's wrong with the output at block i + i + break + end +end + + +if useHTC + fid=fopen('aecOutMatlabC.pcm','w');fwrite(fid,int16(emicrophone),'short');fclose(fid); + %fid=fopen('farendFFT.txt','w');fwrite(fid,int16(afTheFarEnd(:)),'short');fclose(fid); + %fid=fopen('farendFFTreal.txt','w');fwrite(fid,int16(imag(fFar(:))),'short');fclose(fid); + %fid=fopen('farendFFTimag.txt','w');fwrite(fid,int16(real(fFar(:))),'short');fclose(fid); + %fid=fopen('nearendFFT.txt','w');fwrite(fid,int16(afmicrophone(:)),'short');fclose(fid); + %fid=fopen('nearendFFTreal.txt','w');fwrite(fid,int16(real(fNear(:))),'short');fclose(fid); + %fid=fopen('nearendFFTimag.txt','w');fwrite(fid,int16(imag(fNear(:))),'short');fclose(fid); +end +if useHTC + %spclab(setupStruct.samplingfreq,xFar,yNear,emicrophone) +else + spclab(setupStruct.samplingfreq,xFar,yNear,emicrophone,yNearSpeech) +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/simEnvironment.m b/src/modules/audio_processing/aecm/main/matlab/matlab/simEnvironment.m new file mode 100644 index 0000000000..3ebe701dfd --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/simEnvironment.m @@ -0,0 +1,15 @@ +speakerType = 'fm'; +%for k=2:5 +%for k=[2 4 5] +for k=3 + scenario = int2str(k); + fprintf('Current scenario: %d\n',k) + mainProgram + %saveFile = [speakerType, '_s_',scenario,'_delayEst_v2_vad_man.wav']; + %wavwrite(emic,fs,nbits,saveFile); + %saveFile = ['P:\Engineering_share\BjornV\AECM\',speakerType, '_s_',scenario,'_delayEst_v2_vad_man.pcm']; + %saveFile = [speakerType, '_s_',scenario,'_adaptMu_adaptGamma_withVar_gammFilt_HSt.pcm']; + saveFile = ['scenario_',scenario,'_090417_backupH_nlp.pcm']; + fid=fopen(saveFile,'w');fwrite(fid,int16(emicrophone),'short');fclose(fid); + %pause +end diff --git a/src/modules/audio_processing/aecm/main/matlab/matlab/updateSettings.m b/src/modules/audio_processing/aecm/main/matlab/matlab/updateSettings.m new file mode 100644 index 0000000000..c805f1d09f --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/matlab/updateSettings.m @@ -0,0 +1,94 @@ +function [setupStructNew, delayStructNew] = updateSettings(microphone, TheFarEnd, setupStruct, delayStruct); + +% other, constants +setupStruct.hsupport1 = setupStruct.support/2 + 1; +setupStruct.factor = 2 / setupStruct.oversampling; +setupStruct.updatel = setupStruct.support/setupStruct.oversampling; +setupStruct.estLen = round(setupStruct.avtime * setupStruct.samplingfreq/setupStruct.updatel); + +% compute some constants +blockLen = setupStruct.support/setupStruct.oversampling; +delayStruct.maxDelayb = floor(setupStruct.samplingfreq*delayStruct.maxDelay/setupStruct.updatel); % in blocks + +%input +tlength = min([length(microphone),length(TheFarEnd)]); +updateno = floor(tlength/setupStruct.updatel); +setupStruct.tlength = setupStruct.updatel*updateno; +setupStruct.updateno = updateno - setupStruct.oversampling + 1; + +% signal length +n = floor(min([length(TheFarEnd), length(microphone)])/setupStruct.support)*setupStruct.support; +setupStruct.nb = n/blockLen - setupStruct.oversampling + 1; % in blocks + +setupStruct.win = sqrt([0 ; hanning(setupStruct.support-1)]); + +% Construct filterbank in Bark-scale + +K = setupStruct.subBandLength; %Something's wrong when K even +erbs = 21.4*log10(0.00437*setupStruct.samplingfreq/2+1); +fe = (10.^((0:K)'*erbs/K/21.4)-1)/0.00437; +setupStruct.centerFreq = fe; +H = diag(ones(1,K-1))+diag(ones(1,K-2),-1); +Hinv = inv(H); +aty = 2*Hinv(end,:)*fe(2:end-1); +boundary = aty - (setupStruct.samplingfreq/2 + fe(end-1))/2; +if rem(K,2) + x1 = min([fe(2)/2, -boundary]); +else + x1 = max([0, boundary]); +end +%x1 +g = fe(2:end-1); +g(1) = g(1) - x1/2; +x = 2*Hinv*g; +x = [x1;x]; +%figure(42), clf +xy = zeros((K+1)*4,1); +yy = zeros((K+1)*4,1); +xy(1:4) = [fe(1) fe(1) x(1) x(1)]'; +yy(1:4) = [0 1 1 0]'/x(1); +for kk=2:K + xy((kk-1)*4+(1:4)) = [x(kk-1) x(kk-1) x(kk) x(kk)]'; + yy((kk-1)*4+(1:4)) = [0 1 1 0]'/(x(kk)-x(kk-1)); +end +xy(end-3:end) = [x(K) x(K) fe(end) fe(end)]'; +yy(end-3:end) = [0 1 1 0]'/(fe(end)*2-2*x(K)); +%plot(xy,yy,'LineWidth',2) +%fill(xy,yy,'y') + +x = [0;x]; +xk = x*setupStruct.hsupport1/setupStruct.samplingfreq*2; +%setupStruct.erbBoundaries = xk; +numInBand = zeros(length(xk),1); +xh = (0:setupStruct.hsupport1-1); + +for kk=1:length(xk) + if (kk==length(xk)) + numInBand(kk) = length(find(xh>=xk(kk))); + else + numInBand(kk) = length(intersect(find(xh>=xk(kk)),find(xh<xk(kk+1)))); + end +end +setupStruct.numInBand = numInBand; + +setupStructNew = setupStruct; + +delayStructNew = struct(... + 'sxAll2',zeros(setupStructNew.hsupport1,setupStructNew.nb),... + 'syAll2',zeros(setupStructNew.hsupport1,setupStructNew.nb),... + 'z200',zeros(5,setupStructNew.hsupport1),... + 'z500',zeros(5,delayStruct.maxDelayb+1),... + 'bxspectrum',uint32(zeros(setupStructNew.nb,1)),... + 'byspectrum',uint32(zeros(setupStructNew.nb,1)),... + 'bandfirst',delayStruct.bandfirst,'bandlast',delayStruct.bandlast,... + 'bxhist',uint32(zeros(delayStruct.maxDelayb+1,1)),... + 'bcount',zeros(1+delayStruct.maxDelayb,setupStructNew.nb),... + 'fout',zeros(1+delayStruct.maxDelayb,setupStructNew.nb),... + 'new',zeros(1+delayStruct.maxDelayb,setupStructNew.nb),... + 'smlength',delayStruct.smlength,... + 'maxDelay', delayStruct.maxDelay,... + 'maxDelayb', delayStruct.maxDelayb,... + 'oneGoodEstimate', 0,... + 'delayAdjust', 0,... + 'delayNew',zeros(setupStructNew.nb,1),... + 'delay',zeros(setupStructNew.nb,1)); diff --git a/src/modules/audio_processing/aecm/main/matlab/waitbar_j.m b/src/modules/audio_processing/aecm/main/matlab/waitbar_j.m new file mode 100644 index 0000000000..50b9ccf309 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/matlab/waitbar_j.m @@ -0,0 +1,234 @@ +function fout = waitbar_j(x,whichbar, varargin) +%WAITBAR Display wait bar. +% H = WAITBAR(X,'title', property, value, property, value, ...) +% creates and displays a waitbar of fractional length X. The +% handle to the waitbar figure is returned in H. +% X should be between 0 and 1. Optional arguments property and +% value allow to set corresponding waitbar figure properties. +% Property can also be an action keyword 'CreateCancelBtn', in +% which case a cancel button will be added to the figure, and +% the passed value string will be executed upon clicking on the +% cancel button or the close figure button. +% +% WAITBAR(X) will set the length of the bar in the most recently +% created waitbar window to the fractional length X. +% +% WAITBAR(X,H) will set the length of the bar in waitbar H +% to the fractional length X. +% +% WAITBAR(X,H,'updated title') will update the title text in +% the waitbar figure, in addition to setting the fractional +% length to X. +% +% WAITBAR is typically used inside a FOR loop that performs a +% lengthy computation. A sample usage is shown below: +% +% h = waitbar(0,'Please wait...'); +% for i=1:100, +% % computation here % +% waitbar(i/100,h) +% end +% close(h) + +% Clay M. Thompson 11-9-92 +% Vlad Kolesnikov 06-7-99 +% Copyright 1984-2001 The MathWorks, Inc. +% $Revision: 1.22 $ $Date: 2001/04/15 12:03:29 $ + +if nargin>=2 + if ischar(whichbar) + type=2; %we are initializing + name=whichbar; + elseif isnumeric(whichbar) + type=1; %we are updating, given a handle + f=whichbar; + else + error(['Input arguments of type ' class(whichbar) ' not valid.']) + end +elseif nargin==1 + f = findobj(allchild(0),'flat','Tag','TMWWaitbar'); + + if isempty(f) + type=2; + name='Waitbar'; + else + type=1; + f=f(1); + end +else + error('Input arguments not valid.'); +end + +x = max(0,min(100*x,100)); + +switch type + case 1, % waitbar(x) update + p = findobj(f,'Type','patch'); + l = findobj(f,'Type','line'); + if isempty(f) | isempty(p) | isempty(l), + error('Couldn''t find waitbar handles.'); + end + xpatch = get(p,'XData'); + xpatch = [0 x x 0]; + set(p,'XData',xpatch) + xline = get(l,'XData'); + set(l,'XData',xline); + + if nargin>2, + % Update waitbar title: + hAxes = findobj(f,'type','axes'); + hTitle = get(hAxes,'title'); + set(hTitle,'string',varargin{1}); + end + + case 2, % waitbar(x,name) initialize + vertMargin = 0; + if nargin > 2, + % we have optional arguments: property-value pairs + if rem (nargin, 2 ) ~= 0 + error( 'Optional initialization arguments must be passed in pairs' ); + end + end + + oldRootUnits = get(0,'Units'); + + set(0, 'Units', 'points'); + screenSize = get(0,'ScreenSize'); + + axFontSize=get(0,'FactoryAxesFontSize'); + + pointsPerPixel = 72/get(0,'ScreenPixelsPerInch'); + + width = 360 * pointsPerPixel; + height = 75 * pointsPerPixel; + pos = [screenSize(3)/2-width/2 screenSize(4)/2-height/2 width height]; + +%pos= [501.75 589.5 393.75 52.5]; + f = figure(... + 'Units', 'points', ... + 'BusyAction', 'queue', ... + 'Position', pos, ... + 'Resize','on', ... + 'CreateFcn','', ... + 'NumberTitle','off', ... + 'IntegerHandle','off', ... + 'MenuBar', 'none', ... + 'Tag','TMWWaitbar',... + 'Interruptible', 'off', ... + 'Visible','on'); + + %%%%%%%%%%%%%%%%%%%%% + % set figure properties as passed to the fcn + % pay special attention to the 'cancel' request + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + if nargin > 2, + propList = varargin(1:2:end); + valueList = varargin(2:2:end); + cancelBtnCreated = 0; + for ii = 1:length( propList ) + try + if strcmp(lower(propList{ii}), 'createcancelbtn' ) & ~cancelBtnCreated + cancelBtnHeight = 23 * pointsPerPixel; + cancelBtnWidth = 60 * pointsPerPixel; + newPos = pos; + vertMargin = vertMargin + cancelBtnHeight; + newPos(4) = newPos(4)+vertMargin; + callbackFcn = [valueList{ii}]; + set( f, 'Position', newPos, 'CloseRequestFcn', callbackFcn ); + cancelButt = uicontrol('Parent',f, ... + 'Units','points', ... + 'Callback',callbackFcn, ... + 'ButtonDownFcn', callbackFcn, ... + 'Enable','on', ... + 'Interruptible','off', ... + 'Position', [pos(3)-cancelBtnWidth*1.4, 7, ... + cancelBtnWidth, cancelBtnHeight], ... + 'String','Cancel', ... + 'Tag','TMWWaitbarCancelButton'); + cancelBtnCreated = 1; + else + % simply set the prop/value pair of the figure + set( f, propList{ii}, valueList{ii}); + end + catch + disp ( ['Warning: could not set property ''' propList{ii} ''' with value ''' num2str(valueList{ii}) '''' ] ); + end + end + end + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + + colormap([]); + + axNorm=[.05 .3 .9 .2]; + % axNorm=[1 1 1 1]; + axPos=axNorm.*[pos(3:4),pos(3:4)] + [0 vertMargin 0 0]; + + h = axes('XLim',[0 100],... + 'YLim',[0 1],... + 'Box','on', ... + 'Units','Points',... + 'FontSize', axFontSize,... + 'Position',axPos,... + 'XTickMode','manual',... + 'YTickMode','manual',... + 'XTick',[],... + 'YTick',[],... + 'XTickLabelMode','manual',... + 'XTickLabel',[],... + 'YTickLabelMode','manual',... + 'YTickLabel',[]); + + tHandle=title(name); + tHandle=get(h,'title'); + oldTitleUnits=get(tHandle,'Units'); + set(tHandle,... + 'Units', 'points',... + 'String', name); + + tExtent=get(tHandle,'Extent'); + set(tHandle,'Units',oldTitleUnits); + + titleHeight=tExtent(4)+axPos(2)+axPos(4)+5; + if titleHeight>pos(4) + pos(4)=titleHeight; + pos(2)=screenSize(4)/2-pos(4)/2; + figPosDirty=logical(1); + else + figPosDirty=logical(0); + end + + if tExtent(3)>pos(3)*1.10; + pos(3)=min(tExtent(3)*1.10,screenSize(3)); + pos(1)=screenSize(3)/2-pos(3)/2; + + axPos([1,3])=axNorm([1,3])*pos(3); + set(h,'Position',axPos); + + figPosDirty=logical(1); + end + + if figPosDirty + set(f,'Position',pos); + end + + xpatch = [0 x x 0]; + ypatch = [0 0 1 1]; + xline = [100 0 0 100 100]; + yline = [0 0 1 1 0]; + + p = patch(xpatch,ypatch,'r','EdgeColor','r','EraseMode','none'); + l = line(xline,yline,'EraseMode','none'); + set(l,'Color',get(gca,'XColor')); + + + set(f,'HandleVisibility','callback','visible','on', 'resize','off'); + + set(0, 'Units', oldRootUnits); +end % case +drawnow; + +if nargout==1, + fout = f; +end diff --git a/src/modules/audio_processing/aecm/main/source/Android.mk b/src/modules/audio_processing/aecm/main/source/Android.mk new file mode 100644 index 0000000000..7ed9f3616a --- /dev/null +++ b/src/modules/audio_processing/aecm/main/source/Android.mk @@ -0,0 +1,55 @@ +# 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. + +LOCAL_PATH := $(call my-dir) + +include $(CLEAR_VARS) + +LOCAL_ARM_MODE := arm +LOCAL_MODULE_CLASS := STATIC_LIBRARIES +LOCAL_MODULE := libwebrtc_aecm +LOCAL_MODULE_TAGS := optional +LOCAL_GENERATED_SOURCES := +LOCAL_SRC_FILES := echo_control_mobile.c \ + aecm_core.c + +# Flags passed to both C and C++ files. +MY_CFLAGS := +MY_CFLAGS_C := +MY_DEFS := '-DNO_TCMALLOC' \ + '-DNO_HEAPCHECKER' \ + '-DWEBRTC_TARGET_PC' \ + '-DWEBRTC_LINUX' \ + '-DWEBRTC_THREAD_RR' +ifeq ($(TARGET_ARCH),arm) +MY_DEFS += \ + '-DWEBRTC_ANDROID' \ + '-DANDROID' +endif +LOCAL_CFLAGS := $(MY_CFLAGS_C) $(MY_CFLAGS) $(MY_DEFS) + +# Include paths placed before CFLAGS/CPPFLAGS +LOCAL_C_INCLUDES := $(LOCAL_PATH)/../../../../.. \ + $(LOCAL_PATH)/../interface \ + $(LOCAL_PATH)/../../../utility \ + $(LOCAL_PATH)/../../../../../common_audio/signal_processing_library/main/interface + +# Flags passed to only C++ (and not C) files. +LOCAL_CPPFLAGS := + +LOCAL_LDFLAGS := + +LOCAL_STATIC_LIBRARIES := + +LOCAL_SHARED_LIBRARIES := libcutils \ + libdl \ + libstlport +LOCAL_ADDITIONAL_DEPENDENCIES := + +include external/stlport/libstlport.mk +include $(BUILD_STATIC_LIBRARY) diff --git a/src/modules/audio_processing/aecm/main/source/aecm.gyp b/src/modules/audio_processing/aecm/main/source/aecm.gyp new file mode 100644 index 0000000000..a535d2b294 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/source/aecm.gyp @@ -0,0 +1,43 @@ +# 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. + +{ + 'includes': [ + '../../../../../common_settings.gypi', + ], + 'targets': [ + { + 'target_name': 'aecm', + 'type': '<(library)', + 'dependencies': [ + '../../../../../common_audio/signal_processing_library/main/source/spl.gyp:spl', + '../../../utility/util.gyp:apm_util' + ], + 'include_dirs': [ + '../interface', + ], + 'direct_dependent_settings': { + 'include_dirs': [ + '../interface', + ], + }, + 'sources': [ + '../interface/echo_control_mobile.h', + 'echo_control_mobile.c', + 'aecm_core.c', + 'aecm_core.h', + ], + }, + ], +} + +# Local Variables: +# tab-width:2 +# indent-tabs-mode:nil +# End: +# vim: set expandtab tabstop=2 shiftwidth=2: diff --git a/src/modules/audio_processing/aecm/main/source/aecm_core.c b/src/modules/audio_processing/aecm/main/source/aecm_core.c new file mode 100644 index 0000000000..f17f1bf237 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/source/aecm_core.c @@ -0,0 +1,2534 @@ +/* + * 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 <stdlib.h> + +#include "aecm_core.h" +#include "ring_buffer.h" +#include "echo_control_mobile.h" +#include "typedefs.h" + +// TODO(bjornv): Will be removed in final version. +//#include <stdio.h> + +#ifdef ARM_WINM_LOG +#include <stdio.h> +#include <windows.h> +#endif + +// BANDLAST - BANDFIRST must be < 32 +#define BANDFIRST 12 // Only bit BANDFIRST through bit BANDLAST are processed +#define BANDLAST 43 + +#ifdef ARM_WINM +#define WebRtcSpl_AddSatW32(a,b) _AddSatInt(a,b) +#define WebRtcSpl_SubSatW32(a,b) _SubSatInt(a,b) +#endif +// 16 instructions on most risc machines for 32-bit bitcount ! + +#ifdef AEC_DEBUG +FILE *dfile; +FILE *testfile; +#endif + +#ifdef AECM_SHORT + +// Square root of Hanning window in Q14 +static const WebRtc_Word16 kSqrtHanning[] = +{ + 0, 804, 1606, 2404, 3196, 3981, 4756, 5520, + 6270, 7005, 7723, 8423, 9102, 9760, 10394, 11003, + 11585, 12140, 12665, 13160, 13623, 14053, 14449, 14811, + 15137, 15426, 15679, 15893, 16069, 16207, 16305, 16364, + 16384 +}; + +#else + +// Square root of Hanning window in Q14 +static const WebRtc_Word16 kSqrtHanning[] = {0, 399, 798, 1196, 1594, 1990, 2386, 2780, 3172, + 3562, 3951, 4337, 4720, 5101, 5478, 5853, 6224, 6591, 6954, 7313, 7668, 8019, 8364, + 8705, 9040, 9370, 9695, 10013, 10326, 10633, 10933, 11227, 11514, 11795, 12068, 12335, + 12594, 12845, 13089, 13325, 13553, 13773, 13985, 14189, 14384, 14571, 14749, 14918, + 15079, 15231, 15373, 15506, 15631, 15746, 15851, 15947, 16034, 16111, 16179, 16237, + 16286, 16325, 16354, 16373, 16384}; + +#endif + +//Q15 alpha = 0.99439986968132 const Factor for magnitude approximation +static const WebRtc_UWord16 kAlpha1 = 32584; +//Q15 beta = 0.12967166976970 const Factor for magnitude approximation +static const WebRtc_UWord16 kBeta1 = 4249; +//Q15 alpha = 0.94234827210087 const Factor for magnitude approximation +static const WebRtc_UWord16 kAlpha2 = 30879; +//Q15 beta = 0.33787806009150 const Factor for magnitude approximation +static const WebRtc_UWord16 kBeta2 = 11072; +//Q15 alpha = 0.82247698684306 const Factor for magnitude approximation +static const WebRtc_UWord16 kAlpha3 = 26951; +//Q15 beta = 0.57762063060713 const Factor for magnitude approximation +static const WebRtc_UWord16 kBeta3 = 18927; + +// Initialization table for echo channel in 8 kHz +static const WebRtc_Word16 kChannelStored8kHz[PART_LEN1] = { + 2040, 1815, 1590, 1498, 1405, 1395, 1385, 1418, + 1451, 1506, 1562, 1644, 1726, 1804, 1882, 1918, + 1953, 1982, 2010, 2025, 2040, 2034, 2027, 2021, + 2014, 1997, 1980, 1925, 1869, 1800, 1732, 1683, + 1635, 1604, 1572, 1545, 1517, 1481, 1444, 1405, + 1367, 1331, 1294, 1270, 1245, 1239, 1233, 1247, + 1260, 1282, 1303, 1338, 1373, 1407, 1441, 1470, + 1499, 1524, 1549, 1565, 1582, 1601, 1621, 1649, + 1676 +}; + +// Initialization table for echo channel in 16 kHz +static const WebRtc_Word16 kChannelStored16kHz[PART_LEN1] = { + 2040, 1590, 1405, 1385, 1451, 1562, 1726, 1882, + 1953, 2010, 2040, 2027, 2014, 1980, 1869, 1732, + 1635, 1572, 1517, 1444, 1367, 1294, 1245, 1233, + 1260, 1303, 1373, 1441, 1499, 1549, 1582, 1621, + 1676, 1741, 1802, 1861, 1921, 1983, 2040, 2102, + 2170, 2265, 2375, 2515, 2651, 2781, 2922, 3075, + 3253, 3471, 3738, 3976, 4151, 4258, 4308, 4288, + 4270, 4253, 4237, 4179, 4086, 3947, 3757, 3484, + 3153 +}; + +#ifdef ARM_WINM_LOG +HANDLE logFile = NULL; +#endif + +static void WebRtcAecm_ComfortNoise(AecmCore_t* const aecm, const WebRtc_UWord16 * const dfa, + WebRtc_Word16 * const outReal, + WebRtc_Word16 * const outImag, + const WebRtc_Word16 * const lambda); + +static __inline WebRtc_UWord32 WebRtcAecm_SetBit(WebRtc_UWord32 in, WebRtc_Word32 pos) +{ + WebRtc_UWord32 mask, out; + + mask = WEBRTC_SPL_SHIFT_W32(1, pos); + out = (in | mask); + + return out; +} + +// WebRtcAecm_Hisser(...) +// +// This function compares the binary vector specvec with all rows of the binary matrix specmat +// and counts per row the number of times they have the same value. +// Input: +// - specvec : binary "vector" that is stored in a long +// - specmat : binary "matrix" that is stored as a vector of long +// Output: +// - bcount : "Vector" stored as a long, containing for each row the number of times +// the matrix row and the input vector have the same value +// +// +void WebRtcAecm_Hisser(const WebRtc_UWord32 specvec, const WebRtc_UWord32 * const specmat, + WebRtc_UWord32 * const bcount) +{ + int n; + WebRtc_UWord32 a, b; + register WebRtc_UWord32 tmp; + + a = specvec; + // compare binary vector specvec with all rows of the binary matrix specmat + for (n = 0; n < MAX_DELAY; n++) + { + b = specmat[n]; + a = (specvec ^ b); + // 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; + + bcount[n] = tmp; + } +} + +// WebRtcAecm_BSpectrum(...) +// +// Computes the binary spectrum by comparing the input spectrum with a threshold spectrum. +// +// Input: +// - spectrum : Spectrum of which the binary spectrum should be calculated. +// - thresvec : Threshold spectrum with which the input spectrum is compared. +// Return: +// - out : Binary spectrum +// +WebRtc_UWord32 WebRtcAecm_BSpectrum(const WebRtc_UWord16 * const spectrum, + const WebRtc_UWord16 * const thresvec) +{ + int k; + WebRtc_UWord32 out; + + out = 0; + for (k = BANDFIRST; k <= BANDLAST; k++) + { + if (spectrum[k] > thresvec[k]) + { + out = WebRtcAecm_SetBit(out, k - BANDFIRST); + } + } + + return out; +} + +// WebRtcAecm_MedianEstimator(...) +// +// Calculates the median recursively. +// +// Input: +// - newVal : new additional value +// - medianVec : vector with current medians +// - factor : factor for smoothing +// +// Output: +// - medianVec : vector with updated median +// +int WebRtcAecm_MedianEstimator(const WebRtc_UWord16 newVal, WebRtc_UWord16 * const medianVec, + const int factor) +{ + WebRtc_Word32 median; + WebRtc_Word32 diff; + + median = (WebRtc_Word32)medianVec[0]; + + //median = median + ((newVal-median)>>factor); + diff = (WebRtc_Word32)newVal - median; + diff = WEBRTC_SPL_SHIFT_W32(diff, -factor); + median = median + diff; + + medianVec[0] = (WebRtc_UWord16)median; + + return 0; +} + +int WebRtcAecm_CreateCore(AecmCore_t **aecmInst) +{ + AecmCore_t *aecm = malloc(sizeof(AecmCore_t)); + *aecmInst = aecm; + if (aecm == NULL) + { + return -1; + } + + if (WebRtcApm_CreateBuffer(&aecm->farFrameBuf, FRAME_LEN + PART_LEN) == -1) + { + WebRtcAecm_FreeCore(aecm); + aecm = NULL; + return -1; + } + + if (WebRtcApm_CreateBuffer(&aecm->nearNoisyFrameBuf, FRAME_LEN + PART_LEN) == -1) + { + WebRtcAecm_FreeCore(aecm); + aecm = NULL; + return -1; + } + + if (WebRtcApm_CreateBuffer(&aecm->nearCleanFrameBuf, FRAME_LEN + PART_LEN) == -1) + { + WebRtcAecm_FreeCore(aecm); + aecm = NULL; + return -1; + } + + if (WebRtcApm_CreateBuffer(&aecm->outFrameBuf, FRAME_LEN + PART_LEN) == -1) + { + WebRtcAecm_FreeCore(aecm); + aecm = NULL; + return -1; + } + + return 0; +} + +// WebRtcAecm_InitCore(...) +// +// This function initializes the AECM instant created with WebRtcAecm_CreateCore(...) +// Input: +// - aecm : Pointer to the Echo Suppression instance +// - samplingFreq : Sampling Frequency +// +// Output: +// - aecm : Initialized instance +// +// Return value : 0 - Ok +// -1 - Error +// +int WebRtcAecm_InitCore(AecmCore_t * const aecm, int samplingFreq) +{ + int retVal = 0; + WebRtc_Word16 i; + WebRtc_Word16 tmp16; + + if (samplingFreq != 8000 && samplingFreq != 16000) + { + samplingFreq = 8000; + retVal = -1; + } + // sanity check of sampling frequency + aecm->mult = (WebRtc_Word16)samplingFreq / 8000; + + aecm->farBufWritePos = 0; + aecm->farBufReadPos = 0; + aecm->knownDelay = 0; + aecm->lastKnownDelay = 0; + + WebRtcApm_InitBuffer(aecm->farFrameBuf); + WebRtcApm_InitBuffer(aecm->nearNoisyFrameBuf); + WebRtcApm_InitBuffer(aecm->nearCleanFrameBuf); + WebRtcApm_InitBuffer(aecm->outFrameBuf); + + memset(aecm->xBuf, 0, sizeof(aecm->xBuf)); + memset(aecm->dBufClean, 0, sizeof(aecm->dBufClean)); + memset(aecm->dBufNoisy, 0, sizeof(aecm->dBufNoisy)); + memset(aecm->outBuf, 0, sizeof(WebRtc_Word16) * PART_LEN); + + aecm->seed = 666; + aecm->totCount = 0; + + memset(aecm->xfaHistory, 0, sizeof(WebRtc_UWord16) * (PART_LEN1) * MAX_DELAY); + + aecm->delHistoryPos = MAX_DELAY; + + memset(aecm->medianYlogspec, 0, sizeof(WebRtc_UWord16) * PART_LEN1); + memset(aecm->medianXlogspec, 0, sizeof(WebRtc_UWord16) * PART_LEN1); + memset(aecm->medianBCount, 0, sizeof(WebRtc_UWord16) * MAX_DELAY); + memset(aecm->bxHistory, 0, sizeof(aecm->bxHistory)); + + // Initialize to reasonable values + aecm->currentDelay = 8; + aecm->previousDelay = 8; + aecm->delayAdjust = 0; + + aecm->nlpFlag = 1; + aecm->fixedDelay = -1; + + memset(aecm->xfaQDomainBuf, 0, sizeof(WebRtc_Word16) * MAX_DELAY); + aecm->dfaCleanQDomain = 0; + aecm->dfaCleanQDomainOld = 0; + aecm->dfaNoisyQDomain = 0; + aecm->dfaNoisyQDomainOld = 0; + + memset(aecm->nearLogEnergy, 0, sizeof(WebRtc_Word16) * MAX_BUF_LEN); + memset(aecm->farLogEnergy, 0, sizeof(WebRtc_Word16) * MAX_BUF_LEN); + memset(aecm->echoAdaptLogEnergy, 0, sizeof(WebRtc_Word16) * MAX_BUF_LEN); + memset(aecm->echoStoredLogEnergy, 0, sizeof(WebRtc_Word16) * MAX_BUF_LEN); + + // Initialize the echo channels with a stored shape. + if (samplingFreq == 8000) + { + memcpy(aecm->channelAdapt16, kChannelStored8kHz, sizeof(WebRtc_Word16) * PART_LEN1); + } + else + { + memcpy(aecm->channelAdapt16, kChannelStored16kHz, sizeof(WebRtc_Word16) * PART_LEN1); + } + memcpy(aecm->channelStored, aecm->channelAdapt16, sizeof(WebRtc_Word16) * PART_LEN1); + for (i = 0; i < PART_LEN1; i++) + { + aecm->channelAdapt32[i] = WEBRTC_SPL_LSHIFT_W32( + (WebRtc_Word32)(aecm->channelAdapt16[i]), 16); + } + + memset(aecm->echoFilt, 0, sizeof(WebRtc_Word32) * PART_LEN1); + memset(aecm->nearFilt, 0, sizeof(WebRtc_Word16) * PART_LEN1); + aecm->noiseEstCtr = 0; + + aecm->cngMode = AecmTrue; + + // Increase the noise Q domain with increasing frequency, to correspond to the + // expected energy levels. + // Also shape the initial noise level with this consideration. +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + for (i = 0; i < PART_LEN1; i++) + { + if (i < PART_LEN1 >> 2) + { + aecm->noiseEstQDomain[i] = 10; + tmp16 = PART_LEN1 - i; + aecm->noiseEst[i] = (tmp16 * tmp16) << 4; + } else if (i < PART_LEN1 >> 1) + { + aecm->noiseEstQDomain[i] = 11; + tmp16 = PART_LEN1 - i; + aecm->noiseEst[i] = ((tmp16 * tmp16) << 4) << 1; + } else + { + aecm->noiseEstQDomain[i] = 12; + aecm->noiseEst[i] = aecm->noiseEst[(PART_LEN1 >> 1) - 1] << 1; + } + } +#else + for (i = 0; i < PART_LEN1 >> 2; i++) + { + aecm->noiseEstQDomain[i] = 10; + tmp16 = PART_LEN1 - i; + aecm->noiseEst[i] = (tmp16 * tmp16) << 4; + } + for (; i < PART_LEN1 >> 1; i++) + { + aecm->noiseEstQDomain[i] = 11; + tmp16 = PART_LEN1 - i; + aecm->noiseEst[i] = ((tmp16 * tmp16) << 4) << 1; + } + for (; i < PART_LEN1; i++) + { + aecm->noiseEstQDomain[i] = 12; + aecm->noiseEst[i] = aecm->noiseEst[(PART_LEN1 >> 1) - 1] << 1; + } +#endif + + aecm->mseAdaptOld = 1000; + aecm->mseStoredOld = 1000; + aecm->mseThreshold = WEBRTC_SPL_WORD32_MAX; + + aecm->farEnergyMin = WEBRTC_SPL_WORD16_MAX; + aecm->farEnergyMax = WEBRTC_SPL_WORD16_MIN; + aecm->farEnergyMaxMin = 0; + aecm->farEnergyVAD = FAR_ENERGY_MIN; // This prevents false speech detection at the + // beginning. + aecm->farEnergyMSE = 0; + aecm->currentVADValue = 0; + aecm->vadUpdateCount = 0; + aecm->firstVAD = 1; + + aecm->delayCount = 0; + aecm->newDelayCorrData = 0; + aecm->lastDelayUpdateCount = 0; + memset(aecm->delayCorrelation, 0, sizeof(WebRtc_Word16) * ((CORR_MAX << 1) + 1)); + + aecm->startupState = 0; + aecm->mseChannelCount = 0; + aecm->supGain = SUPGAIN_DEFAULT; + aecm->supGainOld = SUPGAIN_DEFAULT; + aecm->delayOffsetFlag = 0; + + memset(aecm->delayHistogram, 0, sizeof(aecm->delayHistogram)); + aecm->delayVadCount = 0; + aecm->maxDelayHistIdx = 0; + aecm->lastMinPos = 0; + + aecm->supGainErrParamA = SUPGAIN_ERROR_PARAM_A; + aecm->supGainErrParamD = SUPGAIN_ERROR_PARAM_D; + aecm->supGainErrParamDiffAB = SUPGAIN_ERROR_PARAM_A - SUPGAIN_ERROR_PARAM_B; + aecm->supGainErrParamDiffBD = SUPGAIN_ERROR_PARAM_B - SUPGAIN_ERROR_PARAM_D; + + return 0; +} + +int WebRtcAecm_Control(AecmCore_t *aecm, int delay, int nlpFlag, int delayOffsetFlag) +{ + aecm->nlpFlag = nlpFlag; + aecm->fixedDelay = delay; + aecm->delayOffsetFlag = delayOffsetFlag; + + return 0; +} + +// WebRtcAecm_GetNewDelPos(...) +// +// Moves the pointer to the next entry. Returns to zero if max position reached. +// +// Input: +// - aecm : Pointer to the AECM instance +// Return: +// - pos : New position in the history. +// +// +WebRtc_Word16 WebRtcAecm_GetNewDelPos(AecmCore_t * const aecm) +{ + WebRtc_Word16 pos; + + pos = aecm->delHistoryPos; + pos++; + if (pos >= MAX_DELAY) + { + pos = 0; + } + aecm->delHistoryPos = pos; + + return pos; +} + +// WebRtcAecm_EstimateDelay(...) +// +// Estimate the delay of the echo signal. +// +// Inputs: +// - aecm : Pointer to the AECM instance +// - farSpec : Delayed farend magnitude spectrum +// - nearSpec : Nearend magnitude spectrum +// - stages : Q-domain of xxFIX and yyFIX (without dynamic Q-domain) +// - xfaQ : normalization factor, i.e., Q-domain before FFT +// Return: +// - delay : Estimated delay +// +WebRtc_Word16 WebRtcAecm_EstimateDelay(AecmCore_t * const aecm, + const WebRtc_UWord16 * const farSpec, + const WebRtc_UWord16 * const nearSpec, + const WebRtc_Word16 xfaQ) +{ + WebRtc_UWord32 bxspectrum, byspectrum; + WebRtc_UWord32 bcount[MAX_DELAY]; + + int i, res; + + WebRtc_UWord16 xmean[PART_LEN1], ymean[PART_LEN1]; + WebRtc_UWord16 dtmp1; + WebRtc_Word16 fcount[MAX_DELAY]; + + //WebRtc_Word16 res; + WebRtc_Word16 histpos; + WebRtc_Word16 maxHistLvl; + WebRtc_UWord16 *state; + WebRtc_Word16 minpos = -1; + + enum + { + kVadCountThreshold = 25 + }; + enum + { + kMaxHistogram = 600 + }; + + histpos = WebRtcAecm_GetNewDelPos(aecm); + + for (i = 0; i < PART_LEN1; i++) + { + aecm->xfaHistory[i][histpos] = farSpec[i]; + + state = &(aecm->medianXlogspec[i]); + res = WebRtcAecm_MedianEstimator(farSpec[i], state, 6); + + state = &(aecm->medianYlogspec[i]); + res = WebRtcAecm_MedianEstimator(nearSpec[i], state, 6); + + // Mean: + // FLOAT: + // ymean = dtmp2/MAX_DELAY + // + // FIX: + // input: dtmp2FIX in Q0 + // output: ymeanFIX in Q8 + // 20 = 1/MAX_DELAY in Q13 = 1/MAX_DELAY * 2^13 + xmean[i] = (aecm->medianXlogspec[i]); + ymean[i] = (aecm->medianYlogspec[i]); + + } + // Update Q-domain buffer + aecm->xfaQDomainBuf[histpos] = xfaQ; + + // Get binary spectra + // FLOAT: + // bxspectrum = bspectrum(xlogspec, xmean); + // + // FIX: + // input: xlogspecFIX,ylogspecFIX in Q8 + // xmeanFIX, ymeanFIX in Q8 + // output: unsigned long bxspectrum, byspectrum in Q0 + bxspectrum = WebRtcAecm_BSpectrum(farSpec, xmean); + byspectrum = WebRtcAecm_BSpectrum(nearSpec, ymean); + + // Shift binary spectrum history + memmove(&(aecm->bxHistory[1]), &(aecm->bxHistory[0]), + (MAX_DELAY - 1) * sizeof(WebRtc_UWord32)); + + aecm->bxHistory[0] = bxspectrum; + + // Compare with delayed spectra + WebRtcAecm_Hisser(byspectrum, aecm->bxHistory, bcount); + + for (i = 0; i < MAX_DELAY; i++) + { + // Update sum + // bcount is constrained to [0, 32], meaning we can smooth with a factor up to 2^11. + dtmp1 = (WebRtc_UWord16)bcount[i]; + dtmp1 = WEBRTC_SPL_LSHIFT_W16(dtmp1, 9); + state = &(aecm->medianBCount[i]); + res = WebRtcAecm_MedianEstimator(dtmp1, state, 9); + fcount[i] = (aecm->medianBCount[i]); + } + + // Find minimum + minpos = WebRtcSpl_MinIndexW16(fcount, MAX_DELAY); + + // 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 (aecm->currentVADValue == 1) + { + if (aecm->delayVadCount >= kVadCountThreshold) + { + // Increment the histogram at the current minimum position. + if (aecm->delayHistogram[minpos] < kMaxHistogram) + { + aecm->delayHistogram[minpos] += 3; + } + +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + // Decrement the entire histogram. + for (i = 0; i < MAX_DELAY; i++) + { + if (aecm->delayHistogram[i] > 0) + { + aecm->delayHistogram[i]--; + } + } + + // Select the histogram index corresponding to the maximum bin as the delay. + maxHistLvl = 0; + aecm->maxDelayHistIdx = 0; + for (i = 0; i < MAX_DELAY; i++) + { + if (aecm->delayHistogram[i] > maxHistLvl) + { + maxHistLvl = aecm->delayHistogram[i]; + aecm->maxDelayHistIdx = i; + } + } +#else + maxHistLvl = 0; + aecm->maxDelayHistIdx = 0; + + for (i = 0; i < MAX_DELAY; i++) + { + WebRtc_Word16 tempVar = aecm->delayHistogram[i]; + + // Decrement the entire histogram. + if (tempVar > 0) + { + tempVar--; + aecm->delayHistogram[i] = tempVar; + + // Select the histogram index corresponding to the maximum bin as the delay. + if (tempVar > maxHistLvl) + { + maxHistLvl = tempVar; + aecm->maxDelayHistIdx = i; + } + } + } +#endif + } else + { + aecm->delayVadCount++; + } + } else + { + aecm->delayVadCount = 0; + } + + return aecm->maxDelayHistIdx; +} + +int WebRtcAecm_FreeCore(AecmCore_t *aecm) +{ + if (aecm == NULL) + { + return -1; + } + + WebRtcApm_FreeBuffer(aecm->farFrameBuf); + WebRtcApm_FreeBuffer(aecm->nearNoisyFrameBuf); + WebRtcApm_FreeBuffer(aecm->nearCleanFrameBuf); + WebRtcApm_FreeBuffer(aecm->outFrameBuf); + + free(aecm); + + return 0; +} + +void WebRtcAecm_ProcessFrame(AecmCore_t * const aecm, const WebRtc_Word16 * const farend, + const WebRtc_Word16 * const nearendNoisy, + const WebRtc_Word16 * const nearendClean, + WebRtc_Word16 * const out) +{ + WebRtc_Word16 farBlock[PART_LEN]; + WebRtc_Word16 nearNoisyBlock[PART_LEN]; + WebRtc_Word16 nearCleanBlock[PART_LEN]; + WebRtc_Word16 outBlock[PART_LEN]; + WebRtc_Word16 farFrame[FRAME_LEN]; + int size = 0; + + // Buffer the current frame. + // Fetch an older one corresponding to the delay. + WebRtcAecm_BufferFarFrame(aecm, farend, FRAME_LEN); + WebRtcAecm_FetchFarFrame(aecm, farFrame, FRAME_LEN, aecm->knownDelay); + + // Buffer the synchronized far and near frames, + // to pass the smaller blocks individually. + WebRtcApm_WriteBuffer(aecm->farFrameBuf, farFrame, FRAME_LEN); + WebRtcApm_WriteBuffer(aecm->nearNoisyFrameBuf, nearendNoisy, FRAME_LEN); + if (nearendClean != NULL) + { + WebRtcApm_WriteBuffer(aecm->nearCleanFrameBuf, nearendClean, FRAME_LEN); + } + + // Process as many blocks as possible. + while (WebRtcApm_get_buffer_size(aecm->farFrameBuf) >= PART_LEN) + { + WebRtcApm_ReadBuffer(aecm->farFrameBuf, farBlock, PART_LEN); + WebRtcApm_ReadBuffer(aecm->nearNoisyFrameBuf, nearNoisyBlock, PART_LEN); + if (nearendClean != NULL) + { + WebRtcApm_ReadBuffer(aecm->nearCleanFrameBuf, nearCleanBlock, PART_LEN); + WebRtcAecm_ProcessBlock(aecm, farBlock, nearNoisyBlock, nearCleanBlock, outBlock); + } else + { + WebRtcAecm_ProcessBlock(aecm, farBlock, nearNoisyBlock, NULL, outBlock); + } + + WebRtcApm_WriteBuffer(aecm->outFrameBuf, outBlock, PART_LEN); + } + + // Stuff the out buffer if we have less than a frame to output. + // This should only happen for the first frame. + size = WebRtcApm_get_buffer_size(aecm->outFrameBuf); + if (size < FRAME_LEN) + { + WebRtcApm_StuffBuffer(aecm->outFrameBuf, FRAME_LEN - size); + } + + // Obtain an output frame. + WebRtcApm_ReadBuffer(aecm->outFrameBuf, out, FRAME_LEN); +} + +// WebRtcAecm_AsymFilt(...) +// +// Performs asymmetric filtering. +// +// Inputs: +// - filtOld : Previous filtered value. +// - inVal : New input value. +// - stepSizePos : Step size when we have a positive contribution. +// - stepSizeNeg : Step size when we have a negative contribution. +// +// Output: +// +// Return: - Filtered value. +// +WebRtc_Word16 WebRtcAecm_AsymFilt(const WebRtc_Word16 filtOld, const WebRtc_Word16 inVal, + const WebRtc_Word16 stepSizePos, + const WebRtc_Word16 stepSizeNeg) +{ + WebRtc_Word16 retVal; + + if ((filtOld == WEBRTC_SPL_WORD16_MAX) | (filtOld == WEBRTC_SPL_WORD16_MIN)) + { + return inVal; + } + retVal = filtOld; + if (filtOld > inVal) + { + retVal -= WEBRTC_SPL_RSHIFT_W16(filtOld - inVal, stepSizeNeg); + } else + { + retVal += WEBRTC_SPL_RSHIFT_W16(inVal - filtOld, stepSizePos); + } + + return retVal; +} + +// WebRtcAecm_CalcEnergies(...) +// +// This function calculates the log of energies for nearend, farend and estimated +// echoes. There is also an update of energy decision levels, i.e. internl VAD. +// +// +// @param aecm [i/o] Handle of the AECM instance. +// @param delayDiff [in] Delay position in farend buffer. +// @param nearEner [in] Near end energy for current block (Q[aecm->dfaQDomain]). +// @param echoEst [i/o] Estimated echo +// (Q[aecm->xfaQDomain[delayDiff]+RESOLUTION_CHANNEL16]). +// +void WebRtcAecm_CalcEnergies(AecmCore_t * const aecm, const WebRtc_Word16 delayDiff, + const WebRtc_UWord32 nearEner, WebRtc_Word32 * const echoEst) +{ + // Local variables + WebRtc_UWord32 tmpAdapt, tmpStored, tmpFar; + + int i; + + WebRtc_Word16 zeros, frac; + WebRtc_Word16 tmp16; + WebRtc_Word16 increase_max_shifts = 4; + WebRtc_Word16 decrease_max_shifts = 11; + WebRtc_Word16 increase_min_shifts = 11; + WebRtc_Word16 decrease_min_shifts = 3; + + // Get log of near end energy and store in buffer + + // Shift buffer + memmove(aecm->nearLogEnergy + 1, aecm->nearLogEnergy, + sizeof(WebRtc_Word16) * (MAX_BUF_LEN - 1)); + + // Logarithm of integrated magnitude spectrum (nearEner) + if (nearEner) + { + zeros = WebRtcSpl_NormU32(nearEner); + frac = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_U32( + (WEBRTC_SPL_LSHIFT_U32(nearEner, zeros) & 0x7FFFFFFF), + 23); + // log2 in Q8 + aecm->nearLogEnergy[0] = WEBRTC_SPL_LSHIFT_W16((31 - zeros), 8) + frac; + aecm->nearLogEnergy[0] -= WEBRTC_SPL_LSHIFT_W16(aecm->dfaNoisyQDomain, 8); + } else + { + aecm->nearLogEnergy[0] = 0; + } + aecm->nearLogEnergy[0] += WEBRTC_SPL_LSHIFT_W16(PART_LEN_SHIFT, 7); + // END: Get log of near end energy + + // Get energy for the delayed far end signal and estimated + // echo using both stored and adapted channels. + tmpAdapt = 0; + tmpStored = 0; + tmpFar = 0; + + for (i = 0; i < PART_LEN1; i++) + { + // Get estimated echo energies for adaptive channel and stored channel + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], + aecm->xfaHistory[i][delayDiff]); + tmpFar += (WebRtc_UWord32)(aecm->xfaHistory[i][delayDiff]); + tmpAdapt += WEBRTC_SPL_UMUL_16_16(aecm->channelAdapt16[i], + aecm->xfaHistory[i][delayDiff]); + tmpStored += (WebRtc_UWord32)echoEst[i]; + } + // Shift buffers + memmove(aecm->farLogEnergy + 1, aecm->farLogEnergy, + sizeof(WebRtc_Word16) * (MAX_BUF_LEN - 1)); + memmove(aecm->echoAdaptLogEnergy + 1, aecm->echoAdaptLogEnergy, + sizeof(WebRtc_Word16) * (MAX_BUF_LEN - 1)); + memmove(aecm->echoStoredLogEnergy + 1, aecm->echoStoredLogEnergy, + sizeof(WebRtc_Word16) * (MAX_BUF_LEN - 1)); + + // Logarithm of delayed far end energy + if (tmpFar) + { + zeros = WebRtcSpl_NormU32(tmpFar); + frac = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_U32((WEBRTC_SPL_LSHIFT_U32(tmpFar, zeros) + & 0x7FFFFFFF), 23); + // log2 in Q8 + aecm->farLogEnergy[0] = WEBRTC_SPL_LSHIFT_W16((31 - zeros), 8) + frac; + aecm->farLogEnergy[0] -= WEBRTC_SPL_LSHIFT_W16(aecm->xfaQDomainBuf[delayDiff], 8); + } else + { + aecm->farLogEnergy[0] = 0; + } + aecm->farLogEnergy[0] += WEBRTC_SPL_LSHIFT_W16(PART_LEN_SHIFT, 7); + + // Logarithm of estimated echo energy through adapted channel + if (tmpAdapt) + { + zeros = WebRtcSpl_NormU32(tmpAdapt); + frac = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_U32((WEBRTC_SPL_LSHIFT_U32(tmpAdapt, zeros) + & 0x7FFFFFFF), 23); + //log2 in Q8 + aecm->echoAdaptLogEnergy[0] = WEBRTC_SPL_LSHIFT_W16((31 - zeros), 8) + frac; + aecm->echoAdaptLogEnergy[0] + -= WEBRTC_SPL_LSHIFT_W16(RESOLUTION_CHANNEL16 + aecm->xfaQDomainBuf[delayDiff], 8); + } else + { + aecm->echoAdaptLogEnergy[0] = 0; + } + aecm->echoAdaptLogEnergy[0] += WEBRTC_SPL_LSHIFT_W16(PART_LEN_SHIFT, 7); + + // Logarithm of estimated echo energy through stored channel + if (tmpStored) + { + zeros = WebRtcSpl_NormU32(tmpStored); + frac = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_U32((WEBRTC_SPL_LSHIFT_U32(tmpStored, zeros) + & 0x7FFFFFFF), 23); + //log2 in Q8 + aecm->echoStoredLogEnergy[0] = WEBRTC_SPL_LSHIFT_W16((31 - zeros), 8) + frac; + aecm->echoStoredLogEnergy[0] + -= WEBRTC_SPL_LSHIFT_W16(RESOLUTION_CHANNEL16 + aecm->xfaQDomainBuf[delayDiff], 8); + } else + { + aecm->echoStoredLogEnergy[0] = 0; + } + aecm->echoStoredLogEnergy[0] += WEBRTC_SPL_LSHIFT_W16(PART_LEN_SHIFT, 7); + + // Update farend energy levels (min, max, vad, mse) + if (aecm->farLogEnergy[0] > FAR_ENERGY_MIN) + { + if (aecm->startupState == 0) + { + increase_max_shifts = 2; + decrease_min_shifts = 2; + increase_min_shifts = 8; + } + + aecm->farEnergyMin = WebRtcAecm_AsymFilt(aecm->farEnergyMin, aecm->farLogEnergy[0], + increase_min_shifts, decrease_min_shifts); + aecm->farEnergyMax = WebRtcAecm_AsymFilt(aecm->farEnergyMax, aecm->farLogEnergy[0], + increase_max_shifts, decrease_max_shifts); + aecm->farEnergyMaxMin = (aecm->farEnergyMax - aecm->farEnergyMin); + + // Dynamic VAD region size + tmp16 = 2560 - aecm->farEnergyMin; + if (tmp16 > 0) + { + tmp16 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(tmp16, FAR_ENERGY_VAD_REGION, 9); + } else + { + tmp16 = 0; + } + tmp16 += FAR_ENERGY_VAD_REGION; + + if ((aecm->startupState == 0) | (aecm->vadUpdateCount > 1024)) + { + // In startup phase or VAD update halted + aecm->farEnergyVAD = aecm->farEnergyMin + tmp16; + } else + { + if (aecm->farEnergyVAD > aecm->farLogEnergy[0]) + { + aecm->farEnergyVAD += WEBRTC_SPL_RSHIFT_W16(aecm->farLogEnergy[0] + tmp16 + - aecm->farEnergyVAD, 6); + aecm->vadUpdateCount = 0; + } else + { + aecm->vadUpdateCount++; + } + } + // Put MSE threshold higher than VAD + aecm->farEnergyMSE = aecm->farEnergyVAD + (1 << 8); + } + + // Update VAD variables + if (aecm->farLogEnergy[0] > aecm->farEnergyVAD) + { + if ((aecm->startupState == 0) | (aecm->farEnergyMaxMin > FAR_ENERGY_DIFF)) + { + // We are in startup or have significant dynamics in input speech level + aecm->currentVADValue = 1; + } + } else + { + aecm->currentVADValue = 0; + } + if ((aecm->currentVADValue) && (aecm->firstVAD)) + { + aecm->firstVAD = 0; + if (aecm->echoAdaptLogEnergy[0] > aecm->nearLogEnergy[0]) + { + // The estimated echo has higher energy than the near end signal. This means that + // the initialization was too aggressive. Scale down by a factor 8 + for (i = 0; i < PART_LEN1; i++) + { + aecm->channelAdapt16[i] >>= 3; + } + // Compensate the adapted echo energy level accordingly. + aecm->echoAdaptLogEnergy[0] -= (3 << 8); + aecm->firstVAD = 1; + } + } + // END: Energies of delayed far, echo estimates + // TODO(bjornv): Will be removed in final version. +#ifdef VAD_DATA + fwrite(&(aecm->currentVADValue), sizeof(WebRtc_Word16), 1, aecm->vad_file); + fwrite(&(aecm->currentDelay), sizeof(WebRtc_Word16), 1, aecm->delay_file); + fwrite(&(aecm->farLogEnergy[0]), sizeof(WebRtc_Word16), 1, aecm->far_cur_file); + fwrite(&(aecm->farEnergyMin), sizeof(WebRtc_Word16), 1, aecm->far_min_file); + fwrite(&(aecm->farEnergyMax), sizeof(WebRtc_Word16), 1, aecm->far_max_file); + fwrite(&(aecm->farEnergyVAD), sizeof(WebRtc_Word16), 1, aecm->far_vad_file); +#endif +} + +// WebRtcAecm_CalcStepSize(...) +// +// This function calculates the step size used in channel estimation +// +// +// @param aecm [in] Handle of the AECM instance. +// @param mu [out] (Return value) Stepsize in log2(), i.e. number of shifts. +// +// +WebRtc_Word16 WebRtcAecm_CalcStepSize(AecmCore_t * const aecm) +{ + + WebRtc_Word32 tmp32; + WebRtc_Word16 tmp16; + WebRtc_Word16 mu; + + // Here we calculate the step size mu used in the + // following NLMS based Channel estimation algorithm + mu = MU_MAX; + if (!aecm->currentVADValue) + { + // Far end energy level too low, no channel update + mu = 0; + } else if (aecm->startupState > 0) + { + if (aecm->farEnergyMin >= aecm->farEnergyMax) + { + mu = MU_MIN; + } else + { + tmp16 = (aecm->farLogEnergy[0] - aecm->farEnergyMin); + tmp32 = WEBRTC_SPL_MUL_16_16(tmp16, MU_DIFF); + tmp32 = WebRtcSpl_DivW32W16(tmp32, aecm->farEnergyMaxMin); + mu = MU_MIN - 1 - (WebRtc_Word16)(tmp32); + // The -1 is an alternative to rounding. This way we get a larger + // stepsize, so we in some sense compensate for truncation in NLMS + } + if (mu < MU_MAX) + { + mu = MU_MAX; // Equivalent with maximum step size of 2^-MU_MAX + } + } + // END: Update step size + + return mu; +} + +// WebRtcAecm_UpdateChannel(...) +// +// This function performs channel estimation. NLMS and decision on channel storage. +// +// +// @param aecm [i/o] Handle of the AECM instance. +// @param dfa [in] Absolute value of the nearend signal (Q[aecm->dfaQDomain]) +// @param delayDiff [in] Delay position in farend buffer. +// @param mu [in] NLMS step size. +// @param echoEst [i/o] Estimated echo +// (Q[aecm->xfaQDomain[delayDiff]+RESOLUTION_CHANNEL16]). +// +void WebRtcAecm_UpdateChannel(AecmCore_t * const aecm, const WebRtc_UWord16 * const dfa, + const WebRtc_Word16 delayDiff, const WebRtc_Word16 mu, + WebRtc_Word32 * const echoEst) +{ + + WebRtc_UWord32 tmpU32no1, tmpU32no2; + WebRtc_Word32 tmp32no1, tmp32no2; + WebRtc_Word32 mseStored; + WebRtc_Word32 mseAdapt; + + int i; + + WebRtc_Word16 zerosFar, zerosNum, zerosCh, zerosDfa; + WebRtc_Word16 shiftChFar, shiftNum, shift2ResChan; + WebRtc_Word16 tmp16no1; + WebRtc_Word16 xfaQ, dfaQ; + + // This is the channel estimation algorithm. It is base on NLMS but has a variable step + // length, which was calculated above. + if (mu) + { + for (i = 0; i < PART_LEN1; i++) + { + // Determine norm of channel and farend to make sure we don't get overflow in + // multiplication + zerosCh = WebRtcSpl_NormU32(aecm->channelAdapt32[i]); + zerosFar = WebRtcSpl_NormU32((WebRtc_UWord32)aecm->xfaHistory[i][delayDiff]); + if (zerosCh + zerosFar > 31) + { + // Multiplication is safe + tmpU32no1 = WEBRTC_SPL_UMUL_32_16(aecm->channelAdapt32[i], + aecm->xfaHistory[i][delayDiff]); + shiftChFar = 0; + } else + { + // We need to shift down before multiplication + shiftChFar = 32 - zerosCh - zerosFar; + tmpU32no1 + = WEBRTC_SPL_UMUL_32_16(WEBRTC_SPL_RSHIFT_W32(aecm->channelAdapt32[i], + shiftChFar), + aecm->xfaHistory[i][delayDiff]); + } + // Determine Q-domain of numerator + zerosNum = WebRtcSpl_NormU32(tmpU32no1); + if (dfa[i]) + { + zerosDfa = WebRtcSpl_NormU32((WebRtc_UWord32)dfa[i]); + } else + { + zerosDfa = 32; + } + tmp16no1 = zerosDfa - 2 + aecm->dfaNoisyQDomain - RESOLUTION_CHANNEL32 + - aecm->xfaQDomainBuf[delayDiff] + shiftChFar; + if (zerosNum > tmp16no1 + 1) + { + xfaQ = tmp16no1; + dfaQ = zerosDfa - 2; + } else + { + xfaQ = zerosNum - 2; + dfaQ = RESOLUTION_CHANNEL32 + aecm->xfaQDomainBuf[delayDiff] + - aecm->dfaNoisyQDomain - shiftChFar + xfaQ; + } + // Add in the same Q-domain + tmpU32no1 = WEBRTC_SPL_SHIFT_W32(tmpU32no1, xfaQ); + tmpU32no2 = WEBRTC_SPL_SHIFT_W32((WebRtc_UWord32)dfa[i], dfaQ); + tmp32no1 = (WebRtc_Word32)tmpU32no2 - (WebRtc_Word32)tmpU32no1; + zerosNum = WebRtcSpl_NormW32(tmp32no1); + if ((tmp32no1) && (aecm->xfaHistory[i][delayDiff] > (CHANNEL_VAD + << aecm->xfaQDomainBuf[delayDiff]))) + { + // + // Update is needed + // + // This is what we would like to compute + // + // tmp32no1 = dfa[i] - (aecm->channelAdapt[i] * aecm->xfaHistory[i][delayDiff]) + // tmp32norm = (i + 1) + // aecm->channelAdapt[i] += (2^mu) * tmp32no1 + // / (tmp32norm * aecm->xfaHistory[i][delayDiff]) + // + + // Make sure we don't get overflow in multiplication. + if (zerosNum + zerosFar > 31) + { + if (tmp32no1 > 0) + { + tmp32no2 = (WebRtc_Word32)WEBRTC_SPL_UMUL_32_16(tmp32no1, + aecm->xfaHistory[i][delayDiff]); + } else + { + tmp32no2 = -(WebRtc_Word32)WEBRTC_SPL_UMUL_32_16(-tmp32no1, + aecm->xfaHistory[i][delayDiff]); + } + shiftNum = 0; + } else + { + shiftNum = 32 - (zerosNum + zerosFar); + if (tmp32no1 > 0) + { + tmp32no2 = (WebRtc_Word32)WEBRTC_SPL_UMUL_32_16( + WEBRTC_SPL_RSHIFT_W32(tmp32no1, shiftNum), + aecm->xfaHistory[i][delayDiff]); + } else + { + tmp32no2 = -(WebRtc_Word32)WEBRTC_SPL_UMUL_32_16( + WEBRTC_SPL_RSHIFT_W32(-tmp32no1, shiftNum), + aecm->xfaHistory[i][delayDiff]); + } + } + // Normalize with respect to frequency bin + tmp32no2 = WebRtcSpl_DivW32W16(tmp32no2, i + 1); + // Make sure we are in the right Q-domain + shift2ResChan = shiftNum + shiftChFar - xfaQ - mu - ((30 - zerosFar) << 1); + if (WebRtcSpl_NormW32(tmp32no2) < shift2ResChan) + { + tmp32no2 = WEBRTC_SPL_WORD32_MAX; + } else + { + tmp32no2 = WEBRTC_SPL_SHIFT_W32(tmp32no2, shift2ResChan); + } + aecm->channelAdapt32[i] = WEBRTC_SPL_ADD_SAT_W32(aecm->channelAdapt32[i], + tmp32no2); + if (aecm->channelAdapt32[i] < 0) + { + // We can never have negative channel gain + aecm->channelAdapt32[i] = 0; + } + aecm->channelAdapt16[i] + = (WebRtc_Word16)WEBRTC_SPL_RSHIFT_W32(aecm->channelAdapt32[i], 16); + } + } + } + // END: Adaptive channel update + + // Determine if we should store or restore the channel + if ((aecm->startupState == 0) & (aecm->currentVADValue)) + { + // During startup we store the channel every block. + memcpy(aecm->channelStored, aecm->channelAdapt16, sizeof(WebRtc_Word16) * PART_LEN1); + // TODO(bjornv): Will be removed in final version. +#ifdef STORE_CHANNEL_DATA + fwrite(aecm->channelStored, sizeof(WebRtc_Word16), PART_LEN1, aecm->channel_file_init); +#endif + // Recalculate echo estimate +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + for (i = 0; i < PART_LEN1; i++) + { + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], + aecm->xfaHistory[i][delayDiff]); + } +#else + for (i = 0; i < PART_LEN; ) //assume PART_LEN is 4's multiples + + { + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], + aecm->xfaHistory[i][delayDiff]); + i++; + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], + aecm->xfaHistory[i][delayDiff]); + i++; + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], + aecm->xfaHistory[i][delayDiff]); + i++; + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], + aecm->xfaHistory[i][delayDiff]); + i++; + } + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], + aecm->xfaHistory[i][delayDiff]); +#endif + } else + { + if (aecm->farLogEnergy[0] < aecm->farEnergyMSE) + { + aecm->mseChannelCount = 0; + aecm->delayCount = 0; + } else + { + aecm->mseChannelCount++; + aecm->delayCount++; + } + // Enough data for validation. Store channel if we can. + if (aecm->mseChannelCount >= (MIN_MSE_COUNT + 10)) + { + // We have enough data. + // Calculate MSE of "Adapt" and "Stored" versions. + // It is actually not MSE, but average absolute error. + mseStored = 0; + mseAdapt = 0; + for (i = 0; i < MIN_MSE_COUNT; i++) + { + tmp32no1 = ((WebRtc_Word32)aecm->echoStoredLogEnergy[i] + - (WebRtc_Word32)aecm->nearLogEnergy[i]); + tmp32no2 = WEBRTC_SPL_ABS_W32(tmp32no1); + mseStored += tmp32no2; + + tmp32no1 = ((WebRtc_Word32)aecm->echoAdaptLogEnergy[i] + - (WebRtc_Word32)aecm->nearLogEnergy[i]); + tmp32no2 = WEBRTC_SPL_ABS_W32(tmp32no1); + mseAdapt += tmp32no2; + } + if (((mseStored << MSE_RESOLUTION) < (MIN_MSE_DIFF * mseAdapt)) + & ((aecm->mseStoredOld << MSE_RESOLUTION) < (MIN_MSE_DIFF + * aecm->mseAdaptOld))) + { + // The stored channel has a significantly lower MSE than the adaptive one for + // two consecutive calculations. Reset the adaptive channel. + memcpy(aecm->channelAdapt16, aecm->channelStored, + sizeof(WebRtc_Word16) * PART_LEN1); + // Restore the W32 channel +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + for (i = 0; i < PART_LEN1; i++) + { + aecm->channelAdapt32[i] + = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)aecm->channelStored[i], 16); + } +#else + for (i = 0; i < PART_LEN; ) //assume PART_LEN is 4's multiples + + { + aecm->channelAdapt32[i] = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)aecm->channelStored[i], 16); + i++; + aecm->channelAdapt32[i] = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)aecm->channelStored[i], 16); + i++; + aecm->channelAdapt32[i] = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)aecm->channelStored[i], 16); + i++; + aecm->channelAdapt32[i] = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)aecm->channelStored[i], 16); + i++; + } + aecm->channelAdapt32[i] = WEBRTC_SPL_LSHIFT_W32((WebRtc_Word32)aecm->channelStored[i], 16); +#endif + + } else if (((MIN_MSE_DIFF * mseStored) > (mseAdapt << MSE_RESOLUTION)) & (mseAdapt + < aecm->mseThreshold) & (aecm->mseAdaptOld < aecm->mseThreshold)) + { + // The adaptive channel has a significantly lower MSE than the stored one. + // The MSE for the adaptive channel has also been low for two consecutive + // calculations. Store the adaptive channel. + memcpy(aecm->channelStored, aecm->channelAdapt16, + sizeof(WebRtc_Word16) * PART_LEN1); + // TODO(bjornv): Will be removed in final version. +#ifdef STORE_CHANNEL_DATA + fwrite(aecm->channelStored, sizeof(WebRtc_Word16), PART_LEN1, + aecm->channel_file); +#endif +// Recalculate echo estimate +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + for (i = 0; i < PART_LEN1; i++) + { + echoEst[i] + = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], aecm->xfaHistory[i][delayDiff]); + } +#else + for (i = 0; i < PART_LEN; ) //assume PART_LEN is 4's multiples + + { + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], aecm->xfaHistory[i][delayDiff]); + i++; + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], aecm->xfaHistory[i][delayDiff]); + i++; + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], aecm->xfaHistory[i][delayDiff]); + i++; + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], aecm->xfaHistory[i][delayDiff]); + i++; + } + echoEst[i] = WEBRTC_SPL_MUL_16_U16(aecm->channelStored[i], aecm->xfaHistory[i][delayDiff]); +#endif + // Update threshold + if (aecm->mseThreshold == WEBRTC_SPL_WORD32_MAX) + { + aecm->mseThreshold = (mseAdapt + aecm->mseAdaptOld); + } else + { + aecm->mseThreshold += WEBRTC_SPL_MUL_16_16_RSFT(mseAdapt + - WEBRTC_SPL_MUL_16_16_RSFT(aecm->mseThreshold, 5, 3), 205, 8); + } + + } + + // Reset counter + aecm->mseChannelCount = 0; + + // Store the MSE values. + aecm->mseStoredOld = mseStored; + aecm->mseAdaptOld = mseAdapt; + } + } + // END: Determine if we should store or reset channel estimate. +} + +// WebRtcAecm_CalcSuppressionGain(...) +// +// This function calculates the suppression gain that is used in the Wiener filter. +// +// +// @param aecm [i/n] Handle of the AECM instance. +// @param supGain [out] (Return value) Suppression gain with which to scale the noise +// level (Q14). +// +// +WebRtc_Word16 WebRtcAecm_CalcSuppressionGain(AecmCore_t * const aecm) +{ + WebRtc_Word32 tmp32no1; + + WebRtc_Word16 supGain; + WebRtc_Word16 tmp16no1; + WebRtc_Word16 dE = 0; + + // Determine suppression gain used in the Wiener filter. The gain is based on a mix of far + // end energy and echo estimation error. + supGain = SUPGAIN_DEFAULT; + // Adjust for the far end signal level. A low signal level indicates no far end signal, + // hence we set the suppression gain to 0 + if (!aecm->currentVADValue) + { + supGain = 0; + } else + { + // Adjust for possible double talk. If we have large variations in estimation error we + // likely have double talk (or poor channel). + tmp16no1 = (aecm->nearLogEnergy[0] - aecm->echoStoredLogEnergy[0] - ENERGY_DEV_OFFSET); + dE = WEBRTC_SPL_ABS_W16(tmp16no1); + + if (dE < ENERGY_DEV_TOL) + { + // Likely no double talk. The better estimation, the more we can suppress signal. + // Update counters + if (dE < SUPGAIN_EPC_DT) + { + tmp32no1 = WEBRTC_SPL_MUL_16_16(aecm->supGainErrParamDiffAB, dE); + tmp32no1 += (SUPGAIN_EPC_DT >> 1); + tmp16no1 = (WebRtc_Word16)WebRtcSpl_DivW32W16(tmp32no1, SUPGAIN_EPC_DT); + supGain = aecm->supGainErrParamA - tmp16no1; + } else + { + tmp32no1 = WEBRTC_SPL_MUL_16_16(aecm->supGainErrParamDiffBD, + (ENERGY_DEV_TOL - dE)); + tmp32no1 += ((ENERGY_DEV_TOL - SUPGAIN_EPC_DT) >> 1); + tmp16no1 = (WebRtc_Word16)WebRtcSpl_DivW32W16(tmp32no1, (ENERGY_DEV_TOL + - SUPGAIN_EPC_DT)); + supGain = aecm->supGainErrParamD + tmp16no1; + } + } else + { + // Likely in double talk. Use default value + supGain = aecm->supGainErrParamD; + } + } + + if (supGain > aecm->supGainOld) + { + tmp16no1 = supGain; + } else + { + tmp16no1 = aecm->supGainOld; + } + aecm->supGainOld = supGain; + if (tmp16no1 < aecm->supGain) + { + aecm->supGain += (WebRtc_Word16)((tmp16no1 - aecm->supGain) >> 4); + } else + { + aecm->supGain += (WebRtc_Word16)((tmp16no1 - aecm->supGain) >> 4); + } + + // END: Update suppression gain + + return aecm->supGain; +} + +// WebRtcAecm_DelayCompensation(...) +// +// Secondary delay estimation that can be used as a backup or for validation. This function is +// still under construction and not activated in current version. +// +// +// @param aecm [i/o] Handle of the AECM instance. +// +// +void WebRtcAecm_DelayCompensation(AecmCore_t * const aecm) +{ + int i, j; + WebRtc_Word32 delayMeanEcho[CORR_BUF_LEN]; + WebRtc_Word32 delayMeanNear[CORR_BUF_LEN]; + WebRtc_Word16 sumBitPattern, bitPatternEcho, bitPatternNear, maxPos, maxValue, + maxValueLeft, maxValueRight; + + // Check delay (calculate the delay offset (if we can)). + if ((aecm->startupState > 0) & (aecm->delayCount >= CORR_MAX_BUF) & aecm->delayOffsetFlag) + { + // Calculate mean values + for (i = 0; i < CORR_BUF_LEN; i++) + { + delayMeanEcho[i] = 0; + delayMeanNear[i] = 0; +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + for (j = 0; j < CORR_WIDTH; j++) + { + delayMeanEcho[i] += (WebRtc_Word32)aecm->echoStoredLogEnergy[i + j]; + delayMeanNear[i] += (WebRtc_Word32)aecm->nearLogEnergy[i + j]; + } +#else + for (j = 0; j < CORR_WIDTH -1; ) + { + delayMeanEcho[i] += (WebRtc_Word32)aecm->echoStoredLogEnergy[i + j]; + delayMeanNear[i] += (WebRtc_Word32)aecm->nearLogEnergy[i + j]; + j++; + delayMeanEcho[i] += (WebRtc_Word32)aecm->echoStoredLogEnergy[i + j]; + delayMeanNear[i] += (WebRtc_Word32)aecm->nearLogEnergy[i + j]; + j++; + } + delayMeanEcho[i] += (WebRtc_Word32)aecm->echoStoredLogEnergy[i + j]; + delayMeanNear[i] += (WebRtc_Word32)aecm->nearLogEnergy[i + j]; +#endif + } + // Calculate correlation values + for (i = 0; i < CORR_BUF_LEN; i++) + { + sumBitPattern = 0; +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + for (j = 0; j < CORR_WIDTH; j++) + { + bitPatternEcho = (WebRtc_Word16)((WebRtc_Word32)aecm->echoStoredLogEnergy[i + + j] * CORR_WIDTH > delayMeanEcho[i]); + bitPatternNear = (WebRtc_Word16)((WebRtc_Word32)aecm->nearLogEnergy[CORR_MAX + + j] * CORR_WIDTH > delayMeanNear[CORR_MAX]); + sumBitPattern += !(bitPatternEcho ^ bitPatternNear); + } +#else + for (j = 0; j < CORR_WIDTH -1; ) + { + bitPatternEcho = (WebRtc_Word16)((WebRtc_Word32)aecm->echoStoredLogEnergy[i + + j] * CORR_WIDTH > delayMeanEcho[i]); + bitPatternNear = (WebRtc_Word16)((WebRtc_Word32)aecm->nearLogEnergy[CORR_MAX + + j] * CORR_WIDTH > delayMeanNear[CORR_MAX]); + sumBitPattern += !(bitPatternEcho ^ bitPatternNear); + j++; + bitPatternEcho = (WebRtc_Word16)((WebRtc_Word32)aecm->echoStoredLogEnergy[i + + j] * CORR_WIDTH > delayMeanEcho[i]); + bitPatternNear = (WebRtc_Word16)((WebRtc_Word32)aecm->nearLogEnergy[CORR_MAX + + j] * CORR_WIDTH > delayMeanNear[CORR_MAX]); + sumBitPattern += !(bitPatternEcho ^ bitPatternNear); + j++; + } + bitPatternEcho = (WebRtc_Word16)((WebRtc_Word32)aecm->echoStoredLogEnergy[i + j] + * CORR_WIDTH > delayMeanEcho[i]); + bitPatternNear = (WebRtc_Word16)((WebRtc_Word32)aecm->nearLogEnergy[CORR_MAX + j] + * CORR_WIDTH > delayMeanNear[CORR_MAX]); + sumBitPattern += !(bitPatternEcho ^ bitPatternNear); +#endif + aecm->delayCorrelation[i] = sumBitPattern; + } + aecm->newDelayCorrData = 1; // Indicate we have new correlation data to evaluate + } + if ((aecm->startupState == 2) & (aecm->lastDelayUpdateCount > (CORR_WIDTH << 1)) + & aecm->newDelayCorrData) + { + // Find maximum value and maximum position as well as values on the sides. + maxPos = 0; + maxValue = aecm->delayCorrelation[0]; + maxValueLeft = maxValue; + maxValueRight = aecm->delayCorrelation[CORR_DEV]; + for (i = 1; i < CORR_BUF_LEN; i++) + { + if (aecm->delayCorrelation[i] > maxValue) + { + maxValue = aecm->delayCorrelation[i]; + maxPos = i; + if (maxPos < CORR_DEV) + { + maxValueLeft = aecm->delayCorrelation[0]; + maxValueRight = aecm->delayCorrelation[i + CORR_DEV]; + } else if (maxPos > (CORR_MAX << 1) - CORR_DEV) + { + maxValueLeft = aecm->delayCorrelation[i - CORR_DEV]; + maxValueRight = aecm->delayCorrelation[(CORR_MAX << 1)]; + } else + { + maxValueLeft = aecm->delayCorrelation[i - CORR_DEV]; + maxValueRight = aecm->delayCorrelation[i + CORR_DEV]; + } + } + } + if ((maxPos > 0) & (maxPos < (CORR_MAX << 1))) + { + // Avoid maximum at boundaries. The maximum peak has to be higher than + // CORR_MAX_LEVEL. It also has to be sharp, i.e. the value CORR_DEV bins off should + // be CORR_MAX_LOW lower than the maximum. + if ((maxValue > CORR_MAX_LEVEL) & (maxValueLeft < maxValue - CORR_MAX_LOW) + & (maxValueRight < maxValue - CORR_MAX_LOW)) + { + aecm->delayAdjust += CORR_MAX - maxPos; + aecm->newDelayCorrData = 0; + aecm->lastDelayUpdateCount = 0; + } + } + } + // END: "Check delay" +} + +void WebRtcAecm_ProcessBlock(AecmCore_t * const aecm, const WebRtc_Word16 * const farend, + const WebRtc_Word16 * const nearendNoisy, + const WebRtc_Word16 * const nearendClean, + WebRtc_Word16 * const output) +{ + int i, j; + + WebRtc_UWord32 xfaSum; + WebRtc_UWord32 dfaNoisySum; + WebRtc_UWord32 echoEst32Gained; + WebRtc_UWord32 tmpU32; + + WebRtc_Word32 tmp32no1; + WebRtc_Word32 tmp32no2; + WebRtc_Word32 echoEst32[PART_LEN1]; + + WebRtc_UWord16 xfa[PART_LEN1]; + WebRtc_UWord16 dfaNoisy[PART_LEN1]; + WebRtc_UWord16 dfaClean[PART_LEN1]; + WebRtc_UWord16* ptrDfaClean = dfaClean; + + int outCFFT; + + WebRtc_Word16 fft[PART_LEN4]; +#if (defined ARM_WINM) || (defined ARM9E_GCC) || (defined ANDROID_AECOPT) + WebRtc_Word16 postFft[PART_LEN4]; +#else + WebRtc_Word16 postFft[PART_LEN2]; +#endif + WebRtc_Word16 dfwReal[PART_LEN1]; + WebRtc_Word16 dfwImag[PART_LEN1]; + WebRtc_Word16 xfwReal[PART_LEN1]; + WebRtc_Word16 xfwImag[PART_LEN1]; + WebRtc_Word16 efwReal[PART_LEN1]; + WebRtc_Word16 efwImag[PART_LEN1]; + WebRtc_Word16 hnl[PART_LEN1]; + WebRtc_Word16 numPosCoef; + WebRtc_Word16 nlpGain; + WebRtc_Word16 delay, diff, diffMinusOne; + WebRtc_Word16 tmp16no1; + WebRtc_Word16 tmp16no2; +#ifdef AECM_WITH_ABS_APPROX + WebRtc_Word16 maxValue; + WebRtc_Word16 minValue; +#endif + WebRtc_Word16 mu; + WebRtc_Word16 supGain; + WebRtc_Word16 zeros32, zeros16; + WebRtc_Word16 zerosDBufNoisy, zerosDBufClean, zerosXBuf; + WebRtc_Word16 resolutionDiff, qDomainDiff; + +#ifdef ARM_WINM_LOG_ + DWORD temp; + static int flag0 = 0; + __int64 freq, start, end, diff__; + unsigned int milliseconds; +#endif + +#ifdef AECM_WITH_ABS_APPROX + WebRtc_UWord16 alpha, beta; +#endif + + // Determine startup state. There are three states: + // (0) the first CONV_LEN blocks + // (1) another CONV_LEN blocks + // (2) the rest + + if (aecm->startupState < 2) + { + aecm->startupState = (aecm->totCount >= CONV_LEN) + (aecm->totCount >= CONV_LEN2); + } + // END: Determine startup state + + // Buffer near and far end signals + memcpy(aecm->xBuf + PART_LEN, farend, sizeof(WebRtc_Word16) * PART_LEN); + memcpy(aecm->dBufNoisy + PART_LEN, nearendNoisy, sizeof(WebRtc_Word16) * PART_LEN); + if (nearendClean != NULL) + { + memcpy(aecm->dBufClean + PART_LEN, nearendClean, sizeof(WebRtc_Word16) * PART_LEN); + } + // TODO(bjornv): Will be removed in final version. +#ifdef VAD_DATA + fwrite(aecm->xBuf, sizeof(WebRtc_Word16), PART_LEN, aecm->far_file); +#endif + +#ifdef AECM_DYNAMIC_Q + tmp16no1 = WebRtcSpl_MaxAbsValueW16(aecm->dBufNoisy, PART_LEN2); + tmp16no2 = WebRtcSpl_MaxAbsValueW16(aecm->xBuf, PART_LEN2); + zerosDBufNoisy = WebRtcSpl_NormW16(tmp16no1); + zerosXBuf = WebRtcSpl_NormW16(tmp16no2); +#else + zerosDBufNoisy = 0; + zerosXBuf = 0; +#endif + aecm->dfaNoisyQDomainOld = aecm->dfaNoisyQDomain; + aecm->dfaNoisyQDomain = zerosDBufNoisy; + + if (nearendClean != NULL) + { +#ifdef AECM_DYNAMIC_Q + tmp16no1 = WebRtcSpl_MaxAbsValueW16(aecm->dBufClean, PART_LEN2); + zerosDBufClean = WebRtcSpl_NormW16(tmp16no1); +#else + zerosDBufClean = 0; +#endif + aecm->dfaCleanQDomainOld = aecm->dfaCleanQDomain; + aecm->dfaCleanQDomain = zerosDBufClean; + } else + { + zerosDBufClean = zerosDBufNoisy; + aecm->dfaCleanQDomainOld = aecm->dfaNoisyQDomainOld; + aecm->dfaCleanQDomain = aecm->dfaNoisyQDomain; + } + +#ifdef ARM_WINM_LOG_ + // measure tick start + QueryPerformanceFrequency((LARGE_INTEGER*)&freq); + QueryPerformanceCounter((LARGE_INTEGER*)&start); +#endif + + // FFT of noisy near end signal + for (i = 0; i < PART_LEN; i++) + { + j = WEBRTC_SPL_LSHIFT_W32(i, 1); + // Window near end + fft[j] = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT((aecm->dBufNoisy[i] + << zerosDBufNoisy), kSqrtHanning[i], 14); + fft[PART_LEN2 + j] = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT( + (aecm->dBufNoisy[PART_LEN + i] << zerosDBufNoisy), + kSqrtHanning[PART_LEN - i], 14); + // Inserting zeros in imaginary parts + fft[j + 1] = 0; + fft[PART_LEN2 + j + 1] = 0; + } + + // Fourier transformation of near end signal. + // The result is scaled with 1/PART_LEN2, that is, the result is in Q(-6) for PART_LEN = 32 + +#if (defined ARM_WINM) || (defined ARM9E_GCC) || (defined ANDROID_AECOPT) + outCFFT = WebRtcSpl_ComplexFFT2(fft, postFft, PART_LEN_SHIFT, 1); + + // The imaginary part has to switch sign + for(i = 1; i < PART_LEN2-1;) + { + postFft[i] = -postFft[i]; + i += 2; + postFft[i] = -postFft[i]; + i += 2; + } +#else + WebRtcSpl_ComplexBitReverse(fft, PART_LEN_SHIFT); + outCFFT = WebRtcSpl_ComplexFFT(fft, PART_LEN_SHIFT, 1); + + // Take only the first PART_LEN2 samples + for (i = 0; i < PART_LEN2; i++) + { + postFft[i] = fft[i]; + } + // The imaginary part has to switch sign + for (i = 1; i < PART_LEN2;) + { + postFft[i] = -postFft[i]; + i += 2; + } +#endif + + // Extract imaginary and real part, calculate the magnitude for all frequency bins + dfwImag[0] = 0; + dfwImag[PART_LEN] = 0; + dfwReal[0] = postFft[0]; +#if (defined ARM_WINM) || (defined ARM9E_GCC) || (defined ANDROID_AECOPT) + dfwReal[PART_LEN] = postFft[PART_LEN2]; +#else + dfwReal[PART_LEN] = fft[PART_LEN2]; +#endif + dfaNoisy[0] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(dfwReal[0]); + dfaNoisy[PART_LEN] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(dfwReal[PART_LEN]); + dfaNoisySum = (WebRtc_UWord32)(dfaNoisy[0]); + dfaNoisySum += (WebRtc_UWord32)(dfaNoisy[PART_LEN]); + + for (i = 1; i < PART_LEN; i++) + { + j = WEBRTC_SPL_LSHIFT_W32(i, 1); + dfwReal[i] = postFft[j]; + dfwImag[i] = postFft[j + 1]; + + if (dfwReal[i] == 0 || dfwImag[i] == 0) + { + dfaNoisy[i] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(dfwReal[i] + dfwImag[i]); + } else + { + // Approximation for magnitude of complex fft output + // magn = sqrt(real^2 + imag^2) + // magn ~= alpha * max(|imag|,|real|) + beta * min(|imag|,|real|) + // + // The parameters alpha and beta are stored in Q15 + + tmp16no1 = WEBRTC_SPL_ABS_W16(postFft[j]); + tmp16no2 = WEBRTC_SPL_ABS_W16(postFft[j + 1]); + +#ifdef AECM_WITH_ABS_APPROX + if(tmp16no1 > tmp16no2) + { + maxValue = tmp16no1; + minValue = tmp16no2; + } else + { + maxValue = tmp16no2; + minValue = tmp16no1; + } + + // Magnitude in Q-6 + if ((maxValue >> 2) > minValue) + { + alpha = kAlpha1; + beta = kBeta1; + } else if ((maxValue >> 1) > minValue) + { + alpha = kAlpha2; + beta = kBeta2; + } else + { + alpha = kAlpha3; + beta = kBeta3; + } + tmp16no1 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(maxValue, alpha, 15); + tmp16no2 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(minValue, beta, 15); + dfaNoisy[i] = (WebRtc_UWord16)tmp16no1 + (WebRtc_UWord16)tmp16no2; +#else + tmp32no1 = WEBRTC_SPL_MUL_16_16(tmp16no1, tmp16no1); + tmp32no2 = WEBRTC_SPL_MUL_16_16(tmp16no2, tmp16no2); + tmp32no2 = WEBRTC_SPL_ADD_SAT_W32(tmp32no1, tmp32no2); + tmp32no1 = WebRtcSpl_Sqrt(tmp32no2); + dfaNoisy[i] = (WebRtc_UWord16)tmp32no1; +#endif + } + dfaNoisySum += (WebRtc_UWord32)dfaNoisy[i]; + } + // END: FFT of noisy near end signal + + if (nearendClean == NULL) + { + ptrDfaClean = dfaNoisy; + } else + { + // FFT of clean near end signal + for (i = 0; i < PART_LEN; i++) + { + j = WEBRTC_SPL_LSHIFT_W32(i, 1); + // Window near end + fft[j] + = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT((aecm->dBufClean[i] << zerosDBufClean), kSqrtHanning[i], 14); + fft[PART_LEN2 + j] + = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT((aecm->dBufClean[PART_LEN + i] << zerosDBufClean), kSqrtHanning[PART_LEN - i], 14); + // Inserting zeros in imaginary parts + fft[j + 1] = 0; + fft[PART_LEN2 + j + 1] = 0; + } + + // Fourier transformation of near end signal. + // The result is scaled with 1/PART_LEN2, that is, in Q(-6) for PART_LEN = 32 + +#if (defined ARM_WINM) || (defined ARM9E_GCC) || (defined ANDROID_AECOPT) + outCFFT = WebRtcSpl_ComplexFFT2(fft, postFft, PART_LEN_SHIFT, 1); + + // The imaginary part has to switch sign + for(i = 1; i < PART_LEN2-1;) + { + postFft[i] = -postFft[i]; + i += 2; + postFft[i] = -postFft[i]; + i += 2; + } +#else + WebRtcSpl_ComplexBitReverse(fft, PART_LEN_SHIFT); + outCFFT = WebRtcSpl_ComplexFFT(fft, PART_LEN_SHIFT, 1); + + // Take only the first PART_LEN2 samples + for (i = 0; i < PART_LEN2; i++) + { + postFft[i] = fft[i]; + } + // The imaginary part has to switch sign + for (i = 1; i < PART_LEN2;) + { + postFft[i] = -postFft[i]; + i += 2; + } +#endif + + // Extract imaginary and real part, calculate the magnitude for all frequency bins + dfwImag[0] = 0; + dfwImag[PART_LEN] = 0; + dfwReal[0] = postFft[0]; +#if (defined ARM_WINM) || (defined ARM9E_GCC) || (defined ANDROID_AECOPT) + dfwReal[PART_LEN] = postFft[PART_LEN2]; +#else + dfwReal[PART_LEN] = fft[PART_LEN2]; +#endif + dfaClean[0] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(dfwReal[0]); + dfaClean[PART_LEN] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(dfwReal[PART_LEN]); + + for (i = 1; i < PART_LEN; i++) + { + j = WEBRTC_SPL_LSHIFT_W32(i, 1); + dfwReal[i] = postFft[j]; + dfwImag[i] = postFft[j + 1]; + + if (dfwReal[i] == 0 || dfwImag[i] == 0) + { + dfaClean[i] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(dfwReal[i] + dfwImag[i]); + } else + { + // Approximation for magnitude of complex fft output + // magn = sqrt(real^2 + imag^2) + // magn ~= alpha * max(|imag|,|real|) + beta * min(|imag|,|real|) + // + // The parameters alpha and beta are stored in Q15 + + tmp16no1 = WEBRTC_SPL_ABS_W16(postFft[j]); + tmp16no2 = WEBRTC_SPL_ABS_W16(postFft[j + 1]); + +#ifdef AECM_WITH_ABS_APPROX + if(tmp16no1 > tmp16no2) + { + maxValue = tmp16no1; + minValue = tmp16no2; + } else + { + maxValue = tmp16no2; + minValue = tmp16no1; + } + + // Magnitude in Q-6 + if ((maxValue >> 2) > minValue) + { + alpha = kAlpha1; + beta = kBeta1; + } else if ((maxValue >> 1) > minValue) + { + alpha = kAlpha2; + beta = kBeta2; + } else + { + alpha = kAlpha3; + beta = kBeta3; + } + tmp16no1 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(maxValue, alpha, 15); + tmp16no2 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(minValue, beta, 15); + dfaClean[i] = (WebRtc_UWord16)tmp16no1 + (WebRtc_UWord16)tmp16no2; +#else + tmp32no1 = WEBRTC_SPL_MUL_16_16(tmp16no1, tmp16no1); + tmp32no2 = WEBRTC_SPL_MUL_16_16(tmp16no2, tmp16no2); + tmp32no2 = WEBRTC_SPL_ADD_SAT_W32(tmp32no1, tmp32no2); + tmp32no1 = WebRtcSpl_Sqrt(tmp32no2); + dfaClean[i] = (WebRtc_UWord16)tmp32no1; +#endif + } + } + } + // END: FFT of clean near end signal + + // FFT of far end signal + for (i = 0; i < PART_LEN; i++) + { + j = WEBRTC_SPL_LSHIFT_W32(i, 1); + // Window farend + fft[j] + = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT((aecm->xBuf[i] << zerosXBuf), kSqrtHanning[i], 14); + fft[PART_LEN2 + j] + = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT((aecm->xBuf[PART_LEN + i] << zerosXBuf), kSqrtHanning[PART_LEN - i], 14); + // Inserting zeros in imaginary parts + fft[j + 1] = 0; + fft[PART_LEN2 + j + 1] = 0; + } + // Fourier transformation of far end signal. + // The result is scaled with 1/PART_LEN2, that is the result is in Q(-6) for PART_LEN = 32 +#if (defined ARM_WINM) || (defined ARM9E_GCC) || (defined ANDROID_AECOPT) + outCFFT = WebRtcSpl_ComplexFFT2(fft, postFft, PART_LEN_SHIFT, 1); + + // The imaginary part has to switch sign + for(i = 1; i < PART_LEN2-1;) + { + postFft[i] = -postFft[i]; + i += 2; + postFft[i] = -postFft[i]; + i += 2; + } +#else + WebRtcSpl_ComplexBitReverse(fft, PART_LEN_SHIFT); + outCFFT = WebRtcSpl_ComplexFFT(fft, PART_LEN_SHIFT, 1); + + // Take only the first PART_LEN2 samples + for (i = 0; i < PART_LEN2; i++) + { + postFft[i] = fft[i]; + } + // The imaginary part has to switch sign + for (i = 1; i < PART_LEN2;) + { + postFft[i] = -postFft[i]; + i += 2; + } +#endif + + // Extract imaginary and real part, calculate the magnitude for all frequency bins + xfwImag[0] = 0; + xfwImag[PART_LEN] = 0; + xfwReal[0] = postFft[0]; +#if (defined ARM_WINM) || (defined ARM9E_GCC) || (defined ANDROID_AECOPT) + xfwReal[PART_LEN] = postFft[PART_LEN2]; +#else + xfwReal[PART_LEN] = fft[PART_LEN2]; +#endif + xfa[0] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(xfwReal[0]); + xfa[PART_LEN] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(xfwReal[PART_LEN]); + xfaSum = (WebRtc_UWord32)(xfa[0]) + (WebRtc_UWord32)(xfa[PART_LEN]); + + for (i = 1; i < PART_LEN; i++) + { + j = WEBRTC_SPL_LSHIFT_W32(i,1); + xfwReal[i] = postFft[j]; + xfwImag[i] = postFft[j + 1]; + + if (xfwReal[i] == 0 || xfwImag[i] == 0) + { + xfa[i] = (WebRtc_UWord16)WEBRTC_SPL_ABS_W16(xfwReal[i] + xfwImag[i]); + } else + { + // Approximation for magnitude of complex fft output + // magn = sqrt(real^2 + imag^2) + // magn ~= alpha * max(|imag|,|real|) + beta * min(|imag|,|real|) + // + // The parameters alpha and beta are stored in Q15 + + tmp16no1 = WEBRTC_SPL_ABS_W16(postFft[j]); + tmp16no2 = WEBRTC_SPL_ABS_W16(postFft[j + 1]); + +#ifdef AECM_WITH_ABS_APPROX + if(tmp16no1 > xfwImag[i]) + { + maxValue = tmp16no1; + minValue = tmp16no2; + } else + { + maxValue = tmp16no2; + minValue = tmp16no1; + } + // Magnitude in Q-6 + if ((maxValue >> 2) > minValue) + { + alpha = kAlpha1; + beta = kBeta1; + } else if ((maxValue >> 1) > minValue) + { + alpha = kAlpha2; + beta = kBeta2; + } else + { + alpha = kAlpha3; + beta = kBeta3; + } + tmp16no1 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(maxValue, alpha, 15); + tmp16no2 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(minValue, beta, 15); + xfa[i] = (WebRtc_UWord16)tmp16no1 + (WebRtc_UWord16)tmp16no2; +#else + tmp32no1 = WEBRTC_SPL_MUL_16_16(tmp16no1, tmp16no1); + tmp32no2 = WEBRTC_SPL_MUL_16_16(tmp16no2, tmp16no2); + tmp32no2 = WEBRTC_SPL_ADD_SAT_W32(tmp32no1, tmp32no2); + tmp32no1 = WebRtcSpl_Sqrt(tmp32no2); + xfa[i] = (WebRtc_UWord16)tmp32no1; +#endif + } + xfaSum += (WebRtc_UWord32)xfa[i]; + } + +#ifdef ARM_WINM_LOG_ + // measure tick end + QueryPerformanceCounter((LARGE_INTEGER*)&end); + diff__ = ((end - start) * 1000) / (freq/1000); + milliseconds = (unsigned int)(diff__ & 0xffffffff); + WriteFile (logFile, &milliseconds, sizeof(unsigned int), &temp, NULL); +#endif + // END: FFT of far end signal + + // Get the delay + + // Fixed delay estimation + // input: dfaFIX, xfaFIX in Q-stages + // output: delay in Q0 + // + // comment on the fixed point accuracy of estimate_delayFIX + // -> due to rounding the fixed point variables xfa and dfa contain a lot more zeros + // than the corresponding floating point variables this results in big differences + // between the floating point and the fixed point logarithmic spectra for small values +#ifdef ARM_WINM_LOG_ + // measure tick start + QueryPerformanceCounter((LARGE_INTEGER*)&start); +#endif + + // Save far-end history and estimate delay + delay = WebRtcAecm_EstimateDelay(aecm, xfa, dfaNoisy, zerosXBuf); + + if (aecm->fixedDelay >= 0) + { + // Use fixed delay + delay = aecm->fixedDelay; + } + + aecm->currentDelay = delay; + + if ((aecm->delayOffsetFlag) & (aecm->startupState > 0)) // If delay compensation is on + { + // If the delay estimate changed from previous block, update the offset + if ((aecm->currentDelay != aecm->previousDelay) & !aecm->currentDelay + & !aecm->previousDelay) + { + aecm->delayAdjust += (aecm->currentDelay - aecm->previousDelay); + } + // Compensate with the offset estimate + aecm->currentDelay -= aecm->delayAdjust; + aecm->previousDelay = delay; + } + + diff = aecm->delHistoryPos - aecm->currentDelay; + if (diff < 0) + { + diff = diff + MAX_DELAY; + } + +#ifdef ARM_WINM_LOG_ + // measure tick end + QueryPerformanceCounter((LARGE_INTEGER*)&end); + diff__ = ((end - start) * 1000) / (freq/1000); + milliseconds = (unsigned int)(diff__ & 0xffffffff); + WriteFile (logFile, &milliseconds, sizeof(unsigned int), &temp, NULL); +#endif + + // END: Get the delay + +#ifdef ARM_WINM_LOG_ + // measure tick start + QueryPerformanceCounter((LARGE_INTEGER*)&start); +#endif + // Calculate log(energy) and update energy threshold levels + WebRtcAecm_CalcEnergies(aecm, diff, dfaNoisySum, echoEst32); + + // Calculate stepsize + mu = WebRtcAecm_CalcStepSize(aecm); + + // Update counters + aecm->totCount++; + aecm->lastDelayUpdateCount++; + + // This is the channel estimation algorithm. + // It is base on NLMS but has a variable step length, which was calculated above. + WebRtcAecm_UpdateChannel(aecm, dfaNoisy, diff, mu, echoEst32); + WebRtcAecm_DelayCompensation(aecm); + supGain = WebRtcAecm_CalcSuppressionGain(aecm); + +#ifdef ARM_WINM_LOG_ + // measure tick end + QueryPerformanceCounter((LARGE_INTEGER*)&end); + diff__ = ((end - start) * 1000) / (freq/1000); + milliseconds = (unsigned int)(diff__ & 0xffffffff); + WriteFile (logFile, &milliseconds, sizeof(unsigned int), &temp, NULL); +#endif + +#ifdef ARM_WINM_LOG_ + // measure tick start + QueryPerformanceCounter((LARGE_INTEGER*)&start); +#endif + + // Calculate Wiener filter hnl[] + numPosCoef = 0; + diffMinusOne = diff - 1; + if (diff == 0) + { + diffMinusOne = MAX_DELAY; + } + for (i = 0; i < PART_LEN1; i++) + { + // Far end signal through channel estimate in Q8 + // How much can we shift right to preserve resolution + tmp32no1 = echoEst32[i] - aecm->echoFilt[i]; + aecm->echoFilt[i] += WEBRTC_SPL_RSHIFT_W32(WEBRTC_SPL_MUL_32_16(tmp32no1, 50), 8); + + zeros32 = WebRtcSpl_NormW32(aecm->echoFilt[i]) + 1; + zeros16 = WebRtcSpl_NormW16(supGain) + 1; + if (zeros32 + zeros16 > 16) + { + // Multiplication is safe + // Result in Q(RESOLUTION_CHANNEL+RESOLUTION_SUPGAIN+aecm->xfaQDomainBuf[diff]) + echoEst32Gained = WEBRTC_SPL_UMUL_32_16((WebRtc_UWord32)aecm->echoFilt[i], + (WebRtc_UWord16)supGain); + resolutionDiff = 14 - RESOLUTION_CHANNEL16 - RESOLUTION_SUPGAIN; + resolutionDiff += (aecm->dfaCleanQDomain - aecm->xfaQDomainBuf[diff]); + } else + { + tmp16no1 = 17 - zeros32 - zeros16; + resolutionDiff = 14 + tmp16no1 - RESOLUTION_CHANNEL16 - RESOLUTION_SUPGAIN; + resolutionDiff += (aecm->dfaCleanQDomain - aecm->xfaQDomainBuf[diff]); + if (zeros32 > tmp16no1) + { + echoEst32Gained = WEBRTC_SPL_UMUL_32_16((WebRtc_UWord32)aecm->echoFilt[i], + (WebRtc_UWord16)WEBRTC_SPL_RSHIFT_W16(supGain, + tmp16no1)); // Q-(RESOLUTION_CHANNEL+RESOLUTION_SUPGAIN-16) + } else + { + // Result in Q-(RESOLUTION_CHANNEL+RESOLUTION_SUPGAIN-16) + echoEst32Gained = WEBRTC_SPL_UMUL_32_16( + (WebRtc_UWord32)WEBRTC_SPL_RSHIFT_W32(aecm->echoFilt[i], tmp16no1), + (WebRtc_UWord16)supGain); + } + } + + zeros16 = WebRtcSpl_NormW16(aecm->nearFilt[i]); + if ((zeros16 < (aecm->dfaCleanQDomain - aecm->dfaCleanQDomainOld)) + & (aecm->nearFilt[i])) + { + tmp16no1 = WEBRTC_SPL_SHIFT_W16(aecm->nearFilt[i], zeros16); + qDomainDiff = zeros16 - aecm->dfaCleanQDomain + aecm->dfaCleanQDomainOld; + } else + { + tmp16no1 = WEBRTC_SPL_SHIFT_W16(aecm->nearFilt[i], aecm->dfaCleanQDomain + - aecm->dfaCleanQDomainOld); + qDomainDiff = 0; + } + tmp16no2 = WEBRTC_SPL_SHIFT_W16(ptrDfaClean[i], qDomainDiff); + tmp16no2 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(tmp16no2 - tmp16no1, 1, 4); + tmp16no2 += tmp16no1; + zeros16 = WebRtcSpl_NormW16(tmp16no2); + if ((tmp16no2) & (-qDomainDiff > zeros16)) + { + aecm->nearFilt[i] = WEBRTC_SPL_WORD16_MAX; + } else + { + aecm->nearFilt[i] = WEBRTC_SPL_SHIFT_W16(tmp16no2, -qDomainDiff); + } + + // Wiener filter coefficients, resulting hnl in Q14 + if (echoEst32Gained == 0) + { + hnl[i] = ONE_Q14; + } else if (aecm->nearFilt[i] == 0) + { + hnl[i] = 0; + } else + { + // Multiply the suppression gain + // Rounding + echoEst32Gained += (WebRtc_UWord32)(aecm->nearFilt[i] >> 1); + tmpU32 = WebRtcSpl_DivU32U16(echoEst32Gained, (WebRtc_UWord16)aecm->nearFilt[i]); + + // Current resolution is + // Q-(RESOLUTION_CHANNEL + RESOLUTION_SUPGAIN - max(0, 17 - zeros16 - zeros32)) + // Make sure we are in Q14 + tmp32no1 = (WebRtc_Word32)WEBRTC_SPL_SHIFT_W32(tmpU32, resolutionDiff); + if (tmp32no1 > ONE_Q14) + { + hnl[i] = 0; + } else if (tmp32no1 < 0) + { + hnl[i] = ONE_Q14; + } else + { + // 1-echoEst/dfa +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + hnl[i] = ONE_Q14 - (WebRtc_Word16)tmp32no1; + if (hnl[i] < 0) + { + hnl[i] = 0; + } +#else + hnl[i] = ((ONE_Q14 - (WebRtc_Word16)tmp32no1) > 0) ? (ONE_Q14 - (WebRtc_Word16)tmp32no1) : 0; +#endif + } + } + if (hnl[i]) + { + numPosCoef++; + } + } + +#ifdef ARM_WINM_LOG_ + // measure tick end + QueryPerformanceCounter((LARGE_INTEGER*)&end); + diff__ = ((end - start) * 1000) / (freq/1000); + milliseconds = (unsigned int)(diff__ & 0xffffffff); + WriteFile (logFile, &milliseconds, sizeof(unsigned int), &temp, NULL); +#endif + +#ifdef ARM_WINM_LOG_ + // measure tick start + QueryPerformanceCounter((LARGE_INTEGER*)&start); +#endif + + // Calculate NLP gain, result is in Q14 + for (i = 0; i < PART_LEN1; i++) + { + if (aecm->nlpFlag) + { + // Truncate values close to zero and one. + if (hnl[i] > NLP_COMP_HIGH) + { + hnl[i] = ONE_Q14; + } else if (hnl[i] < NLP_COMP_LOW) + { + hnl[i] = 0; + } + + // Remove outliers + if (numPosCoef < 3) + { + nlpGain = 0; + } else + { + nlpGain = ONE_Q14; + } + // NLP + if ((hnl[i] == ONE_Q14) && (nlpGain == ONE_Q14)) + { + hnl[i] = ONE_Q14; + } else + { + hnl[i] = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(hnl[i], nlpGain, 14); + } + } + + // multiply with Wiener coefficients + efwReal[i] = (WebRtc_Word16)(WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(dfwReal[i], hnl[i], + 14)); + efwImag[i] = (WebRtc_Word16)(WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND(dfwImag[i], hnl[i], + 14)); + } + + if (aecm->cngMode == AecmTrue) + { + WebRtcAecm_ComfortNoise(aecm, ptrDfaClean, efwReal, efwImag, hnl); + } + +#ifdef ARM_WINM_LOG_ + // measure tick end + QueryPerformanceCounter((LARGE_INTEGER*)&end); + diff__ = ((end - start) * 1000) / (freq/1000); + milliseconds = (unsigned int)(diff__ & 0xffffffff); + WriteFile (logFile, &milliseconds, sizeof(unsigned int), &temp, NULL); +#endif + +#ifdef ARM_WINM_LOG_ + // measure tick start + QueryPerformanceCounter((LARGE_INTEGER*)&start); +#endif + + // Synthesis + for (i = 1; i < PART_LEN; i++) + { + j = WEBRTC_SPL_LSHIFT_W32(i, 1); + fft[j] = efwReal[i]; + + // mirrored data, even + fft[PART_LEN4 - j] = efwReal[i]; + fft[j + 1] = -efwImag[i]; + + //mirrored data, odd + fft[PART_LEN4 - (j - 1)] = efwImag[i]; + } + fft[0] = efwReal[0]; + fft[1] = -efwImag[0]; + + fft[PART_LEN2] = efwReal[PART_LEN]; + fft[PART_LEN2 + 1] = -efwImag[PART_LEN]; + +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + // inverse FFT, result should be scaled with outCFFT + WebRtcSpl_ComplexBitReverse(fft, PART_LEN_SHIFT); + outCFFT = WebRtcSpl_ComplexIFFT(fft, PART_LEN_SHIFT, 1); + + //take only the real values and scale with outCFFT + for (i = 0; i < PART_LEN2; i++) + { + j = WEBRTC_SPL_LSHIFT_W32(i, 1); + fft[i] = fft[j]; + } +#else + outCFFT = WebRtcSpl_ComplexIFFT2(fft, postFft, PART_LEN_SHIFT, 1); + + //take only the real values and scale with outCFFT + for(i = 0, j = 0; i < PART_LEN2;) + { + fft[i] = postFft[j]; + i += 1; + j += 2; + fft[i] = postFft[j]; + i += 1; + j += 2; + } +#endif + + for (i = 0; i < PART_LEN; i++) + { + fft[i] = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT_WITH_ROUND( + fft[i], + kSqrtHanning[i], + 14); + tmp32no1 = WEBRTC_SPL_SHIFT_W32((WebRtc_Word32)fft[i], + outCFFT - aecm->dfaCleanQDomain); + fft[i] = (WebRtc_Word16)WEBRTC_SPL_SAT(WEBRTC_SPL_WORD16_MAX, + tmp32no1 + aecm->outBuf[i], + WEBRTC_SPL_WORD16_MIN); + output[i] = fft[i]; + + tmp32no1 = WEBRTC_SPL_MUL_16_16_RSFT( + fft[PART_LEN + i], + kSqrtHanning[PART_LEN - i], + 14); + tmp32no1 = WEBRTC_SPL_SHIFT_W32(tmp32no1, + outCFFT - aecm->dfaCleanQDomain); + aecm->outBuf[i] = (WebRtc_Word16)WEBRTC_SPL_SAT( + WEBRTC_SPL_WORD16_MAX, + tmp32no1, + WEBRTC_SPL_WORD16_MIN); + } + +#ifdef ARM_WINM_LOG_ + // measure tick end + QueryPerformanceCounter((LARGE_INTEGER*)&end); + diff__ = ((end - start) * 1000) / (freq/1000); + milliseconds = (unsigned int)(diff__ & 0xffffffff); + WriteFile (logFile, &milliseconds, sizeof(unsigned int), &temp, NULL); +#endif + // Copy the current block to the old position (outBuf is shifted elsewhere) + memcpy(aecm->xBuf, aecm->xBuf + PART_LEN, sizeof(WebRtc_Word16) * PART_LEN); + memcpy(aecm->dBufNoisy, aecm->dBufNoisy + PART_LEN, sizeof(WebRtc_Word16) * PART_LEN); + if (nearendClean != NULL) + { + memcpy(aecm->dBufClean, aecm->dBufClean + PART_LEN, sizeof(WebRtc_Word16) * PART_LEN); + } +} + +// Generate comfort noise and add to output signal. +// +// \param[in] aecm Handle of the AECM instance. +// \param[in] dfa Absolute value of the nearend signal (Q[aecm->dfaQDomain]). +// \param[in,out] outReal Real part of the output signal (Q[aecm->dfaQDomain]). +// \param[in,out] outImag Imaginary part of the output signal (Q[aecm->dfaQDomain]). +// \param[in] lambda Suppression gain with which to scale the noise level (Q14). +// +static void WebRtcAecm_ComfortNoise(AecmCore_t * const aecm, const WebRtc_UWord16 * const dfa, + WebRtc_Word16 * const outReal, + WebRtc_Word16 * const outImag, + const WebRtc_Word16 * const lambda) +{ + WebRtc_Word16 i; + WebRtc_Word16 tmp16; + WebRtc_Word32 tmp32; + + WebRtc_Word16 randW16[PART_LEN]; + WebRtc_Word16 uReal[PART_LEN1]; + WebRtc_Word16 uImag[PART_LEN1]; + WebRtc_Word32 outLShift32[PART_LEN1]; + WebRtc_Word16 noiseRShift16[PART_LEN1]; + + WebRtc_Word16 shiftFromNearToNoise[PART_LEN1]; + WebRtc_Word16 minTrackShift; + WebRtc_Word32 upper32; + WebRtc_Word32 lower32; + + if (aecm->noiseEstCtr < 100) + { + // Track the minimum more quickly initially. + aecm->noiseEstCtr++; + minTrackShift = 7; + } else + { + minTrackShift = 9; + } + + // Estimate noise power. + for (i = 0; i < PART_LEN1; i++) + { + shiftFromNearToNoise[i] = aecm->noiseEstQDomain[i] - aecm->dfaCleanQDomain; + + // Shift to the noise domain. + tmp32 = (WebRtc_Word32)dfa[i]; + outLShift32[i] = WEBRTC_SPL_SHIFT_W32(tmp32, shiftFromNearToNoise[i]); + + if (outLShift32[i] < aecm->noiseEst[i]) + { + // Track the minimum. + aecm->noiseEst[i] += ((outLShift32[i] - aecm->noiseEst[i]) >> minTrackShift); + } else + { + // Ramp slowly upwards until we hit the minimum again. + + // Avoid overflow. + if (aecm->noiseEst[i] < 2146435583) + { + // Store the fractional portion. + upper32 = (aecm->noiseEst[i] & 0xffff0000) >> 16; + lower32 = aecm->noiseEst[i] & 0x0000ffff; + upper32 = ((upper32 * 2049) >> 11); + lower32 = ((lower32 * 2049) >> 11); + aecm->noiseEst[i] = WEBRTC_SPL_ADD_SAT_W32(upper32 << 16, lower32); + } + } + } + + for (i = 0; i < PART_LEN1; i++) + { + tmp32 = WEBRTC_SPL_SHIFT_W32(aecm->noiseEst[i], -shiftFromNearToNoise[i]); + if (tmp32 > 32767) + { + tmp32 = 32767; + aecm->noiseEst[i] = WEBRTC_SPL_SHIFT_W32(tmp32, shiftFromNearToNoise[i]); + } + noiseRShift16[i] = (WebRtc_Word16)tmp32; + + tmp16 = ONE_Q14 - lambda[i]; + noiseRShift16[i] + = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(tmp16, noiseRShift16[i], 14); + } + + // Generate a uniform random array on [0 2^15-1]. + WebRtcSpl_RandUArray(randW16, PART_LEN, &aecm->seed); + + // Generate noise according to estimated energy. + uReal[0] = 0; // Reject LF noise. + uImag[0] = 0; + for (i = 1; i < PART_LEN1; i++) + { + // Get a random index for the cos and sin tables over [0 359]. + tmp16 = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(359, randW16[i - 1], 15); + + // Tables are in Q13. + uReal[i] = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(noiseRShift16[i], + WebRtcSpl_kCosTable[tmp16], 13); + uImag[i] = (WebRtc_Word16)WEBRTC_SPL_MUL_16_16_RSFT(-noiseRShift16[i], + WebRtcSpl_kSinTable[tmp16], 13); + } + uImag[PART_LEN] = 0; + +#if (!defined ARM_WINM) && (!defined ARM9E_GCC) && (!defined ANDROID_AECOPT) + for (i = 0; i < PART_LEN1; i++) + { + outReal[i] = WEBRTC_SPL_ADD_SAT_W16(outReal[i], uReal[i]); + outImag[i] = WEBRTC_SPL_ADD_SAT_W16(outImag[i], uImag[i]); + } +#else + for (i = 0; i < PART_LEN1 -1; ) + { + outReal[i] = WEBRTC_SPL_ADD_SAT_W16(outReal[i], uReal[i]); + outImag[i] = WEBRTC_SPL_ADD_SAT_W16(outImag[i], uImag[i]); + i++; + + outReal[i] = WEBRTC_SPL_ADD_SAT_W16(outReal[i], uReal[i]); + outImag[i] = WEBRTC_SPL_ADD_SAT_W16(outImag[i], uImag[i]); + i++; + } + outReal[i] = WEBRTC_SPL_ADD_SAT_W16(outReal[i], uReal[i]); + outImag[i] = WEBRTC_SPL_ADD_SAT_W16(outImag[i], uImag[i]); +#endif +} + +void WebRtcAecm_BufferFarFrame(AecmCore_t * const aecm, const WebRtc_Word16 * const farend, + const int farLen) +{ + int writeLen = farLen, writePos = 0; + + // Check if the write position must be wrapped + while (aecm->farBufWritePos + writeLen > FAR_BUF_LEN) + { + // Write to remaining buffer space before wrapping + writeLen = FAR_BUF_LEN - aecm->farBufWritePos; + memcpy(aecm->farBuf + aecm->farBufWritePos, farend + writePos, + sizeof(WebRtc_Word16) * writeLen); + aecm->farBufWritePos = 0; + writePos = writeLen; + writeLen = farLen - writeLen; + } + + memcpy(aecm->farBuf + aecm->farBufWritePos, farend + writePos, + sizeof(WebRtc_Word16) * writeLen); + aecm->farBufWritePos += writeLen; +} + +void WebRtcAecm_FetchFarFrame(AecmCore_t * const aecm, WebRtc_Word16 * const farend, + const int farLen, const int knownDelay) +{ + int readLen = farLen; + int readPos = 0; + int delayChange = knownDelay - aecm->lastKnownDelay; + + aecm->farBufReadPos -= delayChange; + + // Check if delay forces a read position wrap + while (aecm->farBufReadPos < 0) + { + aecm->farBufReadPos += FAR_BUF_LEN; + } + while (aecm->farBufReadPos > FAR_BUF_LEN - 1) + { + aecm->farBufReadPos -= FAR_BUF_LEN; + } + + aecm->lastKnownDelay = knownDelay; + + // Check if read position must be wrapped + while (aecm->farBufReadPos + readLen > FAR_BUF_LEN) + { + + // Read from remaining buffer space before wrapping + readLen = FAR_BUF_LEN - aecm->farBufReadPos; + memcpy(farend + readPos, aecm->farBuf + aecm->farBufReadPos, + sizeof(WebRtc_Word16) * readLen); + aecm->farBufReadPos = 0; + readPos = readLen; + readLen = farLen - readLen; + } + memcpy(farend + readPos, aecm->farBuf + aecm->farBufReadPos, + sizeof(WebRtc_Word16) * readLen); + aecm->farBufReadPos += readLen; +} diff --git a/src/modules/audio_processing/aecm/main/source/aecm_core.h b/src/modules/audio_processing/aecm/main/source/aecm_core.h new file mode 100644 index 0000000000..5defbe46c7 --- /dev/null +++ b/src/modules/audio_processing/aecm/main/source/aecm_core.h @@ -0,0 +1,338 @@ +/* + * 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 echo control (suppression) with fft routines in fixed-point + +#ifndef WEBRTC_MODULES_AUDIO_PROCESSING_AECM_MAIN_SOURCE_AECM_CORE_H_ +#define WEBRTC_MODULES_AUDIO_PROCESSING_AECM_MAIN_SOURCE_AECM_CORE_H_ + +#define AECM_DYNAMIC_Q // turn on/off dynamic Q-domain +//#define AECM_WITH_ABS_APPROX +//#define AECM_SHORT // for 32 sample partition length (otherwise 64) + +// TODO(bjornv): These defines will be removed in final version. +//#define STORE_CHANNEL_DATA +//#define VAD_DATA + +#include "typedefs.h" +#include "signal_processing_library.h" +// TODO(bjornv): Will be removed in final version. +//#include <stdio.h> + +// Algorithm parameters + +#define FRAME_LEN 80 // Total frame length, 10 ms +#ifdef AECM_SHORT + +#define PART_LEN 32 // Length of partition +#define PART_LEN_SHIFT 6 // Length of (PART_LEN * 2) in base 2 + +#else + +#define PART_LEN 64 // Length of partition +#define PART_LEN_SHIFT 7 // Length of (PART_LEN * 2) in base 2 + +#endif + +#define PART_LEN1 (PART_LEN + 1) // Unique fft coefficients +#define PART_LEN2 (PART_LEN << 1) // Length of partition * 2 +#define PART_LEN4 (PART_LEN << 2) // Length of partition * 4 +#define FAR_BUF_LEN PART_LEN4 // Length of buffers +#define MAX_DELAY 100 + +// Counter parameters +#ifdef AECM_SHORT + +#define CONV_LEN 1024 // Convergence length used at startup +#else + +#define CONV_LEN 512 // Convergence length used at startup +#endif + +#define CONV_LEN2 (CONV_LEN << 1) // Convergence length * 2 used at startup +// Energy parameters +#define MAX_BUF_LEN 64 // History length of energy signals + +#define FAR_ENERGY_MIN 1025 // Lowest Far energy level: At least 2 in energy +#define FAR_ENERGY_DIFF 929 // Allowed difference between max and min + +#define ENERGY_DEV_OFFSET 0 // The energy error offset in Q8 +#define ENERGY_DEV_TOL 400 // The energy estimation tolerance in Q8 +#define FAR_ENERGY_VAD_REGION 230 // Far VAD tolerance region +// Stepsize parameters +#define MU_MIN 10 // Min stepsize 2^-MU_MIN (far end energy dependent) +#define MU_MAX 1 // Max stepsize 2^-MU_MAX (far end energy dependent) +#define MU_DIFF 9 // MU_MIN - MU_MAX +// Channel parameters +#define MIN_MSE_COUNT 20 // Min number of consecutive blocks with enough far end + // energy to compare channel estimates +#define MIN_MSE_DIFF 29 // The ratio between adapted and stored channel to + // accept a new storage (0.8 in Q-MSE_RESOLUTION) +#define MSE_RESOLUTION 5 // MSE parameter resolution +#define RESOLUTION_CHANNEL16 12 // W16 Channel in Q-RESOLUTION_CHANNEL16 +#define RESOLUTION_CHANNEL32 28 // W32 Channel in Q-RESOLUTION_CHANNEL +#define CHANNEL_VAD 16 // Minimum energy in frequency band to update channel +// Suppression gain parameters: SUPGAIN_ parameters in Q-(RESOLUTION_SUPGAIN) +#define RESOLUTION_SUPGAIN 8 // Channel in Q-(RESOLUTION_SUPGAIN) +#define SUPGAIN_DEFAULT (1 << RESOLUTION_SUPGAIN) // Default suppression gain +#define SUPGAIN_ERROR_PARAM_A 3072 // Estimation error parameter (Maximum gain) (8 in Q8) +#define SUPGAIN_ERROR_PARAM_B 1536 // Estimation error parameter (Gain before going down) +#define SUPGAIN_ERROR_PARAM_D SUPGAIN_DEFAULT // Estimation error parameter + // (Should be the same as Default) (1 in Q8) +#define SUPGAIN_EPC_DT 200 // = SUPGAIN_ERROR_PARAM_C * ENERGY_DEV_TOL +// Defines for "check delay estimation" +#define CORR_WIDTH 31 // Number of samples to correlate over. +#define CORR_MAX 16 // Maximum correlation offset +#define CORR_MAX_BUF 63 +#define CORR_DEV 4 +#define CORR_MAX_LEVEL 20 +#define CORR_MAX_LOW 4 +#define CORR_BUF_LEN (CORR_MAX << 1) + 1 +// Note that CORR_WIDTH + 2*CORR_MAX <= MAX_BUF_LEN + +#define ONE_Q14 (1 << 14) + +// NLP defines +#define NLP_COMP_LOW 3277 // 0.2 in Q14 +#define NLP_COMP_HIGH ONE_Q14 // 1 in Q14 + +typedef struct +{ + int farBufWritePos; + int farBufReadPos; + int knownDelay; + int lastKnownDelay; + int firstVAD; // Parameter to control poorly initialized channels + + void *farFrameBuf; + void *nearNoisyFrameBuf; + void *nearCleanFrameBuf; + void *outFrameBuf; + + WebRtc_Word16 xBuf[PART_LEN2]; // farend + WebRtc_Word16 dBufClean[PART_LEN2]; // nearend + WebRtc_Word16 dBufNoisy[PART_LEN2]; // nearend + WebRtc_Word16 outBuf[PART_LEN]; + + WebRtc_Word16 farBuf[FAR_BUF_LEN]; + + WebRtc_Word16 mult; + WebRtc_UWord32 seed; + + // Delay estimation variables + WebRtc_UWord16 medianYlogspec[PART_LEN1]; + WebRtc_UWord16 medianXlogspec[PART_LEN1]; + WebRtc_UWord16 medianBCount[MAX_DELAY]; + WebRtc_UWord16 xfaHistory[PART_LEN1][MAX_DELAY]; + WebRtc_Word16 delHistoryPos; + WebRtc_UWord32 bxHistory[MAX_DELAY]; + WebRtc_UWord16 currentDelay; + WebRtc_UWord16 previousDelay; + WebRtc_Word16 delayAdjust; + + WebRtc_Word16 nlpFlag; + WebRtc_Word16 fixedDelay; + + WebRtc_UWord32 totCount; + + WebRtc_Word16 xfaQDomainBuf[MAX_DELAY]; + WebRtc_Word16 dfaCleanQDomain; + WebRtc_Word16 dfaCleanQDomainOld; + WebRtc_Word16 dfaNoisyQDomain; + WebRtc_Word16 dfaNoisyQDomainOld; + + WebRtc_Word16 nearLogEnergy[MAX_BUF_LEN]; + WebRtc_Word16 farLogEnergy[MAX_BUF_LEN]; + WebRtc_Word16 echoAdaptLogEnergy[MAX_BUF_LEN]; + WebRtc_Word16 echoStoredLogEnergy[MAX_BUF_LEN]; + + WebRtc_Word16 channelAdapt16[PART_LEN1]; + WebRtc_Word32 channelAdapt32[PART_LEN1]; + WebRtc_Word16 channelStored[PART_LEN1]; + WebRtc_Word32 echoFilt[PART_LEN1]; + WebRtc_Word16 nearFilt[PART_LEN1]; + WebRtc_Word32 noiseEst[PART_LEN1]; + WebRtc_Word16 noiseEstQDomain[PART_LEN1]; + WebRtc_Word16 noiseEstCtr; + WebRtc_Word16 cngMode; + + WebRtc_Word32 mseAdaptOld; + WebRtc_Word32 mseStoredOld; + WebRtc_Word32 mseThreshold; + + WebRtc_Word16 farEnergyMin; + WebRtc_Word16 farEnergyMax; + WebRtc_Word16 farEnergyMaxMin; + WebRtc_Word16 farEnergyVAD; + WebRtc_Word16 farEnergyMSE; + WebRtc_Word16 currentVADValue; + WebRtc_Word16 vadUpdateCount; + + WebRtc_Word16 delayHistogram[MAX_DELAY]; + WebRtc_Word16 delayVadCount; + WebRtc_Word16 maxDelayHistIdx; + WebRtc_Word16 lastMinPos; + + WebRtc_Word16 startupState; + WebRtc_Word16 mseChannelCount; + WebRtc_Word16 delayCount; + WebRtc_Word16 newDelayCorrData; + WebRtc_Word16 lastDelayUpdateCount; + WebRtc_Word16 delayCorrelation[CORR_BUF_LEN]; + WebRtc_Word16 supGain; + WebRtc_Word16 supGainOld; + WebRtc_Word16 delayOffsetFlag; + + WebRtc_Word16 supGainErrParamA; + WebRtc_Word16 supGainErrParamD; + WebRtc_Word16 supGainErrParamDiffAB; + WebRtc_Word16 supGainErrParamDiffBD; + + // TODO(bjornv): Will be removed after final version has been committed. +#ifdef VAD_DATA + FILE *vad_file; + FILE *delay_file; + FILE *far_file; + FILE *far_cur_file; + FILE *far_min_file; + FILE *far_max_file; + FILE *far_vad_file; +#endif + + // TODO(bjornv): Will be removed after final version has been committed. +#ifdef STORE_CHANNEL_DATA + FILE *channel_file; + FILE *channel_file_init; +#endif + +#ifdef AEC_DEBUG + FILE *farFile; + FILE *nearFile; + FILE *outFile; +#endif +} AecmCore_t; + +/////////////////////////////////////////////////////////////////////////////////////////////// +// WebRtcAecm_CreateCore(...) +// +// Allocates the memory needed by the AECM. The memory needs to be +// initialized separately using the WebRtcAecm_InitCore() function. +// +// Input: +// - aecm : Instance that should be created +// +// Output: +// - aecm : Created instance +// +// Return value : 0 - Ok +// -1 - Error +// +int WebRtcAecm_CreateCore(AecmCore_t **aecm); + +/////////////////////////////////////////////////////////////////////////////////////////////// +// WebRtcAecm_InitCore(...) +// +// This function initializes the AECM instant created with WebRtcAecm_CreateCore(...) +// Input: +// - aecm : Pointer to the AECM instance +// - samplingFreq : Sampling Frequency +// +// Output: +// - aecm : Initialized instance +// +// Return value : 0 - Ok +// -1 - Error +// +int WebRtcAecm_InitCore(AecmCore_t * const aecm, int samplingFreq); + +/////////////////////////////////////////////////////////////////////////////////////////////// +// WebRtcAecm_FreeCore(...) +// +// This function releases the memory allocated by WebRtcAecm_CreateCore() +// Input: +// - aecm : Pointer to the AECM instance +// +// Return value : 0 - Ok +// -1 - Error +// 11001-11016: Error +// +int WebRtcAecm_FreeCore(AecmCore_t *aecm); + +int WebRtcAecm_Control(AecmCore_t *aecm, int delay, int nlpFlag, int delayOffsetFlag); + +/////////////////////////////////////////////////////////////////////////////////////////////// +// WebRtcAecm_ProcessFrame(...) +// +// This function processes frames and sends blocks to WebRtcAecm_ProcessBlock(...) +// +// Inputs: +// - aecm : Pointer to the AECM instance +// - farend : In buffer containing one frame of echo signal +// - nearendNoisy : In buffer containing one frame of nearend+echo signal without NS +// - nearendClean : In buffer containing one frame of nearend+echo signal with NS +// +// Output: +// - out : Out buffer, one frame of nearend signal : +// +// +void WebRtcAecm_ProcessFrame(AecmCore_t * const aecm, const WebRtc_Word16 * const farend, + const WebRtc_Word16 * const nearendNoisy, + const WebRtc_Word16 * const nearendClean, + WebRtc_Word16 * const out); + +/////////////////////////////////////////////////////////////////////////////////////////////// +// WebRtcAecm_ProcessBlock(...) +// +// This function is called for every block within one frame +// This function is called by WebRtcAecm_ProcessFrame(...) +// +// Inputs: +// - aecm : Pointer to the AECM instance +// - farend : In buffer containing one block of echo signal +// - nearendNoisy : In buffer containing one frame of nearend+echo signal without NS +// - nearendClean : In buffer containing one frame of nearend+echo signal with NS +// +// Output: +// - out : Out buffer, one block of nearend signal : +// +// +void WebRtcAecm_ProcessBlock(AecmCore_t * const aecm, const WebRtc_Word16 * const farend, + const WebRtc_Word16 * const nearendNoisy, + const WebRtc_Word16 * const noisyClean, + WebRtc_Word16 * const out); + +/////////////////////////////////////////////////////////////////////////////////////////////// +// WebRtcAecm_BufferFarFrame() +// +// Inserts a frame of data into farend buffer. +// +// Inputs: +// - aecm : Pointer to the AECM instance +// - farend : In buffer containing one frame of farend signal +// - farLen : Length of frame +// +void WebRtcAecm_BufferFarFrame(AecmCore_t * const aecm, const WebRtc_Word16 * const farend, + const int farLen); + +/////////////////////////////////////////////////////////////////////////////////////////////// +// WebRtcAecm_FetchFarFrame() +// +// Read the farend buffer to account for known delay +// +// Inputs: +// - aecm : Pointer to the AECM instance +// - farend : In buffer containing one frame of farend signal +// - farLen : Length of frame +// - knownDelay : known delay +// +void WebRtcAecm_FetchFarFrame(AecmCore_t * const aecm, WebRtc_Word16 * const farend, + const int farLen, const int knownDelay); + +#endif diff --git a/src/modules/audio_processing/aecm/main/source/echo_control_mobile.c b/src/modules/audio_processing/aecm/main/source/echo_control_mobile.c new file mode 100644 index 0000000000..f9d84f0c4b --- /dev/null +++ b/src/modules/audio_processing/aecm/main/source/echo_control_mobile.c @@ -0,0 +1,733 @@ +/* + * 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 <stdlib.h> +//#include <string.h> + +#include "echo_control_mobile.h" +#include "aecm_core.h" +#include "ring_buffer.h" +#ifdef AEC_DEBUG +#include <stdio.h> +#endif +#ifdef MAC_IPHONE_PRINT +#include <time.h> +#include <stdio.h> +#elif defined ARM_WINM_LOG +#include "windows.h" +extern HANDLE logFile; +#endif + +#define BUF_SIZE_FRAMES 50 // buffer size (frames) +// Maximum length of resampled signal. Must be an integer multiple of frames +// (ceil(1/(1 + MIN_SKEW)*2) + 1)*FRAME_LEN +// The factor of 2 handles wb, and the + 1 is as a safety margin +#define MAX_RESAMP_LEN (5 * FRAME_LEN) + +static const int kBufSizeSamp = BUF_SIZE_FRAMES * FRAME_LEN; // buffer size (samples) +static const int kSampMsNb = 8; // samples per ms in nb +// Target suppression levels for nlp modes +// log{0.001, 0.00001, 0.00000001} +static const int kInitCheck = 42; + +typedef struct +{ + int sampFreq; + int scSampFreq; + short bufSizeStart; + int knownDelay; + + // Stores the last frame added to the farend buffer + short farendOld[2][FRAME_LEN]; + short initFlag; // indicates if AEC has been initialized + + // Variables used for averaging far end buffer size + short counter; + short sum; + short firstVal; + short checkBufSizeCtr; + + // Variables used for delay shifts + short msInSndCardBuf; + short filtDelay; + int timeForDelayChange; + int ECstartup; + int checkBuffSize; + int delayChange; + short lastDelayDiff; + + WebRtc_Word16 echoMode; + +#ifdef AEC_DEBUG + FILE *bufFile; + FILE *delayFile; + FILE *preCompFile; + FILE *postCompFile; +#endif // AEC_DEBUG + // Structures + void *farendBuf; + + int lastError; + + AecmCore_t *aecmCore; +} aecmob_t; + +// Estimates delay to set the position of the farend buffer read pointer +// (controlled by knownDelay) +static int WebRtcAecm_EstBufDelay(aecmob_t *aecmInst, short msInSndCardBuf); + +// Stuffs the farend buffer if the estimated delay is too large +static int WebRtcAecm_DelayComp(aecmob_t *aecmInst); + +WebRtc_Word32 WebRtcAecm_Create(void **aecmInst) +{ + aecmob_t *aecm; + if (aecmInst == NULL) + { + return -1; + } + + aecm = malloc(sizeof(aecmob_t)); + *aecmInst = aecm; + if (aecm == NULL) + { + return -1; + } + + if (WebRtcAecm_CreateCore(&aecm->aecmCore) == -1) + { + WebRtcAecm_Free(aecm); + aecm = NULL; + return -1; + } + + if (WebRtcApm_CreateBuffer(&aecm->farendBuf, kBufSizeSamp) == -1) + { + WebRtcAecm_Free(aecm); + aecm = NULL; + return -1; + } + + aecm->initFlag = 0; + aecm->lastError = 0; + +#ifdef AEC_DEBUG + aecm->aecmCore->farFile = fopen("aecFar.pcm","wb"); + aecm->aecmCore->nearFile = fopen("aecNear.pcm","wb"); + aecm->aecmCore->outFile = fopen("aecOut.pcm","wb"); + //aecm->aecmCore->outLpFile = fopen("aecOutLp.pcm","wb"); + + aecm->bufFile = fopen("aecBuf.dat", "wb"); + aecm->delayFile = fopen("aecDelay.dat", "wb"); + aecm->preCompFile = fopen("preComp.pcm", "wb"); + aecm->postCompFile = fopen("postComp.pcm", "wb"); +#endif // AEC_DEBUG + return 0; +} + +WebRtc_Word32 WebRtcAecm_Free(void *aecmInst) +{ + aecmob_t *aecm = aecmInst; + + if (aecm == NULL) + { + return -1; + } + +#ifdef AEC_DEBUG + fclose(aecm->aecmCore->farFile); + fclose(aecm->aecmCore->nearFile); + fclose(aecm->aecmCore->outFile); + //fclose(aecm->aecmCore->outLpFile); + + fclose(aecm->bufFile); + fclose(aecm->delayFile); + fclose(aecm->preCompFile); + fclose(aecm->postCompFile); +#endif // AEC_DEBUG + WebRtcAecm_FreeCore(aecm->aecmCore); + WebRtcApm_FreeBuffer(aecm->farendBuf); + free(aecm); + + return 0; +} + +WebRtc_Word32 WebRtcAecm_Init(void *aecmInst, WebRtc_Word32 sampFreq, WebRtc_Word32 scSampFreq) +{ + aecmob_t *aecm = aecmInst; + AecmConfig aecConfig; + + if (aecm == NULL) + { + return -1; + } + + if (sampFreq != 8000 && sampFreq != 16000) + { + aecm->lastError = AECM_BAD_PARAMETER_ERROR; + return -1; + } + aecm->sampFreq = sampFreq; + + if (scSampFreq < 1 || scSampFreq > 96000) + { + aecm->lastError = AECM_BAD_PARAMETER_ERROR; + return -1; + } + aecm->scSampFreq = scSampFreq; + + // Initialize AECM core + if (WebRtcAecm_InitCore(aecm->aecmCore, aecm->sampFreq) == -1) + { + aecm->lastError = AECM_UNSPECIFIED_ERROR; + return -1; + } + + // Initialize farend buffer + if (WebRtcApm_InitBuffer(aecm->farendBuf) == -1) + { + aecm->lastError = AECM_UNSPECIFIED_ERROR; + return -1; + } + + aecm->initFlag = kInitCheck; // indicates that initialization has been done + + aecm->delayChange = 1; + + aecm->sum = 0; + aecm->counter = 0; + aecm->checkBuffSize = 1; + aecm->firstVal = 0; + + aecm->ECstartup = 1; + aecm->bufSizeStart = 0; + aecm->checkBufSizeCtr = 0; + aecm->filtDelay = 0; + aecm->timeForDelayChange = 0; + aecm->knownDelay = 0; + aecm->lastDelayDiff = 0; + + memset(&aecm->farendOld[0][0], 0, 160); + + // Default settings. + aecConfig.cngMode = AecmTrue; + aecConfig.echoMode = 3; + + if (WebRtcAecm_set_config(aecm, aecConfig) == -1) + { + aecm->lastError = AECM_UNSPECIFIED_ERROR; + return -1; + } + + return 0; +} + +WebRtc_Word32 WebRtcAecm_BufferFarend(void *aecmInst, const WebRtc_Word16 *farend, + WebRtc_Word16 nrOfSamples) +{ + aecmob_t *aecm = aecmInst; + WebRtc_Word32 retVal = 0; + + if (aecm == NULL) + { + return -1; + } + + if (farend == NULL) + { + aecm->lastError = AECM_NULL_POINTER_ERROR; + return -1; + } + + if (aecm->initFlag != kInitCheck) + { + aecm->lastError = AECM_UNINITIALIZED_ERROR; + return -1; + } + + if (nrOfSamples != 80 && nrOfSamples != 160) + { + aecm->lastError = AECM_BAD_PARAMETER_ERROR; + return -1; + } + + // TODO: Is this really a good idea? + if (!aecm->ECstartup) + { + WebRtcAecm_DelayComp(aecm); + } + + WebRtcApm_WriteBuffer(aecm->farendBuf, farend, nrOfSamples); + + return retVal; +} + +WebRtc_Word32 WebRtcAecm_Process(void *aecmInst, const WebRtc_Word16 *nearendNoisy, + const WebRtc_Word16 *nearendClean, WebRtc_Word16 *out, + WebRtc_Word16 nrOfSamples, WebRtc_Word16 msInSndCardBuf) +{ + aecmob_t *aecm = aecmInst; + WebRtc_Word32 retVal = 0; + short i; + short farend[FRAME_LEN]; + short nmbrOfFilledBuffers; + short nBlocks10ms; + short nFrames; +#ifdef AEC_DEBUG + short msInAECBuf; +#endif + +#ifdef ARM_WINM_LOG + __int64 freq, start, end, diff; + unsigned int milliseconds; + DWORD temp; +#elif defined MAC_IPHONE_PRINT + // double endtime = 0, starttime = 0; + struct timeval starttime; + struct timeval endtime; + static long int timeused = 0; + static int timecount = 0; +#endif + + if (aecm == NULL) + { + return -1; + } + + if (nearendNoisy == NULL) + { + aecm->lastError = AECM_NULL_POINTER_ERROR; + return -1; + } + + if (out == NULL) + { + aecm->lastError = AECM_NULL_POINTER_ERROR; + return -1; + } + + if (aecm->initFlag != kInitCheck) + { + aecm->lastError = AECM_UNINITIALIZED_ERROR; + return -1; + } + + if (nrOfSamples != 80 && nrOfSamples != 160) + { + aecm->lastError = AECM_BAD_PARAMETER_ERROR; + return -1; + } + + if (msInSndCardBuf < 0) + { + msInSndCardBuf = 0; + aecm->lastError = AECM_BAD_PARAMETER_WARNING; + retVal = -1; + } else if (msInSndCardBuf > 500) + { + msInSndCardBuf = 500; + aecm->lastError = AECM_BAD_PARAMETER_WARNING; + retVal = -1; + } + msInSndCardBuf += 10; + aecm->msInSndCardBuf = msInSndCardBuf; + + nFrames = nrOfSamples / FRAME_LEN; + nBlocks10ms = nFrames / aecm->aecmCore->mult; + + if (aecm->ECstartup) + { + if (nearendClean == NULL) + { + memcpy(out, nearendNoisy, sizeof(short) * nrOfSamples); + } else + { + memcpy(out, nearendClean, sizeof(short) * nrOfSamples); + } + + nmbrOfFilledBuffers = WebRtcApm_get_buffer_size(aecm->farendBuf) / FRAME_LEN; + // The AECM is in the start up mode + // AECM is disabled until the soundcard buffer and farend buffers are OK + + // Mechanism to ensure that the soundcard buffer is reasonably stable. + if (aecm->checkBuffSize) + { + aecm->checkBufSizeCtr++; + // Before we fill up the far end buffer we require the amount of data on the + // sound card to be stable (+/-8 ms) compared to the first value. This + // comparison is made during the following 4 consecutive frames. If it seems + // to be stable then we start to fill up the far end buffer. + + if (aecm->counter == 0) + { + aecm->firstVal = aecm->msInSndCardBuf; + aecm->sum = 0; + } + + if (abs(aecm->firstVal - aecm->msInSndCardBuf) + < WEBRTC_SPL_MAX(0.2 * aecm->msInSndCardBuf, kSampMsNb)) + { + aecm->sum += aecm->msInSndCardBuf; + aecm->counter++; + } else + { + aecm->counter = 0; + } + + if (aecm->counter * nBlocks10ms >= 6) + { + // The farend buffer size is determined in blocks of 80 samples + // Use 75% of the average value of the soundcard buffer + aecm->bufSizeStart + = WEBRTC_SPL_MIN((3 * aecm->sum + * aecm->aecmCore->mult) / (aecm->counter * 40), BUF_SIZE_FRAMES); + // buffersize has now been determined + aecm->checkBuffSize = 0; + } + + if (aecm->checkBufSizeCtr * nBlocks10ms > 50) + { + // for really bad sound cards, don't disable echocanceller for more than 0.5 sec + aecm->bufSizeStart = WEBRTC_SPL_MIN((3 * aecm->msInSndCardBuf + * aecm->aecmCore->mult) / 40, BUF_SIZE_FRAMES); + aecm->checkBuffSize = 0; + } + } + + // if checkBuffSize changed in the if-statement above + if (!aecm->checkBuffSize) + { + // soundcard buffer is now reasonably stable + // When the far end buffer is filled with approximately the same amount of + // data as the amount on the sound card we end the start up phase and start + // to cancel echoes. + + if (nmbrOfFilledBuffers == aecm->bufSizeStart) + { + aecm->ECstartup = 0; // Enable the AECM + } else if (nmbrOfFilledBuffers > aecm->bufSizeStart) + { + WebRtcApm_FlushBuffer( + aecm->farendBuf, + WebRtcApm_get_buffer_size(aecm->farendBuf) + - aecm->bufSizeStart * FRAME_LEN); + aecm->ECstartup = 0; + } + } + + } else + { + // AECM is enabled + + // Note only 1 block supported for nb and 2 blocks for wb + for (i = 0; i < nFrames; i++) + { + nmbrOfFilledBuffers = WebRtcApm_get_buffer_size(aecm->farendBuf) / FRAME_LEN; + + // Check that there is data in the far end buffer + if (nmbrOfFilledBuffers > 0) + { + // Get the next 80 samples from the farend buffer + WebRtcApm_ReadBuffer(aecm->farendBuf, farend, FRAME_LEN); + + // Always store the last frame for use when we run out of data + memcpy(&(aecm->farendOld[i][0]), farend, FRAME_LEN * sizeof(short)); + } else + { + // We have no data so we use the last played frame + memcpy(farend, &(aecm->farendOld[i][0]), FRAME_LEN * sizeof(short)); + } + + // Call buffer delay estimator when all data is extracted, + // i,e. i = 0 for NB and i = 1 for WB + if ((i == 0 && aecm->sampFreq == 8000) || (i == 1 && aecm->sampFreq == 16000)) + { + WebRtcAecm_EstBufDelay(aecm, aecm->msInSndCardBuf); + } + +#ifdef ARM_WINM_LOG + // measure tick start + QueryPerformanceFrequency((LARGE_INTEGER*)&freq); + QueryPerformanceCounter((LARGE_INTEGER*)&start); +#elif defined MAC_IPHONE_PRINT + // starttime = clock()/(double)CLOCKS_PER_SEC; + gettimeofday(&starttime, NULL); +#endif + // Call the AECM + /*WebRtcAecm_ProcessFrame(aecm->aecmCore, farend, &nearend[FRAME_LEN * i], + &out[FRAME_LEN * i], aecm->knownDelay);*/ + if (nearendClean == NULL) + { + WebRtcAecm_ProcessFrame(aecm->aecmCore, farend, &nearendNoisy[FRAME_LEN * i], + NULL, &out[FRAME_LEN * i]); + } else + { + WebRtcAecm_ProcessFrame(aecm->aecmCore, farend, &nearendNoisy[FRAME_LEN * i], + &nearendClean[FRAME_LEN * i], &out[FRAME_LEN * i]); + } + +#ifdef ARM_WINM_LOG + + // measure tick end + QueryPerformanceCounter((LARGE_INTEGER*)&end); + + if(end > start) + { + diff = ((end - start) * 1000) / (freq/1000); + milliseconds = (unsigned int)(diff & 0xffffffff); + WriteFile (logFile, &milliseconds, sizeof(unsigned int), &temp, NULL); + } +#elif defined MAC_IPHONE_PRINT + // endtime = clock()/(double)CLOCKS_PER_SEC; + // printf("%f\n", endtime - starttime); + + gettimeofday(&endtime, NULL); + + if( endtime.tv_usec > starttime.tv_usec) + { + timeused += endtime.tv_usec - starttime.tv_usec; + } else + { + timeused += endtime.tv_usec + 1000000 - starttime.tv_usec; + } + + if(++timecount == 1000) + { + timecount = 0; + printf("AEC: %ld\n", timeused); + timeused = 0; + } +#endif + + } + } + +#ifdef AEC_DEBUG + msInAECBuf = WebRtcApm_get_buffer_size(aecm->farendBuf) / (kSampMsNb*aecm->aecmCore->mult); + fwrite(&msInAECBuf, 2, 1, aecm->bufFile); + fwrite(&(aecm->knownDelay), sizeof(aecm->knownDelay), 1, aecm->delayFile); +#endif + + return retVal; +} + +WebRtc_Word32 WebRtcAecm_set_config(void *aecmInst, AecmConfig config) +{ + aecmob_t *aecm = aecmInst; + + if (aecm == NULL) + { + return -1; + } + + if (aecm->initFlag != kInitCheck) + { + aecm->lastError = AECM_UNINITIALIZED_ERROR; + return -1; + } + + if (config.cngMode != AecmFalse && config.cngMode != AecmTrue) + { + aecm->lastError = AECM_BAD_PARAMETER_ERROR; + return -1; + } + aecm->aecmCore->cngMode = config.cngMode; + + if (config.echoMode < 0 || config.echoMode > 4) + { + aecm->lastError = AECM_BAD_PARAMETER_ERROR; + return -1; + } + aecm->echoMode = config.echoMode; + + if (aecm->echoMode == 0) + { + aecm->aecmCore->supGain = SUPGAIN_DEFAULT >> 3; + aecm->aecmCore->supGainOld = SUPGAIN_DEFAULT >> 3; + aecm->aecmCore->supGainErrParamA = SUPGAIN_ERROR_PARAM_A >> 3; + aecm->aecmCore->supGainErrParamD = SUPGAIN_ERROR_PARAM_D >> 3; + aecm->aecmCore->supGainErrParamDiffAB = (SUPGAIN_ERROR_PARAM_A >> 3) + - (SUPGAIN_ERROR_PARAM_B >> 3); + aecm->aecmCore->supGainErrParamDiffBD = (SUPGAIN_ERROR_PARAM_B >> 3) + - (SUPGAIN_ERROR_PARAM_D >> 3); + } else if (aecm->echoMode == 1) + { + aecm->aecmCore->supGain = SUPGAIN_DEFAULT >> 2; + aecm->aecmCore->supGainOld = SUPGAIN_DEFAULT >> 2; + aecm->aecmCore->supGainErrParamA = SUPGAIN_ERROR_PARAM_A >> 2; + aecm->aecmCore->supGainErrParamD = SUPGAIN_ERROR_PARAM_D >> 2; + aecm->aecmCore->supGainErrParamDiffAB = (SUPGAIN_ERROR_PARAM_A >> 2) + - (SUPGAIN_ERROR_PARAM_B >> 2); + aecm->aecmCore->supGainErrParamDiffBD = (SUPGAIN_ERROR_PARAM_B >> 2) + - (SUPGAIN_ERROR_PARAM_D >> 2); + } else if (aecm->echoMode == 2) + { + aecm->aecmCore->supGain = SUPGAIN_DEFAULT >> 1; + aecm->aecmCore->supGainOld = SUPGAIN_DEFAULT >> 1; + aecm->aecmCore->supGainErrParamA = SUPGAIN_ERROR_PARAM_A >> 1; + aecm->aecmCore->supGainErrParamD = SUPGAIN_ERROR_PARAM_D >> 1; + aecm->aecmCore->supGainErrParamDiffAB = (SUPGAIN_ERROR_PARAM_A >> 1) + - (SUPGAIN_ERROR_PARAM_B >> 1); + aecm->aecmCore->supGainErrParamDiffBD = (SUPGAIN_ERROR_PARAM_B >> 1) + - (SUPGAIN_ERROR_PARAM_D >> 1); + } else if (aecm->echoMode == 3) + { + aecm->aecmCore->supGain = SUPGAIN_DEFAULT; + aecm->aecmCore->supGainOld = SUPGAIN_DEFAULT; + aecm->aecmCore->supGainErrParamA = SUPGAIN_ERROR_PARAM_A; + aecm->aecmCore->supGainErrParamD = SUPGAIN_ERROR_PARAM_D; + aecm->aecmCore->supGainErrParamDiffAB = SUPGAIN_ERROR_PARAM_A - SUPGAIN_ERROR_PARAM_B; + aecm->aecmCore->supGainErrParamDiffBD = SUPGAIN_ERROR_PARAM_B - SUPGAIN_ERROR_PARAM_D; + } else if (aecm->echoMode == 4) + { + aecm->aecmCore->supGain = SUPGAIN_DEFAULT << 1; + aecm->aecmCore->supGainOld = SUPGAIN_DEFAULT << 1; + aecm->aecmCore->supGainErrParamA = SUPGAIN_ERROR_PARAM_A << 1; + aecm->aecmCore->supGainErrParamD = SUPGAIN_ERROR_PARAM_D << 1; + aecm->aecmCore->supGainErrParamDiffAB = (SUPGAIN_ERROR_PARAM_A << 1) + - (SUPGAIN_ERROR_PARAM_B << 1); + aecm->aecmCore->supGainErrParamDiffBD = (SUPGAIN_ERROR_PARAM_B << 1) + - (SUPGAIN_ERROR_PARAM_D << 1); + } + + return 0; +} + +WebRtc_Word32 WebRtcAecm_get_config(void *aecmInst, AecmConfig *config) +{ + aecmob_t *aecm = aecmInst; + + if (aecm == NULL) + { + return -1; + } + + if (config == NULL) + { + aecm->lastError = AECM_NULL_POINTER_ERROR; + return -1; + } + + if (aecm->initFlag != kInitCheck) + { + aecm->lastError = AECM_UNINITIALIZED_ERROR; + return -1; + } + + config->cngMode = aecm->aecmCore->cngMode; + config->echoMode = aecm->echoMode; + + return 0; +} + +WebRtc_Word32 WebRtcAecm_get_version(WebRtc_Word8 *versionStr, WebRtc_Word16 len) +{ + const char version[] = "AECM 1.2.0"; + const short versionLen = (short)strlen(version) + 1; // +1 for null-termination + + if (versionStr == NULL) + { + return -1; + } + + if (versionLen > len) + { + return -1; + } + + strncpy(versionStr, version, versionLen); + return 0; +} + +WebRtc_Word32 WebRtcAecm_get_error_code(void *aecmInst) +{ + aecmob_t *aecm = aecmInst; + + if (aecm == NULL) + { + return -1; + } + + return aecm->lastError; +} + +static int WebRtcAecm_EstBufDelay(aecmob_t *aecm, short msInSndCardBuf) +{ + short delayNew, nSampFar, nSampSndCard; + short diff; + + nSampFar = WebRtcApm_get_buffer_size(aecm->farendBuf); + nSampSndCard = msInSndCardBuf * kSampMsNb * aecm->aecmCore->mult; + + delayNew = nSampSndCard - nSampFar; + + if (delayNew < FRAME_LEN) + { + WebRtcApm_FlushBuffer(aecm->farendBuf, FRAME_LEN); + delayNew += FRAME_LEN; + } + + aecm->filtDelay = WEBRTC_SPL_MAX(0, (8 * aecm->filtDelay + 2 * delayNew) / 10); + + diff = aecm->filtDelay - aecm->knownDelay; + if (diff > 224) + { + if (aecm->lastDelayDiff < 96) + { + aecm->timeForDelayChange = 0; + } else + { + aecm->timeForDelayChange++; + } + } else if (diff < 96 && aecm->knownDelay > 0) + { + if (aecm->lastDelayDiff > 224) + { + aecm->timeForDelayChange = 0; + } else + { + aecm->timeForDelayChange++; + } + } else + { + aecm->timeForDelayChange = 0; + } + aecm->lastDelayDiff = diff; + + if (aecm->timeForDelayChange > 25) + { + aecm->knownDelay = WEBRTC_SPL_MAX((int)aecm->filtDelay - 160, 0); + } + return 0; +} + +static int WebRtcAecm_DelayComp(aecmob_t *aecm) +{ + int nSampFar, nSampSndCard, delayNew, nSampAdd; + const int maxStuffSamp = 10 * FRAME_LEN; + + nSampFar = WebRtcApm_get_buffer_size(aecm->farendBuf); + nSampSndCard = aecm->msInSndCardBuf * kSampMsNb * aecm->aecmCore->mult; + delayNew = nSampSndCard - nSampFar; + + if (delayNew > FAR_BUF_LEN - FRAME_LEN * aecm->aecmCore->mult) + { + // The difference of the buffer sizes is larger than the maximum + // allowed known delay. Compensate by stuffing the buffer. + nSampAdd = (int)(WEBRTC_SPL_MAX(((nSampSndCard >> 1) - nSampFar), + FRAME_LEN)); + nSampAdd = WEBRTC_SPL_MIN(nSampAdd, maxStuffSamp); + + WebRtcApm_StuffBuffer(aecm->farendBuf, nSampAdd); + aecm->delayChange = 1; // the delay needs to be updated + } + + return 0; +} |