MFCC語音特征值提取算法

部落客簡介

部落客是一名大二學生，主攻人工智能研究。感謝讓我們在CSDN相遇，部落客緻力于在這裡分享關于人工智能，c++，Python，爬蟲等方面知識的分享。 如果有需要的小夥伴可以關注部落客，部落客會繼續更新的，如果有錯誤之處，大家可以指正。

專欄簡介：   本專欄主要研究python在人工智能方面的應用，涉及算法，案例實踐。包括一些常用的資料處理算法，也會介紹很多的Python第三方庫。如果需要，點選這裡   訂閱專欄。

給大家分享一個我很喜歡的一句話：“每天多努力一點，不為别的，隻為日後，能夠多一些選擇，選擇舒心的日子，選擇自己喜歡的人！”

目錄

 背景引入

MFCC語音特征值提取算法簡介

語音信号分幀

計算MFCC系數

對語音信号進行預加重

對語音信号進行短時傅裡葉變換

定義濾波器組

計算MFCC系數

 背景引入

特征值提取，在模式識别領域是很常見的一種算法和手段。特征值看起來好像很陌生，其實在我們日常生活中也很常見。我們使用的身份認證，ID，都可以視為不同系統下的特征值。

MFCC在語音識别領域就是一組特征向量，它通過對語音信号（頻譜包絡與細節）進行編碼運算來得到。MFCC有39個系數，其中包括13個靜态系數，13個一階差分系數，以及13個二階差分系數。差分系數用來描述動态特征，也就是聲學特征在相鄰幀間的變化情況。這些系數都是通過離散餘弦變換（Discrete Cosine Transform，DCT）計算而來。

MFCC語音特征值提取算法簡介

MFCC意為梅爾頻率倒譜系數，顧名思義，MFCC語音特征提取包含兩個關鍵步驟；将語音信号轉化為梅爾頻率，然後進行倒譜分析。梅爾頻譜是一個可用來代表短期音頻的頻譜，梅爾刻度（Mel Scale）則是一種基于人耳對等距的音高變化的感官判斷而确定的非線性頻率刻度。梅爾頻率和正常的頻率f之間的關系：

當梅爾刻度均勻分布，則對應的頻率之間的距離會越來越大。梅爾刻度的濾波器組在低頻部分的分辨率高，跟人耳的聽覺特性比較相符，這也是梅爾刻度的實體意義。在梅爾頻域内，人對音調的感覺度為線性關系，如果兩段語音的梅爾頻率相差兩倍，則人耳聽起來兩者的音調也相差兩倍。

轉化為梅爾頻率時，首先對時域信号進行離散傅裡葉變換，将信号轉換到頻域，然後再利用梅爾刻度的濾波器組對頻域信号進行切分，使每個頻率段對應一個數值。倒譜（Cepstrum）通過對一個時域信号進行傅裡葉變換後取對數，并再次進行反傅裡葉變換（Inverse Fast Fourier Transform，IFFT）得到。倒譜可分為複倒譜（Complex Cepstrum），實倒譜（Real Cepstrum）和功率倒譜（Power Cepstrum）。倒譜分析可用于信号分解，也就是将兩個信号的卷積轉化為兩個信号的相加。

MFCC的實體含義，簡而言之，可了解為語音信号的能量在不同頻率範圍的分布。

人的發聲過程可以看作是肺裡的氣流通過聲帶這個線性系統。如果用e（t）表示輸入聲音的音高，h(t)表示聲帶的響應（也即我們需要擷取的語音特征），那麼聽到的語音信号x(t)即為二者的卷積：

x（t)為時域信号，對其進行離散傅裡葉變換後可得到頻域信号X（K），亦即頻譜：

時域信号的卷積在頻域内則可表示為二者的乘積：

通常，在頻域分析中我們隻關注頻譜的能量而忽略其相位資訊，即：

對頻譜進行對數運算：

然後進行反傅裡葉計算：

c(n)即為倒譜系數，已經和原始的時域信号x(t)不一樣。并且時域信号的卷積關系已經轉化為頻域信号的線性相加關系。

語音信号分幀

