快捷導(dǎo)航

基于MATLAB和Python實(shí)現(xiàn)MFCC特征參數(shù)提取

更新時(shí)間：2019年08月13日 10:47:08 作者：趙至柔

這篇文章主要介紹了基于MATLAB和Python實(shí)現(xiàn)MFCC特征參數(shù)提取，文中通過示例代碼介紹的非常詳細(xì)，對(duì)大家的學(xué)習(xí)或者工作具有一定的參考學(xué)習(xí)價(jià)值，需要的朋友們下面隨著小編來一起學(xué)習(xí)學(xué)習(xí)吧

1、MFCC概述

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

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

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

2、 MFCC特征參數(shù)提取過程詳解

(1)預(yù)處理

預(yù)處理包括預(yù)加重、分幀、加窗函數(shù)。

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

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

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

后，W(n)形式如下：

(2)FFT

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

(3)譜線能量

對(duì)語音信號(hào)的頻譜取模平方得到語音信號(hào)的譜線能量。

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

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

三角濾波器的頻率響應(yīng)定義為:

對(duì)頻譜進(jìn)行平滑化，并消除諧波的作用，突顯原先語音的共振峰。（因此一段語音的音調(diào)或音高，是不會(huì)呈現(xiàn)在MFCC 參數(shù)內(nèi)，換句話說，以MFCC 為特征的語音辨識(shí)系統(tǒng)，并不會(huì)受到輸入語音的音調(diào)不同而有所影響）此外，還可以降低運(yùn)算量。

計(jì)算每個(gè)濾波器組輸出的對(duì)數(shù)能量為：

(5)計(jì)算DCT倒譜

經(jīng)離散余弦變換（DCT）得到MFCC系數(shù) :

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

3、MATLAB實(shí)現(xiàn)方法

注：在提取MFCC參數(shù)之前需要加載并使用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濾波器的階數(shù)為24，fft變換的長(zhǎng)度為256，采樣頻率為8000Hz 
%歸一化mel濾波器組系數(shù) 
bank=full(bank); 
bank=bank/max(bank(:)); 
% DCT系數(shù),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);
%語音信號(hào)分幀 
xx=enframe(x,256,80);%對(duì)x 256點(diǎn)分為一幀  
%計(jì)算每幀的MFCC參數(shù) 
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');

結(jié)果：

4、Python實(shí)現(xiàn)方法

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

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

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)
#預(yù)加重
#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)  #轉(zhuǎn)置的原因是分幀函數(shù)enframe的輸出矩陣是幀數(shù)*幀長(zhǎng)
frameNum =X.shape[0] #返回矩陣行數(shù)18，獲取幀數(shù)
#print(frameNum)
for i in range(frameNum):
  y=X[i,:]
  #fft
  yf = np.abs(np.fft.fft(y)) 
  #print(yf.shape)
  #譜線能量
  yf = yf**2
  #梅爾濾波器系數(shù)
  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) # 數(shù)值穩(wěn)定性
  filter_banks = 10 * np.log10(filter_banks) # dB 
  filter_banks -= (np.mean(filter_banks, axis=0) + 1e-8)
  #print(filter_banks)
  #DCT系數(shù)
  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()

結(jié)果: