基于倒譜法和線性預測法估計基音頻率（MATLAB和Python）

倒譜法基音檢測在python中實現(xiàn)

一幀信號的基音頻率估計

wlen = 256 inc = 128 pitch = [] x1, Fs = librosa.load("a9.wav",sr=None) plt.subplot(2,1,1) # plt.plot(x1) # 畫一段語音波形 signal = enframe(x1, wlen, inc) # 取一幀 framedata = signal[15] plt.plot(framedata) f_b = pitch_cep(framedata, Fs) plt.show()

一幀信號估計出來的基音頻率是296.29Hz，與matlab求出的296.29Hz是一致的。

一段語音引號的基音頻率估計

def pitch_cep(x, Fs):'''用倒譜法求基音頻率:param x: 分幀后的數(shù)據(jù):param Fs: 采樣率:return: f_b該幀的基音頻率'''x2 = fftpack.fft(x)amp = 20 * np.log10(abs(x2) + 0.0000001)x3 = abs(fftpack.fft(amp))x3[0:27] = 0x3[115:256] = 0y = max(x3)x = x3.tolist().index(y)f_b = Fs / (x - 1.0)# print('基音頻率:', f_b, 'Hz')return f_b wlen = 256 inc = 128 pitch = [] x1, Fs = librosa.load("a9.wav",sr=None) plt.subplot(2,1,1) plt.plot(x1) signal = enframe(x1, wlen, inc) # 取一幀 for i in signal:framedata = if_b = pitch_cep(framedata, Fs)pitch.append(f_b) plt.subplot(2,1,2) plt.plot(pitch) plt.show()

和MATLAB（下圖）作對比，畫出的結果是一致的。

LPC估計基音頻率

線性預測分析是通過矩陣的特殊性質(zhì)來解包含p個未知數(shù)的p個線性方程。自相關解法是的原理是在整個時間范圍內(nèi)使誤差最小，即設$s(n)$在$0?n?N?1$以外的值都是零，等同于假設了$s(n)$經(jīng)過了有限長度的矩形窗、海寧窗或者漢寧窗，就可以用p個方程來解有p個未知數(shù)的方程組了。

通常$s(n)$的加窗自相關函數(shù)定義為：
$$r(j)=\sum _{n=0}^{N?1}s(n)s(n?j),(1?j?p)\text{\hspace{1em}(1)}$$
由于$?(j,i)$等效于$r(j?i)$,由于自相關是偶函數(shù)，所以有：
$$?(j,i)=r(|j?i|)\text{\hspace{1em}(2)}$$
因此式
$$\phi (j,0)=\sum _{i=1}^p{a}_i\phi (j,i),(1?j?p)\text{\hspace{1em}(3)}$$
可以表示為：
$$r(j)=\sum _{i=1}^p{a}_{i}r(|j?i|),(1?j?p)\text{\hspace{1em}(4)}$$
最小均方誤差改寫為：
$$E=r(0)?\sum _{i=1}^p{a}_{i}r(i)\text{\hspace{1em}(5)}$$
展開式(5)可得到方程組：
$$?\begin{array}{ccccc}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{array}??\begin{array}{c}{a}_1\\ a_2\\ a_3\\ ...\\ a_p\end{array}?=?\begin{array}{c}r(1)\\ r(2)\\ r(3)\\ ...\\ r(p)\end{array}?\text{\hspace{1em}(6)}$$
其中左邊為相關函數(shù)的矩陣，以對角線為對稱，其主對角線以及和主對角線平行的任何一條斜線上所有的元素相等。這種矩陣稱為托普利茲（Toepliz）矩陣，而這種方程稱為Yule-Walker方程。這種矩陣方程采用遞歸方法求解，基本思想就是遞歸解法分布進行，常用的是萊文遜-杜賓（Levinson - Durbin）算法。

該算法的計算過程為：

• 當$i=0$時，$E_0=r(0)$，$a_0=1$;
• 對于第$i$次遞歸$(i=1,2,\dots ,p)$

$$k_i=\frac{1}{E_{i?1}}[r(i)?\sum _{j=1}^{i?1}{a}_j^{i?1}r(j?i)]\text{\hspace{1em}(7)}$$

$$i{a}_i^{(i)}=k_i\text{\hspace{1em}(8)}$$

? 對于$j=1到i?1$
$$a_j^{(i)}=a_j^{(i?1)}?a_{i?j}^{(i?1)}\text{\hspace{1em}(9)}$$

