---

自适应回声消除原理

声学回声是指扬声器播出的声音在接受者听到的同时，也通过多种路径被麦克风拾取到。多路径反射的结果产生了不同延时的回声，包括直接回声和间接回声。

直接回声是指由扬声器播出的声音未经任何反射直接进入麦克风。这种回声的延时最短 ，它同远端说话者的语音能量，扬声器与麦克风之间的距离、角度 ,扬声器的播放音量，麦克风的拾取灵敏度等因素直接相关；

间接回声是指由扬声器播出的声音经过不同的路径 (如房屋或房屋内的任何物体 )的一次或多次反射后进入麦克风所产生的回声的集合。房屋内的任何物体的任何变动都会改变回声的通道。因此，这种回声的特点是多路径的、时变的。

自适应回声消除的基本思想是估计回音路径的特征参数，产生一个模拟的回音路径，得出模拟回音信号，从接收信号中减去该信号，实现回音抵消。其关键就是得到回声路径的冲击响应$\hat{h}(n)$，由于回音路径通常是未知的和时变的，所以一般采用自适应滤波器来模拟回音路径。自适应回音消除的显著特点是实时跟踪，实时性强。

图中$y(n)$代表来自远端的信号 , $r(n)$是经过回声通道而产生的回声，$x(n)$是近端的语音信号。D端是近端麦克风，麦克风采集到的房间叠加的回声和近端说话人的语音。对回声消除器来说，接收到的远端信号作为一个参考信号，回声消除器根据参考信号由自适应滤波器产生回声的估计值$\hat{r}(n)$，将$\hat{r}(n)$从近端带有回声的语音信号减去，就得到近端传送出去的信号 。在理想情况下，经过回声消除器处理后，残留的回声误差$e(n)=r(n)-\hat{r}(n)$将为0，从而实现回音消除。

性能指标

• 收敛速度：滤波器的收敛速度越快越好，使正常通话开始后，通话者很快就感觉不到明显的回波存在。
• 稳态残留回波(稳定性)：即当滤波器收敛达到稳态后的回波输出量，实际中总是希望该参数越小越好。
• 算法复杂度：良好的算法应该在保持收敛速度的同时尽量降低计算复杂度，同时也能减少功耗

ITU-T G.168对各种回音抵消器产品在包括以上两个主要指标在内的各种指标规定了必须达到的标准。

自适应滤波器结构

图中所示滤波器的输入是$X(n)={x(n),x(n-1),...x(n-N+1)}^T$，$z^{-1}$滤波器z域模型的延迟单元（零状态），滤波器的权重系数是$h(n)={h_1(n),h_2(n),...,h_N(n)}^T$，$d(n)$为期望输出信号，$\hat{d}(n)$为滤波器的实际输出，也称估计值，$\hat{d}(n)=\sum_{i=1}^Nx(n-i+1)h_i(n)$。$e(n)$是误差，$e(n)=d(n)-\hat{d(n)}$，由误差经过一定的自适应滤波算法来调整滤波系数，使得滤波器的实际输出接近期望输出信号。

传统的IIR和FIR滤波器在处理输入信号的过程中滤波器的参数固定，当环境发生变化时，滤波器无法实现原先设定的目标。自适应滤波器能够根据自身的状态和环境变化调整滤波器的权重。

自适应滤波器类型。可以分为两大类：非线性自适应滤波器、线性自适应滤波器。非线性自适应滤波器包括基于神经网络的自适应滤波器及Volterra滤波器。非线性自适应滤波器信号处理能力更强，但计算复杂度较高。所以实践中，线性自适应滤波器使用较多。

主要分为两类FIR滤波器、IIR滤波器。

1. FIR滤波器时非递归系统，即当前输出样本仅是过去和现在输入样本的函数，其系统冲激响应h(n)是一个有限长序列。具有很好的线性相位，无相位失真，稳定性较好。
2. IIR滤波器时递归系统，即当前输出样本是过去输出和过去输入样本的函数，其系统冲激h(n)是一个无限长序列。IIR系统的相频特性是非线性的，稳定性不能保证。好处是实现阶数较低，计算量较少。

由于IIR存在稳定性问题，因此一般采用FIR。

自适应滤波器算法按照不同的优化准则，常见自适应滤波算法有：递推最小二乘算法(RLS)，最小均方误差算法(LMS)，归一化均方误差算法(NLMS)，快速精确最小均方误差算法，子带滤波，频域的自适应滤波等等。

全带自适应稀疏算法

谱减法

本文研究的背景噪声是以汽车噪声和风声为主要对象，此类噪声的特点为加性的、局部平稳的、且与语音信号统计独立。谱减法根据噪声的加性特点，通过噪声能量估计和增益计算得到噪声的功率谱，然后从带噪语音的功率谱中减去估计出的噪声得到较为纯净的语音。

由于谱减法没有考虑人耳听觉对语音频谱分布的幅度较为敏感的特点。因此采用谱减法进行滤波后，会给原信号带来噪声，使语音质量变差，并且会影响到其他处理，如语音编码等。用实际的语音信号检测谱减法的滤波效果，如图下所示，可以看出，经过谱减法后噪声得到了一定的消除，但频率较低的信号处滤波效果并不是很理想。

最小均方算法(LMS)

1959年由Widrow和Hoff提出了最小均方（Least Mean Square，LMS）算法，LMS基于维纳滤波理论，采用瞬时值估计梯度矢量的算法，通过最小化误差信号的能量来更新自适应滤波器权值系数。

设计一个N 阶滤波器，它的参数为$w(n)$，则滤波器输出为

$y(n)=\sum_{i=0}^{N-1} w_{i}(n) x(n-i)=w^{T}(n)x(n)=x^{T}(n) w(n)$

期望输出为$d(n)$，则误差信号可以定义为：

$e(n)=d(n)-y(n)=d(n)-w^{T}(n) x(n)$

我们的目标就是将误差$e(n)$最小化，采用最小均方误差（MMSE ）准则，最小化目标函数：$J(w)$

$J(w)=E\left\{|e(n)|^{2}\right\}=E\left\{\left|d(n)-w^{T}(n) x(n)\right|^{2}\right\}$

计算目标函数$J(w)$对$w$的导数，令导数为0：

$\nabla (J(w)) = \frac{{\partial [\mathop e\nolimits^2 (n)]}}{{\partial w(n)}} =  - 2e(n)X(n)$

则滤波器系数的更新公式可以写为：

$w(n+1)=w(n)-\frac{1}{2}\mu [-\nabla (J(w))]=w(n)+\mu e(n)X(n)$

上式中的$\mu$为步长因子。$\mu$值越大，算法收敛越快，但稳态误差也越大；$\mu$值越小，算法收敛越慢，但稳态误差也越小。为保证算法稳态收敛，应使$\mu$在以下范围取值：$0 < \mu  < \frac{2}{{\sum\limits_{i = 1}^N {x(i)_{}^2} }}$

• 优点：算法简单，易于实现，算法复杂度低(LMS<RLS)，能够抑制旁瓣效应
• 缺点：
  • 收敛速率较慢(LMS<RLS)，因为LMS滤波器系数更新是逐点的（每来一个新的$x(n)$和$d(n)$，滤波器系数就更新一次）,每一次采样点梯度的估计对于真实梯度会存在误差，导致滤波器系数的每次更新不会严格按照真实梯度方向更新，而是有一定的偏差
  • 跟踪性能较差，并且随着滤波器阶数(步长参数)升高，系统的稳定性下降
  • LMS要求不同时刻的输入向量$x(n)$线性无关——LMS 的独立性假设。如果输入信号存在相关性，会导致前一次迭代产生的梯度噪声传播到下一次迭代，造成误差的反复传播，收敛速度变慢，跟踪性能变差。

所以，理论上，LMS 算法对白噪声的效果最好。为了降低输入信号的相关性，出现了一类“解相关LMS”算法，这里就不展开讲述了。

function [e, y, w] = myLMS(d, x, mu, M)
% Inputs:
% d  - 麦克风语音
% x  - 远端语音
% mu - 步长，0.05
% M  - 滤波器阶数，也称为抽头数
%
% Outputs:
% e - 输出误差，近端语音估计
% y  - 输出系数，远端回声估计
% w - 滤波器参数

    d_length = length(d);
    if (d_length <= M)  
        print('error: 信号长度小于滤波器阶数！');
        return; 
    end
    if (d_length ~= length(x))  
        print('error: 输入信号和参考信号长度不同！');
        return; 
    end

    xx = zeros(M,1);
    w1 = zeros(M,1);    % 滤波器权重
    y = zeros(d_length,1);  % 远端回声估计
    e = zeros(d_length,1);  % 近端语音估计

    % for n = 1:Ns
    %     xx = [xx(2:M);x(n)];    % 纵向拼接
    %     xx(2~40)+x(1)-->xx(2~40)+x(2)-->xx(2~40)+x(3)...
    %     y(n) = w1' * xx;        % 远端回声估计(40,1)'*(40,1)=1; (73113,1)
    %     e(n) = d(n) - y(n);     % 近端语音估计
    %     w1 = w1 + mu * e(n) * xx;   % (40,1)
    %     w(:,n) = w1;        % (40, 73113)
    % end
    % 和上面类似
    for n = M:d_length
        xx = x(n:-1:n-M+1);    % 纵向拼接  (40~1)-->(41~2)-->(42~3)....
        y(n) = w1' * xx;        % 远端回声估计 (40,1)'*(40,1)=1; (73113,1)
        e(n) = d(n) - y(n);     % 近端语音估计
        w1 = w1 + mu * e(n) * xx;   % (40,1)
        w(:,n) = w1;        % (40, 73113)
    end
end

MATLAB实现LMS

import numpy as np

def lms(x, d, N = 4, mu = 0.05):
    L = min(len(x),len(d))
    h = np.zeros(N)
    e = np.zeros(L-N)
    for n in range(L-N):
        x_n = x[n:n+N][::-1]
        d_n = d[n] 
        y_n = np.dot(h, x_n.T)
        e_n = d_n - y_n
        h = h + mu * e_n * x_n
        e[n] = e_n
    return e

python实现LMS

分块 LMS算法(Block LMS)

为了进一步提高LMS算法的稳定性，在LMS基础上提出了块处理LMS（Block LMS）。LMS算法是每来一个采样点就调整一次滤波器权值；而BLMS算法是每K采样点才对滤波器的权值更新一次。

输入数据序列经过串一并变换器被分成若干个长度为L的数据块。然后将这些数据块一次一块的送入长度为M的滤波器中，在收集到每一块数据样值后，进行滤波器抽头权值的更新，使得滤波器的自适应过程一块一块地进行，而不是像时域传统自适应滤波算法那样一个值一个值地进行。

