WebRTC的AEC算法主要包括以下几个重要的模块:
- 回声延迟估计模块
- NLMS算法(见前一章)
- NLP非线性滤波 4 .CNG(comfort nosie generation),我的github上将噪声产生功能去掉了,因为ASR原来就是需要纯净语音的. 此外,像speex之类的算法还包括双讲检测(DT,double talk) WebRTC MATLAB and c code github Address
- 回声延迟估计 回声延迟长短对回声抵消器的性能有较大影响，过长的的滤波器抽头会带来较大的延迟，并且语音信号是短时平稳信号，过长的滤波器抽头也不适合短时平稳信号特征。 基于相关的延迟算法。该算法的主要思想是： 设1表示有说话音，0表无说话声（或者很弱的说话声），参考端（远端）信号和接收端信号可能的组合方式如下： webrtc默认(1, 0)和（0，1）是不可能发生的。设在时间间隔p上，即,频带q,,输入信号加窗后的功率谱用表示，其角标表示其加了窗函数。对每个频带的功率谱设定一个门限,
如果, 则; 如果, 则; 同理,对于信号,加窗信号功率谱和门限如果, 则; 如果, 则; 考虑到实际处理的方便，在webrtc的c代码中，将经过fft变换后的频域功率谱分为32个子带，这样每个特定子带的值可以用1个比特来表示，总共需要32个比特，只用一个32位数据类型就可以表示。
- NLMS归一化最小均放只适应算法 见前一章
NLP(非线性滤波) webrtc采用了维纳滤波器，此处只给出传递函数的表达式，设估计的语音信号的功率谱为, 噪声的功率谱为, 则滤波器的传递函数为
对于webRTC的AEC需要注意两点:
延迟要小,因为算法默认滤波器长度是分为12块，每块64点，按照8000采样率，也就是12×8ms=96ms的数据，超过这个长度就处理不了了。
- 延迟抖动要小，因为算法是默认10块也计算一次参考数据的位置（即滤波器能量最大的那一块），所以如果抖动很大的话，找参考数据不准确的，这样回声就消除不掉了。
实现代码分析
- 几个名词
RERL-residual_echo_return_loss ERL-echo return loss ERLE echo return loss enhancement NLP non-linear processing 前三个是性能评估的参考标准.
初始化和数据读入
%near is micphone captured signal
fid=fopen('near.pcm', 'rb'); % Load far end
ssin=fread(fid,inf,'float32');
fclose(fid);
%far is speaker played music
fid=fopen('far.pcm', 'rb'); % Load fnear end
rrin=fread(fid,inf,'float32');
fclose(fid);
然后对一些变量赋初值
fs=16000;
NLPon=1; % NLP on
M = 16; % Number of partitions
N = 64; % Partition length
L = M*N; % Filter length
VADtd=48;
alp = 0.15; % Power estimation factor
alc = 0.1; % Coherence estimation factor
step = 0.1875;%0.1875; % Downward step size
上述初始化中，M=16和最新的WebRTC代码并不一致，且最新的WebRTC中支持aec3最新一代算法。
len=length(ssin);
NN=len;
Nb=floor(NN/N)-M;
for kk=1:Nb
pos = N * (kk-1) + start;
可以看出Nb是麦克风采集到的数据块数-16（分区数），这是因为第一次输入了16块，所以这里减掉了16。pos是每一次添加一块时的地址指针。
%far is speaker played music
xk = rrin(pos:pos+N-1);
%near is micphone captured signal
dk = ssin(pos:pos+N-1);
xk和dk是读取到的64个点，这里是时域信号。
功率计算
%----------------------- far end signal process
xx = [xo;xk];
xo = xk;
tmp = fft(xx);
XX = tmp(1:N+1);
dd = [do;dk]; % Overlap
do = dk;
tmp = fft(dd); % Frequency domain
DD = tmp(1:N+1);
将xk和上一次的数据结合在一起，做FFT变换，由于两次组合在一起，那么得到的应该是N=128点，这里可以知道得到的谱分辨率是$n*fs/N$,$fs$前面设置过了，是16k，则得到的每一个bin的频谱分辨率是16000/128=125Hz。这里XX和DD取了前65个点，这是因为N点FFT变换后频谱是关于N/2对称的，对称的原因是奈奎斯特采样定理，如果$fs=16000Hz$,那么要求采样到的信号的截止频率必然小于等于$fs/2=8000Hz$,对于实信号，N/2~N,实际上表示的是$-fs/2$ ~ $0$之间的频率。第一个点是直流分量，所以取65个点。和上一帧64个点信号合并在一起的另一个原因是使用重叠保（overlap-save）留法将循环卷积变成线性卷积，这里做的FFT变换，就是为了减少时域里做卷积的计算量的。 计算远端信号功率谱
% ------------------------far end Power estimation
pn0 = (1 - alp) * pn0 + alp * real(XX.* conj(XX));
pn = pn0;
平滑功率谱，上一次的功率谱占85%（alp=0.15），后面的频域共轭相乘等于功率是有帕斯瓦尔定理支撑的。pn0是65*1的矩阵。
滤波
XFm(:,1) = XX;
首先将远端信号频谱赋给XFm，XFm是65*16的矩阵，16就是前面初始化的M值，这里将XX给第一列，其2~16列对应的是之前的输入频谱。
for mm=0:(M-1)
m=mm+1;
YFb(:,m) = XFm(:,m) .* WFb(:,m);
end
YFb，WFb以及XFm都是65*16的矩阵，WFb是自适应滤波器的频谱表示，XFm是原始的speaker数据，上式的意义对应于插图中的$\hat{y}$的频域值，变换到时域后就可以得到$y$的估计值$\hat{y}$.
yfk = sum(YFb,2);
tmp = [yfk ; flipud(conj(yfk(2:N)))];
ykt = real(ifft(tmp));
ykfb = ykt(end-N+1:end);
首先yfk是65_1的矩阵，sum求和就是将估计的频谱按行求和，也就是yfk包含了最近16个块的远端频谱估计信息，这样，只要近端麦克采集到的信号里有这16个块包含的远端信号，那么就可以消掉，从这里也可以看出来，容许的延迟差 在16_64/16=64ms，也就是说，如果麦克风采集到的speaker信号滞后speaker播放超过64ms，那么这种情况是无法消掉的，当然，延迟差越小越好。 flipud(conj(yfk(2:N))是因为前面计算频谱时利用奈奎斯特定理，也即实数的FFT结果按N/2对称，所以这里为了得到正确的ifft变换结果，先把频谱不全到$fs$. ykfb就是.后面再看WFb是如何跟新。
误差估计
ekfb = dk - ykfb;
dk是麦克风采集到的信号，ykfb是$\hat{y}$,这样得到的是误差信号，理想情况下，那么得到的误差信号就是需要的人声信号，而完全滤除 掉了speaker信号（远端信号）。
erfb(pos:pos+N-1) = ekfb;
tmp = fft([zm;ekfb]); % FD version for cancelling part (overlap-save)
Ek = tmp(1:N+1);
erfb是近端信号数组长度×1矩阵，存放的是全部样本对应的误差信号，这个保存仅仅是为了plot用的。 然后补了64个零，然后做FFT，Ek是误差信号FFT的结果。
自适应调节
Ek2 = Ek ./(pn + 0.001); % Normalized error
pn是当前帧远端信号功率谱，Ek是误差信号频谱。Ek2是归一化误差频谱。NLMS公式要求。
absEf = max(abs(Ek2), threshold);
absEf = ones(N+1,1)*threshold./absEf;
Ek2 = Ek2.*absEf;
max的作用是为了防止归一化后误差频谱过小，最终得到的Ek2是一个限幅矩阵，如果该点的值比门限大，则取门限，如果该点的值比门限小，则保持不变。
mEk = mufb.*Ek2;
mufb是步长，对于16000情况，步长取了0.8.NLMS公式。
PP = conj(XFm).*(ones(M,1) * mEk')';
tmp = [PP ; flipud(conj(PP(2:N,:)))];
IFPP = real(ifft(tmp));
PH = IFPP(1:N,:);
tmp = fft([PH;zeros(N,M)]);
FPH = tmp(1:N+1,:);
WFb = WFb + FPH;
PP是将远端信号的共轭乘以误差信号频谱，这一项用于调节步长，NLMS（步长=参考信号×步长×误差）的可变步长就提现在这里。PH是频域到时域的变换值。这和前面频域到时域的变换原理一样。WFb是权中系数的更新。
if mod(kk, 10*mult) == 0
WFbEn = sum(real(WFb.*conj(WFb)));
%WFbEn = sum(abs(WFb));
[tmp, dIdx] = max(WFbEn);
WFbD = sum(abs(WFb(:, dIdx)),2);
%WFbD = WFbD / (mean(WFbD) + 1e-10);
WFbD = min(max(WFbD, 0.5), 4);
end
dIdxV(kk) = dIdx;
上述的作用是更新dIdx和dIdxV。这里的更新并不是每一次都更新，一来是为了稳定，而来也是变相的减少计算量，提高实时性。就算是每一次都更新dIdx，WebRTC计算速度评估的结果也是很满意的。WFb是权重向量的频谱表示，WFbEn是权重向量按列求和，得到的是16_1的矩阵。这样得到的是16个块对权重的累加和。这样的区分度比直接累加和要大。 [tmp, dIdx] = max(WFbEn);作用就是找到16个块中权重累加和最大值及其对应的索引。 WFbD首先计算了权重最大那个块dIdx的列，然后将其按行求和，得到的就是65_1矩阵。WFbD不能低于0.5也不能高于4，算法中并未使用到，plot性能分析时用到。 最后把索引值dIdx存放到dIdxV(kk)中，这样每来一帧，就会有一个最大索引值放到dIdxV向量中。
功率谱密度和相关性计算
NLP
这里的NLP(Non-linear processing)的意思。
ee = [eo;ekfb];
eo = ekfb;
window = wins;
上述作用是将上次的误差和ekfb组合，其中eo可以理解为error old。eo也确实保存了上一次的误差。window是简单将窗函数赋值。
tmp = fft(xx.*window);
xf = tmp(1:N+1);
tmp = fft(dd.*window);
df = tmp(1:N+1);
tmp = fft(ee.*window);
ef = tmp(1:N+1);
上述代码是把xx，dd，ee加窗后做FFT变换，并且取了$fs/2$的频谱部分存放到xf，df以及ef中。加窗的目的是为了减小频谱泄露，提高谱分辨率。
xfwm(:,1) = xf;
xf = xfwm(:,dIdx);
%fprintf(1,'%d: %f\n', kk, xf(4));
dfm(:,1) = df;
将xf存放到xfwm的第一列，xfwm是65*16的矩阵，df做类似处理。 然后把dIdx指向的那一列传给xf，dIdx是之前计算的权重矩阵能量最大的那块的索引，这样从xfwm矩阵里把真正要处理近端信号对应的远端信号（speaker，参考信号）获取到。
Se = gamma*Se + (1-gamma)*real(ef.*conj(ef));
Sd = gamma*Sd + (1-gamma)*real(df.*conj(df));
Sx = gamma*Sx + (1 - gamma)*real(xf.*conj(xf));
计算ef，df和xf的平滑功率谱，gamma这里取值是0.92.相对于8k信号取值略大。它们都是65*1的矩阵，包括了直流分量的能力，剩下的64点是$fs/2$及以下的频点能量。
Sxd = gamma*Sxd + (1 - gamma)*xf.*conj(df);
Sed = gamma*Sed + (1-gamma)*ef.*conj(df);
计算互功率谱，这里计算了远端信号和近端信号功率谱，如果该值较大，则说明它们的相关性较强，说明近端信号采集到了远端信号的概率很大，这时需要进行噪声抑制，同样的如果误差信号和近端信号功率谱较大，则说明估计的$\hat{y}$是比较准的，误差信号里剩余的远端信号较少，需要进行噪声抑制的概率就小。它们也都是65*1矩阵，对应频点的bin。但是上述获得的是绝对值而非相对值，门限设置不容易，需要一个归一化的过程。归一化的过程可以通过求互相关得到。
cohed = real(Sed.*conj(Sed))./(Se.*Sd + 1e-10);
cohedAvg(kk) = mean(cohed(echoBandRange));
cohxd = real(Sxd.*conj(Sxd))./(Sx.*Sd + 1e-10);
如上，计算误差信号和近端信号的互相关，1e-10是为了防止除0报错。cohed越大，表示回声越小，cohxd越大，表示回声越大，这里就可以设置一个统一的门限评判上下限了。
cohedMean = mean(cohed(echoBandRange));
计算设置的echoBandRange里频点的均值，如果echoBandRange设置的过大，则低音部分如鼓点声则不易消掉。
hnled = min(1 - cohxd, cohed);
这里的作用就是把最小值赋值给hnled，该值越大，则说明越不需要消回声。之所以取二者判断，是为了最大可能性的消除掉回声。
[hnlSort, hnlSortIdx] = sort(1-cohxd(echoBandRange));
[xSort, xSortIdx] = sort(Sx);
对1-cohxd（echoBandRange）进行升序排序，同样对Sx也进行升序排序。
hnlSortQ = mean(1 - cohxd(echoBandRange));
对远端和近端信号的带内互相关求均值。hnlSortQ表示远端和近端不相关性的平均值，其值越大约没有回声，与cohed含义一致。
[hnlSort2, hnlSortIdx2] = sort(hnled(echoBandRange));
对hnled进行升序排序。
hnlQuant = 0.75;
hnlQuantLow = 0.5;
qIdx = floor(hnlQuant*length(hnlSort2));
qIdxLow = floor(hnlQuantLow*length(hnlSort2));
hnlPrefAvg = hnlSort2(qIdx);
hnlPrefAvgLow = hnlSort2(qIdxLow);
这里主要取了两个值，一个值取在了排序后的3/4处，一个值取在了排序后的1/2处。
if cohedMean
>
0.98
&
hnlSortQ
>
0.9
suppState = 0;
elseif cohedMean
<
0.95 | hnlSortQ
<
0.8
suppState = 1;
end
如果误差和近端信号的互相关均值大于0.98，且远端和近端频带内的互不相关大于0.9，则说明不需要进行回声抑制，将suppState标志设置成0，如果误差和近端信号小于0.95或者远端和近端频带内信号不相关性小于0.8则需要进行印制，在这个范围之外的，回声抑制标志保持前一帧的状态。
if hnlSortQ
<
cohxdLocalMin
&
hnlSortQ
<
0.75
cohxdLocalMin = hnlSortQ;
end
cohxdLocalMin的初始值是1，表示远端和近端完全不相关，这里判断计算得到的远端近端不相关性是否小于前一次的不相关性，如果小于且不相关性小于0.75，则更新cohxdLocalMin。
if cohxdLocalMin == 1
ovrd = 3;
hnled = 1-cohxd;
hnlPrefAvg = hnlSortQ;
hnlPrefAvgLow = hnlSortQ;
end
如果cohxdLocalMin=1，则说明要么是发现远端和近端完全不相关，要么就是cohxdLocalMin一直没有更新，既然不相关性很大，那么也说明有回声的可能性小，那么使用较小的ovrd（over-driven）值，和较大的hnled（65*1）值。
if suppState == 0
hnled = cohed;
hnlPrefAvg = cohedMean;
hnlPrefAvgLow = cohedMean;
end
如果suppState==0，则说明不需要进行回声消除，直接用误差近端相关性修正误差信号，hnl的两个均值读取cohed的均值，在这种情况下hnled的值接近于1.
if hnlPrefAvgLow
<
hnlLocalMin
&
hnlPrefAvgLow
<
0.6
hnlLocalMin = hnlPrefAvgLow;
hnlMin = hnlPrefAvgLow;
hnlNewMin = 1;
hnlMinCtr = 0;
if hnlMinCtr == 0
hnlMinCtr = hnlMinCtr + 1;
else
hnlMinCtr = 0;
hnlMin = hnlLocalMin;
SeLocalMin = SeQ;
SdLocalMin = SdQ;
SeLocalAvg = 0;
minCtr = 0;
ovrd = max(log(0.0001)/log(hnlMin), 2);
divergeFact = hnlLocalMin;
end
end
hnlLocalMin是hnl系数的最小值，在满足这条判断的情况下发现了更小的值，需要对其进行更新，同时表标志设置成1，计数清0，这种情况下回声的概率较大，所以把ovrd设大以增强抑制能力。
if hnlMinCtr == 2
hnlNewMin = 0;
hnlMinCtr = 0;
ovrd = max(log(0.00000001)/(log(hnlMin + 1e-10) + 1e-10), 5);
end
hnlMinCtr==2，说明之前有满足<0.6的块使得hnlMinCtr=2，然后其下一块又没有满足<0.6的条件，进而更新了ovrd值。默认是和3比较取最大值，这里调节成5是为了增加抑制效果，实际情况还需要针对系统调试。
hnlLocalMin = min(hnlLocalMin + 0.0008/mult, 1);
cohxdLocalMin = min(cohxdLocalMin + 0.0004/mult, 1);
跟新上述两个值，mult是$fs/8000$,对于16kHz，就是2.就是0.0004和0.0002的差异。
if ovrd
<
ovrdSm
ovrdSm = 0.99*ovrdSm + 0.01*ovrd;
else
ovrdSm = 0.9*ovrdSm + 0.1*ovrd;
end
平滑更新ovrdSm，上述结果倾向于保持ovrdSm是一个较大的值，这个较大的值是为了尽量抑制回声。
ekEn = sum(Se);
dkEn = sum(Sd);
按行求和，物理意义分别是误差能量和近端信号能量。
发散处理
if divergeState == 0
if ekEn
>
dkEn
ef = df;
divergeState = 1;
end
else
if ekEn*1.05
<
dkEn
divergeState = 0;
else
ef = df;
end
end
如果不进行发散处理，误差能量大于近端能力，则用近端频谱更新误差频谱并把发散状态设置成1，如果误差能量的1.05倍小于近端能量，则发散处理标志设置成0.
if ekEn
>
dkEn*19.95
WFb=zeros(N+1,M); % Block-based FD NLMS
end
如果误差能量大于近端能量×19.95则，将权重系数矩阵设置成0.
ekEnV(kk) = ekEn;
dkEnV(kk) = dkEn;
相应能量存放在相应的向量中。
hnlLocalMinV(kk) = hnlLocalMin;
cohxdLocalMinV(kk) = cohxdLocalMin;
hnlMinV(kk) = hnlMin;
上述三个向量保存对应值。
平滑滤波器系数和抑制指数
wCurve = [0; aggrFact*sqrt(linspace(0,1,N))' + 0.1];
权重系数曲线生成，线性均匀分布。
hnled = weight.*min(hnlPrefAvg, hnled) + (1 - weight).*hnled;
使用权重平滑hnled，wCurve分布是让65点中频率越高的点本次hnled占的比重越大，反之则以前的平滑结果所占比重大。
od = ovrdSm*(sqrt(linspace(0,1,N+1))' + 1);
产生65*1的曲线，用作更新hnled的幂指数。
hnled = hnled.^(od.*sshift);
od作为幂指数跟新hnled。
输出回声消除后的信号
hnl = hnled;
ef = ef.*(hnl);
用hnl系数与误差频谱相乘，即频域卷积，就是将误差信号通过了传递函数为hnnl的滤波器。
ovrdV(kk) = ovrdSm;
hnledAvg(kk) = 1-mean(1-cohed(echoBandRange));
hnlxdAvg(kk) = 1-mean(cohxd(echoBandRange));
hnlSortQV(kk) = hnlPrefAvgLow;
hnlPrefAvgV(kk) = hnlPrefAvg;
相关值的暂存 没有开启舒适噪声产生，则Fmix=ef。
% Overlap and add in time domain for smoothness
tmp = [Fmix ; flipud(conj(Fmix(2:N)))];
mixw = wins.*real(ifft(tmp));
mola = mbuf(end-N+1:end) + mixw(1:N);
mbuf = mixw;
ercn(pos:pos+N-1) = mola;
则使用重叠想加法获得时域平滑信号。
XFm(:,2:end) = XFm(:,1:end-1);
YFm(:,2:end) = YFm(:,1:end-1);
xfwm(:,2:end) = xfwm(:,1:end-1);
dfm(:,2:end) = dfm(:,1:end-1);
全体后移一个块，进入下一块迭代。 执行完了之后，如果想听回声消除后信号的声音，使用如下命令：sound(10*ercn,16000)其中16000是输入信号的频率。WebRTC MATLAB and c code github Address
c调用实例见:
https://github.com/shichaog/WebRTC-audio-processing/blob/master/src/webrtc_audio_processing.cc
speex AEC
代码见上述库里的speex_aec.m文件. 首先需要自己读取文件并设置相关的初始值 给出自己的一个样例
fid=fopen('near.pcm', 'rb'); % Load far end
ssin=fread(fid,inf,'float32');
fid=fopen('far.pcm', 'rb'); % Load fnear end
rrin=fread(fid,inf,'float32');
ssin=ssin(1:4096*200);
rrin=rrin(1:4096*200);
Fs=16000;
filter_length=4096;
frame_size=128;
speex_mdf_out = speex_mdf(Fs, rrin, ssin, filter_length, frame_size);
执行完之后，需要播放出来听：
sound(speex_mdf_out.e,16000)