纯C实现LPC语音共振峰提取：杜宾递推与牛顿求根

简介：面向语音信号处理学习者与嵌入式开发者的 LPC 语音共振峰提取 C 语言工程，围绕线性预测编码、杜宾递推与牛顿迭代求解预测系数，结合分帧加汉明窗、端点检测等环节，演示从语音样本到共振峰特征的完整处理链路。压缩包共 72 个文件，约 274KB，包含 C 源码、工程与编译配置文件、日志及映射输出、语音样本和说明文本，便于在 CCS 类集成环境中打开、编译并对照调试。已有 385 人学习，适合课程设计、语音识别入门实验及算法验证场景。读者可据此理解 LPC 模型建立、预测系数估计、窗函数处理与共振峰提取之间的衔接，并借助现成工程结构复现实验、修改参数观察结果，为后续语音特征分析或嵌入式语音应用开发提供可运行的参考基础。

1. 从一段 16 位 PCM 说起：LPC 法语音共振峰提取到底能给你什么

手头是一批 16 位单声道录音，要求每 10 毫秒给出 F1、F2、F3 三个数值，用来画元音轨迹。环境是一台只能跑纯 C 的机器，装不了带一堆依赖的语音工具箱，于是这套 LPC 法语音共振峰提取的 C 语言实现就被逼出来了。整条链路不复杂：读 wav、预加重、分帧、加汉明窗、算自相关、用杜宾算法递推 LPC 系数、再用牛顿算法在复平面求根，最后把根换算成频率和带宽。全篇不到一千行，只依赖 libm，编译出来一个可执行文件丢到哪都能跑。

它适合两类人：一类是被 Python 语音库的安装环境折磨过、想拿一份纯 C 实现做嵌入或移植的开发者；另一类是要弄明白"LPC 系数凭什么能变出共振峰"、想自己改阶数和窗函数做实验的人。如果你只想调包出结果，这份代码反而显得啰嗦；但你要真想在 C 里把这条链路走通，它的每个函数都值得抠一遍。

2. 工程骨架与信号流：从 wav 读入到共振峰输出的五级流水线

整套代码按数据流切成五段，每段一个独立函数，互不依赖全局状态。这种切法不是出于架构洁癖，而是因为调试时你一定会想单独抓某一级的中间结果——把第 3 级的 LPC 系数 dump 出来和成熟工具对一遍，比对着最终频率瞎猜快得多。

2.1 目录结构与各文件职责

工程目录保持扁平，六个源文件加一个 Makefile，刻意不引入任何第三方库：

lpc_formant/ ├── Makefile ├── src/ │ ├── wav_io.c # wav 读写，只处理 44 字节标准头的 16bit PCM │ ├── preproc.c # 去直流、预加重、分帧、加窗 │ ├── lpc.c # 自相关 + 杜宾递推，输出 a[1..p] 和预测误差 │ ├── root.c # 复系数多项式求值、牛顿迭代求根 │ ├── formant.c # 根 -> (频率, 带宽)，筛选与排序 │ └── main.c # 串流程、打印 CSV └── data/ └── vowel_a.wav # 一段自录或自合成的测试元音

wav_io.c只认最常见的 16 位 PCM，不去实现 RIFF 的完整解析。遇到fmt块长度不是 16、或者带LIST块的文件，直接用fseek跳到data块再读，这两行处理能免掉后面一大半的"读出来全是噪声"的排查。

preproc.c里预加重系数固定 0.97，这是 8 kHz 时代留下的经验值，16 kHz 用着也没问题；分帧长度和帧移做成运行时参数从命令行传，方便你后面扫参。

lpc.c是整条链路的数学核心，杜宾递推就三十行，但符号约定一旦搞错，后面求根求出来的频率会整体偏移。root.c和formant.c拆开，是因为求根是整个工程里最容易被数值问题咬的地方，单独一个文件方便加日志。

2.2 预加重、分帧与汉明窗：三个不能随手改的参数

语音信号经过声门后高频衰减大约 6 dB/倍频程，预加重就是把这个斜率抬平，让 LPC 不至于把有限阶数全浪费在拟合低频包络上。代码只有一行，但省掉它，高阶极点会明显往低频挤。

