LPC共振峰提取：从全极点模型到杜宾递推与牛顿求根

简介：面向语音信号处理学习者与嵌入式C开发者，这份工程围绕线性预测编码（LPC）实现语音共振峰提取，覆盖语音分帧、汉明窗加权、端点检测、Durbin递推与牛顿迭代求解预测系数、由预测多项式求根定位共振峰等关键环节，可用于语音识别、语音合成与通信系统的教学实验和原型验证。压缩包共72个文件、约274KB，包含多组C源码、工程配置与编译映射及日志文件、数据表文件、wav语音样本和说明文档，便于从预处理到特征提取完整梳理工程结构。目前已有385人学习，适合作为课程设计、算法复现与二次开发的参考，也可帮助读者理解LPC参数估计、加窗降边界效应和共振峰判定的实现细节。

1. 峰值不等于共振峰：为什么 LPC 要绕道去解多项式

做语音的人几乎都干过同一件事：把一帧语音的频谱画出来，在包络上找那几个鼓起来的包，标成 F1、F2、F3。这个做法在教科书示意图上永远正确，在真实数据上却经常出问题——两个靠得很近的共振峰会融合成一个宽包，带宽很窄的极点又会被激励谱的倾斜推离真实位置，基频高的说话人谐波间距大，能把一整个共振峰藏进去。

所以成熟的 LPC 共振峰提取实现从不做峰值搜索。用 C语言 落地时它走的是反方向：先用杜宾算法把一帧语音拟合成 p 阶全极点滤波器，再把分母多项式 A(z) 的根逐个解出来，只有落在特定环形区域里的复根，才有资格被换算成频率和带宽。多绕这两步，换来的是可解释、可复现、帧间连续的输出——而不是一条随基频抖动的曲线。

2. 全极点模型：LPC 拟合的到底是声道的哪一段

2.1 一根管子为什么能用全极点函数描述

声门产生准周期脉冲串，气流经过咽腔、口腔，最后由唇端辐射出去。这套管道可以近似成一串截面积突变的声学管，每一处突变产生一次反射，多个反射叠加形成谐振。集总参数近似下，每一对谐振对应传递函数的一个共轭极点对，而整条管道就是这些极点对的级联——注意是级联，不是并联，所以分母是各因式的乘积，写成 z 域就是：

H(z) = G / (1 + a1*z^-1 + a2*z^-2 + ... + ap*z^-p) = G / A(z)

这里没有任何零点。物理上唯一的零点来自唇端辐射，那是一个近似一阶微分 (1 - μz⁻¹) 的高通效应，而它会在预处理阶段被预加重提前抵消掉——这也是预加重存在的第一理由，后面还会再说。

提示：全极点假设成立的前提是"没有额外零点"。鼻音因为鼻腔分支引入了真正的零点，摩擦音的湍流源不在声门，所以 LPC 在这两类音上的拟合误差天生就大。工程上一般只对浊音段做共振峰提取，清音段直接跳过。

2.2 预测误差最小化怎么变成 Yule-Walker 方程

LPC 的思路是反过来问：如果声道真是个全极点系统，那么当前样本 s(n) 应该能由前面 p 个样本线性预测出来。设预测值为

ŝ(n) = -a1*s(n-1) - a2*s(n-2) - ... - ap*s(n-p)

误差 e(n) = s(n) - ŝ(n) = a0s(n) + a1s(n-1) + ... + ap*s(n-p)，其中 a0 = 1。把这串误差的平方和对所有样本求和，对每个 aₖ 求偏导并令其为零，就得到一组线性方程：

Σ(k=1..p) aₖ * R(|i-k|) = -R(i), i = 1..p R(k) = Σ(n=k..N-1) s(n) * s(n-k)

这叫 Yule-Walker 方程，也叫正规方程。它的系数矩阵是一个托普利茨矩阵——每条对角线上的元素都相同，因为 R 只跟下标差有关。这个结构非常关键：一般矩阵求逆要 O(p³)，托普利茨矩阵可以降到 O(p²)，而且能保证解出来的滤波器是最小相位的。杜宾算法吃的就是这个结构红利。