語音信号屬于準穩态信号，這也意味着，在一定的時間内，信号會保持穩定。這個時常對于我們人類來說，很短，一般隻有10ms~30ms。在這一區間（即幀）内，可将語音信号看成穩态信号，隻有穩态信号才能進行信号處理。

信号分幀一般會涉及到一個加窗的操作，即将原始信号與一個窗函數相乘。我們用計算機處理信号的時候，一般不會取無限長的信号，而是會取其中間的一段信号，這将會減少工作量，也會加快程式法分析的時間。

無限長的信号被截斷後，其頻譜會發生畸變，進而導緻頻譜能量洩露。 為了減少這種能量洩露，我們可采用不同的截取函數對信号進行截斷。執行截斷操作的函數稱為窗函數，簡稱為窗。常用的窗函數有矩形窗，三角窗，漢明（Hamming）窗及漢甯窗等。

漢甯窗也叫升餘弦窗，是很有用的窗函數。如果測試測試信号有多個頻率分量，頻譜表現非常複雜，測試目的更多在于關注頻率點而非能量大小，則用漢甯窗。漢甯窗主瓣加寬并降低，旁瓣則顯著減小，從減少洩漏的觀點出發，漢甯窗明顯優于矩形窗。但漢甯窗主瓣加寬，相當于分析帶寬加寬，頻率分辨率下降，他與矩形窗相比，洩露以及波動較小，選擇性則相應較高。

漢明窗是用來權重餘弦形成的錐形窗，也稱之為改進的升餘弦窗，隻是權重系數不同，其旁瓣更小，但其旁瓣衰減速度比漢甯窗要慢。漢明窗是以著名的美國數學家理查德·衛斯理·漢明（Richard Wesley Hamming）的名字來命名：

下面的代碼就是用python來生成漢明窗和漢甯窗：

import matplotlib.pyplot as plt
import scipy #信号處理工具包
plt.figure(figsize=(6,2))
plt.plot(scipy.hanning(512),"b--",label="Hanning") #繪制漢甯窗
plt.plot(scipy.hamming(512),"r--",label="Hamming") #繪制漢明窗
plt.title("Demo Hanning & Hamming Window")
plt.legend()
plt.show()
           

除了scipy子產品可以實作漢甯窗和漢明窗，我們也可以用NumPy來實作漢甯窗和漢明窗。示例代碼如下：

