【MFCC梅爾倒頻譜參數】淺談語音特征參數MFCC提取過程

MFCC提取過程

• • 一、概述
  • 二、提取過程
  • • 提取語音信号
    • 預加重
    • 分幀
    • 加窗
    • 快速傅裡葉變換
    • 梅爾濾波器組
    • 對數能量
    • 離散餘弦變換（DCT）
    • 動态差分參數的提取（包括一階差分和二階差分）

一、概述

• 在語音識别（Speech Recognition）和話者識别（Speaker

  Recognition）方面，最常用到的語音特征就是梅爾倒譜系數（Mel-scale Frequency Cepstral Coefficients，簡稱MFCC）。根據人耳聽覺機理的研究發現，人耳對不同頻率的聲波有不同的聽覺敏感度。從200Hz到5000Hz的語音信号對語音的清晰度影響對大。兩個響度不等的聲音作用于人耳時，則響度較高的頻率成分的存在會影響到對響度較低的頻率成分的感受，使其變得不易察覺，這種現象稱為掩蔽效應。由于頻率較低的聲音在内耳蝸基底膜上行波傳遞的距離大于頻率較高的聲音，故一般來說，低音容易掩蔽高音，而高音掩蔽低音較困難。在低頻處的聲音掩蔽的臨界帶寬較高頻要小。是以，人們從低頻到高頻這一段頻帶内按臨界帶寬的大小由密到疏安排一組帶通濾波器，對輸入信号進行濾波。将每個帶通濾波器輸出的信号能量作為信号的基本特征，對此特征經過進一步處理後就可以作為語音的輸入特征。由于這種特征不依賴于信号的性質，對輸入信号不做任何的假設和限制，又利用了聽覺模型的研究成果。是以，這種參數比基于聲道模型的LPCC相比具有更好的魯邦性，更符合人耳的聽覺特性，而且當信噪比降低時仍然具有較好的識别性能。

• 梅爾倒譜系數（Mel-scale Frequency Cepstral

  Coefficients，簡稱MFCC）是在Mel标度頻率域提取出來的倒譜參數，Mel标度描述了人耳頻率的非線性特性，它與頻率的關系可用下式近似表示：

  

  下圖是梅爾頻率關系圖

  

• 倒譜(Cepstrum): 信号的傅裡葉變化經過對數運算後在進行逆傅裡葉變換得到的譜

  

• 對濾波器組的輸出使用離散餘弦變換(DCT)去除相關性的倒譜系數c(n):

  

二、提取過程

上圖是MFCC的流程圖。

以下的代碼叙述都是以Python語言作為例子

提取語音信号

1. 提取語音信号

   一般python提取音頻信号可以使用标準庫scipy拿到語音信号值。

def read_wav(audio_path):
    """
    使用scipy庫下的方法讀取音頻檔案wav
    時間序列y / 采樣率sr = 音頻時長
    
    :param audio_path: 音頻路徑
    :return: 傳回時間坐标，時間序列y（資料類型和聲道數由檔案本身決定）和 采樣率sr(赫茲)
    """
    from scipy.io.wavfile import read
    from numpy import arange
    sr, y = read(filename=audio_path)  # 讀取音頻檔案，傳回音頻采樣率和時間序列
    x = arange(0, len(y)/sr, 1/sr)  # 總秒數 = 總采樣點/采樣頻率
    return x, y, sr
           

預加重

1. 預加重

   預加重處理其實是将語音信号通過一個高通濾波器：

   Y(n) = X(n) - μX(n-1)

   y(n) 指輸入的是離散信号

   H(z) = 1 - μz^-1

   H(z) 指輸入的是連續信号

   上式中μ的值介于0.9-1.0之間，我們通常取0.97。預加重的目的是提升高頻部分，使信号的頻譜變得平坦，保持在低頻到高頻的整個頻帶中，能用同樣的信噪比求頻譜。同時，也是為了消除發生過程中聲帶和嘴唇的效應，來補償語音信号受到發音系統所抑制的高頻部分，也為了突出高頻的共振峰。

def pre_emphasis(sinal, coefficient=0.97):
	import numpy
	return numpy.append(signal[0], signal[1:] - coefficient * signal[:-1])
           

分幀

