★ FEATURED ARTICLE 简介面向语音信号处理初学者与嵌入式开发者的LPC共振峰提取C语言实现工程围绕线性预测编码这一经典方法解决从语音采样到共振峰频率估计的完整链路问题。工程将语音分帧、汉明窗加窗、端点检测、杜宾递推求解预测系数以及牛顿算法迭代优化等环节分别落地为独立模块并附readme说明与实测语音样本便于对照调试与算法验证。压缩包共72个文件、约274KB以C源码、工程配置文件、编译中间产物与日志为主另有少量文本数据和一段wav语音可直观看到各算法模块的构建过程与运行结果。目前已有385人学习下载适合课程设计、语音识别入门实验或算法移植参考读者可从中获得可直接编译运行的LPC分析框架、分帧加窗与预测系数求解的具体写法以及将共振峰提取嵌入实际工程的排错思路为后续做语音特征分析打下基础。1. 把一段 16 kHz 录音逐帧变成 F1、F2、F3要给一段 16 kHz 录音逐帧报出 F1、F2、F3用脚本调库两行就完了可一旦要嵌进采集卡上位机或者塞进一颗跑不了操作系统的芯片就只剩 C 语言这一条路。这个项目拆的就是这条路预加重、分帧、汉明窗、自相关杜宾递推拿到 LPC 系数再用牛顿迭代在复平面上求极点映射成共振峰频率和带宽。整套代码只依赖 libm不拉 FFT 库也不绑任何音频框架。适合两类人准备把算法搬进 C 工程的和想弄明白 LPC 到底在算什么的人。如果只是离线跑批出报告现成工具够用不必自己写。2. 工程骨架与处理链路从裸 PCM 到一张 formant.csv拿到这类项目我一般先不看算法实现先看数据怎么流动。因为共振峰提取这个活儿算法写错了顶多是数值难看数据格式和采样率对不上出来的曲线能离谱到让你怀疑人生——所有频率按比例缩放F1 跑到 F2 的位置还特别稳定看不出是错的。2.1 六个阶段与它们之间的数据形状整条链路是严格串行的前一阶段的输出就是后一阶段的输入中间没有分支。把这个表记住后面看代码就不会迷路。阶段输入输出关键参数预加重整段 float 采样点高频抬升后的采样点系数 a0.97分帧整段信号 帧号一帧 512 点帧长 512 / 帧移 256加窗512 点原始帧512 点加窗帧汉明窗自相关512 点r[0..p] 共 13 个值阶数 p12杜宾递推r[0..p]a[1..p] 预测误差 Ep求根与判决a[]F1/F2/F3 各自带宽收敛阈值、带宽上限预加重和分帧的先后顺序容易搞反。常见做法是先对整段信号做一遍预加重再切帧这样帧与帧的边界上是连续的如果先切帧再逐帧预加重每帧第一个点拿不到前一个采样边界处会引入一个极小但周期性的跳变做频谱分析时表现为帧率的边带。我在早期的实现里就吃过这个亏F2 的曲线上一直叠着一条 62.5 Hz 的纹波16000/256排查了半天才想起来是分帧顺序的问题。杜宾递推输出的是预测系数 a不是共振峰。这一点在新手代码里非常常见拿到 a[1..p] 直接当特征送进分类器或者干脆对 a 做谱峰搜索。a 是极点多项式的系数物理意义是用前 p 个采样点线性预测当前点的最优权重它和共振峰之间隔着一个求根或者一次频谱求值的距离。2.2 目录结构与编译一个 Makefile 搞定工程按职责拆文件每个 .c 只做一件事方便你后面把某一块单独抠出来复用。lpc-formant/ ├── Makefile ├── include/ │ ├── preprocess.h # 预加重、分帧、加窗 │ ├── lpc.h # 自相关、杜宾递推 │ └── formant.h # 求根、极点筛选、输出结构体 ├── src/ │ ├── main.c # 命令行解析 帧循环 │ ├── preprocess.c │ ├── autocorr.c │ ├── durbin.c │ ├── polyroot.c # 牛顿迭代求复根 │ ├── formant.c # 频率/带宽映射与平滑 │ └── io_pcm.c # 裸 PCM 读取、CSV 写出 ├── test/ │ └── gen_synth.c # 合成已知共振峰的测试信号 └── data/ └── sample16k.pcm构建用最朴素的 Makefile不引入 CMake 那一套。理由很简单这东西最终大概率要交叉编译到 ARM 上Makefile 改两行就能换工具链CMake 得重写 toolchain 文件。CC gcc CFLAGS -O2 -stdc11 -Wall -Wextra -Iinclude LDLIBS -lm SRC $(wildcard src/*.c) bin/formant: $(SRC) mkdir -p bin $(CC) $(CFLAGS) -o $ $(SRC) $(LDLIBS)编译和运行make ./bin/formant -i data/sample16k.pcm -r 16000 -p 12 -o out/formant.csv命令行参数不多但每一个都能让结果翻车我按踩坑频率排一下参数含义默认常见误用-i输入裸 PCM 路径无直接喂 .wav44 字节文件头被当成 22 个采样点-r采样率 Hz16000填 8000 却给 16k 数据所有频率减半-pLPC 阶数12女声/儿童语音还用 8F3 直接被抹平-o输出 CSVstdout重定向到不存在的目录静默失败关于-stdc11和-lm牛顿迭代那段用了complex.h里的double complex这是 C99 起就有的标准特性但M_PI不是标准 C 的一部分得在头文件里自己定义或者编译时加-D_GNU_SOURCE。我习惯在preprocess.h里写死一行#define PI 3.14159265358979323846交叉编译到某些精简 libc 上时这个宏反而更省心。2.3 输入数据的三个硬约定采样率、位深、字节序工程只吃裸 PCM不解析 wav 头。这不是偷懒是故意的现场采集设备导出的经常就是裸流而 wav 头里的字段各家写法不一有的factchunk 长度都是错的。约定三条单声道、16 bit 有符号、小端。三条里任何一条不满足代码不会报错只会输出垃圾。/* io_pcm.c裸 PCM - float顺手做去直流 */ int pcm_load_f32(const char *path, float **out, size_t *n) { FILE *fp fopen(path, rb); if (!fp) return -1; fseek(fp, 0, SEEK_END); long bytes ftell(fp); fseek(fp, 0, SEEK_SET); if (bytes 0 || (bytes 1)) { fclose(fp); return -2; } /* 16 bit 必须偶数长度 */ size_t cnt (size_t)bytes / 2; int16_t *raw (int16_t *)malloc(cnt * sizeof(int16_t)); float *buf (float *)malloc(cnt * sizeof(float)); if (!raw || !buf) { free(raw); free(buf); fclose(fp); return -3; } if (fread(raw, sizeof(int16_t), cnt, fp) ! cnt) { free(raw); free(buf); fclose(fp); return -4; } fclose(fp); /* 去直流录音设备偏置会让 r[0] 虚高杜宾递推第一级就跑偏 */ double mean 0.0; for (size_t i 0; i cnt; i) mean raw[i]; mean / (double)cnt; for (size_t i 0; i cnt; i) buf[i] (float)((double)raw[i] - mean) / 32768.0f; free(raw); *out buf; *n cnt; return 0; }去直流这一步看起来可有可无实测影响不小。直流分量全部落在 r[0] 上会让杜宾递推的第一级反射系数 k1 接近 1后面的系数被压得极小极点集体往 z1 附近挤最后求出来的 F1 永远卡在几十赫兹。如果设备本身有几十个 LSB 的偏置不去直流基本等于没法用。字节序的问题在 x86 上永远不会暴露一上 ARM 大端或者某些 DSP 就全乱。稳妥写法是不用fread直接读int16_t而是读成uint8_t数组自己拼(int16_t)((uint16_t)lo | ((uint16_t)hi 8))。这段代码我留在io_pcm.c的注释里需要的时候打开就行。3. 前端预处理预加重、分帧与汉明窗预处理这三个步骤写得对不对直接决定后面杜宾递推拿到的是不是人声。很多跑通了但结果不对的案例问题都出在这一层而不是求根算法。3.1 预加重0.97 这个系数到底在补什么语音的频谱天然是往下掉的高频能量比低频低一二十个 dB这是声门脉冲激励的包络决定的。LPC 是个全极点模型它倾向于把能量大的低频拟得很准高频段几乎不管。预加重就是用一个一阶高通把频谱抬平让各频段在最小二乘意义下有差不多的权重。$$y[n] x[n] - a \cdot x[n-1], \quad a \approx 0.9 \sim 0.98$$/* preprocess.c原地预加重从后往前避免覆盖未使用的 x[i-1] */ void pre_emphasis(float *x, size_t n, float a) { if (n 0) return; for (size_t i n; i-- 1; ) /* i 从 n-1 递减到 1 */ x[i] x[i] - a * x[i - 1]; x[0] x[0] * (1.0f - a); /* 首点没有前文按同增益缩放避免爆点 */ }参数说明a取 0.97 是最常见的经验值对应 16 kHz 采样时大约 200 Hz 以下的轻度衰减采样率降到 8 kHz 时习惯取 0.95因为高频段本来就窄压太狠会把 F3 削掉。a0等价于不处理调试时可以拿它做对照。首点x[0] * (1-a)这个处理不是必须的但不做的话整段信号的第一个采样点会保持原幅值和后面被衰减过的点不在一个量级上做自相关时相当于在起点插了一个冲激会给 r[k] 引入一点高频纹波。单帧内影响很小可如果你要拿这段代码做整段信号的能量统计这个尖峰会被当成一次瞬态。还有一个容易忽略的点如果输入音频已经做过高频提升有些采集前端自带再预加重一次等于压了两次F3 的位置会明显偏高。遇到结果系统性偏高的案例先确认一下录音链路里有没有做过类似处理。3.2 分帧参数帧长 512、帧移 256 是怎么定下来的帧长得覆盖至少两到三个基音周期否则这一帧里连一个完整的周期都凑不齐自相关算出来的 r[k] 没有物理意义。男声基频最低到 70 Hz 左右周期 14 ms三个周期 42 ms女声和童声高得多不需要那么长。反过来帧太长声道形状在帧内已经变了LPC 假设的短时平稳不成立共振峰会被平均掉。折中下来32 ms 是个大家都接受的数。16 kHz 采样就是 512 点8 kHz 是 256 点。采样率帧长(点)时长帧移推荐阶数 p800025632 ms1288 ~ 101600051232 ms25612 ~ 142205070432 ms35214 ~ 1644100102423 ms51216 ~ 20帧移取帧长的一半50% 重叠是为了让帧与帧之间的共振峰轨迹平滑一些。想再密一点可以取 1/4代价是计算量翻倍、相邻帧高度相关做平滑时反而容易把真实的过渡抹掉。/* preprocess.c取第 idx 帧并加汉明窗 */ void frame_window(const float *x, size_t n, int flen, int shift, int idx, float *buf) { size_t start (size_t)idx * (size_t)shift; for (int i 0; i flen; i) { size_t k start (size_t)i; buf[i] (k n) ? x[k] : 0.0f; /* 尾帧补零绝不越界读 */ } /* 汉明窗0.54 - 0.46cos端点不为零主瓣略宽于汉宁窗 */ for (int i 0; i flen; i) buf[i] * 0.54f - 0.46f * cosf(2.0f * PI * (float)i / (float)(flen - 1)); }窗函数的选择上汉明窗和汉宁窗在这里差别不大都能把帧两端的突变压下去。汉明窗第一旁瓣 -43 dB比汉宁窗低代价是主瓣稍宽频率分辨率略差。共振峰本身就有一两百赫兹的宽度这点分辨率损失无所谓所以工程里基本都用汉明窗。矩形窗千万别用帧边界的不连续会在频谱上拖出一片旁瓣杜宾递推会把这些旁瓣当成真实谱峰去拟合最后多出一堆带宽极窄的假极点。尾帧补零之后这一帧的有效样本数不足能量统计会偏低。稳妥做法是记录有效长度能量低于门限就直接跳过不要把它当成静音帧参与平滑否则末尾几帧的共振峰会莫名其妙地往平均值上掉。3.3 自相关与加窗带来的偏置自相关只需要算到 p 阶不用算全长。r[k] Σ x[n]x[n-k]k 从 0 到 p13 个值就够杜宾递推用了。/* autocorr.c返回 0 表示成功-1 表示静音帧 */ int autocorr(const float *x, int flen, int p, double *r) { double e 0.0; for (int i 0; i flen; i) e (double)x[i] * (double)x[i]; if (e 1e-6) return -1; /* 能量门限静音帧交给上层丢帧 */ for (int k 0; k p; k) { double s 0.0; for (int i k; i flen; i) s (double)x[i] * (double)x[i - k]; r[k] s; } return 0; }参数说明p就是 LPC 阶数也是自相关要算的阶数两者必须一致e的门限 1e-6 是按 16 bit 归一化到 ±1 之后的经验值如果你的信号幅度不在这个量级门限要跟着改否则要么全都判成静音要么静音帧混进来做除法。加窗会带来一个理论上的偏置加窗后的自相关等于原始自相关乘上窗函数自相关的三角包络。严格做法是除以这个包络来补偿但那会放大高阶项的噪声。实际工程里几乎都不补因为杜宾递推里每一步的反射系数是个比值k acc / E整体缩放会同时作用在分子分母上被约掉了。真正受影响的是预测误差 E 的绝对值如果你打算拿 E 去做浊音/清音判决就得注意它比真实值偏小一点阈值要相应下调。用double累加而不是float在 512 点这个规模下差异不明显可一旦采样率上到 44.1 kHz、帧长过千float 的 24 位有效位就不够用了log 出来的共振峰轨迹会出现细碎的锯齿。4. 杜宾递推接牛顿迭代LPC 系数与极点这一章是全篇最硬的部分。前面准备的数据到这里变成一个 12 阶多项式我们要做的就是把这个多项式的根找出来再把根翻译成频率。4.1 杜宾递推逐行拆一次循环拿到 p 个系数杜宾递推Levinson-Durbin解决的是这个正规方程给定自相关 r[0..p]求使预测误差最小的系数 a[1..p]。直接解 Toeplitz 矩阵是 O(p³)杜宾利用矩阵的对称结构把它降到 O(p²)而且中间产物——反射系数 k——本身就是很好的稳定性指标。递推公式i 从 1 到 p$$k_i \frac{r[i] - \sum_{j1}^{i-1} a_j^{(i-1)} r[i-j]}{E_{i-1}}$$$$a_j^{(i)} a_j^{(i-1)} - k_i \cdot a_{i-j}^{(i-1)}, \quad a_i^{(i)} k_i, \quad E_i (1 - k_i^2) E_{i-1}$$/* durbin.ca[0..p-1] 对应 a_1..a_p误差写回 *err */ int durbin(const double *r, int p, double *a, double *err) { double tmp[MAX_ORDER]; double e r[0]; if (e 1e-12) return -1; /* 能量为零无法递推 */ for (int i 0; i p; i) a[i] 0.0; for (int i 0; i p; i) { /* i 对应公式里的 i1 */ double acc r[i 1]; for (int j 0; j i; j) acc - a[j] * r[i - j]; double k acc / e; if (fabs(k) 1.0) return -2; /* 反射系数越界该帧不稳定 */ tmp[i] k; /* a_{i1}^{(i1)} k */ for (int j 0; j i; j) tmp[j] a[j] - k * a[i - 1 - j]; for (int j 0; j i; j) a[j] tmp[j]; e * (1.0 - k * k); } *err e; return 0; }逐段说明acc算的是公式里的分子注意r[i-j]的下标j0 时取 r[i]ji-1 时取 r[1]都不要碰到 r[0]因为 r[0] 已经在e里了。tmp[j] a[j] - k*a[i-1-j]这一行是递推的核心反向索引i-1-j是杜宾递推区别于普通最小二乘的地方它保证了每一步算出的多项式都是最小相位的。最后e * (1-k*k)更新预测误差同时天然给出了稳定性判据只要所有 |k| 1极点必定在单位圆内。符号约定必须和求根那一步对齐。这里得到的多项式是 $A(z) 1 a_1 z^{-1} \cdots a_p z^{-p}$极点满足 $z^p a_1 z^{p-1} \cdots a_p 0$。如果你的代码里递推公式写成了acc a[j]*r[i-j]那系数符号是反的求根时要用 $z^p - a_1 z^{p-1} - \cdots$。这个坑非常隐蔽症状是所有极点跑到单位圆外、带宽算出负值。我第一次移植的时候就栽在这里对着公式核了两小时。4.2 阶数、带宽扩展与稳定性判断阶数 p 的经验公式是p ≈ fs(Hz)/1000 2~416 kHz 取 12 到 14。阶数不够F3 及以上的极点根本不存在你只能在 12 个根里挑三个F3 经常被 F2 的裙边顶替阶数过高多出来的自由度会去拟合声门脉冲的细结构产生一堆带宽极窄的伪极点而且高阶多项式求根对系数扰动极其敏感。带宽扩展是另一个常用手段。它把多项式系数按 γ 的幂次缩放/* 带宽扩展极点整体向圆心收拢抑制过尖的谱峰 */ for (int i 0; i p; i) a[i] * powf(0.99f, (float)(i 1));当 γ0.99 时一个原本 r0.995带宽约 25 Hz的极点会被拉到 r≈0.985带宽变成约 76 Hz。人耳能分辨的共振峰带宽很少低于 50 Hz太窄的极点基本是数值噪声或者声门谐波压一下能显著减少后面筛选环节的误判。代价是 F1、F2 的位置会有几十赫兹的偏移如果下游要做共振峰轨迹的精确比较γ 尽量取 0.995 以上或者干脆不做。反射系数k[1..p]是最好的单帧体检指标。全都在 (-1, 1) 内说明这帧的自相关是个合法的正定序列出现 |k| ≥ 1多半是能量太低、预加重过头或者帧内有强冲激。我的处理是直接丢帧返回 -2交给上层的平滑逻辑去插值补一个绝不让它进入求根那一步——因为不稳定的多项式牛顿迭代会收敛到单位圆外的根算出来的带宽是负数你在 CSV 里会看到 -300 这种数。4.3 牛顿迭代求复根初值、收敛判据和降阶多项式求根的方法不少这里选牛顿法是因为它不需要求导公式之外的额外信息迭代式简单单根附近二阶收敛12 阶的多项式通常 10 次以内就收敛。缺点是初值敏感、可能不收敛、需要降阶后逐个求。/* polyroot.c对 z^p a[0]z^(p-1) ... a[p-1] 求全部复根 */ int lpc_roots(const double *a, int p, double complex *roots) { double ca[MAX_ORDER]; memcpy(ca, a, sizeof(double) * (size_t)p); int deg p, found 0; while (deg 0) { /* 初值单位圆内均匀撒点半径 0.9 避开 z1 附近的病态区 */ double ang 2.0 * PI * (double)found / (double)deg; double complex z 0.9 * cexp(I * ang); int ok 0; for (int it 0; it 60; it) { double complex f 1.0, df 0.0; for (int k 0; k deg; k) { /* Horner同时求 P 与 P */ df df * z f; f f * z ca[k]; } if (cabs(f) 1e-12) { ok 1; break; } if (cabs(df) 1e-15) break; /* 导数太小换初值 */ double complex dz f / df; z - dz; if (cabs(dz) 1e-10) { ok 1; break; } } if (!ok) return -1; /* 本帧放弃上层补插值 */ roots[found] z; /* 综合除法降阶P p a hrefhttps://download.csdn.net/download/j3210231992/8392611 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p
阅读完成 · 觉得有帮助？