/* preproc.c: 去直流 + 预加重，原地处理 */ void pre_emphasis(double *x, int n, double alpha) { /* 先去直流：一阶高通，截止极低，主要清掉录音设备的直流偏置 */ double mean = 0.0; for (int i = 0; i < n; ++i) mean += x[i]; mean /= n; for (int i = 0; i < n; ++i) x[i] -= mean; /* 预加重必须倒序，否则 x[i-1] 已经被覆盖 */ for (int i = n - 1; i > 0; --i) x[i] = x[i] - alpha * x[i - 1]; x[0] = x[0] * (1.0 - alpha); /* 首点近似处理，避免能量突变 */ }

alpha是预加重系数，常用 0.9～0.97，取值越大高频补偿越强，超过 0.98 会把量化噪声一起抬起来。倒序循环这一步是血泪经验：正序写会让每个采样点都用上已经被改过的前一点，等价于把滤波器变成了递归结构，高频会被过度放大，听起来像加了齿音。

分帧紧接着做。16 kHz 采样下，帧长取 25 ms 即 400 点，帧移取 10 ms 即 160 点。帧长太短，自相关估计的方差大，共振峰会跳；太长，跨越音素边界，两个元音的极点糊在一起。

/* preproc.c: 生成一帧并加汉明窗 */ int make_frame(const double *x, long total, int start, int flen, const double *win, double *out) { if (start + flen > total) return -1; /* 尾部不足一帧直接丢弃 */ for (int i = 0; i < flen; ++i) out[i] = x[start + i] * win[i]; /* 加窗在此一步完成 */ return 0; } /* 汉明窗系数表，启动时算一次，后面每帧复用 */ void init_hamming(double *win, int flen) { for (int i = 0; i < flen; ++i) win[i] = 0.54 - 0.46 * cos(2.0 * M_PI * i / (flen - 1)); }

汉明窗的作用是把帧两端的突变压下去，否则自相关序列会带上矩形窗的 sinc 旁瓣，表现为共振峰附近出现假的极点。公式里分母是flen - 1而不是flen，两者差别很小，但和 MATLAB 的hamming对齐时就看这个细节。

提示：窗函数表一定要在循环外算好。每帧重算一次cos，在 400 点帧长下要多花几十万次三角函数调用，批量处理几千个文件时差距肉眼可见。

参数上我一般这样起步：

阶数这条最容易被忽略。经验公式是p ≈ fs / 1000 + 2，16 kHz 下取 18 左右，但求根运算量随阶数平方增长，实际工程里 14～16 已经够把前三个共振峰分出来，再高只是多出一堆需要筛掉的伪极点。

2.3 自相关与杜宾递推：LPC 系数是怎么被"推"出来的

加完窗的帧送进lpc.c。第一步算自相关r[k] = Σ x[n]·x[n-k]，第二步用杜宾算法解 Yule-Walker 方程。为什么不用高斯消元直接解？因为 Yule-Walker 的系数矩阵是托普利茨矩阵，杜宾递推把 O(p³) 压到 O(p²)，而且递推过程中每一级的反射系数k顺带就是稳定性的判据——|k| < 1时保证所有极点都在单位圆内，这比解完方程再去验根省事太多。

/* lpc.c: 自相关，顺手做数值保护 */ void autocorr(const double *x, int n, int order, double *r) { for (int k = 0; k <= order; ++k) { double s = 0.0; for (int i = k; i < n; ++i) s += x[i] * x[i - k]; r[k] = s; } /* 白噪声校正：给 r[0] 加一个极小的相对量，缓解病态 */ r[0] = r[0] * (1.0 + 1e-9) + 1e-12; }

自相关算完，r[0]就是这一帧的能量。静音帧的r[0]接近 0，后面杜宾递推第一步就要除以它，直接算出 inf。所以工程上一定在进递推之前先判一次能量门限，低于门限的帧标记为"无共振峰"跳过去，比让它在浮点里发散再事后补救干净得多。

