线性预测编码（LPC）笔记

概念：一个语音的抽样能够用过去若干个语音抽样（模板）的线性组合来逼近。通过使实际语音抽样和线性预测抽样之间差值的平方和达到最小，能够决定唯一的一组预测系数。用于语音分析与合成，可估计许多语音基本参数：基音、共振峰、频谱、声道截面积等。

语音信号是由一个激励信号e(k)经过一个时变的全极点滤波器产生的。生成语音信号s(k)表示为：
$$s(n)=\sum _{i=1}^p{a}_{i}s(n-i)+e(n)$$
其中，激励信号e(n)要么是浊音语音的脉冲序列，要么是无声声音的随机噪声。P是滤波器的阶数，$a_i$是滤波器的系数。LPC就是在已知s(n)的情况下获取$a_i$。
$$\tilde{s}(n)=\sum _{i=1}^p{a}_{i}s(n-i)$$

$$e(n)=s(n)-\tilde{s}(n)=s(n)-\sum _{i=1}^p{a}_{i}s(n-i)$$

求取$a_i$最常用的一个方法就是最小化真实信号与预测信号之间的均方误差（MSE）：
$$E=\sum _n{e}^2(n)=\sum _n[s(n)-\tilde{s}(n){]}^2=\sum _n[s(n)-\sum _{i=1}^p{a}_{i}s(n-i){]}^2$$
求$a_i$为多少的情况下，E可取最小值，对$a_i$求偏导：
$$\frac{\partial E}{\partial a_j}=0,j=1,2,...,p$$
可得出：
$$\sum _{n}s(n)s(n-j)=\sum _{i=1}^p{a}_i\sum _{n}s(n-i)s(n-j)$$
根据自相关函数定义Φ：
$${\varphi }_n(i,j)=\sum _{m}s(n-i)s(n-j)$$
所以式子可写成：
$$\varphi (0,j)=\sum _{i=1}^p{a}_i\varphi (i,j)j=1,2,...,p$$
该式表示p个方程构成的方程组，未知数为p个。求解该方程组，就可以得到系统的线性预测系数。

系统的最小均方误差就可以表示为：(sjb，求大佬教这步是怎么推出来的)
$$E=\sum _n[s(n){]}^2-\sum _{i=1}^p{a}_i\sum _{n}s(n)s(n-i)={\varphi }_n(0,0)-\sum _{i=1}^p{a}_i{\varphi }_n(0,i)$$

线性预测之Levinson-Durbin算法

求解这个方程组则用durbin算法进行求解。

拓普利兹矩阵

把Φ变成拓普利兹矩阵R,即第(i,j)个元素为R[i-j]：
$$\sum _{i=1}^p{a}_{i}R[|i-j|]=R[j],1\le j\le p$$
方程的矩阵形式为：
$$Ra=r$$
a和r分别是元素为$a_i$和r[i]=R[i]的列向量。

上式中方程组可做如下变换：
$$R[j]-\sum _{i=1}^p{a}_{i}R[|i-j|]=0,j=1,2,...,p$$

写成矩阵形式为：

$$\left[\begin{array}{ccccc}R[1]& R[0]& R[1]& \cdots & R[p-1]\\ R[2]& R[1]& R[0]& \cdots & R[p-2]\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ R[p]& R[p-1]& R[p-2]& \cdots & R[0]\end{array}\right]\left[\begin{array}{c}1\\ -{\alpha }_1^{(p)}\\ -{\alpha }_2^{(p)}\\ ⋮\\ -{\alpha }_p^{(p)}\end{array}\right]=\left[\begin{array}{c}0\\ 0\\ ⋮\\ 0\end{array}\right]$$

而根据最小均方误差公式，对于一个p阶最佳预测器，其最小均方预测误差为：

$$R[0]-\sum _{i=1}^p{a}_{i}R[i]={\epsilon }^{(p)}$$

写成矩阵的形式为：

