主动噪声控制

主动噪声控制，又称为有源噪声控制，是一种通过发出与外界噪声相位相反且振幅相同的声波来抵消噪音的技术。ANC 的原理是物理学中的声波相消理论。通过产生与原始噪声幅值相等、相位相反的次级声波（反噪声）来抵消不需要的噪声，是传统被动噪声控制的重要补充。

前言

声波相消理论

当两列声波同时在同一媒介中传播并且在某一处相遇时，任意点上的振动将是两者引起振动的叠加，该现象就是声波干涉。如果两列声波具有相同频率及固定相位差，当其传播到同一位置时，若出现同相振动则会产生相长干涉，若出现反相振动则会弧线相消干涉。

其中初级噪声是原始噪声信号，次级噪声是扬声器生成的噪声，用来抵消原始噪声。在系统正常工作的情况下，扬声器的输出信号与原始噪声信号在特定的空间内实现干涉，要想形成稳定的干涉时，两列波需要满足如下条件

传播方向相同
相位差保持恒定
振动频率相同

它的数学表达式为

原始噪声 $p_1(t)=A\sin(\omega t+\phi_1)$
反噪声 $p_2(t)=A\sin(\omega t+\phi_2)$ ，其中为了消除噪声，则 $\phi_2=\phi_1+180^\circ$
合成声音为 $p(t)=p_1(t)+p_2(t)=0$

数字滤波器

数字滤波器 | LuosBlog

数字傅里叶分析

傅里叶分析的本质是，任何满足狄利克雷条件的信号，都可以分解为一系列不同频率、不同幅度、不同相位的正弦或余弦信号的线性叠加

时域视角：信号随时间变化的波形，描述信号在什么时候是什么样
频域视角：信号由哪些频率成分组成，每个频率成分的强度和相位，描述信号由什么组成

一个方波可以分解为基波和无数奇次谐波的叠加，叠加的谐波越多，波形越接近方波。在进行傅里叶分析时，需要先将模拟信号转换为数字信号，再进行傅里叶分析

采样：以固定时间间隔对连续信号进行取值，得到离散时间信号。为了保证无失真地恢复原信号，采样频率必须大于信号最高频率地两倍，否则回发生混叠失真
量化：将连续幅值地采样值转换为有限精度地二进制数，这会引入量化噪声
截断：计算机只能处理有限长度的信号，因此需要将无限长的离散时间信号截断为有限长度的序列

傅里叶变换：傅里叶变换 | LuosBlog

ANC 系统基本结构

控制结构分类如下

前馈控制：产生次级噪声之前就通过传感器测量初级噪声的频率以获取参考信号
反馈控制：不需要测得参考信号就产生次级噪声进行相消干涉
混合控制

前馈控制

前馈 ANC 控制通过在目标噪声源处放置参考传声器或非声学传感器的方式来直接获取参考信号，误差传感器测得的残余噪声信号连同传感器获取的参考信号均作为控制器的输入，产生并调节次级声源信号，驱动次级扬声器发出次级噪声，与初级声源产生的噪声进行相消干涉，最终使得误差传感器处声压值最小。

前馈 ANC 控制有很强的鲁棒性，不仅适用于窄带噪声信号，也可用于处理宽带噪声信号。但是有些特定的场景可能不适合安装参考传声器。

参考传声器或参考麦克风采集噪声信号 $x(n)$
ANC 控制器处理 $x(n)$ 产生反相信号 $y(n)$
次级声源播放 $y(n)$ 以抵消原始噪声
误差麦克风检测残余噪声 $e(n)$ 用于更新控制器的参数

反馈控制

反馈型 ANC 系统中没有传感器来测得参考输入信号，仅通过误差传感器获取经相消干涉后的残余噪声并将其送入到反馈控制器，进而达到调节次级声源的目的，使其发出与初级噪声幅值相等相位相反的次级噪声。

反馈 ANC 系统的结构简单，不存在单反馈问题，因具有一定的主动阻尼而可以有效抑制系统的暂态信号，但是易受外界干扰而导致控制器发散，故系统稳定性不佳，此外，受限于参考信号的估计精度与硬件系统的时延，反馈 ANC 系统能实现的降噪频段相对较窄且降噪时易出现中高频段噪声不降反增的“水床效应”，对宽带噪声的处理能力较差。

误差麦克风同时采集原始噪声和系统输出的反相声波
控制器根据混合信号计算出反相输出
次级声源播放反相声波，形成负反馈回路

混合式 ANC

结合前馈与反馈 ANC 系统的有点，同时使用参考麦克风和误差麦克风。前馈 ANC 提供主要降噪效果，反馈 ANC 修正系统误差和环境变化，提升鲁棒性。

宽带与窄带 ANC 系统

宽带 ANC 系统通常采用参考传声器来构造参考噪声信号，进而实现消除较宽频段内的有效噪声，基本结构类似于前馈 ANC 系统。当参考传声器距离次级声源较近时，声反馈现象较为严重，因此在其中引入声反馈消除技术是十分必要

针对由发动机、压缩机、风扇等旋转机械部件运动而产生的周期性窄频带噪声，通过转速传感器、振动加速度传声器等非声学传感器的使用可辅助构造出较为精准的窄带参考信号，这类以非声学传感器来获取参考信号的 ANC 系统一般称之为窄带 ANC 系统。由于非声学传感器的使用，窄带 ANC 系统不存在声反馈问题，同时也避免了参考传声器所存在的老化和非线性问题，此外，窄带阶次噪声通常具有较强的周期特性，由转速传感器等所构造的参考信号可有效满足系统对目标噪声的连续追踪需求，因此，窄带 ANC 系统通常具有较优良的窄带阶次分量降噪特性。

在实际生活中，大多数目标噪声均具有宽窄带混合的特性，采用由宽带 ANC 子系统和窄带 ANC 子系统所构成的宽窄带混合 ANC 系统往往能实现更为理想的控制效果。但是整体的结构会比较复杂，同时采用了声学传声器与非声学传感器来构造参考噪声信号，次级声源发出的抵消信号也是由宽、窄带 ANC 部分的输出信号求和得出

上述 ANC 类型，宽带与窄带的区别只是在于针对的噪声频率范围，宽带处理宽频率范围内的噪声，需要比较复杂的自适应算法，而窄带 ANC 针对特定频率的噪声，算法较为简单，计算量小。

LMS 算法

LMS 算法是自适应滤波领域中基础且广泛应用的算法，核心是通过最陡下降法随机近似，实时调整滤波器权重，使输出信号域期望信号的均方误差最小。LMS 的计算量小，实现简单，实时性强。

LMS 基于线性自适应滤波器，核心是用输入信号的线性组合逼近期望信号，LMS 中有如下定义

输入向量 $X(n)=[x(n),x(n-1),\cdots,x(n-L+1)]$ 长度为 $L$ 的滑动窗口，其中 $L$ 为滤波器的阶数
权重向量 $w(n)=[w_0(n),w_1(n),\cdots,w_{L-1}(n)]$ 是待优化的滤波系数
滤波器输出 $y(n)=w^T(n)X(n)=\sum_i^{L-1}w_i(n)x(n-i)$
期望信号 $d(n)$ ，是目标信号，例如在降噪中是需要抵消的噪声信号，这个信号是通过麦克风采集得到的
误差信号 $e(n)=d(n)+y(n)$ ，即输出与期望的差值，即 ANC 的残余噪声，这个信号应当是误差接收麦克风所接收到的信号。

LMS 的优化目标是最小化均方误差，即

J(n)=E(e^2(n))=E\Big[(d(n)+w^T(n)X(n))^2\Big]

利用最陡下降法通过沿误差梯度的反方向更新权值，使 $J(n)$ 逐步最小化

w(n+1)=w(n)-\mu\nabla J(n)

其中

$\mu$ 是步长因子，用于控制收敛速度与稳定性。它的大小很大程度上决定了算法的收敛性和稳态性能，当 $\mu$ 越大，算法收敛越快，但稳态误差也越大，当 $\mu$ 越小，则算法收敛慢，但稳态误差也越小，为了保证算法稳态收敛，应道限制 $0<\mu<\frac{2}{\lambda_{\max}}$ ，其中 $\lambda_{\max}$ 是输入信号的自相关矩阵 $R=E[x(n)x^T(n)]$ 的最大特征值。或者咋工程中也有简化 $0<\mu<\frac{1}{tr(R)}$ ，即自相关矩阵的迹。
$\nabla J(n)$ 是均方误差对权值的梯度，展开为 $\nabla J(n)=-2E[e(n)x(n)]$

由于最陡下降法需要统计期望，无法实时计算，所以在 LMS 用瞬时梯度代替期望梯度，实现实时权值更新

\hat{\nabla J(n)}=2e(n)x(n)

带入最陡下降公式，得到 LMS 权值更新的核心公式

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

通常将 $2\mu$ 合并为新的步长 $\mu^\prime$ ，即 $w(n+1)=w(n)-\mu^\prime e(n)x(n)$ 。用当前时刻的瞬时误差和输入信号，实时修正权值，让滤波器逐步逼近最优解。

LMS 的算法步骤如下

初始化：设定滤波器阶数 $L$ 和步长因子 $\mu$ ，将权值向量初始化 $w(0)=0$ ，输入信号缓存 $X(t)$ 初始化全为 0
在每个时刻 $t$ $t$
- 更新输入缓存：利用滑动窗口，将新样本加入，丢弃最旧的样本 $X(t)=[x(t),x(t-1),\cdots,x(t-L+1)]$
- 计算滤波器输出 $y(t)=w^T(t)X(t)$
- 得到误差信号 $e(t)=d(t)+y(t)$
- 更新权值向量 $w(t+1)=w(t)-\mu e(t)X(t)$
迭代直到预设次数或者信号处理结束

NLMS 算法

NLMS 算法是 LMS 算法的归一化版本，在 LMS 算法中，步长是一个稳定的常数，在 NLMS 算法中，步长是一个随时间变化的量，定义为

\mu(n)=\frac{\alpha}{\Vert X(n)\Vert^2+c}

其中定义 $0<\alpha<2$ 以满足收敛条件，其中 $c$ 是一个很小的常数，防止除零异常，其余的步骤与 LMS 一致

VSS-LMS 算法

VSS-LMS 算法即变步长最小均方算法，是解决了传统 LMS 算法中收敛速度和稳态误差之间难以调和的根本矛盾。VSS-LMS 算法使用一个能动态计算的时变参数 $\mu(n)$ ，取代了标准 LMS 算法中固定步长因子 $\mu$ ，从而实现性能的平衡。VSS-LMS 算法遵循一个原则：误差越大，生成的步长 $\mu(n)$ 也越大，常见的函数形式如下