2.3 阶数：p ≈ fs/1000 + 2 这条经验线的由来

p 阶多项式有 p 个极点，实数系数的多项式根成共轭对出现，所以理论上有 p/2 个谐振。但语音频带内真正需要描述的只有 F1 到 F5 这四五个，剩下的阶数是留给频谱倾斜、声门脉冲形状这些"非谐振"成分的。

采样率越高，可描述的频率范围越宽，需要的阶数也越多。常见的经验是 p ≈ fs/1000 + 2，8kHz 采样取 12 阶，10kHz 取 14 阶，16kHz 取 18 到 20 阶。

阶数不足的后果非常直观：F1 和 F2 靠得近时，两个极点会被强行拟合进一个宽极点里，结果只报出一个峰，或者两个频率都被拉到中间。阶数过高的后果更隐蔽——多余的自由度会去拟合基频的谐波结构，冒出一堆带宽极窄（低于 20Hz）的伪极点，看起来像共振峰，实际上每次基频变化它就换个位置。这就是为什么后面筛选规则里必须有一条带宽下限。

3. 前端处理链：预加重、汉明窗与自相关

3.1 预加重：先把 -6dB/oct 的谱倾斜扳平

语音频谱有一个天然的下倾趋势，低频能量比高频高出一大截，大致是每倍频程掉 6dB。这个倾斜来自声门脉冲的谱形和唇端辐射的微分效应。如果直接把它喂给 LPC，最小均方误差准则会把绝大部分建模能力花在低频段，因为那里的误差绝对值最大，高频共振峰基本被放弃。

预加重就是一级一阶高通：y(n) = x(n) - μ·x(n-1)。μ 通常取 0.9 到 0.97，8kHz 采样常用 0.95，采样率越高越接近 1。它恰好抵消了唇端辐射带来的那个零点，同时把高频抬起来，让各频段的误差量级接近。

注意：预加重是双刃剑。它在抬高高频语音的同时也抬高了高频噪声，信噪比差的录音做完预加重可能更糟。实际调试时如果发现 F3、F4 抖得厉害，把 μ 从 0.95 降到 0.9 往往是第一个该试的动作。

3.2 汉明窗：旁瓣 -43dB 换来的到底是什么

分帧之后必须加窗，否则帧边界处的突变会在频谱上产生严重的泄漏。矩形窗的主瓣最窄（2π/N），频率分辨率最好，但第一旁瓣只有 -13dB，意味着一个强低频成分会把自己的"影子"铺到很远的高频去。语音的动态范围轻松超过 30dB，F1 完全可能把 F4 淹没在泄漏里。

汉明窗把第一旁瓣压到 -43dB，代价是主瓣展宽到 4π/N，频率分辨率减半。这个代价在 LPC 里是可以接受的，因为最终频率精度不是由窗的主瓣宽度决定的，而是由极点位置决定的——窗只影响自相关序列，不影响求根的精度。

帧长的选择取决于基音周期。窗内至少要覆盖两到三个基音周期，LPC 才能把激励和声道分开，否则拟合出来的模型会跟着单个脉冲跑。成人基频 80 到 300Hz，对应周期 3 到 12ms，所以帧长取 20 到 30ms 比较稳妥：8kHz 采样下就是 160 到 240 点，常用 200 点配 80 点帧移。

还有一个容易被忽略的约束：帧长 N 必须远大于阶数 p。一般要求 N ≥ 4p，8kHz 下 p=12、N=200，余量充足。如果为了提时间分辨率把帧长压到 80 点还硬用 12 阶，自相关矩阵接近奇异，杜宾递推中途就会崩。

3.3 自相关：为什么只算到 p 就够了

Yule-Walker 方程里只用到 R(0) 到 R(p)，共 p+1 个值，多算一个都是浪费。直接按定义算是 O(N·p)，8kHz 下 N=200、p=12，一次 2400 次乘加，嵌入式设备上也毫无压力；如果帧长很长（比如 1024 点做离线分析），走 FFT 转成功率谱再做逆变换会更快，复杂度 O(N log N)，但短帧反而不划算。

