简介：这套 LPC 法语音共振峰提取 C 语言工程面向语音信号处理初学者、嵌入式开发者与高校实验教学场景，围绕线性预测编码、杜宾递推、牛顿迭代、语音分帧与汉明加窗等环节，演示从端点检测到共振峰参数估计的完整实现路径。压缩包共 72 个文件，约 274KB，以 4 个 C 源文件及配套 pjt、lkf、map、out、obj 等工程与编译文件为主体，另有 txt 说明、wav 语音样本和数据库索引类文件，可用于在 CCS3.3 环境中复现、调试与对照分析。已有 385 人学习，适合需要理解 LPC 求解流程、比较杜宾算法与牛顿算法差异、观察窗函数和工程配置细节的读者。工程保留了分模块的源码与构建记录，便于按语音预处理、系数求解、根求解和特征标记逐段研读，也为课程设计、算法验证和二次开发提供可运行的参考底稿。

1. 峰值最高的那几个点，从来不是共振峰

很多人第一次做 LPC 共振峰提取，都以为谱包络上最高的几个峰就是 F1、F2。真正运行起来才发现：杜宾算法递推出来的系数很稳定，包络画出来也很漂亮；可一旦用牛顿算法求根，频点就开始乱跳，静音帧里还能凭空冒出五六个"共振峰"。问题不在算法本身，而在这条链路上每一环都在做取舍：汉明窗要压掉多少旁瓣、阶数该留几个自由度、求根时先给谁当初值。用 C 语言把整条流程写一遍，这些取舍就躲不掉了，也没有现成库替你处理。理解清楚这些取舍，也就理解了"人耳如何从一段波形里认出元音"的可计算版本。

2. 预加重、分帧与汉明窗：为什么分析前要先"毁掉"一点信号

共振峰提取的输入本来很干净：一段 16 位、16 kHz 的单声道波形。但这条流水线的前三步——预加重、分帧、加窗——每一步都在主动丢弃信息。丢得有道理，因为不丢的话，LPC 的最小二乘会把误差预算花在错误的地方，你想要的 F3、F4 就会被低频淹掉。这一章把这三次"破坏"的动机和代价拆开讲。

2.1 预加重：把 -6 dB/oct 的谱斜坡拉平

声门脉冲串的频谱大致按 -12 dB/oct 衰减，唇端辐射近似一个微分器，贡献 +6 dB/oct，两者抵消后，录音里还剩大约 -6 dB/oct 的整体倾斜。这个倾斜对 LPC 是灾难性的：低频能量比 3 kHz 以上高出二三十 dB，最小二乘求解时误差几乎全部集中在低频段，高频那几对极点就退化成噪声拟合出来的伪结构。

预加重用一个一阶高通把斜坡抬平，y[n] = x[n] - a*x[n-1]，零点在 z = a。系数怎么选有个粗糙的经验：8 kHz 采样常用 0.95，16 kHz 常用 0.97。系数越大，高频抬得越狠；抬过头会把量化噪声和房间噪声一起放大成假峰，我见过 0.99 的参数把底噪抬成一条"第六共振峰"。常见的实现思路是先看一眼语谱图，如果 4 kHz 以上明显发灰，说明抬过头了。

#include <math.h> #define FS 16000 /* 8 kHz 语料要改这里 */ #define FRAME_LEN (FS * 25 / 1000) /* 25 ms -> 400 点 */ #define FRAME_SHIFT (FS * 10 / 1000) /* 10 ms -> 160 点 */ #define PRE_EMPH 0.97f /* 预加重：逆序遍历，避免 x[i-1] 被提前覆盖 */ static void pre_emphasize(float *x, int n, float a) { for (int i = n - 1; i > 0; --i) x[i] -= a * x[i - 1]; /* x[0] 保留原值；若想严格等价于 x[-1] = 0，可令 x[0] *= (1.0f - a) */ } /* 汉明窗表：分母用 N-1 而不是 N，保证首尾对称 */ static void make_hamming(float *w, int n) { for (int i = 0; i < n; ++i) w[i] = 0.54f - 0.46f * cosf(2.0f * (float)M_PI * i / (n - 1)); } /* 取第 k 帧并加窗，结果写入 out[0..FRAME_LEN-1] */ static void get_frame(const float *x, int k, const float *w, float *out) { int start = k * FRAME_SHIFT; for (int i = 0; i < FRAME_LEN; ++i) out[i] = x[start + i] * w[i]; }

