Python语音基础操作--3.5线性预测分析

《语音信号处理试验教程》（梁瑞宇等）的代码主要是Matlab实现的，现在Python比较热门，所以把这个项目大部分内容写成了Python实现，大部分是手动写的。使用CSDN博客查看帮助文件：

Python语音基础操作–2.1语音录制，播放，读取
Python语音基础操作–2.2语音编辑
Python语音基础操作–2.3声强与响度
Python语音基础操作–2.4语音信号生成
Python语音基础操作–3.1语音分帧与加窗
Python语音基础操作–3.2短时时域分析
Python语音基础操作–3.3短时频域分析
Python语音基础操作–3.4倒谱分析与MFCC系数
Python语音基础操作–4.1语音端点检测
Python语音基础操作–4.2基音周期检测
Python语音基础操作–4.3共振峰估计
Python语音基础操作–5.1自适应滤波
Python语音基础操作–5.2谱减法
Python语音基础操作–5.4小波分解
Python语音基础操作–6.1PCM编码
Python语音基础操作–6.2LPC编码
Python语音基础操作–6.3ADPCM编码
Python语音基础操作–7.1帧合并
Python语音基础操作–7.2LPC的语音合成
Python语音基础操作–10.1基于动态时间规整(DTW)的孤立字语音识别试验
Python语音基础操作–10.2隐马尔科夫模型的孤立字识别
Python语音基础操作–11.1矢量量化(VQ)的说话人情感识别
Python语音基础操作–11.2基于GMM的说话人识别模型
Python语音基础操作–12.1基于KNN的情感识别
Python语音基础操作–12.2基于神经网络的情感识别
Python语音基础操作–12.3基于支持向量机SVM的语音情感识别
Python语音基础操作–12.4基于LDA，PCA的语音情感识别

代码可在Github上下载：busyyang/python_sound_open

一个简化的语音模型可以看做一个数字滤波器：
H(z)=S(z)U(z)=G1−∑i=1paizi−1H(z)=\frac{S(z)}{U(z)}=\frac{G}{1-\sum\limits_{i=1}^pa_iz_i^{-1}}H(z)=U(z)S(z)=1−i=1∑paizi−1G

这里s(n)s(n)s(n)和u(n)u(n)u(n)之间的关系可以用差分方程来表示：
s(n)=∑i=1pais(n−1)+Gu(n)s(n)=\sum\limits_{i=1}^pa_is(n-1)+Gu(n)s(n)=i=1∑pais(n−1)+Gu(n)

系统s^(n)=∑i=1pais(n−1)\hat s(n)=\sum\limits_{i=1}^pa_is(n-1)s^(n)=i=1∑pais(n−1)是线性预测器。s^(n)\hat s(n)s^(n)是s(n)s(n)s(n)的估计值。由过去的ppp个值线性组合得到的。ai(i=1,2,...,p)a_i(i=1,2,...,p)ai(i=1,2,...,p)是线性预测系数，预测误差表示为：
e(n)=s(n)−s^(n)=s(n)−∑i=1pais(n−i)e(n)=s(n)-\hat s(n)=s(n)-\sum\limits_{i=1}^pa_is(n-i)e(n)=s(n)−s^(n)=s(n)−i=1∑pais(n−i)

预测的二次方误差为：
E=∑ne2(n)=∑n[s(n)−s^(n)]2=∑n[s(n)−∑i=1pais(n−i)]2E=\sum_ne^2(n)=\sum_n[s(n)-\hat s(n)]^2=\sum_n[s(n)-\sum\limits_{i=1}^pa_is(n-i)]^2E=n∑e2(n)=n∑[s(n)−s^(n)]2=n∑[s(n)−i=1∑pais(n−i)]2

显然预测最准确的话，那么误差EEE要最小。要让参数aia_iai的取值让EEE最小化，要:
∂E∂aj=0,(1⩽j⩽p)\frac{\partial E}{\partial a_j}=0,(1\leqslant j \leqslant p)∂aj∂E=0,(1⩽j⩽p)