$$\left[\begin{array}{ccccc}R[0]& R[1]& R[2]& \cdots & R[p]\end{array}\right]\left[\begin{array}{c}1\\ -{\alpha }_1^{(p)}\\ -{\alpha }_2^{(p)}\\ ⋮\\ -{\alpha }_p^{(p)}\end{array}\right]=\left[\begin{array}{c}{\epsilon }^{(p)}\end{array}\right]$$

再把上面两个矩阵合在一起，写成一个p+1个方程方程组满足p各未知预测器系数和对应的未知均方预测误差${\epsilon }^p$，新的矩阵方程为

$$\left[\begin{array}{ccccc}R[0]& R[1]& R[2]& \cdots & R[p]\\ R[1]& R[0]& R[1]& \cdots & R[p-1]\\ R[2]& R[1]& R[0]& \cdots & R[p-2]\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ R[p]& R[p-1]& R[p-2]& \cdots & R[0]\end{array}\right]\left[\begin{array}{c}1\\ -{\alpha }_1^{(p)}\\ -{\alpha }_2^{(p)}\\ ⋮\\ -{\alpha }_p^{(p)}\end{array}\right]=\left[\begin{array}{c}{\epsilon }^{(p)}\\ 0\\ 0\\ ⋮\\ 0\end{array}\right]$$

上式中构造出的(p+1)x(p+1)阶矩阵也是拓普利兹矩阵。该方程组用Levinson-Durbin算法递归求解。该算法通过在每次迭代中顺序地结合一个新的相关值来实现，并且根据新的相关值和已获得的预测器就能解出下一个高一阶的预测器。

对于任意阶数i，上式中的方程组都可以以矩阵形式表示：
$$R^i{\alpha }^i=e^i$$
我们希望第i阶的解是如何由第i-1阶的解导出的。换而言之是给定${\alpha }^{(i-1)}$,即$R^{(i-1)}{\alpha }^{(i-1)}=e^{(i-1)}$的解，求出$R^{(i)}{\alpha }^{(i)}=e^{(i)}$的解。

durbin算法

步骤：

1、先将矩阵方程$R^{(i-1)}{\alpha }^{(i-1)}=e^{(i-1)}$的扩展形式写为：