具体算法如下：

累积$L$个采样点之后，我们再更新一次滤波器系数

$\begin{array}{l} w(n+1)=w(n)+\mu X(n) e(n) \\ w(n+2)=w(n+1)+\mu X(n+1) e(n+1) \\ \mathbf{. .} \\ w(n+L)=w(n+L-1)+\mu X(n+L-1) e(n+L-1) \end{array}$

上式全部相加，中间的式子可以左右相消，就可以得出Block LMS滤波器权重的更新公式：

$w(n+L)=w(n)+\mu \sum_{i=0}^{L-1}X(n+i) e(n+i)$

其中，$\mu$为步长因子，为了推导方便，我们通常令$L=N$，更新一次block LMS相当于LMS更新了$L$倍。我们用一个新的时间索引$k$来代表Block LMS系数的更新，其中$kL=n$，上式可以写为：

$w(k+1)=w(k)+\mu \sum_{i=0}^{L-1} X(k L+i) e(k L+i)$

• 优点：运算量就比LMS的运算量要小的多
• 缺点：收敛速度却与LMS算法相同

% 参考：https://github.com/Morales97/ASP_P2/blob/4a6a17f6fe/algorithms/block_lms.m
function [e, y, w] = myBlock_LMS(d,x,mu, M,L)
    % 输入：
    % d - 麦克风语音
    % x - 远端语音
    % mu - 步长
    % M - 滤波器阶数
    % L - 块大小
    %
    % Outputs:
    % e - 输出误差，近端语音估计
    % y - 输出系数，远端回声估计
    % w - 滤波器参数

    d_length = length(d); 
    K = floor(d_length / L);    % 块的数量，确保整数
    y = zeros(d_length, 1);
    w = zeros(M,K+1);
    e = zeros(d_length,1);
    x_ = [zeros(M-1,1); x];

    % 根据"块"进行循环
    for k=1:K
        block_sum = 0;
        % 求一个块的和sum
        for i = 1 + L*(k-1):L*k
            X = x_(i:i+M-1);
            y(i) = w(:,k)' * X;  % 滤波器输出

            e(i) = d(i) - y(i);
            block_sum = block_sum + X * e(i);   % 权重更新
        end

        w(:,k+1) = w(:,k) + mu * block_sum;
    end

    % 将指针移动一步，使第n行对应于时间n
    w=w(:,2:K+1);
end

Block_LMS MATLAB代码实现

归一化最小均方算法(NLMS)

归一化最小均方(Normalized Least Mean Squares，NLMS)算法是改进的LMS算法，对于较大的输入，会导致梯度噪声的放大，因此需要用输入向量的平方范数进行归一化。根据原LMS算法中误差信号与远端输入信号的乘积，对远端输入信号的平方(功率)进行归一化处理，将固定步长因子的LMS算法变为根据输入信号时变的变步长NLMS算法，具体算法如下：

$\mu(n)=\frac{1}{x^T(n)x(n)}$

$w(n+1)=w(n)+2 \mu(n) x(n) e(n)$

其中$w(n)$为滤波器的系数，$e(n)$为误差信号

• 优点：改善了LMS算法收敛速度慢的缺点。计算简单、更高的精度
• 缺点：输入相关信号时，收敛速率明显下降