∂E∂aj={∑n[s(n)−∑i=1pais(n−i)]2}′=∑n{2(s(n)−∑i=1pais(n−i))[s(n)−∑i=1pais(n−i)]′}=∑n{2(s(n)−∑i=1pais(n−i))(1−s(n−j))}=2∑n{s(n)−∑i=1paix(n−i)+s(n)s(n−j)−∑i=1pais(n−i)s(n−j)}=2∑ns(n)s(n−j)−2∑i=1pai∑ns(n−i)s(n−j)\begin{array}{ll} \frac{\partial E}{\partial a_j}&=\{\sum_n[s(n)-\sum\limits_{i=1}^pa_is(n-i)]^2\}'\\ &=\sum_n\{2(s(n)-\sum\limits_{i=1}^pa_is(n-i))[s(n)-\sum\limits_{i=1}^pa_is(n-i)]'\}\\ &=\sum_n\{2(s(n)-\sum\limits_{i=1}^pa_is(n-i))(1-s(n-j))\}\\ &=2\sum_n\{s(n)-\sum_{i=1}^pa_ix(n-i)+s(n)s(n-j)-\sum\limits_{i=1}^pa_is(n-i)s(n-j)\}\\ &=2\sum_ns(n)s(n-j)-2\sum\limits_{i=1}^pa_i\sum_ns(n-i)s(n-j) \end{array}∂aj∂E={∑n[s(n)−i=1∑pais(n−i)]2}′=∑n{2(s(n)−i=1∑pais(n−i))[s(n)−i=1∑pais(n−i)]′}=∑n{2(s(n)−i=1∑pais(n−i))(1−s(n−j))}=2∑n{s(n)−∑i=1paix(n−i)+s(n)s(n−j)−i=1∑pais(n−i)s(n−j)}=2∑ns(n)s(n−j)−2i=1∑pai∑ns(n−i)s(n−j)

如果定义：ϕ(j,i)=∑ns(n−i)s(n−j)\phi(j,i)=\sum_ns(n-i)s(n-j)ϕ(j,i)=∑ns(n−i)s(n−j)，那么∂E∂aj=0\frac{\partial E}{\partial a_j}=0∂aj∂E=0就可以写成：
ϕ(j,0)=∑i=1paiϕ(j,i)(求解方程)\phi(j,0)=\sum\limits_{i=1}^pa_i\phi(j,i) \tag{求解方程}ϕ(j,0)=i=1∑paiϕ(j,i)(求解方程)

带入误差方程可以有：
E=ϕ(0,0)−∑i=1paiϕ(0,i)(最小误差方程)E=\phi(0,0)-\sum_{i=1}^pa_i\phi(0,i)\tag{最小误差方程}E=ϕ(0,0)−i=1∑paiϕ(0,i)(最小误差方程)

所以，最小误差是由一个固定分量ϕ(0,0)\phi(0,0)ϕ(0,0)和依赖于预测系数的分量∑i=1paiϕ(0,i)\sum_{i=1}^pa_i\phi(0,i)∑i=1paiϕ(0,i)组成的。

线性预测的自相关解法

自相关解法是s(n)s(n)s(n)在0⩽n⩽N−10 \leqslant n \leqslant N-10⩽n⩽N−1以外的值都是零，等同于假设了s(n)s(n)s(n)经过了有限长度的窗，就可以用p个方程来解有p个未知数的方程组了。定义：
r(j)=∑n=0N−1s(n)s(n−j),1⩽j⩽pr(j)=\sum_{n=0}^{N-1}s(n)s(n-j),1\leqslant j \leqslant pr(j)=n=0∑N−1s(n)s(n−j),1⩽j⩽p

那么ϕ(j,i)\phi(j,i)ϕ(j,i)等效于r(j−i)r(j-i)r(j−i),由于自相关是偶函数，所以有：
ϕ(j,i)=r(∣j−i∣)\phi(j,i)=r(|j-i|)ϕ(j,i)=r(∣j−i∣)

所以求解方程为：
r(j)=∑i=1pair(∣j−i∣),1⩽j⩽pr(j)=\sum_{i=1}^pa_ir(|j-i|),1\leqslant j \leqslant pr(j)=i=1∑pair(∣j−i∣),1⩽j⩽p