/* 一帧语音的前端处理：预加重 -> 加窗 -> 自相关 * x : 输入样本（已按帧移切好） * n : 帧长，8kHz 下取 200（25ms） * p : LPC 阶数 * mu : 预加重系数，8kHz 常用 0.95 * r : 输出自相关，长度 p+1 * 返回 r[0]，即帧能量，可用于清浊/静音判决 */ double frame_autocorr(const double *x, int n, int p, double mu, double *r) { double s[256]; /* 生产代码应由调用方传入缓冲区，避免逐帧 malloc */ int i, k; /* 1) 预加重：y(n) = x(n) - mu*x(n-1)，抵消唇端辐射的一阶零点 */ s[0] = x[0]; for (i = 1; i < n; i++) s[i] = x[i] - mu * x[i - 1]; /* 2) 汉明窗：w(n) = 0.54 - 0.46*cos(2*pi*n/(N-1))，旁瓣压到 -43dB */ for (i = 0; i < n; i++) s[i] *= 0.54 - 0.46 * cos(2.0 * M_PI * i / (n - 1)); /* 3) 自相关：只需 0..p 这 p+1 个值 */ for (k = 0; k <= p; k++) { double acc = 0.0; for (i = k; i < n; i++) acc += s[i] * s[i - k]; r[k] = acc; } return r[0]; }

这段代码里三个参数各司其职：mu控制谱倾斜补偿的强度，n决定了频率分辨率和基音周期覆盖度，p决定模型能容纳多少个谐振。返回值 r[0] 是加窗后的帧能量，后面做清浊判决时用得上——r[0] 低于某个阈值就直接判定为静音，根本不用进杜宾。

4. 杜宾递推与牛顿求根：从自相关到共振峰频率的两段接力

4.1 杜宾递推：O(p²) 和 |kᵢ|<1 的稳定性保证

正规方程直接解要高斯消元，O(p³)。p=12 时是 1728 次运算，其实也不多，但问题不在速度而在数值条件：托普利茨矩阵的条件数可能很差，消元过程中的舍入误差会让解出来的滤波器不稳定。

Levinson-Durbin 递推换了个思路，从一阶开始逐级往上加：已经算好了 i-1 阶的最优预测器，现在加第 i 阶，新的系数只需要在旧系数基础上做一次"旋转"，不需要从头解方程。每一级的中间量 kᵢ 有个明确的物理身份——第 i 级格型滤波器的反射系数。

这个身份带来两个好处。第一，复杂度降到 O(p²)，p=12 时只有 78 次乘加。第二，|kᵢ| < 1 是格型结构稳定的充要条件，而自相关矩阵的半正定性保证了递推出来的 kᵢ 必然满足这个条件。所以自相关法 + 杜宾递推，天然产出最小相位滤波器，不需要额外检查。

另一个可用的中间量是预测误差能量 Eᵢ = (1 - kᵢ²)·Eᵢ₋₁。它单调递减，如果中途算出 Eᵢ ≤ 0 或者 |kᵢ| ≥ 1，说明输入数据出了问题（帧太短、定点溢出、自相关不是真正半正定），这时候与其硬算下去，不如直接丢弃这一帧。

/* Levinson-Durbin 递推：自相关 -> LPC 系数 * r : 自相关，r[0..p] * p : 阶数 * a : 输出，a[0]=1，a[1..p] 为预测系数，对应 A(z)=1+sum(a[k]*z^-k) * 返回最终预测误差 E_p；中途异常返回 -1 */ double durbin(const double *r, int p, double *a) { double tmp[32]; double E = r[0]; int i, j; if (E <= 0.0) return -1.0; /* 静音帧，无模型可拟合 */ for (i = 0; i <= p; i++) a[i] = 0.0; a[0] = 1.0; for (i = 1; i <= p; i++) { double acc = r[i]; for (j = 1; j < i; j++) acc += a[j] * r[i - j]; double k = -acc / E; /* 第 i 级反射系数 */ if (fabs(k) >= 1.0) return -1.0; /* 非最小相位，数值异常 */ /* a_j^(i) = a_j^(i-1) + k * a_{i-j}^(i-1)，先写临时数组避免覆盖 */ for (j = 1; j < i; j++) tmp[j] = a[j] + k * a[i - j]; for (j = 1; j < i; j++) a[j] = tmp[j]; a[i] = k; E *= (1.0 - k * k); /* 误差能量单调递减 */ } return E; }