3. 核心代码逐段拆解：杜宾算法与 LPC 系数的落地实现

杜宾递推的教科书版本只有五行公式，但落在 C 里要处理三件事：符号约定、中间数组的复用、以及递推过程中的数值崩溃。三件事里任何一件出错，最终表现都是"频率算出来了，但和预期差几百赫兹"，这种错误最难查。

3.1 符号约定先定死，再谈实现

LPC 有两套完全相反的符号约定。一套写成x[n] = Σ a_j·x[n-j] + e[n]，另一套写成x[n] + Σ a_j·x[n-j] = e[n]。前者算出来的a是预测系数，后者算出来的是预测误差滤波器的系数，两者每一项都差一个负号。求根的时候，多项式是A(z) = 1 + a_1 z^-1 + … + a_p z^-p，符号一错，根的相位就整体镜像，频率照样能算出来，但 F1 和 F2 会对调。

这套代码统一采用第二种：a[0]恒等于 1，a[1..p]是误差滤波器系数，多项式按A(z)直接构造。所有推导和 MATLAB 的levinson输出对齐，方便你把中间结果导出来逐项比对。

3.2 杜宾递推的完整 C 实现

/* lpc.c: Durbin-Levinson 递推 * 入参 r[0..p] 自相关，p 阶数 * 出参 a[0..p]，a[0] = 1.0；*err 为最终预测误差 * 返回 0 成功，-1 数值异常（能量塌陷或反射系数越界） */ int durbin(const double *r, int p, double *a, double *err) { double tmp[MAX_ORDER + 1]; double e = r[0]; if (e <= 1e-12) return -1; /* 静音帧，拒绝递推 */ for (int i = 0; i <= p; ++i) a[i] = 0.0; a[0] = 1.0; for (int i = 1; i <= p; ++i) { /* 1) 求第 i 级反射系数 k */ double k = -r[i]; for (int j = 1; j < i; ++j) k -= a[j] * r[i - j]; k /= e; /* 2) 反射系数越界，说明自相关矩阵不正定，提前收手 */ if (k <= -1.0 || k >= 1.0) return -1; /* 3) 更新系数，注意必须先写 tmp 再整体拷回 */ tmp[i] = k; for (int j = 1; j < i; ++j) tmp[j] = a[j] + k * a[i - j]; for (int j = 1; j <= i; ++j) a[j] = tmp[j]; /* 4) 更新预测误差 */ e *= (1.0 - k * k); if (e <= 0.0) return -1; } *err = e; return 0; }

这段逻辑分四步走。第一步算反射系数，注意k的初值是-r[i]，负号来自前面定死的符号约定，改成+r[i]会得到一套"看起来也能跑"的系数，但求根结果的实部虚部关系是反的。第二步的区间检查不是可选项：当帧内含强烈周期性（比如录音里有明显嗡声）时，自相关序列可能出现非正定，k冲出(-1, 1)，接下来的递推会在几步内把系数放大到 1e10 量级。

第三步的tmp数组是关键。系数更新要用到上一轮的a[i-j]，如果直接原地改a[j]，当j走到i-j那一项时读到的已经是新值，结果全错。用tmp暂存再整体拷回是标准做法，代价是一次p长度的内存复制。

第四步的e是预测误差，它应该单调递减且始终为正。一旦e变负，说明前面的数值已经崩了，这时候继续算只会得到一堆无意义的极点，不如直接返回 -1 让调用方跳过这一帧。

参数上，p的取值直接决定能分辨出多少个共振峰。理论上p个极点最多对应p/2个共振峰，要稳定拿到 F1～F3 至少需要 8 阶，实际留余量取 12～16。MAX_ORDER定义成 32，够覆盖到 16 kHz 采样下的高阶实验，再往上栈上的tmp数组就要改堆分配了。

3.3 系数稳定性检查与常见翻车点

杜宾递推自带稳定性保证：只要每级|k| < 1，得到的A(z)所有零点都在单位圆内。但这是数学性质，不是浮点性质。实际跑起来有两种情况会让这条保证失效。

