本文内容借鑒于: Speech Processing for Machine Learning: Filter banks, Mel-Frequency Cepstral Coefficients (MFCCs) and What’s In-Between | Haytham Fayek
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秒的波形。現在的音樂檔案采樣率通常是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用的更多。