Kwong-Johnston VSS-LMS
- 更新公式为 $\mu(n)=\alpha \mu(n-1)+\gamma e^2(n)$
- 其中 $0<\alpha<1$ 为平滑因子，其中 $\gamma>0$ 为步长增益，控制步长对误差变化的敏感度
- 通常取值 $\alpha\in[0.95,0.99]$ ，步长增益 $\gamma\in[1e^{-6},1e^{-3}]$
Aboulnasr-Mayyas VSS-LMS
- 更新公式为 $p(n)=\beta p(n-1)+(1-\beta)e(n)e(n-1)$ 和 $\mu(n)=\alpha \mu(n-1)+\gamma p^2(n-1)$
- 通常取值 $\alpha\in[0.95,0.99]$ ，步长增益 $\gamma\in[1e^{-6},1e^{-3}]$ ，其中 $0<\beta<1$ 为误差自相关平滑因子，典型值为 $0.9\sim0.95$
Sigmoid 型 VSS-LMS
- 更新公式为 $\mu(n)=\mu_{min}+\frac{\mu_{max}-\mu_{min}}{1+\exp(-\beta(\vert e(n)\vert-\delta))}$
- 其中 $\beta$ 为 Sigmoid 函数斜率，越大则步长过渡越陡峭，其中 $\delta$ 为误差阈值，当误差 $\vert e(n)\vert>\delta$ 时步长接近 $\mu_{max}$ ，否则接近 $\mu_{min}$
- 通常取值 $\beta\in[1,20]$ ，误差阈值 $\delta\in[0.1\sigma_e,\sigma_e]$ ，其中 $\sigma_e$ 为稳态误差的标准差，可以先运行 LMS 测量标准差
基于误差-输入互相关的 VSS-LMS
- 更新公式为 $\mu(n)=\mu_{min}+(\mu_{max}-\mu_{min})(1-\exp(-\beta\vert e(n)x(n)\vert))$
- 其中控制系数通常取值 $\beta\in[1e^{-3},1e^{-1}]$
误差绝对值 VSS-LMS
- 步长更新公式为 $\mu(n+1)=\alpha\mu(n)+\gamma \vert e(n)\vert$
- 通常取值 $\alpha\in[0.95,0.99]$ ，步长增益 $\gamma\in[1e^{-6},1e^{-3}]$
峰度控制 VSS-LMS
- 步长更新公式 $\mu(n)=\mu_{max}\Big(1-\exp\big(-\alpha\mathrm{Kurt}(e(n))\big)\Big)$
- 其中 $\mathrm{Kurt}(X)=E\Big[(\frac{X-\mu}{\sigma})^4\Big]$ 为峰度函数，衡量的是随机变量分布相对于标准正态分布的偏离程度
- 其中峰度控制系数 $\alpha\in[0.1,1.0]$
归一化 VSS-LMS
- 主要在于权重更新公式 $w(n+1)=w(n)+\frac{\mu(n)e(n)\Vert x(n)\Vert}{\epsilon+\Vert x\Vert_2}$ ，其中 $\epsilon$ 为防止除零的小正数
- 更新步长的公式可以采用上述的方法来更新步长

无论采用哪种更新规律，在实际使用中都会加入饱和限制，即

\mu(n)=\min(\mu_{max},\max(\mu_{min},\mu(n)))

FxLMS

FxLMS 是 LMS 算法针对 ANC 系统中次级通道延迟与幅度衰减的专用改进版本，解决了基础的 LMS 因忽略从扬声器到误差麦克风声的传播路径导致的收敛慢、降噪效果差甚至是发散的问题，是工业上的标准算法。FxLMS 的核心逻辑就是先给参考信号预滤波，再用滤波后的信号更新权值，让自适应滤波器的输出能精准抵消经过次级通道传播后的噪声

基础 LMS 的误差信号为 $e(n)=d(n)+y(n)$ ，在 ANC 中，扬声器输出的反相声波需经过次级通道才能到达误差麦克风，所以在 ANC 中的实际误差信号为 $e(n)=d(n)+s^T(n)Y(n)$ ，其中 $Y(n)=[y(n),y(n-1),\cdots]$ 。基础的 LMS 未考虑 $s(n)$ 的卷积效应，用为滤波的参考信号来更新权值，导致权值更新方向错误，收敛极慢甚至无法降噪。注意这里的 $s(n)$ 就是声波在经过扬声器—功放—声波传播—误差麦克风—ADC这条次级通道的时域脉冲响应，是未知的，需要用 $\hat{s}(n)$ 来近似估计它。需要注意的是，算法里所用到的 $d(n)$ 实际上是麦克风采集到的噪声，它与源噪声肯定会有一些区别的，它就是经过了 $s(n)$ 变换之后的源噪声。FxLMS 中的定义如下

参考麦克风信号 $x(n)$ ，构造向量 $X(n)=[x(n),x(n-1),\cdots,x(n-L+1)]$ ，阶数为 $L$
扬声器输出信号 $y(n)=w^T(n)X(n)$ ，构造向量 $Y(n)=[y(n),y(n-1),\cdots,y(n-L-1)]$
误差麦克风接收误差为 $e(n)=s^T(n)D(n)+s^T(n)Y(n)$ ，其中 $D(n)$ 为噪声源向量，表示从噪声源传播到误差麦克风的过程
利用当前次级路径模型生成滤波参考信号 $x^\prime(n)=\hat{s}^T(n)X(n)$ ，构造向量 $X^\prime(n)=[x^\prime(n),x^\prime(n-1),\cdots,x^\prime(n-L-1)]$
更新权重 $w(n+1)=w(n)-\mu e(n)X^\prime(n)$

为了保证算法稳态收敛，应道限制 $0<\mu<\frac{2}{\lambda_{\max}}$ ，其中 $\lambda_{\max}$ 是输入信号的自相关矩阵 $R=E[x^\prime(n)(x^{\prime})^T(n)]$ 的最大特征值。或者咋工程中也有简化 $0<\mu<\frac{1}{tr(R)}$ ，即自相关矩阵的迹。

另外 FxLMS 的性能高度依赖次级通道估计 $\hat{s}(n)$ 的精度，次级通道估计 $\hat{s}(n)$ 的长度为 $M$ ，得到 $\hat{s}(n)$ 的方法如下

离线辨识：系统启动前，用白噪声激励扬声器，采集误差麦克风信号，通过系统辨识得到 $\hat{s}(n)$
在线辨识：实时更新 $\hat{s}(n)$ ，适合动态环境
在工程中使用短长度 FIR 滤波器近似 $\hat{s}(n)$ ，平衡精度与计算量

次级通道估计在线辨识

离线辨识有个致命缺陷，即环境变化时，会导致真实的 $s(n)$ 漂移，离线 $\hat{s}(n)$ 失效，所以动态场景下需要在线辨识，实时更新 $\hat{s}(n)$ ，需要注意的是，初始化的 $\hat{s}(n)$ 不能直接初始化为 0，建议是使用离线辨识先更新一次，设定初值，或者也可以设定为一个较小的初始值。

无辅助噪声的在线辨识：由于部分场景不允许注入辅助噪声，可以利用主滤波器的输出 $y(n)$ $y (n)$ 作为辨识激励
- 扬声器驱动信号为 $u(n)=y(n)=w^T(n)X(n)$
- 误差麦克风信号为 $e(n)=d(n)+s^TY(n)$ ，其中 $Y(n)=[y(n),y(n-1),\cdots,y(n-L+1)]$
- 定义辨识误差 $f(n)=e(n)-\hat{s}^TY(n)$ ，当 $\hat{s}\rightarrow s$ ，则 $f(n)\rightarrow d(n)$
- 构造辅助变量向量 $Z(n)=[y(n-\Delta),y(n-\Delta-1),\cdots,y(n-\Delta-L+1)]$
- 次级路径模型更新 $\hat{s}(n+1)=\hat{s}(n)+\mu_s f(n)Z(n)$
- 利用当前次级路径模型生成滤波参考信号 $x^\prime(n)=\hat{s}^T(n)X(n)$ ，构造向量 $X^\prime(n)=[x^\prime(n),x^\prime(n-1),\cdots,x^\prime(n-L-1)]$
- 控制滤波器更新 $w(n+1)=w(n)-\mu_we(n)X^\prime(n)$
有辅助噪声的在线辨识：扬声器的输出为主滤波器的输出 $y(n)=w^T(n)X(n)$ $y (n) = w^{T} (n) X (n)$ 与辅助噪声 $v(n)$ $v (n)$ 的叠加
- 误差麦克风信号为 $e(n)=d(n)+s^TY(n)+s^TV(n)$ ，其中 $Y(n)=[y(n),y(n-1),\cdots,y(n-L+1)],V(n)=[v(n),v(n-1),\cdots,v(n-L+1)]$
- 定义辨识误差 $f(n)=e(n)-\hat{s}^TY(n)-\hat{s}^TV(n)$ ，当 $\hat{s}\rightarrow s$ ，则 $f(n)\rightarrow d(n)$
- 以辅助噪声向量 $V(n)$ 作为参考输入，以 $f(n)$ 为期望响应，更新 $\hat{s}(n+1)=\hat{s}(n)+\mu_s f(n)V(n)$
- 由当前次级路径模型生成滤波参考信号 $x^\prime(n)=\hat{s}^T(n)X(n)$ ，形成向量 $X^\prime(n)=[x^\prime(n),x^\prime(n-1),\cdots,x^\prime(n-L+1)]$
- 控制滤波器系数更新 $w(n+1)=w(n)-\mu_we(n)X^\prime(n)$

次级通道估计离线辨识

离线辨识是在系统启动前完成的辨识，适合环境固定的场景，它的精度高、计算简单。它的核心流程如下

准备激励信号：用已知的、宽带、能量均匀的信号激励扬声器，覆盖 ANC 的目标频段
同时同步采集扬声器输入激励 $x(n)$ 和麦克风输出 $y(n)$
用算法从 $\{x(n),y(n)\}$ 中反推出次级通道脉冲响应 $\hat{s}(n)$

其中的激励信号直接决定辨识精度，在工程中常用的有 3 种

白噪声：频谱平坦，覆盖全频段，实现简单
线性扫频信号：能量集中，抗干扰性强，频率分辨率高
伪随机序列：确定性，抗干扰，可重复，频谱接近白噪声

在工程简化中，优先使用白噪声，代码容易实现且满足多数场合需求。

离线辨识有两种方法：即时域辨识和频域辨识