最小误差方程为：
E=r(0)−∑i=1pair(i)E=r(0)-\sum_{i=1}^pa_ir(i)E=r(0)−i=1∑pair(i)

将ϕ(j,i)=r(∣j−i∣)\phi(j,i)=r(|j-i|)ϕ(j,i)=r(∣j−i∣)展开有：
[r(0)r(1)r(2)...r(p−1)r(1)r(0)r(1)...r(p−2)r(2)r(1)r(0)...r(p−3)...............r(p−1)r(p−2)r(p−3)...r(0)][a1a2a3...ap]=[r(1)r(2)r(3)...r(p)]\begin{bmatrix}r(0)&r(1)&r(2)&...&r(p-1)\\r(1)&r(0)&r(1)&...&r(p-2)\\r(2)&r(1)&r(0)&...&r(p-3)\\...&...&...&...&...\\r(p-1)&r(p-2)&r(p-3)&...&r(0)\end{bmatrix}\begin{bmatrix}a_1\\a_2\\a_3\\...\\a_p\end{bmatrix}=\begin{bmatrix}r(1)\\r(2)\\r(3)\\...\\r(p)\end{bmatrix}⎣⎢⎢⎢⎢⎡r(0)r(1)r(2)...r(p−1)r(1)r(0)r(1)...r(p−2)r(2)r(1)r(0)...r(p−3)...............r(p−1)r(p−2)r(p−3)...r(0)⎦⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎡a1a2a3...ap⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎡r(1)r(2)r(3)...r(p)⎦⎥⎥⎥⎥⎤

这种沿着主对角对称，并任何一条与主对角线平行的线上的值都相等的矩阵叫Toeplitz矩阵，这种方程叫Yule-Walker方程。可以采用递归方法求解。
计算的步骤：

当i=0i=0i=0时，E0=r(0),a0=1E_0=r(0),a_0=1E0=r(0),a0=1;
对于第iii次递归(i=1,2,...,pi=1,2,...,pi=1,2,...,p):
ki=1Ei−1[r(i)−∑j=1i−1aji−1r(j−i)]k_i=\frac{1}{E_{i-1}}[r(i)-\sum_{j=1}^{i-1}a_j^{i-1}r(j-i)]ki=Ei−11[r(i)−j=1∑i−1aji−1r(j−i)]

iai(i)=kiia_i^{(i)}=k_iiai(i)=ki

对于j=1j=1j=1到i−1i-1i−1:
aj(i)=aji−1−kiai−j(i−1)a_j^{(i)}=a_j^{i-1}-k_ia_{i-j}^{(i-1)}aj(i)=aji−1−kiai−j(i−1)

Ei=(1−ki2)Ei−1E_i=(1-k_i^2)E_{i-1}Ei=(1−ki2)Ei−1

增益G为：

G=EpG=\sqrt{E_p}G=Ep

递归得到：
Ep=r(0)∏i=1p(1−ki2)E_p=r(0)\prod_{i=1}^p(1-k_i^2)Ep=r(0)i=1∏p(1−ki2)

其他参数

预测误差及其自相关

预测误差为：
e(n)=s(n)−∑i=1pais(n−i)e(n)=s(n)-\sum_{i=1}^pa_is(n-i)e(n)=s(n)−i=1∑pais(n−i)

预测误差的自相关为：
Re(m)=∑n=0N−1−me(n)e(n+m)R_e(m)=\sum_{n=0}^{N-1-m}e(n)e(n+m)Re(m)=n=0∑N−1−me(n)e(n+m)

反射系数和声道面积

反射系数kik_iki在低速率语音编码/语音合成/语音识别和说话人识别等许多领域都是非常重要的特征。
aj(i)=aji−1−kiai−j(i−1)ai−j(i)=ai−ji−1−kiaj(i−1)\begin{array}{ll} a_j^{(i)}=a_j^{i-1}-k_ia_{i-j}^{(i-1)}\\ a_{i-j}^{(i)}=a_{i-j}^{i-1}-k_ia_{j}^{(i-1)} \end{array}aj(i)=aji−1−kiai−j(i−1)ai−j(i)=ai−ji−1−kiaj(i−1)

