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



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

2. 梅爾語譜圖具體是如何獲得的？
梅爾語譜圖分為以下幾個步驟。以一段音樂文件為例，詳細(xì)展示每一步的原理和對應(yīng)的Python實(shí)現(xiàn)。

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

import librosa
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rcParams
import matplotlib.ticker as ticker
 
#這是一個畫圖函數(shù)，方便后續(xù)作圖
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函數(shù)中令duration=20來實(shí)現(xiàn)
tmax,tmin = 20,0
t = np.linspace(tmin,tmax,(tmax-tmin)*sr)
personal_plot(t,y[tmin*sr:tmax*sr])



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



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



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

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數(shù)量
 
#未經(jīng)過預(yù)加重的信號頻譜
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)
 

#預(yù)加重之后的信號頻譜
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)


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

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

2.3 分幀(framing)
預(yù)處理完信號之后，要把原信號按時間分成若干個小塊，一塊就叫一幀(frame)。為啥要做這一步？因?yàn)樵盘柛采w的時間太長，用它整個來做FFT，我們只能得到信號頻率和強(qiáng)度的關(guān)系，而失去了時間信息。我們想要得到頻率隨時間變化的關(guān)系，所以將原信號分成若干幀，對每一幀作FFT（又稱為短時FFT，因?yàn)槲覀冎蝗×艘恍《螘r間)，然后將得到的結(jié)果按照時間順序拼接起來。這就是語譜圖(spectrogram)的原理。

下面定義幾個變量:

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

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

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

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

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

具體代碼如下:

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 #不足的部分補(bǔ)零
pad_y = np.append(emphasized_y,np.zeros(pad_frame))
signal_len = signal_length+pad_frame
2.4 加窗(window)
分幀完畢之后，對每一幀加一個窗函數(shù)，以獲得較好的旁瓣下降幅度。通常使用hamming window。

為啥要加窗？要注意，即使我們什么都不加，在分幀的這個過程中也相當(dāng)于給信號加了矩形窗，學(xué)過離散濾波器設(shè)計的人應(yīng)該知道，矩形窗的頻譜有很大的旁瓣，時域中將窗函數(shù)和原函數(shù)相乘，相當(dāng)于頻域的卷積，矩形窗函數(shù)和原函數(shù)卷積之后，由于旁瓣很大，會造成原信號和加窗之后的對應(yīng)部分的頻譜相差很大，這就是頻譜泄露。hamming window有較小的旁瓣，造成的spectral leakage也就較小。代碼實(shí)現(xiàn)如下：首先定義indices變量，這個變量生成每幀所對應(yīng)的sample的索引。np.tile函數(shù)可以使得array從行或者列擴(kuò)展。然后定義frames，對應(yīng)信號在每一幀的值。frames共有1999行，1102列，分別對應(yīng)一共有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中已經(jīng)獲得了frames變量，其每一行對應(yīng)每一幀，所以我們分別對每一行做FFT。由于每一行是1102個點(diǎn)的信號，所以可以選擇1024點(diǎn)FFT（FFT點(diǎn)數(shù)比原信號點(diǎn)數(shù)少會降低頻率分辨率frequency resolution，但這里相差很小，所以可以忽略）。將得到的FFT變換取其magnitude，并進(jìn)行平方再除以對應(yīng)的FFT點(diǎn)數(shù)，即可得到功率譜。到這一步我們其實(shí)已經(jīng)得到了spectrogram, 只需要用plt.imshow畫出其dB值對應(yīng)的熱力圖即可，代碼和結(jié)果如下:

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)
較后一步是將梅爾濾波器運(yùn)用到上一步得到的pow_frames上。所謂梅爾濾波器組是一個等高的三角濾波器組，每個濾波器的起始點(diǎn)在上一個濾波器的中點(diǎn)處。其對應(yīng)的頻率在梅爾尺度上是線性的，因此稱之為梅爾濾波器組。每個濾波器對應(yīng)的頻率可以將較大頻率（下圖中是4000，我們這里是22.05k)用上文中提到的公式轉(zhuǎn)換成梅爾頻率，在梅爾尺度上線性分成若干個頻段，再轉(zhuǎn)換回實(shí)際頻率尺度即可。實(shí)際操作時，將每個濾波器分別和功率譜pow_frames進(jìn)行點(diǎn)乘，獲得的結(jié)果即為該頻帶上的能量(energy)。這里我們的pow_frames是一個(1999,513)的矩陣(這里可能有人疑問513咋來的？我們剛剛做的不是1024點(diǎn)FFT嗎？這里注意因?yàn)槲覀冇昧藃fft，只用了非負(fù)的那一半頻率，所以是1024/2+1個點(diǎn)），梅爾濾波器fbank是一個(mel_N, 513)的矩陣，其中mel_N代表對應(yīng)的梅爾濾波器個數(shù)，這個值不能太大，因?yàn)檫@里我們一共只有513個點(diǎn)，如果mel_N取得太大，會導(dǎo)致前面幾個濾波器的長度都是0 (因?yàn)榈皖l的梅爾濾波器特別窄)。我們只要將這兩個矩陣相乘pow_frames*fbank.T即可得到mel-spectrogram，結(jié)果是一個(1999, 40)的矩陣，每一行是一幀，每一列代表對應(yīng)的梅爾頻帶的能量。具體梅爾濾波器的圖例和計算公式以及對應(yīng)代碼如下:





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

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