时域辨识：
- 麦克风输出激励信号 $u(n)$
- 误差麦克风接收噪声误差 $e(n)=s^TV(n)$ ，其中 $V(n)=[v(n),v(n-1),\cdots]$
- 计算估计麦克风噪声误差 $y(n)=\hat{s}^TV(n)$
- 计算误差 $e_{off}=e(n)-y(n)$
- 更新 $\hat{s}(n+1)=\hat{s}(n)+\mu e(n)X(n)$ ，为了保证离线辨识时能稳定不发散，可以使用归一化更新的方式，即 $\hat{s}(n+1)=\hat{s}(n)+\frac{\mu}{\Vert X(n)\Vert^2+\delta}e(n)X(n)$
频域辨识：
- 麦克风输出激励信号 $u(n)$ ，误差麦克风采集信号为 $e(n)=s(n)*u(n)+v(n)$ ，其中 $v(n)$ 为测量噪声，与 $u(n)$ 相互独立
- 对上述采集信号做傅里叶变换，得到频域公式 $Y(f)=S(f)U(f)+V(f)$ ，通过有限长的记录来估计 $S(f)$
- 假设输入端 $U(f)$ 不受噪声污染，输出端存在独立噪声 $V(f)$ ，此时对 $S(f)$ 的无偏估计为 $\hat{S}(f)=\frac{P_{UY}(f)}{P_{UU}(f)}=\frac{E[U(f)U^*(f)]}{E[U(f)Y^*(f)]}$ ，即激励与响应的互功率谱密度除以激励的自功率谱密度
- 利用 Welch 法进行谱估计：将长数据分为 $M$ $M$ 个长度为 $N$ $N$ 的重叠段，每段加窗 $w(n)$ $w (n)$ ，计算修正周期图，再对多段平均
  - 第 $i$ 段数据的傅里叶变换 $U_i(k)=\sum_{n=0}^{N-1}u_i(n)w(n)e^{-j\frac{2\pi kn}{N}}$ 和 $Y_i(k)=\sum_{n=0}^{N-1}y_i(n)w(n)e^{-j\frac{2\pi kn}{N}}$
  - 计算平均自谱和互谱 $\hat{P}_{UU}(k)=\frac{1}{M}\sum_{i=1}^M\vert U_i(k)\vert^2$ 和 $\hat{P}_{UY}(k)=\frac{1}{M}\sum_{i=1}^MU_i^*(k)Y_i(k)$
  - 求比值得到频率响应估计 $\hat{S}(k)=\frac{\sum_{i=1}^MU_i^*(k)Y_i(k)}{\sum_{i=1}^M\vert U_i(k)\vert^2}$
- 得到的频率响应估计进行逆傅里叶变换 $\hat{s}_{raw}(n)=IFFT(\hat{S}(k))$ ，取实部 $\hat{s}(n)=Re(\hat{s}_{raw}(n))$
- 此时得到的 $\hat{s}(n)$ 长度等于 $N$ ，但有效部分通常远短于此，观察幅值衰减，选取显著非零的长度 $L$ 进行矩形或平顶窗截断，得到 $\hat{s}_f=\hat{s}(n)w_t(n)$

FxNLMS

FxNLMS 算法就是再 FxLMS 算法上添加步长归一化函数即可，即

\mu(n)=\frac{\alpha}{\Vert X(n)\Vert^2+c}

与 NLMS 算法不同的是，这里需要用到的 $X(n)$ 应该是由当前次级路径模型生成的滤波参考信号 $X^\prime(n)$ ，即

\mu(n)=\frac{\alpha}{\Vert X^\prime(n)\Vert^2+c}

其中定义 $0<\alpha<2$ 以满足收敛条件，其中 $c$ 是一个很小的常数，防止除零异常，其余的步骤与 FxLMS 一致

VFxLMS

VFxLMS 即变步长滤波的 FxLMS 算法，是 FxLMS 算法的核心改进，解决固定步长难以同时兼顾收敛速度与稳态误差的矛盾，在主动噪声控制中广泛应用。它的核心优势就是：初始阶段大步长快速收敛，接近最优解时使用小步长保证低残余噪声，适配动态环境。

在 FxLMS 中，固定步长的选择需要注意

大步长：收敛快，但稳态过量均方误差大，甚至可能导致系统不稳定
小步长：稳态 EMSE 小，但收敛慢，无法跟踪动态噪声变化
实际使用中，步长需与信号统计特性实时匹配，固定步长无法自适应调整

VFxLMS 在 FxLMS 的基础上，将固定的步长替换为时变步长 $\mu(n)$ ，根据系统的状态实时调整步长更新规则。VFxLMS 集成 FxLMS 的次级通道补偿机制，仅修改权重更新的步长项

对于变步长的 $\mu(n)$ 的设计规则如下

基于误差信号的变步长：误差大时增大步长，误差小时减小步长，自适应误差调节步长
- 误差平方型 $\mu(n)=\mu_0\vert e(n)\vert^p$ ，其中 $p=1,2$ ，用于控制步长变化灵敏度，当 $p=2$ 时对误差更敏感
- 改进双曲线正切型 $\mu(n)=\mu_{\max} \tanh(\alpha\vert e(n)\vert^2+\beta)$ 。其中参数 $\alpha$ 控制步长增长速率， $\beta$ 调整阈值， $\mu_{\max}$ 是最大允许步长
归一化变步长：步长与滤波后的参考信号功率成反比，抵消信号功率波动对收敛的影响
- $\mu(n)=\frac{\mu_0}{\Vert X^\prime(n)\Vert^2+\delta}$
- 其中 $\delta$ 为极小的一个正数，防止分母为 0
误差与信号联合归一化：结合误差能量与参考信号功率，进一步提升鲁棒性
- $\mu(n)=\frac{\mu_0\vert e(n)\vert}{\Vert X^\prime(n)\Vert^2+E[e(^2(n)]+\delta}$
- 其中 $E[e^2(n)]$ 是误差信号的指数移动平均功率，用 $E[e^2(n)]=\lambda E[e^2(n-1)]+(1-\lambda)e^2(n)$ 计算，其中的 $\lambda$ 为遗忘因子，一般选择 $0.9\sim0.99$
基于特征值约束的变步长：步长上限由滤波参考信号的自相关矩阵的最大特征值 $\lambda_{\max}$ $λ_{m a x}$ 决定，保证稳定性
- $0<\mu(n)<\frac{2}{\lambda_{\max}}$

需要注意的是，VFxLMS 的性能依赖于 $\hat{s}(n)$ 的精度，需要先通过离线辨识获得高精度次级通道估计，动态场景需要在线更新 $\hat{s}(n)$ ，这个与 FxLMS 中的一致。

RLS

RLS 算法通过递归更新滤波器的系数 $w(n)$ ，它的代价函数为所有历史误差的加权平方和，如下

J(n)=\sum_{i=0}^n\lambda^{n-i}\Vert e(i)\Vert_2

其中的 $0<\lambda<1$ ，表示历史数据的权重随着时间推移越来越小，但是通常取 $\lambda\in[0.95,0.999]$ 之间，跟踪快变系统用 $0.95\sim0.98$ ，对于稳定系统使用 $0.995$ 以上。其中的误差 $e(n)$ 定义如下

e(n)=d(n)+w^TX(n)

令代价函数对 $w$ 的梯度为 0，可得到确定性正规方程

\frac{\partial J(n)}{\partial w(n)}=2\sum_{i=0}^n\lambda^{n-i}X(i)e(i)=2\sum_{i=0}^n\lambda^{n-i}X(i)(d(i)-w^TX(i))=0

由于 $w^TX(i)$ 是一个标量，所以在运算中可以进行换位，将 $X(i)w^TX(i)$ 换为 $X(i)X^T(i)w$ ，变换前后值不变。所以可以得到

\sum_{i=0}^n\lambda^{n-i}X(i)d(i)=\sum_{i=0}^n\lambda^{n-i}X(i)X^T(i)w\\\Downarrow\\R(n)w=r(n)\\\Downarrow\\w=R^{-1}(n)r(n)

但是直接对 $R(n)$ 求逆的话，运算量会比较大，所以可以将 $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)

对上述第一个公式左右两边求逆矩阵，利用 Woodbury 公式 $(A+BCD)^{-1}=A^{-1}-A^{-1}B(C^{-1}+DA^{-1}B)^{-1}DA^{-1}$ ，带入 $A=\lambda R(n-1),B=X(n),C=1,D=X^T(n)$ ，且令 $P(n)=R^{-1}(n)$ 可以得到

P(n)=\lambda^{-1}P(n-1)-\frac{\lambda^{-2}P(n-1)X(n)X^T(n)P(n-1)}{1+\lambda^{-1}X^TP(n-1)X(n)}

得到 $P(n)$ 之后，就需要找 $w(n)$ 满足 $R(n)w(n)=r(n)$ ，也就是权重需要满足这个方程。对于上一时刻的解 $w(n-1)$ 来说，满足 $R(n-1)w(n-1)=r(n-1)$ ，而 $R$ 和 $r$ 都会更新，所以将 $w(n-1)$ 带入方程的左边就会得到

R(n)w(n-1)=\lambda R(n-1)w(n-1)+X(n)X^T(n)w(n-1)

由于 $R(n-1)w(n-1)=r(n-1)$ ，所以带入之后可以得到

R(n)w(n-1)=\lambda r(n-1)+X(n)X^T(n)w(n-1)

对比之前的 $r(n)=\lambda r(n)+X(n)d(n)$ 可以发现它们之间差了一项 $X(n)(d(n)-X^T(n)w(n-1))=X(n)e(n)$ ，即

R(n)w(n-1)=r(n)-X(n)e(n)

所以此时就希望能找到一个 $\Delta w$ 使得 $w(n)=w(n-1)+\Delta w$ 能满足 $R(n)w(n)=r(n)$ ，带入之后，利用上述的公式，两边消去 $r(n)$ 即可得到

R(n)(w(n-1)+\Delta w)=r(n)\\\Downarrow\\R(n)\Delta w=X(n)e(n)

解出修正量

\Delta w=R^{-1}(n)X(n)e(n)=P(n)X(n)e(n)

此时引入一个增益向量来简化计算，带入上述的 $P(n)$ 的值，求解可以得到

K(n)=P(n)X(n)=\frac{P(n-1)X(n)}{\lambda+X^T(n)P(n-1)X(n)}

带入得到 $P(n)$ 的更新公式中，简化计算，得到

P(n)=\lambda^{-1}\Big[P(n-1)-K(n)X^T(n)P(n-1)\Big]

更新权向量

w(n)=w(n-1)-K(n)e(n)

初始时刻设置 $w(0)=0$ 和 $P(0)=\delta^{-1} I$ ，其中的 $\delta$ 的取值表示初始时刻对输入相关矩阵的估计，如果信号功率大则 $\delta$ 就取稍大一点。

FRLS

FRLS 算法即快速递归最小二乘法，是标准 RLS 算法的一种高效实现。在标准 RLS 算法中，权向量更新公式如下

w(n)=w(n-1)-K(n)e(n)

其中的增益向量 $K(n)=P(n)X(n)$ 需要较大计算量，输入信号 $X(n)$ 具有位移结构，利用这一结构，可以通过引入低阶的前向预测和后向预测来递推增益向量，从而避免直接更新更新 $L\times L$ 的逆矩阵 $P(n)$ 。定义先验 Kalman 增益向量

