这一篇在干嘛？
会”学习”的滤波器：LMS/NLMS/符号 LMS、变换域 LMS、RLS 与卡尔曼算法，以及 PCA/ICA 等现代信号处理技术。 原书代码为 VHDL，本篇所有代码已改写为 Verilog。
- 第8章
- 自适应系统
- 8.1 自适应系统的应用
- 8.2 最优估计技术
-
- 信号的性质
-
- 成本函数的定义
- 最优维纳估计
- 例 8.1 二抽头 FIR 滤波器干扰的消除
- 8.3 Widrow-Hoff 最小二乘法算法
- 算法 8.2 Widrow-Hoff LMS 算法
- 例8.3 步长的边界
- 例8.4 标准化LMS
- 8.4 变换域 LMS 算法
- 8.5 LMS算法的实现
- 例 8.5 2 抽头自适应 LMS FIR 滤波器
-
- 采用 FFT 的模块变换
-
- 延迟的 LMS 算法
-
- 先行 DLMS 算法
- 例 8.6 2 抽头流水线自适应 LMS FIR 滤波器
- 例8.7 Signum LMS滤波器中的误差底线
- 8.6 递归最小二乘法算法
- 算法8.8 RLS算法
- 例8.9 RLS学习曲线
- 算法8.10 快速卡尔曼RLS算法
- 8.7 LMS 和 RLS 的参数比较
- 8.8 主成分分析(PCA)
- 例8.11 ECG数据降维
- (1) 直接方法
- (3) 神经网络学习方法
- 例8.12
- 例 8.14 基于 Sanger GHA 的 PCA 设计
- 8.9 独立成分分析(ICA)
- 例8.15 采用EASI算法的ICA设计
- 8.10 语音和音频信号编码
- 例8.16 A律设计
- 例8.17 ADPCM设计
- 8.11 练习
- (e) 对于三阶系统重复(a)~(d)。
- 8.7 假定期望信号如下：
- 8.8 假定期望信号如下：
自测一下
LMS 算法的更新公式里 μ 起什么作用？
μ 是步长：太大会发散（超过输入功率倒数的 2 倍），太小收敛慢。NLMS 用输入功率归一化 μ，对输入幅度不敏感，实际设计首选。
符号-LMS（SLMS）省在哪？
把误差或数据的乘法换成”取符号后加减”，乘法器变成加减器。代价是收敛略慢、稳态误差稍大，适合资源紧张且精度要求不高的场景。
RLS 与 LMS 的本质差别？
RLS 用逆相关矩阵递推，每步都精确最小化最小二乘代价，收敛快一个数量级且对特征值扩散不敏感；代价是 O(N^2) 运算量与数值稳定性问题。
变换域 LMS 为什么能加速收敛？
对输入做正交变换（如 DFT）把各频带去相关，等效把大的特征值扩散度压平，再用归一化步长——各模式可独立用大步长，收敛显著加快。
自适应系统：让滤波器自己”学会”工作
到目前为止，我们遇到的滤波器都有一个共同点：系数一旦设计好就固定不变。这在很多场合是够用的，但现实中的信号往往没那么”听话”——语音处理、通信、雷达、声呐、地震勘探、生物医学等领域的信号，其统计特性会随时间缓慢变化。如果还用固定系数的滤波器，性能就会逐渐变差。解决办法就是自适应数字滤波器（Adaptive Digital Filter, ADF）：它能根据输入信号实时调整自己的系数，让滤波器始终”逼近”当前条件下的最优解。前提是参数变化相对采样频率足够慢，这样我们才有时间不断计算”更好的”系数估计值并加以调整。
从结构上看，任何前面学过的滤波器结构（FIR 或 IIR）原则上都可以当自适应滤波器用，但选择时要注意：
- FIR 滤波器：图3-1 的直接形式最有利，因为所有系数可以同时、并行地更新，而自适应算法恰恰需要每一拍都更新全部系数。
- IIR 滤波器：图4-12 的网格（lattice）结构是较好的选择，因为它的定点舍入误差灵敏度低，而且系数稳定性控制简单。
不过实践表明，FIR 滤波器在实际应用中远比 IIR 成功（IIR 的极点会带来稳定性隐患），所以本章主要讨论自适应 FIR 滤波器的高效实现。理论上，FIR 自适应算法应收敛于由维纳-霍夫方程给出的最优非递归估计器；我们还会顺带提到最优递归估计器（卡尔曼滤波器），并从计算复杂度、算法稳定性、初始收敛速度、收敛一致性以及对附加噪声的鲁棒性等角度比较不同方案。本章最后还会介绍更高级的话题：用主成分分析（PCA）和独立成分分析（ICA）做多通道信号的盲源分离（BSS），以及从 ADPCM 到 MP3 的语音处理。
自适应系统如今已是成熟的 DSP 分支，早在 20 世纪 80 年代中期就有专著出版[268273]，更新的结果见教科书[274276]；期刊（如 IEEE Transactions on Signal Processing）则持续关注 LMS 稳定性及其各种变体等基础问题。
8.1 自适应系统的应用
自适应滤波器的应用五花八门，但仔细一看，几乎全都能归入下面 4 种配置之一：
- 干扰消除：把噪声”减掉”
- 预测：用过去猜现在
- 反演模拟：求一个未知系统的”逆”
- 辨识：求一个未知系统的”拷贝”
为了讲清这些配置，统一使用如下记号：
- = 自适应系统的输入
- = 自适应系统的输出
- = （自适应系统的）期望响应
- = 估计误差
8.1.1 干扰消除
这是应用最广的一类。输入信号里既混着有用信息，也掺杂干扰——比如随机白噪声，或者 50/60 Hz 的电力线交流噪声。图 8-1 给出了配置：传感器收到的信号 （信号+干扰）与自适应滤波器对基准信号 的响应 相减，得到误差信号 ，它同时就是整个系统的输出。关键在于：收敛之后， 会变成干扰的”加性逆”（大小相等、符号相反），于是干扰就从输入中被抵消了，而有用信号基本不动。
图 8-1 用于干扰消除的基本配置
后面我们会用一个电力线交流噪声的完整例子详细研究这一配置。另一个常见应用是电话系统中的自适应回声消除；干扰消除还被用于天线阵列（即所谓波束形成器），可以自适应地抑制来自未知方向的噪声干扰。
8.1.2 预测
预测的任务是给出随机信号当前值的最佳估计（通常在最小均方误差意义下）。显然，只有当输入信号与白噪声”长得不一样”（即有相关性可利用）时，预测才可能成功。图 8-2 说明其原理：期望信号 经一个延迟后作为自适应滤波器的输入，同一信号 本身还充当期望响应来计算误差 。换句话说，滤波器在用”过去”预测”现在”。
图 8-2 用于预测的框图
预测编码已成功用于图像和语音处理：不直接对信号编码，只编码预测误差，因为误差通常小得多、更省存储和带宽。其他应用包括功率谱建模、数据压缩、谱增强与事件检测[269]。
8.1.3 反演模拟
反演模拟的任务是构造未知时变被控对象（plant）的逆模型。典型例子是通信均衡：信道（多径传播）把信号”搅乱”了，自适应滤波器学习去近似这个信道的反演，把信号”复原”。图 8-3 的配置是：输入 先进入未知被控对象，其输出 作为自适应滤波器的输入；延迟后的原输入 用作期望响应来计算误差 并调整系数。收敛后，自适应滤波器的传递函数就近似等于未知被控对象传递函数的倒数。
图 8-3 阐述反演系统建模的示意图
除了通信均衡，反演模拟还用于提高被窄带噪声污染系统的信噪比（S/N）、自适应控制系统，以及语音分析中的去卷积和数字滤波器设计[269]。
8.1.4 系统辨识
系统辨识与反演模拟正好相反：目标是让自适应滤波器直接复制未知系统。结构如图 8-4：同一时间序列 同时送入自适应滤波器和未知被控对象，未知对象的输出 作为期望响应。收敛后 以最小均方意义逼近 。如果自适应滤波器的阶数与未知对象匹配，且输入 是广义稳态（WSS）的，则自适应滤波器的系数将收敛到与未知对象相同的值。实际中未知对象的输出往往带观测噪声、结构也不完全匹配，性能会偏离理想，但这不妨碍它成为评估自适应算法性能的”标准实验台”——后面对比 LMS 与 RLS 时就用这个配置。
图 8-4 系统辨识的基本配置
系统辨识还用于生物学建模、社会与商业系统仿真、自适应控制、滤波器设计和地球物理学[269]；在地震勘探中用于建立分层地球模型，解释地表信号的复杂性[268]。
8.2 最优估计技术
信号的性质：算法能工作的前提
任何自适应算法都有收敛条件，而条件是建立在对输入信号统计特性的假设上的。从概率角度看，输入信号可视为随机变量向量。两个基本假设：
（1）各态历经性（ergodicity）：用一段信号算出的统计量，应等于”所有可能信号”的统计平均。也就是说，时间平均可以代替集合平均，例如均值
和方差
（2）广义稳态（WSS）：均值、方差不随时间变化，且自相关函数
只依赖于时间差 。特别地：
即 就是 WSS 过程的平均功率，后面确定步长边界时会反复用到它。
一些高级算法还需要更高阶的统计量，最常用的是四阶统计”峰态”（kurtosis）。这里采用 MATLAB 的定义（不做 -3 标准化）：
这个定义有个好性质：对信号做直流偏移或比例缩放都不改变峰态值，即
展开计算一个例子：求 上均匀分布信号的峰态。其概率密度函数（PDF）在区间内恒为 ，均值为 。先算四阶矩和二阶矩：
于是峰态：
注意 在计算中约掉了——这正是式(8-3)所说的”比例缩放不影响峰态”。
峰态是区分信号”形状”的利器。以 MATLAB 标准，高斯分布的峰态是 3。峰态高于 3、分布带尖峰脉冲特性的信号叫超高斯信号：例如长度为 、只含 个非零脉冲的信号，其 ，脉冲越稀疏峰态越大。峰态低于 3 的叫次高斯信号，如正弦、余弦、三角波和二进制 ±1 “掷硬币”信号。表 8-1 汇总了典型分布的三个统计量。
表 8-1 典型分布的均值、方差和峰态
| 分布 | 均值 | 方差 | 峰态 |
| 掷硬币分布±1 | 0.0 | 1.0 | 1.0 |
| 正弦/余弦,周期>8 | 0.0 | 0.5 | 1.5 |
| 三角分布[0,10] | 5.0 | 35.0 | 1.8 |
| 均匀分布[-0.5,0.5] | 0.0 | 1/12 | 1.8 |
| 0、1、-1、0、1、-1、... | 0.0 | 0.5 | 2.0 |
| 高斯/标准分布 | 0.0 | 1.0 | 3.0 |
| 瑞利噪声 | 1.25 | 0.43 | 3.3 |
| 拉普拉斯噪声 | 0.0 | 2.0 | 5.8 |
| 长笛 | 0.0 | 1.0 | 3.59 |
| 钢琴 | 0.0 | 1.0 | 3.94 |
| 男声 | 0.0 | 1.0 | 16.4 |
| 女声 | 0.0 | 1.0 | 46.1 |
| I脉冲/长度L | 0.0 | 1.0 | L/I |
可以看到真实信号（男声 16.4、女声 46.1）峰态极大——它们是很”尖”的超高斯信号，这一性质在后面盲源分离章节会直接利用。
成本函数的定义：给误差”打分”
所有自适应算法的核心问题都是：误差 该怎么”加权”成一个人人比得高低的目标？定义
其中 是待估计的随机变量， 是自适应滤波器给出的估计值。最常用的成本函数是最小均方（Least-Mean-Square, LMS）型：
但这不是唯一选择。图 8-5 给出了三种误差成本函数：除了二次函数外，还有绝对值误差和非线性阈值型。如果某个误差水平以内都”可以接受”，用阈值型函数可以显著降低算法计算量——有趣的是，Widrow 最初提出的自适应算法用的恰恰就是阈值函数处理误差。
图 8-5 3 种可能的误差成本函数
不过，二次误差函数有一个决定性优势：它是系数的光滑二次函数，可以借助连续信号领域早已成熟的维纳-霍夫关系构造随机梯度算法。下一小节回顾维纳估计，它直接引出 Widrow 等人提出、至今应用最广的 LMS 算法[277, 278]。
最优维纳估计：理论上的”满分答案”
自适应 FIR 滤波器的输出是一个卷积和：
我们的目标是调整系数 使成本函数 最小。用向量记号写卷积更方便：
其中 ， 都是 向量， 表示转置（复数情形为厄米共轭转置）。对矩阵 ，转置沿主对角线”镜像”：。
把式(8-8)代入误差定义：
于是均方误差展开为：
注意 是系数的二次函数——把它画出来是一个永远非负的”凹超抛物面”，两个系数的情形如图 8-6 所示。自适应的过程，就是沿着这个碗状曲面”往下走”，走到碗底。为此对 微分并令梯度为零：
图 8-6 两个分量情况下的误差成本函数，成本函数的最小值在 和 处
假定滤波器权重向量 与信号向量 统计独立（不相关），就得到：
即最优系数向量：
为书写简洁，定义两个统计量。一是输入的 自相关矩阵：
它沿对角线对称，具有托普利兹（每条对角线上元素相等）结构。二是期望信号与基准信号之间的 互相关向量：
于是维纳解可以写成极简形式：
这就是著名的维纳-霍夫方程[267]。式(8-12)有唯一解的条件是 存在，即自相关矩阵非奇异（行列式非零）；幸运的是，对 WSS 信号可以证明 必非奇异[268]。最优估计处的”碗底”残差为：
其中 是 的方差。
例 8.1 二抽头 FIR 滤波器消除电力线干扰
来看一个把上述理论算到底的完整例子。观测到的通信信号由三部分组成：
- 曼彻斯特编码的传感器信号 ，幅值 （数据，图 8-7(a)）；
- 附加的高斯白噪声 （图 8-7(b)）；
- 60 Hz 电力线交流噪声（美国电网频率，中国为 50 Hz），幅值 （图 8-7(c)）。
采样频率取交流噪声频率的 4 倍，即 Hz。观测信号为：
自适应滤波器可用的基准信号 （图 8-7(d)）取自同一个电力线干扰的参考源，只是带一个常数相位偏移 ：
(a) 数据
(b) 噪声
(c) 交流噪声
(d) 基准信号
图 8-7 在例 8.1 的电力线交流噪声中用到的信号
二抽头滤波器的输出为：
第一步：算基准信号的自相关。延迟 0 和 1：
（余弦平方的平均值是 1/2。）而
（相位差 90° 的正交项平均为零——这解释了为何 会是对角阵。）
第二步：算互相关。数据项 和噪声项 与确定性余弦不相关，平均后消失，只剩交流噪声贡献：
第三步：代入维纳-霍夫方程，解 线性方程组：
逐项核对：，，即最优系数为 、，恰好是图 8-6 抛物面碗底的位置。
第四步：看效果。图 8-8(a) 是三种信号的叠加（曼彻斯特编码的 5 位有效信号、60 Hz 交流噪声、白噪声），干扰把数据几乎完全淹没了；图 8-8(b) 是系统输出 ——交流噪声被干净地消掉，数据信号清晰可见。这就是”两根抽头，胜过一台昂贵的陷波器”的效果。
(a) d[n]=数据+噪声+交流噪声
(b) 系统输出 e[n]
图 8-8 用最优维纳估计方法消除曼彻斯特编码的数据信号中的 60Hz 电力线干扰
用 Verilog 看 LMS 抽头长什么样
本节素材原书没有给出具体 VHDL 代码清单，但为了后续 FPGA 实现做铺垫，这里把书中描述的自适应 FIR 抽头（乘累加 + 系数更新）改写成可综合的 Verilog HDL 参考实现（单抽头，全并行结构，系数更新使用移位近似除法避免浮点）：
// 原书 VHDL 改写为 Verilog：单抽头自适应 FIR（LMS 权值更新）
module lms_tap #(
parameter W = 16, // 数据位宽（Q15 定点）
parameter MU_SHIFT = 3 // 步长 mu = 2^(-MU_SHIFT)，避免乘除法
)(
input wire clk,
input wire reset,
input wire signed [W-1:0] x_in, // 基准信号输入
input wire signed [W-1:0] d_in, // 期望响应
output wire signed [2*W-1:0] y_out, // 滤波器输出 y = f*x
output reg signed [W-1:0] f // 本抽头系数
);
wire signed [2*W-1:0] prod = x_in * f; // y 分量
wire signed [W-1:0] e = d_in - prod[2*W-2:W-1]; // 误差 e = d - y
// 系数更新: f <= f + mu * e * x, 用算术右移实现乘 mu
wire signed [2*W-1:0] grad = e * x_in;
always @(posedge clk) begin
if (reset)
f <= {W{1'b0}};
else
f <= f + grad[2*W-1:MU_SHIFT]; // 近似 mu*e*x
end
assign y_out = prod;
endmodule要点有三：乘累加给出 ； 得到误差；系数更新 用右移（）实现，这样 FPGA 里一个乘法器都不用加。多抽头时把该模块级联、把 一路传给每个抽头即可。
8.3 Widrow-Hoff 最小二乘法（LMS）算法
为什么不直接用维纳公式
维纳-霍夫方程虽然给出”满分答案”，直接计算却不现实，原因有三：
- 统计量难获得：构造 和 本身计算量就大，而且要多长的数据才算”足够”用于估计统计量？事先并不知道。
- 矩阵求逆太贵：即使有了 ，还要求逆。滤波器阶数一大，求逆非常耗时。
- 精度问题：即便算出了 ，中间步骤繁多，特别是用定点运算时，累积误差可能让结果不可信。
Widrow-Hoff 最小均方（LMS）自适应算法[277]巧妙地绕开了所有这三个障碍：它既不需要严格度量相关函数，也不做任何矩阵求逆，只用输入信号的实时测量值工作——代价是精度受单次样本统计波动的影响。
算法推导：用瞬时梯度代替真实梯度
LMS 是最速下降法的一种实现：下一个系数向量等于当前系数向量加上一个正比于负梯度的修正量：
参数 是学习因子（步长），控制稳定性与收敛速度。问题是真实的梯度 需要 ，而期望我们算不起。LMS 的”粗但有效”的招数是：用单次误差平方 的梯度去估计均方误差的梯度。真实梯度为：
估计梯度为：
由式(8-9)， 对系数的偏导恰为 ，故：
代回式(8-14)， 与 2 正好抵消，得到 LMS 的核心更新公式：
算法 8.2：Widrow-Hoff LMS——调整 个系数的全部步骤：
- 初始化 向量 。
- 接收一对新采样 ，并把 移入基准信号向量 。
- 计算 FIR 输出： (8-19)
- 计算误差： (8-20)
- 更新系数： (8-21)
- 回到第 2 步。
注意一个漂亮的性质：虽然算法名叫”均方”，但它不需要平方、不求均值、不微分——每拍只多一次乘加循环，因此极易实现（MATLAB 代码见[276]，C 代码见[279]，PDSP 汇编见[280]）。
用例 8.1 相同的配置对不同步长 仿真（自适应滤波器 1 秒后启动），结果如图 8-9：左侧是系统输出 ，右侧是系数的调整轨迹。可以看到系数最终逼近例 8.1 算出的最优值 、。
图 8-9 运用 LMS 算法在 3 个不同的步长 下对电力线干扰消除的仿真结果，(左侧)系统输出 e[n]，(右侧)滤波器系数
稳定性：步长 的上限
当输入是 WSS 信号时（这正是维纳估计要求的条件），LMS 的梯度估计是无偏的，权向量的期望值收敛到维纳解(8-12)。只要 ，算法就从任意初始值稳定收敛。图 8-10 用另一种方式展示收敛：把系数轨迹投影到 平面上，并叠加误差等高线。注意 LMS 走的是锯齿形路线而非垂直于等高线的真实梯度方向——因为它用的是带噪声的瞬时梯度估计。
图 8-10 电力线干扰例题中的收敛说明(采用 2D 的等高线图 )
确定 的经典方法借助 的特征值。求解特征方程：
（ 为 单位矩阵）得到 个特征值 ，满足：
特征值分析给出的平均意义稳定条件是：
但要清醒：这个条件只保证均值收敛，并不保证有限的均方误差——系数 可能虽然趋于 ，方差却可能发散。所以实践中需要更严格的限制（见下文）。
收敛速度与特征值比（EVR）
图 8-9 中学习曲线呈指数衰减不是偶然。通过特征值分析可以把系数分解成若干互不相关的”模式”，模式数等于自由度（独立分量）数，本例中即系数个数。第 个模式的时间常数为：
最慢的模式对应最小特征值：
把稳定上限(8-25)代入，得到一个重要结论：
也就是说：自相关矩阵的特征值比（EigenValue Ratio, EVR） 越大，LMS 收敛越慢。文献[275]有符合该结论的仿真，8.3.1 节将专门讨论。
还有一个工程”经验法则”[276]：实际选 时，取理论值的 1/10 以下更保险。理论依据是 Horowitz 和 Senne[281]以及 Feuer 和 Weinstein[282]（方法不同、结论一致）给出的严格条件（假定输入向量各元素统计独立）：
这两个条件无法解析求解，但可以证明它们非常接近下面这个实用上极方便的边界：
它的妙处在于：由式(8-24)， 的迹就是基准信号的总平均输入功率——这个量在硬件里用一个累加器就能估计出来，根本不需要算特征值！
例 8.3 步长的边界
继续用例 8.1 的数据，两种方法各算一遍。先按式(8-25)求特征值：
于是经典边界给出：
而按更严格的实用边界(8-31)（，）：
图 8-11 的仿真验证了理论的预言： 时算法发散， 时收敛。同时还揭示了一个细节：即使 收敛更快，系数也在最优值附近”来回波动”（噪声驱动下不精细收敛）。想要平滑、稳定的逼近，必须把 再取小些——这正是”经验法则取理论值 1/10”的由来。
图 8-11 基于最大步长值的 LMS 算法的电力线干扰消除的仿真结果：(左侧)系统输出 e[n]，(右侧)滤波器系数
关于条件(8-29)(8-30)的假设——输入各元素统计独立——要区分场合：如果输入来自 个独立传感器的天线阵列，假设成立；但对抽头延迟线结构的 ADF，各输入显然高度相关，此时 Butterweck[268]指出长滤波器的稳定边界可放宽为：
即比式(8-31)分母上放宽 3 倍。但它只适用于长滤波器，所以实际使用少于式(8-31)。
8.3.1 学习曲线
学习曲线——误差函数 随迭代次数的变化——是比较不同算法和配置性能的标准工具。这里用系统辨识配置（图 8-12）研究 LMS 对两个参数的敏感性：待辨识系统的特征值比（EVR）和信噪比（S/N）。自适应滤波器长度 ，与”未知”系统等长；附加噪声取两个水平：高噪声 dB，低噪声 dB（相当于 8 位量化）。
图 8-12 LMS 学习曲线的系统辨识结构
要制造不同的 EVR，需要一个有色噪声源：用 的高斯白噪声通过一个三抽头对称 FIR 滤波器 （一阶 IIR 生成一阶马尔可夫过程的方法见练习 8.10）。FIR 的好处是功率容易标准化：要求系数满足 ，使输出与输入功率相同，即：
用这样的滤波器可为 生成不同的 EVR（图 8-13）；表 8-2 列出 时四个 10 的幂次 EVR 对应的滤波器系数。
图 8-13 不同系统规模 L 的三抽头滤波器系统的特征值比
表 8-2 四种不同噪声形状的 FIR 滤波器对于 L=16 生成 10 的幂形式的特征值比
| 序号 | 脉冲响应 | EVR |
| 1 | $0+1z^{-1}+0.0z^{-2}$ | 1 |
| 2 | $0.247665+0.936656z^{-1}+0.247665z^{-2}$ | 10 |
| 3 | $0.577582+0.576887z^{-1}+0.577582z^{-2}$ | 100 |
| 4 | $0.432663+0.790952z^{-1}+0.432663z^{-2}$ | 1000 |
对高斯白噪声源， 是对角阵，所有特征值相同（，EVR=1）；其他 EVR 可用 MATLAB 验证（练习 8.9）。“未知”系统取一个奇对称 FIR，系数为 （图 8-14(a)）。按式(8-31)：
为稳妥取一半：。学习曲线（系数误差）用标准化对数式计算：
结果见图 8-14：EVR=1 时（图 8-14(b)），200 次迭代后自适应滤波器已无误差地学会了未知系统的系数。而学习曲线（50 次以上求均值，图 8-14(c)(d)）显示 LMS 对 EVR 极其敏感：EVR 一大就需要多得多的迭代。遗憾的是现实信号 EVR 普遍很高——例如语音信号的 EVR 达 1874[283]。好消息是（图 8-14(c)）LMS 在高噪声环境下依然能正常自适应，鲁棒性不错。
(a) 未知系统
(b) EVR=1
(c) S/N = -10dB
(d) S/N=-48dB
图 8-14 基于图 8-12 所示的系统辨识结构的 LMS 算法的学习曲线：(a) “未知”系统的脉冲响应 ，(b) 随时间学习的系数，(c) 对于较大系统噪声的 50 次以上平均学习曲线，(d) 对于较小系统噪声的 50 次以上平均学习曲线
8.3.2 标准化 LMS（NLMS）
LMS 用的是固定步长 ，其稳定上限 依赖信号功率 。问题在于：这个统计量可能随时间变化，用固定 就不可能始终贴近最优。NLMS 的思路是每拍现场估计功率、现场算步长。当前数据向量的瞬时功率估计为：
代入稳定边界即得”标准化”步长：
如果担心分母偶尔变得很小导致 过大（爆炸），加一个小常数 保护：
出于安全，一般不直接用 ，而用它的一半，如 。
例 8.4 标准化 LMS：仍用图 8-12 的系统辨识结构，但输入改为带干扰的脉幅调制（PAM）信号（图 8-15(a)）。传统 LMS 需要先测出 ，得 并固定使用；NLMS 则按瞬时功率 动态调整。看图 8-15(c) 的 曲线：基准信号幅度大时步长自动变小，幅度小时步长自动变大——于是图 8-15(b) 中系数在小信号区间能迈出更大的学习步子。50 次以上的平均学习曲线见图 8-15(d)。虽然受干扰 PAM 的 EVR 超过 600，NLMS 在收敛性上依然明显优于固定步长的 LMS。
(a) 自适应滤波器和”未知”系统的基准信号输入
(b) NLMS 随时间学习的系数
(c) 用于 LMS 和 NLMS 的步长
(d) 50 次以上学习曲线的均值
图 8-15 采用图 8-12 所示的系统辨识结构的标准化 LMS 算法的学习曲线
不过式(8-41)的功率估计只是当前向量的一幅”快照”，瞬时值可能偏小，带来过大的 。更稳的做法是用递归平均，让功率估计带更长的记忆：
其中 小于 1 且接近 1。对图 8-15 那类剧烈起伏的信号， 要小心选： 太小，NLMS 的行为就退化回原始 LMS（练习 8.14）。
8.4 变换域 LMS 算法
把系数调整搬到变换域中做，有两个动机：
- 省计算：借助 FFT 的快速循环卷积，以”模块更新”方式批量计算滤波输出和系数调整，摊薄每次运算的成本[284]。
- 提速度：寻找一种变换让自适应滤波器的各个模式”解耦”（互不相关），等价于把 EVR 压下来，从而加快收敛[283, 285]。
8.4.1 快速卷积技术
FIR 滤波可以用 FFT 快速循环卷积实现；对自适应滤波器，它带来面向模块（block）的处理方式。模块大小原则上任意，但通常取滤波器长度的两倍——这样系数更新的时延不至于太长，计算量上也往往是好选择。流程是：一块 个输入值经变换后与滤波器系数卷积，产生 个新输出 ；据此算出 个误差 ；再在变换域中用变换后的输入序列更新系数。下面用 走一遍。一个块内计算 3 个输出：
这恰好是长度 6 的循环卷积：
误差信号依次为：
梯度也按块处理：
每个系数的更新把块内三个瞬时梯度累加起来：
仔细看，这又是一个循环卷积——只是这次 序列是逆序出现的：
时域逆序在傅立叶域里对应取 的共轭变换，所以实现时用共轭 FFT 即可。系数按块更新：
全部步骤汇总于图 8-16：两次 FFT（正向变换输入、共轭变换误差）、频域相乘、一次逆 FFT 得输出与梯度、频域更新系数。
图 8-16 采用 FFT 的快速变换域滤波方法
稳定性方面，模块化更新引入的延迟不可忽略。Feuer[286]指出：每 步才更新一次系数时，步长必须降低为：
与单样本 LMS 的边界(8-31)相比，二者非常接近：只有当模块很大（）时 的差异才显著；取 时式(8-46)退化为式(8-31)。但要记住：BLMS 的时间常数是以”块”为单位的——块处理 LMS 的最大时间常数是普通 LMS 的 倍。换句话说，快速卷积省了计算量，却用收敛速度付出了代价：天下没有免费的午餐。
8.4.2 应用正交变换
在 8.3.1 节我们已经看到，LMS 算法的收敛速度对输入信号自相关矩阵的特征值比（Eigen Value Ratio, EVR）高度敏感。EVR 越大，误差曲面的”峡谷”越窄，收敛越慢。遗憾的是，许多实际信号的 EVR 都很高——例如语音信号的 EVR 可高达 1874。如果直接用时域 LMS 去处理这样的信号，滤波器会”学”得非常慢。
变换域算法提供了一条出路：它把输入信号解耦成一组在统计上更加独立的”模式”，等效于把误差曲面”拉圆”了。理论上最优的变换是 Karhunen-Loève 变换（KLT），它能完全对角化自相关矩阵，但 KLT 需要事先知道信号统计特性，不是一种可实时执行的变换。次优而实用的选择是离散余弦变换（DCT）、快速傅立叶变换（FFT），以及沃尔什、阿达马、哈尔等正交变换，它们都能显著降低 EVR，且大多有快速算法。
现在尝试用这个概念加速 8.3.1 节的系统辨识实验：自适应滤波器要学习一个 16 抽头 FIR 滤波器的脉冲响应。做法是——输入向量 先经过一个 的正交变换 （这里是 DCT）变成 ，滤波器系数 在变换域中做 LMS 更新；为了让读者能观察到”原域中”的滤波器长什么样，可以在收敛后再做一次 IDCT（实际应用中 IDCT 只需算一次，不必每个采样点都算）。下面的 MATLAB 代码演示了变换域 DCT-LMS 算法：
for k = L:Iterations % adapt over full length
x = [xin;x(1:L-1)]; % get new sample
din = g'*x + n(k); % "unknown" filter output + AWGN
z = dct(x); % LxL orthogonal transform
y = f' * z; % transformed filter output
err = din-y; % error: primary - reference
f = f + err*mu.*z; % update weight vector
fi = idct(f); % filter in original domain
J(k-L+1) = J(k-L+1) + sum((fi-g).^2); % Learning curve end变换 对特征值分布的改变可以用下式描述：
其中上标 表示转置共轭。
还有一个关键问题没解决：变换后 个”模式”（即变换后的输入 ）的功率在统计上各不相同，原来在原域里选好的步长 可能既不再保证稳定，也不再带来快速收敛。Lee 和 Un 的仿真表明：如果变换域中不做功率标准化，变换域 LMS 的收敛性相比时域 LMS 并没有实质改进。道理很直观——按稳定性条件（8-31），每个频谱分量的稳定性边界取决于它自己的功率，所以应当为这 个分量分别计算各自的步长上限：
也就是只用变换分量的功率 来缩放步长。这其实就是功率标准化：功率大的分量用小步长，功率小的分量用大步长，思路与 8.3 节讨论的归一化 LMS（式 8-42）类似。功率标准化之后，变换 对特征值分布的总效果变为：
其中 是一个对角矩阵，作用是把 的对角元素全部归一化为 1。上面 MATLAB 代码中的 err*mu.*z 正是逐分量（.* 表示逐元素相乘）的功率标准化更新。
图 8-17 给出了 的 FIR 滤波器在 4 个不同特征值比下每个频谱分量的最优步长。图 (a) 中输入是纯高斯白噪声，EVR=1，所有频谱分量功率相等，因此步长也几乎一致；其余几幅图中滤波器对不同频段做了放大或衰减，对应分量的步长就被调低或调高——这相当于对步长做了”噪声整形”。
图 8-17(a) EVR=1 时 DCT-LMS 各频谱分量的最优步长
图 8-17(b) EVR=10 时 DCT-LMS 各频谱分量的最优步长
图 8-17(c) EVR=100 时 DCT-LMS 各频谱分量的最优步长
图 8-17(d) EVR=1000 时 DCT-LMS 各频谱分量的最优步长
图 8-18 展示了功率标准化 DCT-LMS 的学习效果：即便 EVR 高达 1000，学习曲线仍然保持收敛；只是在约 -48dB 的误差底线附近，高 EVR 与低 EVR 情形在达到最低误差的水平和一致性上略有差别。可以说，变换加功率标准化把”坏输入”变得可用了。
图 8-18 4 个不同特征值比的 DCT 变换域 LMS 算法学习曲线（50 次循环的均值）
为实时应用挑选变换时，还要考虑计算复杂度。DCT、DST 这类实时变换优于 FFT 这类复杂变换；有快速算法的变换优于没有快速算法的变换；而哈尔、阿达马这类整数变换完全不需要乘法，是硬件实现的最佳选择。最后别忘了 RLS 算法（稍后讨论）也是一种替代方案：它比 LMS 复杂，但比变换域方法更有效，收敛速度可与基于 KLT 的 LMS 相媲美。
8.5 LMS 算法的实现
前面讨论的 LMS 算法都可以用 HDL 直接实现，但硬件化之后必须额外回答两个问题：一是有限字长带来的量化效应是否可以接受；二是如何用流水线提高吞吐量，同时保证自适应滤波器仍然稳定收敛。下面逐一解决。
8.5.1 量化效应
在动手写硬件之前，应先确认参数和数据落在”绿灯区”。最简单的验证办法是把软件仿真从全精度换成目标整数精度。图 8-19 给出了 8 位整型数据、步长 分别取 1/4、1/8、1/16 的仿真结果。注意两点矛盾： 不能太小，否则系数更新方程（8-21）中由 缩放后的梯度 幅度过小，算法无法在有限次迭代中收敛到最优值 、；而稳定性又给 设了上限。 越小收敛问题越多。工程上的解决办法是给系统增加额外的小数位，让”小 “不再等于”小数值”——这就是把 取成 、用移位实现缩放的原因。
图 8-19 采用整型数据的 LMS 算法电力线干扰消除仿真结果：(左)误差输出 e[n]，(右)滤波器系数
8.5.2 LMS 算法的 FPGA 设计
图 8-20 用信号流图给出了 LMS 算法的一种可行硬件结构。从资源角度看：需要缩放 ，且除滤波部分的 个乘法器外，系数更新路径还要 个乘法器，共 2L 个通用乘法器——工作量是第 3 章例 3.1 可编程 FIR 滤波器的两倍多。这正是 LMS 的硬件代价：自适应能力是”买”来的。
图 8-20 LMS 算法的信号流图
下面用两个例子把结构落到 FPLD 上。
例 8.5：2 抽头自适应 LMS FIR 滤波器
原书 VHDL 改写为 Verilog。下面的设计有两个系数 、，步长 （通过右移实现：先对误差除 2，再取乘积高 8 位隐含除 2，共除 4）。数据/系数宽度 8 位，乘法结果 16 位。
// This is a generic LMS FIR filter generator
// It uses W1 bit data/coefficients bits
module fir_lms #(
parameter L = 2, // Filter length
parameter W1 = 8, // Input bit width
parameter W2 = 16 // Multiplier bit width 2*W1
)(
input wire clk, // System clock
input wire reset, // Asynchronous reset
input wire signed [W1-1:0] x_in, // System input
input wire signed [W1-1:0] d_in, // Reference input
output wire signed [W1-1:0] f0_out, // 1. filter coefficient
output wire signed [W1-1:0] f1_out, // 2. filter coefficient
output wire signed [W2-1:0] y_out, // System output
output wire signed [W2-1:0] e_out // Error signal
);
reg signed [W1-1:0] d;
reg signed [W1-1:0] x [0:L-1]; // Data array
reg signed [W1-1:0] f [0:L-1]; // Coefficient array
reg signed [W2-1:0] p [0:L-1]; // Product array
reg signed [W2-1:0] xemu [0:L-1]; // Update product array
wire signed [W2-1:0] e, sxtd, sxty, y;
wire signed [W1-1:0] emu;
// 16 bit sign extension for input d
assign sxtd = {{(W2-W1){d_in[W1-1]}}, d_in};
// Store data or coefficients
integer i;
always @(posedge clk or posedge reset) begin
if (reset) begin
d <= 0;
x[0] <= 0; x[1] <= 0;
f[0] <= 0; f[1] <= 0;
end else begin
d <= d_in;
x[0] <= x_in;
x[1] <= x[0];
// take upper half: implicit divide by 2
f[0] <= f[0] + xemu[0][W2-1:W1];
f[1] <= f[1] + xemu[1][W2-1:W1];
end
end
// Multiply, add, and error computation (1 pipeline stage via Mul block)
always @(posedge clk or posedge reset) begin
if (reset) begin
p[0] <= 0; p[1] <= 0;
xemu[0] <= 0; xemu[1] <= 0;
y <= 0; e <= 0;
end else begin
p[0] <= f[0] * x[0];
p[1] <= f[1] * x[1];
xemu[0] <= emu * x[0];
xemu[1] <= emu * x[1];
y <= p[0] + p[1]; // Compute ADF output
e <= sxtd - sxty; // error: d - y
end
end
assign emu = e[W1:1]; // e*mu: divide by 2 here and
// divide by 2 in update => mu=1/4
// Scale y by 128 because x is fraction
assign sxty = {{(W2-W1+1){y[W2-1]}}, y[W2-1:W1-1]};
assign y_out = sxty; // Monitor some test signals
assign e_out = e;
assign f0_out = f[0];
assign f1_out = f[1];
endmodule这一设计就是图 8-20 结构的直接实现：抽头延迟线的每个抽头输出乘以对应系数，结果相加得到 ；误差 经 缩放后与 相乘，用于更新系数。仿真结果如图 8-21 所示：滤波器在大约 20 步（1 μs 处）收敛到最优值 、。本设计消耗 50 个 LE 和 4 个嵌入式乘法器，TimeQuest 缓慢 85C 模型下时序性能 。
图 8-21 采用 LMS 算法的电力线干扰消除的仿真结果
这个例子同时暴露了标准 LMS 的短板：系数更新路径要在同一个时钟周期内完成两次乘法和若干次加法，组合逻辑太长， 上不去。下一节用流水线解决。
8.5.3 流水线 LMS 滤波器
从图 8-20 可以看出，原始 LMS 的更新路径很长，8 位数据和系数时性能就已相当低，因此人们提出了许多提高吞吐量的方案。先算流水线级数的最优值：使用嵌入式乘法器（自带一级输出寄存器，见第 2 章图 2-15），加法器树需要额外 级流水线，误差计算再加一级，因此最大吞吐量对应：
这里假定 是 2 的幂常数，缩放不需要额外流水线级；若采用归一化 LMS， 不再是常数，还需要按位宽增加若干级。
难点在于：LMS 含反馈（类似 IIR），不能像 FIR 那样随意插入寄存器——必须保证流水线化的滤波器仍收敛到与非流水线版本相同的系数。主流思路有四类：
- 延迟的 LMS 算法（DLMS）
- 流水线 LMS 的先行变换
- 转置形式的 LMS 滤波器
- 采用 FFT 的模块变换（模块变换已在 8.4.2 讨论过）
1. 延迟的 LMS 算法
DLMS（Delayed LMS）的核心假设是：把系数更新延迟几个采样周期后，误差梯度变化不大，即 。文献已证明：只要延迟 小于滤波器长度 ，这一近似成立，且基本不降低收敛速度。Long 最初的 DLMS 只流水线化加法器树，假定乘法与系数更新可在一个时钟周期内完成（适合 DSP 处理器）；而 FPGA 实现中乘法器和系数更新本身也需要流水线。设在滤波计算路径引入延迟 、系数更新路径引入延迟 ，则 LMS 算法 8.2 变为：
2. 先行 DLMS 算法
对 的长滤波器，延迟更新对收敛影响不大；但短滤波器就必须消除这个影响。参照第 4 章 IIR 流水线的时域交叉思想，可以只在系数计算中做一次先行（look-ahead）。从 DLMS 更新方程出发，只在系数计算中流水线化：
而真正的 LMS 误差应是：
按 Poltmann 的构想计算修正项 ，用 抵消 DLMS 引入的误差变化：
括号中的差可以递归展开：
于是修正后的误差为：
代价是修正项需要额外 次乘法，开销可能过大，因此有文献建议”放宽”修正项的要求，但仍需额外乘法器。
另一条路是用先行原则直接消除系数更新延迟（式 8-50）：
式中的求和相当于对最近 个梯度值取移动均值，使更新更平滑。它的优点是不需要额外乘法——移动均值甚至可以用一阶 CIC 滤波器（见第 5 章图 5-15）实现，运算量降到一加一减。类似 Poltmann 构想的改进方法还可见文献[294~296]。
8.5.4 转置形式的 LMS 滤波器
先行计算能平滑 DLMS 的延迟效应，但通常有额外成本。更巧的办法是直接把 FIR 结构从直接型换成转置型（见第 3 章图 3-3）：转置结构里加法器天然分布在寄存器之间，可以完全消除加法器树的延迟，式(8-49) 的流水线级数要求随之降低 级。对线性时不变系统，直接型与转置型由同一卷积方程描述；但系数时变时，需要把第 个系数从 改记为 。于是梯度估计（式 8-17）变成：
如果假定系数更新得相当慢，即 （式 8-53），梯度就近似回到熟悉的形式：
系数更新方程则为：
Jones 研究了转置型自适应滤波器的学习特性：与原始 LMS 相比，收敛速度稍慢，且稳定性要求更小的 。也就是说，转置型用一点收敛性能换来了更高的吞吐量。
8.5.5 DLMS 算法的设计
现在把例 8.5 的 LMS 滤波器流水线化。按式(8-49)，采用嵌入式乘法器时最优流水线级数为：
（若用 位软乘法器搭流水线 LE，还需要额外 级，见文献[34]。）图 8-22 显示了 8 位精度、延迟为 4 的 DLMS 的 MATLAB 仿真：与原始 LMS 相比，自适应过程中系数出现明显”过摆”，但最终仍收敛到正确值。
图 8-22(a) DLMS（延迟 4、8 位精度）的误差输出 e[n]
图 8-22(b) DLMS（延迟 4、8 位精度）的滤波器系数轨迹
例 8.6：2 抽头流水线自适应 LMS FIR 滤波器
原书 VHDL 改写为 Verilog。该设计两个系数 、，步长 ，每个乘法器及两个加法器各带 1 级流水线延迟（总延迟 4）。
// This is a generic DLMS FIR filter generator
// It uses W1 bit data/coefficients bits
module fir4dlms #(
parameter L = 2, // Filter taps
parameter W1 = 8, // Input bit width
parameter W2 = 16 // Multiplier bit width 2*W1
)(
input wire clk, // System clock
input wire reset, // Asynchronous reset
input wire signed [W1-1:0] x_in, // System input
input wire signed [W1-1:0] d_in, // Reference input
output wire signed [W1-1:0] f0_out, // 1. filter coefficient
output wire signed [W1-1:0] f1_out, // 2. filter coefficient
output wire signed [W2-1:0] y_out, // System output
output wire signed [W2-1:0] e_out // Error signal
);
reg signed [W1-1:0] f [0:L-1]; // Coefficient array
reg signed [W1-1:0] x [0:4]; // Data array (extra pipeline delay)
reg signed [W1-1:0] d [0:2]; // Reference array (extra delay)
reg signed [W2-1:0] p [0:L-1]; // Product array
reg signed [W2-1:0] xemu [0:L-1]; // Update product array
reg signed [W2-1:0] y, e;
wire signed [W2-1:0] sxtd, sxty;
wire signed [W1-1:0] emu;
wire signed [W2-1:0] y_add, e_add;
// make d a 16 bit number
assign sxtd = {{(W2-W1){d[2][W1-1]}}, d[2]};
// Store these data or coefficients in registers
integer k, i;
always @(posedge clk or posedge reset) begin
if (reset) begin
for (k = 0; k <= 2; k = k + 1) d[k] <= 0;
for (k = 0; k <= 4; k = k + 1) x[k] <= 0;
for (k = 0; k <= 1; k = k + 1) f[k] <= 0;
end else begin
d[0] <= d_in; // Shift register for desired data
d[1] <= d[0];
d[2] <= d[1];
x[0] <= x_in; // Shift register for data
x[1] <= x[0];
x[2] <= x[1];
x[3] <= x[2];
x[4] <= x[3];
// take upper half: implicit divide by 2
f[0] <= f[0] + xemu[0][W2-1:W1];
f[1] <= f[1] + xemu[1][W2-1:W1];
end
end
assign emu = e[W1:1]; // e*mu: divide by 2 here and
// divide by 2 in update => mu=1/4
// Store products in registers: 1 pipeline stage for mul & add
always @(posedge clk or posedge reset) begin
if (reset) begin
for (i = 0; i < L; i = i + 1) begin
p[i] <= 0;
xemu[i] <= 0;
end
y <= 0; e <= 0;
end else begin
for (i = 0; i < L; i = i + 1) begin
p[i] <= f[i] * x[i];
xemu[i] <= emu * x[i+3];
end
y <= y_add; // Compute ADF output: log(L) adds
e <= sxtd - sxty; // e*mu divide by 2 and 2
end
end
assign y_add = p[0] + p[1];
// Scale y by 128 because x is fraction
assign sxty = {{(W2-W1+1){y[W2-1]}}, y[W2-1:W1-1]};
assign y_out = sxty; // Monitor some test signals
assign e_out = e;
assign f0_out = f[0];
assign f1_out = f[1];
endmodule注意 x 与 d 的移位寄存器比滤波本身多出几级延迟，用于让误差与”当时参与滤波的”数据对齐——这正是 DLMS 的本质。仿真结果如图 8-23 所示：滤波器约 30 步（2 μs 处）收敛到 、，但过程中有”过摆”。本设计消耗 106 个 LE、4 个嵌入式乘法器，（TimeQuest 缓慢 85C 模型）。
图 8-23 采用延迟为 4 的 DLMS 算法的电力线干扰消除的仿真结果
对比例 8.5：流水线 LMS 的速度提升达 4 倍（69.26 MHz → 261.57 MHz），代价是 LE 数量翻倍（50 → 106），外加基准信号 与输入 的附加延迟。为了保证稳定性， 的选择必须考虑总流水线延迟的限制。
8.5.6 应用 Signum 函数的 LMS 设计
即使流水线化了，短滤波器的 LMS 成本依然很高——瓶颈在通用乘法器数量。滤波部分无法削减，能下手的是系数更新路径。由于稳定性原因，实际选的 通常远小于 ，因此可以用”符号”近似代替全精度值，有三种简化方案：
- 只用输入 的符号（符号数据函数）
- 只用误差 的符号（符号误差函数）
- 两者符号都用（符号-符号函数）
三种方案的系数更新方程统一写成：
图 8-24 的仿真结果很有启发：符号数据函数的结果与全精度几乎一致——这并不奇怪，因为本例输入 本来就近似二值的，符号运算没有引入过多量化损失；而符号误差函数完全不同：对误差的符号化实质上改变了系统的时间常数，虽然从输出 看长期存在波动，但约 2.5 s 后最终达到正确值；符号-符号算法收敛比符号误差更快，但输出同样有波动。结论是：符号简化能省掉 个更新乘法器，但必须针对具体应用仔细评估稳定性与时间常数——对某些信号，全精度 LMS 收敛而符号算法不收敛。此外还要确保整型量化本身不改变期望的系统属性。
图 8-24(a) 符号数据算法：误差输出 e[n]
图 8-24(b) 符号数据算法：滤波器系数
图 8-24(c) 符号误差算法：误差输出 e[n]
图 8-24(d) 符号误差算法：滤波器系数
图 8-24(e) 符号-符号算法：误差输出 e[n]
图 8-24(f) 符号-符号算法：滤波器系数
符号函数还要考虑一个问题：误差底线。
例 8.7：Signum LMS 滤波器中的误差底线
在 8.3.1 的系统辨识结构中使用符号类 ADF 算法，能达到什么误差底线？直观预期：符号化损失精度，误差底线应高于全精度 LMS，学习速度也稍慢。对 50 次迭代均值、两种 EVR 的仿真（图 8-25）证实了这一点：符号数据算法虽然自适应稍慢，但仍达到了 -60dB 的误差底线（与全精度相同的最大误差）；符号误差算法与符号-符号算法延迟更大，且只能达到约 -40dB 的误差底线。这样的误差对一些应用可以接受，对另一些则不行。
图 8-25(a) 有符号数据算法在低 EVR 下的学习曲线
图 8-25(b) 有符号数据算法在高 EVR 下的学习曲线
图 8-25(c) 有符号误差算法的学习曲线
图 8-25(d) 符号-符号 LMS 算法的学习曲线
符号-符号算法无论在软件还是硬件实现上都颇具吸引力，已被用作自适应差分脉冲编码调制（ADPCM）的 ITU（国际电信同盟）标准。不过在硬件中其实不必真正实现”符号×符号”：与 相乘只是一个常数缩放，而且任选一种单符号方案就已经省掉了图 8-20 结构中系数更新所需的 个乘法器。
8.6 递归最小二乘法算法
LMS 的哲学是”随机梯度”：系数一步一步地朝维纳-霍夫最优解挪动。递归最小二乘（Recursive Least Square, RLS）走另一条路：每一对新数据 到来时，直接迭代更新自相关矩阵 和互相关向量 的估计，然后重解维纳-霍夫方程 。每个新输入 放入长度 的数据阵列 中： 加到 ， 加到 ，以此类推。数学上只需计算叉积 并累加。递归方程为：
再用时间递归形式的维纳-霍夫方程：
严格地说，相关量应按求和次数 缩放才是真估计，但两者缩放倍数相同、在算法中互相抵消，于是系数更新就是：
这种”强制版”RLS 计算量巨大（矩阵求逆约需 次运算），但揭示了基本思想，编程很快。例如 MATLAB 中长为 的 RLS 内层循环：
x = [xin;x(1:L-1)]; % get new sample
y = f' * x; % filter output
err = din - y; % error: reference - filter output
Rxx = Rxx + x*x'; % update the autocorrelation matrix
rdx = rdx + din .* x; % update the cross-correlation vector
f = Rxx^(-1) * rdx; % compute filter coefficients其中 Rxx 是 矩阵，rdx 是 向量，通常初始化 。还有一个初始化陷阱：最初 次迭代时 只有少数非零项，是奇异矩阵，没有逆。三种应对办法：
- 等到 再求逆。
- 用伪逆 （求解超定线性方程组的标准结果）。
- 用 初始化 ，高信噪比取小 、低信噪比取大 。
第三种最常用：计算方便，还能用 控制”初始学习速度”。图 8-26 展示了初始化的影响（上半为 4000 次迭代全程，下半为最初 100 次迭代）：高信噪比（-48dB）时选较大的初始化值能得到快速收敛；低信噪比（-10dB）时要用较小初始化值，否则初始迭代会出现大误差——是否可接受取决于具体应用。
图 8-26(a) 高信噪比 -48dB、δ=1000 的学习曲线
图 8-26(b) 高信噪比 -48dB、δ=1 的学习曲线
图 8-26(c) 低信噪比 -10dB、δ=1 的学习曲线
图 8-26(d) 低信噪比 -10dB、δ=1/1000 的学习曲线
更吸引人的方法根本不求 本身，而是直接对逆矩阵 做时间递归。把维纳方程取 时刻并代入式(8-58)：
再由式(8-57)把 换成 ：
两边乘 并整理，把 移到左边：
其中先验误差定义为：
**卡尔曼增益（Kalman gain）**向量定义为：
式(8-62)的形状与 LMS 更新一致，但”步长×输入”换成了 ——一个随数据自动调节的增益向量。剩下的任务是不经矩阵求逆地递归更新 ，这里用矩阵求逆引理：
对所有维数一致（且 非奇异）的矩阵成立。代入 、、、，得到：
借助卡尔曼增益 可以写得更简洁：
初始化取 ，其中 对高信噪比取大的正常数、对低信噪比取小的正常数。这一递归的计算量与 成正比——相比 的直接求逆，对大 是根本性的降低。图 8-27 给出了 RLS 干扰消除的基本结构小结。
图 8-27 采用 RLS 算法进行干扰消除的基本结构
8.6.1 有限记忆的 RLS 算法
由式(8-57)(8-58)可以看出，到目前为止的 RLS 具有无限记忆：系数是零时刻以来所有输入的函数。实践中常引入遗忘因子，让近期数据更重要，同时保证更新过程不会溢出。做法是把代价函数从平方和换成指数加权平方和：
决定有效记忆： 即无限记忆； 时有效记忆长度为 个数据点。
算法 8.8：RLS 算法
按如下步骤调整自适应滤波器的 个系数：
(1) 初始化 和 。
(2) 接受新输入对 ，并把输入信号移入基准信号向量 。
(3) 计算 FIR 滤波器输出：
(4) 计算先验误差：
(5) 计算卡尔曼增益因子：
(6) 更新滤波器系数：
(7) 更新逆自相关矩阵：
之后回到步骤(2)。
RLS 每个输入采样需要 次乘法和 次加/减法，计算量远高于 LMS。但它的好处在下面的例子中一目了然：收敛快，而且根本不需要选步长 ——而选择 恰恰是保证 LMS 稳定时的难点。
例 8.9：RLS 学习曲线
用系统辨识结构比较 RLS 与 LMS 的收敛性（配置同图 8-12）：自适应滤波器长度 ，与”未知”系统同长；未知系统的附加噪声 -48dB（相当于 8 位量化）。LMS 的收敛速度取决于 EVR（见式 8-28），为生成不同 EVR，把 的高斯白噪声通过表 8-2 的 FIR 滤波器，并将系数标准化为 以保持信号功率不变。“未知系统”的脉冲响应是奇对称滤波器，系数为 （图 8-28(a)）。LMS 的最大步长为：
为保稳定取 。DCT-LMS 对每个系数采用功率标准化（见图 8-17）。结果（图 8-29）表明：随 EVR 增大，RLS 收敛明显快于 LMS；DCT-LMS 也比 LMS 快，某些情形与 RLS 一样快，但在残差水平和收敛一致性上不如 RLS——EVR=1 时 DCT-LMS 达到约 50dB 的水平，EVR=100 时只有 40dB；而 RLS 在各情形下都收敛到系统噪声底线。
图 8-28(a) “未知系统”的脉冲响应
图 8-28(b) EVR=100 的 RLS 系数学习曲线
图 8-29(a) EVR=1 时 LMS、DCT-LMS 与 RLS 的学习曲线
图 8-29(b) EVR=10 时 LMS、DCT-LMS 与 RLS 的学习曲线
图 8-29(c) EVR=100 时 LMS、DCT-LMS 与 RLS 的学习曲线
图 8-29(d) EVR=1000 时 LMS、DCT-LMS 与 RLS 的学习曲线
8.6.2 快速 RLS 算法的卡尔曼实现
前面接触到的 RLS 类算法都有一个共同的麻烦：为了更新滤波器系数，需要维护并求逆一个 的自相关矩阵，计算量是 。当滤波器长度 达到几十、上百时，这个代价在 FPGA 上非常可观。Ljung 等人提出的最小二乘 FIR 快速卡尔曼（Kalman）算法换了一条思路：把”一步前向预测”和”一步后向预测”这两个概念推到舞台中央。在所有递归的一维 Levinson-Durbin 类型算法中使用前向和后向预测系数之后，卡尔曼增益向量的更新只需要 次运算——这正是”快速”二字的来历。
一步前向预测器
所谓 阶一步前向预测器，就是用当前时刻之前的 个最近采样值 去估计当前值 ，如图 8-30 所示。
图 8-30 阶线性前向预测
预测得好不好，用后验误差来量化：
上标 表示”前向”（forward）预测误差，下标 表示预测器的阶数（也就是系数个数）。本节其余部分略去下标 ：如果没有特别说明，向量的长度就是 。同样也可以定义先验误差——它使用上一次迭代（ 时刻）的滤波器系数来计算：
为什么要同时保留两种误差？因为后验误差用的是”刚更新完”的系数，反映的是此刻预测器的真实水平；而先验误差在递推链里出现得更早，是系数更新公式中实际可用的量。快速卡尔曼算法正是在这两种误差之间来回衔接。
要使 的最小二乘误差达到最小，令误差平方对系数的偏导为零：
这又导出了一个 自相关矩阵方程，但右侧与维纳–霍夫方程不同：
成本函数的最小值（也就是最小的总平方误差）为：
其中 。
这里的关键依据是 Levinson-Durbin 算法：它能以递归方式求解式(8-72)的最小二乘误差最小值，而不需要显式计算矩阵逆。为了更新预测器系数，要用与式(8-69)相同的卡尔曼增益因子：
利用分块矩阵求逆扩展增益向量
接下来的问题是：卡尔曼增益 本身如何迭代更新？观察一个事实：相邻两次迭代之间，增益向量只有第一个和最后一个元素不同。基于这一点，可以借助长度为 的下一个卡尔曼增益来更新，由式(8-68)得：
为了计算逆矩阵 ，可以利用分块矩阵求逆的一个著名定理：如果 是非奇异矩阵，那么
初学者不必死记这个庞大的公式，只需理解它的作用：把一个 矩阵的求逆，化简成若干个”小块乘法”加上一个已经知道的 逆矩阵 的复用。做如下关联：
于是得到中间结果：
请注意这些结果的形式： 恰好就是新的前向预测系数 ，而那个标量正好是最小的总平方前向误差 ——分块求逆自动”变出”了我们已经熟悉的量。把中间结果代回式(8-77)，可将 改写为：
重新整理式(8-77)，就得到用前向预测系数表达的增益更新：
到这里还没有得到闭合的递归—— 和 还只是”半成品”。
一步后向预测器
为什么只靠前向预测不够？仔细看上式：新向量由 （前 个元素）和 （最后一个元素）组成，而 需要从 和 中”解”出来，这个解需要用到后向预测系数。因此，卡尔曼增益向量的迭代更新除了前向预测系数外，还需要一步后向预测器的系数。
一步后向预测器反过来：用当前的 个采样 去估计 步以前的值 。其后验误差函数是：
所有向量的大小都是 。线性后向预测器的结构如图 8-31 所示。
图 8-31 阶线性后向预测
后向预测器的先验误差是：
后向预测器最小二乘系数的求解方程与前向情况完全对称，由下式给出：
总平方误差的最小值变成：
其中 。更新后向预测器系数同样要借助式(8-69)中的卡尔曼增益因子：
现在再次为扩展的卡尔曼增益向量寻找 Levinson-Durbin 类型的递归方程，只是这次使用后向预测系数：
为了求解矩阵的逆，需要像式(8-79)一样定义 分块矩阵 ，只是这次的分块方式不同：块 需要是非奇异矩阵，于是有另一种形式：
进行如下关联：
由此得到如下中间结果：
可以看到，这次分块求逆自动产生了后向预测系数和最小总平方后向误差。利用这些中间结果就得到：
用后向预测系数重新整理式(8-83)，写成：
两条路（前向系数版、后向系数版）都把 拆成了 和 ，把两个表达式联立，就能解出 ——这正是算法8.10中倒数第二个方程的来历。
补上误差更新方程
目前唯一缺少的是总平方误差最小值 、 的迭代更新方程，由下式给出：
直观理解：每来一个新采样，最小总平方误差要么增加、要么减少一个小的修正量，修正量是后验误差与先验误差的乘积——当预测器已收敛时两种误差几乎相等且趋近于零， 就基本保持稳定。
算法8.10 快速卡尔曼 RLS 算法
预开窗的快速卡尔曼 RLS 算法按如下步骤调整自适应滤波器的 个滤波器系数：
(1) 初始化 ，且 （ 是一个很小的正数，防止初始时除以零）。
(2) 接收新的输入采样对 。
(3) 按顺序更新 、 和 ：
(4) 将输入信号 在基准信号向量 中移位，并通过如下两个方程更新自适应滤波器系数：
接下来重复步骤(2)。
统计一下计算量就会发现：步骤(3)需要两次除法、 次乘法和 次加/减法；步骤(4)中的系数更新还需要额外的 次乘法和加/减法。所以总计算量是 次乘法、 次加/减法和两次除法。对比直接形式的 RLS 的 次乘法，当 较大时（比如 ： 对 ），“快速”名副其实。
8.6.3 快速后验卡尔曼 RLS 算法
仔细观察算法 8.10 就会发现，Ljung 等人提出的原始快速卡尔曼算法主要以先验误差方程为基础。Carayannis 等人提出的快速后验误差时序技术（Fast A posteriori Error Sequential Technique, FAEST）更大程度上采用了后验误差。这一算法探讨了快速卡尔曼算法中不同参数的更多迭代性质，把计算量又降低了 次乘法。原始快速卡尔曼和 FAEST 采用大体相同的思想：将卡尔曼增益的长度增加 1，以及使用前向预测 和后向预测 。FAEST 还引入了遗忘因子 （，用于削弱旧数据的影响）。下述代码清单给出了 MATLAB 中 FAEST 算法的内层循环：
%********** FAEST Update of k, a, and b
ef=xin - a'*x; % a priori forward prediction error
ediva=ef/(rho*af); % a priori forward error/minimal error
ke(1)=-ediva; % extended Kalman gain vector update
ke(2:1+1)=k - ediva*a;% split the l+1 length vector
epsf=ef*psi; % a posteriori forward error
a=a+epsf*k; % update forward coefficients
k=ke(1:1) + ke(1+1).*b; % Kalman gain vector update
eb=-rho*alphab*ke(11); % a priori backward error
alphaf=rho*alphaf+ef*epsf; % forward minimal error
alpha=alpha+ke(1+1)*eb+ediva*ef; % prediction crosspower
psi=1.0/alpha; % psi makes it a 2 div algorithm
epsb=eb*psi; % a posteriori backward error update
alphab=rho*alphab+eb*epsb; % minimum backward error
b=b-k*epsb; % update backward prediction coefficients
x=[xin;x(1:1-1)]; % shift new value into filter taps
%********** Time updating of the LS FIR filter
e=din-f'*x; % error: reference - filter output
eps=-e*psi; % a posteriori error of adaptive filter
f=f+w*eps; % coefficient update从代码可以看出 FAEST 的核心技巧：变量 psi（即 ）被保存下来重复使用，使每次迭代只需要两次除法；先验量与后验量通过 psi 相互转换，链条环环相扣。总计算量（不考虑指数权重 的乘法）是两次除法、 次乘法和 次加/减法。
8.7 LMS 和 RLS 的参数比较
最后用表 8-3 把本章介绍的算法放在一起比较。表中根据计算复杂性对基本随机梯度（Stochastic Gradient, SG）方法（如有符号 LMS（SLMS）、标准化 LMS（NLMS）或使用 FFT 的分块 LMS（BLMS）算法）进行了比较；接下来是变换域算法（TDLMS），但不包括功率标准化，也就是变换域中 的标准化；RLS 一族则包括直接形式、（快速）卡尔曼算法、网格算法和 FAEST 算法。通常情况下，网格算法（本章没有讨论）需要大量的除法和求平方根运算，此时建议使用对数系统（请参阅第 2 章）[299]。
表 8-3 长度 的自适应滤波器的 LMS 和 RLS 算法的复杂性比较。其中 TDLMS 没有标准化。如果在 TDLMS 算法中使用标准化，就要增加 次乘法、 次加/减法和 次除法
| 算法 | 实现方式 | 计算量 | ||
| 乘法 | 加/减法 | 除法 | ||
| SG | LMS | $2L$ | $2L$ | - |
| SLMS | $L$ | $2L$ | - | |
| NLMS | $2L+1$ | $2L+2$ | 1 | |
| BLMS(FFT) | $10\log_{2}(L)+8$ | $15\log_{2}(L)+30$ | - | |
| SG TDLMS | 阿达马 | $2L$ | $4L-2$ | - |
| 哈尔 | $2L$ | $2L+2\log_{2}(L)$ | - | |
| DCT | $2L+\frac{3}{2}L\log_{2}(L)+L$ | $2L+\frac{3}{2}L\log_{2}(L)$ | - | |
| DFT | $2L+\frac{3}{2}L\log_{2}(L)$ | $2L+\frac{3}{2}L\log_{2}(L)$ | - | |
| KLT | $2L+L^{2}+L$ | $2L+2L$ | - | |
| RLS | 直接形式 | $2L^{2}+4L$ | $2L^{2}+2L-2$ | 2 |
| 快速卡尔曼 | $10L+2$ | $9L+2$ | 2 | |
| 网格 | $8L$ | $8L$ | 6L | |
| FAEST | $7L+8$ | $7L+4$ | 2 | |
表 8-3 中的数据基于第 6 章讨论过的 DCT 和 DFT 及其基于快速 DIF 或 DIT 算法的实现方式。对于长度为 8 和 16 的 DCT 或 DFT，已经开发了使用更少次运算的更有效的（Winograd 类型）算法。例如，长度为 8 的 DCT（请参阅图 6-23）使用了 12 次乘法，这样 DCT 变换域算法可以由 次乘法实现，而 FAEST 算法是 次乘法。但要注意：这些计算没有考虑功率标准化，而功率标准化是所有 TDLMS 都必须遵循的（否则与标准 LMS 相比就没有快速收敛性 [287, 288]）。
在硬件实现层面还要记住一点：除法比乘法贵得多。当功率标准化因子可以预先确定时，就可以用硬连线的缩放运算实现除法。FAEST 只需要两次除法，与 ADF 长度无关，这是它对硬件友好的重要原因。
从收敛性能看：例 8.9 已经比较过 RLS 和 LMS 的自适应速度——RLS 类型算法的适应速度比 LMS 快很多，但 LMS 算法可以用变换域算法（如 DCT-LMS）从根本上进行改进。而且通常与 LMS 或 TDLMS 相比，RLS 算法的误差底线和误差一致性更佳。代价是：没有除法操作就不能实现 RLS 类型的算法，而这经常需要更大的整体系统位宽，至少是分数表示方法，甚或是浮点数表示方法 [299]。另一方面，LMS 算法仅用例 8.5 中的几位就可以实现。简单说：LMS 省资源但收敛慢，RLS 收敛快但吃算力和位宽，选型时要在两者之间权衡。
8.8 主成分分析（PCA）
迄今为止讨论的信号处理方法都是单输入的，并且绝大部分模型信号都经过精心定义。接下来讨论更复杂的问题：多输入多输出，而且精确的信号模型不可知。一个典型例子是测量和评估用于心脏心电活动测量的心电图（Electrocardiograph, ECG）。一般从身体的不同部位进行 3 到 12 次测量：在手和腿之间有 3 个双极和 3 个单极导联，围绕心脏还有 6 个。健康的 ECG 信号 [300] 由规则的脉冲构成，每个脉冲由前波 P（大约 80ms）、主峰 QRS（80ms～120ms）和后波 T（160ms）混合而成。典型的幅度是 、、、 以及 ，请参阅图 8-32。
图 8-32 ECG 示例信号
真正的难点在于：从体表测到的是高达 12 路信号的混合，要把其中真实的单个信号（如果要检测胎儿或双胞胎心率，则是 2 或 3 个 [301]）分离出来。对于胎心检测，主要工作是确保胎心率保持在每分钟 110 次（beat per minute, bpm）到 180 次的期望区间内。
选择的工具是 Kahunen-Loève 变换（尽管相当古老），它又称为主成分分析（Principle Component Analysis, PCA）。思路分三步：在输入信号之间构建互相关矩阵，计算特征值和特征向量，然后使用特征值和特征向量转换输入信号。写成公式的形式为：
其中 是输入的均值向量。设 为第 个特征向量， 为 对应的特征值，则第 个主成分由下式得到：
特征值的分布可以很好地说明输入数据中有”多少”个真正的成分：能量大的信号对应大的特征值，噪声对应接近零的小特征值，通常把特征值较小的分量剔除掉。最终，从最小均方意义上得到一个或多个代表输入信号的信号。主成分捕捉了绝大部分信号能量，同时所有输出信号彼此正交——这也是 PCA 常常被称作”白化”操作的原因。下面构建一个胎心监测 ECG 模型来实践它。
例8.11 ECG 数据降维
使用一个简单的 ECG 模型，其中妈妈和胎儿的 ECG 以脉冲序列的形式建模。采样频率设置为 1kHz。妈妈的 ECG 周期长度为 14（约 71bpm），胎儿的周期长度为 8（125bpm）。假设采用 4 通道测量，得到两个信号的叠加，并且信号的极性（正负）也应该确定下来。输入的混合信号以随机混合矩阵的形式建模，并假设母婴 ECG 幅值之间的比例是 10:1。下述 MATLAB 脚本给出了完整的监测系统的实现：
MatLab script ECG PCA example
clear all; close all;
T1=14; % Mother period length at 1 kHz
T2=8; % Fetus period length at 1 kHz
T=T1*T2*2; % Total samples in simulation
s=zeros(2,T);
s(1,:)=upsample(ones(1,T/T1),T1,T1/2); % Mother ECG
s(2,:)=upsample(ones(1,T/T2),T2,T2/2)/10; % Fetus ECG
subplot(3,1,1);L=50;t=1:L; % Plot signals without noise
plot(t,s(1,1:L),'k-x',t,s(2,1:L),'k-o');
title('Undisturbed mother (x) and fetus (o) signals');
axis([0 L -.2 1.2]);ylabel('s[n]')
rng(0); % Initialize random number to start
A=[-1 -1;-1 1;1 -1;1 1];
A=A+(rand(4,2)-.5); % Mixing matrix with polarity change
x=A*s+rand(4,T)/100; % Apply mixing and add Gauss noise
subplot(3,1,2); % Plot measured signals with noise
plot(t,x(1,1:L),'k-x',t,x(2,1:L),'k-o',...
t,x(3,1:L),'k-+',t,x(4,1:L),'k-*');
title('The 4 measured ECG signals (with noise)')
axis([0 L -1 2]);ylabel('x[n]')
Cxx = cov(x', 1) % Calculate the covariance matrix
disp('Eigenvectors and eigenvalues of covariance matrix:')
[E, D] = eig(Cxx)
%apply reconstruction mixing to the x data i.e. do PCA
y=(D^- .5 * E'*x);
subplot(3,1,3);
plot(t,y(4,1:L),'k-x',t,y(3,1:L),'k-o');
title('Mother (x) and fetus (o) signals after PCA')
grid
axis([0 L -.5 4.5]);xlabel('Sample n');ylabel('y[n]')
print -deps ecg1.eps;print -djpeg ecg1.jpg脚本的流程是：先生成并绘制两个测试信号；然后用混合矩阵生成四个测量信号，并向信号中添加高斯噪声；随后基于互相关矩阵 执行 PCA。MATLAB 计算出的互相关矩阵为：
0.0310 0.0264 -0.0276 -0.0635
0.0264 0.0244 -0.0264 -0.0551
-0.0276 -0.0264 0.0291 0.0583
-0.0635 -0.0551 0.0583 0.1308特征值以方阵对角元素的形式给出，并由小到大排列， 为：
对应的特征向量矩阵 为：
| -0.6078 | -0.5971 | -0.3604 | -0.3798 |
| 0.5737 | -0.6215 | 0.4152 | -0.3350 |
| 0.4797 | -0.2568 | -0.7594 | 0.3567 |
| -0.2672 | -0.4374 | 0.3479 | 0.7850 |
如期望的一样，只得到两个非零特征值——也就是说，4 路测量信号里其实只藏着两个主成分：妈妈和胎儿。
图 8-33 中的仿真结果首先给出了两个输入信号，然后是四组测量信号，第三幅图给出了两个重构的 ECG 信号。注意，信号的周期（也就是 8 和 14）保持良好，只有幅值与原始数据相比拥有不同的比例因子。PCA 完成了从 4 路信号到两路信号的降维工作，同时把两个信号分离了出来。
四组 ECG 测量信号（含噪声）
PCA 后妈妈（x）和胎儿（o）的信号
图 8-33 ECG 示例的 PCA 仿真结果
8.8.1 主成分分析的计算
PCA 的理论很优雅，但在 FPGA 上落地时必须回答：特征值和特征向量到底怎么算？这部分需要大量的算术运算。互相关矩阵本身可以通过下式迭代计算：
这是一个一阶低通式的递推， 越小历史数据权重越大。根据噪声水平，上式可能需要很多个周期的循环；由于噪声较低，两个输入采样在大约 50～100 个采样周期后应该就能得到足够稳定的值。由于互相关矩阵 是对称的，计算可以简化一点：对称矩阵的特征值和特征向量都是实数，且特征向量相互正交、可以归一化为单位长度，即 。不过从本质上讲，特征值和特征向量的计算仍然比互相关矩阵的计算更复杂。根据文献 [302, 303]，主要有三种方法：
(1) 直接方法
(2) v.Mises 功率法
(3) 神经网络学习方法
这三种方法在硬件需求、算法编码、收敛循环次数以及需要考虑的具体事项等方面有本质上的不同。下面逐一研究。
直接方法。为了计算特征值和特征向量，首先计算特征多项式来确定特征值，即 ，多项式的零点就是特征值 ；为了计算特征向量，需要求解矩阵系统 。该算法不需要任何迭代，也不会出现收敛异常，只需要少数几个能够在微处理器或 MATLAB 中计算的矩阵等式。但是，即使对于少量的信号，矩阵系统也需要相当大的硬件资源 [304]——所以它适合”离线算好、硬件照搬”，不适合在线自适应。
功率法。只需要少量迭代的第二类方法归功于 v.Mises [302, 303]。该方法从一个随机开始向量 出发，反复左乘 ，计算向量序列 ：
直觉上可以这样理解：任何向量都可以分解成各特征向量的叠加，反复乘以 之后，最大特征值对应的分量被放大的倍数最多，几轮下来就”一枝独秀”了。功率方法并不需要计算互相关矩阵的高次幂，因为：
也就是说每步只需一次矩阵–向量乘法。通过恰当的初始化，向量将收敛到主特征向量 。唯一的例外是：如果 “碰巧”已经是另一个特征向量（例如 ），那么迭代只会放大它自己，永远得不到 。下面的例子展示了该方法收敛得有多快。
例8.12
使用初始向量 时，例 8.11 的 相关矩阵的 v.Mises 迭代如下：
通过 v=v./sqrt(sum(v.^2))' 把向量标准化为长度 1，得到：
与六位小数的精确特征向量一致。仅仅三次迭代后，误差的近似精度 就从 5.1 位依次增加到 10.9 和 16.8 位——大约每迭代一次，精度翻一倍。
特征值的计算可以利用相邻两次迭代对应分量的比值完成：，含 4 位有效数字。道理很简单：收敛后每乘一次 ，向量就缩放一个 的倍数。
由于特征向量相互正交，后续特征值的计算可以采取类似的方法。在 v.Mises 提出的方法中，用随机向量 计算 ，再通过 把已经找到的方向从初始向量中”减掉”，然后正常进行功率迭代。
例8.13
例 8.12 的 相关矩阵已经生成了 。由初始化向量 ，得到 ，随后迭代从向量 开始。对于第二个特征向量，算法的执行过程如下：
对于 ，位精度依次达到 5.4、14.2 和 23.0 位。真实的特征向量 最终需要通过 标准化为长度 1：
所有的 6 个小数数字都是正确的。特征值的计算可以通过 完成。
现在，功率法的优点很明显了：算法只需要少量迭代，硬件只需要实现一个”矩阵–向量乘”等式。主要不足在于向量的初始化——如果”碰巧”用较小特征值对应的特征向量进行初始化，第一次运行将无法得到最大的特征值和相关的特征向量。另一个代价是：最终向量的标准化需要开方、平方根和除法运算，成本可能很高。考虑到特征向量乘以任意常数后仍然是特征向量（任意缩放不改变方向），建议只使用向量的最大分量进行标准化，省掉开方。FPGA 上对功率法的实现通常基于高效的矩阵操作来完成，例如 QR 分解 [305]。
神经网络学习方法。第三种方法以类似于自适应滤波器 LMS 算法的方式逐步调整向量系数。基本算法使用 Hebbian 学习规则——“输入输出乘积越大，权重增加越多”，即向权重增加 ，其中 。Oja 和 Karhunen [306] 指出 将收敛到系统的特征向量。为了避免权重无限制增长，算法需要正确的标准化，由此引出下述方程（被称作 Oja 学习规则 [306]），对于第一特征向量为：
括号里的 项就是隐含的归一化：它防止 无限增长。上式也称为随机梯度上升（Stochastic Gradient Ascent, SGA）算法——LMS 做的是”下降”（减小误差），这里做的是”上升”（增大投影能量）。
后来，Sanger 在他的硕士论文中对该方法做了一点小的修正，使其能够按顺序找出多个特征向量 [307]。Sanger 的方法可以归纳为：
该方法名为通用 Hebbian 算法（Generalized Hebbian Algorithm, GHA）。关键在于求和上限是 ：学习第 个成分时，把前 个成分的贡献都减掉，相当于在”去掉了前面主成分的残差”上继续找方向，这就保证了各 收敛到不同的（相互正交的）特征向量。对于双通道系统，图 8-34 和图 8-35 分别以 SIMULINK 模型的形式给出了第一和第二特征向量的学习系统。输出 就是来自 PCA 分析的期望信号。
图 8-34 采用 Sanger 的 GHA 算法的混合矩阵以及第一主成分学习系统
图 8-35 采用 Sanger 的 GHA 算法的第二主成分学习系统
8.8.2 Sanger GHA PCA 的实现
从所需硬件资源的角度看，GHA 是三种方法中最好的 [308]，下面讨论一个完整的设计示例。
例 8.14 基于 Sanger GHA 的 PCA 设计
假设拥有两个源信号：二进制数字信号 （每位采样 5 次）和正弦波输入 （周期为 6）。这两个信号通过如下混合矩阵组合在一起，生成系统输入 和 ，如图 8-34 所示：
经过一个初始学习过程，期望 PCA 系统收敛并把原始信号分离出来。整个过程约 4K 个时钟周期达到收敛。由于需要学习两个特征向量，首先从 开始：经过 1000 个时钟周期后，学习速率 逐渐减小为零；接下来学习第二个特征向量 ，再经过 1000 个时钟周期使 减小到零，得到一个稳定的输出向量。学习速率 及特征向量值如图 8-36(a) 所示。为了验证学习结果，直接计算互相关矩阵及对应的特征向量，得到：
与 GHA 学习得到的数据比较：两个特征向量之间只有符号发生变化，数值完全一致（符号差异无关紧要，因为 与 代表同一方向）。
(a) 总体学习速率和特征向量
(b) 收敛后的输入、混合和输出信号
图 8-36 Sanger GHA 的 SIMULINK 仿真结果
GHA 的仿真结果如图 8-36(b) 所示，从采样点 4000 开始，该点处特征向量已经稳定。可以看到两个输入信号 ，随后是混合输出 及输出信号 ——输入信号被很好地还原了出来。
原书给出的 VHDL 设计代码可以等价地改写为如下可综合的 Verilog HDL 实现（内部数据格式为 16.16 有符号定点数，即 32 位中高 16 位为整数部分含符号、低 16 位为小数部分）：
// ============================================================
// Sanger GHA PCA (原书 VHDL 改写为 Verilog, 可综合风格)
// 数据格式: Q16.16 有符号定点数 (SFIXED(15 DOWNTO -16) 等价)
// ============================================================
module pca (
input wire clk, // 系统时钟
input wire reset, // 系统复位
input wire [31:0] s1_in, // 1. 信号输入 (Q16.16)
input wire [31:0] s2_in, // 2. 信号输入 (Q16.16)
input wire [31:0] mu1_in, // 学习速率 1. 主成分
input wire [31:0] mu2_in, // 学习速率 2. 主成分
output wire [31:0] x1_out, // 混合 1. 输出
output wire [31:0] x2_out, // 混合 2. 输出
output wire [31:0] w11_out, // 特征向量 [1,1]
output wire [31:0] w12_out, // 特征向量 [1,2]
output wire [31:0] w21_out, // 特征向量 [2,1]
output wire [31:0] w22_out, // 特征向量 [2,2]
output reg [31:0] y1_out, // 1. 主成分输出
output reg [31:0] y2_out // 2. 主成分输出
);
// 混合矩阵常量 a_k,i (Q16.16), 采用截断(不回绕)策略
localparam signed [31:0] A11 = 32'sd49152; // 0.75 = 3/4
localparam signed [31:0] A12 = 32'sd98304; // 1.5
localparam signed [31:0] A21 = 32'sd32768; // 0.5 = 1/2
localparam signed [31:0] A22 = 32'sd21846; // 1/3 ≈ 0.333333
localparam signed [31:0] INI = 32'sd32768; // 权重初值 0.5 (非零随机初值)
// 定点乘法: Q16.16 * Q16.16 = Q32.32, 算术右移 16 位截断回 Q16.16
function signed [31:0] mul_trunc;
input signed [31:0] a;
input signed [31:0] b;
reg signed [63:0] p;
begin
p = a * b; // 64 位乘积, 不回绕
mul_trunc = p >>> 16; // 截断低位, 取高 32 位
end
endfunction
// 输入位重定义为有符号 Q16.16 (不消耗任何逻辑资源)
wire signed [31:0] s1 = $signed(s1_in);
wire signed [31:0] s2 = $signed(s2_in);
wire signed [31:0] mu1 = $signed(mu1_in);
wire signed [31:0] mu2 = $signed(mu2_in);
reg signed [31:0] x1, x2;
reg signed [31:0] w11, w12, w21, w22;
reg signed [31:0] y1, y2, h11, h12; // 中间变量(阻塞赋值)
always @(posedge clk or posedge reset) begin
if (reset) begin // 复位/初始化所有寄存器
x1 <= 32'sd0; x2 <= 32'sd0;
w11 <= INI; w12 <= INI;
w21 <= INI; w22 <= INI;
y1_out <= 32'sd0; y2_out <= 32'sd0;
end else begin
// ---------- 混合矩阵 ----------
x1 <= mul_trunc(A11, s1) + mul_trunc(A12, s2);
x2 <= mul_trunc(A21, s1) - mul_trunc(A22, s2);
// ---------- 第一主成分及特征向量 (Oja/SGA 规则) ----------
y1 = mul_trunc(x1, w11) + mul_trunc(x2, w12); // y1 = w1'x
h11 = mul_trunc(w11, y1); // h11 = w11*y1
w11 <= w11 + mul_trunc(mul_trunc(mu1, x1 - h11), y1);
h12 = mul_trunc(w12, y1); // h12 = w12*y1
w12 <= w12 + mul_trunc(mul_trunc(mu1, x2 - h12), y1);
// ---------- 第二主成分及特征向量 (Sanger GHA 规则) ----------
y2 = mul_trunc(x1, w21) + mul_trunc(x2, w22); // y2 = w2'x
w21 <= w21 + mul_trunc(mul_trunc(mu2, x1 - h11 - mul_trunc(w21, y2)), y2);
w22 <= w22 + mul_trunc(mul_trunc(mu2, x2 - h12 - mul_trunc(w22, y2)), y2);
// ---------- 输出寄存 (便于测试时序性能) ----------
y1_out <= y1;
y2_out <= y2;
end
end
// 内部数据重新定义为 32 位 SLV 输出, 便于监控
assign x1_out = x1;
assign x2_out = x2;
assign w11_out = w11;
assign w12_out = w12;
assign w21_out = w21;
assign w22_out = w22;
endmodule对于第一主成分，上述设计采用了图 8-34 的方案（标准 Oja/SGA 规则）；对于第二主成分，采用了图 8-35 的方案（Sanger 修正：减去前一个成分的贡献 ）。设计上值得注意的是代码的组织方式：
- 指定 I/O 之后，将混合矩阵 定义为常值（
localparam）； - 输入信号以 Q16.16 有符号格式重新定义，该过程不需要任何资源，只是把同样的 32 位比特重新”看待”为有符号数；
- 单个时钟
always块执行整个 PCA 算法：复位时寄存器 清零，四个权重寄存器 初始化为非零值（此处选 0.5）——权重初值不能为零，否则 ，学习永远无法启动； - 复位语句之后，先执行混合矩阵，随后的代码依次计算第一主成分和第二主成分；
always块还包括对 的输出寄存，以便测试设计的时序性能；- 最后，为了便于监控，把一些内部数据（、）分配给输出引脚。
该设计用到了 2447 个 LE、180 个嵌入式乘法器，使用 TimeQuest 缓慢 85C 模型的时序电路性能 。可以看到，即便是 这么小的 PCA，嵌入式乘法器的消耗也不低——这正是学习型算法的硬件代价。
图 8-36 的 GHA 仿真结果由 HDL 进行了二次生成。首先，通过图 8-37 可以看到学习速率和特征向量的收敛轨迹；接下来，图 8-38 给出了输入信号和复原的主成分。从这一点上，仿真过程中得到了稳定的特征向量。
图 8-37 采用 HDL 的 Sanger GHA 的整体仿真结果：学习速率和特征向量
图 8-38 采用 HDL 的 Sanger GHA 的整体仿真结果：输入信号以及收敛后的输出主成分
最后要指出 GHA 的定位：它是迄今为止我们讨论过的最慢的 PCA 算法——每个时钟节拍只做很小一步梯度调整，收敛需要数千个周期。有没有更快的？遗憾的是，诸如 APEX、PAST 或 OPAST [309-311] 等最快的算法都是非 PCA 的子空间方法（Subspace Method）：它们不计算特征向量本身，而是计算一组旋转基向量。比如 LMS 算法，绝大多数此类算法都使用一种 RLS 类型替换缓慢的学习速率，但总体上需要相当数量的矩阵运算和除法操作，从而降低了硬件的时序性能。换句话说：收敛速度与硬件友好度在 PCA 领域同样是一对需要权衡的矛盾——这与本章开头 LMS 与 RLS 的取舍形成了漂亮的呼应。
8.9 独立成分分析(ICA)
上一节的 PCA 能把输入数据空间的维数降下来，但如果想让 PCA 直接把混合信号拆成一个个独立的源信号，往往会失败。原因在于：PCA 从均方意义上捕捉的是”功率最大”的分量，它只关心二阶统计量（相关性），而”独立”这个要求比”不相关”强得多。于是经常出现这样的情形：PCA 给出的结果恰好是我们想要的，但只要条件稍一变化，它就失效了。
典型的应用场景是盲源分离(Blind Source Separation, BSS)：在一段混合信号里，设法把未知的独立源恢复出来。所谓”盲”，并不是说信号测不到，而是说我们事先不知道源的准确参数、形状或周期——就像上节的 ECG 例子，心跳周期本来就不是先验知识。PCA 对系统条件非常敏感：当两个信号源的功率非常接近、混合矩阵的输入顺序被颠倒，或者数字信号的幅值从别的值变成 1 时，PCA 就可能分不开了。图 8-39(a) 是两个幅值都为 1 的信号经过 PCA 分离的仿真结果，正弦波和数字信号已经无法明确区分。
(a) 采用 PCA 方法
(b) 采用 Herault 和 Jutten 的 ICA 方法
图 8-39 幅值均为 1 的两个信号的 SIMULINK 仿真结果
Herault 和 Jutten 的开创性工作给出了解决这类问题的思路。这一系列方法后来被统称为独立成分分析(Independent Component Analysis, ICA)。他们的关键观察是：引入高阶统计量，许多 PCA 无法完成的 BSS 任务就能迎刃而解。高阶统计量无法用线性运算得到，必须借助非线性函数——ICA 中常用的几种非线性函数如图 8-41 所示。图 8-40 画的是一个双通道系统，其方程为
其中 和 是非线性函数，典型取法是一个为 ，另一个为 。Herault 和 Jutten 还发现，取更简单的 甚至 ，同样能分离信号源。图 8-39(b) 给出幅值均为 1 的数字信号与正弦信号双通道系统的实验结果：开始有一段短暂的学习时间，大约 50 个时钟周期后源信号就分离得比较好了。仿真中前 50 个时钟周期的学习速率为 1/16，随后 50 个时钟周期学习速率逐渐减小到零（冻结权重）。
图 8-40 双通道 Herault 和 Jutten 系统
(b) 低复杂性近似
图 8-41 ICA 中采用的典型非线性方法
尽管 Herault 和 Jutten 给出了一些成功案例，局限性也很快暴露出来：其一，该方法无法推广到多于两通道的系统；其二， 直接和 耦合，增益依赖性太强。随后出现了大量改进算法，成果汇集在教科书、专题报告和工具箱中。有意思的是，互信息最小化、非高斯性最大化和最大似然估计这三大理论殊途同归，得到相同的解决方案。下面简要梳理主流的 ICA 算法。
8.9.1 白噪声化和正交化
许多 ICA 算法要求先对信号做白噪声化预处理：即通过某种变换，使互相关矩阵 除对角元素外全部为零。利用 PCA 方法即可完成白噪声化：
Cxx = cov(x', 1);
[E, D] = eig(covarianceMatrix);
V=D^- .5 * E';
z=V*x;变换后 域的相关矩阵 的非对角元素全为零。不过从上一节可以看到，完整地做 PCA 代价不小。如果只需要白噪声化，有比 PCA 轻得多的做法。一种是与功率法类似的在线学习方式，通过如下两个矩阵方程实现：
另一种颇具吸引力的选择是经典的 Gram-Schmidt 正交化方法，它同样能产生白噪声的互相关矩阵。对三个输入信号 求正交输出 的 MATLAB 代码如下：
% Step K=1
z1 = x1 ./ sqrt(x1'*x1);
z2 = x2 - (z1'*x2)*z1;
v3 = x3 - (z1'*x3)*z1;
% Step K=2
z2 = z2 ./ sqrt(z2'*z2);
z3 = z3 - (z2'*z3)*z2;
% Step K=3
z3 = z3 ./ sqrt(z3'*z3);思路是逐个处理信号：先把第一个信号归一化，再从第二个信号中减去它在第一个信号方向上的投影，如此递推，保证每个新信号都与之前所有信号正交。教科书中通常在最后再做一次标准化，但如果只需要”正交”而不需要”单位长度”，归一化步骤可以省略，硬件上就少开几次平方根。
8.9.2 独立成分分析算法
两类经典 ICA 算法都可以看成 PCA 的非线性扩展，由此得到非线性 PCA 算法，MATLAB 实现为：
yk = W' * z(:,k);
H = f(yk,n1) * (z(:,k)' - f(yk',n1) * W);
W = W + mu * H;非线性函数 用了两次；为方便起见，输入（已白噪声化的）信号按行向量存放。学习速率 通常取得非常小，例如 。
另一类算法同样需要先白噪声化输入信号，称为自然梯度算法(Natural Gradient Algorithm, NGA)，实现为：
yk = W * z(:,k);
H = I - f(yk, n1) * yk';
W = W + mu * H * W;其中 是单位矩阵。与非线性能 PCA 相比，NGA 的复杂度更低——自然梯度相当于在更新式中乘了 ，把搜索方向”扭”到了统计意义上的最速下降方向。
前面这些算法都要先做白噪声化，代价较高。Cardoso 和 Laheld 首次证明：白噪声化可以与 ICA 的高阶统计非线性处理合并在一起完成，这就是基于独立分量分析的等变量自适应分离算法(Equivariant Adaptive Separation via Independence, EASI)：
yk = B * x(:,k);
H = I - yk * yk' + f(yk,n1)*yk' - yk*f(yk,n1);
B = B + mu * H * B;注意：Cardoso 和 Laheld 原文使用 tanh 非线性，而这里直接用 而不是白噪声化后的 ，也就是说 EASI 把”去相关”和”独立性学习”合并进了同一次矩阵更新中，硬件实现自然省了一整级预处理电路。
8.9.3 EASI ICA 算法的实现
在硬件实现上，EASI 是最具吸引力的 ICA 算法。下面设计一个 的 EASI 系统，把例 8.14 中的正弦和数字信号分离开。
例8.15 采用EASI算法的ICA设计
先搭建 SIMULINK 模型生成 HDL 测试数据。仍用位宽为 5 个采样周期长度的二进制数字信号 和周期为 6 的正弦波 ，二者经过混合矩阵 (8-95) 组合成两个系统输入 、。EASI 算法需要更新矩阵 B 和 H 的全部四个系数。对 ， 的各元素为
其中 表示对 施加非线性函数。矩阵 B 的更新方程为
到这里，构建 EASI 系统的全部方程已经齐了，完整的 SIMULINK 系统如图 8-42 所示。图 8-43(a) 的仿真结果中，前 100 个时钟周期学习速率为 1/16，随后 100 个时钟周期逐步减小到零，最后 100 个时钟周期保持为零。与前面讨论的 PCA 或 Herault-Jutten ICA 系统相比，EASI 系统的鲁棒性明显更好：切换输入顺序、调整幅值大小，仍能良好收敛；甚至可以把 tanh 换成更简单的形式，例如 （带饱和的线性函数）或正负号函数（参见图 8-41），在允许多一点残留噪声的前提下依然能分离信号，见图 8-43(b)。
图 8-42 双通道 EASI 系统
(a) 使用 tanh 非线性化
(b) 使用 sign 非线性化
图 8-43 EASI 的 SIMULINK 仿真结果
基于 EASI 算法的 ICA 设计原书以 VHDL 给出，下面将其改写为等价的可综合 Verilog HDL（内部采用 16.16 有符号定点数，即 32 位数据的高 16 位为整数部分、低 16 位为小数部分）：
// 原书 VHDL 改写为 Verilog：双通道 EASI ICA
module ica (
input wire clk, // 系统时钟
input wire reset, // 系统复位
input wire [31:0] s1_in, // 第 1 路信号输入
input wire [31:0] s2_in, // 第 2 路信号输入
input wire [31:0] mu_in, // 学习速率
output wire [31:0] x1_out, // 混合 1. 输出
output wire [31:0] x2_out, // 混合 2. 输出
output wire [31:0] B11_out, // 去混合系数 1,1
output wire [31:0] B12_out, // 去混合系数 1,2
output wire [31:0] B21_out, // 去混合系数 2,1
output wire [31:0] B22_out, // 去混合系数 2,2
output reg [31:0] y1_out, // 第 1. 成分输出
output reg [31:0] y2_out // 第 2. 成分输出
);
// Q16.16 有符号定点类型
typedef signed [31:0] fix_t;
function automatic signed [31:0] mul16_16; // (a*b)>>16，Q16.16 乘法
input signed [31:0] a, b;
begin
mul16_16 = (a * b) >>> 16;
end
endfunction
// 混合矩阵常数
localparam signed [31:0] a11 = 32'sd49152; // 0.75
localparam signed [31:0] a12 = 32'sd98304; // 1.5
localparam signed [31:0] a21 = 32'sd32768; // 0.5
localparam signed [31:0] a22 = 32'sd21845; // 0.333333
localparam signed [31:0] one = 32'sd65536; // 1.0
localparam signed [31:0] negone = -32'sd65536; // -1.0
reg signed [31:0] x1 = 0, x2 = 0;
reg signed [31:0] B11 = 32'sd65536, B12 = 0; // B = I
reg signed [31:0] B21 = 0, B22 = 32'sd65536;
reg signed [31:0] mu_r = 0;
always @(posedge clk or posedge reset) begin
if (reset) begin
x1 <= 0; x2 <= 0;
B11 <= one; B12 <= 0; B21 <= 0; B22 <= one;
mu_r <= 0;
end else begin
mu_r <= $signed(mu_in);
// ---- 混合矩阵 ----
x1 <= mul16_16(a11, $signed(s1_in)) + mul16_16(a12, $signed(s2_in));
x2 <= mul16_16(a21, $signed(s1_in)) - mul16_16(a22, $signed(s2_in));
end
end
// ---- EASI 核：计算 y、H、ΔB 并更新 B ----
reg signed [31:0] y1, y2, f1, f2, H11, H12, H21, H22;
reg signed [31:0] DB11, DB12, DB21, DB22;
always @(posedge clk or posedge reset) begin
if (reset) begin
y1_out <= 0; y2_out <= 0;
end else begin
// 先计算新的 y 分量
y1 = mul16_16(x1, B11) + mul16_16(x2, B12);
y2 = mul16_16(x1, B21) + mul16_16(x2, B22);
// tanh 逼近：f(y)=y, |y|<=1 处饱和
f1 = y1;
if (y1 > one) f1 = one;
else if (y1 < negone) f1 = negone;
f2 = y2;
if (y2 > one) f2 = one;
else if (y2 < negone) f2 = negone;
// 计算 H 矩阵
H11 = one - mul16_16(y1, y1);
H12 = mul16_16(f1, y2) - mul16_16(y1, y2) - mul16_16(y1, f2);
H21 = mul16_16(f2, y1) - mul16_16(y2, y1) - mul16_16(y2, f1);
H22 = one - mul16_16(y2, y2);
// 计算 ΔB
DB11 = mul16_16(B11, H11) + mul16_16(H12, B21);
DB12 = mul16_16(B12, H11) + mul16_16(H12, B22);
DB21 = mul16_16(B11, H21) + mul16_16(H22, B21);
DB22 = mul16_16(B12, H21) + mul16_16(H22, B22);
// 更新矩阵 B（寄存器）
B11 <= B11 + mul16_16(mu_r, DB11);
B12 <= B12 + mul16_16(mu_r, DB12);
B21 <= B21 + mul16_16(mu_r, DB21);
B22 <= B22 + mul16_16(mu_r, DB22);
// 寄存 y 输出
y1_out <= y1;
y2_out <= y2;
end
end
assign x1_out = x1;
assign x2_out = x2;
assign B11_out = B11;
assign B12_out = B12;
assign B21_out = B21;
assign B22_out = B22;
endmodule该设计采用图 8-42 的方案，数据格式为 32 位，内部是 16.16 有符号小数。指定 I/O 之后，先把混合矩阵 定义为常值，再按 Q16.16 格式重新解释输入信号。主体是两个时钟过程：第一级完成混合矩阵运算并寄存 ，第二级执行 EASI 算法——先算 分量（含基于 的 tanh 逼近），接着依次计算 、 和更新后的 。原书中所有定点算术通过 fixed_wrap/fixed_truncate 选项实现饱和与截断，Verilog 版用乘后右移 16 位（>>> 16）完成同样的 Q16.16 运算，代价最小。为了便于测试时序性能， 的结果保存在寄存器中；同时把一些内部数据引到输出端口方便监测。
该设计使用了 2275 个 LE、172 个嵌入式乘法器，用 TimeQuest 缓慢 85C 模型测得时序电路性能 。
现在用 HDL 重新生成图 8-43 所示的 EASI ICA 仿真，结果见图 8-44。前 100 个时钟周期（010μs）学习速率保持不变，随后的 100 个时钟周期（10μs20μs）逐步减小为零，最后 100 个时钟周期（20μs~30μs）保持为零，从那时起信号被分离成两个原始信号。
图 8-44 HDL 中的 EASI 仿真结果：学习(0-10μs)、离散学习速率(10μs20μs)、稳定输出(20μs30μs)
把 EASI 算法与 GHA PCA 设计相比，硬件资源需求和时序性能基本相同，但 EASI 的鲁棒性更好，可以分离比 PCA 更多的信号。
8.9.4 备选BSS算法
Karhunen 等人通过双非线性把更高阶统计量引入 EASI，得到的算法可以用以下 MATLAB 代码描述：
其中图 8-41 所示的非线性形式可以通过参数 n1 和 n2 选择。
EASI 中”把学习和白噪声化结合起来生成健壮算法”这一技巧，同样可以套用到其他 ICA 算法上。例如非线性子空间学习算法按此改进后变为：
该算法同样不需要输入信号先做白噪声化。
可选的 ICA 算法还有很多，再叠加不同的非线性选择，空间更大。例如 FastICA 在微处理器或 MATLAB 上非常流行，因为它收敛快、鲁棒性好；但从硬件实现角度讲，其对称正交化的成本非常高。
还有几种非 ICA 算法被成功用于 BSS 问题，比如 SOBI 和 AMUSE。它们只用二阶统计量、不用非线性，但需要计算带时延的（互）相关。实现简单、仿真效果良好的一个代表是多未知信号提取算法(Algorithm for Multiple Unknown Signals Extraction, AMUSE)，处理过程如下：
% Calculate the covariance matrix with delay
czd = corr(d(z, 1);
czd=(czd+czd')/2; % make it symmetric
% Calculate the eigenvalues and eigenvectors of cov. matrix
[Ed, Dd] = eig(czd);
%apply 2. reconstruction mixing to the y data
yy=real(Dd^- .5 * Ed' * z);其中 z 是预白噪声化的混合信号 x，corr 用于创建 与 的相关矩阵，绝大多数仿真中取 即可稳定工作：
function [C] = corr(x, d)
[r l]=size(z); u=z(:,1:l-d); v=z(:,1+d:l);
for k=1:r
for n=1:r
C(k,n) = xcorr(u(k,:),v(n,:),0,'biased');
end
end总的来说，在 BSS 这类任务中 ICA 提供了更大的选择空间，但也面临不少挑战。选择 ICA 意味着在性能和实现上都要付出较大成本，而且效果依赖于所选非线性与信号峰态的匹配程度。回忆一下，这类算法采用 Hebbian 学习：
符号 取决于信号是次高斯还是超高斯。技术信号通常是次高斯信号，而声音等自然信号通常是超高斯信号（参见表 8-1）。第二个限制是：ICA 无法以确定的阶数给出各成分——输出顺序是随机的；并且与 PCA 不同，ICA 要求源信号与输出信号数量必须相同。
不同算法之间的比较也很有挑战性，涉及学习次数、计算量和残差等多个维度。除了可见的度量之外，很多算法还需要测量混合矩阵与重构矩阵的乘积 。由于输出阶数随机，这个乘积不可能精确等于单位矩阵；当然，乘积结果的每一行、每一列中也只应出现一个占优的非零值，否则说明发生了去相关失败。
8.10 语音和音频信号编码
语音信号的压缩多年来一直是 DSP 领域研究最热的课题之一，涌现出大量复杂而高度自适应的语音和音频压缩算法。这些算法大多需要跨学科知识——不仅来自 DSP，还来自心理声学（这本身就独立成一个研究方向）。道理很直白：对语音编码每提高几个百分点，乘以电话系统中数百万的用户，就是巨大的经济效益。于是各类编码器层出不穷，也定义了多种电信标准，覆盖从高质量音频编码到语音合成的范围。表 8-4 给出了常用方法和标准的总览。
表 8-4 语音和音频压缩算法的比较
| 算法 | PCM | A律 | ADPCM | LPC-10 | MPEG 1 层次1和2 | MPEG 1 层次3 |
|---|---|---|---|---|---|---|
| 年份 | 1948 | 1972 | 1984 | 1984 | 1993 | 1993 |
| 标准 | G.711 | G.722 | FS1015 | MPEG | MPEG | |
| 类型 | 波编码器 | 波编码器 | 波编码器 | 音码器 | 部分波段 | 部分波段+变换 |
| 压缩率 | 1:1 | 1:2 | 1:4 | 1:25 | 1:11 | 1:11 |
| 复杂性 | 低 | 低 | 中 | 高 | 高 | 高 |
| 语音质量 | 高 | 高 | 优 | 优 | 优 | 高 |
| 音频质量 | 高 | 优 | 优 | 低 | 优 | 高 |
最基础的是线性脉冲编码调制(PCM)：以 1216 位精度、8kHz 频率对高品质语音采样，得到 96Kb/s128Kb/s 的数据流。ITU-T 制定的第一个波形压缩方案标准 G.711 采用了对数型幅值压缩：每次采样时，一个 13 或 14 位的语音信号可以压缩到 8 位而品质没有明显变化。与 MP3 等许多其他方案类似，G.711 依赖所谓的心理声学效应——编码方案利用了听觉系统的某些感知特性。下面先来了解 G.711。
8.10.1 A律和 律编码
G.711 利用的心理声学事实是：与较低声级相比，高音量（高分贝）声级的幅值变化更不容易被察觉；换句话说，听觉系统在低分贝时对量化噪声更敏感。G.711 巧妙地反向利用这一点：对低幅值声级用更高的精度（更小的量化噪声）编码，对高幅值声级则允许较大噪声。标准中集成了两种编码方法。 律编码采用的非线性变换为
通常取 。另一种是 A 律编码，在欧洲更流行：
G.711 中取 A=87.56。图 8-45 对两种方法进行了比较：虽然输入数据不同，两条压缩曲线却非常接近。注意 A 律编码覆盖 （即 12 位数据范围），而 律编码覆盖 （约 13 位的正数范围）。
图 8-45 A 律和 律编码的比较
从实现角度看，A 律有一个诱人的优点——可以把它看作一种”短浮点”格式。输出为 8 位，各位布局见表 8-5：第 1 位是符号位，接下来 3 位描述有符号幅值格式下左侧零（即对数”指数”），最后 4 位是幅值的最高有效位（“尾数”）。A 律采用 7 段，4 个幅值位的基本计算方法是
段 2 到 7 各包含 16 个输出值，段 1 包含 32 个值。表 8-5 给出了全部 7 种左零样本情形的完整编码。
表 8-5 采用短浮点格式的 A 律编码(s=符号位, a、b、c、…为各个数据位)
| 段号 | 线性 PCM 输入 | 8 位 A 律编码输出 |
|---|---|---|
| 7 | slabcedfghij | s111abcd |
| 6 | s01abcdedfghi | s110abcd |
| 5 | s00labcedfgh | s101abcd |
| 4 | s000labcdefgh | s100abcd |
| 3 | s00001abcdefg | s011abcd |
| 2 | s000001abcdef | s010abcd |
| 1 | s000000abcdef | s00abcde |
借助直方图，从熵编码的角度看，与原始语音数据相比，A 律或 律编码能更充分地利用可用的信道容量。用简单的 MATLAB 仿真即可验证，如图 8-46 所示：图 8-46(a) 为变化很小或无变化的原始信号，图 8-46(b) 为编码/译码后的 A 律信号——两种信号的波形接近，但直方图看起来完全不同（编码后的分布更加均匀，熵更高）。
图 8-46 原始语音信号和 A 律编码信号的直方图
下面通过一个小的 HDL 设计实现 G.711 标准。
例8.16 A律设计
G.711 编码器和译码器原书以 VHDL 给出，下面改写为等价的可综合 Verilog HDL（采用无流水寄存器的纯组合译码，输入为 13 位符号幅值格式）：
// 原书 VHDL 改写为 Verilog：G.711 A 律编解码器
// A ~= 87.56; |x| <= 4095, 即 12 位幅值加 1 位符号
// mu ~= 255; |x| <= 8160, 即 14 位
module g711alaw #(parameter WIDTH = 13) (
input wire clk, // 系统时钟
input wire reset, // 异步复位
input wire [12:0] x_in, // 系统输入(符号幅值, 非2的补码!)
output reg [7:0] enc, // 编码器输出
output reg [12:0] dec, // 译码器输出
output wire [12:0] err // 结果误差
);
wire s = x_in[WIDTH-1]; // 符号位
wire [12:0] abs_x = {1'b0, x_in[WIDTH-2:0]};
// 迷你浮点格式编码器
always @* begin
case (abs_x)
13'd0 : enc = {s, 7'b0000000}; // +0 或 -0
13'd1 : enc = {s, 2'b00, abs_x[5:1]}; // 段 1
13'd2 : enc = {s, 3'b010, abs_x[5:2]}; // 段 2
13'd3 : enc = {s, 3'b011, abs_x[6:3]}; // 段 3
13'd4 : enc = {s, 3'b100, abs_x[7:4]}; // 段 4
13'd5 : enc = {s, 3'b101, abs_x[8:5]}; // 段 5
13'd6 : enc = {s, 3'b110, abs_x[9:6]}; // 段 6
default : enc = {s, 3'b111, abs_x[10:7]}; // 段 7
endcase
end
// 迷你浮点格式译码器
always @* begin
case (enc[6:4])
3'd0, 3'd1 : dec = {s, 6'b000000, enc[4:0], 1'b1};
3'd2 : dec = {s, 5'b00000, enc[3:0], 2'b10};
3'd3 : dec = {s, 4'b0000, enc[3:0], 3'b100};
3'd4 : dec = {s, 3'b000, enc[3:0], 4'b1000};
3'd5 : dec = {s, 2'b00, enc[3:0], 5'b10000};
3'd6 : dec = {s, 1'b0, enc[3:0], 6'b100000};
default : dec = {s, enc[3:0], 7'b1000000};
endcase
end
// 误差监测
assign err = (dec > x_in) ? (dec - x_in) : (x_in - dec);
endmodule编码器和译码器都按表 8-5 的编码方案实现。由于未使用寄存器，输入、编码和译码数据之间不存在延迟。从图 8-47 的仿真结果可以看到：输入值较小时量化误差小，幅值较大时量化误差也更大——这正是对数量化的本意。由于没有寄存器，时序电路性能也无法测量。
图 8-47 A 律编解码器的 VHDL 仿真结果
本书学习资料提供了语音样本，为便于比较，分别以 SpeechPCM16bit.wav、Speech_A_LAW8bit.wav 和 SpeechPCM8bit.wav 编码。从听感上讲，8 位线性 PCM 的噪声水平远高于 8 位 A 律编码信号，后者几乎保持了原始信号的质量。注意，由于 MATLAB 使用的 WAV 文件格式是 16 位 8kHz（128kbps）的 PCM 数据，资料中所有文件的实际大小是一样的。
8.10.2 线性和自适应 PCM 编码
另一类编码器基于这样一个事实：音频和语音信号通常是缓变的，不像图像那样有尖锐的边沿。因此相邻帧之间相关性很大。如果对”当前帧与上一帧之间的差异”进行编码，需要编码的幅值会变小，编码负担也随之降低。实用上，还需要一个能充分利用幅值新形态的”优秀”量化器。可以设计一种只允许极小变化的系统：对这些微小变化进行累加，这样只需转换或存储增量操作的符号；但如果出现较大的幅值变化，增量方法可能需要若干时钟周期才能跟上。为了防止误差过大，需要在奈奎斯特速率之上施加一个过采样率，比如 4 或 16。这种方法通常被称为增量调制(Delta Modulation, DM)。
另一种方法不用过采样，而是采用一种位数很少、但量化步长能随输入信号自适应调整的量化器，并结合之前若干次采样（即预测方案，如图 8-2 所示）。输入 变化大时，差分信号变大，步长随之变大；信号平稳时，步长自动缩小到比固定小步长 DM 更精细的水平。由于它非常类似 DM、但量化步长具有自适应性，这类系统被称为自适应差分脉冲编码调制(Adaptive Differential Pulse Code Modulation, ADPCM)。完整的编码器和译码器如图 8-48 所示——译码器实际上是编码器的一个内嵌副本（预测器 + 步长自适应部分），因此收发两端不需要传输额外的状态信息。
(a) 编码器
图8-48 ADPCM
ITU 有多种 ADPCM 标准，例如 G.726、G.727 以及使用 ADPCM 的 G.722。与 MP3 或 LPC-10 相比，ADPCM 更容易实现，且比 A 律编码系统获得更好的编码增益。下面以 IMA 倡导的 4 位 ADPCM 作为示例系统，因为与 ITU 标准相比，这款单抽头预测器更容易实现。
例8.17 ADPCM设计
IMA 编码器和译码器原书以 VHDL 给出，下面改写为等价的可综合 Verilog HDL：
// 原书 VHDL 改写为 Verilog：IMA ADPCM 编码器
module adpcm (
input wire clk, // 系统时钟
input wire reset, // 异步复位
input wire signed [15:0] x_in, // 编码器输入
output wire [3:0] y_out, // 4 位 ADPCM 码字
output wire signed [15:0] p_out, // 预测器/译码器输出
output reg p_underflow, // 预测器下溢标志
output reg p_overflow, // 预测器上溢标志
output wire signed [7:0] i_out, // 步长表索引
output reg i_underflow, // 索引下溢标志
output reg i_overflow, // 索引上溢标志
output wire signed [15:0] err, // 系统误差
output wire [14:0] sz_out, // 步长
output wire s_out // 符号位
);
// 步长变化表
function signed [7:0] indexTable;
input [2:0] d;
begin
case (d)
3'd0: indexTable = -1; 3'd1: indexTable = -1;
3'd2: indexTable = -1; 3'd3: indexTable = -1;
3'd4: indexTable = 2; 3'd5: indexTable = 4;
3'd6: indexTable = 6; default: indexTable = 8;
endcase
end
endfunction
// 量化步长查找表(89 项)
function [14:0] stepsizeTable;
input [6:0] i;
begin
case (i)
7'd0 : stepsizeTable = 15'd7;
7'd1 : stepsizeTable = 15'd8;
7'd2 : stepsizeTable = 15'd9;
7'd3 : stepsizeTable = 15'd10;
7'd4 : stepsizeTable = 15'd11;
7'd5 : stepsizeTable = 15'd12;
7'd6 : stepsizeTable = 15'd13;
7'd7 : stepsizeTable = 15'd14;
7'd8 : stepsizeTable = 15'd16;
7'd9 : stepsizeTable = 15'd17;
7'd10: stepsizeTable = 15'd19;
7'd11: stepsizeTable = 15'd21;
7'd12: stepsizeTable = 15'd23;
7'd13: stepsizeTable = 15'd25;
7'd14: stepsizeTable = 15'd28;
7'd15: stepsizeTable = 15'd31;
7'd16: stepsizeTable = 15'd34;
7'd17: stepsizeTable = 15'd37;
7'd18: stepsizeTable = 15'd41;
7'd19: stepsizeTable = 15'd45;
7'd20: stepsizeTable = 15'd50;
7'd21: stepsizeTable = 15'd55;
7'd22: stepsizeTable = 15'd60;
7'd23: stepsizeTable = 15'd66;
7'd24: stepsizeTable = 15'd73;
7'd25: stepsizeTable = 15'd80;
7'd26: stepsizeTable = 15'd88;
7'd27: stepsizeTable = 15'd97;
7'd28: stepsizeTable = 15'd107;
7'd29: stepsizeTable = 15'd118;
7'd30: stepsizeTable = 15'd130;
7'd31: stepsizeTable = 15'd143;
7'd32: stepsizeTable = 15'd157;
7'd33: stepsizeTable = 15'd173;
7'd34: stepsizeTable = 15'd190;
7'd35: stepsizeTable = 15'd209;
7'd36: stepsizeTable = 15'd230;
7'd37: stepsizeTable = 15'd253;
7'd38: stepsizeTable = 15'd279;
7'd39: stepsizeTable = 15'd307;
7'd40: stepsizeTable = 15'd337;
7'd41: stepsizeTable = 15'd371;
7'd42: stepsizeTable = 15'd408;
7'd43: stepsizeTable = 15'd449;
7'd44: stepsizeTable = 15'd494;
7'd45: stepsizeTable = 15'd544;
7'd46: stepsizeTable = 15'd598;
7'd47: stepsizeTable = 15'd658;
7'd48: stepsizeTable = 15'd724;
7'd49: stepsizeTable = 15'd796;
7'd50: stepsizeTable = 15'd876;
7'd51: stepsizeTable = 15'd963;
7'd52: stepsizeTable = 15'd1060;
7'd53: stepsizeTable = 15'd1166;
7'd54: stepsizeTable = 15'd1282;
7'd55: stepsizeTable = 15'd1411;
7'd56: stepsizeTable = 15'd1552;
7'd57: stepsizeTable = 15'd1707;
7'd58: stepsizeTable = 15'd1878;
7'd59: stepsizeTable = 15'd2066;
7'd60: stepsizeTable = 15'd2272;
7'd61: stepsizeTable = 15'd2499;
7'd62: stepsizeTable = 15'd2749;
7'd63: stepsizeTable = 15'd3024;
7'd64: stepsizeTable = 15'd3327;
7'd65: stepsizeTable = 15'd3660;
7'd66: stepsizeTable = 15'd4026;
7'd67: stepsizeTable = 15'd4428;
7'd68: stepsizeTable = 15'd4871;
7'd69: stepsizeTable = 15'd5358;
7'd70: stepsizeTable = 15'd5894;
7'd71: stepsizeTable = 15'd6484;
7'd72: stepsizeTable = 15'd7132;
7'd73: stepsizeTable = 15'd7845;
7'd74: stepsizeTable = 15'd8630;
7'd75: stepsizeTable = 15'd9493;
7'd76: stepsizeTable = 15'd10442;
7'd77: stepsizeTable = 15'd11487;
7'd78: stepsizeTable = 15'd12635;
7'd79: stepsizeTable = 15'd13899;
7'd80: stepsizeTable = 15'd15289;
7'd81: stepsizeTable = 15'd16818;
7'd82: stepsizeTable = 15'd18500;
7'd83: stepsizeTable = 15'd20350;
7'd84: stepsizeTable = 15'd22385;
7'd85: stepsizeTable = 15'd24623;
7'd86: stepsizeTable = 15'd27086;
7'd87: stepsizeTable = 15'd29794;
default: stepsizeTable = 15'd32767;
endcase
end
endfunction
reg signed [15:0] va = 0, va_d = 0; // 当前输入及其延迟
reg sign = 0; // 当前 ADPCM 符号位
reg [3:0] sdelta = 0; // 带符号的 ADPCM 输出
reg [14:0] step = 15'd7; // 步长
reg signed [15:0] valpred = 0; // 预测输出值
reg signed [7:0] index = 0; // 当前步长变化索引
integer k;
reg signed [16:0] diff; // 输入值 - 预测值
reg signed [16:0] p; // 下一个 valpred
reg signed [7:0] i; // 下一个索引
reg [2:0] delta; // 当前 ADPCM 绝对输出
reg [14:0] tStep;
reg signed [15:0] vpdiff; // 对 valpred 的当前修正量
// 输入寄存器
always @(posedge clk or posedge reset) begin
if (reset) begin
va <= 0; va_d <= 0;
end else begin
va <= x_in;
va_d <= va; // 延迟一拍用于误差比较
end
end
// 编码器状态机
always @(posedge clk or posedge reset) begin
if (reset) begin
valpred <= 0; step <= 15'd7; index <= 0;
end else begin
// 状态 1：计算与预测值之差
diff = va - valpred;
sign <= 1'b0;
if (diff < 0) begin
sign <= 1'b1; // 负差值置符号位
diff = -diff; // 量化用绝对值
end
// 状态 2+3：量化(除法)与逆量化
// delta = floor(diff*4/step); vpdiff = floor((delta+0.5)*step/4)
delta = 0;
tStep = step;
vpdiff = tStep[14:3] + tStep[2]; // 近似 tStep/8
if (diff >= tStep) begin
delta = 4; diff = diff - tStep; vpdiff = vpdiff + tStep;
end
tStep = tStep >> 1;
if (diff >= tStep) begin
delta = delta + 2; diff = diff - tStep; vpdiff = vpdiff + tStep;
end
tStep = tStep >> 1;
if (diff >= tStep) begin
delta = delta + 1; diff = diff - tStep; vpdiff = vpdiff + tStep;
end
// 状态 4：按逆量化值更新预测采样
if (sign) p = valpred - vpdiff;
else p = valpred + vpdiff;
// 状态 5：预测值限幅，防止 16 位溢出
p_overflow <= 1'b0;
p_underflow <= 1'b0;
if (p > 17'sd32767) begin
p = 17'sd32767; p_overflow <= 1'b1;
end
if (p < -17'sd32768) begin
p = -17'sd32768; p_underflow <= 1'b1;
end
valpred <= p;
// 状态 6：更新步长与步长表索引
i_overflow <= 1'b0;
i_underflow <= 1'b0;
i = index + indexTable(delta);
if (i < 0) begin // 索引范围 [0...88]
i = 0; i_underflow <= 1'b1;
end
if (i > 88) begin
i = 88; i_overflow <= 1'b1;
end
step <= stepsizeTable(i[6:0]);
index <= i;
// 输出码字 = 符号位 + 幅值
sdelta <= sign ? {1'b1, delta} : {1'b0, delta};
end
end
assign y_out = sdelta; // 监测信号
assign p_out = valpred;
assign i_out = index;
assign sz_out = step;
assign s_out = sign;
assign err = va_d - valpred;
endmodule编码器和译码器都采用图 8-48 的方案；译码器可以看成编码器的后续步骤，故只展示编码器。HDL 先定义 I/O 端口，随后是两张查找表：步长变化表（量化输出越大，步长索引增长越快，从 -1 到 +8）和步长量化表（89 项，从最小步长 7 一直到大步长约 3K）。处理流程分六个状态：第 1 步计算新值与预测值的差；第 2、3 步用除法完成量化并计算逆量化值；第 4 步用逆量化值更新预测采样；第 5 步对预测值限幅以防溢出；第 6 步更新步长，并把索引和步长存入寄存器。编码输出与符号组合成 4 位码字 sdelta。原书 VHDL 中基于整数的查找表在 Verilog 中改用 case 函数实现，综合后即映射为组合逻辑或块 RAM。
图 8-49 的仿真中，输入是一个叠加常值 1000 的斜坡信号。可以看到步长 sz_out 在斜坡段先增大，输入变为恒值后又逐渐减小；斜坡阶段的误差较大，恒值输入阶段误差就非常小了——这正是自适应步长想要的行为。
图 8-49 ADPCM CODEC 的 VHDL 仿真结果
本设计使用了 531 个 LE，未使用嵌入式乘法器，用 TimeQuest 缓慢 85C 模型测得时序电路性能 。
本书学习资料提供了以下语音样本：
- SpeechPCM16bit.wav：16 位编码的数字语音信号。
- Speech_PCM4bit.wav：4 位线性 PCM。
- Speech_APCM4.wav：4 位 ADPCM 编码的语音信号。
- Speech_PCM4Lloyd.wav：4 位 Lloyd 编码的语音信号。
所有文件都保存为 WAV 格式。Lloyd 方法基于所用位数和提供的训练数据计算编码簿，使量化噪声最小（参见图 8-50），在 MATLAB 中通过以下代码计算：
[partition, codebook, distor, rel_distor] = ...
lloyds(training_set, ini_codebook, [], 3);
%% encode + decode
[indx, quant, distor] = quantiz(sig, partition, codebook);
(c) 优化的量化编码分割和编码簿
图 8-50 优化的 Lloyd 量化器
听感上，4 位线性 PCM 的噪声水平比较高；4 位 Lloyd 方法好一点，但仍高于 4 位 ADPCM 信号，后者几乎保持了原始信号的质量。同样注意，由于 MATLAB 的 WAV 文件格式是 16 位 8kHz（128kbps）PCM 数据，资料提供的所有文件实际大小是一样的。
8.10.3 模型化编码：LPC-10e 方法
前面讨论的方法对语音和音频都适用。但如果语音压缩比要超过 4:1，就必须放弃波形编码，改用一种”重建声道、只对声道参数编码”的编码器——编码对象变成音高、共振峰、短时谱以及声道转换函数等参数。语音的生成可以粗略建模为两类：有声语音（有元音，激励为脉冲序列，周期即音高）和无声语音（无元音，激励为随机噪声），声道用滤波器形状建模，如图 8-51 所示。这类编码系统称为线性预测编码器(Linear Predictive Coder)，其中 LPC-10e 是一款非常流行且功能强大的语音编码器。这类系统通常被称作声码器：US DOD 标准 FS1015 在 2400bps 时能获得非常好的效果，即使在 800bps 条件下也成功通过了测试。低于 800bps 时只有语音合成系统才能获得更高压缩比。LPC 声码器的不足是：对音乐等通用音频压缩，效果不尽如人意。
图 8-51 简化的语音生成模型
以 2400bps 传输速率为例，流行的 LPC-10e 编码方案的流程是：以 8kHz 采样输入信号，切成互不重叠的 180 个样本一组，即 44.44 帧/秒。在 2400bps 速率下，每帧只有 54 位可用于编码所有参数。声道滤波器模型采用来自网格 FIR 滤波器的 10 个滤波器系数；音高通过均振幅差函数的自相关估计计算：
如果 是周期为 的函数， 将在 处出现凹点，这个位置就是语音信号的音高周期。LPC-10e 用 7 位编码音高。滤波器系数方面，较低阶的偏相关(Partial Correlation, PARCOR)系数分配更多位，末端的网格滤波器系数分配更少位，总计 位；滤波器增益因子另需 5 位；同步占 1 位。每帧合计 位。LPC-10e 看似复杂，但已有 7.8kbps 和 2.4kbps 的实现案例，例如《ADSP 应用手册：卷 2》。
8.10.4 MPEG 音频编码方法
近年来，MPEG 标准中的音频编码器取得了巨大成功。MPEG1 的第 3 层编码器通常被称为 MP3，是当今网络上使用最频繁的音频格式。由于 MPEG 的主要目标是高品质音频而非单纯压缩语音，其压缩率通常不如 LPC-10e 声码器。复杂性方面的另一个因素是：要保证 CD 音质，音频采样率为 44.1kHz，而语音只需 8kHz。即便如此，11:1 的压缩率下 MP3 依然能生成品质优秀的信号。MPEG 编码器的主要处理单元是一个余弦调制的、32 通道等间隔排列的、长度为 513 的滤波器组，通带滤波器的脉冲响应为
原型滤波器具有 -90dB 的阻带衰减，等间隔滤波器组的通道带宽约为 700Hz。第 3 层编码器进一步改进了频率分解：对每个通道的 12 到 36 次采样使用改进的 DCT，把频率分辨率提高到约 40Hz；此外在最后的压缩环节增加了一个非均匀霍夫曼编码。
MPEG 对所有音频信号都要压缩，但只能依靠基本的心理声学效应，没有声码器系统可用。MPEG 编码器用到的最重要的心理声学效应是频率遮蔽：一个大幅值正弦波会掩盖临近频率的其他信号——听觉系统无法分辨这些被遮蔽的频率分量是否存在。既然听不到，就没有必要把它们编码进去。一般假设遮蔽范围是高频每倍频程 12dB 衰减、低频每倍频程 24dB 衰减。通过长度为 512（MPEG 1+2）或 1024（MPEG 3）的 FFT 确定遮蔽门限，再用它控制信号的量化。图 8-52 给出了 MP3 编码器的结构描述。
图 8-52 又名 MP3 的第 3 个层次的 MPEG1。灰色代表相对于第 1 个层次的改进
MP3 编码器在设计上需要 32 个长度为 513 的滤波器、DCT 和霍夫曼编码器，实现颇具挑战性；对实时嵌入式系统而言过于复杂，因此考虑第 1 层的 MPEG 1 时需要更实际的方案。
8.11 练习
这一节是第 8 章”自适应系统”的配套练习。与其把它当成一份题单，不如把它看成一条由浅入深的学习路径：前面几题练的是自适应滤波器的”数学地基”（功率、自相关、维纳解），中间几题练”数值手感”（用 MATLAB 把统计量算出来），最后几题回到本书的主题——把算法真正放到 FPGA 上跑起来，并评估速度和资源。做题之前，建议先熟悉 Quartus II 的基本流程（可参考 1.4.3 节的案例研究），并注意本书统一使用 Cyclone IV E 系列的 EP4CE115F29C7 器件来评估综合结果；部分练习还会与 Cyclone II 系列的 EP2C35F672C6 器件对比，体会两代器件在嵌入式乘法器和时序性能上的差异。
相关信号的基础题（8.1～8.3）
这三道题练的是自适应系统中最常见的几个统计量：功率（方差）、自相关函数和互相关函数。它们为什么重要？因为后面所有的维纳解、LMS 算法的收敛速度，全部由这些统计量决定。不用会手算这些量，就看不懂为什么自适应滤波器”该收敛到哪”以及”收敛得快不快”。
8.1 纯余弦信号。 给定
(a) 求功率（即方差 ）。对余弦信号逐点求平方再取平均：，其中的交流项在一个周期内平均为零，只剩下直流项，所以
注意这个结果与初相 无关——功率只看振幅。
(b) 求自相关函数 。利用积化和差公式 ，和频项 在 上平均为零，得到
再次注意：自相关函数与初相 无关。这就是说，自相关”看不见”相位，只看见波形的重叠程度。
(c) 由上式直接读出： 是 的余弦，周期为 。信号是周期的，它的自相关也是同周期的——这对后面判断互相关矩阵是否奇异很有用。
8.2 余弦加噪声。 给定
其中 是方差为 的高斯白噪声。
(a) 信号与噪声相互独立，功率直接相加：
(b) 白噪声的自相关是冲激函数 （只在自己身上相关，隔一点就无关），余弦部分沿用 8.1 的结论，于是
这个结果值得多看一眼：在 处自相关有一个”尖峰”（噪声贡献），在 处只有余弦的摆动。白噪声的存在抬高了 ，也就是抬高了自相关矩阵对角线上的元素——这正是白噪声通常让特征值比变差、让 LMS 收敛变慢的原因。
(c) 除 处的冲激外， 的摆动周期仍是 。
8.3 两个不同频率的余弦。 给定 ，。
(a) 互相关函数 。积化和差后得到 。只要 ，这两项都是随 摆动的交流项，在足够长的观测（公共周期）上平均为零。
(b) 因此，只有当 时 才不为零（此时它就是 8.1(b) 的形式）；一般情况下 。结论是：不同频率的确定性信号互不相关。这就是练习 8.6、8.8 中” 某些分量为零”的根源。
维纳解的手算题（8.4、8.5）
这两道题练习二阶系统的”精确解流程”：由统计量构造 和 ，求逆矩阵，得到维纳权重，再算最小误差和特征值比。这条流程在 FPGA 实现中不会真的执行（FPGA 用的是 LMS 这类迭代算法），但它是检验仿真结果是否收敛到正确位置的”标尺”。
8.4 数值题。 已知
(a) 求逆。二阶对称矩阵求逆有固定套路：行列式 ，逆矩阵等于伴随矩阵除以行列式：
(b) 维纳权重 ：
结果非常整齐——题目是特意设计的。
(c) 最优（维纳）权重误差就是最小均方误差：
直觉上：期望信号功率是 20，最优滤波器”解释”掉了 14，剩下的 6 是任何线性滤波器都消不掉的残差。
(d) 特征值。对二阶对称矩阵 形式（本题 ），特征值为 ，即 ，。特征值比
EVR 越接近 1，LMS 收敛越快、越稳；EVR 很大时，收敛速度被最小特征值拖住。这里 EVR = 3 属于”好条件”的情形。
8.5 符号题。 统计量改为
(a) 沿用 (a) 的套路：，
(b) 最优权重误差（最小均方误差）为
也可以写成通用形式 ，代入本题符号即 ；若题目要求”以 为变量的函数形式”，则是 在最优点的取值。
(c) 最优权重本身：
(d) 当 （输入样本互不相关，即”白化”后的理想情形）时，，上式退化成
两个系数完全解耦，各自等于互相关除以输入功率。这解释了一个关键事实：输入越”白”，各系数的收敛越独立、越快；后面用变换（DCT 等）做去相关的练习，动机正在于此。
构造统计量的综合题（8.6～8.8）
这三道题的套路是：先由给定的 和 算出 、、，再走一遍 8.4 的流程。难度在于第一步——要熟练运用三角恒等式和”不同频率信号不相关”的性质。
8.6 展开示例。 ，，即 、。以二阶系统为例演示计算方法：
(a) 输入自相关。 含频率 和 两个成分，两者互不相关，各自贡献自相关。余弦（或正弦）功率为振幅平方的一半：，故 。对 ，正弦成分贡献 ； 成分贡献 ；合计 。于是
互相关 ： 中， 与同频 的相关为零（差一个 90° 相位），与 频率项因频率不同也为零，故 ；而 （不同频项为零）。。
(b)~(d) 之后按 8.4 的流程求 、、、特征值 。这道题的启示是：基准信号里混着的 频率分量只会抬高 的对角线（增大 EVR、拖慢收敛），却对 没有任何贡献——它纯属”干扰”。
(e) 三阶系统把矩阵扩成 （需要 三条对角线），求逆不再有二阶公式，建议交给 MATLAB 的 inv/eig 完成，把主要精力放在构造统计量上。
8.7 含噪期望信号。 （噪声方差为 1），。噪声与 独立，所以 与 8.6 的正弦成分结论一致：，；而 ； 只由正弦决定：，。注意噪声”只进 、不进 和 “——它抬高了误差曲面下限 ，但不改变最优权重的位置。这正是噪声抵消类应用的理论基础：权重照样收敛到正确值，只是残差里多了一份消不掉的噪声功率。
8.8 谐波干扰。 （二倍频），（基频 + 二倍频）。由”不同频率不相关”， 只与 中的二倍频成分相关，且相关值就是 的量级（符号取决于 中该分量的符号）。可以预期最优权重主要用来”重建” ，而一倍频分量在 中几乎不出现。这道题模拟的是谐波对消场景：自适应滤波器学会从混合信号中挑出并抵消目标谐波。
仿真与数值实验题（8.9～8.17）
这批题全部用 C 或 MATLAB 完成，重点是把前几章的理论”跑”出来。
8.9 用 FIR 滤波器制造病态输入。 用表 8-2 的 4 个 FIR 滤波器过滤白噪声（10 000 个样本），再用 xcorr 估计自相关、toeplitz 组装自相关矩阵、eig 算特征值，得到 时的 EVR。一个典型 MATLAB 片段思路：
w = randn(1,10000); % 白噪声
y = filter(b, 1, w); % FIR 滤波（b 为表 8-2 的系数）
r = xcorr(y, L-1, 'biased'); % 自相关估计
R = toeplitz(r(L:end)); % 组装 L 阶 Toeplitz 矩阵
evr = max(eig(R))/min(eig(R)); % 特征值比预期结论：滤波器把白噪声”染色”得越厉害（频谱起伏越大），EVR 越大，而且随 增大有变差的趋势——矩阵越大，越容易”看到”输入频谱中的深谷。
8.10 一阶 IIR（一阶马尔可夫过程）。 用单极点 的 IIR 滤波器重复 8.9，并把仿真 EVR 与理论值
比较。 时该值急剧增大：一阶低通把噪声压到极低频，输入高度相关，自相关矩阵接近奇异。这是理解”DCT/变换域自适应为什么有用”的反面教材——输入越不白，越需要先去相关。
8.11 特征向量与 DCT 基。 对 EVR = 1000 的 FIR 情形，算出 自相关矩阵的特征向量，与 DCT 基向量逐列比较。预期两者高度相似（图形几乎重合）：这说明 DCT 近似实现了 KLT（Karhunen-Loéve 变换）——而 KLT 才是使各系数完全解耦的”最优”变换。这就是第 8 章正文里变换域 LMS 的合法性来源。
8.12 比较各种变换。 对式(8-48)的变换域方法计算功率归一化后变换域自相关矩阵的 EVR，，比较恒等变换、DCT、Hadamard、Haar、KLT。预期排序大致为：恒等（EVR = 1000）最差；Haar 较差（基太简陋）；Hadamard 居中；DCT 接近最优；KLT 理论最优（对角化后 EVR = 1）。(f) 的小问让你构造这些变换的级联，体会”再补一级去相关”的收益递减。
8.13 随极点扫参。 让 在 0.5～0.95 取 10 个值重复 8.12，画出各变换的 EVR 随 的曲线。预期：恒等变换的 EVR 随 接近 1 而爆炸式增长，DCT/KLT 几乎保持平坦——这正是”变换域 LMS 对强相关输入鲁棒”的定量证据。
8.14 非平稳功率估计。 用 C/MATLAB 重构图 8-15 的仿真，比较式(8-41)（普通平均）与式(8-44)（遗忘因子平均）在 和 下的表现。 越大，历史权重越大：对缓慢变化的信号平滑好，对突变的信号跟踪慢。 跟踪快但抖动大， 平滑但滞后——这就是非平稳环境下”平滑 vs 跟踪”的经典折中。
8.15～8.17 仿真复现。 8.15、8.16 分别复现例 8.1、8.3 的仿真（采样率取 次/秒，即对 和 每周期采 个点），。8.15 的 (d)(e) 还要求对 手算精确维纳解，用来对照 LMS 的收敛终点。8.17 复现例 8.6 的 DLMS 仿真，流水线级数 ，观察”延迟误差”随 增大对收敛的影响： 越大，有效步长必须越小、稳态误差越大，但硬件时钟越快——这组数据是后面 FPGA 设计的”预期答案”。
FPGA 设计与资源评估题（8.18～8.22）
最后四题回到硬件：改滤波器长度或流水线级数，完成”功能编译 → 功能仿真 → 对照 MATLAB 结果 → 时序/资源评估”的完整闭环。虽然原书给出的是 VHDL 代码，建议改用等价的 Verilog HDL 描述实现，两者综合后资源占用基本一致，Verilog 对初学者更常见。评估时注意三点：
- ：用 TimeQuest 静态时序分析，选缓慢 85C 模型（最坏工艺/电压/温度角），报告设计能跑的最高时钟。流水线级数增加（8.20～8.22 中 ）会显著提高 ，因为组合逻辑被寄存器切短了。
- LE（逻辑单元）：L 从 2 加到 4（8.18、8.19），乘加结构数量线性增加，LE 数随之上升。
- 嵌入式乘法器和 M9K：例 8.5、8.6 的设计把 、 等关键乘法映射到嵌入式乘法器（18×18 硬核），把延迟线等存储放进 M9K 块 RAM。报告值应显示 LE 占用下降、硬乘法器占用等于并行乘法个数——这正是”用专用资源换通用逻辑”的验证。
8.20～8.22 分别对应”只流水 （）”、“流水 和 （）”、“再在 后加一级且 不流水（）“三种修改。做完后与 8.17 的 MATLAB 仿真逐项对照：功能仿真中误差信号 的收敛轨迹应与对应 值的 DLMS 仿真一致（收敛步长需按 折减）；若不一致，多半是流水线寄存器位置与延迟补偿逻辑不匹配。最后把 、LE、乘法器、M9K 在 EP4CE115F29C7（Cyclone IV E）与 EP2C35F672C6（Cyclone II）上各报一遍，体会较新一代器件在同类设计上的性能与资源优势。
小结
这组练习覆盖了自适应系统设计的四个层次：(1) 统计量手算——弄清 、、维纳解、EVR 之间的因果链；(2) 构造统计量——把具体信号变成矩阵，理解”频率不同即不相关""噪声只抬高 “这类结构性结论；(3) 数值实验——用仿真验证染色输入、变换去相关、非平稳跟踪、DLMS 延迟等理论预言；(4) 硬件落地——在 FPGA 上验证流水线换速度、专用资源换面积的工程规律。建议按 8.1 → 8.4 → 8.6 → 8.9 → 8.12 → 8.18 的顺序做，前一层结果是后一层的对照标准，做完后再回头看第 8 章正文，很多公式会豁然开朗。