基于倒譜法和線性預(yù)測(cè)法估計(jì)基音頻率（MATLAB和Python）

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

一幀信號(hào)的基音頻率估計(jì)

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

一幀信號(hào)估計(jì)出來的基音頻率是296.29Hz，與matlab求出的296.29Hz是一致的。

一段語(yǔ)音引號(hào)的基音頻率估計(jì)

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（下圖）作對(duì)比，畫出的結(jié)果是一致的。

LPC估計(jì)基音頻率

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

通常$s(n)$的加窗自相關(guā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)$,由于自相關(guān)是偶函數(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)}$$
其中左邊為相關(guān)函數(shù)的矩陣，以對(duì)角線為對(duì)稱，其主對(duì)角線以及和主對(duì)角線平行的任何一條斜線上所有的元素相等。這種矩陣稱為托普利茲（Toepliz）矩陣，而這種方程稱為Yule-Walker方程。這種矩陣方程采用遞歸方法求解，基本思想就是遞歸解法分布進(jìn)行，常用的是萊文遜-杜賓（Levinson - Durbin）算法。

該算法的計(jì)算過程為：

• 當(dāng)$i=0$時(shí)，$E_0=r(0)$，$a_0=1$;
• 對(duì)于第$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)}$$

? 對(duì)于$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一幀語(yǔ)音的基音頻率

估計(jì)一幀語(yǔ)音的基音頻率

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

%用LPC法計(jì)算基音頻率 clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); wlen=256; inc=128; % 給出幀長(zhǎng)和幀移 N=length(x1); % 信號(hào)長(zhǎng)度 time=(0:N-1)/Fs; % 計(jì)算出信號(hào)的時(shí)間刻度 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

估計(jì)一段語(yǔ)音的基音頻率

%求出一段語(yǔ)音的基音頻率 clc; close all; clear all; [x1,Fs] = audioread('voice/a9.wav'); subplot(2,1,1); plot(x1); wlen=256; inc=128; % 給出幀長(zhǎng)和幀移 N=length(x1); % 信號(hào)長(zhǎng)度 time=(0:N-1)/Fs; % 計(jì)算出信號(hào)的時(shí)間刻度 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); %自相關(guān)法 % f_b=pitch_admf(framedata,Fs); %平均幅度差f_b=pitch_lpc(framedata,Fs); %LPCx(i)=f_b; end subplot(2,1,2); plot(x);

python實(shí)現(xiàn)

估計(jì)一幀語(yǔ)音的基音頻率

def pitch_lpc(s, p):'''此函數(shù)用LPC法求基音頻率:param s: 一幀數(shù)據(jù):param p: 預(yù)測(cè)階數(shù):return: ar：預(yù)測(cè)系數(shù)'''n = len(s)# 計(jì)算自相關(guān)函數(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開始，遞歸計(jì)算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) # 畫一段語(yǔ)音波形 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坐標(biāo):\n', x) print('基音頻率:\n', f_b) plt.plot(ar) plt.plot(x, y, 'r')

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

估計(jì)一段語(yǔ)音的基音頻率

def pitch_lpc(s, p):'''此函數(shù)用LPC法求基音頻率:param s: 一幀數(shù)據(jù):param p: 預(yù)測(cè)階數(shù):return: ar：預(yù)測(cè)系數(shù)'''n = len(s)# 計(jì)算自相關(guān)函數(shù)Rp = np.zeros(p) # 給一個(gè)預(yù)測(cè)階數(shù)的零矩陣for i in range(p): # 求自相關(guān)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開始，遞歸計(jì)算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] # 求得預(yù)測(cè)系數(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) #找到最大值對(duì)應(yīng)的坐標(biāo)print('Fs=',Fs)print("x:",x)f_b = Fs / (x) # 計(jì)算基頻return f_bwlen = 256 inc = 128 pitch = [] x1, Fs = librosa.load("a9.wav", sr=None) plt.subplot(2,1,1) plt.plot(x1) # 畫一段語(yǔ)音波形 signal = enframe(x1, wlen, inc) # 求一段語(yǔ)音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)一語(yǔ)音的基音頻率相比，兩個(gè)圖基本是一致的。

總結(jié)

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

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