K^\prime=P(n-1)X(n)

这与标准 RLS 算法中的 Kalman 增益不同。定义转换因子

\theta(n)=\lambda+X^TK^\prime(n)

是一个标量，用于归一化 Kalman 增益。定义前向预测误差和后向预测误差

e_f(n)=x(n)-a^T(n-1)X(n-1)\\e_b(n)=x(n-L)-b^T(n-1)X(n)

其中 $a(n)$ 是前向预测器系数向量，用于用过去的 $L$ 个样本预测当前样本 $x(n)$ ，而 $b(n)$ 为后向预测器系数向量，用于用当前的 $L$ 个样本预测过去的样本 $x(n-L)$ 。另外定义前向预测误差能量和后向预测误差能量

E_f(n)=\lambda E_f(n-1)+\frac{e_f^2(n)}{\theta(n-1)}\\E_b(n)=\lambda E_b(n-1)+\frac{e_b^2(n)}{\theta(n)}

FRLS 算法的核心是阶数更新和时间更新的结合，定义扩展 Kalman 增益向量 $K_{L+1}(n)$ ，它是 $L+1$ 阶的 Kalman 增益向量，可以通过 $L$ 阶的 Kalman 增益向量和前向预测器递推得到

K_{L+1}(n)=\begin{bmatrix}0\\K_L(n-1)\end{bmatrix}+\begin{bmatrix}1\\-a(n-1)\end{bmatrix}\frac{e_f(n)}{E_f(n-1)}\\K_{L+1}(n)=\begin{bmatrix}K_L(n)\\j(n)\end{bmatrix}

将 $K_L(n)$ 分块就是需要的 $L$ 阶先验 Kalman 增益向量，上述的 $j(n)$ 为标量。前向和后向预测器系数可以通过卡尔曼增益和预测误差递归更新

前向预测器更新 $a(n)=a(n-1)+K_L(n-1)\frac{e_f(n)}{\theta(n-1)}$
后向预测器更新 $b(n)=b(n-1)+K_L(n)\frac{e_b(n)}{\theta(n)}$

上述的转换因子更新公式为 $\theta(n)=\theta(n-1)+\frac{e_f^2(n)}{E_f(n-1)}-e_b(n)j(n)$ 。整理上述公式可以得到 FRLS 算法的完整流程

初始化
$w(n)=0\\a(0)=0\\b(0)=0\\K(0)=0\\\theta(0)=\lambda\\E_f(0)=\delta\\E_b(0)=\delta\\\lambda\in(0,1]$
计算前向预测误差 $e_f(n)=x(n)-a^T(n-1)X(n-1)$
更新前向预测误差能量 $E_f(n)=\lambda E_f(n-1)+\frac{e_f^2(n)}{\theta(n-1)}$
前向预测器更新 $a(n)=a(n-1)+K_L(n-1)\frac{e_f(n)}{\theta(n-1)}$
计算后向预测误差 $e_b(n)=x(n-L)-b^T(n-1)X(n)$
更新扩展 Kalman 增益向量 $\begin{bmatrix}K_L(n)\\j(n)\end{bmatrix}=\begin{bmatrix}0\\K_L(n-1)\end{bmatrix}+\begin{bmatrix}1\\-a(n-1)\end{bmatrix}\frac{e_f(n)}{E_f(n-1)}$
更新转换因子 $\theta(n)=\theta(n-1)+\frac{e_f^2(n)}{E_f(n-1)}-e_b(n)j(n)$
更新后向预测误差能量 $E_b(n)=\lambda E_b(n-1)+\frac{e_b^2(n)}{\theta(n)}$
后向预测器更新 $b(n)=b(n-1)+K_L(n)\frac{e_b(n)}{\theta(n)}$
得到误差信号 $e(n)=d(n)-w^T(n-1)X(n)$
更新滤波器参数 $w(n)=w(n-1)-K_L(n)\frac{e(n)}{\theta(n)}$

FxRLS

FxRLS 算法是递归最小二乘 RLS 算法上存在次级路径系统中的扩展，继承了 RLS 算法快速收敛、稳态误差小的优点，同时引入滤波操作解决了标准 RLS 在次级路径存在时的稳定性问题和收敛性能下降问题

FxRLS 相较于 RLS 的改进点在于：在进行权值更新前，先将原始参考信号 $x(n)$ 通过估计的次级路径 $\hat{s}(n)$ 进行滤波，得到滤波后的参考信号 $x^\prime(n)$ ，然后用 $x^\prime(n)$ 代替 $x(n)$ 进行 RLS 更新。剩下的部分与 RLS 算法一致，将响应的定义代入即可，更新公式如下

初始化 $w(0)=0$ 和 $P(0)=\delta^{-1}I$ ，其中的 $\delta$ 为很小的正数
滤波后的参考信号 $x^\prime(n)=\hat{s}(n)X(n)$ ，构造向量 $X^\prime(n)$
麦克风输出信号 $y(n)=w^TX(n)$ ，构造向量 $Y(n)$
获得误差信号 $e(n)=s^T(n)D(n)+s^T(n)Y(n)$
计算增益向量 $K(n)=P(n)X^\prime(n)=\frac{P(n-1)X^\prime(n)}{\lambda+(X^\prime)^T(n)P(n-1)X^\prime(n)}$
计算协方差矩阵 $P(n)=\lambda^{-1}(P(n-1)-K(n)(X^\prime)^T(n)P(n-1))$
更新权值向量 $w(n)=w(n-1)-K(n)e(n)$

次级通道模型辨识的方法与 FxLMS 中的辨识方法一致

FxKalman

FxKalman 算法建立在 Kalman 滤波的最优状态估计框架之上，核心是将自适应滤波器的系数向量 $w$ 当做系统的状态，用参考信号的滤波 $X^\prime(n)$ 构建观测方程，再通过 Kalman 实时最优估计这个状态。相比于 FxLMS 这类随机梯度算法，FxKalman 的收敛速度不依赖参考信号的自相关矩阵特征值散布，对宽带噪声有极快的收敛能力和极低的稳态失调。更新公式如下

参考麦克风采集的参考信号 $x(n)$ ，构造向量 $X(n)$
计算扬声器输出信号 $y(n)=w^T(n)X(n)$ ，构造向量 $Y(n)$
误差麦克风采集的残余信号 $e(n)=s^T(n)D(n)+s^T(n)Y(n)$
计算次级路径估计滤波信号 $x^\prime(n)=\hat{s}^T(n)X(n)$ ，构造向量 $X^\prime(n)$
更新预测 $w(n\vert n-1)=w(n-1\vert n-1)$ 和 $P(n\vert n-1)=P(n-1\vert n-1)+Q$
计算 Kalman 增益 $K(n)=\frac{P(n\vert n-1)X^\prime(n)}{R+(X^\prime)^T(n)P(n\vert n-1)X^\prime(n)}$
测量更新 $w(n\vert n)=w(n\vert n-1)+K(n)e(n)$ 和 $P(n\vert n)=\Big[I-K(n)(X^\prime)^T(n)\Big]P(n\vert n-1)$

其中 $Q$ 为过程噪声协方差矩阵，越大则表示算法对系数变化的跟踪能力越强，但稳态失调也越大；其中 $R$ 为观测噪声方差，决定了算法对测量误差的信任程度，它越小则卡尔曼增益越大，收敛越快，但可能对测量噪声更敏感。参数初始化 $P(0\vert0)=\delta^{-1}I$ ，其中的 $\delta$ 为小正数。

APA

APA（Affine Projection Algorithm，仿射投影算法）是一种经典的自适应滤波算法，是为了解决 LMS 和 NLMS 算法在处理相关输入信号时收敛速度显著下降的问题的，它可以看作时 NLMS 算法的自然推广。NLMS 算法每次迭代时只使用 1 个当前输入向量进行权值更新，而 APA 算法每次迭代同时使用 $K$ 个最近输入的向量进行权值更新。它的收敛速度相较于 LMS 和 NLMS 快得多，计算复杂度比 RLS 低得多

从几何角度来看，K 阶仿射投影就是将当前权向量正交投影到由 K 个超平面相交形成的 K-L 维仿射子空间上。当 $K=1$ 时算法就是 NLMS 算法，在 NLMS 算法中，权重步长更新公式如下

\mu(n)=\frac{\alpha}{\Vert X(n)\Vert^2+c}

是当前权重向量在单个超平面上的正交投影，这个超平面由当前的参考信号向量定义，满足后验误差为 0 的约束。在 APA 算法中，同时施加 $K$ 个约束条件，这些约束条件共同定义了一个仿射子空间，更新后的权重向量是当前权向量在这个仿射子空间上的正交投影

在 APA 算法中，定义了输入矩阵 $M$ ，如下

M(n)=\begin{bmatrix}X(n)&X(n-1)&\cdots&X(n-K+1)\end{bmatrix}\\X(n)=\begin{bmatrix}x(n)&x(n-1)&\cdots&x(n-L+1)\end{bmatrix}^T

定义期望信号向量

D(n)=\begin{bmatrix}d(n)&\cdots&d(n-K+1)\end{bmatrix}

定义先验误差向量

e(n)=D(n)-M^T(n)w(n-1)

此时优化问题可以重写如下

\min_{w(n)}\Vert w(n)-w(n-1)\Vert_2^2\\s.t.\quad D(n)=M^T(n)w(n)

此时利用拉格朗日乘数法来求解，构造拉格朗日函数

L(w(n),\lambda)=\Big\Vert w(n)-w(n-1)\Big\Vert_2^2+\lambda^T\Big(D(n)-M^T(n)w(n)\Big)

对 $w(n)$ 求偏导，并令其为 0

\frac{\partial L(w(n),\lambda)}{\partial w(n)}=2\Big(w(n)-w(n-1)\Big)-M(n)\lambda=0\\\Downarrow\\w(n)=w(n-1)+\frac{1}{2}M(n)\lambda

对 $\lambda$ 求偏导，并令其为 0

\frac{\partial L(w(n),\lambda)}{\partial \lambda}=D(n)-M^Tw(n)=0

将 $w(n)=w(n-1)+\frac{1}{2}M(n)\lambda$ 带入上述公式中得到

D(n)-M^T\Big(w(n-1)+\frac{1}{2}M(n)\lambda\Big)=0\\\Downarrow\\\lambda=2\Big(M^T(n)M(n)\Big)^{-1}e(n)

带入 $w(n)=w(n-1)+\frac{1}{2}M(n)\lambda$ 中得到

w(n)=w(n-1)+M(n)\Big(M^T(n)M(n)\Big)^{-1}e(n)

实际中为了防止矩阵 $M^T(n)M(n)$ 奇异，通常加入一个小的正则化参数 $\delta$ ，同时，为了控制收敛速度和稳态误差，引入步长因子 $\mu$ ，即

