简介这套 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.9516 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 ms240 点16 kHz 用 25 ms400 点帧移统一取窗长的 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^-kk 从 1 到 pG 是增益。LPC 做的事就是找一组a_k让x[n]用前 p 个样本线性预测的均方误差最小。注意A(z)的根就是H(z)的极点而且是 A(z) 在单位圆内的根——这个约定很重要后面求根时所有的符号都和它绑定。3.2 自相关法为什么共振峰提取几乎不用协方差法求解a_k有两条经典路线。协方差法不加窗直接用帧内数据构造非 Toeplitz 的协方差矩阵精度更高但解出来的极点可能跑到单位圆外对应的带宽是负数。共振峰提取最怕这个一个不稳定的极点会让整个 LPC 谱的形状失真。自相关法先用加窗后的序列算出自相关r[k]矩阵变成 ToeplitzLevinson-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~1010 kHz 取 10~1216 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^-1A(z) 1 - Σ a_k z^-k在 z 0 处有多重极点直接对它求导做 Newton数值上极不稳定。改用Ã(z) z^p · A(z) z^p - a_1 z^(p-1) - ... - a_p根完全不变首项系数归一为 1Horner 法可以一次遍历同时算出函数值和导数值。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 vv 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; } /* 复 Newtonz - 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的那一半。典型的成年人元音男性大致落在下表范围超出很多的频点基本可以怀疑是伪峰共振峰典型频率范围典型带宽说明F1200~1000 Hz50~150 Hz与舌位高低强相关F2700~2500 Hz60~200 Hz与舌位前后强相关F31500~3500 Hz100~300 Hz受卷舌、儿化影响F43000~4500 Hz150~400 Hz常被阶数不足吞掉F54000~5500 Hz200~500 Hz16 kHz 采样才有意义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 突然交换**本文还有配套的精品资源点击获取