1. 什么是梅爾語譜圖和梅爾倒頻系數？
機器學習的第一步都是要提取出相應的特征(feature)，如果輸入數據是圖片，例如28*28的圖片，那么只需要把每個像素(pixel)作為特征，對應的像素值大小(代表顏色的強度)作為特征值即可。那么在音頻、語音信號處理領域，我們需要將信號轉換成對應的語譜圖(spectrogram)，將語譜圖上的數據作為信號的特征。語譜圖的橫軸x為時間，縱軸y為頻率，(x,y)對應的數值代表在時間x時頻率y的幅值。通常的語譜圖其頻率是線性分布的，但是人耳對頻率的感受是對數的(logarithmic)，即對低頻段的變化敏感，對高頻段的變化遲鈍，所以線性分布的語譜圖顯然在特征提取上會出現“特征不夠有用的情況”，因此梅爾語譜圖應運而生。梅爾語譜圖的縱軸頻率和原頻率經過如下公式互換:



其中f代表原本的頻率，m代表轉換后的梅爾頻率，顯然，當f很大時，m的變化趨于平緩。而梅爾倒頻系數(MFCCs)是在得到梅爾語譜圖之后進行余弦變換(DCT,一種類似于傅里葉變換的線性變換），然后取其中一部分系數即可。

2. 梅爾語譜圖具體是如何獲得的？
梅爾語譜圖分為以下幾個步驟。以一段音樂文件為例，詳細展示每一步的原理和對應的Python實現。

2.1 獲取音頻信號
python可以用librosa庫來讀取音頻文件，但是對于MP3文件，它會自動調用audio_read函數，所以如果是MP3文件，務必保證將ffmpeg.exe的路徑添加到系統環境變量中，不然audio_read函數會出錯。這里我們首先讀取音頻文件，并作出0-20秒的波形?，F在的音樂文件采樣率通常是44.1kHz。用y和sr分別表示信號和采樣率。代碼和圖形如下：

import librosa
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import matplotlib.ticker as ticker
 
#這是一個畫圖函數，方便后續作圖
def personal_plot(x,y):
    plt.figure(dpi=200,figsize=(12,6))
    rcParams['font.family']='Comic Sans MS'
    plt.plot(x,y)
    plt.xlim(x[0],x[-1])
    plt.xlabel('time/s',fontsize=20)
    plt.ylabel('Amplitude',fontsize=20)
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    plt.grid()
 
#注意如果文件名不加路徑，則文件必須存在于python的工作目錄中
y,sr = librosa.load('笑顏.mp3',sr=None)
 
#這里只獲取0-20秒的部分，這里也可以在上一步的load函數中令duration=20來實現
tmax,tmin = 20,0
t = np.linspace(tmin,tmax,(tmax-tmin)*sr)
personal_plot(t,y[tmin*sr:tmax*sr])



2.2 信號預加重(pre-emphasis)
通常來講語音/音頻信號的高頻分量強度較小，低頻分量強度較大，信號預加重就是讓信號通過一個高通濾波器，讓信號的高低頻分量的強度不至于相差太多。在時域中，對信號x[n]作如下操作：



α通常取一個很接近1的值，typical value為0.97或0.95. 從時域公式來看，可能有部分人不懂為啥這是一個高通濾波器，我們從z變換的角度看一下濾波器的transfer function：



可以看出濾波器有一個極點0，和一個零點α。當頻率為0時，z=1, 放大系數為(1-α）。當頻率漸漸增大，放大系數不斷變大，當頻率到pi時，放大系數為(1+α)。離散域中，[0,pi]對應連續域中的[0, fs/2](單位Hz)。其中fs為采樣率，在我們這里是44.1kHz。因此當頻率到22000Hz時，放大系數為(1+α)。下面用兩段代碼和對應的圖像給出一個直觀感受:

alpha = 0.97
emphasized_y = np.append(y[tmin*sr],y[tmin*sr+1:tmax*sr]-alpha*y[tmin*sr:tmax*sr-1])
n = int((tmax-tmin)*sr) #信號一共的sample數量
 
#未經過預加重的信號頻譜
plt.figure(dpi=300,figsize=(7,4))
freq = sr/n*np.linspace(0,n/2,int(n/2)+1)
plt.plot(freq,np.absolute(np.fft.rfft(y[tmin*sr:tmax*sr],n)**2)/n)
plt.xlim(0,5000)
plt.xlabel('Frequency/Hz',fontsize=14)
 

#預加重之后的信號頻譜
plt.figure(dpi=300,figsize=(7,4))
plt.plot(freq,np.absolute(np.fft.rfft(emphasized_y,n)**2)/n)
plt.xlim(0,5000)
plt.xlabel('Frequency/Hz',fontsize=14)


這兩段代碼里用了函數librosa.fft.rfft(y,n)，rfft表示經過fft變換之后只取其中一半（因為另一半對應負頻率，沒有用處）, y對應信號，n對應要做多少點的FFT。我們這里的信號有44.1k*20=882000個點，所以對應的FFT 也做882000點的FFT，每一個點所對應的實際頻率是該點的索引值*fs/n，這是咋得出來的？因為第882000個點應該對應(約等于)fs(或者離散域中的2pi)，所以前面的點根據線性關系一一對應即可。這里只展示0-5000Hz，可以看出，經過預加重之后的信號高頻分量明顯和低頻分量的差距沒那么大了。

這樣預加重的好處有什么?原文提到了三點:(1)就是我們剛剛提到的平衡一下高頻和低頻 (2)避免FFT中的數值問題(也就是高頻值太小出現在分母的時候可能會出問題) (3)或許可以提高SNR。

2.3 分幀(framing)
預處理完信號之后，要把原信號按時間分成若干個小塊，一塊就叫一幀(frame)。為啥要做這一步？因為原信號覆蓋的時間太長，用它整個來做FFT，我們只能得到信號頻率和強度的關系，而失去了時間信息。我們想要得到頻率隨時間變化的關系，所以將原信號分成若干幀，對每一幀作FFT（又稱為短時FFT，因為我們只取了一小段時間)，然后將得到的結果按照時間順序拼接起來。這就是語譜圖(spectrogram)的原理。

下面定義幾個變量:

frame_size: 每一幀的長度。通常取20-40ms。太長會使時間上的分辨率(time resolution)較小，太小會加重運算成本。這里取25ms.

frame_length: 每一幀對應的sample數量。等于fs*frame_size。我們這里是44.1k*0.025=1102.

frame_stride: 相鄰兩幀的間隔。通常間隔必須小于每一幀的長度，即兩幀之間要有重疊，否則我們可能會實去兩幀邊界附近的信息。做特征提取的時候，我們是絕不希望實去有用信息的。 這里取10ms，即有60%的重疊。

frame_step: 相鄰兩幀的sample數量。這里是441.

frame_num: 整個信號所需要的幀數。一般希望所需要的幀數是個整數值，所以這里要對信號補0(zero padding)讓信號的長度正好能分成整數幀。

具體代碼如下:

frame_size, frame_stride = 0.025,0.01
frame_length, frame_step = int(round(sr*frame_size)),int(round(sr*frame_stride))
signal_length = (tmax-tmin)*sr
frame_num = int(np.ceil((signal_length-frame_length)/frame_step))+1 #向上舍入
pad_frame = (frame_num-1)*frame_step+frame_length-signal_length #不足的部分補零
pad_y = np.append(emphasized_y,np.zeros(pad_frame))
signal_len = signal_length+pad_frame
2.4 加窗(window)
分幀完畢之后，對每一幀加一個窗函數，以獲得較好的旁瓣下降幅度。通常使用hamming window。

為啥要加窗？要注意，即使我們什么都不加，在分幀的這個過程中也相當于給信號加了矩形窗，學過離散濾波器設計的人應該知道，矩形窗的頻譜有很大的旁瓣，時域中將窗函數和原函數相乘，相當于頻域的卷積，矩形窗函數和原函數卷積之后，由于旁瓣很大，會造成原信號和加窗之后的對應部分的頻譜相差很大，這就是頻譜泄露。hamming window有較小的旁瓣，造成的spectral leakage也就較小。代碼實現如下：首先定義indices變量，這個變量生成每幀所對應的sample的索引。np.tile函數可以使得array從行或者列擴展。然后定義frames，對應信號在每一幀的值。frames共有1999行，1102列，分別對應一共有1999幀和每一幀有1102個sample。將得到的frames和hamming window直接相乘即可，注意這里不是矩陣乘法。

indices = np.tile(np.arange(0, frame_length), (frame_num, 1)) + np.tile(
    np.arange(0, frame_num * frame_step, frame_step), (frame_length, 1)).T
frames = pad_y[indices] #frame的每一行代表每一幀的sample值
frames *= np.hamming(frame_length) #加hamming window 注意這里不是矩陣乘法
2.5 獲取功率譜
我們在2.4中已經獲得了frames變量，其每一行對應每一幀，所以我們分別對每一行做FFT。由于每一行是1102個點的信號，所以可以選擇1024點FFT（FFT點數比原信號點數少會降低頻率分辨率frequency resolution，但這里相差很小，所以可以忽略）。將得到的FFT變換取其magnitude，并進行平方再除以對應的FFT點數，即可得到功率譜。到這一步我們其實已經得到了spectrogram, 只需要用plt.imshow畫出其dB值對應的熱力圖即可，代碼和結果如下:

NFFT = 1024 #frame_length=1102，所以用1024足夠了
mag_frames = np.absolute(np.fft.rfft(frames,NFFT))
pow_frames = mag_frames**2/NFFT
 
plt.figure(dpi=300,figsize=(12,6))
plt.imshow(20*np.log10(pow_frames[40:].T),cmap=plt.cm.jet,aspect='auto')
plt.yticks([0,128,256,384,512],np.array([0,128,256,384,512])*sr/NFFT)


2.6 梅爾濾波器組(Mel-filter banks)
較后一步是將梅爾濾波器運用到上一步得到的pow_frames上。所謂梅爾濾波器組是一個等高的三角濾波器組，每個濾波器的起始點在上一個濾波器的中點處。其對應的頻率在梅爾尺度上是線性的，因此稱之為梅爾濾波器組。每個濾波器對應的頻率可以將較大頻率（下圖中是4000，我們這里是22.05k)用上文中提到的公式轉換成梅爾頻率，在梅爾尺度上線性分成若干個頻段，再轉換回實際頻率尺度即可。實際操作時，將每個濾波器分別和功率譜pow_frames進行點乘，獲得的結果即為該頻帶上的能量(energy)。這里我們的pow_frames是一個(1999,513)的矩陣(這里可能有人疑問513咋來的？我們剛剛做的不是1024點FFT嗎？這里注意因為我們用了rfft，只用了非負的那一半頻率，所以是1024/2+1個點），梅爾濾波器fbank是一個(mel_N, 513)的矩陣，其中mel_N代表對應的梅爾濾波器個數，這個值不能太大，因為這里我們一共只有513個點，如果mel_N取得太大，會導致前面幾個濾波器的長度都是0 (因為低頻的梅爾濾波器特別窄)。我們只要將這兩個矩陣相乘pow_frames*fbank.T即可得到mel-spectrogram，結果是一個(1999, 40)的矩陣，每一行是一幀，每一列代表對應的梅爾頻帶的能量。具體梅爾濾波器的圖例和計算公式以及對應代碼如下:





其中m代表濾波器的序號，f(m-1)和f(m)、f(m+1)分別對應第m個濾波器的起始點、中間點和結束點。大家一定要注意的一點是，這里的f(m)對應的值不是頻率值，而是對應的sample的索引！比如，我們這里較大頻率是22050 Hz, 所以22050Hz對應的是第513個sample，即頻率f所對應的值是f/fs*NFFT。

代碼中有一段np.where(condition,a,b)，這個函數的功能是檢索b中的元素，當condition滿足的時候，輸出a否則，輸出b中的原元素。這一步的操作是為了將其中的全部0值以一個很小的非負值代替，否則在計算dB的時候，log中出現0會出錯。

#下面定義mel filter
mel_N = 40 #濾波器數量,這個數字若要提高，則NFFT也要相應提高
mel_low, mel_high = 0, (2595*np.log10(1+(sr/2)/700))
mel_freq = np.linspace(mel_low,mel_high,mel_N+2)
hz_freq = (700 * (10**(mel_freq / 2595) - 1))
bins = np.floor((NFFT)*hz_freq/sr) #將頻率轉換成對應的sample位置
fbank = np.zeros((mel_N,int(NFFT/2+1))) #每一行儲存一個梅爾濾波器的數據
for m in range(1, mel_N + 1):
    f_m_minus = int(bins[m - 1])   # left
    f_m = int(bins[m])             # center
    f_m_plus = int(bins[m + 1])    # right
 
    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.matmul(pow_frames, fbank.T)
filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)  # np.finfo(float)是較小正值
filter_banks = 20 * np.log10(filter_banks)  # dB
#filter_banks -= np.mean(filter_banks,axis=1).reshape(-1,1)
plt.figure(dpi=300,figsize=(12,6))
plt.imshow(filter_banks[40:].T, cmap=plt.cm.jet,aspect='auto')
plt.yticks([0,10,20,30,39],[0,1200,3800,9900,22000])