w(n)=w(n-1)+\mu M(n)\Big(M^T(n)M(n)+\delta I\Big)^{-1}e(n)

APA 算法的迭代过程如下

参考麦克风信号 $x(n)$ ，构造输入向量 $X(n)$
更新输入矩阵 $M(n)=\begin{bmatrix}X(n),M(:,1:K-1)\end{bmatrix}$
计算先验误差向量 $e(n)=D(n)-M^T(n)w(n-1)$ ，其中的 $D(n)=\begin{bmatrix}d(n)&\cdots&d(n-K+1)\end{bmatrix}$
计算矩阵逆 $R(n)=M^T(n)M(n)+\delta I$
更新权重向量 $w(n)=w(n-1)+\mu M(n)R^{-1}(n)e(n)$
计算当前滤波器输出 $y(n)=w^T(n)X(n)$

由于整体需要求逆计算，所以在运行中的计算量较大

VSS-APA

对于固定步长 APA 算法，存在一个固有的矛盾，即大步长时收敛速度块，但稳态误差大，甚至可能导致算法发散，而小步长稳态误差小，但收敛速度慢，跟踪能力差

为了推导最优变步长，定义权误差向量

\tilde{w}(n)=w_0-\hat{w}(n)

其中 $w_0$ 为滤波器的最优权向量，将权更新公式带入上述公式

\tilde{w}(n+1)=\tilde{w}(n)-\mu M(n)R^{-1}(n)e(n)

对两边取平方欧几里得范数并求期望

E(\Vert\tilde{w}(n+1)\Vert^2)=E\Big(\Vert\tilde{w}(n)\Vert^2\Big)-2\mu(n)E\Big(e^T(n)R^{-1}(n)M^T(n)\tilde{w}(n)\Big)+\mu^2(n)E\Big(\Vert M(n)R^{-1}e(n)\Vert^2\Big)

为了使权误差向量的范数最小化，对 $\mu(n)$ 求导并令其导数为 0，得到最优步长

\mu(n)=\frac{E\Big(e^T(n)R^{-1}(n)M^T(n)\tilde{w}(n)\Big)}{E\Big(\Vert M(n)R^{-1}e(n)\Vert^2\Big)}

但是由于 $\hat{w}(n)$ 未知，实际应用中需要对这个最优步长进行近似。典型的变步长算法如下

基于误差范数的变步长
- 步长更新公式 $\mu(n)=\mu_{max}\exp(-\alpha\Vert e(n)\Vert^2)$ 或 $\mu(n)=\frac{\mu_{max}}{1+\beta\Vert e(n)\Vert^2}$
- 其中 $\mu_{max}$ 为最大步长， $\alpha,\beta$ 为控制步长衰减速度的参数
基于误差自相关的变步长，提高对噪声的鲁棒性
- 引入 $p(n)=\gamma p(n-1)+(1-\gamma)e(n)e(n-1)$
- 步长更新 $\mu(n)=\mu_{max}(1-\exp(-\beta\vert p(n)\vert^2))$
- 其中的 $\gamma$ 为平滑因子，通常取 $0.9\sim0.99$
基于集员滤波的变步长：只有当误差超过某个预设的阈值时才更新权向量，这样可以大大降低计算复杂度
- 步长更新公式 $\mu(n)=\max(0,1-\frac{\gamma}{\Vert e(n)\Vert})$
- 其中 $\gamma$ 为误差阈值，当 $\Vert e(n)\leq\gamma$ 时 $\mu(n)=0$ ，权向量不更新。步长随着误差增大而增大
指数型变步长：直接将投影误差的范数通过指数函数映射得到步长
- 步长更新公式 $\mu(n)=\mu_{min}+(\mu_{max}-\mu_{min})\exp(-\alpha\Vert e(n)\Vert^2)$

变步长公式需要在权向量更新时进行计算

FxAPA

首先是将期望矩阵 $D(n)$ 替换为 $D_f(n)$ ，另外将普通的输入矩阵 $M(n)$ 替换为 $M_f(n)$ ，如下

D(n)=\begin{bmatrix}d(n)&\cdots&d(n-L+1)\end{bmatrix}\\d_f(n)=\hat{s}^TD(n)\\D_f(n)=\begin{bmatrix}d_f(n)&\cdots&d_f(n-K+1)\end{bmatrix}\\x_f(n)=\hat{s}^TX(n)\\X_f(n)=\begin{bmatrix}x_f(n)\\\vdots\\x_f(n-L+1)\end{bmatrix}\\M_f=\begin{bmatrix}X_f(n)&\cdots&X_f(n-L+1)\end{bmatrix}

另外误差的公式为 $e(n)=s^T\Big(D_f(n)-M^T(n)w(n-1)\Big)$

FAP

FAP（Fast Affine Projection Algorithm，快速仿射投影算法）是标准 APA 算法的高效实现版本，解决了标准 APA 算法中 $K\times K$ 大小矩阵求逆的计算瓶颈，将计算复杂度从 $O(LK^2+K^3)$ 降低到 $O(L+K)$ 。在 FAP 算法中，实现加速的算法的流程如下

残差向量的时间相关性近似

在计算误差向量 $e(n)$ 时，直接使用 APA 中的公式计算需要 $O(LK)$ 次乘法运算，复杂度较高。所以利用残差向量的时间相关性，将其复杂度降低到 $O(L+K)$ 。首先将 $e(n)$ 按第一个元素和剩余 $K-1$ 个元素进行分块，得到

e(n)=\begin{bmatrix}d(n)-X^T(n)w(n-1)\\\bar{D}(n-1)-\bar{M}(n-1)w(n-1)\end{bmatrix}\in\R^{K\times1}

其中 $\bar{M}(n-1)$ 是 $M(n-1)$ 的前 $n-1$ 列，其中的 $\bar{D}(n-1)$ 表示参考向量 $D(n-1)$ 的前 $n-1$ 个元素。对于分块的下半部分，它正好是上一时刻后验残差向量 $e_1(n-1)$ 的前 $n-1$ 个元素。后验残差定义为使用更新后的权值 $w(n)$ 计算的残差

e_1(n)=D(n)-M^T(n)\Big(w(n-1)+\mu M(n)R^{-1}(n)e(n)\Big)\\=D(n)-M^T(n)w(n-1)-\mu M^T(n)M(n)R^{-1}e(n)\\=e(n)-\mu M^T(n)M(n)R^{-1}(n)e(n)\\=\Big(I-\mu M^T(n)M(n)\big(M^T(n)M(n)+\delta I\big)^{-1}\Big)e(n)

对自相关矩阵 $M^T(n)M(n)$ 做特征值分解

M^T(n)M(n)=V(n)\Lambda(n)V^T(n)

其中 $V(n)$ 为正交特征向量矩阵， $\Lambda(n)$ 为特征值对角矩阵。之后定义模态残差向量

e^\prime(n)=V^T(n)e(n)\\e_1^\prime(n)=V^T(n)e_1(n)

将上述模态残差和特征值分解带入上述的后验残差公式中，矩阵运算转化为标量运算，得到每个模态的后验-先验残差关系

er_{1,i,n-1}=\Big(1-\frac{\mu\lambda_{i,n-1}}{\delta+\lambda_{i,n-1}}\Big)er_{i,n-1}

此时分为两种极端情况

当 $\lambda_{i,n-1}\gg\delta$ 时，此时系数 $\approx1-\mu$
当 $\lambda_{i,n-1}\ll\delta$ 时，此时系数 $\approx1$

一般来说， $\delta$ 取值很小，通常为 $10^{-6}\sim10^{-3}$ ，所以将所有模态系数近似为 $1-\mu$ ，得到

e_1(n)\approx(1-\mu)e(n)

将其带入分块展开式，得到

e(n)=\begin{bmatrix}d(n)-X^T(n)w(n-1)\\\bar{e}_1(n-1)\end{bmatrix}=\begin{bmatrix}e_n\\(1-\mu)\bar{e}(n-1)\end{bmatrix}

其中的 $e_n=d(n)-X^T(n)w(n-1)$

权值向量分解和递归

在原始 APA 算法中，权值更新公式 $w(n)=w(n-1)+\mu M(n)R^{-1}(n)e(n)$ 复杂度为 $O(LK)$ ，令 $\varepsilon(n)=R^{-1}(n)e(n)$ ，此时引入辅助权重向量 $\hat{w}(n)$ 和累积向量 $E(n)$ ，将复杂度降低

w(n)=w(0)+\mu\sum_{i=0}^{n-1}M(n-i)\epsilon(n-i)\\M(n-i)\epsilon(n-i)=\sum_{j=0}^{K-1}X(n-i-j)\epsilon_j(n-i)\\w(n)=w(0)+\mu\sum_{i=0}^{n-1}\sum_{j=0}^{K-1}X(n-i-j)\epsilon_j(n-i)

其中 $M(n-i)$ 为 $n-i$ 时刻的输入矩阵，其中 $\epsilon(n-i)$ 为 $n-i$ 时刻的归一化残差向量，其中 $X(n-i-j)$ 为 $n-i-j$ 时刻的输入向量，其中 $\epsilon_j(n-i)$ 为在 $n-i$ 时刻的归一化残差向量 $\epsilon(n-i)$ 的第 $j$ 个元素，其中 $X(n-i-j)\epsilon_j(n-i)$ 表示单个输入向量对权值的贡献。将上述权重更新矩阵分解如下

w(n )=w(0)+\mu\sum_{k=0}^{K-1}X(n-k)\sum_{j=0}^k\epsilon_j(n-k+j)+\mu\sum_{k=K}^{n-1}X(n-K)\sum_{j=0}^{K-1}\epsilon_j(n-k+j)

定义辅助权值向量 $\hat{w}(n)$ ，对应为第二部分求和，即

\hat{w}(n-1)=\hat{w}(0)+\mu\sum_{k=K}^{n-1}X(n-K)\sum_{j=0}^{K-1}\epsilon_j(n-k+j)

包含了所有早于 $n-K$ 时刻的输入对权值的贡献。定义累积归一化残差向量 $E(n)$ ，对应第一部分求和中的内层和

E(n)=\begin{bmatrix}\epsilon_0(n)\\\epsilon_1(n)+\epsilon_0(n-1)\\\cdots\\\epsilon_{K-1}(n)+\epsilon_{K-2}(n-1)+\cdots+\epsilon_0(n-(K-1))\end{bmatrix}

它是最近 $K$ 个归一化残差向量的累积和，而此时第一部分的求和可以表示为 $\mu M(n)E(n)$ ，因此上述权重跟新公式可以简化为

w(n)=\hat{w}(n-1)+\mu M(n)E(n)

此时权重向量的递归更新变为了辅助权重的递归更新 $\hat{w}(n)$ ，对比辅助权重的定义，可以得到

\hat{w}(n)=\hat{w}(n-1)+\mu X(n-(K-1))\sum_{j=0}^{K-1}\epsilon_j(n-K+1+j)\\=\hat{w}(n-1)+\mu X(n-K+1)E_{K-1}(n)