第一种是长帧加高阶数。帧长 400 点、阶数 16 时，自相关矩阵的条件数可以到 1e6 以上，k在最后几级容易出现0.9999这种临界值，递推出来的系数在双精度下勉强能看，换成float就直接发散。

第二种是录音里有低频轰鸣。预加重只压了高频衰减，没管低频，50 Hz 附近的强分量会让r[1]、r[2]异常大，第一级反射系数就顶到 0.99 以上，后面所有极点被这一个强分量带偏。

处理方法我一般两手都上：一是在自相关之后加一个滞后窗（lag window），对r[k]乘一个高斯衰减exp(-0.5·(2π·k·bw)²)，bw取 60～80 Hz，把高阶自相关的权重压下去；二是在杜宾返回 -1 时不要把整帧丢掉，而是降阶重试，p从 16 降到 12 再试一次，通常就过了。全丢掉会让共振峰轨迹出现大片空洞，后续插值反而引入假峰。

注意：降阶重试会改变极点数量，输出前一定要记录本帧实际用的阶数。我曾经在批量跑的时候忘了这一步，后面统计"每帧三个共振峰"时发现有 8% 的帧只有两个，追了半天才发现是降阶路径没打标记。

4. 从 LPC 系数到共振峰：多项式求根与牛顿迭代的工程写法

拿到a[1..p]只是半程。LPC 谱的包络峰对应A(z)在单位圆内的零点，而共振峰就是这些零点的频率和带宽。求根这一步没有闭式解，只能数值迭代，牛顿法是最直接的选择——收敛快，但初值敏感，需要工程上的补丁。

4.1 构造多项式与复数运算底座

A(z) = 1 + a_1 z^-1 + … + a_p z^-p两边乘z^p，得到P(z) = z^p + a_1 z^{p-1} + … + a_p，这是一个首一的p次多项式，系数直接就是a的倒序排列。求它的p个复根即可。

C 里没有内置复数（C99 的complex.h在部分编译器上支持不全），所以先搭一层极简复数结构：

/* root.c: 最小复数工具，只实现求根需要的四个运算 */ typedef struct { double re, im; } cplx; static cplx c_add(cplx a, cplx b) { cplx c = {a.re + b.re, a.im + b.im}; return c; } static cplx c_sub(cplx a, cplx b) { cplx c = {a.re - b.re, a.im - b.im}; return c; } static cplx c_mul(cplx a, cplx b) { cplx c = {a.re * b.re - a.im * b.im, a.re * b.im + a.im * b.re}; return c; } static cplx c_div(cplx a, cplx b) { double d = b.re * b.re + b.im * b.im; if (d < 1e-30) { cplx bad = {0.0, 0.0}; return bad; } /* 分母塌陷保护 */ cplx c = {(a.re * b.re + a.im * b.im) / d, (a.im * b.re - a.re * b.im) / d}; return c; } static double c_abs(cplx a) { return sqrt(a.re * a.re + a.im * a.im); }

c_div里那个1e-30的判据不是形式主义。牛顿迭代中导数值接近 0 时，除法会把根弹到 1e8 以外，下一轮求值直接溢出成nan，而nan参与比较永远返回假，整个筛选逻辑会静默地把所有根都当成有效根输出。

多项式求值用霍纳法：

/* c[0..p] 为降幂系数，c[0] 必须为 1.0 */ static cplx polyval(const double *c, int p, cplx z) { cplx v = {c[0], 0.0}; for (int i = 1; i <= p; ++i) { cplx ci = {c[i], 0.0}; v = c_add(c_mul(v, z), ci); } return v; }

霍纳法把p次求值压到p次乘加，比逐项算幂快一个量级，而且避免了z^p在|z| ≈ 1时的精度损失。

4.2 牛顿迭代求复根：初值、收敛判据与降阶

