梅爾倒譜系數（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的步驟
- 將信號幀化為短幀
- 對於每一幀，計算功率譜的周期圖估計
- 將mel濾波器組應用於功率譜，求濾波器組的能量，將每個濾波器中的能量相加
- 取所有濾波器組能量的對數
- 取對數濾波器組能量的離散余弦變換（DCT）。
- 保持DCT系數2-13，其余部分丟棄
通常還有其他事情要做，有時會將幀能量附加到每個特征向量上。通常還會附加Delta和Delta-Delta 特征。提升也通常應用於最終特征。
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）、預加重（Pre-Emphasis）
對信號應用預加重濾波器以放大高頻，預加重濾波器在以下幾個方面很有用：
- 平衡頻譜，因為高頻通常與較低頻率相比具有較小的幅度，
- 避免在傅里葉變換操作操作過程中出現數值問題
- 改善信號 - 噪聲比（SNR）
- 消除發聲過程中聲帶和嘴唇的效應，來補償語音信號受到發音系統所抑制的高頻部分，也為了突出高頻的共振峰。
預加重處理其實是將語音信號通過一個高通濾波器：
$$y(t) = x(t) - \alpha x(t-1)$$
其中濾波器系數（$\alpha$）的通常為0.95或0.97，這里取pre_emphasis =0.97：
emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])
預加重后的時域信號
題外話：預加重在現代系統中的影響不大，主要是因為除避免在現代FFT實現中不應成為問題的FFT數值問題，大多數預加重濾波器的動機都可以通過均值歸一化來實現（在本文后面討論）。 在現代FFT實現中。
（2）、分幀（Framing）
在預加重之后，我們需要將信號分成短時幀。因此在大多數情況下，語音信號是非平穩的，對整個信號進行傅里葉變換是沒有意義的，因為我們會隨着時間的推移丟失信號的頻率輪廓。語音信號是短時平穩信號。因此我們在短時幀上進行傅里葉變換，通過連接相鄰幀來獲得信號頻率輪廓的良好近似。
將信號幀化為20-40 ms幀。標准是25毫秒 frame_size = 0.025。這意味着8kHz信號的幀長度為0.025 * 8000 = 200個采樣點。幀移通常為10毫秒 frame_stride = 0.01（80個采樣點），為了避免相鄰兩幀的變化過大，因此會讓兩相鄰幀之間有一段重疊區域，通常約為每幀語音的1/2或1/3或50%(+/-10%)，我們設置為15毫秒 overlap=0.015，因此此重疊區域包含了0.015*8000=120個取樣點。第一個語音幀0開始，下一個語音幀從80開始，直到到達語音文件的末尾。如果語音文件沒有划分為偶數個幀，則用零填充它以使其完成。
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)]
（3）、加窗（Window）
將信號分割成幀后，我們再對每個幀乘以一個窗函數，如Hamming窗口。以增加幀左端和右端的連續性。抵消FFT假設（數據是無限的），並減少頻譜泄漏。漢明窗的形式如下：
$$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$$，N是窗口長度，我們這里假設$a=0.46$
frames *= numpy.hamming(frame_length) # frames *= 0.54 - 0.46 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1)) # 內部實現
二、FFT（Fourier-Transform）
由於信號在時域上的變換通常很難看出信號的特性，通常對它做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的幅度(magnitude)
三、功率譜（Power Spectrum）
然后我們使用以下公式計算功率譜（周期圖periodogram），對語音信號的頻譜取模平方（取對數或者去平方，因為頻率不可能為負，負值要舍去）得到語音信號的譜線能量。
$$P = \frac{|FFT(x_i)|^2}{N}$$
其中，$X_i$是信號X的第$i$幀，這可以用以下幾行來實現：
pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2)) # 功率譜
四、濾波器組（Filter Banks）
計算Mel濾波器組，將功率譜通過一組Mel刻度(通常取40個濾波器，nfilt=40)的三角濾波器(triangular filters)來提取頻帶(frequency bands)。
這個Mel濾波器組就像人類的聽覺感知系統（耳朵），人耳只關注某些特定的頻率分量（人的聽覺對頻率是有選擇性的）。它對不同頻率信號的靈敏度是不同的，換言之，它只讓某些頻率的信號通過，而壓根就直接無視它不想感知的某些頻率信號。但是這些濾波器在頻率坐標軸上卻不是統一分布的，在低頻區域有很多的濾波器，他們分布比較密集，但在高頻區域，濾波器的數目就變得比較少，分布很稀疏。因此Mel刻度的目的是模擬人耳對聲音的非線性感知，在較低的頻率下更具辨別力，在較高的頻率下則不具辨別力。我們可以使用以下公式在赫茲（f）和梅爾（m）之間進行轉換：
我們可以用下面的公式，在語音頻率和Mel頻率間轉換
- 從頻率轉換為梅爾刻度的公式為：$f_{mel}=2595*\log _{10}(1+\frac{f}{700})$
- 從梅爾回到頻率：$f = 700 (10^{f_{mel}/2595} - 1)$
定義一個有M個三角濾波器的濾波器組（濾波器的個數和臨界帶的個數相近），M通常取22-40，26是標准，本文取nfilt = 40。濾波器組中的每個濾波器都是三角形的，中心頻率為f(m) ，中心頻率處的響應為1，並向0線性減小，直到達到兩個相鄰濾波器的中心頻率，其中響應為0，各f(m)之間的間隔隨着m值的增大而增寬，如圖所示：
Mel-Scale上的Filter Bank
這可以通過以下等式建模，三角濾波器的頻率響應定義為:
$$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個能量值
nfilt = 40 low_freq_mel = 0 high_freq_mel = (2595 * np.log10(1 + (sample_rate / 2) / 700)) # 將Hz轉換為Mel # 我們要做40個濾波器組，為此需要42個點，這意味着在們需要low_freq_mel和high_freq_mel之間線性間隔40個點 mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2) # 使得Mel scale間距相等 hz_points = (700 * (10 ** (mel_points / 2595) - 1)) # 將Mel轉換回-Hz # bin = sample_rate/NFFT # frequency bin的計算公式 # bins = hz_points/bin=hz_points*NFFT/ sample_rate # 得出每個hz_point中有多少frequency bin bins = 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(bins[m - 1]) # 左 f_m = int(bins[m]) # 中 f_m_plus = int(bins[m + 1]) # 右 for k in range(f_m_minus, f_m): fbank[m - 1, k] = (k - bins[m - 1]) / (bins[m] - bins[m - 1]) for k in range(f_m, f_m_plus): fbank[m - 1, k] = (bins[m + 1] - k) / (bins[m + 1] - bins[m]) 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) # dB
信號的功率譜經過 Filter Bank 后，得到的譜圖為：
信號的頻譜圖
如果經過Mel scale濾波器組是所需的特征，那么我們可以跳過下一步。
五、梅爾頻率倒譜系數（MFCCs）
上一步驟中計算的濾波器組系數是高度相關的，這在某些機器學習算法中可能是有問題的。因此，我們可以應用離散余弦變換（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
六、均值歸一化（Mean Normalization）
如前所述，為了平衡頻譜並改善信噪比（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）的過程。
我們從動機和實現的角度討論了計算Filter Banks和MFCCs的步驟。值得注意的是，計算濾波器組所需的所有步驟都是由語音信號的性質和人類對這些信號的感知所驅動的。相反，計算MFCC所需的額外步驟是由一些機器學習算法的限制所驅動的。需要離散余弦變換（DCT）來對濾波器組系數去相關化，該過程也稱為白化。特別是，當高斯混合模型 - 隱馬爾可夫模型（GMMs-HMMs）非常受歡迎時，MFCC非常受歡迎，MFCC和GMM-HMM共同演化為自動語音識別（ASR）的標准方式2。隨着深度學習在語音系統中的出現，人們可能會質疑MFCC是否仍然是正確的選擇，因為深度神經網絡不太容易受到高度相關輸入的影響，因此離散余弦變換（DCT）不再是必要的步驟。值得注意的是，離散余弦變換（DCT）是一種線性變換，因此在語音信號中丟棄一些高度非線性的信息是不可取的。
質疑傅里葉變換是否是必要的操作是明智的。鑒於傅立葉變換本身也是線性運算，忽略它並嘗試直接從時域中的信號中學習可能是有益的。實際上，最近的一些工作已經嘗試過，並且報告了積極的結果。然而，傅立葉變換操作是很難學習的操作，可能會增加實現相同性能所需的數據量和模型復雜性。此外，在進行短時傅里葉變換（stft）時，我們假設信號在這一短時間內是平穩的，因此傅里葉變換的線性不會構成一個關鍵問題。
如果機器學習算法不易受高度相關輸入的影響，則使用Mel刻度濾波器組。如果機器學習算法易受相關輸入的影響，則使用MFCCs。
用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 import librosa.display 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', cmp="jet") plt.colorbar(format='%+2.0f dB') # 右邊的色度條 plt.title('Beat wavform') plt.show()
參考文獻
MFCC特征參數提取（一）（基於MATLAB和Python實現）