手写一个 LPC 语音 Codec从一帧真实音频看编码到解码的全过程关键词LPC、线性预测、Levinson-Durbin、Yule-Walker、反射系数、全极点合成、有损压缩适用读者想真正算一遍搞清楚 LPC 怎么工作的同学 / 工程师 / 教程作者摘要本文不堆公式而是取一段真实录音里的 9 个采样一步一步把 LPC 语音编码器的每一个中间结果都算出来再原样解回去。你会发现LPC 本身是一个无损的可逆变换真正产生压缩和失真的是最后那一步对残差做 4 bit 粗量化。全文数据前后完全一致可直接复算。1. LPC 到底是什么一句话版语音是短时相关的信号现在的采样很大程度上能用过去几个采样加权预测出来。预测: x_hat[n] a1·x[n-1] a2·x[n-2] a3·x[n-3] 残差: e[n] x[n] - x_hat[n]a预测器描述声道/共振峰e残差是预测剩下的激励又小又像噪声。编码端只传a或其等价形式 K 残差 e解码端用 a 把 e 合成回 x。因为 e 很小粗量化也听不出于是压缩就来了。历史注自相关法来自 Yule(1927)/Walker(1931) 的 Yule-Walker 方程高效解法 Levinson(1947) 递推、Durbin(1959) 应用于 AR 模型反射系数 K 后来被 Itakura-Saito 用于语音PARCOR/LPC。2. 实验素材真实音频里的一帧取一段 16 kHz 单声道语音myvoice.wav挑一帧中等响度的 9 个采样3 个历史 6 个当前帧预测阶数 P3x [0.03650, 0.03955, 0.04218, | 0.04437, 0.04623, 0.04779, 0.04901, 0.04968, 0.04971] └── 3 历史 ──┘ └──────── 6 帧 ────────┘为方便手算选了幅度 ~0.04 的片段数字小、好读换任意其它帧算法完全一样。3. 编码端从音频到两串整数STAGE1 — 原始音频就是上面的x9 点浮点来自.wav的int16 / 32768。STAGE2 — 自相关r[k] Σ x[n]·x[n-k]只算 P1 4 个k0…3r[0] 0.04437²0.04623²...0.04971² 0.01373 (能量) r[1] 0.04437·0.04623 ... 0.04968·0.04971 0.01338 r[2] ... 0.01290 r[3] ... 0.01231 r[0] 加 0.2% 噪声地板: 0.01373 × 1.002 0.01376 (防除零)归一化看相关度r[1]/r[0]0.97, r[2]/r[0]0.94, r[3]/r[0]0.89—— 相邻采样 97% 相关冗余很大这正是 LPC 能压缩的前提。STAGE3 — Levinson-Durbin 递推 → K 和 a逐阶手算E 初值 r[0]每阶乘(1-k²)i0: accr[1]0.01338, k0.01338/0.013760.9500 → K[0]0.95(夹) E0.01376·(1-0.95²)0.00134 i1: accr[2]-a[0]·r[1]0.01290-0.95·0.013380.00019 k0.00019/0.001340.1429 → K[1]0.1429 i2: accr[3]-a[0]·r[2]-a[1]·r[1]-0.00011 k-0.00011/0.00131-0.0809 → K[2]-0.0809 最终: K [0.95, 0.1429, -0.0809] 反射系数 a [0.8258, 0.2087, -0.0809] 预测器(直接型)K[0]0.95表示只用前 1 个采样就能抓 95% 的相关性。|K|1保证滤波器稳定。STAGE4 — 量化 K准备传输K 在 [-1,1]映射到 10 bit 整数 0…1023码 round((K1)/2 × 1023)K[0]0.95 → 码 997 → 反量化 0.9492 → 夹到 ±0.90 → 0.90 K[1]0.1429→ 码 585 → 反量化 0.1437 → 0.1437 K[2]-0.0809→码 470 → 反量化 -0.0811 → -0.0811 Kq [0.90, 0.1437, -0.0811] ← 真正传输的反射系数STAGE5 — Kq → a_dec阶升 step-up解码端也会做同样一步所以编码端也用同一份保证两端预测器一致a_dec refl_to_lpc(Kq) [0.7823, 0.2062, -0.0811]注意 a_dec ≠ raw a因为 K[0] 被夹 0.95→0.90。STAGE6 — 用 a_dec 算残差e_dec[n] x[n] - (0.7823·x[n-1] 0.2062·x[n-2] - 0.0811·x[n-3]) e_dec [0.00618, 0.00603, 0.00589, 0.00569, 0.00524, 0.00461]对比信号幅度 ~0.047残差只剩 ~0.006。信号 RMS 0.04784 → 残差 RMS 0.00563能量压到 1/72≈18.6 dB。注若用未量化的 raw a残差 RMS 仅 0.00352能量压到 1/185 ≈ 22.7 dBK 量化让增益略降代价极小。STAGE7 — 帧增益 rms 残差 4 bit 量化帧 rms 0.00563 残差/rms 的 4bit 码 [10, 10, 9, 9, 9, 9] 重建 rq [0.00751, 0.00751, 0.00451, 0.00451, 0.00451, 0.00451]残差码全是 9~104 bit 范围 0~15中点 7.5说明残差小、靠近 04 bit 足够。STAGE8 — 打包Kq: 3 × 10 bit 30 bit 残差: 6 × 4 bit 24 bit 合计 54 bit 6.75 字节 (原 PCM: 6 采样 × 16 bit 96 bit)4. 传输线路上只跑这两串整数Kq码流: [997, 585, 470] 残差码流: [10, 10, 9, 9, 9, 9]5. 解码端从两串整数反算回音频STAGE9 — 读码恢复参数反算起点读 Kq码 → dq → 夹 → a_dec [0.7823, 0.2062, -0.0811] (与编码端完全相同 ✓) 读残差码 → rq [0.00751, 0.00751, 0.00451, 0.00451, 0.00451, 0.00451]STAGE10 — 残差反量化就是 STAGE7 的逆rq 码/15 × 8 - 4再 × rms。STAGE11 — 全极点合成反算核心LPC 的逆运算用合成出来的过去做反馈s[n] rq[n] 0.7823·s[n-1] 0.2062·s[n-2] - 0.0811·s[n-3]历史用移位寄存器首帧用那 3 个历史采样。逐步算出s [0.04570, 0.04875, 0.04865, 0.04891, 0.04885, 0.04886]与原始x[3..8]只差一点点见下节误差分析。STAGE12 — 转 PCMPCM int16 round(s × 32767), 限幅 [-1, 1] 还原 PCM [1497, 1597, 1594, 1603, 1601, 1601]6. 结果误差与哪一步有损原始 PCM [1454, 1515, 1566, 1606, 1628, 1629] 还原 PCM [1497, 1597, 1594, 1603, 1601, 1601] 还原误差 RMS 0.001308 (原信号 RMS 0.04784 → 相对误差 2.7%)逐 STAGE 定性STAGE操作有损1–3x→r→K,a浮点运算无损4K→10bit量化夹有量化但两端一致→不添重建误差5–6Kq→a_dec、算 e_dec无损7残差→4bit量化有损 ← 误差 0.001308 全来自这里8打包无损9–11读码、合成无损逆运算12float→int16末端格式量化≈0.00003可忽略结论只有 STAGE7残差 4 bit 量化是真正有损的步骤。证明把 STAGE7 的rq换成理想残差e_dec不量化STAGE11 还原误差 0——LPC 分析/合成严格互逆。7. 为什么能压缩为什么有损和 MP3/Opus 什么关系能压缩LPC 去掉了短时相关性把能量压到 1/72本帧。剩下的残差幅度只有原信号的 1/8.5、且像噪声于是对它做粗量化4 bit代价很小。原 16 bit/采样 → 本帧 9.0 bit/采样小帧开销占比大真实 codec 用 frame160、P8 时摊薄头开销可到4.6 bit/采样≈3.5× 压缩全文件相关度 0.9855听感正常。有损压缩来自残差少 bit而少 bit 量化 失真。有损是为换压缩比主动付出的不是算法缺陷。和工业 codec 的关系MP3/AAC/Opus(SILK) 同一思路——去冗余LPC 或 MDCT 对小且像噪声的残差/系数做感知量化。区别只在它们阶数更高、量化更精细、还加了心理声学模型。本文的 mini-codec 就是它的最小可运行内核。8. 结论用真实音频的一帧我们完整走通了编码: x ─自相关→ K ─量化→ Kq码[997,585,470] 残差─量化→ [10,10,9,9,9,9] 解码: Kq码→a_dec, 残差码→rq ─全极点合成→ s ─×32767→ PCMLPC 变换本身无损可逆压缩与失真100% 来自残差 4 bit 量化STAGE7整条链数据前后完全对上可逐数复算。附录 A本帧完整数据流一览x [0.03650, 0.03955, 0.04218, 0.04437, 0.04623, 0.04779, 0.04901, 0.04968, 0.04971] r [0.01376, 0.01338, 0.01290, 0.01231] K [0.95, 0.1429, -0.0809] Kq码 [997, 585, 470] → Kq [0.90, 0.1437, -0.0811] a_dec [0.7823, 0.2062, -0.0811] e_dec [0.00618, 0.00603, 0.00589, 0.00569, 0.00524, 0.00461] rms 0.00563 残差4bit [10, 10, 9, 9, 9, 9] → rq [0.00751, 0.00751, 0.00451, 0.00451, 0.00451, 0.00451] 传输 [997,585,470] [10,10,9,9,9,9] 还原 s [0.04570, 0.04875, 0.04865, 0.04891, 0.04885, 0.04886] 还原 PCM [1497, 1597, 1594, 1603, 1601, 1601] 误差 RMS 0.001308附录 B核心代码片段Python / numpydeflevinson(r,P):anp.zeros(P);Er[0]1e-9;Knp.zeros(P)foriinrange(P):accr[i1]-sum(a[j]*r[i-j]forjinrange(i))knp.clip(acc/E,-0.95,0.95);K[i]kforjinrange(i//2):ta[j];a[j]t-k*a[i-1-j];a[i-1-j]a[i-1-j]-k*tifi1:mi//2;a[m]a[m]-k*a[m]a[i]k;E*(1-k*k)returna,E,Kdefrefl_to_lpc(K):# K - 直接型预测器 aPlen(K);anp.zeros(P1);a[0]1.0forminrange(1,P1):kmK[m-1];newa.copy()foriinrange(1,m):new[i]a[i]-km*a[m-i]new[m]km;anewreturna[1:]# 编码(单帧): K 量化 残差 4bitKqnp.clip(dq_lpc(q_lpc(K)),-0.90,0.90);a_decrefl_to_lpc(Kq)exw[P:]-[sum(a_dec[k]*xw[Pn-1-k]forkinrange(P))forninrange(frame)]rqdq_res(q_res(e/np.sqrt(np.mean(e*e)),4),4)*np.sqrt(np.mean(e*e))# 解码(单帧): 读码 - 合成a_decrefl_to_lpc(np.clip(dq_lpc(qK),-0.90,0.90))sbuflist(history)forninrange(frame):srq[n]sum(a_dec[k]*sbuf[P-1-k]forkinrange(P))sbuf.append(s);sbuf.pop(0)全文数据由myvoice.wav16 kHz 单声道在 P3、frame6 的真实帧start4377逐数复算得到。