【語音-01】Fbank和MFCC介紹-理論和代碼

目錄

1. 簡介
2. Fbank處理過程
3. MFCC
4. fbank與mfcc的標准化
5. fbank與mfcc的比較

一、簡介

Fbank：FilterBank：人耳對聲音頻譜的響應是非線性的，Fbank就是一種前端處理算法，以類似於人耳的方式對音頻進行處理，可以提高語音識別的性能。獲得語音信號的fbank特征的一般步驟是：預加重、分幀、加窗、短時傅里葉變換（STFT）、mel濾波、去均值等。對fbank做離散余弦變換（DCT）即可獲得mfcc特征。

MFCC(Mel-frequency cepstral coefficients):梅爾頻率倒譜系數。梅爾頻率是基於人耳聽覺特性提出來的， 它與Hz頻率成非線性對應關系。梅爾頻率倒譜系數(MFCC)則是利用它們之間的這種關系，計算得到的Hz頻譜特征。主要用於語音數據特征提取和降低運算維度。例如：對於一幀有512維(采樣點)數據，經過MFCC后可以提取出最重要的40維(一般而言)數據同時也達到了降維的目的。

下面看看每個步驟的過程及知識點。

二、Fbank處理過程

2.0 讀取數據

數據下載：PCM編碼的WAV格式

# 導包

import numpy as np

from scipy.io import wavfile

from scipy.fftpack import dct

import matplotlib.pyplot as plt

   

def plot_time(sig, fs):

#繪制時域圖

time = np.arange(0, len(sig)) * (1.0 / fs)

plt.figure(figsize=(20, 5))

plt.plot(time, sig)

plt.xlabel('Time(s)')

plt.ylabel('Amplitude')#振幅

plt.grid()

   

def plot_freq(sig, sample_rate, nfft=512):

#繪制頻域圖

freqs = np.linspace(0, sample_rate/2, nfft//2 + 1)

xf = np.fft.rfft(sig, nfft) / nfft

xfp = 20 * np.log10(np.clip(np.abs(xf), 1e-20, 1e100))#強度

plt.figure(figsize=(20, 5))

plt.plot(freqs, xfp)

plt.xlabel('Freq(hz)')

plt.ylabel('dB')#強度

plt.grid()

 

def plot_spectrogram(spec, ylabel = 'ylabel'):

#繪制二維數組

fig = plt.figure(figsize=(20, 5))

heatmap = plt.pcolor(spec)

fig.colorbar(mappable=heatmap)

plt.xlabel('Time(s)')

plt.ylabel(ylabel)

plt.tight_layout()

plt.show()

   

wav_file = 'OSR_us_000_0010_8k.wav'

fs, sig = wavfile.read(wav_file)

#fs是wav文件的采樣率，signal是wav文件的內容，filename是要讀取的音頻文件的路徑

sig = sig[0: int(10 * fs)]

# 保留前10s數據

plot_time(sig, fs)

plot_freq(sig, fs)

2.1 預處理—預加重（Pre-Emphasis）

定義：預加重即對語音的高頻部分進行加重。

預加重的目的：

1. 平衡頻譜，因為高頻通常與較低頻率相比具有較小的幅度，提升高頻部分，使信號的頻譜變得平坦，保持在低頻到高頻的整個頻帶中，能用同樣的噪聲比（SNR）求頻譜。
2. 也是為了消除發生過程中聲帶和嘴唇的效應，來補償語音信號受到發音系統所抑制的高頻部分，也為了突出高頻的共振峰。

預加重的過程：

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

其中濾波器系數（α）的通常為0.95或0.97

代碼形式：emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])

   

題外話：預加重在現代系統中的影響不大，主要是因為除避免在現代FFT實現中不應成為問題的FFT數值問題，大多數預加重濾波器的動機都可以通過均值歸一化來實現（在本文后面討論）。 在現代FFT實現中。

效果展示

pre_emphasis = 0.97

sig = np.append(sig[0], sig[1:] - pre_emphasis * sig[:-1])

plot_time(sig, fs)

plot_freq(sig, fs)

print(sig.shape,fs)

2.2 預處理—分幀（Framing）

目的：

由於語音信號是一個非平穩態過程，不能用處理平穩信號的信號處理技術對其進行分析處理。語音信號是短時平穩信號。因此我們在短時幀上進行傅里葉變換，通過連接相鄰幀來獲得信號頻率輪廓的良好近似。

過程：

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

1. 短時分析將語音流分為一幀來處理，幀長：10～30ms，20ms常見；
2. 幀移：STRIDE，0～1/2幀長，幀與幀之間的平滑長度；