1. 分幀

   先将N個采樣點集合成一個觀測機關，稱為幀。通常情況下N的值為256或512，涵蓋的時間約為20~30ms左右。為了避免相鄰兩幀的變化過大，是以會讓兩相鄰幀之間有一段重疊區域，此重疊區域包含了M個取樣點，通常M的值約為N的1/2或1/3。通常語音識别所采用語音信号的采樣頻率為8KHz或16KHz，以8KHz來說，若幀長度為256個采樣點，則對應的時間長度是256/8000×1000=32ms。

def audio2frame(signal, frame_length, frame_step, winfunc=lambda x: numpy.ones((x,))):
    """将音頻信号轉化為幀。
	參數含義：
	signal:原始音頻型号
	frame_length:每一幀的長度(這裡指采樣點的長度，即采樣頻率乘以時間間隔)
	frame_step:相鄰幀的間隔（同上定義）
	winfunc:lambda函數，用于生成一個向量
    """
    signal_length = len(signal)  # 信号總長度
    frame_length = int(round(frame_length))  # 以幀幀時間長度
    frame_step = int(round(frame_step))  # 相鄰幀之間的步長
    if signal_length <= frame_length:  # 若信号長度小于一個幀的長度，則幀數定義為1
        frames_num = 1
    else:  # 否則，計算幀的總長度
        frames_num = 1+int(math.ceil((1.0*signal_length-frame_length)/frame_step))
    pad_length = int((frames_num-1)*frame_step+frame_length)  # 所有幀加起來總的鋪平後的長度, 最後一幀長度當滿長度處理,幀與幀之間的相同部分不算
    zeros = numpy.zeros((pad_length-signal_length,))  # 不夠的長度使用0填補，類似于FFT中的擴充數組操作
    pad_signal = numpy.concatenate((signal, zeros))  # numpy的矩陣拼接, 填補後的信号記為pad_signal
    # numpy.tile(x, y) 生成新數組，新數組為複制y次x數組->一行中有y個x的數組； y=[x0, y0]可為數組時->有x0行的x數組，y0列的x數組
    indices = numpy.tile(numpy.arange(0, frame_length), (frames_num, 1))+numpy.tile(numpy.arange(0, frames_num * frame_step, frame_step), (frame_length, 1)).T  # 相當于對所有幀的時間點進行抽取，得到frames_num*frame_length長度的矩陣
    indices = numpy.array(indices, dtype=numpy.int32)  # 将indices轉化為矩陣
    frames = pad_signal[indices]  # 得到幀信号
    win = numpy.tile(winfunc(frame_length), (frames_num, 1))  # window窗函數，這裡預設取1
    return frames*win   # 傳回幀信号矩陣
           

• 下面來詳解以下這個算法, 本人會以表格形式來展現信号，設一個信号長度為11，時間序列值為1-11的信号，并且令幀長度為4，幀間隔為2.

  ∵幀間隔一般是幀長度的1/2或1/3,在這裡取1/2

  (1) 依照算法，123行代碼所得到的值依次是11、4、2

  (2) 接下來通過判斷信号總長度與幀長度大小：

  如果信号長度小于幀長度，說明一個幀就足夠了，幀總數指派為1。否則計算幀。

  

  前面的1的作用是将最後一幀的長度若不夠幀長度的則按滿長度處理，因為除法得到的結果是小數，則進行取整，這會導緻缺少一幀，固先把最後一幀當作滿處理，再通過減一個幀長度，就得到的總幀數。在這裡得到5

  (3) 計算處理後的總長度，就是每幀按滿幀處理，其中幀與幀的重複部分隻取一次,是以在這裡得到值12

  (4) 通過計算處理後的長度與原信号長度的差，拿到需要補0的數目。在這裡長度差為12 - 11 = 1，需要補一個0，其中在這裡先将信号轉為矩陣，再向信号矩陣後面添加0矩陣。

  

  (5) indices的生成代碼:

  第一part: 數組[0, …, 幀長度] 的數組擴充至 [幀總數，1]的矩陣。 在這裡第一部分的矩陣為:

  

  第二部分：數組生成規則是，從0到幀總數*幀步長的範圍内，每隔一個幀步長就取一個數。并将該數組複制至[幀長度，1]的矩陣最後取逆矩陣。在這裡第二部分的矩陣為:

  

  取逆矩陣後為：

  

  最後第一部分與第二部分進行矩陣加法:

  

  在這裡就顯而易見了，每一行代表一個幀序列，一共有5個幀，且每一個幀為4個幀長度，幀的值為信号矩陣的索引值，即signal[0]， signal[1]， …

  (6) 将幀生成二維數組轉化為矩陣

  (7) 将拿到的幀矩陣通過signal的值轉化為幀信号

  

  (8) 窗矩陣生成 (一般選用的窗函數是漢明窗，下面有詳解介紹)

  預設不使用窗函數，在這裡即窗函數為[1, 1, 1, 1]

  再将窗函數複制擴充至[幀總數行，幀長度列]

  在這裡即為

  

  若窗函數為[2, 2, 2, 2]，則意義是将每一幀的信号擴大1倍

  (9) 使用窗矩陣點乘幀信号矩陣，對每一個幀都通過窗函數進行處理。

  (10) 最後傳回處理完的幀信号矩陣

  在這裡，最終傳回的信号是

  