pre_emphasize逆序循环这个细节值得单独说：正序写会把x[i-1]改成预加重之后的值，等于把一阶差分做成了递归滤波，频响完全变样。这类错误不会报错，只会让共振峰整体偏高。make_hamming里分母取n-1而不是n，是为了让窗首尾严格对称；用n会在自相关序列上引入一个微小的线性相位偏斜，低阶时看不出来，阶数上到 16 就能看到 F1 漂十几赫兹。get_frame只做乘法，不加任何判断，因为越界检查应该放在调用层做一次，而不是每帧每个样本判一次。

2.2 分帧与帧移：25 ms 背后是时间-频率的硬折中

声道形状的变化速率有限，10~30 ms 之内可以近似当成时不变系统，这就是短时分析的物理前提。窗长定下来，频率分辨率也就定了：汉明窗的主瓣宽度约4π/N，换算成赫兹是2*FS/N。16 kHz、400 点，主瓣宽度大约 160 Hz。

这个数字决定了你能分到什么程度。成年男性 /ɑ/ 的 F1 约 730 Hz、F2 约 1090 Hz，间距 360 Hz，两个主瓣勉强分得开；但 /u/ 的 F1 只有 300 Hz 左右、F2 约 870 Hz，或者某些发音人 F1/F2 挤在 200 Hz 以内，25 ms 的窗就开始"抹平"，两个峰并成一个。反过来把窗拉到 50 ms，频率分辨率好看了，但一帧里跨过了塞音的爆破，得到的共振峰是把两个音素平均出来的产物，轨迹上全是台阶一样的假跳变。

我一般会这么定：8 kHz 语料用 30 ms（240 点），16 kHz 用 25 ms（400 点），帧移统一取窗长的 40%，也就是 10 ms。帧移小于窗长意味着帧间有重叠，这对后面的轨迹跟踪是好事——相邻帧的极点位置变化足够小，才谈得上连续性约束。

2.3 汉明窗的 0.54/0.46：旁瓣抑制与谱平滑的双重身份

矩形窗的第一旁瓣只有 -13 dB，衰减还慢，一个 300 Hz 的强 F1 能把自己的旁瓣铺到 2 kHz 以上，把本来就弱的 F3 盖住。汉明窗第一旁瓣 -43 dB，代价是主瓣变宽一倍，且远端滚降只有 -6 dB/oct（汉宁窗远端是 -18 dB/oct，但第一旁瓣只到 -31 dB）。共振峰提取更在意近端旁瓣，所以选汉明是常见做法。

更关键的一层作用常被忽略：加窗之后，自相关序列实际上是"信号自相关"与"窗自相关"的卷积，等价于在频域对功率谱做了一次平滑。而谱分裂（一个宽共振峰被拆成两个尖峰）本质上就是谱包络上出现了不该有的高频起伏——平滑正好压住它。换句话说，汉明窗在 LPC 里不只是防泄漏，还承担了正则化的角色。这也是为什么"去掉窗试试"这种实验往往能得到更尖锐、但完全不可用的极点：不是精度变高了，是没人管的谱分裂被放出来了。

3. 自相关与杜宾递推：LPC 究竟在"预测"什么

到这里信号已经被切成一段段 400 点的短时序列。接下来的问题是：给定这 400 个样本，怎么用十几个系数描述产生它们的声道？杜宾算法就是那个把"最小化预测误差"这个优化问题变成一次 O(p²) 递推的东西。它之所以重要，不只是快，而是顺带把稳定性判据一起算出来了。

3.1 全极点模型：声道为什么能写成 1/A(z)

源-滤波器模型把一个音素拆成三部分：声门激励、声道共振、唇端辐射。声道可以看作一根截面渐变的管子，用 p 段无损管级联近似时，从声门到唇端的传递函数天然是有理分式，且分子是常数——也就是全极点。零点的物理来源主要是鼻腔耦合和唇端的辐射细节，能量占比小，多数情况下被吸收进极点里也不影响前三个共振峰的位置。

形式上，H(z) = G / A(z)，其中A(z) = 1 - Σ a_k z^-k（k 从 1 到 p），G 是增益。LPC 做的事就是找一组a_k，让x[n]用前 p 个样本线性预测的均方误差最小。注意A(z)的根就是H(z)的极点，而且是 A(z) 在单位圆内的根——这个约定很重要，后面求根时所有的符号都和它绑定。