import numpy as np
import matplotlib.pyplot as plt
hanWing=np.hanning(512)#定義漢甯窗
hamWin=np.hamming(512) #定義漢明窗
plt.plot(hanWing,'y--',label="Hanning")
plt.plot(hamWin,'b--',label="Hamming")
plt.title("Hamming & Hanning window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.legend()
plt.show()
           

信号加窗，從本質上而言，就是将原始信号與一個窗函數相乘。進行加窗操作之後，我們就可以對信号進行傅裡葉展開。加窗的代價就是，一幀信号的兩端部分将會被消弱。是以在進行信号分幀處理時，幀與幀之間需要有部分重疊。相鄰兩幀重疊後，其起始位置的時間差稱之為幀移，即步長（Stride） 。

以下為簡單的信号加窗操作示意圖：

import numpy as np
import matplotlib.pyplot as plt
import scipy
x=np.linspace(0,10,1000)
originWav=np.sin(x**2) #示例原信号
win=scipy.hamming(1000) #定義一個窗函數，這裡使用的漢明窗
winFrame=originWav*win
#結果可視化
plt.title("Signal Chunk with Hamming Windows")
plt.plot(originWav)
plt.plot(win)
plt.plot(winFrame)
plt.legend()
plt.show()
           

運作程式，其中藍色波形為原信号，橙色波形為窗函數，綠色為加窗操作之後的信号。

假設x為語音信号，w為窗函數，則分幀信号為：

其中，w(n-m)為視窗序列，當n去不同的值時，視窗w(n-m)沿x(m)。是以，w(n-m)是一個“滑動的”視窗。y(n)為短時傅裡葉變換（SIFT）。由于視窗是有限長度的，滿足絕對可和條件，是以這個變幻的前提條件是存在的，這也是信号分幀的理論依據。

以下示例代碼從指定檔案夾讀取一個音頻檔案，然後将該音頻檔案分幀并顯示其中一個分幀信号的波形：

#讀取指定音頻檔案
import matplotlib.pyplot as plt
import numpy as np
import wave #導入波形處理工具包
import os
import soundfile
def audioSignalFrame(signal,nw,inc):
    '''
    signal:原始音頻信号
    nw:每一幀的長度
    inc:相鄰幀的間隔
    '''
    #信号總長度
    signal_length=len(signal)
    #若信号長度小于一個幀的長度，則幀數定義為1
    if signal_length<=nw:
        nf=1
    else:
        nf=int(np.ceil((1.0*signal_length-nw+inc)/inc))
        #所有幀加起來的總的鋪平的長度
    pad_length=int((nf-1)*inc+nw)
    #長度不夠時，使用0填補，類似于FFT中的擴充數組長度
    zeros=np.zeros((pad_length-signal_length,))
    #填補後的信号
    pad_signal=np.concatenate((signal,zeros))
    #相當于對所有幀的時間點進行抽取，得到nf*nw的長度的矩陣
    indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T
    #将indices轉化為矩陣
    indices=np.array(indices,dtype=np.int32)
    #得到幀信号
    frames=pad_signal[indices]
    #窗函數，這裡預設取1
    return frames
def readSignalWave(filename):
    f=wave.open(filename,'rb')
    params=f.getparams()
    nchannels,sampwidth,frammerate,nframes=params[:4]
    #讀取音頻，字元串格式
    strData=f.readframes(nframes)
    #将字元串轉化為int
    waveData=np.fromstring(strData,dtype=np.int16)
    f.close()
    #信号幅值歸一化
    waveData=waveData*1.0/(max(abs(waveData)))
    waveData=np.reshape(waveData,[nframes,nchannels]).T
    return waveData
if __name__=='__main__':
    filepath="./test.wav"
    #dirname=os.listdir(filepath)
    #filename=filepath+dirname[3]
    data=readSignalWave(filepath)
    #初始化每幀長度及幀間隔
    nw=512
    inc=128
    Frame=audioSignalFrame(data[0],nw,inc)
    #顯示原始信号
    plt.plot(data[0])
    plt.title("Original Signal")
    plt.show()

    #顯示第一幀信号
    plt.plot(Frame[0])
    plt.title("First Frame")
    plt.show()
           

上面的代碼中，沒有對信号進行加窗處理，若要執行信号加床操作，隻需将分幀函數稍作修改，

def audioSignalFrame(signal,nw,inc,winfunc):
    '''
    signal:原始音頻信号
    nw:每一幀的長度
    inc:相鄰幀的間隔
    '''
    #信号總長度
    signal_length=len(signal)
    #若信号長度小于一個幀的長度，則幀數定義為1
    if signal_length<=nw:
        nf=1
    else:
        nf=int(np.ceil((1.0*signal_length-nw+inc)/inc))
        #所有幀加起來的總的鋪平的長度
    pad_length=int((nf-1)*inc+nw)
    #長度不夠時，使用0填補，類似于FFT中的擴充數組長度
    zeros=np.zeros((pad_length-signal_length,))
    #填補後的信号
    pad_signal=np.concatenate((signal,zeros))
    #相當于對所有幀的時間點進行抽取，得到nf*nw的長度的矩陣
    indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T
    #将indices轉化為矩陣
    indices=np.array(indices,dtype=np.int32)
    #得到幀信号
    frames=pad_signal[indices]
    #窗函數，這裡預設取1
    win=np.tile(winfunc,(nf,1))

    return frames*win
           

當然，随着函數的改變，主函數中對函數的調用也需要改變，隻需要改變參數即可。除了調用工具包中的漢明窗函數，也可以使用公式來定義。

def hamming(n):

    return 0.54-0.46*cos(2*pi/n*(arange(n)+0.5))

 語音信号在進行分幀之前，一般需要進行一個與加重操作。語音信号的預加重，是為了對語音的高頻部分進行加重，使信号變得平坦，保持在地頻到高頻的整個頻帶中能用同樣的信噪比求頻譜。同時也為了消除發聲過程中聲帶和口唇輻射效應，補償語音信号受到發音系統所抑制的高頻部分，增加語音的高頻分辨率。

我們一般通過一階有限長機關沖激響應（Finite Impulse Response，FIR）高通數字濾波器來實作預加重。FIR濾波器

作為傳遞函數，其中a為預加重系數，0.9<a<1.0。假設t時刻的語音采樣值為x(t)，經過預加重處理後的結果為

.a一般預設取0.95.

信号的預加重處理示例代碼：

def preemphasis(signal,coeff=0.95):
    '''
    signal:要濾波的輸入信号
    coeff:預加重系數。0表示無過濾，預設為0.95
    傳回值：濾波信号
    '''
    return numpy.append(signal[0],signal[1:]-coeff*signal[:-1])
           

計算MFCC系數

由于信号在時域上的變換很難看出特征，是以，我們通常将它轉換為頻域上的能量分布以便于觀察。不同的能量分布，代表不同語音的特征。語音原信号在與窗函數（如漢明窗）相乘後，每幀還必須再經過快速傅裡葉變換以得到頻譜上的能量分布。對語音信号分幀加窗後的各幀的頻譜，然後對頻譜進行取模平方運算後即為語音信号的功率譜。

對信号幅度譜、功率譜以及對數功譜的計算執行個體代碼如下：

import numpy
import logging
def msgspec(frames,NFFT):
    """
    計算幀中每個幀的幅度譜。如果幀為N*D，則輸出N*（NFFT/2+1）
    """
    if numpy.shape(frames)[1]>NFFT:
        logging.warn('frame length (%d)is greater than FFT size(%d),frame will be truncated .Increase NFFT to avoid.',numpy.shape(frames)[1],NFFT)
        complex_spec=numpy.fft.rfft(frames,NFFT)
        return numpy.absolute(complex_spec)
    
def power_spectrum(frames,NFFT):
    return 1.0/NFFT*numpy.square(spectrum_magnitude(frames,NFFT))
def log_power_spectrum(frames,NFFT,norm=1):
    spec_power=power_spectrum(frames,NFFT)
    spec_power[spec_power<1e-30]
    log_spec_power=10*numpy.log10(spec_power)
    if norm:
        return log_spec_power-numpy.max(log_spec_power)
    else:
        return log_spec_power
           

此外，信号的每一幀的音量（即能量），也是語音的特征，而且非常容易計算。是以，通常會再加上一幀的能量，使得每一幀基本的語音特征增加一個次元，包括一個對數能量和倒譜參數。标準的倒譜參數MFCC，隻反映了語音參數的靜态特征，語音參數的動态特征可以用這些靜态特征的差分普來描述。

MFCC的全部組成如下：N維MFCC系數（N/3 MFCC系數+N/3 一階差分系數+N/3二階差分系數）+幀能量。以語音識别中常用的39維MFCC為例，即為：13個靜态系數+13個一階差分系數（Delta系數）+13個二階差分系數（Delta-Delta系數）。其中，差分系數用來描述動态特征，即聲學特征在相鄰幀間的變化情況。

在MFCC計算中還涉及頻率與梅爾刻度之間的轉換，其轉換方式如下：

def hz2mel(hz):

    return 2595*numpy.log10(1+hz/700.0)

同樣，我們也可以推出下列公式：

 Delta系數的計算公式為：

其中，

為Delta系數，從幀t根據靜态系數

 計算而得。N一般取值為2。Delta-Delta（加速度）系數的計算方法相同，但他們是根據Delta而不是靜态系數來進行計算得到的。計算Deltaa系數的示例代碼如下：

def delta(feat,N):
    if N<1:
        raise ValueError('N must be an integer>=1')
    NUMFRAMES=len(feat)
    denominator=2*sum([i**2 for i in range(1,N+1)])
    delta_feat=numpy.pad(feat,((N,N),(0,0)),mode='edge')
    for t in range(NUMFRAMES):
        delta_feat[t]=numpy.dot(numpy.arange(-N,N+1),padded[t:t+2*N+1])/denominator
        return delta_feat
           

 當然除了自己定義函數，也可以直接使用工具包中的API。

對語音信号進行預加重

import numpy as np
import matplotlib.pyplot as plt
from python_speech_features.sigproc import *
from python_speech_features import *
from scipy.fftpack import dct
import scipy.io.wavfile as wav
sample_rate,signal=wav.read('./test.wav')
#保留語音的前3.5秒
signal=signal[0:int(3.5*sample_rate)]
#信号預加重
emphasized_signal=preemphasis(signal,coeff=0.95)
#顯示信号
plt.plot(signal)
plt.title("Original Signal")
plt.plot(emphasized_signal)
plt.title("Preemphasis Signal")
plt.show()
           

上述示例代碼，對信号進行預加重處理的是preemphasis（signal，coeff）函數，除了這個函數，也可以使用以下代碼實作：

 pre_emphasis=0.95

emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1])

