语音信号处理——线性预测编码LPC

语音信号处理二：干净语音的特征提取：

今天的信号与系统，DSP知识点参考 Spoken Language Processing 第5, 6 章

LPC方程的Durbin算法推导：语音信号数字处理（杨行峻，迟惠生）第四章，数字语音处理（Rabiner）第九章

作业是自己实现语音信号的LPC预测算法：输入一段语音信号，选定不同阶数p，在最小二乘准则下，用自相关法估计预测系数$a_i$，对比重建语音和原始语音的时域&短时频谱差别

自相关法可以用普通的矩阵求逆，和Durbin算法做对比。

语音信号的生成模型：激励-滤波模型：

语音信号的激励部分：声门激励

• 声门：选择激励源
• 声带振动：周期性信号
  • 声带打开：产生激励
  • 声带闭合：信号为0
  • 对应：浊音/元音
  • 声带有个开闭的过程
• 声带松弛：白噪声（与频率无关）
  • 时域：高斯分布
  • 短时频谱：均匀分布
    

语音信号的滤波部分：无损声管模型

• 声管：声道各个器官的抽象模型
  
  $$H(z)=\frac{X(z)}{E(z)}=\frac{1}{1-{\sum }_{k=1}^p{a}_k{z}^{-k}}=\frac{1}{A(z)}$$

• 鼻腔的作用：并行的通道

  • 简化：不考虑口鼻同时打开的情况
    

声源是由声带产生的，声带向声道提供激励信号，这种激励可以是周期性的或非周期性的。当声带处于发声状态（振动）时，会产生有声声音（例如，元音）；而当声带处于无声状态时，会产生无声声音（例如，辅音）。声道可以看作是一个滤波器，它可以对来自声带的激励信号频谱进行整形以产生各种声音。

线性预测编码（LPC）：Linear Predictive Coding

LPC编码的基本思想：

“一个语音取样的现在值可以用若干个语音取样过去值的加权线性组合来逼近”，用过去p个样本点预测当前值：
$$\tilde{x}[n]=\sum _{k=1}^p{a}_{k}x[n-k]$$

在线性组合中的加权系数$a_k$称为预测器系数。通过使实际语音抽样和线性预测抽样之间差值的平方和达到最小值，能够决定唯一的一组预测器系数。

预测误差：
$$e[n]=x[n]-\tilde{x}[n]=x[n]-\sum _{k=1}^p{a}_{k}x[n-k]$$

m个语音信号样本片段的周期延拓：$x_m[n]=x[m+n]$

短时预测误差：
$$E_m=\sum _n{e}_m^2[n]=\sum _n{\left(x_m[n]-{\tilde{x}}_m[n]\right)}^2=\sum _n{\left(x_m[n]-\sum _{j=1}^p{a}_j{x}_m[n-j]\right)}^2$$

线性预测编码通过估计共振峰、剔除它们在语音信号中的作用、估计保留的蜂鸣音强度与频率来分析语音信号。剔除共振峰的过程称为逆滤波，经过这个过程剩余的信号称为残余信号（residue）。

描述峰鸣强度与频率、共鸣峰、残余信号的数字可以保存、发送到其它地方。线性预测编码通过逆向的过程合成语音信号：使用蜂鸣参数与残余信号生成源信号、使用共振峰生成表示声道的滤波器，源信号经过滤波器的处理就得到语音信号。
$$J=E\left[e^2(k)\right]=E\left[{\left(s(k)-\sum _{p=1}^P{a}_{p}s(k-p)\right)}^2\right]$$

LPC特点：

• LPC分析/AR模型

• 把声道抽象为一个全极点模型：
  $$H(z)=\frac{X(z)}{E(z)}=\frac{1}{1-{\sum }_{k=1}^p{a}_k{z}^{-k}}=\frac{1}{A(z)}$$

• p：级联声管个数，LPC分析阶数

• 时域表示：
  $$x[n]=\sum _{k=1}^p{a}_{k}x[n-k]+e[n]$$

LPC分析的正交性原理：

• 预测误差与当前样本正交：
  $$<e_m,x_m^i>=\sum _n{e}_m[n]x_m[n-i]=01\le i\le p$$$$\sum _n{x}_m[n-i]x_m[n]=\sum _{j=1}^p{a}_j\sum _n{x}_m[n-i]x_m[n-j]i=1,2,...,p$$