較后，得到的mel-spectrogram如下:



2.7 Mel-spectogram feature
機器學習的時候，每一個音頻段即可用對應的mel-spectogram表示，每一幀對應的某個頻段即為一個feature。因此我們一共獲得了1999*40個feature和對應的值。實際操作中，每個音頻要采用同樣的長度，這樣我們的feature數量才是相同的。通常還要進行歸一化，即每一幀的每個元素要減去該幀的平均值，以保證每一幀的均值均為0.

3. MFCCs原理
得到了梅爾語譜圖，想得到MFCCs就很簡單了。首先，為啥要用MFCCs? 因為2中得到的梅爾譜系數是互相關的，在一些機器學習算法中可能會出問題，因為有些算法假設數據不存在互相關性。因此，可以用DCT變換來壓縮梅爾譜，得到一組不相關的系數。DCT在圖像壓縮領域很常見，大家可以自己查閱相關資料其原理。在語音識別中，得到的梅爾倒頻系數只保存前2-13個，剩下的不用，因為研究表明其他系數代表了系數中高階的變化，在ASR中沒啥用。

當然，更深層次的原因是MFCC是倒譜系數，所謂倒譜系數，就是對log之后的梅爾譜系數進行DCT變換，其實相當于將實際上是頻域的信號當成時域信號強行進行頻域變換，得到的是頻域信號在偽頻域的幅頻相應，前2-13個系數代表的是包絡，因為他們在偽頻域上是低頻信號，所以在前面，后面的系數是偽頻域的高頻信號，代表的是spectral details，在語音識別的時候，對我們幫助更大的是包絡，因為包含了formants等信息。

4. 總結
總的來說，過去在HMM、GMM等模型用的比較火的時候，多將MFCC用于特征提取，因為當時的機器學習算法有相應的不足。如今較熱門的是以神經網絡為代表的深度學習算法，神經網絡內部復雜，在訓練的過程中可以在網絡內部將互相關的問題弱化，也因此DCT變換顯得有些多余，何況還會提高計算量，而且DCT作為一種線性變換，有可能會導致損失信號中一些非線性信息。因此，如今Mel-spectogram用的更多。