加窗

1. 加窗(預設是漢明窗)

   将每一幀乘以漢明窗，以增加幀左端和右端的連續性。假設分幀後的信号為S(n), n=0,1…,N-1, N為幀數目，那麼乘上漢明窗後,W(n)形式如下：

   

   不同的a值會産生不同的漢明窗，一般情況下a取0.46.

def hanming(n, a, N):
    import numpy as np
    return (1 - a) - a * np.cos((2 * np.pi * n) / (N - 1))
           

# 漢明窗展示
from matplotlib.pyplot import plot, show


def hanming(n, a, N):
    import numpy as np
    return (1 - a) - a * np.cos((2 * np.pi * n) / (N - 1))


N = 100
a = 0.46
x = linspace(0, 1, N)	# 頻率為N
y = [hanming(i, a, N) for i in range(N)]
plot(x, y)
show()
           

• N越大，曲線越平滑。随着a的變化，趨勢如下圖.

  

  漢明窗的特性是窗内中間的資料極大程度被保留，而邊緣的兩邊被削弱。是以有上面的幀步長為1/2的幀長度，将窗内邊緣的部分通過下個幀表現出來。

  而為什麼要加窗呢？

  在信号進行中，可以說加窗處理是一個必經的過程，因為我們的計算機隻能處理有限長度的信号，是以原始信号X(t)要以T（采樣時間）截斷，即有限化，成為XT(t)後再進一步處理，這個過程式就是加窗處理，但什麼時候用什麼窗呢？這時我們就要對所需用到的函數窗做一定的了解。在平時，我們用得最多的是矩形窗。實際的信号處理過程中，我們用的矩形窗在邊緣處将信号突然截斷時，窗外時域資訊全部消失，導緻在頻域增加了頻率分量的現象，即頻譜洩漏。避免洩漏的最佳方法是滿足整周期采樣條件，但實際中是不可能做到的。對于非整周期采樣的情況，必須考慮如何減少加窗時造成的洩漏誤差，主要的措施是使用合理的加窗函數，使信号截斷的銳角鈍化，進而使頻譜的擴散減到最少。

  首先介紹一下為什麼要用函數窗：函數窗的主要用于對截斷處的不連續變化進行平滑，減少洩漏。此外，加窗處理還有很多其它的原因，如減少噪聲幹擾、限定測試的持續時間、從頻率接近的信号中分離出幅值不同的信号……

  常見的幾種窗的基本名額：

  

  

  下面是通過生成一個不同長度的餘弦函數，對餘弦函數進行快速傅裡葉變換以及對加了漢明窗的餘弦函數進行快速傅裡葉變換，進行對比。在這裡，長度取值N我分别取值了10、100、500、1000，步長為0.1.(在這裡我并沒有對x坐标進行頻率變換)

  每一行代表的是同一個長度的餘弦函數，第一列是原函數，第二列是原函數的快速傅裡葉變換圖，第三列是加漢明窗後的函數，第四列是加漢明窗後函數的快速傅裡葉變換。

  

  在這裡x坐标軸我沒有做頻率的變換，這是示範一下圖像在加窗前後的傅裡葉變換的差别。

  通過上圖可以看到，在n等于10的時候，圖像已經完全失去了cos的頻率特性，是以在原函數的FFT圖像中，頻譜已經失真了，但是經過加窗後，頻譜竟然有些修複。(cos函數的頻譜圖不清楚的建議去看一下信号與系統這本書或者百度一下)。n=100時，其頻譜圖也有些失真，但是經過加窗後，失真現象也緩和了。到n=500, n=1000時，因為時間序列已經比較長了，是以原函數的FFT的圖像與原FFT圖像盡管有些誤差，但是誤差已經可以算是可以允許了。這說明，漢明窗的确是能緩和頻譜洩露現象。