可以解出：
aji−1=(aj(i)+aj(i)ai−j(i))/(1−ki2)a_{j}^{i-1}=(a_j^{(i)}+a_j^{(i)}a_{i-j}^{(i)})/(1-k_i^2)aji−1=(aj(i)+aj(i)ai−j(i))/(1−ki2)

由线性预测系数可以反解出反射系数。声道可以被模拟成一系列截面积不等的无损声道级联。反射系数反应了声波在各管道界面处的反射量：
ki=Ai+1−AiAi+1+Aik_i=\frac{A_{i+1}-A_i}{A_{i+1}+A_i}ki=Ai+1+AiAi+1−Ai

同样可以反解出：
AiAi+1=1−ki1+ki\frac{A_i}{A_{i+1}}=\frac{1-k_i}{1+k_i}Ai+1Ai=1+ki1−ki

线性预测频谱

一帧信号x(n)x(n)x(n)模型可以转化为一个ppp阶的线性预测模型。当z=ejwz=e^{jw}z=ejw时，能得到线性预测系数的频谱：
H(ejw)=11−∑n=1panz−jwnH(e^{jw})=\frac{1}{1-\sum_{n=1}^pa_nz^{-jwn}}H(ejw)=1−∑n=1panz−jwn1

不考虑激励和辐射时，H(ejw)H(e^{jw})H(ejw)就是X(ejw)X(e^{jw})X(ejw)的包络线。

线性预测倒谱

语音信号的倒谱可以通过对信号做傅里叶变换，取模的对数，再求傅里叶逆变换得到。由于频率响应H(eiω)H(e^{iω})H(eiω)反映声道的频率响应和被分析信号的谱包络，因此用log⁡∣H(eiw)∣\log|H(e^{iw})|log∣H(eiw)∣做傅里叶拟变换得到线性预测倒谱系数（Linear Prediction Cepstrum Coeffient, LPCC）。

更根据同态处理法有：H^(z)=log⁡H(z)\hat H(z)=\log H(z)H^(z)=logH(z),H^(z)\hat H(z)H^(z)可以写成级数展开的形式：
H^(z)=∑n=1+∞h^(n)z−n\hat H(z)=\sum_{n=1}^{+\infin}\hat h(n)z^{-n}H^(z)=n=1∑+∞h^(n)z−n

对上式同时对z−1z^{-1}z−1求导：
∂∂z−111−∑i=1paiz−i=∂∂z−1∑n=1+∞h^(n)z−n\frac{\partial }{\partial z^{-1}}\frac{1}{1-\sum\limits_{i=1}^pa_iz^{-i}}=\frac{\partial }{\partial z^{-1}}\sum\limits_{n=1}^{+\infin}\hat h(n)z^{-n}∂z−1∂1−i=1∑paiz−i1=∂z−1∂n=1∑+∞h^(n)z−n

解为：
(1−∑i=1paiz−1)∑n=1+∞nh^(n)z−n+1=∑n=1+∞iaiz−i+1(1-\sum\limits_{i=1}^pa_iz^{-1})\sum\limits_{n=1}^{+\infin}n\hat h(n)z^{-n+1}=\sum\limits_{n=1}^{+\infin}ia_iz^{-i+1}(1−i=1∑paiz−1)n=1∑+∞nh^(n)z−n+1=n=1∑+∞iaiz−i+1

