基于FDAF的ACE实现

题目说明

• ACE 声学回声消除
• FDAF 频域自适应滤波算法

回音产生原因

声学回音

• 麦克风直接采集到喇叭播放的声音：直接耦合
• 麦克风采集到空间反射造成的回声：远端信号

• 直接耦合：距离很近，基本为接收信号延迟，一般设备会直接滤除。
• 远端信号：成分复杂，为接收信号的多次残留叠加。

线路回音

• 由于2-4线匹配耦合转换引入的线路回音，不做详细介绍。

系统建模

先验条件

1. 麦克风采集信号有延迟
2. 开始对话前有空白时间
3. 对话过程会产生双向通话，产生混叠。

设：
接收信号的人为$A$，发送信号到$B$
接收信号为$x$
经过空间反射的回声为$y$
近端麦克风采集信号为$v$，
那么麦克风采集声音为$d=v+y$
$x$和$y$有很高的相关度
希望找到一个关系来模拟回声路径,设$y‘ = hx$
理想输出为$v=d-y'$
误差$e=y-y'=y-hx$
对于一个固定场景，$h$最终应该是一个较为稳定的场景。

假设$B$先说话，对于$A$短时间处有：
$v=0$，麦克风采集声音为$y$，理想输出为$0$，实际为$e$此时可以对$h$估计；
双向通话中无法确定理想输出$v$，此时无法对h进行估计；
反射回声$y=0$时无法更新。
故$h$的估计建立在对话开始很短一段时间和单项通话过程。

完整的回声消除系统

此处主要针对后两个进行建模处理。

1. 时延估计
2. （线性）回声消除
3. 双讲检测
4. 非线性残余声学回声抑制

算法实现

FDAF算法

设$x$分块$L$，$h$长度$M$，一般取$L=M$

$$\begin{split} &X_k=FFT\begin{bmatrix} x_{k-L}\\ x_k \end{bmatrix}\\ &2M点后端补0，H_k = FFT\begin{bmatrix} h_{k}\\ 0 \end{bmatrix}\\ &Y_k = X_k \odot H_k，取后M点\\ &y_k = IFFT(Y_k)\\ &e_k = d_k - y_k\\ &2M点前端补0，E_k = FFT\begin{bmatrix} 0\\ e_k \end{bmatrix}\\ &\nabla _k = IFFT(X_K^H\odot E_k)，取前M点\\ &H_{k+1} = H_k + \mu FFT\begin{bmatrix} \nabla _{k-L}\\ 0 \end{bmatrix} \end{split}$$

在数值上麦克风期望输出等于$e_k$

双讲检测&确定h更新条件

采用DTD信号相关度，余弦相似性。
当没有近处信号$v$,
麦克风采集信号$d=v+y=y$，和输入信号相似度很大:

$$\begin{split} \frac{x^Td}{\sqrt{x^Tx+d^Td}}=&\frac{x^Ty}{\sqrt{x^Tx+y^Ty}}\\ =&\frac{x^Thx}{\sqrt{x^Tx+(hx)^Thx}}\\ \approx & 1 \end{split}$$

当双向通信时候麦克风采集信号$d$和输入信号$x$相似度很小:

$$\frac{x^Td}{\sqrt{x^Tx+d^Td}}\approx 0$$

设置合适阈值判断是否更新

1. 当$v$，$y$均为$0$，$h$不更新;
2. 当$v=0$，$y \ne$，$h$开始更新;
3. 当$v \ne 0$，$y=0$，$h$不更新.

为了确保1.可以有效计算，修正公式，其中$\epsilon$为一很小的给定常数。

$$\frac{x^Td}{\sqrt{x^Tx+d^Td+\epsilon}}$$

性能分析&对比NLMS算法：

此处选择了一段5s以后开始的双讲混合音频。音频采样率16k，滤波器阶数2k。
从时域上来看还是很难看出性能，但是看时频谱图就可以看出FDAF算法还是很强的。
其实就是把NLMS的batch_size开大了。

音频对比

原始音频

NLMS滤波

DADF滤波

NLMS 时域

FDAF 时域

NLMS 对数

FDAF 对数

NLMS 时谱图

FDAF 时谱图

其他应用

• 经过一定程度修改可用于音频中的背景音乐和人声分离。

不足

• 未考虑滤波器参数更新计算时间。
• 环境建模过于理想，实际的回声传播很复杂

代码参考

实际上并没有做双讲检测，但是效果依然不错。

