一步一步講解和實現ASR中常用的語音特征——FBank和MFCC的提取，包括算法原理、代碼和可視化等。
完整Jupyter Notebook鏈接：https://github.com/Magic-Bubble/SpeechProcessForMachineLearning/blob/master/speech_process.ipynb
文章目錄
- 語音信號的產生
- 準備工作
- 1. 導包
- 2. 繪圖工具
- 3. 數據準備
- 預加重（Pre-Emphasis）
- 分幀（Framing）
- 加窗（Window）
- 快速傅里葉變換（FFT）
- FBank特征（Filter Banks）
- MFCC特征（Mel-frequency Cepstral Coefficients）
- FBank與MFCC比較
- 其他特征
- 標準化
- 總結
- 傳送門
語音信號的產生
語音通常是指人說話的聲音。從生物學的角度來看，是氣流通過聲帶、咽喉、口腔、鼻腔等發出聲音；從信號的角度來看，不同位置的震動頻率不一樣，最后的信號是由基頻和一些諧波構成。
之后被設備接收后（比如麥克風），會通過A/D轉換，將模擬信號轉換為數字信號，一般會有采樣、量化和編碼三個步驟，采樣率要遵循奈奎斯特采樣定律： f s > = 2 f fs >= 2f f s > = 2 f ，比如電話語音的頻率一般在300Hz~3400Hz，所以采用8kHz的采樣率足矣。
下面采用一個30s左右的16比特PCM編碼后的語音wav為例。
準備工作
1. 導包
import
numpy
as
np
from
scipy
.
io
import
wavfile
from
scipy
.
fftpack
import
dct
import
warnings
warnings
.
filterwarnings
(
'ignore'
)
import
matplotlib
.
pyplot
as
plt
%
matplotlib inline
2. 繪圖工具
# 繪制時域圖
def
plot_time
(
signal
,
sample_rate
)
:
time
=
np
.
arange
(
0
,
len
(
signal
)
)
*
(
1.0
/
sample_rate
)
plt
.
figure
(
figsize
=
(
20
,
5
)
)
plt
.
plot
(
time
,
signal
)
plt
.
xlabel
(
'Time(s)'
)
plt
.
ylabel
(
'Amplitude'
)
plt
.
grid
(
)
# 繪制頻域圖
def
plot_freq
(
signal
,
sample_rate
,
fft_size
=
512
)
:
xf
=
np
.
fft
.
rfft
(
signal
,
fft_size
)
/
fft_size
freqs
=
np
.
linspace
(
0
,
sample_rate
/
2
,
fft_size
/
2
+
1
)
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
,
note
)
:
fig
=
plt
.
figure
(
figsize
=
(
20
,
5
)
)
heatmap
=
plt
.
pcolor
(
spec
)
fig
.
colorbar
(
mappable
=
heatmap
)
plt
.
xlabel
(
'Time(s)'
)
plt
.
ylabel
(
note
)
plt
.
tight_layout
(
)
3. 數據準備
sample_rate
,
signal
=
wavfile
.
read
(
'./resources/OSR_us_000_0010_8k.wav'
)
signal
=
signal
[
0
:
int
(
3.5
*
sample_rate
)
]
# Keep the first 3.5 seconds
print
(
'sample rate:'
,
sample_rate
,
', frame length:'
,
len
(
signal
)
)
sample rate: 8000 , frame length: 28000
plot_time
(
signal
,
sample_rate
)
plot_freq
(
signal
,
sample_rate
)
預加重（Pre-Emphasis）
預加重一般是數字語音信號處理的第一步。語音信號往往會有頻譜傾斜（Spectral Tilt）現象，即高頻部分的幅度會比低頻部分的小，預加重在這里就是起到一個平衡頻譜的作用，增大高頻部分的幅度。它使用如下的一階濾波器來實現：
y ( t ) = x ( t ) ? α x ( t ? 1 ) , ???? 0.95 < α < 0.99 y(t) = x(t) - \alpha x(t-1), \ \ \ \ 0.95 < \alpha < 0.99 y ( t ) = x ( t ) ? α x ( t ? 1 ) , ? ? ? ? 0 . 9 5 < α < 0 . 9 9
筆者對這個公式的理解是：信號頻率的高低主要是由信號電平變化的速度所決定，對信號做一階差分時，高頻部分（變化快的地方）差分值大，低頻部分（變化慢的地方）差分值小，達到平衡頻譜的作用。
pre_emphasis
=
0.97
emphasized_signal
=
np
.
append
(
signal
[
0
]
,
signal
[
1
:
]
-
pre_emphasis
*
signal
[
:
-
1
]
)
plot_time
(
emphasized_signal
,
sample_rate
)
plot_freq
(
emphasized_signal
,
sample_rate
)
從下面這個圖來看，確實起到了平衡頻譜的作用。
分幀（Framing）
在預加重之后，需要將信號分成短時幀。做這一步的原因是：信號中的頻率會隨時間變化（不穩定的），一些信號處理算法（比如傅里葉變換）通常希望信號是穩定，也就是說對整個信號進行處理是沒有意義的，因為信號的頻率輪廓會隨著時間的推移而丟失。為了避免這種情況，需要對信號進行分幀處理，認為每一幀之內的信號是短時不變的。一般設置幀長取20ms~40ms，相鄰幀之間50%（+/-10%）的覆蓋。對于ASR而言，通常取幀長為25ms，覆蓋為10ms。
frame_size
,
frame_stride
=
0.025
,
0.01
frame_length
,
frame_step
=
int
(
round
(
frame_size
*
sample_rate
)
)
,
int
(
round
(
frame_stride
*
sample_rate
)
)
signal_length
=
len
(
emphasized_signal
)
num_frames
=
int
(
np
.
ceil
(
np
.
abs
(
signal_length
-
frame_length
)
/
frame_step
)
)
+
1
pad_signal_length
=
(
num_frames
-
1
)
*
frame_step
+
frame_length
z
=
np
.
zeros
(
(
pad_signal_length
-
signal_length
)
)
pad_signal
=
np
.
append
(
emphasized_signal
,
z
)
indices
=
np
.
arange
(
0
,
frame_length
)
.
reshape
(
1
,
-
1
)
+
np
.
arange
(
0
,
num_frames
*
frame_step
,
frame_step
)
.
reshape
(
-
1
,
1
)
frames
=
pad_signal
[
indices
]
print
(
frames
.
shape
)
(349, 200)
加窗（Window）
在分幀之后，通常需要對每幀的信號進行加窗處理。目的是讓幀兩端平滑地衰減，這樣可以降低后續傅里葉變換后旁瓣的強度，取得更高質量的頻譜。常用的窗有：矩形窗、漢明（Hamming）窗、漢寧窗（Hanning），以漢明窗為例，其窗函數為：
w ( n ) = 0.54 ? 0.46 c o s ( 2 π n N ? 1 ) w(n) = 0.54 - 0.46 cos(\frac{2\pi n}{N-1}) w ( n ) = 0 . 5 4 ? 0 . 4 6 c o s ( N ? 1 2 π n ? )
這里的 0 < = n < = N ? 1 0<=n<=N-1 0 < = n < = N ? 1 ， N N N 是窗的寬度。
hamming
=
np
.
hamming
(
frame_length
)
# hamming = 0.54 - 0.46 * np.cos(2 * np.pi * np.arange(0, frame_length) / (frame_length - 1))
plt
.
figure
(
figsize
=
(
20
,
5
)
)
plt
.
plot
(
hamming
)
plt
.
grid
(
)
plt
.
xlim
(
0
,
200
)
plt
.
ylim
(
0
,
1
)
plt
.
xlabel
(
'Samples'
)
plt
.
ylabel
(
'Amplitude'
)
frames
*=
hamming
plot_time
(
frames
[
1
]
,
sample_rate
)
plot_freq
(
frames
[
1
]
,
sample_rate
)
快速傅里葉變換（FFT）
對于每一幀的加窗信號，進行N點FFT變換，也稱短時傅里葉變換（STFT），N通常取256或512，然后用如下的公式計算能量譜：
P = ∣ F F T ( x i ) ∣ 2 N P = \frac{|FFT(x_i)|^2}{N} P = N ∣ F F T ( x i ? ) ∣ 2 ?
NFFT
=
512
mag_frames
=
np
.
absolute
(
np
.
fft
.
rfft
(
frames
,
NFFT
)
)
pow_frames
=
(
(
1.0
/
NFFT
)
*
(
mag_frames
**
2
)
)
print
(
pow_frames
.
shape
)
(349, 257)
plt
.
figure
(
figsize
=
(
20
,
5
)
)
plt
.
plot
(
pow_frames
[
1
]
)
plt
.
grid
(
)
FBank特征（Filter Banks）
經過上面的步驟之后，在能量譜上應用Mel濾波器組，就能提取到FBank特征。
在介紹Mel濾波器組之前，先介紹一下Mel刻度，這是一個能模擬人耳接收聲音規律的刻度，人耳在接收聲音時呈現非線性狀態，對高頻的更不敏感，因此Mel刻度在低頻區分辨度較高，在高頻區分辨度較低，與頻率之間的換算關系為：
m = 2595 l o g 10 ( 1 + f 700 ) m = 2595 log_{10} (1 + \frac{f}{700}) m = 2 5 9 5 l o g 1 0 ? ( 1 + 7 0 0 f ? )
f = 700 ( 1 0 m / 2595 ? 1 ) f = 700(10^{m/2595} - 1) f = 7 0 0 ( 1 0 m / 2 5 9 5 ? 1 )
Mel濾波器組就是一系列的三角形濾波器，通常有40個或80個，在中心頻率點響應值為1，在兩邊的濾波器中心點衰減到0，如下圖：
具體公式可以寫為：
最后在能量譜上應用Mel濾波器組，其公式為：
Y t ( m ) = ∑ k = 1 N H m ( k ) ∣ X t ( k ) ∣ 2 Y_t(m) = \sum_{k=1}^{N} H_m(k)|X_t(k)|^2 Y t ? ( m ) = k = 1 ∑ N ? H m ? ( k ) ∣ X t ? ( k ) ∣ 2
其中，k表示FFT變換后的編號，m表示mel濾波器的編號。
low_freq_mel
=
0
high_freq_mel
=
2595
*
np
.
log10
(
1
+
(
sample_rate
/
2
)
/
700
)
print
(
low_freq_mel
,
high_freq_mel
)
0 2146.06452750619
nfilt
=
40
mel_points
=
np
.
linspace
(
low_freq_mel
,
high_freq_mel
,
nfilt
+
2
)
# 所有的mel中心點，為了方便后面計算mel濾波器組，左右兩邊各補一個中心點
hz_points
=
700
*
(
10
**
(
mel_points
/
2595
)
-
1
)
fbank
=
np
.
zeros
(
(
nfilt
,
int
(
NFFT
/
2
+
1
)
)
)
# 各個mel濾波器在能量譜對應點的取值
bin
=
(
hz_points
/
(
sample_rate
/
2
)
)
*
(
NFFT
/
2
)
# 各個mel濾波器中心點對應FFT的區域編碼，找到有值的位置
for
i
in
range
(
1
,
nfilt
+
1
)
:
left
=
int
(
bin
[
i
-
1
]
)
center
=
int
(
bin
[
i
]
)
right
=
int
(
bin
[
i
+
1
]
)
for
j
in
range
(
left
,
center
)
:
fbank
[
i
-
1
,
j
+
1
]
=
(
j
+
1
-
bin
[
i
-
1
]
)
/
(
bin
[
i
]
-
bin
[
i
-
1
]
)
for
j
in
range
(
center
,
right
)
:
fbank
[
i
-
1
,
j
+
1
]
=
(
bin
[
i
+
1
]
-
(
j
+
1
)
)
/
(
bin
[
i
+
1
]
-
bin
[
i
]
)
print
(
fbank
)
[[0. 0.46952675 0.93905351 … 0. 0. 0. ]
[0. 0. 0. … 0. 0. 0. ]
[0. 0. 0. … 0. 0. 0. ]
…
[0. 0. 0. … 0. 0. 0. ]
[0. 0. 0. … 0. 0. 0. ]
[0. 0. 0. … 0.14650797 0.07325398 0. ]]
filter_banks
=
np
.
dot
(
pow_frames
,
fbank
.
T
)
filter_banks
=
np
.
where
(
filter_banks
==
0
,
np
.
finfo
(
float
)
.
eps
,
filter_banks
)
filter_banks
=
20
*
np
.
log10
(
filter_banks
)
# dB
print
(
filter_banks
.
shape
)
(349, 40)
plot_spectrogram
(
filter_banks
.
T
,
'Filter Banks'
)
PS：“log mel-filter bank outputs”和“FBANK features”說的是同一個東西。
MFCC特征（Mel-frequency Cepstral Coefficients）
前面提取到的FBank特征，往往是高度相關的。因此可以繼續用DCT變換，將這些相關的濾波器組系數進行壓縮。對于ASR來說，通常取2~13維，扔掉的信息里面包含濾波器組系數快速變化部分，這些細節信息在ASR任務上可能沒有幫助。
DCT變換其實是逆傅里葉變換的等價替代：
y t ( n ) = ∑ m = 0 M ? 1 l o g ( Y t ( m ) ) c o s ( n ( m + 0.5 ) π M ) , ????? n = 0 , . . . , J y_t(n) = \sum_{m=0}^{M-1}log(Y_t(m))cos(n(m + 0.5) \frac{\pi}{M}), \ \ \ \ \ n = 0, ..., J y t ? ( n ) = m = 0 ∑ M ? 1 ? l o g ( Y t ? ( m ) ) c o s ( n ( m + 0 . 5 ) M π ? ) , ? ? ? ? ? n = 0 , . . . , J
所以MFCC名字里面有倒譜（Cepstral）。
num_ceps
=
12
mfcc
=
dct
(
filter_banks
,
type
=
2
,
axis
=
1
,
norm
=
'ortho'
)
[
:
,
1
:
(
num_ceps
+
1
)
]
print
(
mfcc
.
shape
)
(349, 12)
plot_spectrogram
(
mfcc
.
T
,
'MFCC Coefficients'
)
一般對于ASR來說，對MFCC進行一個正弦提升（sinusoidal liftering）操作，可以提升在噪聲信號中最后的識別率：
M F C C i ′ = w i M F C C i MFCC'_i = w_i MFCC_i M F C C i ′ ? = w i ? M F C C i ?
w i = D 2 s i n ( π ? i D ) w_i = \frac{D}{2} sin(\frac{\pi * i}{D}) w i ? = 2 D ? s i n ( D π ? i ? )
從公式看，猜測原因可能是對頻譜做一個平滑，如果 D D D 取值較大時，會加重高頻部分，使得噪聲被弱化？
cep_lifter
=
23
(
nframes
,
ncoeff
)
=
mfcc
.
shape
n
=
np
.
arange
(
ncoeff
)
lift
=
1
+
(
cep_lifter
/
2
)
*
np
.
sin
(
np
.
pi
*
n
/
cep_lifter
)
mfcc
*=
lift
plot_spectrogram
(
mfcc
.
T
,
'MFCC Coefficients'
)
FBank與MFCC比較
FBank特征的提取更多的是希望符合聲音信號的本質，擬合人耳接收的特性。而MFCC特征多的那一步則是受限于一些機器學習算法。很早之前MFCC特征和GMMs-HMMs方法結合是ASR的主流。而當一些深度學習方法出來之后，MFCC則不一定是最優選擇，因為神經網絡對高度相關的信息不敏感，而且DCT變換是線性的，會丟失語音信號中原本的一些非線性成分。
還有一些說法是在質疑傅里葉變換的使用，因為傅里葉變換也是線性的。因此也有很多方法，設計模型直接從原始的音頻信號中提取特征，但這種方法會增加模型的復雜度，而且本身傅里葉變換不太容易擬合。同時傅里葉變換是在短時上應用的，可以建設信號在這個短的時間內是靜止的，因此傅里葉變換的線性也不會造成很嚴重的問題。
結論就是：在模型對高相關的信號不敏感時（比如神經網絡），可以用FBank特征；在模型對高相關的信號敏感時（比如GMMs-HMMs），需要用MFCC特征。從目前的趨勢來看，因為神經網絡的逐步發展，FBank特征越來越流行。
其他特征
- PLP（Perceptual Linear Prediction）
另外一種特征，與MFCC相比有一些優勢，具體提取方式見下圖：
- 動態特征
加入表現幀之間變化的特征，用如下公式：
d ( t ) = c ( t + 1 ) ? c ( t ? 1 ) 2 d(t) = \frac{c(t+1) - c(t-1)}{2} d ( t ) = 2 c ( t + 1 ) ? c ( t ? 1 ) ?
一般在ASR中使用的特征（用于GMM相關的系統），是39維的；包括(12維MFCC+1維能量) + delta + delta^2
具體提取過程見下圖：
標準化
其目的是希望減少訓練集與測試集之間的不匹配。有三種操作：
- 去均值 （CMN）
為了均衡頻譜，提升信噪比，可以做一個去均值的操作
filter_banks
-=
(
np
.
mean
(
filter_banks
,
axis
=
0
)
+
1e
-
8
)
plot_spectrogram
(
filter_banks
.
T
,
'Filter Banks'
)
mfcc
-=
(
np
.
mean
(
mfcc
,
axis
=
0
)
+
1e
-
8
)
plot_spectrogram
(
mfcc
.
T
,
'MFCC Coefficients'
)
- 方差歸一（CVN）
除以標準差，從而使得方差為1
- 標準化（CMVN）
y t ( j ) = y t ( j ) ? μ ( y ( j ) ) σ ( y ( j ) ) y_t(j) = \frac{y_t(j) - \mu (y(j))}{\sigma (y(j))} y t ? ( j ) = σ ( y ( j ) ) y t ? ( j ) ? μ ( y ( j ) ) ?
PS：這些操作，還可以針對speaker/channel做；在實時情景下，可以計算moving average。
總結
最后引用文末slide里面的一個總結：
傳送門
Speech Processing for Machine Learning: Filter banks, Mel-Frequency Cepstral Coefficients (MFCCs) and What’s In-Between 一個很優質，講的很清楚的英文博客
Speech Signal Analysis 英國愛丁堡大學一門ASR課程的講義
python_speech_features 一個很成熟的python提取這些特征的包
ASR中常用的語音特征之FBank和MFCC（原理 + Python實現） 個人博客
更多文章、技術交流、商務合作、聯系博主
微信掃碼或搜索：z360901061
微信掃一掃加我為好友
QQ號聯系： 360901061
您的支持是博主寫作最大的動力，如果您喜歡我的文章，感覺我的文章對您有幫助，請用微信掃描下面二維碼支持博主2元、5元、10元、20元等您想捐的金額吧，狠狠點擊下面給點支持吧，站長非常感激您！手機微信長按不能支付解決辦法：請將微信支付二維碼保存到相冊，切換到微信，然后點擊微信右上角掃一掃功能，選擇支付二維碼完成支付。
【本文對您有幫助就好】元
