简介面向语音信号处理学习者的LPC共振峰提取C语言实现工程包基于TI Code Composer Studio 3.3开发环境完整演示线性预测编码方法从分帧加窗、端点检测到杜宾/牛顿算法求解预测系数并提取共振峰的流程适合嵌入式语音处理、语音识别方向的学生与工程师学习参考。压缩包共72个文件约274KB主要包含C源码、CCS工程配置与调试文件.pjt/.out/.map/.obj、数据存储文件.dbf/.fpt/.cdx、说明文档和WAV语音样本。包内按功能划分为db、voicemark、hamming、ndroots等模块便于对照理解各阶段实现。已有383人学习下载。读者可依据readme和源码在CCS3.3中直接编译运行观察不同语音样本的共振峰提取结果尤其适合借助杜宾算法迭代与牛顿算法加速收敛的细节展开深入研究。从元音识别到发声评测用LPC在C语言里稳定提取共振峰手头有个语音评测的小项目需求是把这个元音发得标不标准变成机器可判断的指标。翻了一圈方案最后还是老老实实回来用LPC线性预测编码做共振峰提取。共振峰Formant说白了就是语音谱上那几个能量集中的峰值区域对应声道谐振的特性元音识别、声纹分析、发音质量评估全靠它。基于帧的LPC分析只需要乘加运算计算量可控非常适合在C语言环境下做嵌入式或离线批处理。本文从LPC的原理、C语言算法实现到工程调参完整走一遍适合正在做语音信号处理、想绕开MATLAB直接用C/C落地算法的读者也适合毕业设计里需要真跑通而不是只仿真的同学。LPC的想法很直白语音信号的当前采样点可以用过去若干个采样点的线性组合来近似组合系数就是我们要解的LPC系数。而这个线性预测器本质上是声道传输函数的建模共振峰就藏在这个模型里。整条处理链路是原始音频 → 预加重 → 分帧加窗 → 自相关序列 → 莱文森-杜宾递推求LPC系数 → 共振峰频率提取 → 后处理平滑下面逐个环节拆开了讲。1. 为什么是LPC共振峰提取的选型逻辑共振峰的物理意义要先捋清楚。发一个元音时声带产生周期性的脉冲激励这个激励经过口腔、鼻腔这些腔体的调谐就形成了一系列谐振峰。这些谐振峰在频谱上体现为局部极大值从低频到高频分别叫F1、F2、F3。元音之所以不同主要就是F1和F2的相对位置不一样比如/a/的F1一般在700-1000HzF2在1100-1400Hz附近而/i/的F1低到300Hz左右F2却冲到2200Hz以上。提取这些峰值主要有三条路Fourier频谱直接找峰、倒谱法、LPC包络法。简单频谱找峰的问题是语音信号里不仅有声道谐振的包络还叠着声带激励的精细谐波结构谱峰里混着一堆假峰直接找极值点容易被谐波干扰。倒谱法能把激励和声道分开但计算复杂而且对基频较高的女声和孩子声音容易产生同态分析混叠。LPC的优势在于它是参数化建模全极点模型直接模拟声道传输特性系数本身就是声道的指纹。一段语音帧用LPC系数表征之后共振峰可以通过两种方式读出来——一种是计算预测误差滤波器的根另一种是对LPC谱做峰值搜索。相比直接FFTLPC谱非常平滑匹配的都是声道的整体谐振特性对基频谐波不敏感稳定性好得多。另外一个实际考虑是计算量。共振峰提取要逐帧处理每一帧都要做自相关和递推求解。对于一个30ms的帧10kHz采样率就是300个点LPC阶数12到14阶莱文森-杜宾递推的复杂度大约是O(p²)p是阶数也就是一百多次乘加运算这在单片机上都是可以接受的。2. 预处理三件套预加重、分帧加窗的参数怎么配很多人一上来就急着写递推忽略预处理结果提取出的共振峰乱跳。我在这上面栽过跟头这里把每个环节的参数和理由都交代清楚。2.1 预加重先把高频能量捞回来语音信号的频谱大致按每倍频程6dB衰减高频段能量太弱如果直接分析共振峰检测会偏向低频导致高频共振峰特别是F3以上估计不准确。预加重是一阶高通滤波器y[n] x[n] - α·x[n-1]α一般取0.95到0.97之间。这段代码很简单static double pre_emphasis(double *sig, int n, double alpha) { double prev 0.0; for (int i 0; i n; i) { double cur sig[i]; sig[i] cur - alpha * prev; prev cur; } return 0.0; }实测下来α取0.95时共振峰曲线最平滑取太高会让高频的噪声同时被放大。另外一个细节做共振峰提取时预加重只是分析用不需要做去加重还原因为我们只关心频率位置不关心频谱绝对幅值。2.2 分帧帧长帧移的选择直接影响时间分辨率语音是短时平稳的一般10-30ms内可以认为声道形状基本不变。帧长选多少取决于你关心共振峰的时间变化速度。帧太长共振峰变化被平均掉过渡音比如双元音的轨迹会失真帧太短频率分辨率下降LPC谱不够细腻。我的经验配比采样率Fs 10kHz → 帧长N 240点24ms帧移 80点采样率Fs 16kHz → 帧长N 400点25ms帧移 160点有一个必须注意的点一次LPC分析至少需要覆盖2-3个基音周期。成年男性基频约100-120Hz周期约8-10ms女声约200-220Hz周期约4-5ms。如果帧长只有10ms对于低频男声一帧里只有1个周期LPC预测会偏向激励源而不是声道模型结果直接崩掉。所以帧长宁可偏长也不要为了追求时间分辨率而牺牲稳定性。2.3 加窗汉明窗是默认选择分帧后直接算自相关会带来截断效应频谱出现旁瓣泄漏所以每帧要乘上一个窗函数。汉明窗是共振峰提取里最稳妥的选择公式w[n] 0.54 - 0.46·cos(2πn/(N-1)), 0 ≤ n N汉明窗的好处是主瓣宽度和旁瓣衰减平衡得好。矩形窗旁瓣太高不推荐布莱克曼窗旁瓣低但主瓣宽会把相近的共振峰比如F2和F3靠得近时抹在一块儿。代码static void apply_hamming(double *frame, int n) { for (int i 0; i n; i) { frame[i] * 0.54 - 0.46 * cos(2.0 * M_PI * i / (n - 1)); } }预处理三个环节做完接下来就是核心的LPC系数求解。3. 莱文森-杜宾递推LPC系数求解的C语言核心实现LPC的核心假设是当前信号值可以用前p个值的线性组合近似。s[n] ≈ a[1]·s[n-1] a[2]·s[n-2] ... a[p]·s[n-p]误差是e[n] s[n] - Σ(a[k]·s[n-k])我们要找一组a[k]让误差的能量最小。用自相关法求解最终会落到一个托普利兹Toeplitz方程组R[0]R[1]...R[p-1] [a[1]] [R[1]] R[1]R[0]...R[p-2] [a[2]] [R[2]] ... [ ... ] [ ...] R[p-1]R[p-2]...R[0] [a[p]] [R[p]]R[k]是信号的自相关函数。托普利兹矩阵的特殊结构决定了它可以用莱文森-杜宾Levinson-Durbin递推高效求解复杂度O(p²)比常规高斯消元快得多。3.1 自相关函数计算自相关序列要在整帧范围内做R[k] Σ(s[n]·s[n-k]), k 0, 1, ..., p注意这里s[n]是已经加过窗的帧。一个工程上的关键问题计算R[k]时需要把帧外数据当零还是直接截断自相关法的标准做法是把窗外信号视为零所以R[k]的求和范围是n k到n N-1而不是从0到N-1-k。static void autocorrelation(double *frame, int n, double *r, int p) { for (int k 0; k p; k) { double sum 0.0; for (int i k; i n; i) { sum frame[i] * frame[i - k]; } r[k] sum; } }这里有一个别偷懒的理由如果R[0]非常小接近静音帧后面的递推会除以接近零的数导致数值爆炸。后面我会讲怎么处理。3.2 莱文森-杜宾递推实现我把完整代码贴出来这个可以直接抄typedef struct { double a[LPC_MAX_ORDER 1]; // a[0]恒为1.0预测误差滤波器系数 double e; // 预测误差能量 int order; // 实际阶数 } LPCResult; // 返回0表示成功-1表示数值异常如静音帧 static int levinson_durbin(double *r, int p, LPCResult *lpc) { if (r[0] 1e-10) { lpc-e 0.0; return -1; } double *a lpc-a; double *b (double *)calloc(p 1, sizeof(double)); if (!b) return -1; a[0] 1.0; double error r[0]; lpc-e error; for (int i 1; i p; i) { // 计算反射系数 k_i double sum 0.0; for (int j 1; j i; j) { sum a[j] * r[i - j]; } double ki (r[i] - sum) / error; // 更新预测系数复制旧系数到临时数组 memcpy(b, a, (i 1) * sizeof(double)); for (int j 1; j i; j) { a[j] b[j] - ki * b[i - j]; } a[i] -ki; // 注意符号约定 // 更新误差能量 error * (1.0 - ki * ki); if (error 1e-10) { // 误差能量太小继续递推只会放大数值误差提前终止 free(b); lpc-e error; for (int j i 1; j p; j) a[j] 0.0; return 0; } } lpc-e error; free(b); return 0; }几点说明反射系数ki的绝对值理论上一定小于1但浮点误差可能导致略微越界。如果|ki| 1说明数值不稳定或输入不是平稳信号这时应该强制截断并返回异常而不是让它继续跑下去。符号约定要统一我这里的a[i]是预测误差滤波器的系数构造的是A(z) 1 - Σ(a[k]·z^(-k))后面求共振峰时会用到这个约定。如果你在别的代码里看到符号相反不用慌只是约定不同。误差能量error每步都会乘以(1 - ki²)整体是递减的。如果原始帧能量很小error会迅速跌到接近零。3.3 数值稳定性的两个硬性要求第一变量一律用double。float在自相关累加里会严重损失精度特别是帧长256点以上时累加误差会被莱文森递推放大导致a[k]乱跳。我在STM32F4上试过用float跑同样的输入共振峰结果能偏出100Hz。第二对于R[0]特别小的帧直接跳过。语音里有清音、停顿、静音段这些帧的R[0]可能比正常帧低几个数量级强行递推会算出巨大的无用系数。工程上我设定一个阈值如果帧能量即R[0]低于本段语音最大帧能量的1%就标记为无效帧共振峰输出置为-1或NaN由上层逻辑处理。4. 共振峰读取求根法和峰值搜索法的工程取舍拿到LPC系数a[1..p]之后模型是H(z) 1 / (1 - Σ a[k]·z^(-k))共振峰对应的是分母多项式的复共轭根的极角位置。求出A(z)的根保留上半平面的根极角θ对应的频率是f θ·Fs/(2π)。另一个更省事的方法是直接对A(z)在单位圆上采样得到LPC谱包络然后找局部极大值。两种方法的对比方法优点缺点适用场景多项式的根精确共振峰频率和带宽都能拿到需要复根求解如牛顿法实现复杂需要带宽信息、离线高精度分析LPC谱峰值搜索实现简单FFT结果直接可视化频率分辨率受FFT点数限制只能拿频率拿不到带宽实时处理、嵌入式部署我实际项目里两个都写过结论是如果你只需要F1、F2、F3的中心频率峰值搜索法完全够用而且好调试。求根法适合学术论文级别的分析或者需要共振峰带宽的场景。4.1 峰值搜索法实现要点在单位圆上对A(z)做DFT就得LPC频谱的倒数LPC频谱包络就是1/|A(e^jω)|²。具体步骤对A(z)的系数做N_fft点FFT得到A(k)。求功率谱包络P(k) 1 / |A(k)|²注意幅度归一化不需要太讲究我们只找峰值位置。在感兴趣的频段比如200Hz到Fs/2-500Hz内找局部极大值。对峰值做抛物线插值把频率分辨率从Fs/N_fft细化到真实值。如果不用FFT也可以直接用Goertzel算法在候选频率点上算|A(e^jω)|²省内存但不如FFT直观。抛物线插值公式Δ 0.5·(P[k-1] - P[k1]) / (P[k-1] - 2·P[k] P[k1]) f_peak (k Δ)·Fs / N_fft别小看这一步插值。如果你只用离散FFT的峰值索引直接读频率当N_fft512、Fs10k时频率分辨率接近20Hz共振峰相邻帧之间会看到明显的台阶状跳变。插值之后平滑多了也是后面共振峰轨迹连续的保障之一。4.2 求根法实现要点求A(z)的复根我用最简单的牛顿法// A(z) 1 - a[1]z^-1 - a[2]z^-2 - ... - a[p]z^-p // 多项式按z的正幂次重写: z^p - a[1]z^(p-1) - ... - a[p] 0用初始猜测e^(j·2π·k/p)迭代k取0到p-1每个根收敛后用多项式降阶合成除法去除已找到的根。这个方法对p14以内很稳定但是要处理根在迭代中漂移的情况。实际工程里我很少用求根法除非发论文需要精确带宽。5. 实测效果与参数调优阶数、帧长、采样率的配合方式5.1 LPC阶数怎么定经验公式p Fs/1000 2 到 Fs/1000 4 之间。10kHz采样率用12到14阶16kHz用14到16阶。阶数太低LPC谱包络太粗相邻共振峰分不开阶数太高除了声道谐振之外会把基频谐波也建模进去产生假共振峰。我实测过一组数8kHz采样率的语音用8阶LPC提取元音/a/的共振峰F1位置还能看F2直接和F1黏在一起完全分不开提到14阶后F1、F2、F3都清晰了。继续加到24阶反而在500Hz附近额外冒出一个假峰就是谐波被建模了。5.2 一个典型的完整提取流程代码把前面的片段整合成一个可运行的提取函数#define Fs 16000 #define FRAME_LEN 400 // 25ms #define FRAME_SHIFT 160 // 10ms #define LPC_ORDER 16 typedef struct { double f1, f2, f3; // 实际读到的共振峰频率 int valid; // 1有效 0无效 } FormantFrame; static double *g_fftOut; int extract_formant_frame(short *pcm, int n, FormantFrame *out) { // 1. 转浮点预加重 static double buf[FRAME_LEN]; for (int i 0; i FRAME_LEN i n; i) { buf[i] (double)pcm[i] / 32768.0; } pre_emphasis(buf, FRAME_LEN, 0.95); // 2. 加窗 apply_hamming(buf, FRAME_LEN); // 3. 自相关 莱文森-杜宾 double r[LPC_ORDER 1]; autocorrelation(buf, FRAME_LEN, r, LPC_ORDER); LPCResult lpc; if (levinson_durbin(r, LPC_ORDER, lpc) ! 0) { out-valid 0; return -1; } // 4. LPC谱包络FFT方式 // 用256点FFT计算P(k) 1 / |A(k)|^2参考4.1节 // 5. 峰值搜索插值输出前三个共振峰 // 搜索范围: F1在200-1000HzF2在800-2500HzF3在1500-3500Hz // 具体实现省略就是分段找局部最大值 out-valid 1; return 0; }这里面每一步都有明确的输出可以单独打印出来调试验证。我建议你在本机先用一段录好的元音音频跑一遍把每一步的自相关数组、LPC系数、LPC谱都打印出来和MATLAB的lpc函数结果对一下确认无误再往下做嵌入式移植。5.3 共振峰轨迹平滑即便预处理和插值都做了相邻帧的共振峰偶尔还是会跳变。最常见的是在元音和辅音的过渡段LPC模型会突然把F1和F2交换位置——因为它们俩在那个瞬间靠得太近峰值搜索会选错顺序。平滑手段有两种一是对输出频率做中值滤波窗口3到5帧二是对LPC系数本身做一阶低通平滑a_smooth[n] β·a_smooth[n-1] (1-β)·a_raw[n]β取值0.5到0.7β太大会让共振峰轨迹滞后严重。我在已经标注好元音段落的评测任务里直接对输出频率做5点中值滤波就够用了效果明显但不迟钝。6. 我在工程里踩过的坑共振峰跳变、静音帧误判和边界效应这部分不是泛泛的经验之谈是我踩过的坑记录按影响从大到小排。6.1 共振峰交换问题说白了就是不同共振峰轨迹在某个时间点交叉或靠得极近时逐帧独立搜索算法会跟错人。这个问题在合成语音里非常明显自然语音的过渡段也会有。我的处理方式是给峰值搜索加上跟踪约束初始化时按频率把前三个峰标为F1、F2、F3后续帧搜索时以上一帧位置为中心开搜索窗。比如这一帧的F1搜索范围是上一帧F1±250Hz这个假设基于声道变化不会跳变的物理约束。加了跟踪约束之后跳变率从5%左右降到0.5%以内。6.2 静音帧的R[0]阈值最开始我只判断R[0]是否为0结果安静环境下一段几秒的录音有些帧能量是正常帧的百分之一LPC系数算出来非常夸张共振峰直接输出到几千赫兹的假峰。后来加了相对能量判断当前帧能量 / 本段最大帧能量 0.01直接跳过。这个阈值在安静室内录音环境表现很好噪声大的环境要适当提高到0.02。6.3 首尾帧的自相关窗长第一帧虽然没有前一帧可参考但它的预测误差会偏大因为窗内数据不足。我在处理时会在音频开头补一段50ms的静音全零让前面几帧自然消耗掉再把输出对齐回来。这种做法虽然浪费几十毫秒但避免了开头几帧共振峰的异常波动。6.4 浮点误差累积如果你把 pre_emphasis 和 hamming 合在一起做成函数每次运算都写回原数组长时间运行会看到结果逐渐漂移。我最后改成了在每帧开头的局部double数组里完成所有中间计算只有最终的共振峰结果回传外部。这样既保证了可重入性也让每个帧之间是完全独立的。总结LPC共振峰提取不是什么高深算法但真正做到工程级稳定需要抠的细节一点不少预处理参数配比、自相关数值范围、递推的异常终止、峰值搜索的顺序跟踪、以及静音帧的阈值判断。就我实际项目体验而言这套C语言实现在音频评测、发声矫正设备里跑几十个小时没有出过明显问题处理一帧16kHz的语音大约在零点几毫秒级别不含FFT优化性能完全不是瓶颈。最后留一个小技巧调试共振峰提取时不要只看数字把提取结果叠加到LPC谱包络上画出来检查。代码里临时输出一张PPM格式的灰度图或者直接打印ASCII波形到终端比盯着printf的几百个数字直观得多。语音处理是用耳朵和眼睛一起调的活数字跑出来再准也得先看着像才算数。本文还有配套的精品资源点击获取 用LPC在C语言中实现共振峰提取：从元音识别到发声评测
简介面向语音信号处理学习者的LPC共振峰提取C语言实现工程包基于TI Code Composer Studio 3.3开发环境完整演示线性预测编码方法从分帧加窗、端点检测到杜宾/牛顿算法求解预测系数并提取共振峰的流程适合嵌入式语音处理、语音识别方向的学生与工程师学习参考。压缩包共72个文件约274KB主要包含C源码、CCS工程配置与调试文件.pjt/.out/.map/.obj、数据存储文件.dbf/.fpt/.cdx、说明文档和WAV语音样本。包内按功能划分为db、voicemark、hamming、ndroots等模块便于对照理解各阶段实现。已有383人学习下载。读者可依据readme和源码在CCS3.3中直接编译运行观察不同语音样本的共振峰提取结果尤其适合借助杜宾算法迭代与牛顿算法加速收敛的细节展开深入研究。从元音识别到发声评测用LPC在C语言里稳定提取共振峰手头有个语音评测的小项目需求是把这个元音发得标不标准变成机器可判断的指标。翻了一圈方案最后还是老老实实回来用LPC线性预测编码做共振峰提取。共振峰Formant说白了就是语音谱上那几个能量集中的峰值区域对应声道谐振的特性元音识别、声纹分析、发音质量评估全靠它。基于帧的LPC分析只需要乘加运算计算量可控非常适合在C语言环境下做嵌入式或离线批处理。本文从LPC的原理、C语言算法实现到工程调参完整走一遍适合正在做语音信号处理、想绕开MATLAB直接用C/C落地算法的读者也适合毕业设计里需要真跑通而不是只仿真的同学。LPC的想法很直白语音信号的当前采样点可以用过去若干个采样点的线性组合来近似组合系数就是我们要解的LPC系数。而这个线性预测器本质上是声道传输函数的建模共振峰就藏在这个模型里。整条处理链路是原始音频 → 预加重 → 分帧加窗 → 自相关序列 → 莱文森-杜宾递推求LPC系数 → 共振峰频率提取 → 后处理平滑下面逐个环节拆开了讲。1. 为什么是LPC共振峰提取的选型逻辑共振峰的物理意义要先捋清楚。发一个元音时声带产生周期性的脉冲激励这个激励经过口腔、鼻腔这些腔体的调谐就形成了一系列谐振峰。这些谐振峰在频谱上体现为局部极大值从低频到高频分别叫F1、F2、F3。元音之所以不同主要就是F1和F2的相对位置不一样比如/a/的F1一般在700-1000HzF2在1100-1400Hz附近而/i/的F1低到300Hz左右F2却冲到2200Hz以上。提取这些峰值主要有三条路Fourier频谱直接找峰、倒谱法、LPC包络法。简单频谱找峰的问题是语音信号里不仅有声道谐振的包络还叠着声带激励的精细谐波结构谱峰里混着一堆假峰直接找极值点容易被谐波干扰。倒谱法能把激励和声道分开但计算复杂而且对基频较高的女声和孩子声音容易产生同态分析混叠。LPC的优势在于它是参数化建模全极点模型直接模拟声道传输特性系数本身就是声道的指纹。一段语音帧用LPC系数表征之后共振峰可以通过两种方式读出来——一种是计算预测误差滤波器的根另一种是对LPC谱做峰值搜索。相比直接FFTLPC谱非常平滑匹配的都是声道的整体谐振特性对基频谐波不敏感稳定性好得多。另外一个实际考虑是计算量。共振峰提取要逐帧处理每一帧都要做自相关和递推求解。对于一个30ms的帧10kHz采样率就是300个点LPC阶数12到14阶莱文森-杜宾递推的复杂度大约是O(p²)p是阶数也就是一百多次乘加运算这在单片机上都是可以接受的。2. 预处理三件套预加重、分帧加窗的参数怎么配很多人一上来就急着写递推忽略预处理结果提取出的共振峰乱跳。我在这上面栽过跟头这里把每个环节的参数和理由都交代清楚。2.1 预加重先把高频能量捞回来语音信号的频谱大致按每倍频程6dB衰减高频段能量太弱如果直接分析共振峰检测会偏向低频导致高频共振峰特别是F3以上估计不准确。预加重是一阶高通滤波器y[n] x[n] - α·x[n-1]α一般取0.95到0.97之间。这段代码很简单static double pre_emphasis(double *sig, int n, double alpha) { double prev 0.0; for (int i 0; i n; i) { double cur sig[i]; sig[i] cur - alpha * prev; prev cur; } return 0.0; }实测下来α取0.95时共振峰曲线最平滑取太高会让高频的噪声同时被放大。另外一个细节做共振峰提取时预加重只是分析用不需要做去加重还原因为我们只关心频率位置不关心频谱绝对幅值。2.2 分帧帧长帧移的选择直接影响时间分辨率语音是短时平稳的一般10-30ms内可以认为声道形状基本不变。帧长选多少取决于你关心共振峰的时间变化速度。帧太长共振峰变化被平均掉过渡音比如双元音的轨迹会失真帧太短频率分辨率下降LPC谱不够细腻。我的经验配比采样率Fs 10kHz → 帧长N 240点24ms帧移 80点采样率Fs 16kHz → 帧长N 400点25ms帧移 160点有一个必须注意的点一次LPC分析至少需要覆盖2-3个基音周期。成年男性基频约100-120Hz周期约8-10ms女声约200-220Hz周期约4-5ms。如果帧长只有10ms对于低频男声一帧里只有1个周期LPC预测会偏向激励源而不是声道模型结果直接崩掉。所以帧长宁可偏长也不要为了追求时间分辨率而牺牲稳定性。2.3 加窗汉明窗是默认选择分帧后直接算自相关会带来截断效应频谱出现旁瓣泄漏所以每帧要乘上一个窗函数。汉明窗是共振峰提取里最稳妥的选择公式w[n] 0.54 - 0.46·cos(2πn/(N-1)), 0 ≤ n N汉明窗的好处是主瓣宽度和旁瓣衰减平衡得好。矩形窗旁瓣太高不推荐布莱克曼窗旁瓣低但主瓣宽会把相近的共振峰比如F2和F3靠得近时抹在一块儿。代码static void apply_hamming(double *frame, int n) { for (int i 0; i n; i) { frame[i] * 0.54 - 0.46 * cos(2.0 * M_PI * i / (n - 1)); } }预处理三个环节做完接下来就是核心的LPC系数求解。3. 莱文森-杜宾递推LPC系数求解的C语言核心实现LPC的核心假设是当前信号值可以用前p个值的线性组合近似。s[n] ≈ a[1]·s[n-1] a[2]·s[n-2] ... a[p]·s[n-p]误差是e[n] s[n] - Σ(a[k]·s[n-k])我们要找一组a[k]让误差的能量最小。用自相关法求解最终会落到一个托普利兹Toeplitz方程组R[0]R[1]...R[p-1] [a[1]] [R[1]] R[1]R[0]...R[p-2] [a[2]] [R[2]] ... [ ... ] [ ...] R[p-1]R[p-2]...R[0] [a[p]] [R[p]]R[k]是信号的自相关函数。托普利兹矩阵的特殊结构决定了它可以用莱文森-杜宾Levinson-Durbin递推高效求解复杂度O(p²)比常规高斯消元快得多。3.1 自相关函数计算自相关序列要在整帧范围内做R[k] Σ(s[n]·s[n-k]), k 0, 1, ..., p注意这里s[n]是已经加过窗的帧。一个工程上的关键问题计算R[k]时需要把帧外数据当零还是直接截断自相关法的标准做法是把窗外信号视为零所以R[k]的求和范围是n k到n N-1而不是从0到N-1-k。static void autocorrelation(double *frame, int n, double *r, int p) { for (int k 0; k p; k) { double sum 0.0; for (int i k; i n; i) { sum frame[i] * frame[i - k]; } r[k] sum; } }这里有一个别偷懒的理由如果R[0]非常小接近静音帧后面的递推会除以接近零的数导致数值爆炸。后面我会讲怎么处理。3.2 莱文森-杜宾递推实现我把完整代码贴出来这个可以直接抄typedef struct { double a[LPC_MAX_ORDER 1]; // a[0]恒为1.0预测误差滤波器系数 double e; // 预测误差能量 int order; // 实际阶数 } LPCResult; // 返回0表示成功-1表示数值异常如静音帧 static int levinson_durbin(double *r, int p, LPCResult *lpc) { if (r[0] 1e-10) { lpc-e 0.0; return -1; } double *a lpc-a; double *b (double *)calloc(p 1, sizeof(double)); if (!b) return -1; a[0] 1.0; double error r[0]; lpc-e error; for (int i 1; i p; i) { // 计算反射系数 k_i double sum 0.0; for (int j 1; j i; j) { sum a[j] * r[i - j]; } double ki (r[i] - sum) / error; // 更新预测系数复制旧系数到临时数组 memcpy(b, a, (i 1) * sizeof(double)); for (int j 1; j i; j) { a[j] b[j] - ki * b[i - j]; } a[i] -ki; // 注意符号约定 // 更新误差能量 error * (1.0 - ki * ki); if (error 1e-10) { // 误差能量太小继续递推只会放大数值误差提前终止 free(b); lpc-e error; for (int j i 1; j p; j) a[j] 0.0; return 0; } } lpc-e error; free(b); return 0; }几点说明反射系数ki的绝对值理论上一定小于1但浮点误差可能导致略微越界。如果|ki| 1说明数值不稳定或输入不是平稳信号这时应该强制截断并返回异常而不是让它继续跑下去。符号约定要统一我这里的a[i]是预测误差滤波器的系数构造的是A(z) 1 - Σ(a[k]·z^(-k))后面求共振峰时会用到这个约定。如果你在别的代码里看到符号相反不用慌只是约定不同。误差能量error每步都会乘以(1 - ki²)整体是递减的。如果原始帧能量很小error会迅速跌到接近零。3.3 数值稳定性的两个硬性要求第一变量一律用double。float在自相关累加里会严重损失精度特别是帧长256点以上时累加误差会被莱文森递推放大导致a[k]乱跳。我在STM32F4上试过用float跑同样的输入共振峰结果能偏出100Hz。第二对于R[0]特别小的帧直接跳过。语音里有清音、停顿、静音段这些帧的R[0]可能比正常帧低几个数量级强行递推会算出巨大的无用系数。工程上我设定一个阈值如果帧能量即R[0]低于本段语音最大帧能量的1%就标记为无效帧共振峰输出置为-1或NaN由上层逻辑处理。4. 共振峰读取求根法和峰值搜索法的工程取舍拿到LPC系数a[1..p]之后模型是H(z) 1 / (1 - Σ a[k]·z^(-k))共振峰对应的是分母多项式的复共轭根的极角位置。求出A(z)的根保留上半平面的根极角θ对应的频率是f θ·Fs/(2π)。另一个更省事的方法是直接对A(z)在单位圆上采样得到LPC谱包络然后找局部极大值。两种方法的对比方法优点缺点适用场景多项式的根精确共振峰频率和带宽都能拿到需要复根求解如牛顿法实现复杂需要带宽信息、离线高精度分析LPC谱峰值搜索实现简单FFT结果直接可视化频率分辨率受FFT点数限制只能拿频率拿不到带宽实时处理、嵌入式部署我实际项目里两个都写过结论是如果你只需要F1、F2、F3的中心频率峰值搜索法完全够用而且好调试。求根法适合学术论文级别的分析或者需要共振峰带宽的场景。4.1 峰值搜索法实现要点在单位圆上对A(z)做DFT就得LPC频谱的倒数LPC频谱包络就是1/|A(e^jω)|²。具体步骤对A(z)的系数做N_fft点FFT得到A(k)。求功率谱包络P(k) 1 / |A(k)|²注意幅度归一化不需要太讲究我们只找峰值位置。在感兴趣的频段比如200Hz到Fs/2-500Hz内找局部极大值。对峰值做抛物线插值把频率分辨率从Fs/N_fft细化到真实值。如果不用FFT也可以直接用Goertzel算法在候选频率点上算|A(e^jω)|²省内存但不如FFT直观。抛物线插值公式Δ 0.5·(P[k-1] - P[k1]) / (P[k-1] - 2·P[k] P[k1]) f_peak (k Δ)·Fs / N_fft别小看这一步插值。如果你只用离散FFT的峰值索引直接读频率当N_fft512、Fs10k时频率分辨率接近20Hz共振峰相邻帧之间会看到明显的台阶状跳变。插值之后平滑多了也是后面共振峰轨迹连续的保障之一。4.2 求根法实现要点求A(z)的复根我用最简单的牛顿法// A(z) 1 - a[1]z^-1 - a[2]z^-2 - ... - a[p]z^-p // 多项式按z的正幂次重写: z^p - a[1]z^(p-1) - ... - a[p] 0用初始猜测e^(j·2π·k/p)迭代k取0到p-1每个根收敛后用多项式降阶合成除法去除已找到的根。这个方法对p14以内很稳定但是要处理根在迭代中漂移的情况。实际工程里我很少用求根法除非发论文需要精确带宽。5. 实测效果与参数调优阶数、帧长、采样率的配合方式5.1 LPC阶数怎么定经验公式p Fs/1000 2 到 Fs/1000 4 之间。10kHz采样率用12到14阶16kHz用14到16阶。阶数太低LPC谱包络太粗相邻共振峰分不开阶数太高除了声道谐振之外会把基频谐波也建模进去产生假共振峰。我实测过一组数8kHz采样率的语音用8阶LPC提取元音/a/的共振峰F1位置还能看F2直接和F1黏在一起完全分不开提到14阶后F1、F2、F3都清晰了。继续加到24阶反而在500Hz附近额外冒出一个假峰就是谐波被建模了。5.2 一个典型的完整提取流程代码把前面的片段整合成一个可运行的提取函数#define Fs 16000 #define FRAME_LEN 400 // 25ms #define FRAME_SHIFT 160 // 10ms #define LPC_ORDER 16 typedef struct { double f1, f2, f3; // 实际读到的共振峰频率 int valid; // 1有效 0无效 } FormantFrame; static double *g_fftOut; int extract_formant_frame(short *pcm, int n, FormantFrame *out) { // 1. 转浮点预加重 static double buf[FRAME_LEN]; for (int i 0; i FRAME_LEN i n; i) { buf[i] (double)pcm[i] / 32768.0; } pre_emphasis(buf, FRAME_LEN, 0.95); // 2. 加窗 apply_hamming(buf, FRAME_LEN); // 3. 自相关 莱文森-杜宾 double r[LPC_ORDER 1]; autocorrelation(buf, FRAME_LEN, r, LPC_ORDER); LPCResult lpc; if (levinson_durbin(r, LPC_ORDER, lpc) ! 0) { out-valid 0; return -1; } // 4. LPC谱包络FFT方式 // 用256点FFT计算P(k) 1 / |A(k)|^2参考4.1节 // 5. 峰值搜索插值输出前三个共振峰 // 搜索范围: F1在200-1000HzF2在800-2500HzF3在1500-3500Hz // 具体实现省略就是分段找局部最大值 out-valid 1; return 0; }这里面每一步都有明确的输出可以单独打印出来调试验证。我建议你在本机先用一段录好的元音音频跑一遍把每一步的自相关数组、LPC系数、LPC谱都打印出来和MATLAB的lpc函数结果对一下确认无误再往下做嵌入式移植。5.3 共振峰轨迹平滑即便预处理和插值都做了相邻帧的共振峰偶尔还是会跳变。最常见的是在元音和辅音的过渡段LPC模型会突然把F1和F2交换位置——因为它们俩在那个瞬间靠得太近峰值搜索会选错顺序。平滑手段有两种一是对输出频率做中值滤波窗口3到5帧二是对LPC系数本身做一阶低通平滑a_smooth[n] β·a_smooth[n-1] (1-β)·a_raw[n]β取值0.5到0.7β太大会让共振峰轨迹滞后严重。我在已经标注好元音段落的评测任务里直接对输出频率做5点中值滤波就够用了效果明显但不迟钝。6. 我在工程里踩过的坑共振峰跳变、静音帧误判和边界效应这部分不是泛泛的经验之谈是我踩过的坑记录按影响从大到小排。6.1 共振峰交换问题说白了就是不同共振峰轨迹在某个时间点交叉或靠得极近时逐帧独立搜索算法会跟错人。这个问题在合成语音里非常明显自然语音的过渡段也会有。我的处理方式是给峰值搜索加上跟踪约束初始化时按频率把前三个峰标为F1、F2、F3后续帧搜索时以上一帧位置为中心开搜索窗。比如这一帧的F1搜索范围是上一帧F1±250Hz这个假设基于声道变化不会跳变的物理约束。加了跟踪约束之后跳变率从5%左右降到0.5%以内。6.2 静音帧的R[0]阈值最开始我只判断R[0]是否为0结果安静环境下一段几秒的录音有些帧能量是正常帧的百分之一LPC系数算出来非常夸张共振峰直接输出到几千赫兹的假峰。后来加了相对能量判断当前帧能量 / 本段最大帧能量 0.01直接跳过。这个阈值在安静室内录音环境表现很好噪声大的环境要适当提高到0.02。6.3 首尾帧的自相关窗长第一帧虽然没有前一帧可参考但它的预测误差会偏大因为窗内数据不足。我在处理时会在音频开头补一段50ms的静音全零让前面几帧自然消耗掉再把输出对齐回来。这种做法虽然浪费几十毫秒但避免了开头几帧共振峰的异常波动。6.4 浮点误差累积如果你把 pre_emphasis 和 hamming 合在一起做成函数每次运算都写回原数组长时间运行会看到结果逐渐漂移。我最后改成了在每帧开头的局部double数组里完成所有中间计算只有最终的共振峰结果回传外部。这样既保证了可重入性也让每个帧之间是完全独立的。总结LPC共振峰提取不是什么高深算法但真正做到工程级稳定需要抠的细节一点不少预处理参数配比、自相关数值范围、递推的异常终止、峰值搜索的顺序跟踪、以及静音帧的阈值判断。就我实际项目体验而言这套C语言实现在音频评测、发声矫正设备里跑几十个小时没有出过明显问题处理一帧16kHz的语音大约在零点几毫秒级别不含FFT优化性能完全不是瓶颈。最后留一个小技巧调试共振峰提取时不要只看数字把提取结果叠加到LPC谱包络上画出来检查。代码里临时输出一张PPM格式的灰度图或者直接打印ASCII波形到终端比盯着printf的几百个数字直观得多。语音处理是用耳朵和眼睛一起调的活数字跑出来再准也得先看着像才算数。本文还有配套的精品资源点击获取