3.2 自相关法：为什么共振峰提取几乎不用协方差法

求解a_k有两条经典路线。协方差法不加窗，直接用帧内数据构造非 Toeplitz 的协方差矩阵，精度更高，但解出来的极点可能跑到单位圆外，对应的带宽是负数。共振峰提取最怕这个：一个不稳定的极点会让整个 LPC 谱的形状失真。

自相关法先用加窗后的序列算出自相关r[k]，矩阵变成 Toeplitz，Levinson-Durbin 递推保证解是最小相位的，极点一定在圆内。代价是加窗带来的谱平滑，但对共振峰来说这个代价是划算的——我们要的是包络，不是精确的谱拟合。

#define MAX_ORDER 20 /* 自相关 r[0..p]。注意是有偏估计：求和上限固定为 N-1， 等价于给 r[k] 乘了一个三角形窗，反而让极点估计更稳 */ static void autocorr(const float *x, int n, int p, double *r) { for (int k = 0; k <= p; ++k) { double s = 0.0; for (int i = k; i < n; ++i) s += (double)x[i] * x[i - k]; r[k] = s; /* 未归一化，r[0] 即帧能量 */ } }

内层用double累加不是洁癖：400 个float相乘再相加，误差累积在 1e-5 量级，杜宾递推的分母E会一路除下去，到第 12 阶时反射系数能偏出 1e-3，正好卡在稳定性判据的边界上。代价是慢一点，一帧 400×13 次乘加，现代 CPU 上不到 10 微秒，完全值得。

3.3 Levinson-Durbin：三个不变量和一个副产品

递推的核心是把 p 阶问题拆成 p 个一阶问题，每一步只求一个反射系数k_i。

/* 输入 r[0..p]，输出 a[1..p]（A(z)=1-Σa_k z^-k）和最终预测误差 *err 返回 0 成功，负数表示该帧不可用 */ static int durbin(const double *r, int p, double *a, double *err) { double prev[MAX_ORDER + 1]; double E = r[0]; if (E <= 1e-9) return -1; /* 静音帧，直接拒绝 */ for (int i = 0; i <= p; ++i) a[i] = 0.0; for (int i = 1; i <= p; ++i) { double acc = r[i]; for (int j = 1; j < i; ++j) acc -= a[j] * r[i - j]; double k = acc / E; /* 第 i 个反射系数 */ if (k >= 1.0 || k <= -1.0) return -2; /* 数值失稳，极点出圈 */ for (int j = 0; j <= i; ++j) prev[j] = a[j]; a[i] = k; for (int j = 1; j < i; ++j) a[j] = prev[j] - k * prev[i - j]; /* 对称更新 */ E *= (1.0 - k * k); /* 误差能量单调递减 */ } *err = E; return 0; }

三个不变量可以当自检用：E必须单调递减且恒正；每个|k_i| < 1；更新式里prev[i-j]的下标必须关于 i 对称。任何一条不成立，说明输入有问题而不是算法有问题——最常见的是r[0]太小（静音）或者采样数据里混进了直流偏置（r[k]几乎不衰减）。

副产品比主产品还好用：err / r[0]就是归一化预测误差。浊音帧这个值通常在 0.01~0.1，清音和静音帧接近 1。不必另做一套 VAD，直接用这个比值门限就能把不该出共振峰的帧筛掉，这是整条链路上成本最低的一道保护措施。

3.4 阶数怎么定：每对极点覆盖约 1 kHz

经验规则是p ≈ FS / 1000 + (2~4)。8 kHz 取 8~10，10 kHz 取 10~12，16 kHz 取 14~16。背后的道理是：一个二阶谐振器大致能描述 1 kHz 宽度内的谱包络起伏，而人耳能分辨的共振峰一般到 5 kHz 为止；多出来的 2~4 阶留给谱倾斜、声源细节和非共振的极点，它们最后会被带宽过滤掉。

阶数给高了，谱分裂就会露出来：一个本来 200 Hz 宽的 F1 被两对极点共同拟合，拆成两个相距 150 Hz 的尖峰，求根之后你会看到 F1 和"F1.5"。阶数给低了，F4/F5 直接合并，而且 F1/F2 会系统性地偏高十几赫兹——因为模型被迫用更少的自由度去拟合同样的包络，峰的位置被推着走。这个偏差是系统性的，不是噪声，所以调阶数的时候不能只看单帧，要统计一批稳定元音帧的均值。

