如何使用MATLAB和Python實現MFCC特征參數提取

這篇文章主要介紹如何使用MATLAB和Python實現MFCC特征參數提取，文中介紹的非常詳細，具有一定的參考價值，感興趣的小伙伴們一定要看完！

1、MFCC概述

在語音識別（Speech Recognition）和話者識別（Speaker Recognition）方面，最常用到的語音特征就是梅爾倒譜系數（Mel-scale FrequencyCepstral Coefficients，簡稱MFCC）。根據人耳聽覺機理的研究發現，人耳對不同頻率的聲波有不同的聽覺敏感度。從200Hz到5000Hz的語音信號對語音的清晰度影響較大。兩個響度不等的聲音作用于人耳時，則響度較高的頻率成分的存在會影響到對響度較低的頻率成分的感受，使其變得不易察覺，這種現象稱為掩蔽效應。由于頻率較低的聲音在內耳蝸基底膜上行波傳遞的距離大于頻率較高的聲音，故一般來說，低音容易掩蔽高音，而高音掩蔽低音較困難。在低頻處的聲音掩蔽的臨界帶寬較高頻要小。所以，人們從低頻到高頻這一段頻帶內按臨界帶寬的大小由密到疏安排一組帶通濾波器，對輸入信號進行濾波。將每個帶通濾波器輸出的信號能量作為信號的基本特征，對此特征經過進一步處理后就可以作為語音的輸入特征。由于這種特征不依賴于信號的性質，對輸入信號不做任何的假設和限制，又利用了聽覺模型的研究成果。因此，這種參數比基于聲道模型的LPCC相比具有更好的魯邦性，更符合人耳的聽覺特性，而且當信噪比降低時仍然具有較好的識別性能。

梅爾倒譜系數是在Mel標度頻率域提取出來的倒譜參數，Mel標度描述了人耳頻率的非線性特性，它與頻率的關系可用下式近似表示：

式中f為頻率，單位為Hz。下圖為Mel頻率與線性頻率的關系：

2、 MFCC特征參數提取過程詳解

(1)預處理

預處理包括預加重、分幀、加窗函數。

預加重：預加重的目的是提升高頻部分，使信號的頻譜變得平坦，保持在低頻到高頻的整個頻帶中，能用同樣的信噪比求頻譜。同時，也是為了消除發生過程中聲帶和嘴唇的效應，來補償語音信號受到發音系統所抑制的高頻部分，也為了突出高頻的共振峰。預加重處理其實是將語音信號通過一個高通濾波器：

分幀：先將N個采樣點集合成一個觀測單位，稱為幀。通常情況下N的值為256或512，涵蓋的時間約為20~30ms左右。為了避免相鄰兩幀的變化過大，因此會讓兩相鄰幀之間有一段重疊區域，此重疊區域包含了M個取樣點，通常M的值約為N的1/2或1/3。通常語音識別所采用語音信號的采樣頻率為8KHz或16KHz，以8KHz來說，若幀長度為256個采樣點，則對應的時間長度是256/8000×1000=32ms。

加窗：將每一幀乘以漢明窗，以增加幀左端和右端的連續性。假設分幀后的信號為S(n), n=0,1…,N-1, N為幀的大小，那么乘上漢明窗

后 ，W(n)形式如下：

(2)FFT

由于信號在時域上的變換通常很難看出信號的特性，所以通常將它轉換為頻域上的能量分布來觀察，不同的能量分布，就能代表不同語音的特性。所以在乘上漢明窗后，每幀還必須再經過快速傅里葉變換以得到在頻譜上的能量分布。對分幀加窗后的各幀信號進行快速傅里葉變換得到各幀的頻譜。

(3)譜線能量

對語音信號的頻譜取模平方得到語音信號的譜線能量。

(4)計算通過Mel濾波器的能量

將能量譜通過一組Mel尺度的三角形濾波器組，定義一個有M個濾波器的濾波器組（濾波器的個數和臨界帶的個數相近），采用的濾波器為三角濾波器，中心頻率為f(m) 。M通常取22-26。各f(m)之間的間隔隨著m值的減小而縮小，隨著m值的增大而增寬，如圖所示：

三角濾波器的頻率響應定義為:

對頻譜進行平滑化，并消除諧波的作用，突顯原先語音的共振峰。（因此一段語音的音調或音高，是不會呈現在MFCC 參數內，換句話說，以MFCC 為特征的語音辨識系統，并不會受到輸入語音的音調不同而有所影響）此外，還可以降低運算量。

計算每個濾波器組輸出的對數能量為 ：

(5)計算DCT倒譜

經離散余弦變換（DCT）得到MFCC系數 :

將上述的對數能量帶入離散余弦變換，求出L階的Mel參數。L階指MFCC系數階數，通常取12-16。這里M是三角濾波器個數。

3、MATLAB實現方法

注：在提取MFCC參數之前需要加載并使用VOICEBOX工具包

