語音信號處理——線性預測編碼LPC「建議收藏」

語音信號處理二：干凈語音的特征提?。?/h2>

今天的信號與系統，DSP知識點參考 Spoken Language Processing 第5, 6 章

LPC方程的Durbin算法推導：語音信號數字處理（楊行峻，遲惠生）第四章，數字語音處理（Rabiner）第九章

作業是自己實現語音信號的LPC預測算法：輸入一段語音信號，選定不同階數p，在最小二乘準則下，用自相關法估計預測系數 a i a_i ai，對比重建語音和原始語音的時域&短時頻譜差別

自相關法可以用普通的矩陣求逆，和Durbin算法做對比。

語音信號的生成模型：激勵-濾波模型：

語音信號的激勵部分：聲門激勵

• 聲門：選擇激勵源
• 聲帶振動：周期性信號
  • 聲帶打開：產生激勵
  • 聲帶閉合：信號為0
  • 對應：濁音/元音
  • 聲帶有個開閉的過程
• 聲帶松弛：白噪聲（與頻率無關）
  • 時域：高斯分布
  • 短時頻譜：均勻分布
    

語音信號的濾波部分：無損聲管模型

• 聲管：聲道各個器官的抽象模型
  
  H ( z ) = X ( z ) E ( z ) = 1 1 - ∑ k = 1 p a k z - k = 1 A ( z ) H(z)=\frac{X(z)}{E(z)}=\frac{1}{1-\sum_{k=1}^pa_kz^{-k}}=\frac{1}{A(z)} H(z)=E(z)X(z)=1-∑k=1pakz-k1=A(z)1

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

  • 簡化：不考慮口鼻同時打開的情況
    

聲源是由聲帶產生的，聲帶向聲道提供激勵信號，這種激勵可以是周期性的或非周期性的。當聲帶處于發聲狀態（振動）時，會產生有聲聲音（例如，元音）；而當聲帶處于無聲狀態時，會產生無聲聲音（例如，輔音）。聲道可以看作是一個濾波器，它可以對來自聲帶的激勵信號頻譜進行整形以產生各種聲音。

線性預測編碼（LPC）：Linear Predictive Coding

LPC編碼的基本思想：

“一個語音取樣的現在值可以用若干個語音取樣過去值的加權線性組合來逼近”，用過去p個樣本點預測當前值：
x ~ [ n ] = ∑ k = 1 p a k x [ n - k ] \widetilde{x}[n]=\sum_{k=1}^pa_kx[n-k] x
[n]=k=1∑pakx[n-k]

在線性組合中的加權系數 a k a_k ak稱為預測器系數。通過使實際語音抽樣和線性預測抽樣之間差值的平方和達到最小值，能夠決定唯一的一組預測器系數。

預測誤差：
e [ n ] = x [ n ] - x ~ [ n ] = x [ n ] - ∑ k = 1 p a k x [ n - k ] e[n]=x[n]-\widetilde{x}[n] = x[n] – \sum_{k=1}^pa_kx[n-k] e[n]=x[n]-x
[n]=x[n]-k=1∑pakx[n-k]

m個語音信號樣本片段的周期延拓： x m [ n ] = x [ m + n ] x_m[n] = x[m+n] xm[n]=x[m+n]

短時預測誤差：
E m = ∑ n e m 2 [ n ] = ∑ n ( x m [ n ] - x ~ m [ n ] ) 2 = ∑ n ( x m [ n ] - ∑ j = 1 p a j x m [ n - j ] ) 2 E_m=\sum_ne_m^2[n] = \sum_n\left(x_m[n]-\widetilde{x}_ m[n]\right)^2=\sum_n\left(x_m[n]-\sum_{j=1}^pa_jx_m[n-j]\right)^2 Em=n∑em2[n]=n∑(xm[n]-x
m[n])2=n∑(xm[n]-j=1∑pajxm[n-j])2

線性預測編碼通過估計共振峰、剔除它們在語音信號中的作用、估計保留的蜂鳴音強度與頻率來分析語音信號。剔除共振峰的過程稱為逆濾波，經過這個過程剩余的信號稱為殘余信號（residue）。

描述峰鳴強度與頻率、共鳴峰、殘余信號的數字可以保存、發送到其它地方。線性預測編碼通過逆向的過程合成語音信號：使用蜂鳴參數與殘余信號生成源信號、使用共振峰生成表示聲道的濾波器，源信號經過濾波器的處理就得到語音信號。
J = E [ e 2 ( k ) ] = E [ ( s ( k ) - ∑ p = 1 P a p s ( k - p ) ) 2 ] 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=E[e2(k)]=E(s(k)-p=1∑Paps(k-p))2