令上市两边的zzz各次幂的系数分别相等，得到h^(n)\hat h(n)h^(n)和aia_iai之间的递推关系：
{h^(1)=a1h^(n)=an+∑i=1n−1(1−in)aih^(n−i)1<n⩽ph^(n)=∑i=1n−1(1−in)aih^(n−i)n>p\left \{\begin{array}{ll} \hat h(1)=a_1&\\ \hat h(n)=a_n+\sum_{i=1}^{n-1}(1-\frac{i}{n})a_i\hat h(n-i)&1<n\leqslant p\\ \hat h(n)=\sum_{i=1}^{n-1}(1-\frac{i}{n})a_i\hat h(n-i)&n>p \end{array}\right.⎩⎨⎧h^(1)=a1h^(n)=an+∑i=1n−1(1−ni)aih^(n−i)h^(n)=∑i=1n−1(1−ni)aih^(n−i)1<n⩽pn>p

# lpc.py
import numpy as npdef lpc_coeff(s, p):""":param s: 一帧数据:param p: 线性预测的阶数:return:"""n = len(s)# 计算自相关函数Rp = np.zeros(p)for i in range(p):Rp[i] = np.sum(np.multiply(s[i + 1:n], s[:n - i - 1]))Rp0 = np.matmul(s, s.T)Ep = np.zeros((p, 1))k = np.zeros((p, 1))a = np.zeros((p, p))# 处理i=0的情况Ep0 = Rp0k[0] = Rp[0] / Rp0a[0, 0] = k[0]Ep[0] = (1 - k[0] * k[0]) * Ep0# i=1开始，递归计算if p > 1:for i in range(1, p):k[i] = (Rp[i] - np.sum(np.multiply(a[:i, i - 1], Rp[i - 1::-1]))) / Ep[i - 1]a[i, i] = k[i]Ep[i] = (1 - k[i] * k[i]) * Ep[i - 1]for j in range(i - 1, -1, -1):a[j, i] = a[j, i - 1] - k[i] * a[i - j - 1, i - 1]ar = np.zeros(p + 1)ar[0] = 1ar[1:] = -a[:, p - 1]G = np.sqrt(Ep[p - 1])return ar, Gdef lpcff(ar, npp=None):""":param ar: 线性预测系数:param npp: FFT阶数:return:"""p1 = ar.shape[0]if npp is None:npp = p1 - 1ff = 1 / np.fft.fft(ar, 2 * npp + 2)return ff[:len(ff) // 2]def lpc_lpccm(ar, n_lpc, n_lpcc):lpcc = np.zeros(n_lpcc)lpcc[0] = ar[0]  # 计算n=1的lpccfor n in range(1, n_lpc):  # 计算n=2,,p的lpcclpcc[n] = ar[n]for l in range(n - 1):lpcc[n] += ar[l] * lpcc[n - 1] * (n - l) / nfor n in range(n_lpc, n_lpcc):  # 计算n>p的lpcclpcc[n] = 0for l in range(n_lpc):lpcc[n] += ar[l] * lpcc[n - 1] * (n - l) / nreturn -lpcc

from chapter2_基础.soundBase import *
from chapter3_分析实验.lpc import *
from scipy.signal import lfilterdata, fs = soundBase('C3_5_y.wav').audioread()
L = 240
x = data[8000:8000 + L]
x = (x - np.mean(x)) / np.std(x)
p = 12
ar, G = lpc_coeff(x, p)
b = np.zeros(p + 1)
b[0] = 1
b[1:] = -ar[1:]
est_x = lfilter(b, 1, x)
plt.subplot(2, 1, 1)
plt.plot(x)
plt.subplot(2, 1, 2)
plt.plot(est_x)
plt.savefig('images/lpc.png')
plt.close()

from chapter2_基础.soundBase import *
from chapter3_分析实验.lpc import *data, fs = soundBase('C3_5_y.wav').audioread()
L = 240
p = 12
x = data[8000:8000 + L]
ar, G = lpc_coeff(x, p)
nfft = 512
W2 = nfft // 2
m = np.array([i for i in range(W2)])
Y = np.fft.fft(x, nfft)
Y1 = lpcff(ar, W2)
plt.subplot(2, 1, 1)
plt.plot(x)
plt.subplot(2, 1, 2)
plt.plot(m, 20 * np.log(np.abs(Y[m])))
plt.plot(m, 20 * np.log(np.abs(Y1[m])))
plt.savefig('images/lpcff.png')
plt.close()

from chapter2_基础.soundBase import *
from chapter3_分析实验.lpc import *
from chapter3_分析实验.倒谱计算 import *data, fs = soundBase('C3_5_y.wav').audioread()
L = 240
p = 12
x = data[8000:8000 + L]
ar, G = lpc_coeff(x, p)lpcc1 = lpc_lpccm(ar, p, p)
lpcc2 = rcceps(ar)
plt.subplot(2, 1, 1)
plt.plot(lpcc1)
plt.subplot(2, 1, 2)
plt.plot(lpcc2)
plt.show()