• 相关系数：${\varphi }_m[i,j]={\sum }_n{x}_m[n-i]x_m[n-j]$

• Yule-Walker Equations：${\sum }_{j=1}^p{a}_j{\varphi }_m[i,j]={\varphi }_m[i,0]i=1,2,...,p$

• 预测误差：
  $$E_m=\sum _n{x}_m^2[n]-\sum _{j=1}^p{a}_j\sum _n{x}_m[n]x_m[n-j]=\varphi [0,0]-\sum _{j=1}^p{a}_j\varphi [0,j]$$

• 预测误差的能量归一化：
  $$e_m[n]=G{u}_m[n]\sum _n{u}_m^2[n]=1E_m=\sum _n{e}_m^2[n]=G^2\sum _n{u}_m^2[n]=G^2$$

LPC方程求解

• Yule-Walker Equations ${\sum }_{j=1}^p{a}_j{\varphi }_m[i,j]={\varphi }_m[i,0]i=1,2,...,p$
  • 本质：矩阵求逆，但是存在有效解法
• 方差矩阵法（Spoken Language Processing , Section 6.3.2.1）
  
• 自相关法
  • 加窗 $x_m[n]=x[m+n]w[n]$
  • 预测误差 $E_m={\sum }_{n=0}^{N+p-1}{e}_m^2[n]$
  • 自相关
    $${\varphi }_m[i,j]=\sum _{n=0}^{N+p-1}{x}_m[n-i]x_m[n-j]=\sum _{n=0}^{N-1-(i-j)}{x}_m[n]x_m[n+i-j]$$$${\varphi }_m[i,j]=R_m[i-j]$$$$R_m[k]=\sum _{n=0}^{N-1-k}{x}_m[n]x_m[n+k]$$

语音生成模型

参考线性预测编码
语音生成模型：

LPC正是基于这个模型的语音生成技术。在该模型中，语音信号是由一个激励信号 $e(k)$经过一个时变的全极点滤波器产生。全极点滤波器的系数取决于所产生的特定声音的声道形状。激励信号$e_k$要么是浊音语音的脉冲序列，要么是无声声音的随机噪声。生成语音信号 $s(k)$可以表示为：
$s(k)={\sum }_{p=1}^P{a}_{p}s(k-p)+e(k)$
其中， P 是滤波器的阶数， $a_p$ 是滤波器的系数。LPC就是在已知 $s(k)$ 的情况下获取 $a_p$ .

求取 $a_p$ 最常用的一个方法就是最小化真实信号与预测信号之间的均方误差（Mean Squared Error, MSE）。MSE函数可以表示为
$$J=E\left[e^2(k)\right]=E\left[{\left(s(k)-\sum _{p=1}^P{a}_{p}s(k-p)\right)}^2\right]$$
然后，计算 J 关于每个滤波器系数的偏导，并令其值等于0，可得（3）：
$$\frac{\partial J}{\partial a_p}=0(3)$$

通过对（3）计算，可以得到（4）：
$$\sum _{u=1}^P{a}_{u}E\left[s(k-p)s(k-u)\right]=E\left[s(k)s(k-u)\right],1\le u\le p(4)$$

其中， $1\le p\le P$ 。用数值 1,2,…,P 分别替换（4）中的变量 p ，我们可以得到 P 个关于滤波器系数的线性方程组，求解该线性方程组，即可得到滤波器系数的解。求解该方程组最常用高效的方法是Levinson-Durbin算法。

Matlab参考：

1. https://wenku.baidu.com/view/25a086fc19e8b8f67c1cb938.html
2. https://blog.csdn.net/kaixinshier/article/details/72142889

上述MSE期望也可以写作：$e(n)=x(n)-\tilde{x}(n)=x(n)-{\sum }_{i=1}^p{a}_{i}x(n-i)$
对$a_i$求偏导可得：
$$\sum _{n}x(n)x(n-j)=\sum _{i=1}^p{a}_i\sum _{n}x(n-i)x(n-j)$$$$E=\sum _n[x(n){]}^2-\sum _{i=1}^p{a}_i\sum _{n}x(n)x(n-i)$$
写成自相关形式（Yule-Walker方程）：
$$R(j)=-\sum _{i=1}^p{a}_{i}R(j-i)1\le j\le p$$