注意更新系数时用了临时数组。这一步不能原地做，因为 a[j] 的新值依赖 a[i-j] 的旧值，当 j 和 i-j 交叉时会互相污染。这是手写杜宾最容易出错的地方，符号约定（k 前面那个负号）也是另一个——如果 A(z) 定义成 1 - Σaₖz⁻ᵏ，符号要整体反过来，混用会导致求出来的根全部跑到单位圆外。

4.2 牛顿迭代：为什么不能在幂基上直接开根

拿到 a[1..p] 之后，要解的是 P(z) = z^p + a₁z^(p-1) + ... + aₚ = 0。这是把 A(z) 乘以 z^p 得到的，去掉了 z=0 的平凡根，剩下的根和 A(z) 完全相同。

p=12 时这是个 12 次多项式，没有解析解。可选路线有几条：一是 QR 算法求全部特征值，鲁棒但要写一大坨 Householder 变换和隐式位移；二是把多项式因式分解成二次因式（Bairstow 法），实现简单但收敛慢；三是复数牛顿法，只需要一个初值加 Horner 求值，二十行代码。

牛顿法的迭代式是 z ← z - P(z)/P'(z)。复平面上它的收敛是二次的，只要初值落在根的吸引域内，五六次迭代就能到 1e-9 精度。麻烦在于"只要"这两个字——牛顿法只保证局部收敛，初值选偏了会跑到无穷远或者卡在重根附近。

工程上的解法是在半径 0.95 的圆上均匀分布 p 个初值，每个根至少会被其中一个吸引到。这个策略没有严格的理论保证，但实践中对 LPC 多项式足够可靠，因为 LPC 多项式的根本来就集中在单位圆内侧，分布相对均匀。更稳的做法是先用上一帧的极点做初值——语音的共振峰在 10ms 内变化很小，这个"热启动"能把迭代次数压到两三次，顺带保证帧间频率连续。

typedef struct { double r, i; } cplx; /* 不用 complex.h，方便移植到嵌入式 */ static cplx c_add(cplx a, cplx b) { return (cplx){a.r + b.r, a.i + b.i}; } static cplx c_sub(cplx a, cplx b) { return (cplx){a.r - b.r, a.i - b.i}; } static cplx c_mul(cplx a, cplx b) { return (cplx){a.r*b.r - a.i*b.i, a.r*b.i + a.i*b.r}; } static cplx c_div(cplx a, cplx b) { double d = b.r*b.r + b.i*b.i; return (cplx){(a.r*b.r + a.i*b.i)/d, (a.i*b.r - a.r*b.i)/d}; } static double c_abs(cplx a) { return sqrt(a.r*a.r + a.i*a.i); } /* Horner 同时求 P(z) 与 P'(z)，一次遍历 O(p) * P(z) = z^p + a[1]*z^(p-1) + ... + a[p] */ static cplx poly_eval(const double *a, int p, cplx z, cplx *dp) { cplx P = {1.0, 0.0}, D = {0.0, 0.0}; for (int k = 1; k <= p; k++) { D = c_add(c_mul(D, z), P); /* D = D*z + P */ P = c_add(c_mul(P, z), (cplx){a[k], 0.0}); } *dp = D; return P; } /* 复牛顿迭代求单根，成功返回迭代次数，失败返回负值 */ static int newton_root(const double *a, int p, cplx z0, cplx *out) { cplx z = z0, dp, f; for (int it = 0; it < 30; it++) { f = poly_eval(a, p, z, &dp); if (c_abs(dp) < 1e-12) return -1; /* 导数为零，退化成重根 */ cplx dz = c_div(f, dp); z = c_sub(z, dz); if (c_abs(dz) < 1e-9) { *out = z; return it; } } *out = z; return -2; /* 30 次未收敛，丢弃 */ }