#下面定義mel filter
mel_N = 40 #濾波器數(shù)量,這個數(shù)字若要提高，則NFFT也要相應(yīng)提高
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) #將頻率轉(zhuǎn)換成對應(yīng)的sample位置
fbank = np.zeros((mel_N,int(NFFT/2+1))) #每一行儲存一個梅爾濾波器的數(shù)據(jù)
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
機(jī)器學(xué)習(xí)的時候，每一個音頻段即可用對應(yīng)的mel-spectogram表示，每一幀對應(yīng)的某個頻段即為一個feature。因此我們一共獲得了1999*40個feature和對應(yīng)的值。實(shí)際操作中，每個音頻要采用同樣的長度，這樣我們的feature數(shù)量才是相同的。通常還要進(jìn)行歸一化，即每一幀的每個元素要減去該幀的平均值，以保證每一幀的均值均為0.

3. MFCCs原理
得到了梅爾語譜圖，想得到MFCCs就很簡單了。首先，為啥要用MFCCs? 因?yàn)?中得到的梅爾譜系數(shù)是互相關(guān)的，在一些機(jī)器學(xué)習(xí)算法中可能會出問題，因?yàn)橛行┧惴僭O(shè)數(shù)據(jù)不存在互相關(guān)性。因此，可以用DCT變換來壓縮梅爾譜，得到一組不相關(guān)的系數(shù)。DCT在圖像壓縮領(lǐng)域很常見，大家可以自己查閱相關(guān)資料其原理。在語音識別中，得到的梅爾倒頻系數(shù)只保存前2-13個，剩下的不用，因?yàn)檠芯勘砻髌渌禂?shù)代表了系數(shù)中高階的變化，在ASR中沒啥用。

當(dāng)然，更深層次的原因是MFCC是倒譜系數(shù)，所謂倒譜系數(shù)，就是對log之后的梅爾譜系數(shù)進(jìn)行DCT變換，其實(shí)相當(dāng)于將實(shí)際上是頻域的信號當(dāng)成時域信號強(qiáng)行進(jìn)行頻域變換，得到的是頻域信號在偽頻域的幅頻相應(yīng)，前2-13個系數(shù)代表的是包絡(luò)，因?yàn)樗麄冊趥晤l域上是低頻信號，所以在前面，后面的系數(shù)是偽頻域的高頻信號，代表的是spectral details，在語音識別的時候，對我們幫助更大的是包絡(luò)，因?yàn)榘薴ormants等信息。

4. 總結(jié)
總的來說，過去在HMM、GMM等模型用的比較火的時候，多將MFCC用于特征提取，因?yàn)楫?dāng)時的機(jī)器學(xué)習(xí)算法有相應(yīng)的不足。如今較熱門的是以神經(jīng)網(wǎng)絡(luò)為代表的深度學(xué)習(xí)算法，神經(jīng)網(wǎng)絡(luò)內(nèi)部復(fù)雜，在訓(xùn)練的過程中可以在網(wǎng)絡(luò)內(nèi)部將互相關(guān)的問題弱化，也因此DCT變換顯得有些多余，何況還會提高計算量，而且DCT作為一種線性變換，有可能會導(dǎo)致?lián)p失信號中一些非線性信息。因此，如今Mel-spectogram用的更多。