語音信號的梅爾頻率倒譜系數(MFCC)的原理講解及python實現

梅爾倒譜系數（MFCC）

　　梅爾倒譜系數（Mel-scale FrequencyCepstral Coefficients，簡稱MFCC）。依據人的聽覺實驗結果來分析語音的頻譜，

MFCC分析依據的聽覺機理有兩個

第一Mel scale：人耳感知的聲音頻率和聲音的實際頻率並非線性的，有下面公式

$$f_{mel}=2595*\log _{10}(1+\frac{f}{700})$$

$$f = 700 (10^{f_{mel}/2595} - 1)$$

式中$f_{mel}$是以梅爾(Mel)爲單位的感知頻域（簡稱梅爾頻域），$f$是以$Hz$爲單位的實際語音頻率。$f_{mel}$與$f$的關係曲線以下圖所示，若能將語音信號的頻域變換爲感知頻域中，能更好的模擬聽覺過程的處理。

第二臨界帶（Critical Band）：把進入人耳的聲音頻率用臨界帶進行劃分，將語音在頻域上就被劃分紅一系列的頻率羣，組成了濾波器組，即Mel濾波器組。

　　研究代表，人耳對不一樣頻率的聲波有不一樣的聽覺敏感度。從200Hz到5000Hz的語音信號對語音的清晰度影響較大。兩個響度不等的聲音做用於人耳時，則響度較高的頻率成分的存在會影響到對響度較低的頻率成分的感覺，使其變得不易察覺，這種現象稱爲掩蔽效應。

　　因爲頻率較低的聲音（低音）在內耳蝸基底膜上行波傳遞距離大於頻率較高的聲音（高音），所以低音容易掩蔽高音。低音掩蔽的臨界帶寬較高頻要小。因此，人們從低頻到高頻這一段頻帶內按臨界帶寬的大小由密到疏安排一組帶通濾波器，對輸入信號進行濾波。將每一個帶通濾波器輸出的信號能量做爲信號的基本特徵，對此特徵通過進一步處理後就能夠做爲語音的輸入特徵。因爲這種特徵不依賴於信號的性質，對輸入信號不作任何的假設和限制，又利用了聽覺模型的研究成果。所以，這種參數比基於聲道模型的LPCC相比具備更好的魯棒性，更符合人耳的聽覺特性，並且當信噪比下降時仍然具備較好的識別性能。

求MFCC的步驟

1. 將信號幀化爲短幀
2. 對於每一個幀，計算每幀語音的功率譜（週期圖估計）
   
3. 將mel濾波器組應用於功率譜，求濾波器組的能量，將每一個濾波器中的能量相加
4. 取全部濾波器組能量的對數
5. DCT變換

MFCC的提取過程

預處理

預處理包括預加劇、分幀、加窗函數。假設咱們的語音信號採樣頻率爲8000Hz，語音數據在這裏獲取

import numpy
import scipy.io.wavfile
from scipy.fftpack import dct

sample_rate, signal = scipy.io.wavfile.read('OSR_us_000_0010_8k.wav') 
signal = signal[0:int(3.5 * sample_rate)]  # 咱們只取前3.5s

 時域中的語音信號

預加劇

預加劇濾波器在如下幾個方面頗有用：

1. 平衡頻譜，由於高頻一般與較低頻率相比具備較小的幅度，
2. 避免在傅里葉變換操做期間的數值問題和
3. 改善信號 - 噪聲比（SNR）
4. 消除發聲過程當中聲帶和嘴脣的效應，來補償語音信號受到發音系統所抑制的高頻部分，也爲了突出高頻的共振峯。

預加劇處理實際上是將語音信號經過一個高通濾波器：

$$y(t) = x(t) - \alpha x(t-1)$$

其中濾波器係數（α）的一般爲0.95或0.97，這裏取pre_emphasis =0.97：

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

 預加劇後的時域信號

分幀

　　在預加劇以後，咱們須要將信號分紅短時幀。所以在大多數狀況下，語音信號是非平穩的，對整個信號進行傅里葉變換是沒有意義的，由於咱們會隨着時間的推移丟失信號的頻率輪廓。語音信號是短時平穩信號。所以咱們在短時幀上進行傅里葉變換，經過鏈接相鄰幀來得到信號頻率輪廓的良好近似。