4. 牛顿法求复根：从 LPC 系数到 F1~F5 的最后一步

拿到A(z)的系数只是有了一个隐式的谱包络，真正要输出的是"这一帧有五个共振峰，频率分别是多少、带宽分别是多少"。这一步在数学上是求一个 p 次多项式的复根，在工程上是跟数值误差要求较高。

4.1 为什么不能直接在 LPC 谱上找峰

最省事的做法是算S(ω) = 1/|A(e^{jω})|²的离散谱，然后找局部极大值。这条路在 F1/F2 相距较远时能用，一遇到靠近的就崩：两个极点在频率轴上互相"推"，谱峰的频率不等于极点角度，峰的数量也不等于极点对数。更麻烦的是极点之间的肩部会被误判成峰，你没法用峰值高度区分它和真共振峰。

求根的价值在于它把每个极点单独拿出来，顺便给了带宽。带宽是判断真伪的关键指标：真实共振峰的 -3 dB 带宽通常在 50~400 Hz，而谱分裂产生的伪峰、或者拟合噪声产生的极点，带宽往往超过 500 Hz，甚至极点半径接近 1 导致带宽只有几赫兹——两种情况都不是共振峰。

4.2 复 Newton 迭代：用正幂形式，别碰 z^-1

A(z) = 1 - Σ a_k z^-k在 z = 0 处有多重极点，直接对它求导做 Newton，数值上极不稳定。改用Ã(z) = z^p · A(z) = z^p - a_1 z^(p-1) - ... - a_p，根完全不变，首项系数归一为 1，Horner 法可以一次遍历同时算出函数值和导数值。

typedef struct { double re, im; } cplx; /* Horner：Ã(z) = 1*z^p + c[1]*z^(p-1) + ... + c[p]，其中 c[k] = -a[k] */ static cplx poly_eval(const double *c, int p, cplx z) { cplx v = { 1.0, 0.0 }; for (int i = 1; i <= p; ++i) { double re = v.re * z.re - v.im * z.im + c[i]; double im = v.re * z.im + v.im * z.re; v.re = re; v.im = im; } return v; } /* 同步算导数：d = d*z + v，v = v*z + c[i] */ static cplx poly_eval_deriv(const double *c, int p, cplx z) { cplx v = { 1.0, 0.0 }, d = { 0.0, 0.0 }; for (int i = 1; i <= p; ++i) { double dre = d.re * z.re - d.im * z.im + v.re; double dim = d.re * z.im + d.im * z.re + v.im; double vre = v.re * z.re - v.im * z.im + c[i]; double vim = v.re * z.im + v.im * z.re; d.re = dre; d.im = dim; v.re = vre; v.im = vim; } return d; } /* 复 Newton：z <- z - Ã(z)/Ã'(z)，收敛返回 1 */ static int newton_root(const double *c, int p, cplx *z) { for (int it = 0; it < 50; ++it) { cplx v = poly_eval(c, p, *z); cplx d = poly_eval_deriv(c, p, *z); double den = d.re * d.re + d.im * d.im; if (den < 1e-18) return 0; /* 导数为零，驻点 */ double dr = -(v.re * d.re + v.im * d.im) / den; double di = -(v.im * d.re - v.re * d.im) / den; z->re += dr; z->im += di; if (dr * dr + di * di < 1e-14) return 1; /* |Δz| < 1e-7 */ } return 0; }

迭代上限 50 次看着宽松，实际绝大多数根五六次就收敛了。真正决定成败的是初值：Ã(z)有 p 个根，Newton 只会收敛到离初值最近的那个。我一般会在单位圆内、半径 0.92 的位置均匀撒 128~256 个角度，算|Ã|找局部极小点当初值。半径取 0.92 而不是 1.0，是因为极点都在圆内，初值落在圆内更不容易被相邻的根抢走；角度密集一些是为了保证两个靠得很近的极点各有各的初值。

4.3 降阶与共轭配对：找到一根就"除掉"它

不复用初值、不对多项式降阶，你会反复收敛到同一个根，然后误以为这一帧只有两个共振峰。标准做法是每找到一根z0，用二次因子z² - 2Re(z0)z + |z0|²对Ã做综合除法，把阶数减 2，再重新粗搜初值。这样每个根只会被找到一次。

