回声消除（AEC）初体验之matlab仿真

从代码传送门查看本博客所使用测试的源码。

一个普通的音频系统结构如下图所示。

假设结构中所有的增益都是1，则可以得到以下关系式。

由上式可以看出来spk发出的信号，会被麦克风采集到，然后又通过codec芯片重新播放出来，从而形成回声。为抑制这种现象提出回声消除算法。

经过回声消除后的理想效果如以下公式，即从MIC传回codec芯片的信号仅仅包含信号源A的声音信号。

回声消除算法功能模块主要有两个：一个是检测声源信号的检测器；另一个是用来生成参考信号的滤波器。

自适应滤波算法用于自动跟踪输入的参考信号，根据MICin与REFin信号之间的相关性，生成随着输入信号变化的滤波器参数，从而生成预测的回声信号。现在主流的自适应滤波算法是最小均方（LMS-least mean square）算法，以及基于LMS算法衍生出来的一系列归一化最小均方（NLMS）算法、块最小均方误差（BLMS）算法、归一化块处理最小均方（NBLMS）算法、Filter-X LMS算法，另外还有递推最小二乘（RLS）算法、仿射投影滤波器（affine projection filter）算法等等。

app中算法相关的最佳参数搭配通过Filter.plot文件计算，如我想计算RLS滤波算法的参数length、ff可以执行如下代码。

clear all; 
[x,Fs] = audioread('handel.wav');  
[y,Fs1]=audioread('handel_echo.wav');


long=[2 4 8 16 32 64 128];
ff=(0.9:0.002:1);
co_lms=zeros(length(long),length(mu));
co_rls=zeros(length(long),length(ff));
t0=clock;
StartMatlabPool_fun(2);
for i=1:length(long)
    for j=1:length(ff)
       ha =dsp.RLSFilter('Length',long(i),'ForgettingFactor',ff(j));
        [res,e] = ha(x,y);
        co_rls(i,j)=corr(x,res);
    end
