Speex AEC代码分析1——MDF公式1和2_for (i=0;iframe_size;i++) st->last_y[st->fram

 

微信扫一扫,关注公众号“音频算法与工程实践”

Speex AEC算法分两部分：自适应滤波和最优步长因子的计算，分别基于以下两篇论文：

[1] J.-S. Soo and K. Pang, Multidelayblock frequency domain adaptive filter, 1990

[2] Jean-Marc Valin, On Adjusting theLearning Rate in Frequency Domain Echo Cancellation With Double-Talk, 2007

网上已经有博客讲解Speex AEC的代码流程，这种做法容易陷入代码运算细节而忽略对算法原理的整体把握。本系列文章先根据论文公式讲解代码实现，然后梳理AEC算法的前因后果，最后再回顾这两篇论文。

官网speex-1.1.11.1.tar.gz是浮点AEC，代码与论文基本一致。speex-1.1.12.tar.gz的AEC代码支持浮点和定点版本，代码新增加了一些其他处理。因此本文选择speex-1.1.11.1.tar.gz讲解，MDF论文公式如下图。

公式中N是指滤波器系数的长度。

公式中FFT长度N'，论文要求N'>=2N/M，2是overlap-save方法要求使用两帧数据。

公式中delay blocks数量M=取整(滤波器长度/帧长)。

因此AEC初始化函数speex_echo_state_init()计算代码如下：frame_size、window_size、M分别对应公式中的N、N'、M。

1. /** Creates a new echo canceller state */

2. SpeexEchoState*speex_echo_state_init

3. (int frame_size, int filter_length)

4. {

5. SpeexEchoState *st = (SpeexEchoState *)

6. speex_alloc(sizeof(SpeexEchoState));

7. st->frame_size= frame_size;//160

8. st->window_size= 2*frame_size;//320

9. st->M =

10. (filter_length+st->frame_size-1)/frame_size;

11. //(8*160+160-1)/160=8

12. return st;

13. }

测试代码testecho.c如下，设置参数frame_size=160、filter_length=8*160、window_size=320、M=8，对应公式中的N=8*160、N' =320、M=8。

1. #define NN 160

2. st = speex_echo_state_init(NN,8*NN);

公式1是将连续的2帧音频数据从时域转换到频域，公式2是更新前M-1帧频域缓存数据，因此公式2应该在公式1的FFT计算之前完成。

对应AEC处理函数speex_echo_cancel()代码如下。

1. /** Performs echo cancellation on a frame */

2. void speex_echo_cancel(SpeexEchoState *st,

3. short*ref, short *echo,

4. short *out, float *Yout)

5. {

6. int i,j;

7. int N,M;

8. N =st->window_size;//320

9. M =st->M;//8

10. /* Copyinput data to buffer */

11. for(i=0;i<st->frame_size;i++)

12. {

13. st->x[i] = st->x[i+st->frame_size];

14. st->x[i+st->frame_size] = echo[i];

15. st->d[i] = st->d[i+st->frame_size];

16. st->d[i+st->frame_size] = ref[i];

17. }

18. /* Shiftmemory:

19. this could be optimized eventually*/

20. for(i=0;i<N*(M-1);i++)

21. st->X[i]=st->X[i+N];

22. /* Copynew echo frame */

23. for(i=0;i<N;i++)

24. st->X[(M-1)*N+i]=st->x[i];

25. /* Convertx (echo input) to

26. frequency domain */

27. spx_drft_forward(st->fft_lookup,

28. &st->X[(M-1)*N]);

29. }

代码说明：

1、在现在的技术交流中，ref一般指远端信号，echo指麦克风接收到的回声信号。但是V1.1.11.1代码含义刚好相反：变量echo是送给喇叭的声音(远端信号)，也是自适应滤波器的输入信号。变量ref是麦克风采集的声音，也是自适应滤波器的期望输出信号。V1.1.12已经修改了这两个变量命名。

2、代码中变量N=window_size=320对应公式中N'，与公式中N的含义不同。

3、spx_drft_forward输入输出buf都是同一个buf，频域数据丢弃两个0元素i[0]和i[N/2]，因此时域和频域所需的buf长度也相同。最终的顺序是r[0], r[1],i[1],r[2],i[2],  … , r[N/2-1],i[N/2-1], r[N/2]。

4、下列代码实现公式2，

1. /* Shiftmemory:

2. this could be optimized eventually*/

3. for(i=0;i<N*(M-1);i++)

4. st->X[i]=st->X[i+N];

变量d缓存的数据在下一帧不会再用，因此这2行代码多余。计算公式4时直接使用变量ref。

1. st->d[i] = st->d[i+st->frame_size];

2. st->d[i+st->frame_size] = ref[i];

其余代码都是实现公式1。