function [en, H] = FDAF(d,x,mu,mu_unconst, M, DTDbegin)
% 参数:
% d                输入信号(麦克风语音)
% x                远端语音
% mu                约束FDAF的步长
% M                 滤波器阶数

x_new = zeros(M,1);     % 将新块的M个样本初始化为0
x_old = zeros(M,1);     % 将旧块的M个样本初始化为0

AdaptStart = 2*M;       % 在获得2M个样本块后开始自适应
H = zeros(2*M,1);       % 将2M个滤波器权重初始化为0
d_block = zeros(M,1);   % 将期望信号的M个样本初始化为0

power_alpha = 0.5;        % 常数以更新每个frequency bin的功率
power = zeros(2*M,1);   % 将每个bin的平均功率初始化为0
d_length = length(d);             % 输入序列的长度
en = [];                       % 误差向量
window_save_first_M = [ones(1,M), zeros(1,M)]';  % 设置向量以提取前M个元素 (2M,1)

for k = 1:d_length
    x_new = [x_new(2:end); x(k)];         % 开始的输入信号块 (2M,1)
    d_block = [d_block(2:end); d(k)];     % 开始的期望信号快 (M,1)
    if mod(k,M)==0                        % If iteration == block length, 
        x_blocks = [x_old; x_new];        % 2M样本的输入信号样本块 (2M,1)
        x_old = x_new;
        if k >= AdaptStart                % 频域自适应滤波器

            % 将参考信号转换到频域
            Xk = fft(x_blocks);     % (2M,1)

            % 计算滤波器估计信号
            Yk = Xk.*H;                  % 输入和权重向量的乘积(2M,1)*(2M,1)=(2M,1)
            temp = real(ifft(Yk));            % IFFT 输出的实部 (2M,1)
            yk = temp(M+1:2*M);               % 抛弃前M个元素，保留后M个元素 (M,1)
            % 计算误差信号
            error = d_block-yk;              % 误差信号块 (M,1)
            Ek = fft([zeros(1,M),error']');   % 在FFT之前插入零块以形成2M块(2M,1)

            % 更新信号功率估算
            power = (power_alpha.*power) + (1 - power_alpha).*(abs(Xk).^2); % (2M,1)
            % 计算频域中的梯度和权重更新
            if k <DTDbegin
                gradient = real(ifft((1./power).*conj(Xk).* Ek));   % (2M,1)
                gradient = gradient.*window_save_first_M;   % 去除后一个数据块，并且补零 (2M,1)
                H = H + mu.*(fft(gradient));    % 权重是频域的 (2M,1)
            else
                gradient = conj(Xk).* Ek;   %  (2M,1)
                H = H + mu_unconst.*gradient;    % (2M,1)
            end
            en = [en; error]; 
        end
    end
end

clear all;

[x, ~] = audioread('ref.wav');
[d, Fs] = audioread('mic_sig_spe_5.wav');
M =2001;
mu = 1;
mu_unconst = 0;
N=length(x); 
DTDbegin = 50000;

[en, W] = FDAF(d,x,mu,mu_unconst, M, DTDbegin);
audiowrite('output_FDAF.wav',en,Fs)
% 
figure(1)
subplot(3,1,1)
plot(d)
xlabel('time (samples)'); 
ylabel('d(n)');
grid on
axis([0 N -1 1]);
subplot(3,1,2)
plot(x)
xlabel('time (samples)'); 
ylabel('x(n)');
grid on
axis([0 N -1 1]);
subplot(3,1,3)
plot(en)
xlabel('time (samples)'); 
ylabel('e(n)');
grid on
axis([0 N -1 1]);


figure(2);
ee= 20 * log10(abs(en));
figure(2)
plot(ee)
xlabel('time (samples)'); 
ylabel('e(n)/dB');
grid on

s = spectrogram(d, 256);
t = linspace(0, 4*16000/Fs, size(s,1));
f = linspace(0, Fs/2, size(s,2));
figure(3);
subplot(1,3,1)
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('dn/Freqency');
colorbar;

s = spectrogram(d, 256);
t = linspace(0, 4*16000/Fs, size(s,1));
f = linspace(0, Fs/2, size(s,2));
subplot(1,3,2)
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('xn/Freqency');
colorbar;

s = spectrogram(en, 256);
t = linspace(0, 4*16000/Fs, size(s,1));
f = linspace(0, Fs/2, size(s,2));
subplot(1,3,3)
imagesc(t, f, 20*log10((abs(s))));xlabel('Samples'); ylabel('en/Freqency');
colorbar;