代碼：

def framing(frame_len_s, frame_shift_s, fs, sig):

"""

分幀，主要是計算對應下標

:param frame_len_s: 幀長，s

:param frame_shift_s: 幀移，s

:param fs: 采樣率，hz

:param sig: 信號

:return: 二維list，一個元素為一幀信號

"""

sig_n = len(sig)

frame_len_n, frame_shift_n = int(round(fs * frame_len_s)), int(round(fs * frame_shift_s))

num_frame = int(np.ceil(float(sig_n - frame_len_n) / frame_shift_n) + 1)

pad_num = frame_shift_n * (num_frame - 1) + frame_len_n - sig_n # 待補0的個數

pad_zero = np.zeros(int(pad_num)) # 補0

pad_sig = np.append(sig, pad_zero)

 

# 計算下標

# 每個幀的內部下標

frame_inner_index = np.arange(0, frame_len_n)

# 分幀后的信號每個幀的起始下標

frame_index = np.arange(0, num_frame) * frame_shift_n

 

# 復制每個幀的內部下標，信號有多少幀，就復制多少個，在行方向上進行復制

frame_inner_index_extend = np.tile(frame_inner_index, (num_frame, 1))

# 各幀起始下標擴展維度，便於后續相加

frame_index_extend = np.expand_dims(frame_index, 1)

# 分幀后各幀的下標，二維數組，一個元素為一幀的下標

each_frame_index = frame_inner_index_extend + frame_index_extend

each_frame_index = each_frame_index.astype(np.int, copy=False)

 

frame_sig = pad_sig[each_frame_index]

return frame_sig

   

frame_len_s = 0.025

frame_shift_s = 0.01

frame_sig = framing(frame_len_s, frame_shift_s, fs, sig)

2.3 預處理—加窗（Window）

語音在長范圍內是不停變動的，沒有固定的特性無法做處理，所以將每一幀代入窗函數，窗外的值設定為0，其目的是消除各個幀兩端可能會造成的信號不連續性。

常用的窗函數有方窗、漢明窗和漢寧窗等，根據窗函數的頻域特性，常采用漢明窗。

目的：

1. 使全局更加連續，避免出現吉布斯效應；
2. 加窗時，原本沒有周期性的語音信號呈現出周期函數的部分特征。

過程：

將信號分割成幀后，我們再對每個幀乘以一個窗函數，如Hamming窗口。以增加幀左端和右端的連續性。抵消FFT假設（數據是無限的），並減少頻譜泄漏。漢明窗的形式如下：

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

效果展示：

#frames *= numpy.hamming(frame_length)

# frames *= 0.54 - 0.46 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1)) # 內部實現

# 加窗

window = np.hamming(int(round(frame_len_s * fs)))

plt.figure(figsize=(20, 5))

plt.plot(window)

plt.grid()

plt.xlim(0, 200)

plt.ylim(0, 1)

plt.xlabel('Samples')

plt.ylabel('Amplitude')

frame_sig *= window

plot_time(frame_sig, fs)

plot_freq(frame_sig.reshape(-1,), fs) #需要先變成一維數據

print(frame_sig.reshape(-1,).shape)

   

   

2.4 傅里葉變換FFT（Fourier-Transform）

目的

由於信號在時域上的變換通常很難看出信號的特性，所以通常將它轉換為頻域上的能量分布來觀察，不同的能量分布，就能代表不同語音的特性。所以在乘上漢明窗后，每幀還必須再經過快速傅里葉變換以得到在頻譜上的能量分布。對分幀加窗后的各幀信號進行快速傅里葉變換得到各幀的頻譜。並對語音信號的頻譜取模平方得到語音信號的功率譜。設語音信號的DFT為：

過程：

接我們對分幀加窗后的各幀信號進行做一個N點FFT來計算頻譜，也稱為短時傅立葉變換（STFT），其中N通常為256或512，NFFT＝512；

代碼：

mag_frames = numpy.absolute(numpy.fft.rfft(frames, NFFT)) # fft的幅度(magnitude)

   

功率譜（Power Spectrum）

然后我們使用以下公式計算功率譜（周期圖periodogram），對語音信號的頻譜取模平方（取對數或者去平方，因為頻率不可能為負，負值要舍去）得到語音信號的譜線能量。

短時傅里葉變換返回的就是功率譜的值

代碼和效果：

def stft(frame_sig, nfft=512):