$$\left[\begin{array}{ccccc}R[0]& R[1]& R[2]& \cdots & R[i-1]\\ R[1]& R[0]& R[1]& \cdots & R[i-2]\\ R[2]& R[1]& R[0]& \cdots & R[i-3]\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ R[i-1]& R[i-2]& R[i-3]& \cdots & R[0]\end{array}\right]\left[\begin{array}{c}1\\ -{\alpha }_1^{(i-1)}\\ -{\alpha }_2^{(i-1)}\\ ⋮\\ -{\alpha }_{i-1}^{(i-1)}\end{array}\right]=\left[\begin{array}{c}{\epsilon }^{(i-1)}\\ 0\\ 0\\ ⋮\\ 0\end{array}\right]$$
2、将数0附加到向量$a^{(i-1)}$中，并与矩阵$R^{(i)}$相乘，这时得到新的一组i+1个方程：
$$\left[\begin{array}{ccccc}R[0]& R[1]& R[2]& \cdots & R[i]\\ R[1]& R[0]& R[1]& \cdots & R[i-1]\\ R[2]& R[1]& R[0]& \cdots & R[i-2]\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ R[i-1]& R[i-2]& R[i-3]& \cdots & R[1]\\ R[i]& R[i-1]& R[i-2]& \cdots & R[0]\end{array}\right]\left[\begin{array}{c}1\\ -{\alpha }_1^{(i-1)}\\ -{\alpha }_2^{(i-1)}\\ ⋮\\ -{\alpha }_{i-1}^{(i-1)}\\ 0\end{array}\right]=\left[\begin{array}{c}{\epsilon }^{(i-1)}\\ 0\\ 0\\ ⋮\\ 0\\ {\gamma }^{(i-1)}\end{array}\right]$$
上式成立的条件是：
$${\gamma }^{(i-1)}=R[i]-\sum _{j=1}^{i-1}{\alpha }_j^{(i-1)}R[i-j]$$
3、根据拓普利兹矩阵$R^{(i)}$的对称性，方程组可以以相反顺序写出(第一个方程写在最后面，最后一个方程写在最前面，以此类推)，可得:
$$\left[\begin{array}{ccccc}R[0]& R[1]& R[2]& \cdots & R[i]\\ R[1]& R[0]& R[1]& \cdots & R[i-1]\\ R[2]& R[1]& R[0]& \cdots & R[i-2]\\ ⋮& ⋮& ⋮& \cdots & ⋮\\ R[i-1]& R[i-2]& R[i-3]& \cdots & R[1]\\ R[i]& R[i-1]& R[i-2]& \cdots & R[0]\end{array}\right]\left[\begin{array}{c}0\\ -{\alpha }_{i-1}^{(i-1)}\\ -{\alpha }_{i-2}^{(i-1)}\\ ⋮\\ -{\alpha }_1^{(i-1)}\\ 1\end{array}\right]=\left[\begin{array}{c}{\gamma }^{(i-1)}\\ 0\\ 0\\ ⋮\\ 0\\ {\epsilon }^{(i-1)}\end{array}\right]$$
4、将2和3两步中的公式合并可得
$$R^i\cdot \left[\begin{array}{c}\left[\begin{array}{c}1\\ -{\alpha }_1^{(i-1)}\\ -{\alpha }_2^{(i-1)}\\ ⋮\\ -{\alpha }_{i-1}^{(i-1)}\\ 0\end{array}\right]-k_i\left[\begin{array}{c}0\\ -{\alpha }_{i-1}^{(i-1)}\\ -{\alpha }_{i-2}^{(i-1)}\\ ⋮\\ -{\alpha }_1^{(i-1)}\\ 1\end{array}\right]\end{array}\right]=\left[\begin{array}{c}\left[\begin{array}{c}{\epsilon }^{(i-1)}\\ 0\\ 0\\ ⋮\\ 0\\ {\gamma }^{(i-1)}\end{array}\right]-k_i\left[\begin{array}{c}{\gamma }^{(i-1)}\\ 0\\ 0\\ ⋮\\ 0\\ {\epsilon }^{(i-1)}\end{array}\right]\end{array}\right]$$
该式逼近了表达式$R^{(i)}{\alpha }^{(i)}=e^{(i)}$，然后选择${\gamma }^{(i-1)}$使方程右边的向量只有一个非零元素，此时新参数$k_i$必须为:
$$k_i=\frac{{\gamma }^{(i-1)}}{{\epsilon }^{i-1}}=\frac{R[i]-{\sum }_{j=1}^{i-1}{\alpha }_j^{(i-1)}R[i-j]}{{\epsilon }^{i-1}}$$
这样就确保方程右边的向量的最后一个元素为0，则此时第一个元素为:
$${\epsilon }^i={\epsilon }^{i-1}-k_i{\gamma }^{(i-1)}={\epsilon }^{i-1}(1-k_i^2)$$
5、选好${\gamma }^{(i-1)}$后，第i阶的预测系数向量为：
$${\alpha }_i^{(i)}=k_i$$
6、可得系数更新方程组为
$${\alpha }_j^{(i)}={\alpha }_i^{(i-1)}-k_i{\alpha }_{i-j}^{(i-1)},j=1,2,\cdots ,i-1$$

$${\alpha }_i^{(i)}=k_i$$

因此，对于某个特定的阶数p，最佳预测系数为:
$${\alpha }_j={\alpha }_i^{(i)}=k_i,j=1,2,\cdots ,p$$
最终伪代码：

matlab实现

用到lpc函数：

a=lpc(x,n);

这里x为一帧语音信号，n为计算LPC参数的阶数。通常x为240点或256点的数据，n取10-12。

博主表示还是似懂非懂的，略难受。

希望有大牛推荐推荐书籍入个门。