　　將信號幀化爲20-40 ms幀。25毫秒是標準的frame_size = 0.025。這意味着8kHz信號的幀長度爲0.025 * 8000 = 200個採樣點。幀移一般爲10m sframe_stride = 0.01（80個採樣點），爲了不相鄰兩幀的變化過大，所以會讓兩相鄰幀之間有一段重疊區域，此重疊區域包含了120個取樣點，一般約爲每幀語音的1/2或1/3。第一個語音幀0開始，下一個語音幀從160開始，直到到達語音文件的末尾。若是語音文件沒有劃分爲偶數個幀，則用零填充它以使其完成。

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))
# 確保咱們至少有1幀
num_frames = int(numpy.ceil(float(numpy.abs(signal_length - frame_length)) / frame_step))  

pad_signal_length = num_frames * frame_step + frame_length
z = numpy.zeros((pad_signal_length - signal_length))
# 填充信號，確保全部幀的採樣數相等，而不從原始信號中截斷任何採樣
pad_signal = numpy.append(emphasized_signal, z) 

indices = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1)) + numpy.tile(numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T
frames = pad_signal[indices.astype(numpy.int32, copy=False)]

加窗

　　將每一幀乘以漢明窗，以增長幀左端和右端的連續性。抵消fft假設的數據是無限的，並減小頻譜泄漏。

假設分幀後的信號爲S(n), n=0,1…,N-1, N爲窗口長度，那麼乘上漢明窗

$${S}'(n)=S(n)\times W(n)$$

後 ，漢明窗W(n)形式以下：

 $$W\left( n,a \right)=\left( 1-a \right)-a\times \cos \left( \frac{2\pi n}{N-1} \right),0\le n\le N-1$$

frames *= numpy.hamming(frame_length)
# frames *= 0.54 - 0.46 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1))  # 內部實現

FFT

　　因爲信號在時域上的變換一般很難看出信號的特性，因此一般將它轉換爲頻域上的能量分佈來觀察，不一樣的能量分佈，就能表明不一樣語音的特性。對分幀加窗後的各幀信號進行作一個N點FFT來計算頻譜，也稱爲短時傅立葉變換（STFT），其中N一般爲256或512，NFFT＝512；

$$S_i(k)=\sum_{n=1}^{N}s_i(n)e^{-j2\pi kn/N} 1\le k \le K$$

mag_frames = numpy.absolute(numpy.fft.rfft(frames, NFFT))   # fft的幅度

能量譜

　　計算功率譜（週期圖），對語音信號的頻譜取模平方（取對數或者去平方，由於頻率不可能爲負，負值要捨去）獲得語音信號的譜線能量。

$$P = \frac{|FFT(x_i)|^2}{N}$$

其中，$X_i$是信號X的第$i$幀，這能夠用如下幾行來實現：

pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2))  # 功率譜

使用Mel刻度濾波器組過濾

　　計算經過Mel濾波器，將功率譜經過一組Mel刻度(這裏取40個)的三角濾波器組提取帶寬。

　　這個Mel濾波器組就像人類的聽覺感知系統（耳朵），人耳只關注某些特定的頻率份量（人的聽覺對頻率是有選擇性的）。它對不一樣頻率信號的靈敏度是不一樣的，換言之，它只讓某些頻率的信號經過，而壓根就直接無視它不想感知的某些頻率信號。可是這些濾波器在頻率座標軸上卻不是統一分佈的，在低頻區域有不少的濾波器，他們分佈比較密集，但在高頻區域，濾波器的數目就變得比較少，分佈很稀疏。

　　咱們能夠用下面的公式，在語音頻率和Mel頻率間轉換

$$f_{mel}=2595*\log _{10}(1+\frac{f}{700})$$

$$f = 700 (10^{f_{mel}/2595} - 1)$$

　　定義一個有M個三角濾波器的濾波器組（濾波器的個數和臨界帶的個數相近），中心頻率爲f(m) ，中心頻率響應爲1，而且朝向0線性減少，直到它達到響應爲0的兩個相鄰濾波器的中心頻率，以下圖所示。M一般取22-40，26是標準，本文取nfilt = 40。各f(m)之間的間隔隨着m值的增大而增寬，如圖所示：