function [e,y, w] = NLMS(d, x, M)
    % 输入:
    % d  - 麦克风语音
    % x  - 远端语音
    % M  - 滤波器阶数
    %
    % 输出:
    % e - 近端语音估计
    % y - 远端回声估计
    % w - 滤波器参数
    d_length = length(d);
    if (d_length <= M)
        print('error: 信号长度小于滤波器阶数！');
        return; 
    end
    if (d_length ~= length(x))  
        print('error: 输入信号和参考信号长度不同！');
        return; 
    end

    xx = zeros(M,1);
    w1 = zeros(M,1);        % 滤波器权重
    y = zeros(d_length,1);   % 近端语音
    e = zeros(d_length,1);  % 误差

    for n = 1:d_length
        %在输入信号x前补上M-1个0，使输出y与输入具有相同长度
        xx = [xx(2:M);x(n)];    % (39,1)+(1,1)=(40,1)
        y(n) = w1' * xx;        % (40,1)'*(40,1)=1
        mu = 1/(xx'*xx);    % 步长
        e(n) = d(n) - y(n);     % 误差
        w1 = w1 + mu * e(n) * xx;
        w(:,n) = w1;
    end
end

MATLAB实现NLMS

import numpy as np

def nlms(x, d, N = 4, mu = 0.05):
    L = min(len(x),len(d))
    h = np.zeros(N)
    e = np.zeros(L-N)
    for n in range(L-N):
        x_n = x[n:n+N][::-1]
        d_n = d[n] 
        y_n = np.dot(h, x_n.T)
        e_n = d_n - y_n
        h = h + mu * e_n * x_n / (np.dot(x_n,x_n)+1e-5)
        e[n] = e_n
    return e

python代码实现NLMS

后续研发除了归一化块处理LMS（BNLMS）：结合以上NLMS和BLMS两者的特点则有归一化块处理LMS（BNLMS）

功率归一化最小均方(PNLMS)

$\mu(n)=\frac{\alpha}{\beta+\mathbf{x}^{T}(n) \mathbf{x}(n)}, \quad 0<\alpha<2, \quad \beta \geqslant 0$

$w(n+1)=w(n)+2 \mu(n) x(n) e(n)$

PNLMS 在NLMS的基础上引入了$\alpha$和$\beta$，稳定性要好于 NLMS 。$\beta$一般设为一个较小的整数，防止输入数据矢量$x(n)$的内积过小使得$\mu(n)$过大而引起稳定性能下降，一般取0.0001。当$\alpha=1$，$\beta=0$时， PNLMS 退化为 NLMS 。

function [e,y, w] = PNLMS(d, x, M, aplha, beta)
    % 输入:
    % d  - 麦克风语音
    % x  - 远端语音
    % a  - 偏置参数
    % M  - 滤波器阶数
    %
    % 输出:
    % e - 近端语音估计
    % y - 远端回声估计
    % w - 滤波器参数
    d_length = length(d);
    if (d_length <= M)
        print('error: 信号长度小于滤波器阶数！');
        return; 
    end
    if (d_length ~= length(x))  
        print('error: 输入信号和参考信号长度不同！');
        return; 
    end

    xx = zeros(M,1);
    w1 = zeros(M,1);        % 滤波器权重
    y = zeros(d_length,1);   % 近端语音
    e = zeros(d_length,1);  % 误差

    for n = 1:d_length
        %在输入信号x前补上M-1个0，使输出y与输入具有相同长度
        xx = [xx(2:M);x(n)];    % (39,1)+(1,1)=(40,1)
        y(n) = w1' * xx;        % (40,1)'*(40,1)=1
        mu = aplha/(beta + xx'*xx);    % 步长
        e(n) = d(n) - y(n);     % 误差
        w1 = w1 + mu * e(n) * xx;
        w(:,n) = w1;
    end
end

MATLAB实现PNLMS

变步长LMS(VSS LMS)

针对$\mu$值，人们研究了许多变步长LMS 算法（Variable Step-Size LMS），一般是在滤波器工作的开始阶段采用较大的μ值，以加快收敛速度，而在后阶段采用较小的μ值，可以减小稳态误差。这类算法的关键是确定在整个过程中$\mu$值如何变化或$\mu$值在何种条件满足下才改变。

$w(n+1)=w(n)+\mu(n) x(n) e(n)$

$\mu(n+1)=\alpha\mu(n)+\gamma e^2(n)$

式中$\alpha$和$\gamma$为参数。

另外VSS LMS还对$\mu$进行了约束，$\mu_{max}$为上下学习率，$\mu_{min}$为下界学习率，$\mu_{max}$一般设置为$\mu_{max}=\frac{1}{var(x)*M}$，var为求方差，$x$为输入信号(麦克风信号)，$M$为滤波器阶数：

$\mu(n)=\left\{\begin{matrix} \mu_{max} &\mu(n)>\mu_{max}\\ \mu_{min} &\mu(n)<\mu_{min} \end{matrix}\right.$

• 优点：收敛速度快，
• 缺点：算法的稳定性和跟踪能力上较易受输入噪声的影响

% 学习率的上下界
% 以上信息可从论文中获取，这些值是在下述论文中定义的
% R. H. Kwong and E. W. Johnston, "A variable step size LMS algorithm," in IEEE Transactions on Signal Processing, vol. 40, no. 7, pp. 1633-1642, July 1992.

function [e, y, w] = myVSS_LMS(d, x, M)
% Inputs:
% d  - 麦克风语音
% x  - 远端语音
% mu - 步长，0.05
% M  - 滤波器阶数，也称为抽头数
%
% Outputs:
% e - 输出误差
% y  - 输出系数
% w - 滤波器参数

    d_length = length(d);
    if (d_length <= M)  
        print('error: 信号长度小于滤波器阶数！');
        return; 
    end
    if (d_length ~= length(x))  
        print('error: 输入信号和参考信号长度不同！');
        return; 
    end
    
    % 定义VSS-LMS算法的初始参数
    
    input_var = var(x);     % 输入的方差
    % 如果mu_max和mu_min之间的差异不够大，LMS和VSS LMS算法的误差曲线都是相同的
    mu_max = 1/(input_var*M); % 上界=1/(filter_length * input variance)
    mu_min = 0.0004;        % 下界=LMS算法的学习率
    mu_VSS(1)=1;    % VSS算法的mu初始值
    alpha  = 0.97;
    gamma = 4.8e-4;
    
    xx = zeros(M,1);
    w1 = zeros(M,1);    % 滤波器权重
    y = zeros(d_length,1);  % 近端语音
    e = zeros(d_length,1);  % 误差

    for n = 1:d_length
        xx = [xx(2:M);x(n)];    % 纵向拼接 或者[x(n); xx(1:end-1)]
        y(n) = w1' * xx;        % (40,1)'*(40,1)=1; (73113,1)
        e(n) = d(n) - y(n);
        w1 = w1 + mu_VSS(n) * e(n) * xx;   % 更新权重系数 (40,1)
        w(:,n) = w1;        % (40, 73113)
        
        mu_VSS(n+1) = alpha * mu_VSS(n) + gamma * e(n) * e(n) ;% 使用VSS算法更新mu值
       % 检查论文中给出的mu的约束条件
        if (mu_VSS(n+1)>mu_max)
            mu_VSS(n+1)=mu_max; % max
        elseif(mu_VSS(n+1)<mu_min)
            mu_VSS(n+1)= mu_min;
        else
            mu_VSS(n+1) = mu_VSS(n+1) ;
        end
    end
    % 和上面类似
%     for n = M:d_length
%         xx = x(n:-1:n-M+1);    % 纵向拼接  (40~1)-->(41~2)-->(42~3)....
%         y(n) = w1' * xx;        % (40,1)'*(40,1)=1; (73113,1)
%         e(n) = d(n) - y(n);
%         w1 = w1 + mu * e(n) * xx;   % (40,1)
%         w(:,n) = w1;        % (40, 73113)
%     end
end

MATLAB代码实现VSS LMS

递归最小二乘法(RLS)

MMSE适合处理平稳序列，因为MMSE是一个均匀加权的最优化问题，也就是说，每一时刻的误差信号对目标函数的贡献权重是相同的，如果对于非平稳的语音信号效果就不太好了。

RLS重新定义了目标函数：

$J_{n}(w)=\sum_{i=0}^{n} \lambda^{n-i}|e(i)|^{2}=\sum_{i=0}^{n} \lambda^{n-i}\left|d(i)-w^{T}(n) x(i)\right|^{2}$

其中，滤波器系数选用的是$w^T(n)$，而不是$w^T(i)$。因为在自适应更新过程中，滤波器总是变得越来越好，这意味着对于任何的$i<n$，$|d(i)-w^T(n)x(i)|$总是比$|d(i)-w^T(i)x(i)|$小 。

$\lambda$称为遗忘因子($0<\lambda\leq 1$)。对离$n$时刻越近的误差加比较大的权重，遗忘越少，而对离$n$时刻越远的误差加比较小的权重，遗忘越多

• $\lambda=1$：无任何遗忘功能，此时 RLS 退化为 LMS 方法 
• $\lambda->0$：只对当前时刻的误差起作用，而过去时刻的误差完全被遗忘 

我们让目标函数$J_{n}(w)$对$w$求导，令梯度等于0，得到$w$的公式为：

$w(n)=R^{-1}(n)r(n)=\frac{r(n)}{R(n)}=\frac{\sum_{i=0}^{n} \lambda^{n-i} x(i) d(i)}{\sum_{i=0}^{n} \lambda^{n-i} x(i) x^{T}(i)}$

根据$R(n)$和$r(n)$的等式，得其时间递推公式：

$R(n)=\lambda R(n-1)+x(n)x^T(n)$

$r(n)=\lambda r(n-1)+x(n)d(n)$

弊端：需要计算$R(n)$每个点的逆矩阵，很不划算。更优的方法？

令$P(n)=R^{-1}(n)$，$P(n)$的时间递推公式为：

$P(n)=R^{-1}(n)=(\lambda R(n-1)+x(n)x^T(n))^{-1}$

根据矩阵求逆引理：

$(A+B C D)^{-1}=A^{-1}-A^{-1} B\left(C^{-1}+D A^{-1} B\right)^{-1} D A^{-1}$

令：$A=R(n-1) \quad B=x(n) \quad C=\frac{1}{\lambda} \quad D=x^{T}(n)$

可以得到：

$\begin{aligned} P(n) &=\frac{1}{\lambda}\left[P(n-1)-\frac{P(n-1) x(n) x^{T}(n) P(n-1)}{\lambda+x^{T}(n) P(n-1) x(n)}\right] \\ &=\frac{1}{\lambda}\left[P(n-1)-k(n) x^{T}(n) P(n-1)\right] \end{aligned}$

其中，增益向量：

$k(n)=\frac{P(n-1) x(n)}{\lambda+x^{T}(n) P(n-1) x(n)}$

$\begin{aligned} \mathbf{P}(n) \mathbf{x}(n) &=\frac{1}{\lambda}\left[\mathbf{P}(n-1) \mathbf{x}(n)-\mathbf{k}(n) \mathbf{x}^{T}(n) \mathbf{P}(n-1) \mathbf{x}(n)\right] \\ &=\frac{1}{\lambda}\left\{\left[\lambda+\mathbf{x}^{T}(n) \mathbf{P}(n-1) \mathbf{x}(n)\right] \mathbf{k}(n)-\mathbf{k}(n) \mathbf{x}^{T}(n) \mathbf{P}(n-1) \mathbf{x}(n)\right\} \\ &=\mathbf{k}(n) \end{aligned}$

探索$w(n)$的时间递推公式

$\begin{aligned} \mathbf{w}(n)=& \mathbf{R}^{-1}(n) \mathbf{r}(n)=\mathbf{P}(n) \mathbf{r}(n) \\ =& \frac{1}{\lambda}\left[\mathbf{P}(n-1)-\mathbf{k}(n) \mathbf{x}^{T}(n) \mathbf{P}(n-1)\right][\lambda \mathbf{r}(n-1)+\mathbf{x}(n) d(n)] \\ =& \mathbf{P}(n-1) \mathbf{r}(n-1)-\mathbf{k}(n) \mathbf{x}^{T}(n) \mathbf{P}(n-1) \mathbf{r}(n-1) \\ &+\frac{1}{\lambda} d(n)\left[\mathbf{P}(n-1) \mathbf{x}(n)-\mathbf{k}(n) \mathbf{x}^{T}(n) \mathbf{P}(n-1) \mathbf{x}(n)\right] \\ =& \mathbf{w}(n-1)+d(n) \mathbf{k}(n)-\mathbf{k}(n) \mathbf{x}^{T}(n) \mathbf{w}(n-1) \\ =& \mathbf{w}(n-1)+\mathbf{k}(n)\left[d(n)-\mathbf{x}^{T}(n) \mathbf{w}(n-1)\right] \\ =& \mathbf{w}(n-1)+\mathbf{k}(n) \varepsilon(n) \end{aligned}$

这就是标准RLS算法的更新公式，这里的$\varepsilon (n)$为先验估计误差。

标准RLS 算法的执行流程：

1. 初始化： $w(0)=0，p(0)=\delta ^{-1}I$，$\delta$是一个很小的正数，$I$是单位矩阵。
2. 对每一个时刻$n$，重复如下计算

• • 先验误差：$\varepsilon(n)=d(n)-\mathbf{w}^{T}(n-1) \mathbf{x}(n)$
  • 增益向量：$\mathbf{k}(n)=\frac{\mathbf{P}(n-1) \mathbf{x}(n)}{\lambda+\mathbf{x}^{T}(n) \mathbf{P}(n-1) \mathbf{x}(n)}$
  • 逆矩阵更新：$\mathbf{P}(n)=\frac{1}{\lambda}\left[\mathbf{P}(n-1)-\mathbf{k}(n) \mathbf{x}^{T}(n) \mathbf{P}(n-1)\right]$
  • 滤波器更新：$\mathbf{w}(n)=\mathbf{w}(n-1)+\mathbf{k}(n) \varepsilon(n)$

一个广泛的共识是RLS 算法的收敛速度和跟踪性能都优于 LMS 算法，所付出的代价是需要更复杂的计算 。

算法	误差计算	收敛速度	跟踪性能	计算代价
LMS	后验误差	慢	慢	小
RLS	先验误差	快	快	大

• 优点：RLS自适应滤波器提供更快的收敛速度和跟踪性能。
• 缺点：由于RLS 使用了自相关矩阵的逆矩阵的递推，所以，一旦输入信号的自相关矩阵接近奇异时RLS 的收敛速度和跟踪性能会严重恶化 。

function [e, y, w] = myRLS(d, x, lamda, M)
% Inputs:
% d  - 麦克风语音
% x  - 远端语音
% lamda - the weight parameter, 权重
% M  - the number of taps. 滤波器阶数
%
% Outputs:
% e - 大小为Ns的输出误差向量
% y - 输出的近端语音
% w - 滤波器参数

    Ns = length(d);
    if (Ns <= M)  
        print('error: 信号长度小于滤波器阶数！');
        return; 
    end
    if (Ns ~= length(x))
        print('error: 输入信号和参考信号长度不同！');
        return;
    end

    I = eye(M);
    a = 0.01;
    p = a * I;
    
    xx = zeros(M,1);
    w1 = zeros(M,1);        % 滤波器权重
    y = zeros(Ns, 1);        % 近端语音
    e = zeros(Ns, 1);       % 误差
    
    for n = 1:Ns
        %在输入信号x后补上M-1个0，使输出y与输入具有相同长度
        xx = [x(n); xx(1:M-1)];
        k = (p * xx) ./ (lamda + xx' * p * xx);
        y(n) = xx'*w1;
        e(n) = d(n) - y(n);
        w1 = w1 + k * e(n);
        p = (p - k * xx' * p) ./ lamda;
        w(:,n) = w1;
    end
end

MATLAB代码实现RLS

import numpy as np

def rls(x, d, N = 4, lmbd = 0.999, delta = 0.0002):
    L = min(len(x),len(d))
    lmbd_inv = 1/lmbd
    h = np.zeros((N, 1))
    P = np.eye(N)/delta
    e = np.zeros(L-N)
    for n in range(L-N):
        x_n = np.array(x[n:n+N][::-1]).reshape(N, 1)
        d_n = d[n] 
        y_n = np.dot(x_n.T, h)
        e_n = d_n - y_n
        g = np.dot(P, x_n)
        g = g / (lmbd + np.dot(x_n.T, g))
        h = h + e_n * g
        P = lmbd_inv*(P - np.dot(g, np.dot(x_n.T, P)))
        e[n] = e_n
    return e

python代码实现RLS

仿射投影算法(AP)

仿射投影算法 (Affine projection algorithm) 是运算量介于 LMS 和 RLS 之间的一种自适应算法 。

$y(n)=w^H(n)x(n)=X^{T}(n) w^{*}(n)$

其中，$X(n)$是输入向量矩阵，$w^*(n)$是$w(n)$的共轭转置

权重更新方程：$w(n)=w(n-1)+\mu \Delta w_{n}$

滤波器输出：：$y(n)=X^{T}(n) w^{*}(n-1)+\mu X^{T}(n) \Delta w_{n}^{*}$

令滤波器的实际输出作为期望输出$d(n)$，则先验误差向量表示为：$e(n)=y(n)-X^{T}(n) w^{*}(n-1)$

为了简化，我们令$\mu =1$，则：

$X^{T}(n) \Delta w_{n}^{*}=e(n)$

$X^{H}(n) \Delta w_{n}=e^{*}(n)$

我们发现 滤波器系数的更新量$\Delta W^{*}_n(\Delta W_n)$，是由$L$个线性方程组成的线性方程组决定的

线性方程组相关知识补充：

对于线性方程组：$A_{m*n}X_{n*1}=B_{m*1}$

考虑（行/列）满秩的情况，分下面三种情况：

（一）如果$m=n$，则：

$\begin{aligned} &\mathbf{x}=\mathbf{A}^{-1} \mathbf{b}&\text { 唯一解 精确解 } \end{aligned}$

（二）如果$m<n$即方程的个数小于未知数的个数 此时方程组有无穷多个解 。 为了得到唯一解 必须增加约束条件 要求 x 的范数最小 这样得到的解称为最小范数解 。

$\begin{aligned} &\mathbf{x}=\mathbf{A}^{H} \mathbf{(AA^H)^{-1}b}&\text { 最小范数解 } \end{aligned}$

（三）如果$m>n$即方程的个数大于未知数的个数 此时方程组不存在精确解 只存在近似解 。 我们自然希望找到一个使方程组两边的误差平方和为最小的解 即最小二乘解 。

$\mathbf{x}=\left(\mathbf{A}^{H} \mathbf{A}\right)^{-1} \mathbf{A}^{H} \mathbf{b} \quad \text { 最小二乘解 }$

利用线性方程组的结论，重新观察公式：

$公式19：\mathbf{X}^{H}(n) \Delta \mathbf{w}_{n}=\mathbf{e}^{*}(n)$

分为两种情况讨论：

情况1：$1\leq L< N$，方程组为欠定方程，具有唯一的最小范数解

$\Delta \mathbf{w}_{n}=\mathbf{X}(n)\left(\mathbf{X}^{H}(n) \mathbf{X}(n)\right)^{-1} \mathbf{e}^{*}(n)$

此时，权重更新公式变为：

$\mathbf{w}(n)=\mathbf{w}(n-1)+\mu \mathbf{X}(n)\left(\mathbf{X}^{H}(n) \mathbf{X}(n)\right)^{-1} \mathbf{e}^{*}(n)$

如果$L=1$，则退化为 NLMS 算法：

$\mathbf{w}(n)=\mathbf{w}(n-1)+\mu e^{*}(n) \frac{\mathbf{x}(n)}{\mathbf{x}^{H}(n) \mathbf{x}(n)}$

实际中，L 取 2 、 3 就足够了。

情况2：$L\geq N$，方程组为超定方程，其唯一解为最小二乘解：

$\Delta \mathbf{w}_{n}=\left(\mathbf{X}(n) \mathbf{X}^{H}(n)\right)^{-1} \mathbf{X}(n) \mathbf{e}^{*}(n)$

此时，权重更新公式变为：

$\mathbf{w}(n)=\mathbf{w}(n-1)+\mu\left(\mathbf{X}(n) \mathbf{X}^{H}(n)\right)^{-1} \mathbf{X}(n) \mathbf{e}^{*}(n)$

如果我们令：$\mu =1$, $R(n)=\frac{1}{L}X(n)X^H(n)$, $r(n)=\frac{1}{L}X(n)y^*(n)$并结合公式17，可得

$\begin{aligned} \mathbf{w}(n) &=\mathbf{w}(n-1)+\left(\mathbf{X}(n) \mathbf{X}^{H}(n)\right)^{-1} \mathbf{X}(n)\left(\mathbf{y}^{*}(n)-\mathbf{X}^{H}(n) \mathbf{w}(n-1)\right) \\ &=\mathbf{w}(n-1)+\mathbf{R}^{-1}(n) \mathbf{r}(n)-\mathbf{R}^{-1}(n) \mathbf{R}(n) \mathbf{w}(n-1) \end{aligned}$

等价于RLS 算法

• 优点：收敛速度和计算复杂度介于 LMS 和RLS 之间

% 参考：
% https://github.com/rohillarohan/System-Identification-for-Echo-Cancellation
% https://github.com/3arbouch/ActiveNoiseCancelling/blob/c29ccfcd869657a5f58e1ce7755fe08c3a195bd9/ANC/BookExamples/MATLAB/APLMS_AEC_mono.m
function [e,y,w] = myAP(d,x,mu,M,L,psi)
    % input: 
    % d   -- 麦克风信号
    % x   -- 远端语音
    % mu  -- 步长
    % M   -- 滤波器阶数 40
    % L   -- 列数 2
    % psi -- 1e-4
    %
    % Outputs:
    % e - 输出误差 the output error vector of size Ns
    % y  - 输出系数 output coefficients
    % w - 滤波器参数 filter parameters
    
    d_length = length(d);
    if (d_length <= M)  
        print('error: 信号长度小于滤波器阶数！');
        return; 
    end
    if (d_length ~= length(x))  
        print('error: 输入信号和参考信号长度不同！');
        return; 
    end
    
    XAF = zeros(M,L);
    w1 = zeros(M,1);  % 滤波器权重 (40, 1)
    y = zeros(d_length,1);  % 估计的近端语音
    e = zeros(d_length,1);  % 误差
    
    for m = M+L:d_length    % 采样点数
        for k = 1:L % 列数
            XAF(:,k) = x(m-k+1:-1:m-k+1-M+1);   % (40,2)
        end
%         y(m) = (XAF'*w)'*(XAF'*w);        % 不太确定是不是这样(40,2)'*(40,1)=(2,1)
        E = d(m:-1:m-L+1) -XAF'*w1;    % (2,1)-(2,1)
        w1 = w1 + mu*XAF*inv((XAF'*XAF + psi*eye(L)))*E;
        w(:,m) = w1;        % (40, 73113)
        e(m) = E(1)'*E(1);
    end
end

MATLAB代码实现AP

稀疏类自适应算法(PNLMS)

通过对回声路径模型的分析，发现回声能量中较活跃系数均在时域聚集，且比重很小，其数值只有很少不为零的有效值，大多数都是零值或者接近零值，这就是回声路径具有的稀疏特性（基于时域对信号进行处理的一种特性）

比例归一化最小均方 PNLMS

回声脉冲响应具有稀疏特性，其系数绝对值只在短时间内大于零，其余则无限趋近于零。上千阶的滤波器系数中实则只有几百个发挥实际作用，且决定着回声路径估计的准确性。因此，在回声路径具有长脉冲响应的情况下，基本的 LMS 及 NLMS 算法的收敛速度一般不能达到要求，必须研究出一种拥有较快收敛速度的自适应滤波算法，在此背景下，根据回声路径的稀疏性，Duttweiler引入了比例自适应的思想，提出了比例归一化最小均方算法(Proportionate Normalized Least Mean Square, PNLMS)，按比例分配滤波器的权值向量大小，该算法对回声消除的发展具有非常重要的意义。 

PNLMS 算法拥有初始收敛速度的优势，尤其在回声路径具有稀疏性条件下，其收敛速度则更加明显。PNLMS算法采用与滤波器抽头稀疏成正比的可变步长参数来调整算法收敛速度，利用其抽头稀疏的比例值来判断当前权重稀疏所属的活跃状态，根据状态的不同，所分配的步长大小也有所差异，活跃抽头系数分配较大的步长参数，这样可以加速其收敛，而不活跃的抽头系数则相反，通过分配其较小的步长参数来提高算法的稳态误差。每个滤波器抽头被分别赋予了不同估计值，算法的稳态收敛性得到了明显改善。滤波器权值更新公式可表示为：

$\mathbf{w}(n+1)=\mathbf{w}(n)+\frac{\mu \mathbf{G}(n+1)}{\delta+\mathbf{x}^{T}(n) \mathbf{G}(n+1) \mathbf{x}(n)} \cdot e(n) \mathbf{x}(n)$

式中$G(n)$为$N$维的对角增益矩阵，以对滤波器各权值系数分配步长权重，表示为

$\mathbf{G}(n+1)=\operatorname{diag}\left\{g_{1}(n+1), g_{2}(n+1), . ., g_{N}(n+1)\right\}$

其中

$g_{k}(n+1)=\frac{\gamma_{k}(n+1)}{\frac{1}{N} \sum_{i=1}^{N} \gamma_{i}(n+1)}, 1 \leq k \leq N$

式中

$\gamma_{k}(n+1)=\max \left\{\rho \cdot \max \left[\delta_{p},\left|w_{1}(n)\right|,\left|w_{1}(n)\right|, \ldots,\left|w_{N}(n)\right|\right],\left|w_{k}(n)\right|\right\}$

$\delta_p$ 为在滤波器初始阶段调整更新参数的因子，$\rho$的取值将影响算法的整体收敛速度，较大的$\rho$值会使得收敛速度减缓，但也不宜取值过小，一般取$\frac{5}{M}$。若将$\rho$取值为 1，则PNLMS 将变为 NLMS 算法。

优点：使算法对于稀疏的回声路径，在初始阶段拥有快速的收敛速率，与此同时降低了稳态误差，

缺点：

1. 比例步长参数的选择引入了定值，这会导致估计误差累积，最后可能出现算法后期收敛速度慢、不能及时收敛的情况
2. 相比于NLMS算法增加了算法的复杂度
3. 另外在回声路径非稀疏情况下，收敛速度会变得比NLMS算法更慢

改进的PNLMS

PNLMS 算法的收敛速度提升是以回声路径脉冲响应具有高稀疏程度为前提条件的，因此在一些低稀疏度的场景下，其效果甚至没有 NLMS 算法表现好，且随着迭代的进一步进行，其收敛速度会发生显著减缓的情况。为了解决这一问题，P.A.Naylor提出了改进的IPNLMS（Improved IPNLMS，IIPNLMS）算法。结合了PNLMS 和 NLMS 算法的优势，IPNLMS 算法在二者的基础上采取了可调整算法比重的参数$\theta$ ：

$公式1：g_{k}(n+1)=\frac{1-\theta}{2}+(1+\theta) \frac{N\left|w_{k}(n)\right|}{2\|\mathbf{w}(n)\|+\varepsilon}, 1 \leq k \leq N$

$\varepsilon$为一个小的正常数，式中 

$\|\mathbf{w}(n)\|=\sum_{k=1}^{N}\left|w_{k}(n)\right|$

则公式1右半部分第一项可以看成是 IPNLMS 算法的 NLMS 成分，而第二项可以是PNLMS 成分。当$\theta =1$时 IPNLMS 算法退化为 PNLMS 算法，当$\theta =-1$时 IPNLMS 算法退化为 NLMS 算法。IPNLMS 算法内部包含了 NLMS 算法和 PNLMS 算法的成分，且由于调整算法比重参数的存在使得其内部两种算法的成分达到协调，因此使得其收敛速度保持在适中的状态。

PNLMS++

PNLMS ++算法，在每个采样周期内，通过将NLMS算法和PNLMS算法之间进行交替来实现收敛速度方面的提升。但是PNLMS ++算法只在回声路径稀疏或高度非稀疏的情况下表现良好。

CPNLMS

CPNLMS 算法（Composite PNLMS）通过比较误差信号功率和已设阈值大小，来判断算法采用 PNLMS 还是 NLMS 算法

• 缺点：阈值的选取往往和实际环境有关，因此该算法在实际应用中不能通用

$\mu$准则MPNLMS

$\mu$准则MPNLMS算法：为了解决PNLMS算法后程收敛速度慢的问题，将最速下降法应用到PNLMS算法中，使用$\mu$准则来计算比例因子（使用对数函数代替PNLMS算法中的绝对值函数）

• 优点：使得滤波器权值向量的更新更加平衡，提升了PNLMS算法接近稳态时期的收敛速率
• 缺点：复杂性加大

MPNLMS

MPNLMS算法引入最优步长控制矩阵，均衡了滤波器大小系数之间的更新，MPNLMS 中 P 步长的计算方式克服了 PNLMS 过分注重大系数收敛忽略小系数收敛的缺点，修正了 PNLMS 算法收敛后期速度慢的缺点，但

• 缺点：MPNLMS算法的运算过程中包含对数的计算，因此算法的运算复杂度相对较高。

SPNLMS

SPNLMS算法，将MPNLMS滤波器更新过程中的对数函数关系简为两段折线形式的形式。

SPNLMS的P步长计算函数$F(|w_l(k)|)$可以描述为：

$F(|w_l(k)|)=\left\{\begin{matrix} 600*|w_l(k)|&&|w_l(k)\leq 0.005|\\ 3&&其他 \end{matrix}\right.$

• 优点：相对于MPNLMS为减小运法复杂度
• 缺点：但收敛速度有所下降

改进的 MPNLMS算法，采用多个分段函数来近似 MPNLMS 的对数函数，从而降低算法的复杂度。

改进的 SPNLMS算法，通过控制步长控制矩阵迭代的频率降低算法复杂度。收敛速度也所下降

以上对 MPNLMS 的改进算法在减小算法运算复杂度的情况下损害了收敛性能和稳定性能。

稀疏控制(Sparse Control, SC)

有时候回声路径的稀疏程度会根据温度、压力以及房间墙面的吸声系数等等因素而产生变化，这就需要一种能够适应不同稀疏度变化的AEC算法。

稀疏控制比例回声消除算法（SC-PNLMS、SC-MPLNMS、SC-IPNLMS），使用新的稀疏控制方法动态地适应回声路径的稀疏程度，使得算法在稀疏和非稀疏的回声路径条件下都有良好的性能表现，体现了稀疏控制类的自适应滤波算法能够提高算法对回声路径稀疏程度的鲁棒性。

小总结：通过上述对 PNLMS 算法的分析，可知导致 PNLMS 算法整体收敛速度慢的原因主要是大系数和小系数之间收敛的不均衡。尽管很多学者针对此缺陷提出了修正算法，如 PNLMS++、CPNLMS 等，但是 PNLMS 忽略小系数收敛的缺陷并未从根本上得到改善，因此这些改进算法的效果不是十分理想。Deng通过对滤波器系数收敛过程的定量分析，在权系数更新过程中推导出最优步长的计算方式，提出一种改进的算法——MPNLMS 算法。MPNLMS 中 P 步长的计算方式克服了 PNLMS 过分注重大系数收敛忽略小系数收敛的缺点，修正了 PNLMS 收敛后期速度慢的缺陷。

一种新的改进型 PNLMS 算法，因 PNLMS 算法只注重大系数更新，忽视小系数收敛，致使算法在收敛后期速度下降，因此在 P 步长引入的同时必须注意小系数的更新。MPNLMS 算法通过建立 P 步长与当前滤波器权系数的函数关系，在一定程度上解决了 PNLMS 算法后期收敛速度慢的问题，但滤波过程包含了对数运算，不利于系统的实时实现。

​ 本文通过定量分析滤波过程，并考虑到大、小系数的收敛，建立了一种新的 P步长与当前滤波器系数之间的映射关系，降低了算法的运算复杂度。 该改进算法以 PNLMS 算法为基础，通过改变收敛过程来克服 PNLMS 算法的缺陷

子带自适应滤波器(SAF)

在声学回声消除应用中，远端输入语音信号的相关性较高，然而，传统的方法是基于“信号的无关性”假设的，传统的全带LMS 和NLMS 等计算复杂度低的随机梯度算法很难满足系统对收敛速度的要求。

远端语音信号相关性有两层含义：

• 时域上：它表征语音信号相关矩阵特征值的扩散度
• 频域上：它表征远端语音信号的频谱动态范围

一般来说，语音信号相比白色信号，前者明显有更大的频谱动态范围，即更大的信号相关性。因此，可以 通过降低输入信号的相关性来加快算法收敛速度，但是行之有效的一种方法就是是 子带自适应滤波算法，子带结构是基于频域对信号进行的一种处理（节省计算量、提高收敛速度）。

子带自适应滤波器（subband adaptive filter，SAF）：将相关信号通过滤波器组分割成近似无关的各个子带独立信号(子带分割)。然后对子带信号进行多速率抽取来获得采样信号，再进行信号的自适应处理。为研究子带自适应滤波器，首先需要了解 多速率信号抽取系统 和 滤波器组 。

多速率系统

用于 子带自适应滤波器 的多速率抽取系统 有 下采样 和 上采样 两种，主要通过 抽取 和 插值 方法来使系统获得不同采样率。输入信号经过 N 个滤波器分频后的总采样点数是原信号的 N 倍，大幅度提高的采样数增加了计算量。设采样因子是$K$，通过保留信号的$K$倍采样点可将采样速率从$f$减小为$\frac{f}{K}$，这样就降低了自适应算法的计算量。反之，上采样通过在信号相邻点间插入 K-1个 0 则使原采样率从$f$增大为$Kf$。

下采样的实现通过如下所示的抽取器和延时器。

下采样的时域表达为$x_D(n)=x(mK)$，$n$是块序列号，其频域表达式为：

$X_{D}\left(e^{j w}\right)=\frac{1}{K} \sum_{n=0}^{K-1} X\left(e^{j \frac{w-2 \pi n}{K}}\right)$

此过程可理解为将$x(n)$的频率进行$K$倍扩展，并进行$2\pi$的周期延拓，即生成了$x_D(n)$频谱。

上采样通过在相邻采样点间加入$K-1$个零点，可得上采样信号$x_I(n)$表达式为：

$x_{I}(n)=\left\{\begin{array}{cc} x\left(\frac{n}{K}\right) & n=m K, m \in Z \\ 0 & \text { 其它 } \end{array}\right.$

而频域表达式$X_I(e^{jw})$则是将$x(n)$的频率进行$K$倍压缩，再以$2\pi$为周期进行延拓得到：

$X_{I}\left(e^{j w}\right)=X\left(e^{j w K}\right)$

滤波器组

信号子带分割通过 滤波器组 实现。滤波器组是由一系列的带通滤波器组成，主要包含以下环节：1)分析滤波器组，2)抽取，3)插值，4)综合滤波器组

分析滤波器组将数字信号分割后抽取成多个子带信号，经过信号处理后，综合滤波器组再对子带信号进行插值和滤波相加而恢复成原来的信号。

上图是一个N个子带的滤波器组，其中$H_i(z)$和$G_i(z)$分别是分析滤波器组和综合滤波器组的传递函数， $i=0,1,...N-1$，代表第$i$个子带。

分析滤波器组将输入信号$x(n)$分割为子带信号$x_i(k)$，每个子带信号占据输入信号频带的一部分。综合滤波器组将$N$个子带信号重建为输出信号$y(n)$。

传递函数均为带通滤波器，带宽一样，且中心频率是均匀分布的，滤波器覆盖全部频域，无重叠部分。采用 N 通道的均匀滤波器组分解全带信号，每个子带信号$X_i(z)$只包含全带信号$\frac{1}{N}$频带，因此子带信号可以用全带信号的1/ $\frac{1}{N}$倍采样率抽取，且能保留全部的原始信号信息。

如果滤波器组的采样因子恰好等于子带数，即$K=N$，这个滤波器组就称为临界采样滤波器组。临界采样使用 N 个采样子带信号保留了有效地采样率，每个子带信号，$X_{i}(x)$的采样率是全带信号$x(n)$采样率的$\frac{1}{N}$，因此所有子带信号样值数与全带信号样值数相等。

在综合过程中，子带信号$x_i(k)$采用相同的插值因子$K$，然后由综合滤波器组 合并为全带信号。因此，在重建全带信号$y(n)$时原来的采样率就被恢复了。

信号经过分析滤波器组后子带信号表达式为：

$X_{i}(z)=H_{i}(z) X(z)$

采样后的子带信号可以写为：

$X_{i, D}(z)=\frac{1}{N} \sum_{l=0}^{N-1} H_{i}\left(z^{1 / N} W_{N}^{l}\right) X\left(z^{1 / N} W_{N}^{l}\right), i=0,1, \ldots, N-1$

先将采样后的子带信号，$X_i(z)$经过相同因子 N 插值，再经过综合滤波器组合并为全带信号$Y(z)$。综合滤波器组的输出信号可以写为

$\begin{aligned} Y_{i}(z) &=G_{i}(z)\left[\frac{1}{N} \sum_{l=0}^{N-1} H_{i}\left(z W_{N}^{l}\right) X\left(z W_{N}^{l}\right)\right] \\ &=\frac{G_{i}(z)}{N}\left[H_{i}(z) X(z)+\sum_{l=1}^{N-1} H_{i}\left(z W_{N}^{l}\right) X\left(z W_{N}^{l}\right)\right] \\ &=\frac{G_{i}(z)}{N}\left[X_{i}(z)+\sum_{l=1}^{N-1} X_{i}\left(z W_{N}^{l}\right)\right] \end{aligned}$

明显地 ，$Y_i(z)$是 由 期 望 部 分$X_i(z)=H_i(z)X(z)$和其它频移部分$X_{i}\left(z W_{N}^{l}\right)=H_{i}\left(z W_{N}^{l}\right) X\left(z W_{N}^{l}\right), l=1,2, \ldots, N-1$共同组成，$X_{i}\left(z W_{N}^{l}\right)$称为$l$阶混叠部分，它是由于子带信号$X_i(z)$采样率偏移造成的。在实际应用中，滤波器存在过渡带和有限的阻带衰减，因此临界采样存在混叠的问题。正交镜像滤波器组(QMF)是无混叠的滤波器组：当插值信号$Y_i(z)$经过综合滤波器组合并为全带信号$Y(z)$时，混叠部分可以通过子带相互抵消。即使使用无混叠的 QMF 滤波器组，采样子带信号$X_i(z)$依然存在着混叠。当然，混叠信号可以利用理想的分析滤波器组消除，但是在现实中无法实现。

可以看出$X_i(zW_N^l0$是综合滤波器组的输出信号$Y_0(z),...,Y_{N-1}$中共用的部分。通过把共用项组合起来， N 通道的 QMF 滤波器组的输出结果

$\begin{aligned} Y(z) &=\sum_{i=0}^{N-1} Y_{i}(z) \\ &=\sum_{i=0}^{N-1}\left[\frac{1}{N} \sum_{l=0}^{N-1} F_{i}(z) H_{i}\left(z W_{N}^{l}\right) X\left(z W_{N}^{l}\right)\right] \\ &=\sum_{i=0}^{N-1}\left[\frac{1}{N} \sum_{l=0}^{N-1} F_{i}(z) H_{i}\left(z W_{N}^{l}\right)\right] X\left(z W_{N}^{l}\right) \\ &=\sum_{i=0}^{N-1} A_{l}(z) X\left(z W_{N}^{l}\right) \end{aligned}$

其中$A_{l}(z)=(1 / N) \sum_{l=0}^{N-1} F_{i}(z) H_{i}\left(z W_{N}^{l}\right)$。上式表明 QMF 滤波器组的输出信号Y(z)是输入信号 X(z)和$X(zW_N^l)$的加权和，即

$\begin{aligned} Y(z) &=\left[(1 / N) \cdot \sum_{l=0}^{N-1} F_{i}(z) H_{i}(z)\right] X(z) \\ &=T(z) X(z) \end{aligned}$

其中$T(z)=1 / N \cdot \sum_{l=0}^{N-1} F_{i}(z) H_{i}(z)$为滤波器组失真的传递函数。$T(z)$是衡量滤波器输出的幅值和相位失真情况，假如T(z)在频谱上幅值固定且相位是线性的，就称之为完全重建滤波器组。

在实际应用中，可以使用余弦调制滤波器组来实现。学者们对余弦调制滤波器组的理论和设计进行了深入研究。余弦调制滤波器组的分析滤波器组和综合滤波器组都是 N 通道的临界采样滤波器组，它们可以用酉矩阵分别表达为$E(z)$和$R(z)$。分析滤波器组$E(z)$的完全重建系统可以通过仿酉特性实现：

$\mathbf{R}(z)=z^{-K+1} \tilde{\mathbf{E}}(z)$

其中，K-1是多相分量矩阵的阶数。仿酉条件意味着$E^{-1}(z)=\tilde{E}(z)$。综合滤波器组中的多相分量可以用上式实现，即R_{r, i}(z)=z^{-K+1} \tilde{E}_{r, i}(z)。

N 通道的余弦调制滤波器组的分析滤波器组和综合滤波器组都以低通滤波器为原型滤波器，这些低通滤波器的截止频率为$\frac{\pi}{2N}$。余弦调制的分析滤波器表示为

$H_{i}(z)=\alpha_{i} P\left[z W_{2 N}^{(i+0.5)}\right]+\alpha_{i}^{*} P\left[z W_{2 N}^{-(i+0.5)}\right], i=0,1, \ldots, N-1$

其中，$W_{2 N}=e^{-j \pi / N}$是旋转因子，$\alpha_i$定义为 

$\alpha_{i}=\exp \left\{j\left[\theta_{i}-\frac{\pi}{N}(i+0.5)\left(\frac{L-1}{2}\right)\right]\right\}, \theta_{i}=(-1)^{i} \frac{\pi}{4}$

L 是原型滤波器的长度。综合滤波器可以用分析滤波器的时域反转的形式表示

$F_{i}(z)=z^{-L+1} H_{i}\left(z^{-1}\right)$

这样就可以减少原型滤波器的设计。

余弦调制分析滤波器 $H_i(z)$由两个复数滤波器组成，即$P\left[z W_{2 N}^{(i+0.5)}\right]$和$P\left[z W_{2 N}^{-(i+0.5)}\right]$，它们都由 P(z) 的复数调制产生，$P\left[z W_{2 N}^{(i+0.5)}\right]$的脉冲响应对应于$P\left[z W_{2 N}^{-(i+0.5)}\right]$的复共轭。因此，分析滤波器组$H_i(z)$可以用实数脉冲响应表示

$h_{i}(n)=2 p(n) \cos \left[\frac{\pi}{N}(i+0.5)\left(n-\frac{L-1}{2}\right)+\theta_{i}\right]$

类似的，综合滤波器组$F_i(z)$用实数脉冲响应表示为 

$f_{i}(n)=2 p(n) \cos \left[\frac{\pi}{N}(i+0.5)\left(n-\frac{L-1}{2}\right)-\theta_{i}\right]$

通过优化原型滤波器的设计来满足一系列的预定限制，这样就能完全重构信号。不同的限制和优化方法可以推导出两种不同的余弦调制滤波器组类，即伪 QMF 余弦调制滤波器组和仿酉余弦调制滤波器组。

子带自适应算法结构

 在传统的 SAF中，子带自适应算法都是以最小化子带误差信号为目标的，这样基于局部目标函数误差的最小化不一定是全局误差能量最小化。而分析滤波器组在子带切割和综合滤波器组重建全带信号时皆会引入时延，在AEC 应用中，这样的时延会使包含近端语音的全带误差信号传到远端，为了消除时延的影响，无延时子带闭环结构系统以全局误差能量最小化为约束条件来调整滤波器系数。最后，确保自适应滤波算法能够收敛到最佳的滤波器系数。

下图所示是闭环无延时系统的结构图，无延时的子带自适应算法就是通过该结构来实现的。

该结构可以消除信道延迟带来的影响，将子带更新的抽头滤波器系数$W_i(z)$映射到全带自适应滤波器$W(z)$上，全带滤波器的表达式为：

$\hat{\boldsymbol{W}}(z)=\sum_{i=0}^{K-1} \hat{\boldsymbol{W}}_{i}\left(z^{K}\right) z^{i}$

从无延迟子带系统结构图中可以看出，输入信号直接经过实际房间脉冲响应后就可以得到期望信号$d(n)$，而预测信号$\hat{y}(n)$则还要经过分析滤波器组和自适应滤波器，这会导致两者间的时延，引入$\Delta$，作为一个补偿参数，其值的大小可用如下表达式求得：

$公式1：\Delta=\left\lfloor\frac{2 N C-1}{K}\right\rfloor$

式中，$[\ ]$是求整数运算，C 是一个常数。子带自适应滤波器系数矩阵$w(n)$加入补偿因子后公式为：

$公式2：\hat{\boldsymbol{W}}_{i}(z)=\sqrt{K} \boldsymbol{w}_{i}(z) z^{-\frac{K-(i-1)}{K}+\Delta+1}$

$w_i(z)$由$w(n)$经傅里叶变换得到，将$w_i(z)$代入:2后再代入式1后就完成了子带自适应滤器系数矩阵到全带系数矩阵$\hat{W}(n)$的映射，$\hat{W}(n)$与输入信号$x(n)$相乘得到$\hat{y}(n)$，全局误差信号$e(n)$可以由期望信号$d(n)$和预测信号$\hat{y}(n)$相减得到。

• 优点：改善了全带自适应滤波算法在相关信号条件下的收敛速率
• 缺点：
  • 但其稳态误差由于输出时存在的混叠分量而显著升高
  • 当采用正交镜像滤波器组时，虽然可以通过子带系统将混叠部分相互抵消掉，但在现实中却无法实现

归一化子带自适应滤波 NSAF

问题：针对SAF算法稳态误差较高的问题

解决：提出了 基于最小扰动原理提出了归一化的SAF（normalized SAF，NSAF）算法。

优点：由于SAF类算法固有的解相关特性，NSAF在处理相关输入信号时比全带的NLMS收敛速度快，而且计算成本与NLMS不相上下

由上图可知$x(n)$与$d(n)$经过分析滤波器组分割成数量为 N 的子带信号$x_i(n)$与$d_i(n)$，子带输出信号$y_i(n)$由子带输入信号$x_i(n)$通过估计的滤波器$\hat{w}(k)$进行滤波而获得。接着$y_i(n)$和$d_i(n)$分别经过$N$倍抽取得到 $y_{i,D}(k)$和$d_{i,D}(k)$。其中$n$表示原始序列，$k$表示抽取的序列。如图所示，子带误差信号表示为：

$e_{i, D}(k)=d_{i, D}(k)-\mathbf{x}^{\mathrm{T}}(k) \hat{\mathbf{w}}(k), i=0,1, \ldots, N-1$

 NSAF 算法的滤波器权值系数更新迭代公式：

$\hat{\mathbf{w}}(k)=\hat{\mathbf{w}}(k-1)+\mu \sum_{i=0}^{N-1} \frac{e_{i, D}(k) \mathbf{x}_{i}(k)}{\delta+\left\|\mathbf{x}_{i}(k)\right\|_{2}^{2}}$

式中，$\delta$为正则化参数。NSAF 针对相关信号的收敛性能比 NLMS 强的主要原因是 SAF类算法的固有去相关特性。此外，对于 AEC 等高阶的自适应滤波器应用，NSAF 的计算复杂度几乎与 NLMS 相同。另外，当 NSAF 算法只有一个子带时，其本质相当于NLMS 算法。NSAF 算法解决了输入为相关信号的全带自适应滤波算法收敛性能降低的问题。

近几年研究人员为了能够提升AEC算法的收敛性能和稳态性能，在NSAF的基础上结合全带自适应滤波算法的成比例理论，提出了几种改进的NSAF算法，例如不同形式的变步长因子NSAF以及变正则化参数NSAF。为了在识别稀疏回声路径时快速收敛，文献[22,23]将NLMS算法中的成比例思想以类比的方式融合到NSAF算法中，提出了比例NASF（proportionate NASF，PNSAF）算法和μ准则PNSAF（μ-law PNASF，MPNSAF）算法。

因为子带结构中存在混叠分量问题

1. Keermann于1988年利用采样滤波器组技术消除了混叠现象，但是此举增加了算法复杂度。
2. 相邻子带间留安全频带，缺点：引入了空白频带，降低了信号质量。
3. 重叠子滤波器补偿法，缺点：因为交叉项而增加了运算量，还降低了收敛速度。

2004年K. A. Lee 和 W. S. Gan提出了基于最小扰动原理的多带结构式自适应滤波器(Multiband Structured SAF，MSAF)算法，并给出了自适应滤波器抽头系数的更新方程。该结构完全不存在滤波器输出端的混叠分量问题。

多带自适应滤波器

子带自适应滤波器中每个子带单独使用一个子滤波器。该结构会导致输出端产生混叠分量，解决此问题的传统方法多以降低信号的质量或增大稳态误差为代价，Lee 和 Gan在 2004 年提出了一种全新多带结构。不同于子带滤波器在每个子带都使用不同的滤波器，多带结构的每个子带使用相同的全带滤波器，这很好地克服了输出端存在混叠分量的问题。

假设全带自适应滤波器$W(k,z)$长度为 M ， N 个子带自适应滤波器只占原始输出频率的一部分，输出信号为

$y_{i}(n)=\sum_{m=0}^{M-1} w_{m}(k) x_{i}(n-m), i=0,1, \ldots, N-1,$

其中，$\boldsymbol{W}(k, z)=\sum_{m=0}^{M-1} w_{m}(k) \mathbf{z}^{-m}$ 。信号经过采样因子$N$抽取后，第$i$个子带信号为：

$\begin{aligned} y_{i, D}(k) &=\mathrm{y}_{i}(k N) \\ &=\sum_{m=0}^{M-1} w_{m}(k) x_{i}(k N-m) \\ &=\boldsymbol{w}^{T}(k) \boldsymbol{x}_{\mathrm{i}}(k) \end{aligned}$

其中，$\boldsymbol{w}(k)=\left[w_{0}(k), w_{1}(k), \ldots, w_{M-1}(k)\right]^{\mathrm{T}}$是全带自适应滤波器$W(k, z)$的权重矢量，$x_i(k)$表示第$i$个子带输入矩阵：

$\begin{array}{r} \boldsymbol{x}_{\mathrm{i}}(k)=\left[x_{i}(k N), x_{i}(k N-1), \ldots, x_{i}(k N-N+1)\right. \\ \left.x_{i}(k N-N), \ldots, x_{i}(k N-M+1)\right]^{T} \end{array}$

每个独立自适应运算的子带误差信号如下：

$e_{i, D}(k)=\mathrm{d}_{i, D}(k)-\boldsymbol{w}^{\mathrm{T}}(k) \boldsymbol{x}_{i}(k) \text { for } i=0,1, \ldots, N-1$

N 个子带误差信号向量：$\boldsymbol{e}_{\mathrm{D}}(k) \equiv\left[e_{0, \mathrm{D}}(k), e_{\mathrm{l}, \mathrm{D}}(k), \ldots, e_{N-1, \mathrm{D}}(k)\right]^{T}$写出如下形式：

$\boldsymbol{e}_{\mathrm{D}}(k)=\boldsymbol{d}_{\mathrm{D}}(k)-\boldsymbol{X}^{\mathrm{T}}(k) \boldsymbol{w}(k)$

其中，$M*N$子带信号矩阵$X(k)$ 和 $N*1$ 期望响应矢量$d_D(k)$分别为：

$\boldsymbol{X}(k) \equiv\left[\boldsymbol{x}_{0}(k), \boldsymbol{x}_{1}(k), \ldots, \boldsymbol{x}_{N-1}(k)\right]$

$\boldsymbol{d}_{\mathrm{D}}(k) \equiv\left[d_{0, \mathrm{D}}(k), d_{1, \mathrm{D}}(k), \ldots, d_{N-1, \mathrm{D}}(k)\right]^{\mathrm{T}}$

频域分块LMS(FDAF)

在讲频域分块LMS之前建议大家回顾一个时域分块LMS算法

时域分块LMS：

$e(n+m)=d(n+m)-X^{T}(n+m) \mathbf{w}(n)$

$w(k+1)=w(k)+\mu \sum_{i=0}^{L-1} X(k L+i) e(k L+i)$

Block  LMS的误差计算 和 权重更新公式中，

$y(n+m)=X^{T}(n+m) w(n)$：是输入向量与滤波器系数向量的线性卷积

$\hat{\nabla}(k)=-\sum_{i=0}^{L-1} X(k L+i) e(k L+i)$：是误差信号与输入向量的线性相关

对于线性卷积和线性相关的运算量较大，特别是在

当回声路径很长且复杂，并且回声延迟较高时，时域自适应滤波算法中的线性卷积和线性相关运算量较大，导致计算复杂度升高，我们更愿意把这两个信号变换到频域，通过频域相乘的方式来取代时域复杂度相当高的卷积或相关运算。所以我们把时域自适应滤波算法，衍生到频域自适应滤波算法。

预备知识：线性卷积（相关）和圆周卷积（相关）之间的关系

1. 一般的，如果两个有限长序列的长度为$N_1$和$N_2$，且满足$N_1\geq N_2$，则有圆周卷积的后 $N_1-N_2+1$个点，与线性卷积的结果一致。
2. 一般的，如果两个有限长序列的长度为$N_1$和$N_2$，且满足$N_1\geq N_2$，则有圆周相关的前 $N_1-N_2+1$个点，与线性相关的结果一致。
3. 时域中的圆周卷积对应于其离散傅里叶变换的乘积
4. 时域中的圆周相关对应于其离散傅里叶变换共轭谱的乘积

因此我们需要将 分块LMS自适应滤波算法中的线性卷积和线性相关 通过快速傅里叶变换(FFT)转换到频域来实现。这样实现的分块LMS自适应滤波算法称为频域块LMS自适应滤波算法(FDAF，Frequency-Domain Block Least Mean Square Adaptive Filter)。FDAF算法将长度为L的自适应滤波器分成FFT长度的整数倍个子块，对输入信号的每个子块进行频域内的LMS算法。

第一步：计算线性卷积

$y(n+m)=X^{T}(n+m) w(n)$

利用FFT计算线性卷积的方法有：

• 重叠存储法（overlap save method，更常用）
• 重叠相加法（overlap add method）

我们这里以overlap save method为例，为了确保能得到$N$个点的线性卷积输出信号，我们至少要保证有N个点的线性卷积和圆周卷积的结果一致（预备知识）

$N_1-N_2+1=N$

由于$N_1\geq N_2$ (输入信号长度通常大于滤波器的阶数)，且$N_2=N$ (滤波器的阶数为N)，那么要求每次参与运算的输入信号长度$N_1$至少为$2N-1$，为了计算FFT方便，我们令输入信号的长度为：$N_1=2N$，那么我们FFT的长度也为$2N$

为了构造长度为$2N$的数据，我们需要在每个$N$阶滤波器后面$N$补零

要 求线性卷积(预备知识1)，我们就需要 求圆周卷积后$N_1-N_2+1$个点，根据预备知识3，我们只需要求 离散傅里叶变换的乘积 就能得到圆周卷积的结果，接下来我们分别计算 输入信号向量 和 滤波器系数向量 的FFT：

$X(k)=\operatorname{diag}\{F[x(k N-N), \ldots, x(k N-1),| x(k N), \ldots, x(k N+N-1))]\}$

$W(k)=F\left[W^{T}(k), |0, \ldots, 0\right]^{T}$

频域相乘

$\mathbf{Y}(k)=\mathbf{X}(k) \mathbf{W}(k)$

则，$N$点线性卷积输出信号$y(k)=[y(kN),y(kN+1),...,y(kN+N-1)]$，就等于$Y(k)$的傅里叶逆变换的后N个点：$y(k)=last\ N\ samples\ of\ F^{-1} Y(k)$

图: overlap save 切片法

第二步：计算线性相关

$\hat{\nabla}(k)=-\sum_{i=0}^{L-1} X(k L+i) e(k L+i)$

根据预备知识2，4可知，需要求线性相关，我们可以通过获得 圆周相关 来获得。因此我们需要求输入信号的共轭谱与误差信号谱的乘积。

输入信号$X(k)$在上一步处理卷积运算是已经求得。

那么，剩下的工作就是，将误差向量$e(k)=[e(kN),e(kN+1),...,e(kN+N-1)]$也扩展到$2N$长度，因为是求相关，我们需要在误差向量前面补0，然后经过FFT：$E(k)=F\left[0,0, \ldots, 0, e^{T}(k)\right]^{T}$

频域相乘：

$\boldsymbol{\Delta}(k)=X^{*}(k) E(k)$

最后，将梯度向量进行傅里叶逆变换$\vec{\nabla}(k)=F^{-1} \Delta(k)$，取前$N$个点，就是我们求的线性相关。

第三步：滤波器系数更新

$W(k+1)=W(k)+\mu F\left[\vec{\nabla}^{T}(k), 0,0, \ldots, 0\right]^{T}$

注意：

• 第一，滤波器系数直接在频域更新，所以需要将梯度向量再次变换到频域；
• 第二，由于滤波器系数向量后面补了N 个零 ，为了保证结果的正确性梯度向量也需要在后面补 N 个零。

基于重叠存储法的频域块LMS自适应滤波算法的信号流程图

• 优点：
  • 降低了时域自适应滤波器的计算复杂度
  • 提高了收敛速度
• 缺点：
  • 增加了延迟（需要通过频域滤波器延迟期望信号）
  • 增加了内存需求（需要同时存储激励信号和期望信号）

如前所述，回声消除等应用回波路径较长，会导致较大的延迟和存储需求。该缺点可以通过诸如多延迟自适应滤波器之类的方法来克服。在这种方法中，块大小可以小于所需的时域自适应滤波器，并且可以应用每个频点中的自适应滤波器来代替单个系数。因此，可以减轻FDAF的缺点，同时保持降低的计算复杂性和提高的收敛速度。

% 参考：https://github.com/CharlesThaCat/acoustic-interference-cancellation
function [en, yk, W] = myFDAF(d,x,mu,mu_unconst, M, select)
% 参数:
% d                输入信号(麦克风语音)
% x                远端语音
% mu                约束 FDAF的步长
% mu_unconst        不受约束的 FDAF的步长
% M                 滤波器阶数
% select;            选择有约束或无约束FDAF算法
%
% 参考:        
% S. Haykin, Adaptive Filter Theory, 4th Ed, 2002, Prentice Hall
% by Lee, Gan, and Kuo, 2008
% Subband Adaptive Filtering: Theory and Implementation
% Publisher: John Wiley and Sons, Ltd

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

AdaptStart = 2*M;       % 在获得2M个样本块后开始自适应
W = 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)
            % FFT[old block; new block]
            % Old block 包含M个先前的输入样本 (u_old)
            % New 包含M个新的输入样本 (u_new)

            % 计算滤波器估计信号
            Yk = Xk.*W;                  % 输入和权重向量的乘积(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 select == 1
                gradient = real(ifft((1./power).*conj(Xk).* Ek));   % (2M,1)
                gradient = gradient.*window_save_first_M;   % 去除后一个数据块，并且补零 (2M,1)
                W = W + mu.*(fft(gradient));    % 权重是频域的 (2M,1)
            else
                gradient = conj(Xk).* Ek;   %  (2M,1)
                W = W + mu_unconst.*gradient;    % (2M,1)
            end
            
            en = [en; error];             % 更新误差块
        end
    end
end

MATLAB代码实现FDAF

""" frequency domain adaptive filter """

import numpy as np
from numpy.fft import rfft as fft
from numpy.fft import irfft as ifft


def fdaf(x, d, M, mu=0.05, beta=0.9):
    H = np.zeros(M+1,dtype=np.complex)
    norm = np.full(M+1,1e-8)

    window =  np.hanning(M)
    x_old = np.zeros(M)

    num_block = len(x) // M
    e = np.zeros(num_block*M)

    for n in range(num_block):
        x_n = np.concatenate([x_old,x[n*M:(n+1)*M]])
        d_n = d[n*M:(n+1)*M]
        x_old = x[n*M:(n+1)*M]
 
        X_n = np.fft.rfft(x_n)
        y_n = ifft(H*X_n)[M:]
        e_n = d_n-y_n

        e_fft = np.concatenate([np.zeros(M),e_n*window])
        E_n = fft(e_fft)

        norm = beta*norm + (1-beta)*np.abs(X_n)**2
        G = X_n.conj()*E_n/norm
        H = H + mu*G

        h = ifft(H)
        h[M:] = 0
        H = fft(h)

        e[n*M:(n+1)*M] = e_n

    return e

python代码实现FDAF

学习速率如何选取

$W(k+1)=W(k)+2 \mu F\left[\vec{\nabla}^{T}(k), 0,0, \ldots, 0\right]^{T}$

我们在推导频域自适应滤波方法的时候，为了简化问题，将每一个frequency bin 上的学习速率均设为常数$\mu$。

在实际工程应用中，用的更多的一种方法是，对第$m$个 frequency bin ，利用输入信号在这个频点的功率 对学习速率$P_m(k)$进行归一化：

$\mu_{m}(k)=\frac{\mu}{P_{m}(k)}$

频点功率$P_m(k)$通常采用迭代的方式求得：

$P_{m}(k)=\lambda P_{m}(k-1)+(1-\lambda)\left|X_{m}(k)\right|^{2}$

补零与倒序

以 LMS 算法权重更新为例，

有人写的算法是前面补零

function [e,w,y] = q_buling(mu, d_length, M, mic)
    xx = zeros(M,1);
    w_updata = zeros(M,1);    % 滤波器权重
    for n = 1:d_length
        xx = [xx(2:M);mic(n)];    % 纵向拼接
        % xx(2~40)+x(1)-->xx(2~40)+x(2)-->xx(2~40)+x(3)...
        y(n) = w_updata' * xx;        % 远端回声估计(40,1)'*(40,1)=1; (73113,1)
        e(n) = mic(n) - y(n);     % 近端语音估计
        w_updata = w_updata + mu * e(n) * xx;   % (40,1)
        w(:,n) = w_updata;        % (40, 73113)
    end

有人写的算法是前面补零

function [e,w,y] = h_buling(mu, d_length, M, mic)
    xx = zeros(M,1);
    w_updata = zeros(M,1);    % 滤波器权重
    for n = M:d_length
        xx = [mic(n);xx(2:M)];
        y(n) = w_updata' * xx;        % 远端回声估计 (40,1)'*(40,1)=1; (73113,1)
        e(n) = mic(n) - y(n);     % 近端语音估计
        w_updata = w_updata + mu * e(n) * xx;   % (40,1)
        w(:,n) = w_updata;        % (40, 73113)
    end

有人写的算法是信号倒序(更常见)，（为什么要倒序呢？有人能解答一下吗？）

function [e,w,y] = daoxu(mu, d_length, M, mic)
    xx = zeros(M,1);
    w_updata = zeros(M,1);    % 滤波器权重
    for n = M:d_length
        xx = mic(n:-1:n-M+1);    % 纵向拼接  (40~1)-->(41~2)-->(42~3)....
        y(n) = w_updata' * xx;        % 远端回声估计 (40,1)'*(40,1)=1; (73113,1)
        e(n) = mic(n) - y(n);     % 近端语音估计
        w_updata = w_updata + mu * e(n) * xx;   % (40,1)
        w(:,n) = w_updata;        % (40, 73113)
    end

好像没看到正序的，那我自己写一个，下面代码写错了，有空改一下

function [e,w,y] = zhenxu(mu, d_length, M, mic)
    w_updata = zeros(M,1);    % 滤波器权重
    for n = M:d_length
        xx2 = mic(n:-1:n-M+1);    % 纵向拼接  (40~1)-->(41~2)-->(42~3)....
        y(n) = w_updata' * xx2;        % 远端回声估计 (40,1)'*(40,1)=1; (73113,1)
        e(n) = mic(n) - y(n);     % 近端语音估计
        w_updata = w_updata + mu * e(n) * xx2;   % (40,1)
        w(:,n) = w_updata;        % (40, 73113)
    end

最后画个图，看看它们的收敛情况：

clc;clear;
%%
power=2;     % 信噪比
far = wgn (160000,1,power);     % 生成高斯白噪声,作为远端语音 (160000, 1)
far = mapminmax(far, -1, 1);
noise_length = length(far);
audiowrite("noise.wav",far,16000)
mic = audioread("./add_revred.wav");        % 麦克风语音
mic = mic(1:noise_length);

%% 参数
mu = 0.5;   % 步长
M = 40;     % 滤波器阶数
d_length = length(mic);

%%
[e_1,w_1,y_1] = daoxu(mu, d_length, M, mic);
[e_2,w_2,y_2] = q_buling(mu, d_length, M, mic);
[e_3, w_3, y_3] = h_buling(mu, d_length, M, mic);
[e_4,w_4,y_4] = zhenxu(mu, d_length, M, mic);

error1_LMS=10*log10(e_1.^2);
error2_LMS=10*log10(e_2.^2);
error3_LMS=10*log10(e_3.^2);
error4_LMS=10*log10(e_4.^2);

plot(error1_LMS,':b');              % 蓝色
hold on
plot(error2_LMS,'--r');              % 蓝色
plot(error3_LMS,'-.g');              % 蓝色
plot(error4_LMS,"--c");              % 蓝色
legend({'倒序',"前补零","后补零","正序"});   % 图例

axis tight;                         % 使用紧凑的坐标轴
title('LMS算法误差曲线');            % 图标题
xlabel('样本');                     % x轴标签
ylabel('误差/dB');                  % y轴标签
grid on;                            % 网格线

我不了解 前面补零 和 后面补零的区别。为什么前面补零会更好一点呢？为什么需要样本倒序呢？

综上所述

《声学回声消除中子带和块稀疏自适应算法研究_魏丹丹》

远端输入语音信号的相关性较高，且声学回声信道的冲击响应一般只有少量的非零系数，因此是一个稀疏信道。

针对用于声学回声消除的子带和块稀疏算法进行了研究和改进，以达到提高算法跟踪性能和抗冲激鲁棒性的目的。本文的主要贡献如下：

首先，区别于传统文献中子带归一化自适应算法消除回声的方法，

我提出一种用于声学回声消除的新型子带归一化自适应滤波切换算法(LMS-NSAF)。

该算法核心思想是根据语音信号的状态不同，采用 VAD 快慢包络技术切换算法，当输入远端信号的瞬时能量值较大时，使用收敛速度快的子带 NLMS 算法，当输入信号的瞬时能量值较小时，则使用计算复杂度低的权重矢量更新公式，从而使得改进的子带 NLMS 算法在提高收敛性的同时又能降低算法的计算复杂度

基于多带结构的改进型自适应滤波切换算法NLMS-NSAF

首先远端语音信号利用包络法判别有无语音段，

然后将信号状态输出到自适应多带结构算法模块当中。

若语音区输入信号的短时能量较大，则使用收敛速度快的自适应滤波算法(NLMS)；

若语音区输入信号的短时能量较小，需考虑计算量低的算法(NSAF)，

当然，语音在无语音区时算法迭代停止。

对输入语音信号能量高低的判定是通过和阈值比较得到的。在充分考虑语音特性的情况下，切换算法实现了算法在收敛速度的优势，同时完成了同算法复杂度的优化选择。最后达到了提高滤波算法性能、降低运算量的目的。

回声路径是经典的稀疏路径，且语音信号作为远端输入时相关性较强，

总结

回声消除的挑战在于能否快速跟踪回声路径中的变化，同时又对大声干扰（例如双向通话）保持鲁棒性。这两个目标是矛盾的，因为为了快速适应回声路径的变化，系统需要具有快速收敛速度的自适应算法，这又意味着在出现双向通话时很容易发散。

在大多数系统中，实现了双向通话检测器以冻结自适应并防止发散。双向通话检测器的决策阈值通常需要对信号处理环境进行一些调整。这是一个不良的特性，因此，大多数系统都被调整为过于保守。换句话说，由于双向通话决策遗漏比双向通话虚假警报更重要，因此要进行调整以确保不会发生任何丢失的决策。这可能会导致系统产生高频率的错误警报，从而导致系统无法及时调整。

健壮的自适应算法在回声消除系统中的应用使双向通话检测器的保守性降低。由于大多数双向通话遗漏都发生在语音段的开始和偏移处，因此鲁棒的自适应算法试图使滤波器系数的自适应不受干扰，直到可以做出双向通话决定为止。由过去误差信号的电平调节的比例因子控制自适应算法的速率。因此，例如尚未收敛的系统将导致高比例因子。在双向通话的情况下，由于使用了过去的误差信号值，因此其开始会被延迟。

消除回声的两路径方法（AEC）论文进一步讨论了该主题。双向通话检测器的复杂性以及使用两组滤波器系数作为缓解快速收敛和双向通话问题的方法的目的变得更加清晰。双路径方法的问题是需要两个FIR操作，这给系统增加了内存和计算复杂度。鲁棒的自适应算法的应用可以帮助缓解过于保守的双向通话检测器的问题，而不会增加两径方法的负担。

除了居于核心地位的自适应滤波技术外，实际回音抵消技术应用系统中还包括远端信号检测、近端信号检测、舒适噪声产生、残留回波的非线性处理技术等，这些技术也有待改进。这样整个回音抵消器才能实现一个较好的回音抵消效果。

 大家可以根据实际项目中对计算复杂度、 收敛速度的要求 选择合适的自适应算法 解决自身场景的问题 。

参考

Sub-band Acoustic Echo Cancellation

Acoustic Echo Cancellation Using Subband Adaptive Filtering

Subband Adaptive Filtering: Theory and Implementation

《声学回声消除中子带和块稀疏自适应算法研究_魏丹丹》抄袭了(有比较大的嫌疑)《子带自适应滤波器及其应用_胡伟》部分内容，改了之后语句都不通

《车载免提系统降噪算法的研究以及硬件实现》__张雪

频域块LMS自适应滤波算法的研究__田超

【代码】pyaec

【代码】acoustic-interference-cancellation