源代碼：

import numpy as np
import matplotlib.pyplot as plt
from python_speech_features.sigproc import *
from python_speech_features import *
from scipy.fftpack import dct
import scipy.io.wavfile as wav
sample_rate,signal=wav.read('./test.wav')
pre_emphasis=0.95
emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1])
#保留語音的前3.5秒
#signal=signal[0:int(3.5*sample_rate)]
#信号預加重
#emphasized_signal=preemphasis(signal,coeff=0.95)
#顯示信号
plt.plot(signal)
plt.title("Original Signal")
plt.plot(emphasized_signal)
plt.title("Preemphasis Signal")
plt.show()
           

通過上面的程式可知，兩種函數都可以進行預加重處理，可以自行選擇合适的方法。 

對語音信号進行短時傅裡葉變換

在對語音信号進行處理之前，我們需要對不穩定的語音信号進行短時分幀以擷取傅裡葉變換必需的穩定信号。語音處理範圍内的典型幀大小範圍為20ms~40ms，連續幀之間重疊50%左右。是以一般将幀長度設定為25ms。短時傅裡葉變換（Short-Time Fourier Transform，SIFT）在MFCC計算過程中主要用于短時分幀處理後，通過對信号進行時域到頻域的轉換來擷取語音信号的頻譜。

