一、speex回聲消除aec小析：

頻域自適應演算法採用了分塊處理的思想，以減少高階自適應濾波器的計算複雜度，多延遲自適應濾波器（MDF)則更一般可以分M塊來處理，其中塊的長度比自適應濾波器的階數更小。

        後置濾波器係數一直在更新，單講時前置濾波器用後置濾波器係數更新，雙講時前置濾波器係數不更新。

        如果近端遠端的資料線性比較好的話，用speex處理回聲效果還不錯。若是非線性嚴重效果就差了。另外memath認為相比其它的回聲消除演算法，speex對在濾波器更新時做了遠端檢測，回聲估計比較好些。

以前寫的文章，分享下，不當之處，請指正。

二、原始碼解讀：

標記：N=2*frame_size

      M= filter_length/frame_size，其中filter_length為回聲消除拖尾

1、按照公式（1）把近端訊號in[frame_size]直流去除掉，存入快取input[frame_size]。

2、按照公式（2）對遠端與近端訊號分別進行預加重處理，遠端訊號存入x[frame_size]，並對近端延遲的訊號進行frame_size的交疊，以下均按預加重處理後訊號處理。

3、頻域X[M*frame_size]後面M-1個頻域段的歷史資料後移，然後對遠端訊號x使用FFT變換到頻域X[M*frame_size]前面部分。

4、求遠端訊號x時域的能量Sxx和以及其功率譜Xf[frame_size+1]。   

5、前置濾波器foreground與遠端訊號x的輸出再與近端訊號相減得到誤差訊號存入e[N]前半部分，並求其能量方差Sff。

6、求取後置濾波器W各段的係數幅值並進行歸一化處理

，

其中初始化max_sum=1。

更新系數w(i)=w(i)+PHI(i)，PHI由遠端訊號頻域X(i)與誤差E(i)複數相乘並乘以power_1*prop(i)計算得到，其中power_1為自適應步長因子。

7、為防止迴圈卷積，對於第1段以及當前幀數對M-1求餘為j-1的j(j=0...M-1)段作處理。（把W[j*N]段做反傅立葉變換後的後一半部分全部置0，再做傅立葉變換回來）

8、計算遠端訊號與後置濾波器的輸出並反傅立葉變換到時域y，前置、後置濾波器輸出之間的誤差與它們誤差的方差Dbf；計算近端訊號減去後置濾波器輸出差值的方差See。

9、分別計算平均能量差別與平均方差差別Davg1= 0.6*Davg1 + 0.4*(Sff-See); Davg2 = 0.85*Davg2 + 0.15*(Sff-See);Dvar1 =0.36*Dvar1 + 0.16*Sff*Dbf；Dvar2 = 0.7225*Dvar2 + 0.0225*Sff*Dbf；以便判斷是否要更新前置後置濾波器係數。

如果;其中在公式(3)中d1=1，d2=0.5，d3=0.25任一個條件成立，就要更新前置濾波器係數。若是要更新，則把後置濾波器的係數複製給前置濾波器，並平滑前置濾波器的輸出誤差與後置濾波器的輸出誤差得到總的誤差。Davg1，Davg2，Dvar1，Dvar2重置為0。

10、若在公式(4)中d1=4，d2=4，d3=4任一個條件成立，就要更新後置濾波器係數。則把前置濾波器的係數複製給後置濾波器，並把前置濾波器的輸出誤差複製給後置濾波器的輸出誤差，再更新後置濾波器輸出前部分。把See=Sff,Davg1，Davg2，Dvar1，Dvar2重置為0。

11、自適應濾波器的輸出out[frame_size]由近端訊號減去前置濾波器的輸出，並且去預加重處理公式（4）；儲存前置濾波器生成人造回聲訊號並去預加重處理公式（5）。

12、誤差快取區前部分拷貝給後部分，並置前部分為0。

13、計算前置濾波器輸出與前置濾波器的輸出的互相關Sey，後置濾波器輸出的自相關Syy，近端訊號的自相關Sdd；計算分別前置、後置濾波器輸出誤差的頻域資料E與Y，及其功率譜Rf與Yf。

14、若!(Syy>=0 且Sxx>=0且 See >= 0)|| !(Sff < N*1e9且 Syy < N*1e9且 Sxx< N*1e9),則screwed_up=screwed_up+50;否則Sff>Sdd+10000N時，screwed_up=screwed_up+1;否則screwed_up=0；當screwed_up>=50時，重置全部引數。

15、計算遠端訊號的能量Sxx及其功率譜Xf;平滑遠端訊號的功率

16、計算相關係數及標準差

Eht=Rf(i)-Eh(i);

Yht=Yf(i)-Yh(i);

Peyt=Peyt+Eht*Yht;

Pyyt=Pyyt+Yht*Yht;

更新：Eh(i)=(1- N/Fs)*Eh(i)+N/Fs*Rf(i)

Yh(i)=(1- N/Fs)*Yh(i)+N/Fs*Yf(i)

最後Pyyt=sqrt(Pyyt),  Peyt=Peyt/Pyyt;

17、計算相關係數更新率alpha=min(Syy*2*frame_size/Fs,See*0.5*frame_size/Fs)/See; alpha_1=1-alpha。

18、更新相關係數Pey=alpha_1*Pey+ alpha*Peyt; Pyy=alpha_1*Pyy+ alpha*Pyyt;

19、Pyy=max(Pyy,1);Pey=max(Pey,0.005*Pyy);Pey=min(Pey,Pyy)

20、洩漏估計leak_estimate=Pey/Pyy,若leak_estimate>16383,則令leak_estimate=32767。

21、計算錯誤比率RER = (0.0001*Sxx + 3*leak_estimate*Syy)/ See;對RER作上下限限制。

RER=max(RER, Sey*Sey/(1+See*Syy);RER=min(RER,0.5);

22、若sum_adapt>M並且leak_estimate>0.03,則更新自適應濾波器步長因子:先求後置濾波器的回聲消除功率         r=min(leak_estimate*Yf(i),0.5(Rf(i)+1))*0.7+RER*(Rf(i)+1)*0.3，

最優步長為power_1=r/((Rf(i)+1)*(power(i)+10))；否則，adapt_rate=0,若Sxx>1000N時，

adapt_rate=min(Sxx, See)*0.25/See;更新power_1= adapt_rate/(power+10)，

總自適應次數sum_adapt=sum_adapt+adapt_rate。

三、speex回聲消除流程圖：