在Mel-Scale上的濾波器組

這能夠經過如下等式建模，三角濾波器的頻率響應定義爲: 

$$H_m(k) =
\begin{cases}
\hfill 0 \hfill & k < f(m - 1) \\
\\
\hfill \dfrac{k - f(m - 1)}{f(m) - f(m - 1)} \hfill & f(m - 1) \leq k < f(m) \\
\\
\hfill 1 \hfill & k = f(m) \\
\\
\hfill \dfrac{f(m + 1) - k}{f(m + 1) - f(m)} \hfill & f(m) < k \leq f(m + 1) \\
\\
\hfill 0 \hfill & k > f(m + 1) \\
\end{cases}$$

　　對於FFT獲得的幅度譜，分別跟每個濾波器進行頻率相乘累加，獲得的值即爲該幀數據在在該濾波器對應頻段的能量值。若是濾波器的個數爲22，那麼此時應該獲得22個能量值

low_freq_mel = 0
high_freq_mel = (2595 * numpy.log10(1 + (sample_rate / 2) / 700))  # 將Hz轉換爲Mel
mel_points = numpy.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # 使得Mel scale間距相等
hz_points = (700 * (10**(mel_points / 2595) - 1))  # 將Mel轉換爲Hz
bin = numpy.floor((NFFT + 1) * hz_points / sample_rate)

fbank = numpy.zeros((nfilt, int(numpy.floor(NFFT / 2 + 1))))
for m in range(1, nfilt + 1):
    f_m_minus = int(bin[m - 1])   # 左
    f_m = int(bin[m])             # 中
    f_m_plus = int(bin[m + 1])    # 右

    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 = numpy.dot(pow_frames, fbank.T)
filter_banks = numpy.where(filter_banks == 0, numpy.finfo(float).eps, filter_banks)  # 數值穩定性
filter_banks = 20 * numpy.log10(filter_banks)  # dB

信號的功率譜通過濾波器組後，獲得的譜圖爲：

 信號的頻譜圖

若是通過Mel scale濾波器組獲得的頻譜是所需的特徵，那麼咱們能夠跳過均值歸一化。

能量值取對數

因爲人耳對聲音的感知並非線性的，用log這種非線性關係更好描述。取完log之後才能夠進行倒譜分析。

 $$s\left( m \right)=\ln \left( {{\sum\limits_{K=0}^{N-1}{\left| {{X}_{a}}\left( k \right) \right|}}^{2}}{{H}_{m}}\left( k \right) \right),0\le m\le M$$

離散餘弦變換（DCT）

　　前一步驟中計算的濾波器組係數是高度相關的，這在一些機器學習算法中多是有問題的。所以，咱們能夠應用離散餘弦變換（DCT）去相關濾波器組係數併產生濾波器組的壓縮表示。一般，對於自動語音識別（ASR），保留所獲得的倒頻譜系數2-13，其他部分被丟棄; num_ceps = 12。丟棄其餘係數的緣由是它們表明濾波器組係數的快速變化，而且這些精細的細節對自動語音識別（ASR）沒有貢獻。

$$C(n)=\sum\limits_{m=0}^{N-1}{s( m)\cos( \frac{\pi n( m-0.5)}{M} )},n=1,2,...,L$$

L階指MFCC係數階數，一般取2-13。這裏M是三角濾波器個數。

mfcc = dct(filter_banks, type=2, axis=1, norm='ortho')[:, 1 : (num_ceps + 1)] # 保持在2-13

能夠將正弦提高器(Liftering在倒譜域中進行過濾。 注意在譜圖和倒譜圖中分別使用filtering和liftering)應用於MFCC以去強調更高的MFCC，其已經聲稱能夠改善噪聲信號中的語音識別。

(nframes, ncoeff) = mfcc.shape
n = numpy.arange(ncoeff)
lift = 1 + (cep_lifter / 2) * numpy.sin(numpy.pi * n / cep_lifter)
mfcc *= lift

由此產生的MFCC：

 MFCCs

 均值歸一化

如前所述，爲了平衡頻譜並改善信噪比（SNR），咱們能夠簡單地從全部幀中減去每一個係數的平均值。