"""

:param frame_sig: 分幀后的信號

:param nfft: fft點數

:return: 返回分幀信號的功率譜

np.fft.fft vs np.fft.rfft

fft 返回 nfft

rfft 返回 nfft // 2 + 1，即rfft僅返回有效部分

"""

frame_spec = np.fft.rfft(frame_sig, nfft)

# 幅度譜

frame_mag = np.abs(frame_spec)

# 功率譜

frame_pow = (frame_mag ** 2) * 1.0 / nfft

return frame_pow

   

nfft = 512

frame_pow = stft(frame_sig, nfft)

   

plt.figure(figsize=(20, 5))

plt.plot(frame_pow[1])

plt.grid()

2.5 梅爾濾波器組（Filter Banks）

計算Mel濾波器組，將功率譜通過一組Mel刻度(通常取40個濾波器，nfilt=40)的三角濾波器(triangular filters)來提取頻帶(frequency bands)。

這個Mel濾波器組就像人類的聽覺感知系統（耳朵），人耳只關注某些特定的頻率分量（人的聽覺對頻率是有選擇性的）。它對不同頻率信號的靈敏度是不同的，換言之，它只讓某些頻率的信號通過，而壓根就直接無視它不想感知的某些頻率信號。但是這些濾波器在頻率坐標軸上卻不是統一分布的，在低頻區域有很多的濾波器，他們分布比較密集，但在高頻區域，濾波器的數目就變得比較少，分布很稀疏。因此Mel刻度的目的是模擬人耳對聲音的非線性感知，在較低的頻率下更具辨別力，在較高的頻率下則不具辨別力。我們可以使用以下公式在赫茲（f）和梅爾（m）之間進行轉換：

定義一個有M個三角濾波器的濾波器組（濾波器的個數和臨界帶的個數相近）：M通常取22-40，26是標准。濾波器組中的每個濾波器都是三角形的，中心頻率為f(m) ，中心頻率處的響應為1，並向0線性減小，直到達到兩個相鄰濾波器的中心頻率，其中響應為0，各f(m)之間的間隔隨着m值的增大而增寬，如圖所示：

這可以通過以下等式建模，三角濾波器的頻率響應定義為:

三角帶通濾波器有兩個主要目的：對頻譜進行平滑化，並消除諧波的作用，突顯原先語音的共振峰。（因此一段語音的音調或音高，是不會呈現在MFCC 參數內，換句話說，以MFCC 為特征的語音辨識系統，並不會受到輸入語音的音調不同而有所影響）此外，還可以降低運算量。

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

代碼和效果：

def mel_filter(frame_pow, fs, n_filter, nfft):

"""

mel 濾波器系數計算

:param frame_pow: 分幀信號功率譜

:param fs: 采樣率 hz

:param n_filter: 濾波器個數

:param nfft: fft點數

:return: 分幀信號功率譜mel濾波后的值的對數值

mel = 2595 * log10(1 + f/700) # 頻率到mel值映射

f = 700 * (10^(m/2595) - 1 # mel值到頻率映射

上述過程本質上是對頻率f對數化

"""

mel_min = 0 # 最低mel值

mel_max = 2595 * np.log10(1 + fs / 2.0 / 700) # 最高mel值，最大信號頻率為 fs/2

mel_points = np.linspace(mel_min, mel_max, n_filter + 2) # n_filter個mel值均勻分布與最低與最高mel值之間

hz_points = 700 * (10 ** (mel_points / 2595.0) - 1) # mel值對應回頻率點，頻率間隔指數化

filter_edge = np.floor(hz_points * (nfft + 1) / fs) # 對應到fft的點數比例上

 

# 求mel濾波器系數

fbank = np.zeros((n_filter, int(nfft / 2 + 1)))

for m in range(1, 1 + n_filter):

f_left = int(filter_edge[m - 1]) # 左邊界點

f_center = int(filter_edge[m]) # 中心點

f_right = int(filter_edge[m + 1]) # 右邊界點

 

for k in range(f_left, f_center):

fbank[m - 1, k] = (k - f_left) / (f_center - f_left)

for k in range(f_center, f_right):

fbank[m - 1, k] = (f_right - k) / (f_right - f_center)

 

# mel 濾波

# [num_frame, nfft/2 + 1] * [nfft/2 + 1, n_filter] = [num_frame, n_filter]

filter_banks = np.dot(frame_pow, fbank.T)

filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)

# 取對數

filter_banks = 20 * np.log10(filter_banks) # dB

 

return filter_banks

   

# mel 濾波

n_filter = 40 # mel濾波器個數