/* root.c: 牛顿法求一个复根 * c[0..p] 降幂系数，z0 初值，max_iter 上限 * 成功写 *out 返回 0，失败返回 -1 */ int newton_root(const double *c, int p, cplx z0, cplx *out, int max_iter) { double dc[MAX_ORDER]; for (int i = 0; i < p; ++i) dc[i] = c[i] * (double)(p - i); /* 逐项求导，次数降一 */ cplx z = z0; for (int it = 0; it < max_iter; ++it) { cplx f = polyval(c, p, z); cplx df = polyval(dc, p - 1, z); if (c_abs(df) < 1e-12) return -1; /* 导数过小，换初值重来 */ cplx dz = c_div(f, df); z = c_sub(z, dz); if (c_abs(dz) < 1e-9) { *out = z; return 0; } /* 步长足够小 */ if (!isfinite(z.re) || !isfinite(z.im)) return -1; } return -1; }

三个参数决定成败。z0初值：我在单位圆上均匀取p个点，半径取 0.9，即z0 = 0.9·exp(j·2πk/p)，k = 0..p-1。半径小于 1 是因为共振峰极点必然在单位圆内，从外面往里收敛比从圆心往外跑稳定。max_iter取 50，正常情况 5～8 轮就收敛，取 50 是给病态帧留余量。

收敛判据用步长而不是函数值，因为P(z)在根附近的量级可能是 1e-6，用固定阈值判函数值会误判。步长小于 1e-9 时，对应的频率误差已经在 0.001 Hz 量级，远小于物理意义。

真正的坑在"两个初值收敛到同一个根"。复数多项式在单位圆内有p个根，但牛顿法的收敛域不保证互斥，p个初值里可能有三四个收敛到同一个根，剩下的根就漏了。工程上的补法是降阶（deflation）：每求出一个根z_r，就用综合除法把多项式除以(z - z_r)，次数减一，再对降阶后的多项式求下一个根。

/* root.c: 综合除法，c 为降幂系数，zr 为已知根 * 结果写回 c，长度减一，c[0] 保持 1.0 */ void deflate(double *c, int p, cplx zr) { double nr[MAX_ORDER]; nr[0] = 1.0; for (int i = 1; i <= p - 1; ++i) nr[i] = c[i] + zr.re * nr[i - 1]; /* 只保留实部路径，见下文说明 */ for (int i = 0; i <= p - 1; ++i) c[i] = nr[i]; }

严格来说综合除法要处理复数系数，但 LPC 多项式是实系数多项式，复根必然成对共轭出现，所以每求出一个复根，同时把它的共轭也去掉，降两阶、c保持实系数。这一步省掉了复数系数递推，代码短一半，代价是必须先验证求出的根确实成对——只要输入是实系数、durbin的数值没崩，这个前提基本成立。

4.3 由根到频率带宽：筛选与排序

根求出来之后，换算成物理量。单位圆内的复根z = r·exp(jθ)，θ 对应角频率，-ln(r)对应带宽：

/* formant.c: 极点 -> (频率, 带宽)，并做有效性筛选 */ int root_to_formant(const cplx *roots, int p, double fs, Formant *out, int max_out) { int n = 0; for (int i = 0; i < p; ++i) { double mag = c_abs(roots[i]); /* 1) 单位圆外或贴圆的根一律丢弃 */ if (mag >= 1.0 || mag < 1e-6) continue; /* 2) 角频率 -> Hz，注意 atan2 返回 [-pi, pi] */ double f = atan2(roots[i].im, roots[i].re) * fs / (2.0 * M_PI); /* 3) 带宽：极点离单位圆越远，带宽越宽 */ double bw = -log(mag) * fs / M_PI; /* 4) 物理范围过滤 */ if (f < 90.0 || f > fs * 0.5 - 100.0) continue; if (bw > 400.0 || bw < 10.0) continue; if (n < max_out) { out[n].freq = f; out[n].bw = bw; out[n].mag = mag; ++n; } } /* 5) 按频率升序，前三个即 F1/F2/F3 */ 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; }

四个阈值都有来历。f < 90丢掉的是基频和低频噪声对应的极点，成年男声基频可以低到 80 Hz，再往下就纯是直流残留了。f > fs/2 - 100是抗混叠滤波器的过渡带，那里的能量不可信。bw > 400丢掉的是极宽的极点，这种极点不构成谐振，只是频谱斜坡。bw < 10反过来丢掉的是一种假象：带宽小于 10 Hz 的极点通常来自数值误差而非真实谐振，真实语音的 F1 带宽一般在 50～150 Hz。