其中 $E_{K-1}(n)$ 为 $E(n)$ 的最后一个元素，上述辅助权重的更新仅需要 $L$ 次乘法。将辅助权重的递归公式带入到权重的计算公式中，可以得到

w(n)=\hat{w}(n)+\mu \bar{M}(n)\bar{E}(n)

其中 $\bar{M}(n)$ 是 $M(n)$ 的前 $K-1$ 列，其中的 $\bar{E}(n)$ 为 $E(n)$ 的前 $K-1$ 个元素。另外矩阵 $E(n)$ 也需要递归更新，从上述的定义即可知 $E(n)$ 的第 $k$ 个元素等于 $\epsilon_{k}(n)$ 加上 $E(n-1)$ 的第 $k-1$ 个元素，第 0 个元素等于 $\epsilon_0(n)$ ，得到

E(n)=\epsilon(n)+\begin{bmatrix}0\\\bar{E}(n-1)\end{bmatrix}

其中 $\bar{E}(n-1)$ 是 $E(n-1)$ 的前 $K-1$ 个元素，此次更新需要 $K$ 次加法。根据 $w(n)=\hat{w}(n-1)+\mu X(n)E(n)$ ，带入标量残差 $e_n=d(n)-X^T(n)w(n-1)$ ，另外根据辅助权重的递归公式，可以得到残差如下

e_n=d(n)-X^T(n)\Big(\hat{w}(n-1)+\mu \bar{M}(n-1)\bar{E}(n-1)\Big)\\=\hat{e}_n-\mu R_{xx}(n)\bar{E}(n-1)

其中 $\hat{e}_n=d(n)-X^T(n)\hat{w}(n-1)$ 为基于辅助权值的残差。其中 $R_{xx}(n)=\begin{bmatrix}r_{xx,n}(1)&\cdots&r_{xx,n}(K-1)\end{bmatrix}^T$ 为输入信号的自相关向量，其中的 $r_{xx,n}(k)$ 为输入信号的 $k$ 阶自相关。自相关向量 $R_{xx}(n)$ 可以通过滑动窗口递归更新

R_{xx}(n)=R_{xx}(n-1)+x(n)\alpha(n)-x(n-L)\alpha(n-L)

其中的 $\alpha(n)=\begin{bmatrix}x(n-1)&\cdots&x(n-K+1)\end{bmatrix}^T$ ，每次更新需要 $2N$ 次乘法

近似 Toeplitz 矩阵的快速求逆

另外上述公式中用到了 $\epsilon(n)=\Big(M^T(n)M(n)+\delta I\Big)^{-1}e(n)$ ，其中有求逆运算，需要 $O(K^3)$ 计算量，所以可以利用 Toeplitz 矩阵的快速求逆性质和滑动窗口 FRLS 算法来简化算法。首先定义正则化自相关矩阵 $R(n)=M^T(n)M(n)+\delta I$ ，近似为一个 $K\times K$ 对称正定 Toeplitz 矩阵，它的逆矩阵可以通过前后向线性预测器快速分解。定义 $a(n)$ 为 $R(n)$ 对应的最优前向预测器，定义 $b(n)$ 为最优后向预测器，对应预测误差能量分别为 $E_{a}(n)$ 和 $E_b(n)$ ，这些量可由 sliding-window FRLS。之后定义 $R(n)$ 的两个子矩阵，定义 $\bar{R}(n)$ 为 $R(n)$ 左上角 $(K-1)\times(K-1)$ 子矩阵，定义 $\tilde{R}(n)$ 为 $R(n)$ 右下角 $(K-1)\times(K-1)$ 子矩阵。然后给出矩阵逆的两个秩一分解式

R^{-1}(n)=\begin{bmatrix}\bar{R}^{-1}(n)&0\\0^T&0\end{bmatrix}+\frac{1}{E_{b}(n)}b(n)b^T(n)\\R^{-1}(n)=\begin{bmatrix}0&0^T\\0&\tilde{R}^{-1}(n)\end{bmatrix}+\frac{1}{E_{a}(n)}a(n)a^T(n)

可以得知 $R^{-1}(n)$ 可以由一个 $(K-1)$ 阶子矩阵逆加一个秩一修正得到。此时定义

\tilde{\epsilon}(n)=\tilde{R}^{-1}(n)\tilde{e}(n)\\\bar{\epsilon}(n)=\bar{R}^{-1}(n)\bar{e}(n)

其中 $\tilde{e}(n)$ 为 $e(n)$ 的后 $K-1$ 个元素，其中 $\bar{e}(n)$ 为 $e(n)$ 的前 $K-1$ 个元素。将上述矩阵逆的两个秩一分解式分别右乘 $e(n)$ ，得到

R^{-1}(n)e(n)=\begin{bmatrix}\bar{\epsilon}(n)\\0\end{bmatrix}+\frac{1}{E_{b}(n)}b(n)b^T(n)e(n)\\R^{-1}(n)e(n)=\begin{bmatrix}0\\\tilde{\epsilon}(n)\end{bmatrix}+\frac{1}{E_{a}(n)}a(n)a^T(n)e(n)

由于滑动窗口结构 $\tilde{R}(n)=\bar{R}(n-1)$ ，即当前右下角子矩阵等于上一时刻左上角子矩阵。又由于前面有 $\tilde{e}(n)=(1-\mu)\bar{e}(n-1)$ ，所以

\tilde{\epsilon}(n)=\tilde{R}^{-1}(n)\tilde{e}(n)=\bar{R}^{-1}(n-1)(1-\mu)\bar{e}(n-1)\\=(1-\mu)\bar{\epsilon}(n-1)

这个公式让 $\epsilon(n)$ 也能递推，不需要每次直接求逆矩阵了。

整体流程

将上述所有公式整理一下，得到如下的 FAP 算法流程

初始化
- 初始化 $E_{a}(0)=E_b(0)=\delta$
- 初始化 $a(0)=\begin{bmatrix}1&\vec{0}^T\end{bmatrix}\in\R^{K\times1}$
- 初始化 $b(0)=\begin{bmatrix}\vec{0}^T&1\end{bmatrix}\in\R^{K\times1}$
- 初始化 $\hat{w}(0)=\vec{0}\in\R^{L\times1}$
- 初始化 $e(0)=\vec{0}\in\R^{K\times1}$
- 初始化 $\epsilon(0)=\vec{0}\in\R^{K\times1}$
- 初始化 $E(0)=\vec{0}\in\R^{K\times1}$
- 初始化 $R_{xx}(0)=\vec{0}\in\R^{(K-1)\times1}$
利用
- 计算先验前向预测误差 $e_f(n)=x(n)+a^T(n-1)\tilde{X}(n-1)$ 其中 $\tilde{X}=\begin{bmatrix}x(n-1)&\cdots&x(n-K)\end{bmatrix}$
- 更新前向预测器 $a^{(1)}(n)=\begin{bmatrix}1\\\bar{a}(n-1)\end{bmatrix}+\frac{e_f(n)}{E_b(n-1)}\begin{bmatrix}0\\\bar{b}(n-1)\end{bmatrix}$ ，其中 $\bar{a}(n-1)$ 和 $\bar{b}(n-1)$ 为 $a(n-1)$ 和 $b(n-1)$ 的
- 更新前向预测误差能量 $E_a^{(1)}(n)=E_a(n-1)+\frac{e_f^2}{E_{b}(n-1)}$
- 计算先验后向预测误差 $e_b(n)=x(n-K)+b^T(n-1)\tilde{X}(n)$
- 更新后向预测其 $b^{(1)}(n)=\begin{bmatrix}0\\\bar{b}(n-1)\end{bmatrix}+\frac{e_b(n)}{E^{(1)}_a(n)}\begin{bmatrix}0\\\bar{a}(n-1)\end{bmatrix}$
- 更新后向预测误差能量 $E^{(1)}_b(n)=E_b(n-1)+\frac{e_b^2}{E^{(1)}_a(n)}$
- 剔除最旧样本 $a(n)=a^{(1)}(n)-\frac{a_K^{(1)}(n)}{E_b^{(1)}(n)}b^{(1)}(n)$ ，其中 $a_K^{(1)}(n)$ 为 $a^{(1)}(n)$ 的第 $K$ 项
- 剔除最旧样本 $E_a(n)=E^{(1)}_a(n)-\frac{(a_K^{(1)}(n))^2}{E_b^{(1)}(n)}$
- 剔除最旧样本 $b(n)=b^{(1)}(n)-\frac{b^{(1)}_1(n)}{E_a(n)}a(n)$ ，其中 $b_1(n)$ 表示 $b(n)$ 的第 1 项
- 剔除最旧样本 $E_b(n)=E^{(1)}_b(n)-\frac{(b_1^{(1)}(n))^2}{E_a(n)}$
更新输入相关向量 $R_{xx}(n)=R_{xx}(n-1)+x(n)\alpha(n)-x(n-L)\alpha(n-L)$ 其中 $\alpha(n)=\begin{bmatrix}x(n-1)&\cdots&x(n-K+1)\end{bmatrix}$
计算辅助残差 $\hat{e}_n=d(n)-X^T(n)\hat{w}(n-1)$
计算标量残差 $e_n=\hat{e}_n-\mu R^T_{xx}(n)\bar{E}(n-1)$
更新残差向量 $e(n)=\begin{bmatrix}e_n\\(1-\mu)\bar{e}(n-1)\end{bmatrix}$
计算归一化残差 $\epsilon(n)=R^{-1}(n)e(n)=\begin{bmatrix}0\\\tilde{\epsilon}(n)\end{bmatrix}+\frac{1}{E_{a}(n)}a(n)a^T(n)e(n)$
计算 $\begin{bmatrix}\bar{\epsilon}(n)\\0\end{bmatrix}=\epsilon(n)-\frac{1}{E_b(n)}b(n)b^T(n)e(n)$
更新累积向量 $E(n)=\begin{bmatrix}0\\\bar{E}(n-1)\end{bmatrix}+\epsilon(n)$
更新辅助权值 $\hat{w}(n)=\hat{w}(n-1)+\mu X(n-K+1)E_{K-1}(n)$
递推下一时刻的 $\tilde{\epsilon}(n+1)=(1-\mu)\bar{\epsilon}(n)$

精确递归

上述的方法就是近似 $M^T(n)M(n)\approx Toeplitz Matrix$ 进行计算的，后续又有人提出了精确利用输入矩阵的滑动结构做精确递推的方法

定义输入矩阵 $M(n)=\begin{bmatrix}X(n)&X(n-1)&\cdots&X(n-L+1)\end{bmatrix}$ ，可以看出当前时刻的矩阵相比于前一时刻的矩阵，新增了一个新列，且删除了最旧的列。此时定义矩阵 $G(n)$