LPC特點：

• LPC分析/AR模型

• 把聲道抽象為一個全極點模型：
  H ( z ) = X ( z ) E ( z ) = 1 1 - ∑ k = 1 p a k z - k = 1 A ( z ) H(z)=\frac{X(z)}{E(z)}=\frac{1}{1-\sum_{k=1}^pa_kz^{-k}}=\frac{1}{A(z)} H(z)=E(z)X(z)=1-∑k=1pakz-k1=A(z)1

• p：級聯聲管個數，LPC分析階數

• 時域表示：
  x [ n ] = ∑ k = 1 p a k x [ n - k ] + e [ n ] x[n] = \sum_{k=1}^pa_kx[n-k]+e[n] x[n]=k=1∑pakx[n-k]+e[n]

LPC分析的正交性原理：

• 預測誤差與當前樣本正交：
  &lt; e m , x m i &gt; = ∑ n e m [ n ] x m [ n - i ] = 0 1 ≤ i ≤ p &lt;e_m,x_m^i&gt;=\sum_ne_m[n]x_m[n-i] =0 \quad\quad1≤i≤p <em,xmi>=n∑em[n]xm[n-i]=01≤i≤p ∑ n x m [ n - i ] x m [ n ] = ∑ j = 1 p a j ∑ n x m [ n - i ] x m [ n - j ] i = 1 , 2 , . . . , p \sum_nx_m[n-i]x_m[n]=\sum_{j=1}^pa_j\sum_nx_m[n-i]x_m[n-j]\quad\quad i=1,2,…,p n∑xm[n-i]xm[n]=j=1∑pajn∑xm[n-i]xm[n-j]i=1,2,...,p

• 相關系數： φ m [ i , j ] = ∑ n x m [ n - i ] x m [ n - j ] \phi_m[i,j]=\sum_nx_m[n-i]x_m[n-j] φm[i,j]=∑nxm[n-i]xm[n-j]

• Yule-Walker Equations： ∑ j = 1 p a j φ m [ i , j ] = φ m [ i , 0 ] i = 1 , 2 , . . . , p \sum_{j=1}^pa_j\phi_m[i,j]=\phi_m[i,0]\quad i=1,2,…,p ∑j=1pajφm[i,j]=φm[i,0]i=1,2,...,p

• 預測誤差：
  E m = ∑ n x m 2 [ n ] - ∑ j = 1 p a j ∑ n x m [ n ] x m [ n - j ] = φ [ 0 , 0 ] - ∑ j = 1 p a j φ [ 0 , j ] E_m=\sum_nx_m^2[n]-\sum_{j=1}^pa_j\sum_nx_m[n]x_m[n-j]=\phi[0,0]-\sum_{j=1}^pa_j\phi[0,j] Em=n∑xm2[n]-j=1∑pajn∑xm[n]xm[n-j]=φ[0,0]-j=1∑pajφ[0,j]

• 預測誤差的能量歸一化：
  e m [ n ] = G u m [ n ] ∑ n u m 2 [ n ] = 1 E m = ∑ n e m 2 [ n ] = G 2 ∑ n u m 2 [ n ] = G 2 e_m[n]=Gu_m[n]\quad\quad \sum_nu_m^2[n]=1\quad\quad E_m=\sum_ne_m^2[n]=G^2\sum_nu_m^2[n]=G^2 em[n]=Gum[n]n∑um2[n]=1Em=n∑em2[n]=G2n∑um2[n]=G2

LPC方程求解

• Yule-Walker Equations ∑ j = 1 p a j φ m [ i , j ] = φ m [ i , 0 ] i = 1 , 2 , . . . , p \sum_{j=1}^pa_j\phi_m[i,j]=\phi_m[i,0]\quad\quad i=1,2,…,p ∑j=1pajφm[i,j]=φm[i,0]i=1,2,...,p
  • 本質：矩陣求逆，但是存在有效解法
• 方差矩陣法（Spoken Language Processing , Section 6.3.2.1）
  