#對信号進行短時分幀處理
frame_size=0.025 #設定幀長
#計算幀對應采樣數（frame_length)以及步長對應采樣數（frame_step)
frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate
signal_length=len(emphasized_signal) #信号總采樣數
frame_length=int(round(frame_length)) #幀采樣數
frame_step=int(round(frame_step))
#num_frames為總幀數，確定我們至少有一個幀
num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step))
pad_signal_length=num_frames*frame_step+frame_length
z=np.zeros((pad_signal_length-signal_length))
#填充信号以後確定所有的幀的采樣數相等
pad_signal=np.append(emphasized_signal,z)
indices=np.tile(np.arange(0,frame_length),(num_frames,1))+np.tile(np.arange(0,num_frames*frame_step,frame_step),(frame_length,1)).T
frames=pad_signal[indices.astype(np.int32,copy=False)]
           

信号經過短時分幀之後，可通過短時傅裡葉變換得到各種頻譜

NFFT=512
mag_frames=np.absolute(np.fft.rfft(frames,NFFT))
pow_frames=((1.0/NFFT)*((mag_frames)**2))
log_pow_frames=logpowspec(pow_frames,NFFT,norm=1)
#保留語音的前3.5秒
#signal=signal[0:int(3.5*sample_rate)]
#信号預加重
#emphasized_signal=preemphasis(signal,coeff=0.95)
#顯示信号
plt.plot(mag_frames)
plt.title("Mag_Spectrum")
plt.plot(emphasized_signal)
plt.show()
plt.plot(pow_frames)
plt.title("Power_Spectrum")
plt.show()
plt.plot(pow_frames)
plt.title("Log_Power_Spectrum")
plt.show()
           

運作上面的程式，就可以得到處理結果，下面展示原有的所有代碼：

import numpy as np
import matplotlib.pyplot as plt
from python_speech_features.sigproc import *
from python_speech_features import *
from scipy.fftpack import dct
import scipy.io.wavfile as wav
sample_rate,signal=wav.read('./test.wav')
pre_emphasis=0.95
emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1])
#對信号進行短時分幀處理
frame_size=0.025 #設定幀長
frame_stride=0.1
#計算幀對應采樣數（frame_length)以及步長對應采樣數（frame_step)
frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate
signal_length=len(emphasized_signal) #信号總采樣數
frame_length=int(round(frame_length)) #幀采樣數
frame_step=int(round(frame_step))
#num_frames為總幀數，確定我們至少有一個幀
num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step))
pad_signal_length=num_frames*frame_step+frame_length
z=np.zeros((pad_signal_length-signal_length))
#填充信号以後確定所有的幀的采樣數相等
pad_signal=np.append(emphasized_signal,z)
indices=np.tile(np.arange(0,frame_length),(num_frames,1))+np.tile(np.arange(0,num_frames*frame_step,frame_step),(frame_length,1)).T
frames=pad_signal[indices.astype(np.int32,copy=False)]