拆写加权式子，即为Toeplize矩阵：
使用Durbin算法来求解Toeplize矩阵，即可计算出滤波器系数$a_i$。

Matlab中自带lpc函数，数学推导过程看《语音信号数字处理(L.R.Rabiner)》。

Matlab程序：

[x,fs] = audioread('1.wav');   %这里读取的双声道信号 x数据，fs采样率
sound(x,fs);    %播放音频

x1 = x(:,1);	
% figure
% plot(t,x1)
n = 200; % n is the length of a frame.  % n行   
p0 = 50; % p0 is the overlap length.    % 在前面填充p0个0。自相关法 两端都需要加P个零取样值，会造成谱估计失真
xx = buffer(x1, n, p0);  % 列：L/(n-p0)，列理解为帧数

m = 8; % (m)th frame of data.
y = xx((m-1)*n+1:m*n); % select a frame of data.
p = 12; % p is the order of the AR model.	阶数
ar = lpc(y,p); % calculate the coefficients of AR model.就是前面的滤波器系数a_i
est_x = filter([0 -ar(2:end )],1,y); % calculate the predicted signal. 建立语音帧的正则方程
err = y - est_x; % calculate the residual signal.

figure
plot(x1)
title('原始信号');
figure
subplot(2,2,1); 
plot(y,'r');
title('原始一帧');
subplot(2,2,2); 
plot(est_x);
title('lpc预测的一帧');
subplot(2,2,3); 
plot(err,'r');
title('残余信号');

效果图：

参考《语音信号处理》实验3-LPC特征提取：

I = audioread('1.wav');   %读入原始语音
I = I(:,1);
plot(I);
title('原始语音波形');%对指定帧位置进行加窗处理
Q = I';
N = 256;    %窗长
Hamm = hamming(N); %加窗
frame = 60;%需要处理的帧位置
M = Q(((frame - 1) * (N / 2) + 1):((frame - 1) * (N / 2) + N));
Frame = M.*Hamm';   %加窗后的语音帧
[B,F,T] = specgram(I,N,N/2,N);
[m,n] = size(B);
for i = 1:m 
    FTframe1(i) = B(i,frame);
end
% P = input('请输入预测器阶数?=?');
P = 5;    % 预测器阶数   改变不同阶数 观察变化
ai = lpc(Frame,P);  %计算lpc系数
LP = filter( [0 - ai(2:end)],1,Frame); %建立语音帧的正则方程
FFTlp = fft(LP);
E = Frame - LP;     % 预测误差
figure
subplot(2,1,1),plot(1:N,Frame,1:N,LP,'-r');grid;
title('原始语音和预测语音波形 ');
subplot(2,1,2)
plot(E);
grid;
title('预测误差');
% pause

fLength(1:2*N) = [M,zeros(1,N)];
Xm = fft(fLength, 2 * N);
X = Xm .* conj(Xm);
Y = fft(X , 2* N);
Rk = Y(1 : N);
PART = sum(ai(2:P+1) .* Rk(1:P));
G = sqrt(sum(Frame.^2) - PART);
A = (FTframe1 - FFTlp(1:length(F')))./FTframe1;
figure
subplot(2,1,1),plot(F',20*log(abs(FTframe1)), F',(20*log(abs(1 ./A))),'-r');
grid;
xlabel('频率/dB');ylabel('幅度');
title('短时谱');
subplot(2,1,2),plot(F',(20*log(abs(G./A))));grid;
xlabel('频率/dB');ylabel('幅度');
title('LPC谱');
% pause

%求出预测误差的倒谱
pitch = fftshift(rceps(E));
M_pitch = fftshift(rceps(Frame));
figure
subplot(2,1,1),plot(M_pitch);grid;
xlabel('语音帧');ylabel('/dB');
title('原始语音帧倒谱');
subplot(2,1,2),plot(pitch);grid;
xlabel('语音帧');ylabel('/dB');
title('预测误差倒谱');
% pause

%画出语谱图
ai1 = lpc(I,P);   %计算原始语音lpc系数
LP1 = filter([0 - ai(2:end)], 1 ,I); % 建立原始语音的正则方程
figure
subplot(2,1,1);
specgram(I,N,N/2,N);
title('原始语音语谱图');
subplot(2,1,2);
specgram(LP1,N,N/2,N);
title('预测语音语谱图')

效果图：