★ FEATURED ARTICLE 简介面向语音信号处理学习者与嵌入式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ₖ 求偏导并令其为零就得到一组线性方程Σ(k1..p) aₖ * R(|i-k|) -R(i), i 1..p R(k) Σ(nk..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 28kHz 采样取 12 阶10kHz 取 14 阶16kHz 取 18 到 20 阶。采样率常用阶数理论谐振数实际可用共振峰8 kHz10~125~6F1~F410 kHz12~146~7F1~F416 kHz16~208~10F1~F5阶数不足的后果非常直观F1 和 F2 靠得近时两个极点会被强行拟合进一个宽极点里结果只报出一个峰或者两个频率都被拉到中间。阶数过高的后果更隐蔽——多余的自由度会去拟合基频的谐波结构冒出一堆带宽极窄低于 20Hz的伪极点看起来像共振峰实际上每次基频变化它就换个位置。这就是为什么后面筛选规则里必须有一条带宽下限。3. 前端处理链预加重、汉明窗与自相关3.1 预加重先把 -6dB/oct 的谱倾斜扳平语音频谱有一个天然的下倾趋势低频能量比高频高出一大截大致是每倍频程掉 6dB。这个倾斜来自声门脉冲的谱形和唇端辐射的微分效应。如果直接把它喂给 LPC最小均方误差准则会把绝大部分建模能力花在低频段因为那里的误差绝对值最大高频共振峰基本被放弃。预加重就是一级一阶高通y(n) x(n) - μ·x(n-1)。μ 通常取 0.9 到 0.978kHz 采样常用 0.95采样率越高越接近 1。它恰好抵消了唇端辐射带来的那个零点同时把高频抬起来让各频段的误差量级接近。注意预加重是双刃剑。它在抬高高频语音的同时也抬高了高频噪声信噪比差的录音做完预加重可能更糟。实际调试时如果发现 F3、F4 抖得厉害把 μ 从 0.95 降到 0.9 往往是第一个该试的动作。3.2 汉明窗旁瓣 -43dB 换来的到底是什么分帧之后必须加窗否则帧边界处的突变会在频谱上产生严重的泄漏。矩形窗的主瓣最窄2π/N频率分辨率最好但第一旁瓣只有 -13dB意味着一个强低频成分会把自己的影子铺到很远的高频去。语音的动态范围轻松超过 30dBF1 完全可能把 F4 淹没在泄漏里。汉明窗把第一旁瓣压到 -43dB代价是主瓣展宽到 4π/N频率分辨率减半。这个代价在 LPC 里是可以接受的因为最终频率精度不是由窗的主瓣宽度决定的而是由极点位置决定的——窗只影响自相关序列不影响求根的精度。帧长的选择取决于基音周期。窗内至少要覆盖两到三个基音周期LPC 才能把激励和声道分开否则拟合出来的模型会跟着单个脉冲跑。成人基频 80 到 300Hz对应周期 3 到 12ms所以帧长取 20 到 30ms 比较稳妥8kHz 采样下就是 160 到 240 点常用 200 点配 80 点帧移。还有一个容易被忽略的约束帧长 N 必须远大于阶数 p。一般要求 N ≥ 4p8kHz 下 p12、N200余量充足。如果为了提时间分辨率把帧长压到 80 点还硬用 12 阶自相关矩阵接近奇异杜宾递推中途就会崩。3.3 自相关为什么只算到 p 就够了Yule-Walker 方程里只用到 R(0) 到 R(p)共 p1 个值多算一个都是浪费。直接按定义算是 O(N·p)8kHz 下 N200、p12一次 2400 次乘加嵌入式设备上也毫无压力如果帧长很长比如 1024 点做离线分析走 FFT 转成功率谱再做逆变换会更快复杂度 O(N log N)但短帧反而不划算。/* 一帧语音的前端处理预加重 - 加窗 - 自相关 * x : 输入样本已按帧移切好 * n : 帧长8kHz 下取 20025ms * p : LPC 阶数 * mu : 预加重系数8kHz 常用 0.95 * r : 输出自相关长度 p1 * 返回 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 这 p1 个值 */ 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³)。p12 时是 1728 次运算其实也不多但问题不在速度而在数值条件托普利茨矩阵的条件数可能很差消元过程中的舍入误差会让解出来的滤波器不稳定。Levinson-Durbin 递推换了个思路从一阶开始逐级往上加已经算好了 i-1 阶的最优预测器现在加第 i 阶新的系数只需要在旧系数基础上做一次旋转不需要从头解方程。每一级的中间量 kᵢ 有个明确的物理身份——第 i 级格型滤波器的反射系数。这个身份带来两个好处。第一复杂度降到 O(p²)p12 时只有 78 次乘加。第二|kᵢ| 1 是格型结构稳定的充要条件而自相关矩阵的半正定性保证了递推出来的 kᵢ 必然满足这个条件。所以自相关法 杜宾递推天然产出最小相位滤波器不需要额外检查。另一个可用的中间量是预测误差能量 Eᵢ (1 - kᵢ²)·Eᵢ₋₁。它单调递减如果中途算出 Eᵢ ≤ 0 或者 |kᵢ| ≥ 1说明输入数据出了问题帧太短、定点溢出、自相关不是真正半正定这时候与其硬算下去不如直接丢弃这一帧。/* Levinson-Durbin 递推自相关 - LPC 系数 * r : 自相关r[0..p] * p : 阶数 * a : 输出a[0]1a[1..p] 为预测系数对应 A(z)1sum(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 得到的去掉了 z0 的平凡根剩下的根和 A(z) 完全相同。p12 时这是个 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] 永远是当前帧最低的那个共振峰帧间才有本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