/* 用二次因子 z^2 + D*z + C 去除 Ã，商写回 c[]，阶数减 2 */ static void deflate_quad(double *c, int *p, double re, double im) { int n = *p; double D = -2.0 * re; double C = re * re + im * im; double b[MAX_ORDER + 1]; b[0] = c[0]; /* 恒为 1 */ b[1] = c[1] - D * b[0]; for (int k = 2; k <= n - 2; ++k) b[k] = c[k] - D * b[k - 1] - C * b[k - 2]; /* 余项应为 0，调试时可算出来当收敛判据： rem1 = c[n-1] - D*b[n-2] - C*b[n-3]; rem0 = c[n] - C*b[n-2]; */ for (int k = 0; k <= n - 2; ++k) c[k] = b[k]; *p = n - 2; }

实根要单独处理：Ã的实根对应一个不构成谐振的极点（带宽极大），它对共振峰没贡献，但会占用一次迭代。用线性因子z - r降一阶即可，别硬套二次因子。另外要注意，降阶用的是求出来的复根，不是共轭对里随便挑一个——D和C都是实数，所以取re和|im|的结果一样，但如果你在别处缓存了根，记得保证实部来自同一个根。

4.4 从根到频率和带宽：两个换算公式

设极点为z0 = r·e^(jθ)，它对应连续域里的s = -πB + j2πF，经过z = e^(sT)映射后：

• 频率F = θ · FS / (2π)，θ 取atan2(im, re)，范围在 (0, π) 内对应 0 到奈奎斯特频率；
• 带宽B = -ln(r) · FS / π，半径越接近 1，带宽越窄。

下半平面和上半平面的极点共轭成对，物理上是同一个共振峰，只记一次，取im > 0的那一半。典型的成年人元音（男性）大致落在下表范围，超出很多的频点基本可以怀疑是伪峰：

typedef struct { double f, bw; } Formant; /* 输入 c[]：Ã 的系数，c[0]=1，最高阶 *p；输出按频率升序 */ int lpc_formants(double *c, int p, double fs, Formant *out, int max_out) { int n = p, cnt = 0; while (n >= 2) { double theta = coarse_peak_angle(c, n); /* 圆内粗搜，见正文 */ cplx z = { 0.92 * cos(theta), 0.92 * sin(theta) }; if (!newton_root(c, n, &z)) break; /* 不收敛，整帧作废 */ if (fabs(z.im) < 1e-6) { /* 实根，非谐振 */ deflate_linear(c, &n, z.re); continue; } if (z.im > 0.0) { /* 共轭对只记一次 */ double r = hypot(z.re, z.im); if (r < 1.0 && cnt < max_out) { out[cnt].f = atan2(z.im, z.re) * fs / (2.0 * M_PI); out[cnt].bw = -log(r) * fs / M_PI; ++cnt; } } deflate_quad(c, &n, z.re, z.im); } qsort_by_freq(out, cnt); return cnt; }

注意c[]传递进来的是正幂形式的系数，需要先从杜宾的输出转换：c[0] = 1.0; c[i] = -a[i];。这一步的符号搞反，求出来的根会是A(z)的零点的倒数，落在圆外，表现为"所有带宽都是负数"。排序之后还要做一次过滤：带宽大于 500 Hz 的丢掉，相邻两个频率差小于主瓣宽度（2·FS/N）的合并成一个，剩下的才是可以往上送轨迹的 F1~F5。

5. 避坑与排查：五个真实翻过车的点

算法本身不复杂，复杂的是数值边界。下面这几条都是调试日志里反复出现的模式，按"现象 → 原因 → 解决"列出来，方便对照。

5.1 静音帧里冒出五六个共振峰

现象：一段前导静音，每帧都能输出 5 个频点，而且帧与帧之间完全没有相关性。原因：r[0]只有 1e-6 量级，杜宾递推里E一路除下去，反射系数变成随机的，极点散落在圆内各处。求根照常能解出根，只是这些根毫无物理意义。解决：两道门。第一道在自相关之后，r[0] < 阈值直接跳过；第二道在杜宾之后，用err / r[0] > 0.6判为非浊音帧。阈值不要写死绝对值，用整段语音的r[0]中位数乘一个比例，否则换一个录音设备就要重调。

5.2 相邻帧的 F1 和 F2 突然交换

**

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