NFFT=512
mag_frames=np.absolute(np.fft.rfft(frames,NFFT))
pow_frames=((1.0/NFFT)*((mag_frames)**2))
log_pow_frames=logpowspec(pow_frames,NFFT,norm=1)
#保留語音的前3.5秒
#signal=signal[0:int(3.5*sample_rate)]
#信号預加重
#emphasized_signal=preemphasis(signal,coeff=0.95)
#顯示信号
plt.plot(mag_frames)
plt.title("Mag_Spectrum")
plt.plot(emphasized_signal)
plt.show()
plt.plot(pow_frames)
plt.title("Power_Spectrum")
plt.show()
plt.plot(log_pow_frames)
plt.title("Log_Power_Spectrum")
plt.show()
           

                                                        (a)幅度譜

                                                             (b)功率譜

                                                             （c)功率對數譜

音頻檔案使用不同，最終結果也會不同，大家自己使用自己的音頻，注意音頻格式為“.wav”

定義濾波器組

将信号通過一組梅爾刻度的三角形濾波器組，采用的濾波器為三角形濾波器，中心頻率為f(m),m=1,2,3,```````,M,M通常取22~26. 各f(m)之間的間隔随着m值的減少而減少。随着m值的增大而增大。如圖：

 三角形濾波器的頻率響應定義公式：4

對于其他的情況，例如，k<f(m-1)和k>=f(m+1)則為0，當k=f(m)時為1.

定義梅爾刻度的三角形濾波器組的示例代碼為：

low_freq_MEL=0 #将頻率轉換為梅爾刻度
nfilt=40 #窗的數目
#計算m=2595*log10(1+f/700)
high_freq_mel=(2595*np.log10(1+(sample_rate/2)/700))
mel_points=np.linspace(low_freq_MEL,high_freq_mel,nfilt+2) #梅爾刻度的均勻分布
#計算f=700(10**(m/2595)-1)
hz_points=(700*(10**(mel_points/2595)-1))
bin=np.floor((NFFT+1)*hz_points/sample_rate)
fbank=np.zeros((nfilt,int(np.floor(NFFT/2+1))))
#計算三角形濾波器頻率響應
for m in range(1,nfilt+1):
    f_m_minus=int(bin[m-1]) #三角形濾波器左邊頻率f(m-1)
    f_m=int(bin[m]) #三角形濾波器中間頻率fm
    f_m_plus=int(bin[m+1]) #三角形濾波器右邊頻率f(m-1)
    for k in range(f_m_minus,f_m):
        fbank[m-1,k]=(k-bin[m-1])/(bin[m+1]-bin[m])
plt.plot(fbank.T)
plt.show()
           

三角形濾波器有兩個主要功能，其一，對頻譜進行平滑并消除諧波的作用，突顯原先語音的共振峰；其二，用以降低運算量。如圖所示的濾波器組中的每個濾波器在中心頻率處響應為1，并朝着0線性減少，直至達到響應為0的兩個相鄰濾波器的中心頻率。

計算MFCC系數

如果計算出的濾波器組系數高度相關，則在某些機器學習算法中可能會存在問題。我們可用離散餘弦變換對濾波器組系數進行去相關，并産生濾波器組的壓縮表示。濾波器組輸出的對數能量經離散餘弦變換後，即可得到MFCC系數。示例代碼如下：

import numpy as np
import matplotlib.pyplot as plt
from python_speech_features.sigproc import *
from python_speech_features import *
from scipy.fftpack import dct
import scipy.io.wavfile as wav
sample_rate,signal=wav.read('./test.wav')
pre_emphasis=0.95
emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1])
#對信号進行短時分幀處理
frame_size=0.025 #設定幀長
frame_stride=0.1
#計算幀對應采樣數（frame_length)以及步長對應采樣數（frame_step)
frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate
signal_length=len(emphasized_signal) #信号總采樣數
frame_length=int(round(frame_length)) #幀采樣數
frame_step=int(round(frame_step))
#num_frames為總幀數，確定我們至少有一個幀
num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step))
pad_signal_length=num_frames*frame_step+frame_length
z=np.zeros((pad_signal_length-signal_length))
#填充信号以後確定所有的幀的采樣數相等
pad_signal=np.append(emphasized_signal,z)
indices=np.tile(np.arange(0,frame_length),(num_frames,1))+np.tile(np.arange(0,num_frames*frame_step,frame_step),(frame_length,1)).T
frames=pad_signal[indices.astype(np.int32,copy=False)]

NFFT=512
mag_frames=np.absolute(np.fft.rfft(frames,NFFT))
pow_frames=((1.0/NFFT)*((mag_frames)**2))
log_pow_frames=logpowspec(pow_frames,NFFT,norm=1)
#保留語音的前3.5秒
#signal=signal[0:int(3.5*sample_rate)]
#信号預加重
#emphasized_signal=preemphasis(signal,coeff=0.95)
#顯示信号
'''
plt.plot(mag_frames)
plt.title("Mag_Spectrum")
plt.plot(emphasized_signal)
plt.show()
plt.plot(pow_frames)
plt.title("Power_Spectrum")
plt.show()
plt.plot(log_pow_frames)
plt.title("Log_Power_Spectrum")
plt.show()
'''

low_freq_MEL=0 #将頻率轉換為梅爾刻度
nfilt=40 #窗的數目
#計算m=2595*log10(1+f/700)
high_freq_mel=(2595*np.log10(1+(sample_rate/2)/700))
mel_points=np.linspace(low_freq_MEL,high_freq_mel,nfilt+2) #梅爾刻度的均勻分布
#計算f=700(10**(m/2595)-1)
hz_points=(700*(10**(mel_points/2595)-1))
bin=np.floor((NFFT+1)*hz_points/sample_rate)
fbank=np.zeros((nfilt,int(np.floor(NFFT/2+1))))
#計算三角形濾波器頻率響應
for m in range(1,nfilt+1):
    f_m_minus=int(bin[m-1]) #三角形濾波器左邊頻率f(m-1)
    f_m=int(bin[m]) #三角形濾波器中間頻率fm
    f_m_plus=int(bin[m+1]) #三角形濾波器右邊頻率f(m-1)
    for k in range(f_m_minus,f_m):
        fbank[m-1,k]=(k-bin[m-1])/(bin[m+1]-bin[m])
plt.plot(fbank.T)
plt.show()


filter_banks=np.dot(pow_frames,fbank.T)
filter_banks=np.where(filter_banks==0,np.finfo(float).eps,filter_banks)
filter_banks=20*np.log10(filter_banks)
num_ceps=12 #取12個系數
#通過DCT計算MFCC系數
mfcc=dct(filter_banks,type=2,axis=1,norm='ortho')[:,1:(num_ceps+1)]
#對MFCC進行倒譜提升可以改善噪聲信号中的語音識别
(nframes,ncoeff)=mfcc.shape
n=np.arange(ncoeff)
cep_lifter=22 #倒譜濾波系數，定義倒譜所用到的濾波器組内濾波器個數
lift=1+(cep_lifter/2)*np.sin(np.pi*n/cep_lifter)
mfcc*=lift
mfcc-=(np.mean(mfcc,axis=0)+1e-8)
plt.imshow(np.flipud(mfcc.T),cmap=plt.cm.jet,aspect=0.2,extent=[0,mfcc.shape[0],0,mfcc.shape[1]]) #繪制MFCC熱力圖
plt.show()
           

對MFCC進行如下的歸一化操作，運作操作，其相應的熱力圖如下：

filter_banks-=(np.mean(filter_banks,axis=0)+1e-8)
plt.imshow(np.flipud(filter_banks.T),cmap=plt.cm.jet,aspect=0.2,extent=[0,filter_banks.shape[1],0,filter_banks.shape[0]])
plt.show()
           

                                 歸一化的MFCC熱力圖

好了，本篇文章介紹就到此結束了，拜了個拜！