filter_banks = mel_filter(frame_pow, fs, n_filter, nfft)

plot_spectrogram(filter_banks.T, ylabel='Filter Banks')

   

三、MFCC

MFCC實際上就是在Fbank的基礎上增加一個離散余弦變換（DCT）。

3.1 離散余弦變換（DCT）

離散余弦變換經常用於信號處理和圖像處理，用來對信號和圖像進行有損數據壓縮，這是由於離散余弦變換具有很強的"能量集中"特性：大多數的自然信號（包括聲音和圖像）的能量都集中在離散余弦變換后的低頻部分，實際就是對每幀數據在進行一次降維。

上一步驟中計算的濾波器組系數是高度相關的，這在某些機器學習算法中可能是有問題的。因此，我們可以應用離散余弦變換（DCT）對濾波器組系數去相關處理，並產生濾波器組的壓縮表示。通常，對於自動語音識別（ASR），保留所得到的個倒頻譜系數2-13，其余部分被丟棄; 我們這里取 num_ceps = 12。丟棄其他的原因是它們代表了濾波器組系數的快速變化，並且這些精細的細節對自動語音識別（ASR）沒有貢獻。

L階指MFCC系數階數，通常取2-13。這里M是三角濾波器個數。

   

代碼和效果

num_ceps = 12

mfcc = dct(filter_banks, type=2, axis=1, norm='ortho')[:, 1:(num_ceps+1)] # 保持在2-13

plot_spectrogram(mfcc.T, 'MFCC Coefficients')

3.2 動態差分參數的提取

標准的倒譜參數MFCC只反映了語音參數的靜態特性，語音的動態特性可以用這些靜態特征的差分譜來描述。實驗證明：把動、靜態特征結合起來才能有效提高系統的識別性能。

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

MFCC的全部組成其實是由： N維MFCC參數（N/3 MFCC系數+ N/3 一階差分參數+ N/3 二階差分參數）+幀能量（此項可根據需求替換）。

這里的幀能量是指一幀的音量（即能量），也是語音的重要特征，而且非常容易計算。因此，通常再加上一幀的對數能量（定義：一幀內信號的平方和，再取以10為底的對數值，再乘以10）使得每一幀基本的語音特征就多了一維，包括一個對數能量和剩下的倒頻譜參數。另外，解釋下40維是怎么回事，假設離散余弦變換的階數取13，那么經過一階二階差分后就是39維了再加上幀能量總共就是40維，當然這個可以根據實際需要動態調整。

四、fbank與mfcc的標准化

其目的是希望減少訓練集與測試集之間的不匹配。有三種操作：

4.1 去均值 （CMN）

為了均衡頻譜，提升信噪比，可以做一個去均值的操作

# filter_banks去均值

filter_banks -= (np.mean(filter_banks, axis=0) + 1e-8)

plot_spectrogram(filter_banks.T, ylabel='Filter Banks')

#去均值之后的fbank:如下圖所示

#同樣可以對mfcc進行去均值操作。

mfcc -= (np.mean(mfcc, axis=0) + 1e-8)

plot_spectrogram(mfcc.T, 'MFCC Coefficients')

五、fbank與mfcc的比較

fbank特征更多是希望符合聲音信號的本質，擬合人耳的接收特性。DCT是線性變換，會丟失語音信號中原本的一些高度非線性成分。在深度學習之前，受限於算法，mfcc配GMMs-HMMs是ASR的主流做法。當深度學習方法出來之后，由於神經網絡對高度相關的信息不敏感，mfcc不是最優選擇，經過實際驗證，其在神經網絡中的表現也明顯不如fbank。

質疑傅里葉變換是否是必要的操作是明智的。鑒於傅立葉變換本身也是線性運算，忽略它並嘗試直接從時域中的信號中學習可能是有益的。實際上，最近的一些工作已經嘗試過，並且報告了積極的結果。然而，傅立葉變換操作是很難學習的操作，可能會增加實現相同性能所需的數據量和模型復雜性。此外，在進行短時傅里葉變換（stft）時，我們假設信號在這一短時間內是平穩的，因此傅里葉變換的線性不會構成一個關鍵問題。

   

上述代碼：https://github.com/yifanhunter/audio

參考文獻

【1】 https://blog.csdn.net/fengzhonghen/article/details/51722555

【2】MFCC的理解： https://www.cnblogs.com/LXP-Never/p/10918590.html

【3】數據分析圖例（代碼主要來源於此）：https://zhuanlan.zhihu.com/p/130926693

【4】語音： https://zhuanlan.zhihu.com/p/61467187