$$E_i=(i?k_i^2)E_{i?1}\text{\hspace{1em}(9)}$$

• 增益$G$為

$$G=\sqrt{E_p}$$

• 遞歸得到

$$E_p=r(0)\prod _{i=1}^p(1?k_i^2)$$

Matlab一幀語音的基音頻率

估計一幀語音的基音頻率

matlab中有l(wèi)pc函數(shù)直接可以調(diào)用，因此在python中編寫算法去計算，利用matlab中的lpc函數(shù)來驗證算法的正確性。

%用LPC法計算基音頻率 clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); wlen=256; inc=128; % 給出幀長和幀移 N=length(x1); % 信號長度 time=(0:N-1)/Fs; % 計算出信號的時間刻度 signal=enframe(x1,wlen,inc)'; % 分幀 framedata = signal(:,15); subplot(2,1,1) plot(framedata); %LPC [x3,r]=lpc(framedata,256); x3=abs(x3); x3(1:27) = 0;%在話音基頻范圍外的都取零 x3(115:256) = 0; [M,idx] = max(x3); subplot(2,1,2); plot(x3); f_b=Fs/(idx-1);

所求出的基音頻率是285.71Hz

估計一段語音的基音頻率

%求出一段語音的基音頻率 clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); subplot(2,1,1); plot(x1); wlen=256; inc=128; % 給出幀長和幀移 N=length(x1); % 信號長度 time=(0:N-1)/Fs; % 計算出信號的時間刻度 signal=enframe(x1,wlen,inc)'; % 分幀 [n,m]=size(signal); for i=1:mframedata = signal(:,i); % f_b=pitch_cep(framedata,Fs); %倒譜法 % f_b=pitch_cor(framedata,Fs); %自相關法 % f_b=pitch_admf(framedata,Fs); %平均幅度差f_b=pitch_lpc(framedata,Fs); %LPCx(i)=f_b; end subplot(2,1,2); plot(x);

python實現(xiàn)

估計一幀語音的基音頻率

def pitch_lpc(s, p):'''此函數(shù)用LPC法求基音頻率:param s: 一幀數(shù)據(jù):param p: 預測階數(shù):return: ar：預測系數(shù)'''n = len(s)# 計算自相關函數(shù)Rp = np.zeros(p) #創(chuàng)建0矩陣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, Gwlen = 256 inc = 128 pitch = [] x1, Fs = librosa.load("a9.wav",sr=None) # plt.plot(x1) # 畫一段語音波形 signal = enframe(x1, wlen, inc) # 取一幀 framedata = signal[15] ar, G = pitch_lpc(framedata, p=256) ar = abs(ar)ar[0:27] = 0 ar[115:257] = 0 y = max(ar) x = ar.tolist().index(y) f_b = Fs / x print('x坐標:\n', x) print('基音頻率:\n', f_b) plt.plot(ar) plt.plot(x, y, 'r')

所求出的基音頻率是296.29Hz，與matlab求出的有10Hz左右的差距

估計一段語音的基音頻率

def pitch_lpc(s, p):'''此函數(shù)用LPC法求基音頻率:param s: 一幀數(shù)據(jù):param p: 預測階數(shù):return: ar：預測系數(shù)'''n = len(s)# 計算自相關函數(shù)Rp = np.zeros(p) # 給一個預測階數(shù)的零矩陣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] # 求得預測系數(shù)G = np.sqrt(Ep[p - 1]) # 得到遞歸增益Gar = abs(ar)ar[0:27] = 0 #將話音范圍外置零ar[115:257] = 0y = max(ar) #找最大值x = ar.tolist().index(y) #找到最大值對應的坐標print('Fs=',Fs)print("x:",x)f_b = Fs / (x) # 計算基頻return f_bwlen = 256 inc = 128 pitch = [] x1, Fs = librosa.load("a9.wav", sr=None) plt.subplot(2,1,1) plt.plot(x1) # 畫一段語音波形 signal = enframe(x1, wlen, inc) # 求一段語音pitch用 for i in signal:framedata = if_b = pitch_lpc(framedata, 256)pitch.append(f_b) plt.subplot(2,1,2) plt.plot(pitch) plt.show()

與MATLAB畫出的統(tǒng)一語音的基音頻率相比，兩個圖基本是一致的。

總結

以上是生活随笔為你收集整理的基于倒谱法和线性预测法估计基音频率（MATLAB和Python）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。