filter_banks -= (numpy.mean(filter_banks, axis=0) + 1e-8)

均值歸一化濾波器組：

 歸一化濾波器數組

一樣對於MFCC：

mfcc -= (numpy.mean(mfcc, axis=0) + 1e-8)

均值歸一化MFCC：

 標準的MFCC

總結

　　在這篇文章中，咱們探討了計算Mel標度濾波器組和Mel頻率倒譜系數（MFCC）的過程。

　　到目前爲止，計算濾波器組和MFCC的步驟是根據其動機和實現進行討論的。值得注意的是，計算濾波器組所需的全部步驟都是由語音信號的性質和人類對這些信號的感知所驅動的。相反，計算MFCC所需的額外步驟是由一些機器學習算法的限制所驅動的。須要離散餘弦變換（DCT）來去相關濾波器組係數，該過程也稱爲白化。特別是，當高斯混合模型 - 隱馬爾可夫模型（GMMs-HMMs）很是受歡迎時，MFCC很是受歡迎，MFCC和GMM-HMM共同演化爲自動語音識別（ASR）的標準方式²。隨着深度學習在語音系統中的出現，人們可能會質疑MFCC是否仍然是正確的選擇，由於深度神經網絡不太容易受到高度相關輸入的影響，所以離散餘弦變換（DCT）再也不是必要的步驟。值得注意的是，離散餘弦變換（DCT）是一種線性變換，所以在語音信號中丟棄一些高度非線性的信息是不可取的。

　　質疑傅里葉變換是不是必要的操做是明智的。鑑於傅立葉變換自己也是線性運算，忽略它並嘗試直接從時域中的信號中學習多是有益的。實際上，最近的一些工做已經嘗試過，而且報告了積極的結果。然而，傅立葉變換操做是很難學習的操做，可能會增長實現相同性能所需的數據量和模型複雜性。此外，在進行短時傅里葉變換（stft）時，咱們假設信號在這一短期內是平穩的，所以傅里葉變換的線性不會構成一個關鍵問題。

用librosa提取MFCC特徵

　　MFCC特徵是一種在自動語音識別和說話人識別中普遍使用的特徵。在librosa中，提取MFCC特徵只須要一個函數：

librosa.feature.mfcc(y=None, sr=22050, S=None, n_mfcc=20, dct_type=2, norm='ortho', **kwargs)

參數：

• y：音頻時間序列
• sr：音頻的採樣率
• S：np.ndarray，對數功率梅爾頻譜，這個函數既能夠支持時間序列輸入也能夠支持頻譜輸入，都是返回MFCC係數。
• n_mfcc：int>0，要返回的MFCC數，scale濾波器組的個數，一般取22~40，本文取40
• dct_type：None, or {1, 2, 3}  離散餘弦變換（DCT）類型。默認狀況下，使用DCT類型2。
• norm： None or ‘ortho’ 規範。若是dct_type爲2或3，則設置norm =’ortho’使用正交DCT基礎。norm不支持dct_type = 1。

返回：

M：np.ndarray MFCC序列

import librosa

y, sr = librosa.load('./train_nb.wav', sr=16000)
# 提取 MFCC feature
mfccs = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=40)

print(mfccs.shape)        # (40, 65)

繪製頻譜圖

Librosa有顯示頻譜圖波形函數specshow( ):

import matplotlib.pyplot as plt 
import librosa

y, sr = librosa.load('./train_wb.wav', sr=16000)
# 提取 mel spectrogram feature
melspec = librosa.feature.melspectrogram(y, sr, n_fft=1024, hop_length=512, n_mels=128)
logmelspec = librosa.power_to_db(melspec)       # 轉換爲對數刻度
# 繪製 mel 頻譜圖
plt.figure()
librosa.display.specshow(logmelspec, sr=sr, x_axis='time', y_axis='mel')
plt.colorbar(format='%+2.0f dB')        # 右邊的色度條
plt.title('Beat wavform')
plt.show()

參考文獻

MFCC特徵參數提取（一）（基於MATLAB和Python實現）

http://haythamfayek.com/2016/04/21/speech-processing-for-machine-learning.html

Mel Frequency Cepstral Coefficient (MFCC) tutorial

音頻處理庫—librosa的安裝與使用