end
etime(clock,t0)
CloseMatlabPool_fun;
value=max(max(co_rls));
[row, col]=find(value==co_rls)
[L,F]=meshgrid(long,ff);
figure;
surf(L,F,co_rls');
title('Correlation Function of Filter Order and Stepsize');
xlabel('Filter Order'); ylabel('Stepsize'); zlabel('Correlation');

执行结果如下。

>> test_plot
Starting parallel pool (parpool) using the 'local' profile ...
connected to 2 workers.

ans =

  444.8310

Parallel pool using the 'local' profile is shutting down.

row =

     1


col =

    51

根据结果分析long=2(当ff=1时候，相关系数取值已经趋近于1，long取值对其影响已经可以忽略），ff=1

在echo_algo_hallelu.m文件中修改代码如下，然后执行。

clear all;  

[x,Fs] = audioread('handel.wav');  
[y,Fs1]=audioread('handel_echo.wav');

ha =dsp.RLSFilter(2, 'ForgettingFactor', 1);
[res,e] = ha(x,y);

%Regenerated Sound from Echoed Sound
sound(res,Fs);
audiowrite('handle Regenerated.wav',res,Fs)

%Original Signal
subplot(4,1,1); plot(x);
title('Original Signal');
xlabel('Time Index'); ylabel('Signal Value');

%Echoed Signal
subplot(4,1,2); plot(y);
title('Echoed Signal');
xlabel('Time Index'); ylabel('Signal Value');

%Regenerated Signal
subplot(4,1,3); plot(e);
title('Regenerated Signal');
xlabel('Time Index'); ylabel('Signal Value'); 
co=corr(x,res);
subplot(4,1,4); bar(co);

得到如下执行结果

或者在对应app中设置length=2，Fogettingfactor=1得到以下结果。

对于LMS滤波器参数寻优，则将代码改动如下。

clear all;  
[x,Fs] = audioread('handel.wav');  
[y,Fs1]=audioread('handel_echo.wav');


long=[2 4 8 16 32 64 128];
mu=(0.002:0.02:0.8);
ff=(0.9:0.002:1);
co_lms=zeros(length(long),length(mu));
co_rls=zeros(length(long),length(ff));
bar=waitbar(0,'参数寻优进度');
steps=length(long)*length(mu);
StartMatlabPool_fun(4);
for i=1:length(long)
    for j=1:length(mu)
       ha = dsp.LMSFilter('Length',long(i),'StepSize',mu(j));
        [res,e] = ha(x,y);
        co_lms(i,j)=corr(x,res);
    end
    waitbar((i*j)/steps);
end
CloseMatlabPool_fun;
close(bar);
value=max(max(co_lms));
[row, col]=find(value==co_lms);
length=long(row)
stepSize=mu(col)
[L,M]=meshgrid(long,mu);
figure;
surf(L,M,co_lms');
title('Correlation Function of Filter Order and Stepsize');
xlabel('Filter Order'); ylabel('Stepsize'); zlabel('Correlation');

执行结果如下。

test_plot
Starting parallel pool (parpool) using the 'local' profile ...
connected to 4 workers.
Parallel pool using the 'local' profile is shutting down.

length =

     2


stepSize =

    0.0220

同样在app中验证此组参数的效果。

在信噪比大于 60 dB 的理想信道中仿真自适应算法，LMS 和 NLMS 具有很好的收敛特性，回声可以有大比例的衰减。实际通信信道中，绝大多数近端场所（如会场）都有较大环境背景声音，回声比背景声电平只大 0~20 dB，在这样的系统中，LMS 和 NLMS 算法几乎无法收敛。AEC 算法要求 y(n)序列、x(n)序列保持严格同步（延时控制在100mS之内应该就可以，而且应该也可以提前获取固定延时量，然后在算法中做补偿），即有固定的小于预测窗的延迟时间，否则预测无法收敛。

回声检测只需要检测SPK是否有声音输出，如果输出就控制自适应算法模块生成回声预测信号，把MICin中所有关于SPK的部分全部滤掉，不用分近端远端。

滤波器一般用FIR或者IIR都可以。

集中列出app中不同的滤波器最优参数对应的回声消除效果对比。

在计算过程中，matlab通过matlabpool_open指令打开了并行计算模式。

并行计算模式启动打开指令如下（我电脑上最大允许开启四个核，默认的也是四个核）

StartMatlabPool_fun(4);
%%  此间放置需要进行并行计算的代码部分

CloseMatlabPool_fun;

实现函数代码如下

function StartMatlabPool_fun(size)
%====================================================================
%启动 MATLAB 并行运算（ling）
%StartMatlabPool(size)
%====================================================================
%输入参数：
%   [1] size 为 CPU 核心数量
p = gcp('nocreate'); 
if isempty(p)
    poolsize = 0;
else
    poolsize = p.NumWorkers;
end
if poolsize == 0
    if nargin == 0
        parpool('local');
    else
        try
            parpool('local',size);
        catch ce
            parpool;
            fail_p = gcp('nocreate');
            fail_size = fail_p.NumWorkers;
            display(ce.message);
            display(strcat('输入的size不正确，采用的默认配置size=',num2str(fail_size)));
        end
    end
else
    disp('parpool start');
    if poolsize ~= size
        closematlabpool();
        startmatlabpool(size);
    end
end

我对比了使用四核并行计算和两核并行计算之间的速度，以下是加速的效果，可以看出节省了12秒钟≈3%的时间节约，但考虑到Matlab重复运行相同的程序速度回一次比一次快的特点，这12秒钟时间是不是并行计算省出来的真不好说。

>> test_plot
parpool start

ans =

  426.6580

Parallel pool using the 'local' profile is shutting down.
>> test_plot
Starting parallel pool (parpool) using the 'local' profile ...
connected to 2 workers.

ans =

  438.2910

Parallel pool using the 'local' profile is shutting down.