排序之后取前三个，就是这一帧的 F1/F2/F3。但要注意：n有可能小于 3，尤其是清音段和静音段。输出时把缺失的位置填nan而不是填 0，否则画出来的轨迹会在清音处突然掉到 0 Hz，看图的人以为是算法崩了。

5. 避坑与排查：共振峰飘、数值发散、结果对不上的六类现场

这套代码写完到能稳定出结果，中间踩的坑比写代码的时间还长。下面六条是我实际遇到过、并且能复现的，按现象—原因—解决整理，你跑的时候大概率会撞上其中三四个。

5.1 现象：共振峰轨迹每帧剧烈跳变，图画出来像噪声

原因通常是三个叠加的。第一，帧长太短导致自相关估计方差大，10 ms 帧长在 16 kHz 下只有 160 点，r[k]的估计误差能到 20% 以上；第二，阶数不够，p = 10时只有 5 对极点，勉强够 F1/F2，F3 时有时无；第三，没有做预加重，低频能量占绝对主导，高阶极点被挤到高频去凑数。

解决顺序是先把帧长拉到 25 ms、阶数拉到 14、确认预加重那一行代码真的执行了（加个printf打到文件里看前 10 个样点的值），再看跳变是否收敛。如果还跳，检查是不是在对每一帧独立排序之后直接输出——相邻帧的 F1/F2 可能互换，需要在帧间做一次连续性约束，把当前帧的三个频率和上一帧按最小距离配对。

5.2 现象：杜宾递推返回的系数出现 1e8 量级，后面求根全变成 nan

直接原因是反射系数k越界没被拦住。追下去有两条路径：一是输入帧是纯静音，r[0]接近 0，除法放大；二是录音里有强烈的周期性干扰，比如电源嗡声，让r[1]、r[2]异常大，第一级k就超过 0.999。

解决是在durbin入口加能量门限，在每级加|k| < 1的硬检查并且返回错误码，调用方收到 -1 时降阶重试而不是继续用垃圾系数。另外把预加重之后的自相关加一个滞后窗，bw取 60 Hz 左右，能压掉大部分由嗡声引起的病态。

5.3 现象：牛顿迭代收敛到的根重复，输出的共振峰数量少于阶数的一半

牛顿法的收敛域不互斥，p个初值里三四个跑到同一个根是常态。没有降阶的话，求出来的有效极点数量会明显偏少，表现为某些帧只有一两个共振峰。

解决就是 4.2 节说的降阶：每求出一个根就除掉，同时把共轭根一起去掉。加了降阶之后，p = 14的帧基本能稳定拿到 4～5 个候选极点，筛完剩 3 个正好。注意降阶之后多项式的系数会累积误差，降到 6 阶以下就该停了，剩下的根用原多项式加"排除已求根附近"的初值策略补。

5.4 现象：算出来的 F1 永远固定在 300 Hz 附近，不管什么元音

这是符号约定错了的典型表现。durbin里k的初值写成+r[i]而不是-r[i]，得到的a整体反号，求根之后所有极点的相位发生镜像，频率排序也跟着乱。因为它"照样能算出数"，很容易被当成正常结果接受。

排查方法很简单：找一段明确的 /a/ 元音，正常结果应该是 F1 在 700～900 Hz、F2 在 1100～1400 Hz。如果 F1 卡在 300 Hz 附近不动，八成是符号问题。改回来之后记得重新验证，因为降阶代码里也隐含了符号假设。

5.5 现象：单帧结果和成熟工具差 100 Hz 以上，但趋势一致

参数没对齐。常见差异点按影响从大到小排：预加重系数（0.95 vs 0.97）、窗函数（汉明 vs 汉宁，汉宁窗的主瓣更宽，会让相邻共振峰糊在一起）、帧长、以及自相关是否加了滞后窗。任何一个不同，绝对频率都会有几十赫兹的偏移。