G(n)=M^T(n)M(n)+\delta I\\G(n-1)=M^T(n-1)M(n-1)+\delta I

根据输入矩阵的滑动，可以得到 $G(n)$ 矩阵的右下角的 $(L-1)\times(L-1)$ 的子矩阵和 $G(n-1)$ 左上角的 $(L-1)\times(L-1)$ 的子矩阵相同。对矩阵 $G(n)$ 进行分块

G(n)=\begin{bmatrix}a(n)&b^T(n)\\b(n)&C(n-1)\end{bmatrix}\\a(n)=X^T(n)X(n)+\delta\\b(n)=\begin{bmatrix}X^T(n)X(n-1)\\X^T(n)X(n-2)\\\vdots\\X^T(n)X(n-L+1)\end{bmatrix}\\G(n-1)=\begin{bmatrix}C(n-1)&h^T(n)\\h(n)&\gamma(n)\end{bmatrix}

假设已有 $R(n-1)=G^{-1}(n-1)$ ，将逆矩阵分块

R(n-1)=\begin{bmatrix}R_{11}(n-1)&R_{12}(n-1)\\R_{21}(n-1)&R_{22}(n-1)\end{bmatrix}

其中 $R_{11}(n)\in \R^{(K-1)\times(K-1)}$ ，根据分块矩阵逆公式，可以得到

C^{-1}(n-1)=R_{11}(n-1)-R_{12}(n-1)R_{22}^{-1}(n-1)R_{21}(n-1)

由于 $G(n-1)$ 为对称矩阵，则对于它的逆矩阵有 $R_{21}(n-1)=R_{12}^T(n-1)$ ，所以

C^{-1}(n-1)=R_{11}(n-1)-\frac{R_{12}(n-1)R_{12}^T(n-1)}{R_{22}(n-1)}

然后构造 $G^{-1}(n)$ ，此时假设已经得到了 $C^{-1}(n-1)$ ，根据分块矩阵的求逆公式，定义 Schur 补 $\alpha=a-b^TC^{-1}b$ ，则得到

G^{-1}(n)=\begin{bmatrix}\frac{1}{\alpha}&-\frac{1}{\alpha}b^TC^{-1}\\-\frac{1}{\alpha}C^{-1}b&C^{-1}+\frac{1}{\alpha}C^{-1}bb^TC^{-1}\end{bmatrix}

此时就不需要每次都重新求逆，而是利用上一时刻的逆矩阵递推得到下一时刻的逆矩阵。之后计算 $b(n)$ 的递推。由于 $b(n)$ 中的每一项都是输入向量之间的内积，即

X^T(n)X(n-i)=\sum_{m=0}^{K-1}x(n+1-m)x(n-i-m)

这些内积也可以利用滑动窗口递推，定义 $\rho_i(n)=X^T(n)X(n-i)$ ，则

\rho_i(n)=\rho_i(n-1)+x(n)x(n-i)-x(n-K)x(n-i-K)

因此有

b(n)=\begin{bmatrix}\rho_1(n)&\rho_2(n)&\cdots&\rho_{K-1}(n)\end{bmatrix}\\a(n)=\rho_0(n)+\delta

则可计算出 $R(n)$ ，计算 $\epsilon(n)=R(n)e(n)$ ，权值更新

w(n+1)=w(n)+\mu M(n)\epsilon(n)

整体的计算复杂度约为 $O(K^2+K)$

FxFAP

将 FAP 算法应用到 ANC 中，将参考信号 $x(n)$ 和对应的向量替换为如下的 $x_f(n)$ 进行运算，其中矩阵中的每一列

x_f(n)=\hat{s}^TX(n)\\X_f(n)=\begin{bmatrix}x_f(n)\\\vdots\\x_f(n-L+1)\end{bmatrix}\\M=\begin{bmatrix}X_f(n)&X_f(n-1)&\cdots&X_f(n-K+1)\end{bmatrix}\\d_f(n)=s^TD(n)\\D_f=\begin{bmatrix}d_f(n)\\\vdots\\d_f(n-L+1)\end{bmatrix}\\y(n)=X^T(n)w(n-1)\\Y(n)=\begin{bmatrix}y(n)\\\vdots\\y(n-L+1)\end{bmatrix}

另外误差的公式为 $\hat{e}(n)=s^T\Big(D_f(n)-s^TY\Big)$ ，辅助权重的更新公式为 $\hat{w}=\hat{w}+\mu X_f(n-K+1)E(n-K+1)$

APSA

APSA (Affine Projection Sign Algorithm，仿射投影符号算法) 是一种广泛用于自适应滤波领域的算法，结合了仿射投影算法对有色输入信号的快速收敛特性和符号算法对脉冲噪声的鲁棒性，也避免了 APA 复杂的矩阵求逆运算，计算复杂度降低

APSA 基于 L1 范数优化准则，通过最小化后验误差向量的 L1 范数，同时满足权重向量变化量的 L2 范数约束来更新滤波器系数。定义先验误差向量和后验误差向量

e(n)=D(n)-M^T(n)w(n-1)\\e_p(n)=D(n)-M^Tw(n)

APSA 的设计基于一下约束优化问题

\min_{w(n)}\Vert e_p(n)\Vert_1\\s.t.\quad\Vert w(n)-w(n-1)\Vert_2^2\leq\mu^2

引入拉格朗日乘数 $\gamma>0$ ，将约束优化问题转为无约束优化问题

J(w(n))=\Vert e_p(n)\Vert_1+\gamma\Big(\Vert w(n)-w(n-1)\Vert_2^2-\mu^2\Big)

对 $w(n)$ 求偏导并令其为零

\frac{\partial J}{\partial w(n)}=-M(n)\mathrm{sgn}(e_p(n))+2\gamma(w(n)-w(n-1))=0\\\Downarrow\\w(n)=w(n-1)+\frac{1}{2\gamma}M(n)\mathrm{sgn}(e_p(n))

由于后验误差在实际中无法直接获得，通常使用先验误差近似代替。另外为其添加归一化和防止分母为 0，得到 APSA 的权重更新公式

w(n)=w(n-1)+\mu\frac{M(n)\mathrm{sgn}(e(n))}{\Vert M(n)\mathrm{sgn}(e(n))\Vert_2+\delta}

最终滤波器输出为

y(n)=X^T(n)w(n)

宽带 ANC 算法

SAF

SAF 即子带自适应滤波（Subband Adaptive Filtering），通过分析滤波器组将宽带输入信号分解为多个窄带子带信号，对每个子带独立进行自适应滤波处理，再通过综合滤波器组重构输出信号。

定义全带权向量 $w(n)=\begin{bmatrix}w_0(n)&\cdots&w_{L-1}(n)\end{bmatrix}$ ，定义第 $i$ 个子带的输入向量 $X_i(n)=\begin{bmatrix}x_i(n)&x_i(n-1)&\cdots&x_i(n-L+1)\end{bmatrix}$ ，定义第 $i$ 个子带的输出 $y_i(n)=w^T(n)X_i(n)$ ，第 $i$ 个子带的期望信号为 $d_i(n)=d_i(nM)$ ，则第 $i$ 个子带的误差信号为 $e_i(n)=d_i(n)-y_i(n)$ 。具体的 SAF 算法有如下的几种变种

SBLMS

子带 LMS 基于最小扰动原理，权重更新公式为 $w(n+1)=w(n)+\mu\sum_{i=0}^{M-1}X_i(n)e_i(n)$

NSAF

NSAF 基于最小扰动原理，在满足子带误差约束的条件下，使权向量的更新量最小

定义优化问题 $\min_{w(n+1)}\Vert w(n+1)-w(n)\Vert^2$ ，相等条件 $d_i(n)=X_i^T(n)w(n+1)$
使用拉格朗日乘数法求解得到权值更新公式 $w(n+1)=w(n)+\mu\sum_{i=0}^{M-1}\frac{X_i(n)e_i(n)}{X_i^T(n)X_i(n)+\epsilon}$
其中 $\mu$ 为步长因子，取值 $0<\mu<2$

变步长 NSAF

变步 NSAF 解决固定步长 NSAF 收敛速度稳态失调的矛盾，变步长公式如下

基于误差平方的步长控制 $\mu(n)=\mu_{min}+(\mu_{max}-\mu_{min})(1-e^{-\alpha\bar{e}^2(n)})$ $μ (n) = μ_{m i n} + (μ_{m a x} - μ_{m i n}) (1 - e^{- α \overset{e}{ˉ}^{2} (n)})$
- 其中的 $\bar{e}^2(n)=\beta\bar{e}^2(n-1)+(1-\beta)\frac{1}{M}\sum_{i=0}^{M-1}e_i^2(n)$
基于最小均方偏差 $\mu(n)=\max(\mu_{min},\min(\mu_{max},\frac{\sum_{i=0}^{M-1}\frac{\epsilon_i^2(n)}{\Vert X_i(n)\Vert^2}}{\sum_{i=0}^{M-1}\frac{e_i^2(n)}{\Vert X_i(n)\Vert^2}}))$ $μ (n) = max (μ_{m i n}, min (μ_{m a x}, \frac{\sum _{i = 0}^{M - 1} \frac{ϵ _{i}^{2} ( n )}{∥ X _{i} ( n ) ∥ ^{2}}}{\sum _{i = 0}^{M - 1} \frac{e _{i}^{2} ( n )}{∥ X _{i} ( n ) ∥ ^{2}}}))$
- 其中的 $epsilon_i^2(n)=e_i^2(n)-\sigma_i^2(n)$
- 其中 $\sigma_i(n)=\min(\sigma_i(n-1),\gamma\sigma_i(n-1)+(1-\gamma)\vert e_i(n)\vert)$
基于噪声方差估计 $\mu_i(n)=1-\frac{\sigma_i(n)}{e_i(n)}$ $μ_{i} (n) = 1 - \frac{σ _{i} ( n )}{e _{i} ( n )}$
- 其中 $\sigma_i(n)=\min(\sigma_i(n-1),\gamma\sigma_i(n-1)+(1-\gamma)\vert e_i(n)\vert)$
- 其中的 $\gamma$ 为跟踪因子，取值为 $0.99\sim0.999$

FDAF

FDAF 即频域自适应滤波器算法（Frequency-Domain Adaptive Filter），是基于块处理和快速傅里叶变换的自适应滤波算法。LMS 算法的计算复杂度为 $O(L^2)$ ，在 FDAF 通过将时域卷积和梯度计算转换为频域乘法，将复杂度降低到 $O(L\log L)$ ，同时保留了 LMS 算法的简单性和鲁棒性。算法的核心思想如下

