LPC语音共振峰提取C语言实现与参数调优

简介面向语音信号处理学习者与嵌入式音频开发者的LPC共振峰提取C语言实现包围绕线性预测编码的建模流程将语音分帧、汉明窗平滑、端点检测、杜宾递推与牛顿迭代求解预测系数、求根及共振峰估计串联为可编译的工程示例适合课程设计、算法验证与语音识别前端研究参考。压缩包共72个文件、约274KB以C源码、CCS工程与链接配置、编译输出及调试日志为主辅以数据库/索引类中间文件和一个wav语音样本便于对照工程结构理解各算法模块的调用关系与运行结果。现有385人学习可为LPC入门与共振峰提取实践提供较完整的代码骨架和排错线索读者可从中获取从语音样本预处理到共振峰输出的完整实现思路理解杜宾与牛顿算法在LPC系数求解中的差异并借助工程日志定位编译与调试问题。1. 语音共振峰提取为什么我盯着 LPC 这条链路不放LPC 法语音共振峰提取的 C 语言实现最早是我在一个 16 kHz 单声道录音的小工具里被逼着啃下来的。当时频谱图上 F1、F2 本该平稳的元音段逐帧看却像心电图换了几版 FFT 峰值检测都没压住。后来把链路拆成预加重、汉明窗、自相关、杜宾递推、牛顿迭代求根五段才发现问题不在“找峰”而在全极点模型和根筛选没有对齐。这份资源就是一套纯 C 的落地实现不依赖大型信号库按帧读 PCM算 LPC 系数解多项式根再换算成共振峰频率和带宽。适合做语音前端、教学实验、嵌入式低频分析的人也适合想弄明白“为什么阶数和窗长一换结果就翻车”的熟手。2. LPC 共振峰链路拆解从汉明窗到杜宾递推的五个环节先别急着把代码整段抄进工程。共振峰提取不是对一帧信号做一次 FFT 再挑峰值而是用全极点模型去逼近声道响应再从多项式根的幅角和模长里读频率与带宽。链路顺序错一步后面调参就会变成玄学窗加错了自相关偏自相关偏了杜宾递推出来的预测系数符号反符号反了牛顿法求根会收敛到共轭镜像上最后 F1/F2 全部对不上。2.1 预加重、分帧与汉明窗先把短时平稳这件事做扎实语音在 2030 ms 内可以近似短时平稳所以分帧长度通常取 2030 ms帧移取 10 ms 左右。16 kHz 采样率下25 ms 对应 400 点10 ms 帧移对应 160 点。预加重放在分帧之前或之后都可以常见做法是先对整段信号做一阶高通系数取 0.950.9716 kHz 录音我一般用 0.97目的是把高频端抬起来抵消声门激励带来的低频倾斜。#define PI 3.14159265358979323846f #define FRAME_LEN 400 #define FRAME_SHIFT 160 static void pre_emphasis(const float *in, float *out, int n, float alpha) { out[0] in[0]; for (int i 1; i n; i) { out[i] in[i] - alpha * in[i - 1]; } } static void hamming_window(float *x, int n) { for (int i 0; i n; i) { x[i] * 0.54f - 0.46f * cosf(2.0f * PI * i / (n - 1)); } }alpha控制预加重强度太大会把高频噪声一起抬起来太小则低频压不住n是帧长汉明窗分母用n - 1是为了让窗两端不完全落到零减少频谱泄漏。加窗后帧能量会下降但自相关和杜宾递推看的是相对关系不需要把窗后能量归一化回原幅值。2.2 自相关序列与杜宾递推把预测系数解出来自相关序列r[k]用加窗后的帧算累加时用双精度避免 16 位 PCM 转 float 后长帧累加出现精度抖动。杜宾递推Levinson-Durbin负责从r[0..p]解出预测系数a[1..p]和预测误差err。这里最容易翻车的是符号常见定义是x[n] ≈ sum(a[i] * x[n-i])多项式写成A(z) 1 - sum(a[i] * z^{-i})后面求根必须用这个多项式不能把减号写成加号。#define MAX_ORDER 32 int levinson_durbin(const float *r, int p, float *a, float *err) { float tmp[MAX_ORDER 1]; a[0] 1.0f; *err r[0]; if (*err 1e-9f) return -1; for (int i 1; i p; i) { float acc r[i]; for (int j 1; j i; j) { acc - a[j] * r[i - j]; } float k acc / *err; // 反射系数 if (k 1.0f || k -1.0f) return -2; 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]; } *err * (1.0f - k * k); if (*err 1e-12f) return -3; } return 0; }p是 LPC 阶数16 kHz 下常用 1012r[0]太小说明这一帧接近静音直接跳过比硬算更稳反射系数k绝对值接近 1 时递推会数值发散返回失败让上层丢帧或降阶。err是预测误差递推过程中如果变成负数基本可以判定自相关或符号出了问题。2.3 带 Newton 迭代的根求解从多项式根到共振峰频率把A(z)乘以z^p得到P(z) z^p - a[1] z^{p-1} - ... - a[p]求它的复根。牛顿迭代的初值不能全放在实轴上否则容易收敛到同一个实根或共轭镜像。我一般把单位圆内按角度均匀撒初值半径取 0.9角度取2π(k0.5)/pk 0..p-1再逐个迭代。求值和求导都用 Horner导数不是数值差分避免复数步长选不好导致不收敛。typedef struct { float re, im; } complex_t; static complex_t poly_eval_deriv(const float *a, int p, complex_t z, complex_t *deriv) { complex_t acc {1.0f, 0.0f}; complex_t d {0.0f, 0.0f}; for (int i 1; i p; i) { complex_t nd; nd.re d.re * z.re - d.im * z.im acc.re; nd.im d.re * z.im d.im * z.re acc.im; d nd; complex_t na; na.re acc.re * z.re - acc.im * z.im - a[i]; na.im acc.re * z.im acc.im * z.re; acc na; } *deriv d; return acc; } static int newton_root(const float *a, int p, complex_t *z, int max_iter, float eps) { for (int it 0; it max_iter; it) { complex_t d, f poly_eval_deriv(a, p, *z, d); float den d.re * d.re d.im * d.im; if (den 1e-20f) return -1; complex_t step; step.re (f.re * d.re f.im * d.im) / den; step.im (f.im * d.re - f.re * d.im) / den; z-re - step.re; z-im - step.im; if (step.re * step.re step.im * step.im eps * eps) return 0; } return -2; }求出根后频率f atan2(im, re) / (2π) * fs带宽BW -ln(sqrt(re² im²)) / π * fs。只保留虚部大于 0 的根按频率排序再按带宽阈值过滤。eps取1e-6左右max_iter取 30 足够如果某个初值 30 次不收敛直接丢弃不要让失败根污染后面的筛选。2.4 参数表与 C 语言数据结构一份能直接落地的接口约定把配置和输出拆成两个结构体后面换采样率、换窗长只改配置不动算法主体。我习惯把lpc_order、frame_len、frame_shift、pre_emphasis、bandwidth_max都放进配置调试时打印整份配置比只打印一个频率值更容易回溯。typedef struct { int sample_rate; int frame_len; int frame_shift; int lpc_order; float pre_emphasis; float bandwidth_max; int max_formants; } lpc_config_t; typedef struct { float freq_hz; float bandwidth_hz; } formant_t;参数16 kHz 常用值作用改大后的典型影响frame_len400 点25 ms短时分析窗长频率分辨率变细但时间定位变糊frame_shift160 点10 ms帧间步进帧数变少过渡段容易漏lpc_order1012全极点模型阶数根变多虚根和带宽异常根增加pre_emphasis0.97高频抬升高频噪声被放大静音帧更敏感bandwidth_max400500 Hz共振峰带宽上限放宽后噪声根混入收紧后真实 F3 被丢3. 在 C 里搭一条可复现的提取流水线从 PCM 到前三个共振峰这一章把上一章的零件串成能编译、能跑出数字的流水线。输入建议先用 16 kHz、16 bit、单声道 WAV别一上来就啃 48 kHz 双声道。双声道要下混48 kHz 要重采样任何一个环节没对齐共振峰频率都会整体偏移。流水线顺序固定为读 WAV → 转 float 并归一化 → 逐帧预加重 → 加汉明窗 → 自相关 → 杜宾递推 → 牛顿求根 → 筛选排序 → 落盘。3.1 读 WAV 与预处理的 C 实现标准 44 字节头的 PCM WAV 可以先按固定偏移读但工程里更稳的做法是扫描fmt和data块。下面这段只处理最常见的 44 字节头目的是让示例能跑遇到带LIST块的文件固定偏移会读错数据长度。读进来后立刻除以 32768.0f统一到[-1, 1]后面调参才不会因为文件不同而出现数量级差异。#include stdio.h #include stdlib.h #include string.h static int read_wav_pcm16(const char *path, int *fs, float **samples, int *n) { FILE *fp fopen(path, rb); if (!fp) return -1; unsigned char hdr[44]; if (fread(hdr, 1, 44, fp) ! 44) { fclose(fp); return -2; } int channels hdr[22] | (hdr[23] 8); int sample_rate hdr[24] | (hdr[25] 8) | (hdr[26] 16) | (hdr[27] 24); int bits hdr[34] | (hdr[35] 8); int data_bytes hdr[40] | (hdr[41] 8) | (hdr[42] 16) | (hdr[43] 24); if (channels ! 1 || bits ! 16 || data_bytes 0) { fclose(fp); return -3; } int count data_bytes / 2; short *buf (short *)malloc((size_t)count * sizeof(short)); float *out (float *)malloc((size_t)count * sizeof(float)); if (!buf || !out) { free(buf); free(out); fclose(fp); return -4; } if (fread(buf, sizeof(short), (size_t)count, fp) ! (size_t)count) { free(buf); free(out); fclose(fp); return -5; } for (int i 0; i count; i) { out[i] buf[i] / 32768.0f; } free(buf); fclose(fp); *fs sample_rate; *samples out; *n count; return 0; }channels、bits、data_bytes从头部固定偏移取小端拼接data_bytes / 2是 16 bit 采样点数。count为 0 或fread数量不对时返回错误不要继续往下算。这个函数只适合单声道 16 bit PCM遇到 IEEE float WAV 或 24 bit 要先转换否则归一化系数全错。3.2 自相关、杜宾与 LPC 系数自相关用双精度累加k从 0 到p每个k只累加i k..n-1别把边界写错。帧内如果r[0]小于1e-6直接判为静音帧跳过不进入杜宾递推。加窗之后不要再做一次预加重预加重在分帧前做一次就够重复做会把高频抬得过分F2 容易虚高。static void autocorr(const float *x, int n, float *r, int p) { for (int k 0; k p; k) { double sum 0.0; for (int i k; i n; i) { sum (double)x[i] * (double)x[i - k]; } r[k] (float)sum; } } static int frame_lpc(const float *frame, int n, int p, float *a) { float r[MAX_ORDER 1]; float err 0.0f; autocorr(frame, n, r, p); if (r[0] 1e-6f) return 1; // 静音帧 int ret levinson_durbin(r, p, a, err); if (ret ! 0) return 2; // 递推失败 return 0; }autocorr的n是帧长p是阶数r长度至少p 1。frame_lpc返回 1 表示静音返回 2 表示递推失败调用方按帧记录失败原因后面排查时能看出是能量问题还是数值问题。不要把所有失败都吞成 0否则结果表里会出现成片空缺却不知道原因。3.3 多项式求根与共振峰筛选求根之后不要直接取前三个根。先过滤虚部小于等于 0 的共轭根再过滤模长不在[0.5, 1.2]的根然后算频率和带宽按90 Hz f 0.45 * fs、20 Hz BW bandwidth_max筛。同一个根可能被两个初值收敛到频率差小于 20 Hz 时只保留一个。最后按频率升序排取前max_formants个。static int cmp_formant(const void *a, const void *b) { const formant_t *fa (const formant_t *)a; const formant_t *fb (const formant_t *)b; return (fa-freq_hz fb-freq_hz) - (fa-freq_hz fb-freq_hz); } static int roots_to_formants(const float *a, int p, int fs, float bw_max, formant_t *out, int max_out) { int cnt 0; for (int k 0; k p; k) { float theta 2.0f * PI * (k 0.5f) / p; complex_t z {0.9f * cosf(theta), 0.9f * sinf(theta)}; if (newton_root(a, p, z, 30, 1e-6f) ! 0) continue; float mag sqrtf(z.re * z.re z.im * z.im); if (mag 0.5f || mag 1.2f) continue; if (z.im 0.0f) continue; float freq atan2f(z.im, z.re) / (2.0f * PI) * (float)fs; float bw -logf(mag) / PI * (float)fs; if (freq 90.0f || freq 0.45f * (float)fs) continue; if (bw 20.0f || bw bw_max) continue; int dup 0; for (int i 0; i cnt; i) { if (fabsf(out[i].freq_hz - freq) 20.0f) { dup 1; break; } } if (dup) continue; out[cnt].freq_hz freq; out[cnt].bandwidth_hz bw; cnt; if (cnt max_out) break; } qsort(out, (size_t)cnt, sizeof(formant_t), cmp_formant); return cnt; }bw_max是带宽上限16 kHz 下我一般先给 500 Hz再根据结果收紧到 400 Hz。max_out不要超过 5超过之后 F4 以上很容易被噪声根占掉。qsort之后 F1、F2、F3 的顺序才稳定否则同一帧两次运行可能给出不同顺序批处理对比时会误判成频率跳动。3.4 编译、运行与结果对照把上面的函数放进一个main.c加一个逐帧循环for (start 0; start frame_len n; start frame_shift)拷贝帧、预加重、加窗、算 LPC、求根、打印。编译时开-O2 -Wall先保证没有隐式声明和越界警告再谈调参。gcc -O2 -Wall -Wextra -stdc11 main.c -lm -o lpc_formant ./lpc_formant input_16k_mono.wav formants.csv输出 CSV 至少包含frame_index,time_ms,f1,f2,f3,bw1,bw2,bw3。拿一段稳定元音对照16 kHz 下常见范围是/a/ 的 F1 大约 700900 Hz、F2 大约 11001400 Hz/i/ 的 F1 大约 250350 Hz、F2 大约 22002800 Hz/u/ 的 F1 大约 300400 Hz、F2 大约 700900 Hz。不同发音人、不同录音设备会有偏移但同一段稳定段内不应该逐帧大幅跳动。如果 F1 整体偏低 100 Hz 以上先查预加重和窗长如果 F2 整体偏高且带宽很大先查阶数是不是给高了。检查项正常表现异常时先看静音帧返回 1不输出共振峰是否把噪声帧硬算成 F1杜宾递推返回 0err为正反射系数是否越界根的数量每帧 35 个有效根是否全是虚部小于 0 的镜像根F1 稳定性相邻帧差通常小于 100 Hz帧长是否短于 20 ms4. 参数怎么调阶数、窗长、预加重和带宽阈值的边界同样一份 C 代码参数换一组结果可以从“能用”变成“完全不能看”。这一章把四个最常动的参数拆开LPC 阶数、窗长与帧移、预加重系数、带宽阈值。调参顺序建议是先固定帧长和帧移再定阶数最后收紧带宽筛选。反过来先调带宽会把阶数不够导致的合并峰误当成噪声根丢掉。4.1 阶数 p 与采样率、声道长度的关系阶数p决定全极点模型能表示多少个共振峰。经验上每 1 kHz 采样率大约需要 1 个极点16 kHz 常用 10128 kHz 常用 81044.1 kHz 常用 1620。但这不是硬公式声道长度、发音人性别、录音带宽都会影响。阶数太低F1 和 F2 容易合并成一个宽峰阶数太高多项式根变多噪声根和虚根混进来带宽筛选稍微一松就冒出假的 F3。我一般先取p 12跑一遍稳定元音看 F1/F2 是否落在合理范围再看带宽。如果 F1 和 F2 分不开先把p加到 14如果 F3 之后出现大量带宽超过 600 Hz 的根把p降回 10 或 12。不要一上来就上 20 阶C 里求根的计算量翻倍不说牛顿迭代初值也要跟着加密否则漏根更严重。4.2 窗长与帧移时间分辨率与频率分辨率的拉锯窗长决定频率分辨率帧移决定时间分辨率。16 kHz 下 25 ms 窗是 400 点频率分辨率大约 40 Hz30 ms 窗是 480 点频率分辨率大约 33 Hz。分析稳定元音可以用 30 ms分析辅音到元音的过渡段用 2025 ms 更稳。窗太长辅音段和元音段被混在一帧里F1 会先被拉低再弹回窗太短低频分辨率不够F1 容易在相邻帧之间来回跳。帧移一般取窗长的 1/3 到 1/2。25 ms 窗配 10 ms 帧移是常见组合既不会漏掉过渡段也不会让输出帧数爆炸。做批处理时我习惯把帧移固定成 10 ms换窗长只改frame_len这样时间轴不变前后两版结果可以直接按time_ms对齐。如果帧移大于窗长的一半相邻帧重叠太少共振峰轨迹会出现台阶。4.3 共振峰筛选与带宽阈值为什么多出来的根要丢掉LPC 多项式解出来的根不一定都是共振峰。接近单位圆的根对应窄带共振模长小、带宽大的根往往是噪声或模型拟合出来的虚根。带宽公式BW -ln(mag) / π * fs里模长越接近 1带宽越小。16 kHz 下bandwidth_max先给 500 Hz再看结果收紧到 400 Hz通常能把大部分噪声根挡掉。筛选时不要只按“前三个根”取。正确顺序是先过滤负频率镜像根再过滤模长异常根再算频率和带宽再按频率范围过滤最后按频率排序取前几个。F1 一般落在 2001000 HzF2 落在 8003000 HzF3 落在 15004000 Hz但这不是绝对边界只是排查时的参考。如果 F1 跑到 1500 Hz 以上先看是不是把 F2 当成了 F1再看带宽是不是大得离谱。调参动作适用场景风险p从 10 加到 12F1/F2 分不开虚根增多F3 不稳窗长从 20 ms 加到 30 ms稳定元音分析过渡段被抹平pre_emphasis从 0.95 加到 0.97高频弱、F2 不明显静音帧噪声被抬高bandwidth_max从 500 收到 400噪声根多真实高带宽共振峰被丢5. 避坑与排查共振峰乱跳、虚根、数值溢出这三类翻车这一章记录的是我在 C 实现里真实踩过的四类问题。共同点是代码能跑结果也能出但数字不可信。排查时不要只看最终 F1/F2要把中间量打出来r[0]、反射系数、err、根的模长、带宽。很多时候问题不在求根而在进求根之前就已经歪了。5.1 现象F1 在相邻帧跳 300 Hz 以上现象很直观稳定元音段导出 CSVF1 一帧 600 Hz、下一帧 950 Hz再下一帧又回到 700 HzF2 也跟着抖。原因通常有三个窗长太短、阶数太高、根匹配没做。20 ms 窗在 16 kHz 下只有 320 点低频分辨率不够F1 会在两个根之间来回选阶数给到 16 以上多出来的根让筛选边界变得敏感相邻帧各自排序取前三个没有做频率连续性约束F1 和 F2 可能互换。解决分两步。先把窗长固定到 2530 ms阶数降到 12看跳动是否收窄。如果还跳再加一层简单的帧间平滑当前帧的 F1 只接受与上一帧 F1 相差 150 Hz 以内的候选根超出就标记为可疑帧不要直接覆盖。下面这段不是万能平滑只是把候选根按上一帧频率做一次最近邻选择。static int pick_nearest(const formant_t *cand, int cnt, float prev_freq) { int best -1; float best_dist 1e9f; for (int i 0; i cnt; i) { float d fabsf(cand[i].freq_hz - prev_freq); if (d best_dist) { best_dist d; best i; } } if (best 0 best_dist 150.0f) return best; return -1; }prev_freq是上一帧同一阶共振峰的频率150.0f是容忍窗口。超过窗口不要硬拉宁可输出空值也不要把 F2 强行贴到 F1 上。平滑只用于显示或批处理对齐不要拿平滑后的值去反推 LPC 系数。5.2 现象根求解不收敛或返回复数对异常现象是某些帧的newton_root返回 -2或者收敛出来的根实部很大、虚部接近 0换算出的频率要么是 0 Hz要么是奈奎斯特附近。原因常见于三种初值全放在实轴上、多项式符号写反、导数 Horner 写错。初值全实轴时牛顿法容易沿着实轴收敛到实根而共振峰对应的是共轭复根符号写反时P(z)的根会关于单位圆镜像频率看起来对但带宽完全不对。排查时先把a[1..p]打印出来用一个小脚本或手算验证A(z)在某个已知根附近是否接近 0。再把初值改成0.9 * exp(j * 2π(k0.5)/p)确保每个初值都有非零虚部。导数 Horner 的更新顺序是“先d d*z acc再acc acc*z - a[i]”这两步顺序反了导数会差一阶迭代会震荡。若某帧 30 次不收敛直接丢弃该帧不要用上一次的根顶替。5.3 现象杜宾递推出现 NaN 或反射系数越界现象是err变成NaN或者k的绝对值大于等于 1函数返回 -2。原因通常是r[0]太小、自相关累加溢出、或者前置预加重把静音帧抬成了噪声。16 位 PCM 转 float 后如果没归一化r[0]可能到1e9量级err累乘后精度损失如果归一化了但静音帧r[0]只有1e-8除法又会放大噪声。解决是两道门进门先判r[0] 1e-6f直接跳过出门判k 1.0f || k -1.0f返回失败。自相关累加用doubler存float可以但累加过程不要用float。另外预加重系数不要大于 0.98接近 1 时高频提升过猛静音帧的自相关条件数会变差。批量跑的时候把失败帧计数打出来失败率超过 10% 就说明前置处理有问题不是算法本身的问题。5.4 现象不同 WAV 读入幅值差 32768 倍现象是同一套参数A 文件算出来的共振峰正常B 文件 F1 直接跑到 0 Hz 或奈奎斯特CSV 里数值量级也不对。原因多半是 PCM 没有归一化short直接赋给float幅值在[-32768, 32767]自相关和预测误差全部放大另一种情况是 WAV 头不是 44 字节固定偏移读到了错误data_bytes把头部当采样点算进去。解决是读入后立刻除以32768.0f并在读头时校验channels 1 bits 16。如果工程里必须支持带LIST块的 WAV不要用固定偏移改成逐块扫描fmt和data。调试时打印前 10 个采样点的最大值和r[0]归一化后r[0]应该在1e-3到1e2量级如果看到1e8以上基本就是没归一化。提示排查顺序固定为“读入幅值 → 帧能量 → 自相关 → 杜宾返回码 → 根模长 → 带宽”不要一上来改求根迭代次数很多不收敛是前面传进来的系数已经歪了。6. 进阶验证与工程化技巧用合成元音和 LPC 谱包络交叉核对调参调到后面最容易出现“看着像对但不知道对不对”。我后来养成一个习惯每换一次采样率或窗长先跑一段合成元音把已知共振峰频率喂进去看提取结果能不能回到设定值附近。合成不用太复杂用二阶谐振器叠三个共振峰给不同带宽生成 1 秒 16 kHz 单声道 float 序列再走同一条 C 流水线。设定 F1700、F21200、F32500带宽分别给 80、100、150 Hz提取结果如果偏差超过 10%先别怀疑算法查窗长和阶数是不是不匹配。第二个交叉核对是 LPC 谱包络。求出a[1..p]后计算H(e^{jω}) 1 / A(e^{jω})的幅度在 0fs/2 上扫一遍找峰值位置和根换算出来的频率对比。两者应该大致重合如果根给出的 F2 在 1800 Hz而谱包络峰值在 1200 Hz说明求根筛选把某个根丢了或者排序后取错了阶。这个检查不需要额外库几十行 C 就能写调试时比盯着 CSV 猜有用。第三个习惯是批处理落盘时保留中间列frame_index,time_ms,r0,err,k_max,f1,f2,f3,bw1,bw2,bw3,valid。valid标记静音、递推失败、根不足三种状态。后面做统计时先按valid过滤再看频率分布不会把失败帧的 0 值混进均值。从那以后我每次换采样率或换录音设备都先跑合成元音和谱包络这两道检查再去看真实录音。希望帮到你。本文还有配套的精品资源点击获取