快速傅裡葉變換

1. 快速傅裡葉變換

   由于信号在時域上的變換通常很難看出信号的特性，是以通常将它轉換為頻域上的能量分布來觀察，不同的能量分布，就能代表不同語音的特性，這一轉換可以通過傅裡葉變換進行。是以在乘上漢明窗後，每幀還必須再經過快速傅裡葉變換以得到在頻譜上的能量分布，快速傅裡葉變換與傅裡葉變換的差別可以看作是運算量的降低。對分幀加窗後的各幀信号進行快速傅裡葉變換得到各幀的頻譜後，由于傅裡葉變換後的值是含有虛數部分的，是以對語音信号的頻譜取模平方得到語音信号的功率譜。設語音信号的DFT為：

   

   式中x(n)為輸入的語音信号，N表示傅裡葉變換的點數。

   代碼為：

from scipy import fft
result = fft(signal)	# signal是指信号的時間序列
           

梅爾濾波器組

1. 三角帶通濾波器

   将能量譜通過一組Mel尺度的三角形濾波器組，定義一個有M個濾波器的濾波器組（濾波器的個數和臨界 帶的個數相近），采用的濾波器為三角濾波器，中心頻率為f(m),m=1,2,…,M。M通常取22-26。各f(m)之間的間隔随着m值的減小而縮小，随着m值的增大而增寬，如圖所示：

   

   可以從三角波函數可以看出，随着m的增大，幅值越來越小，這從頻域上看，就是顯低頻，越到高頻的部分越衰減的厲害。

   三角濾波器的頻率響應定義為：

   

   三角帶通濾波器有兩個主要目的：

   對頻譜進行平滑化，并消除高次諧波的作用(高次諧波也就是高頻分量)，突顯原先語音的共振峰。（是以一段語音的音調或音高，是不會呈現在 MFCC 參數内，換句話說，以 MFCC 為特征的語音辨識系統，并不會受到輸入語音的音調不同而有所影響） 此外，還可以降低運算量。

對數能量

1. 計算每個濾波器組輸出的對數能量

   

   X(t)是信号序列， Hm是梅爾濾波器組

   對數能量：一幀内信号的平方和，再取以10為底的對數值，再乘以10

   此外，一幀的音量（即能量），也是語音的重要特征。

   注：若要加入其它語音特征以測試識别率，也可以在此階段加入，這些常用的其它語音特征包含音高、過零率以及共振峰等。

離散餘弦變換（DCT）

1. 離散餘弦變換（DCT）得到MFCC系數

   

   将上述的對數能量帶入離散餘弦變換，求出L階的Mel-scale Cepstrum參數。L階指MFCC系數階數，通常取12-16。這裡m是三角濾波器的階數。

動态差分參數的提取（包括一階差分和二階差分）

1. 動态差分參數的提取（包括一階差分和二階差分）

   标準的倒譜參數MFCC隻反映了語音參數的靜态特性，語音的動态特性可以用這些靜态特征的差分譜來描述。實驗證明：把動、靜态特征結合起來才能有效提高系統的識别性能。差分參數的計算可以采用下面的公式：

   

   式中,dt表示第t個一階差分，Ct表示第t個倒譜系數，Q表示倒譜系數的階數，K表示一階導數的時間差，可取1或2。将上式的結果再代入就可以得到二階差分的參數。

參考文獻或部落格：

[1]說話人識别中改進特征提取算法的研究[J].宋樂.太原理工大學

[2]基于時頻分布與MFCC的說話人識别[J]. 金銀燕,于鳳芹,何豔. 計算機系統應用. 2012(04)

[3]語音識别之——mfcc什麼是漢明窗，為什麼加漢明窗[CSDN部落格].gxiaoyaya

[4]語音識别的第一步MFCC特征提取代碼（Python）[ITKeyword社群].斷橋殘雪斷橋殘雪

[5]語音識别之MFCC特征提取[百度文庫].湯旭國

[6]MFCC [百度文庫]