块处理：每积累 $L$ $L$ 个新样本才更新一次滤波器的权值
- 更新公式为 $w(n+1)=w(n)+\mu\sum_{l=0}^{L-1}X(nL+l)e(nL+l)$
- 其中 $n$ 为索引，每块包含 $L$ 个样本，块梯度是对 $L$ 个瞬时梯度的平均，估计更准确，但步长上限需要降低到原来的 $\frac{1}{N}$
频域变换：利用 FFT 将时域卷积转换为频域逐点乘法，降低计算量
重叠保留法：解决频域循环卷积与线性卷积的差异问题，保证输出无混叠失真
- 直接使用 FFT 计算卷积会得到循环卷积结果，而非线性卷积
- 对长度为 $L$ 的滤波器权向量补 $L$ 个 0，得到长度为 $2L$ 的向量
- 每次处理长度为 $2L$ 的输入块，其中前 $L$ 个样本来自上一个块，后 $L$ 个是当前块
- 进行 $2L$ 点 FFT 变换到频域，逐点相乘之后做 IFFT 变换变回时域
- 只保留输出块的后 $N$ 个样本作为有效输出

约束 FDAF

约束 FDAF 通过时域约束保证权向量的长度始终为 $L$ ，避免循环卷积带来的混叠问题

W(n+1)=W(n)+\mu GX^H(n)E(n)

其中的 $G=F\begin{bmatrix}I_L&0_L\\0_L&0_L\end{bmatrix}F^{-1}$ 为约束矩阵，用于强制时域权向量的后 $L$ 个样本为 0，其中的 $F$ 为 $2L$ 点 DFT 矩阵。每次频域更新之后，需要将权向量变换回时域，将后 $N$ 个样本置 0 后再变回频域

w(n)=IFFT(W(n))\\w_{L:end}(n)=0\\W(n)=FFT(w(n))

无约束 FDAF

无约束 FDAF 移除了时域约束，计算更简单，它的权值更新公式如下

W(n+1)=W(n)+\mu X^H(n)E(n)

与上述约束 FDAF 相比，去掉了约束矩阵 $G$

归一化 FDAF

归一化 FDAF 算法为每个频点分配独立的步长，与该频点的输入功率成反比，大幅提升了对相关输入信号的收敛速度

W(n+1)=W(n)+\mu \Lambda^{-1}(n)X^H(n)E(n)

其中 $\Lambda(n)=\mathrm{diag}\{\Phi_{xx,0}(n),\Phi_{xx,1}(n),\cdots,\Phi_{xx,2L-1}(n)\}$ ，其中的 $\Phi_{xx,i}(n)$ 表示第 $i$ 个频点的输入功率估计，通过递归方式更新

\Phi_{xx,i}(n)=\lambda\Phi_{xx,i}(n)+(1-\lambda)\Vert X_i(n)\Vert^2

其中 $\lambda$ 为平滑因子，通常取 $0.9\sim0.999$

非线性 ANC 算法

基于核方法的算法

传统线性 ANC 算法基于线性时不变系统假设，通过自适应线性滤波器生成反相声波抵消噪声，但是在实际使用中，系统普遍存在非线性失真。此时线性滤波器的容量不足，降噪效果就会急剧下降，核方法通过将低维输入空间的非线性问题映射到高维再生核希尔伯特空间，再高维空间中执行线性自适应滤波，等价于原空间的非线性运算，解决了线性 ANC 的非线性建模能力不足问题

再生核希尔伯特空间 RKHS 与核技巧：设输入空间 $X\subseteq R^M$ ，存在非线性映射 $\phi:X\rightarrow H$ ，将输入向量映射到高维甚至无限维的 RKHS 中，任意两个向量的内积可以通过核函数直接计算，无需显式求解映射 $\phi$ ，这就是核技巧，即

k(x,y)=<\phi(x),\phi(y)>_H

其中的核函数需要满足 Mercer 定理，即对于任意有限个输入 $\{x_1,\cdots,x_n\}$ ，核矩阵 $K\in\R^{n\times n}$ ，其中 $K_{ij}=k(x_i,x_j)$ 是半正定的。常用的核函数如下

高斯核函数 $k(x,y)=\exp(-\frac{\Vert x-y\Vert^2}{2\sigma^2})$ ，非线性能力最强
多项式核函数 $k(x,y)=(x\cdot y+c)^d$ ，全局非线性，适合多项式型非线性
Sigmoid 核函数 $k(x,y)=\tanh(ax\cdot y+c)$ ，适合局部非线性

核 LMS 算法

设定在 $n$ 时刻的参考向量为

X(n)=\begin{bmatrix}x(n)&x(n-1)&\cdots&x(n-L+1)\end{bmatrix}

其中 $L$ 为滤波器阶数。将其映射到 RKHS 中 $\phi(X(n))\in H$ ，在 RKHS 中的权重向量为 $w(n)$ ，输出为内积形式

y(n)=<w(n-1),\phi(X(n))>_H

误差信号为

e(n)=d(n)-y(n)

采用随机梯度下降法最小化瞬时平方误差 $\xi(n)=e^2(n)$ ，梯度为

\nabla\xi(n)=-2e(n)\phi(X(n))

此时权重更新公式为

w(n)=w(n-1)+\mu e(n)\phi(X(n))

其中的 $\mu>0$ 为步长参数。根据**表示定理（**表示定理说明在高维或无限维的函数空间中，正则化优化问题的解可以表示为训练样本的有限线性组合），RKHS 中的最优权向量可以表示为所有历史输入映射的线性组合，即

w(n)=\sum_{i=1}^n\alpha(i)\phi(X(i))

其中的 $\alpha(i)$ 为展开系数，将其带入到输出公式中，利用核技巧即可得到

y(n)=\sum_{i=1}^{n-1}\alpha(i)k\Big(X(i),X(n)\Big)

权值更新转化为系数更新，即

\alpha(n)=\mu e(n)

核 FxLMS 算法

在核 LMS 算法中，未考虑到 ANC 系统中次级路径的影响，无法直接应用于 ANC 算法中，而核 FxLMS 算法是 KLMS 算法的扩展，用于解决次级路径的补偿问题

对于参考信号 $X(n)=\begin{bmatrix}x(n)&x(n-1)&\cdots&x(n-L+1)\end{bmatrix}$ ，自适应核滤波器的输出为 $y(n)$ ，通过次级路径输出 $y_f(n)=s(n)*y(n)$ ，其中 $s(n)$ 就是次级路径。得到误差信号 $e(n)=d(n)-y_f(n)$

代价函数为 $\xi(n)=e^2(n)$ ，对权向量 $w(n-1)$ 求梯度

\nabla\xi(n)=-2e(n)\Big(s(n)*\nabla y(n)\Big)

由于 $y(n)=<w(n-1),\phi(X(n))>_H$ ，则 $\nabla y(n)=\phi(X(n))$ ，带入得到

\nabla\xi(n)=-2e(n)\Big(s(n)*\phi(X(n))\Big)

但是上述中直接对高维向量 $\phi(X(n))$ 进行卷积运算不可行，所以采用近似的算法：先对低维参考信号 $X(n)$ 进行次级路径滤波，得到滤波后的参考信号 $X_f=\hat{s}(n)*X(n)$ ，其中 $\hat{s}(n)$ 为次级路径的估计。再将 $X_f(n)$ 映射到 RKHS 得到 $\phi(X_f(n))$ ，作为 $s(n)*\phi(X(n))$ 的近似。剩下的公式与核 LMS 算法相似，如下

\nabla\xi(n)=-2e(n)\Big(\phi(X_f(n))\Big)\\w(n)=w(n-1)+\mu e(n)\phi(X_f(n))\\y(n)=\sum_{i=1}^{n-1}\alpha(i)k\Big(X_f(i),X(n)\Big)\\\alpha(n)=\mu e(n)

Volterra

Volterra 算法是基于 Volterra 级数的非线性系统建模与信号处理的方法，是泰勒级数在时域的推广，能同时描述系统的非线性特性和记忆效应。Volterra 级数公式如下

y(t)=w_0+\int h_1(\tau)x(t-\tau)d\tau+\iint h_2(\tau_1,\tau_2)x(t-\tau_1)x(t-\tau_2)d\tau_1d\tau_2+\cdots

其中的 $h_n(\tau_1,\cdots,\tau_n)$ 称为 $n$ 阶 Volterra 核，是系统的高阶脉冲响应。在数字信号处理中，有离散时间的 Volterra 级数

y(n)=h_0+\sum_{p=1}^P\sum_{\tau_1=0}^{M-1}\cdots\sum_{\tau_p=0}^{M-1}h_p(\tau_1,\cdots,\tau_p)\prod_{j=1}^px(n-\tau_j)

其中的 $P$ 就是 Volterra 级数的最高阶数，而 $M$ 为系统的记忆长度，而 $h_p(\tau_1,\cdots,\tau_p)$ 就是离散时间下 $p$ 阶 Volterra 核。通常在实际使用中一般选择 2-3 阶，而且当阶数越高，所需要的参数数量就越多。对于 2 阶的 Volterra 核，所需要的参数数量为 $C(M,2)$ ，将所有核系数展开为一个权向量，所有输入项展开为一个输入向量，则输出可以表示为

y(t)=w^TX(n)

其中的输入向量为

X(n)=\begin{bmatrix}1\\x(n)\\\vdots\\x(n-M+1)\\x(n)x(n)\\\vdots\\x(n)x(n-M+1)\\x(n-1)x(n-1)\\\vdots\\x(n-M+1)x(n-M+1)\end{bmatrix}\\w=\begin{bmatrix}h_0\\h_1(0)\\\vdots\\h_1(M-1)\\h_2(0,0)\\\vdots\\h_2(0,M-1)\\h_2(1,1)\\\vdots\\h_2(M-1,M-1)\end{bmatrix}

Volterra 算法的核心是估计核函数的系数，常用方法如下

最小二乘 LS

最小化误差的平方和 $J(w)=\sum_{n=0}^{N-1}e^2(n)=\sum_{n=0}^{N-1}\Big(d(n)+w^TX(n)\Big)^2$
得到最优解 $w=-(X^TX)^{-1}X^Td$ ，其中的 $X$ 为输入矩阵， $d$ 为期望输入向量

最小均方误差 LMS

最小化误差的平方和 $J(w)=\sum_{n=0}^{N-1}e^2(n)=\sum_{n=0}^{N-1}\Big(d(n)+w^TX(n)\Big)^2$
基于梯度下降法更新权向量 $w=w-\mu e(n)X(n)$

递归最小二乘 RLS

最小化加权误差平方和 $J(w)=\sum_{n=0}^{N-1}\lambda^{N-n}e^2(n)$
计算增益向量 $K(n)=\frac{P(n-1)X(n)}{\lambda+X^T(n)P(n-1)X(n)}$
计算误差信号 $e(n)=d(n)+w^T(n-1)X(n)$
更新权向量 $w(n)=w(n-1)-K(n)e(n)$
更新逆相关矩阵 $P(n)=\lambda^{-1}\Big(P(n-1)-K(n)X^T(n)P(n-1)\Big)$

前言