Df=5;
fs=8000; 
N=fs/Df;
t=0:1./fs:(N-1)./fs;
x=sin(2*pi*200*t);  
bank=melbankm(24,256,8000,0,0.5,'t');%Mel濾波器的階數為24，fft變換的長度為256，采樣頻率為8000Hz 
%歸一化mel濾波器組系數 
bank=full(bank); 
bank=bank/max(bank(:)); 
% DCT系數,12*p
for k=1:12      
n=0:23; 
dctcoef(k,:)=cos((2*n+1)*k*pi/(2*24)); 
end 
%歸一化倒譜提升窗口 
w=1+6*sin(pi*[1:12]./12); 
%w=w/max(w);
%語音信號分幀 
xx=enframe(x,256,80);%對x 256點分為一幀  
%計算每幀的MFCC參數 
for i=1:size(xx,1) 
y=xx(i,:); 
s=y'.*hamming(256); 
t=abs(fft(s));%fft快速傅立葉變換 
t=t.^2; 
c1=dctcoef*log(bank*t(1:129)); 
c2=c1.*w'; 
end  
plot(c2);title('MFCC');

結果：

4、Python實現方法

import numpy as np 
from scipy import signal
from scipy.fftpack import dct
import pylab as plt

def enframe(wave_data, nw, inc, winfunc):
  '''將音頻信號轉化為幀。
  參數含義：
  wave_data:原始音頻型號
  nw:每一幀的長度(這里指采樣點的長度，即采樣頻率乘以時間間隔)
  inc:相鄰幀的間隔（同上定義）
  '''
  wlen=len(wave_data) #信號總長度
  if wlen<=nw: #若信號長度小于一個幀的長度，則幀數定義為1
    nf=1
  else: #否則，計算幀的總長度
    nf=int(np.ceil((1.0*wlen-nw+inc)/inc))
  pad_length=int((nf-1)*inc+nw) #所有幀加起來總的鋪平后的長度
  zeros=np.zeros((pad_length-wlen,)) #不夠的長度使用0填補，類似于FFT中的擴充數組操作
  pad_signal=np.concatenate((wave_data,zeros)) #填補后的信號記為pad_signal
  indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T #相當于對所有幀的時間點進行抽取，得到nf*nw長度的矩陣
  indices=np.array(indices,dtype=np.int32) #將indices轉化為矩陣
  frames=pad_signal[indices] #得到幀信號
  win=np.tile(winfunc,(nf,1)) #window窗函數，這里默認取1
  return frames*win  #返回幀信號矩陣

Df=5
fs=8000
N=fs/Df
t = np.arange(0,(N-1)/fs,1/fs)   
wave_data=np.sin(2*np.pi*200*t)
#預加重
#b,a = signal.butter(1,1-0.97,'high')
#emphasized_signal = signal.filtfilt(b,a,wave_data)
#歸一化倒譜提升窗口
lifts=[]
for n in range(1,13):
  lift =1 + 6 * np.sin(np.pi * n / 12)
  lifts.append(lift)
#print(lifts)  

#分幀、加窗 
winfunc = signal.hamming(256) 
X=enframe(wave_data, 256, 80, winfunc)  #轉置的原因是分幀函數enframe的輸出矩陣是幀數*幀長
frameNum =X.shape[0] #返回矩陣行數18，獲取幀數
#print(frameNum)
for i in range(frameNum):
  y=X[i,:]
  #fft
  yf = np.abs(np.fft.fft(y)) 
  #print(yf.shape)
  #譜線能量
  yf = yf**2
  #梅爾濾波器系數
  nfilt = 24
  low_freq_mel = 0
  NFFT=256
  high_freq_mel = (2595 * np.log10(1 + (fs / 2) / 700)) # 把 Hz 變成 Mel
  mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2) # 將梅爾刻度等間隔
  hz_points = (700 * (10**(mel_points / 2595) - 1)) # 把 Mel 變成 Hz
  bin = np.floor((NFFT + 1) * hz_points / fs)
  fbank = np.zeros((nfilt, int(np.floor(NFFT / 2 + 1))))
  for m in range(1, nfilt + 1):
    f_m_minus = int(bin[m - 1])  # left
    f_m = int(bin[m])       # center
    f_m_plus = int(bin[m + 1])  # right
    for k in range(f_m_minus, f_m):
      fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1])
    for k in range(f_m, f_m_plus):
      fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m])
  filter_banks = np.dot(yf[0:129], fbank.T)
  filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks) # 數值穩定性
  filter_banks = 10 * np.log10(filter_banks) # dB 
  filter_banks -= (np.mean(filter_banks, axis=0) + 1e-8)
  #print(filter_banks)
  #DCT系數
  num_ceps = 12
  c2 = dct(filter_banks, type=2, axis=-1, norm='ortho')[ 1 : (num_ceps + 1)] # Keep 2-13
  c2 *= lifts
print(c2)
plt.plot(c2)
plt.show()

結果:

由MATLAB和Python繪制出來的波形可以看出二者計算出來的MFCC倒譜系數的基本走勢相同。

以上是“如何使用MATLAB和Python實現MFCC特征參數提取”這篇文章的所有內容，感謝各位的閱讀！希望分享的內容對大家有幫助，更多相關知識，歡迎關注億速云行業資訊頻道！