• 自相關法
  • 加窗 x m [ n ] = x [ m + n ] w [ n ] x_m[n]=x[m+n]w[n] xm[n]=x[m+n]w[n]
  • 預測誤差 E m = ∑ n = 0 N + p - 1 e m 2 [ n ] E_m=\sum_{n=0}^{N+p-1}e_m^2[n] Em=∑n=0N+p-1em2[n]
  • 自相關
    φ m [ i , j ] = ∑ n = 0 N + p - 1 x m [ n - i ] x m [ n - j ] = ∑ n = 0 N - 1 - ( i - j ) x m [ n ] x m [ n + i - j ] \phi_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] φm[i,j]=n=0∑N+p-1xm[n-i]xm[n-j]=n=0∑N-1-(i-j)xm[n]xm[n+i-j] φ m [ i , j ] = R m [ i - j ] \phi_m[i,j]=R_m[i-j] φm[i,j]=Rm[i-j] R m [ k ] = ∑ n = 0 N - 1 - k x m [ n ] x m [ n + k ] R_m[k]=\sum_{n=0}^{N-1-k}x_m[n]x_m[n+k] Rm[k]=n=0∑N-1-kxm[n]xm[n+k]

語音生成模型

參考線性預測編碼
語音生成模型：

LPC正是基于這個模型的語音生成技術。在該模型中，語音信號是由一個激勵信號 e ( k ) e(k) e(k)經過一個時變的全極點濾波器產生。全極點濾波器的系數取決于所產生的特定聲音的聲道形狀。激勵信號 e k e_{k} ek要么是濁音語音的脈沖序列，要么是無聲聲音的隨機噪聲。生成語音信號 s ( k ) s(k) s(k)可以表示為：
s ( k ) = ∑ p = 1 P a p s ( k - p ) + e ( k ) s(k)=\sum_{p=1}^{P}a_{p}s(k-p)+e(k) s(k)=∑p=1Paps(k-p)+e(k)
其中， P 是濾波器的階數， a p a_{p} ap 是濾波器的系數。LPC就是在已知 s ( k ) s(k) s(k) 的情況下獲取 a p a_{p} ap .

求取 a p a_{p} ap 最常用的一個方法就是最小化真實信號與預測信號之間的均方誤差（Mean Squared Error, MSE）。MSE函數可以表示為
J = E [ e 2 ( k ) ] = E [ ( s ( k ) - ∑ p = 1 P a p s ( k - p ) ) 2 ] 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=E[e2(k)]=E(s(k)-p=1∑Paps(k-p))2
然后，計算 J 關于每個濾波器系數的偏導，并令其值等于0，可得（3）：
J a p = 0 ( 3 ) \frac{\partial J}{\partial a_{p}}=0\quad\quad\quad(3) apJ=0(3)

通過對（3）計算，可以得到（4）：
∑ u = 1 P a u E [ s ( k - p ) s ( k - u ) ] = E [ s ( k ) s ( k - u ) ] , 1 ≤ u ≤ p ( 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\leq u\leq p\quad\quad\quad(4) u=1∑PauE[s(k-p)s(k-u)]=E[s(k)s(k-u)],1≤u≤p(4)

其中， 1 ≤ p ≤ P 1\leq p \leq P 1≤p≤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 ) - x ~ ( n ) = x ( n ) - ∑ i = 1 p a i x ( n - i ) e(n)=x(n)-\widetilde{x}(n)=x(n)-\sum_{i=1}^pa_ix(n-i) e(n)=x(n)-x
(n)=x(n)-∑i=1paix(n-i)
對 a i a_i ai求偏導可得：
∑ n x ( n ) x ( n - j ) = ∑ i = 1 p a i ∑ n x ( n - i ) x ( n - j ) \sum_nx(n)x(n-j)=\sum_{i=1}^pa_i\sum_nx(n-i)x(n-j) n∑x(n)x(n-j)=i=1∑pain∑x(n-i)x(n-j) E = ∑ n [ x ( n ) ] 2 - ∑ i = 1 p a i ∑ n x ( n ) x ( n - i ) E=\sum_n[x(n)]^2-\sum_{i=1}^pa_i\sum_nx(n)x(n-i) E=n∑[x(n)]2-i=1∑pain∑x(n)x(n-i)
寫成自相關形式（Yule-Walker方程）：
R ( j ) = - ∑ i = 1 p a i R ( j - i ) 1 ≤ j ≤ p R(j)=-\sum_{i=1}^pa_iR(j-i)\quad\quad\quad1≤j≤p R(j)=-i=1∑paiR(j-i)1≤j≤p

拆寫加權式子，即為Toeplize矩陣：
使用Durbin算法來求解Toeplize矩陣，即可計算出濾波器系數 a i a_i ai。

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特征提?。?/p>

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('預測語音語譜圖')

效果圖：