解决方法不是去猜，而是把中间结果导出来逐级比对。在durbin之后把a[1..p]打成 CSV，用另一套工具算同样的帧，对系数数组。系数对上了频率对不上，那就是求根环节的问题；系数就对不上，往自相关和加窗那边查。这一步能省掉大量盲猜。

5.6 现象：debug 编译结果正常，-O2之后结果全乱

优化器把浮点运算重排了。杜宾递推和霍纳法都是有严格顺序依赖的累加，-ffast-math一开，编译器有权把a + b + c重排成a + (b + c)，误差累积路径变了，临界情况下的k可能从 0.999 变成 1.001，直接触发越界返回。

解决办法是编译选项里不要开-ffast-math，只开-O2，并且给递推和求根函数加volatile或者__attribute__((optimize("no-fast-math")))。gcc 下还可以加-fno-aggressive-loop-optimizations保险。这一条在嵌入式和交叉编译里尤其容易中招，因为工具链的默认选项各不相同。

6. 二次开发：三个改法和一个自检习惯

拿到这份代码之后，绝大多数改造需求落在三个地方，我按改造成本从低到高排。

第一是换窗函数。init_hamming换个公式就行，汉宁窗是0.5 - 0.5·cos(2πi/(N-1))，布莱克曼窗是0.42 - 0.5·cos(2πi/(N-1)) + 0.08·cos(4πi/(N-1))。换完一定要重跑同一段测试音频对比，汉宁窗主瓣宽，相邻共振峰容易合并，通常会让 F2/F3 的平均绝对误差上升 20～40 Hz。想做窗长扫描实验，就把窗函数做成函数指针数组，主循环里遍历，输出文件名带上窗函数名，一次跑完。

第二是调阶数和采样率。阶数不是越高越好，16 kHz 下p = 16已经能稳定拿到 F1～F3，p = 20多出来的极点全是需要筛掉的伪峰，反而增加了筛选阈值调参的工作量。真正需要调阶数的是采样率变化：8 kHz 采样取p = 10，44.1 kHz 采样要取到p = 24才够。改的时候同步改MAX_ORDER，别让栈上的tmp、dc、nr三个数组越界。

第三是加批量评测。最实用的做法是把每帧结果输出成 CSV（frame_idx, time_ms, f1, f2, f3, bw1, n_valid），再用十几行 Python 画轨迹、算帧间跳变率。判断算法稳不稳，看两个数：有效帧占比（n_valid >= 3的帧比例）和 F1 的一阶差分绝对均值。前者低于 85% 说明门限太严或降阶逻辑有问题，后者大于 150 Hz 说明帧间连续性没做好。

# eval_track.py：读 CSV，算两个稳定性指标 import csv, statistics rows = [r for r in csv.DictReader(open("formant.csv")) if r["n_valid"] != "0"] valid_ratio = len(rows) / sum(1 for _ in csv.DictReader(open("formant.csv"))) f1 = [float(r["f1"]) for r in rows if r["f1"] != "nan"] jitter = statistics.mean(abs(a - b) for a, b in zip(f1, f1[1:])) if len(f1) > 1 else 0.0 print(f"valid={valid_ratio:.2%} f1_jitter={jitter:.1f} Hz")

还有一个我后来养成、并且强制执行的验证手段：造测试信号。用两个或三个共振峰做级联合成一段元音，F1=800, F2=1200, F3=2500，采样率 16 kHz，时长 1 秒，存成 wav。任何一次改动——换窗、改阶数、改预加重、换编译器优化级别——都先让这段合成音频过一遍，看输出的三个频率有没有落回 ±50 Hz 内。真实录音没有标准答案，合成信号有。这个自检跑一遍不到一秒，但拦下的问题比你想象的多数倍。

从那以后我每次动 LPC 相关的参数，哪怕只是把预加重从 0.97 改成 0.96，都强制先用合成元音走一遍，确认 F1/F2/F3 还在预期区间，再去碰真实录音。共振峰提取这种东西，出错的时候不会崩，只会安静地给你一堆看起来很正常的错数字。希望帮到你。

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