Horner 那五行是整个求根的核心。它从最高次项开始，每轮做一次"乘 z 加系数"，顺便把导数也递推出来——D 的新值等于 D 的旧值乘 z 再加 P 的旧值，这正好是乘积求导法则的递推形式。一轮循环同时得到 P 和 P'，比分别算两遍省一半乘法。

4.3 极点转频率/带宽：F = θ·fs/2π，B = -ln(r)·fs/π

解出来的根是复数，要变成有物理意义的频率和带宽，靠的是拉普拉斯域和 z 域的映射关系。一个谐振在 s 平面上的位置是 s = -πB + j2πF，其中 F 是中心频率、B 是带宽。映射到 z 平面 z = e^(sT)，T = 1/fs，得到

z = exp(-πB/fs) * exp(j*2πF/fs) = r * exp(jθ)

所以模长 r = exp(-πB/fs)，辐角 θ = 2πF/fs，反解就是：

F = θ * fs / (2π) B = -ln(r) * fs / π

这两个公式说明了一件重要的事：极点的模长决定带宽，辐角决定频率，两者完全解耦。极点越靠近单位圆，r 越接近 1，带宽越窄，谐振越尖锐；极点跑到单位圆上，带宽为零，系统临界稳定；跑到圆外，系统发散。这就是为什么筛选时必须检查 r < 1。

筛选规则通常包括这几条：只取上半平面（虚部为正），因为实系数多项式的根成共轭对，下半平面是重复信息；模长必须严格小于 1；频率落在 90Hz 到 fs/2-50Hz 之间；带宽落在 20Hz 到 700Hz 之间。带宽下限用来干掉谐波伪极点，上限用来干掉那些根本没形成谐振的宽极点。

typedef struct { double freq; double bw; } Formant; /* 从 LPC 系数提取共振峰，返回个数（按频率升序写入 out） */ int lpc_to_formants(const double *a, int p, double fs, Formant *out, int max_out) { int n = 0; for (int k = 0; k < p; k++) { /* 初值撒在半径 0.95 的上半圆，避开实轴附近的根 */ double th = M_PI * (k + 0.5) / p; cplx z0 = { 0.95 * cos(th), 0.95 * sin(th) }; cplx z; if (newton_root(a, p, z0, &z) != 0) continue; /* 未收敛则丢弃 */ if (z.i <= 1e-9) continue; /* 共轭对称，只取上半平面 */ double r = c_abs(z); if (r <= 0.0 || r >= 1.0) continue; /* 必须在单位圆内 */ Formant f; f.freq = atan2(z.i, z.r) * fs / (2.0 * M_PI); f.bw = -log(r) * fs / M_PI; if (f.freq < 90.0 || f.freq > fs * 0.5 - 50.0) continue; if (f.bw < 20.0 || f.bw > 700.0) continue; /* 去重：不同初值可能收敛到同一个根 */ int dup = 0; for (int i = 0; i < n; i++) if (fabs(out[i].freq - f.freq) < 20.0) { dup = 1; break; } if (dup || n >= max_out) continue; out[n++] = f; } /* 插入排序，保证 F1 < F2 < ...，帧间才好比较 */ for (int i = 1; i < n; i++) { Formant t = out[i]; int j = i - 1; while (j >= 0 && out[j].freq > t.freq) { out[j + 1] = out[j]; j--; } out[j + 1] = t; } return n; }

去重那一段不能省。牛顿法从不同初值出发，完全可能收敛到同一个根，尤其是带宽窄、吸引域大的极点。20Hz 的阈值是个经验值，比正常人耳能分辨的频率差小得多，不会误伤真正相邻的两个共振峰——F1 和 F2 在极端情况下能接近到 100Hz 以内，仍然远大于这个阈值。

排序也不是可有可无的装饰。不排序的话，同一段元音相邻两帧输出的 out[0] 可能是 F1，也可能是 F3，下游做平滑或者做轨迹跟踪时会看到剧烈的跳变，误以为算法不稳定。按频率排序之后，out[0] 永远是当前帧最低的那个共振峰，帧间才有

本文还有配套的精品资源，点击获取