python求音頻的梅爾倒譜系數

文章目錄

• 寫在前面
• 一、必備的關于音頻的知識
• • 1.關于采樣率
  • 2.關于波形圖，聲譜圖
• 二、音頻的梅爾倒譜系數的含義和求解步驟
• • 1. 梅爾倒譜系數的含義
  • 2. 梅爾倒譜系數的求解步驟概括
• 三、求梅爾倒譜系數詳細過程
• • 1. 預加重、分幀和加窗處理
  • 2. 用周期圖(periodogram)法來進行功率譜(power spectrum)估計；（短時傅裡葉變換）
  • 3.對功率譜用Mel濾波器組進行濾波，計算每個濾波器裡的能量，對每個濾波器的能量取log
  • 3. 進行離散餘弦變換（DCT）變換，保留DCT的第2-13個系數
• 四、代碼總和

寫在前面

唔

這次實驗是對音頻進行處理，求音頻的梅爾倒譜系數，雖然上學期上過了多媒體課，早就學習過梅爾倒譜系數，傅裡葉變換，DCT變換，可這次實驗中遇到的時候腦海裡還是一片空白。不禁想起一段話：“我蹲在菜市場裡，農民伯伯從我身旁走過，隻聽他發自内心的感歎：真是好菜呀！”。希望大家平時學習的時候探究問題的根源，不要像我一樣淺嘗辄止，到頭來隻能感歎自己真是太菜了。

這次參考了好多篇大佬的部落格，這裡列舉一下（代碼基本來源于第一篇，我隻不過加了一些解讀）：

python 計算MFCC詳細步驟

論文筆記:語音情感識别(四）語音特征之聲譜圖，log梅爾譜，MFCC，deltas

梅爾倒譜系數特征(Mel-frequency cepstral coefficients，MFCC)

梅爾頻率倒譜系數（MFCC） 學習筆記

如何決定要使用多少點來做FFT

詳解離散餘弦變換（DCT）

最後寫一句我耳機裡正在放的歌吧“我該怎樣度過人生的旅途，孤獨或者庸俗”

一、必備的關于音頻的知識

想要看懂下面的實驗，這些知識必不可少。

1.關于采樣率

采樣率即在一秒的音頻時間裡進行采樣的次數，采樣率為20k時即是一秒鐘對音頻進行20000次采樣。

根據香農采樣定理，想要完整地還原聲音，采樣率至少要為音頻中最高頻率的2倍。

人耳的可感覺頻率是20-20khz，是以為真實地還原音頻，CD采用了40+kHz的采樣率。

2.關于波形圖，聲譜圖

(1)波形圖：

又稱振幅圖，是音頻的振幅（或能量）這個次元的圖形表達。波形圖的橫坐标一般為時間，縱坐标一般為dB(即分貝)來表示。

Python 在讀音頻時，可以使用 librosa 子產品或者 scipy 子產品，兩種讀取代碼如下：

代碼如下（示例）：

import librosa
import librosa.display
import matplotlib.pyplot as plt

data1 , sampling_rate1 = librosa.load('chew.wav', sr=22050)#采樣率預設值是22050

from scipy.io import wavfile
sample_rate, signal = wavfile.read('chew.wav')
           

librosa子產品得到波形圖的代碼如下：

#波形圖
    plt.figure()
    librosa.display.waveplot(data1*32767, sr=sampling_rate1)
    plt.show()
           

示例波形圖如下：

注意librosa子產品讀音頻時進行了歸一化處理，想要得到真實的振幅可以将信号乘于32767。

(2)聲譜圖：

聲音信号是一維信号，直覺上隻能看到時域資訊，不能看到頻域資訊。通過傅裡葉變換(FT)可以變換到頻域，但是丢失了時域資訊，無法看到時頻關系。聲譜圖是聲音或其他信号的頻率随時間變化時的頻譜（spectrum）的一種直覺表示。

librosa子產品得到聲譜圖的代碼如下：

#聲譜圖
    D = librosa.amplitude_to_db(np.abs(librosa.stft(data1)), ref=np.max)
    librosa.display.specshow(D, y_axis='linear')
    plt.colorbar(format = '%+2.0f dB')
    plt.title('Linear-frequency power spectrogram of aloe')
    plt.show()
           

示例聲譜圖如下：

二、音頻的梅爾倒譜系數的含義和求解步驟

1. 梅爾倒譜系數的含義

感覺實驗表明，人耳對于聲音信号的感覺聚焦于某一特定頻率區域内，而非在整個頻譜包絡中。耳蝸的濾波作用是在對數頻率尺度進行的，在1000Hz以下為線性，在1000Hz以上為對數，這就使得人耳對低頻比高頻更敏感。

心理實體學研究表明，人類對語音信号頻率内容的感覺遵循一種主觀上定義的非線性尺度，該非線性标度可被稱為“Mel”标度。

一般來說，聲音的頻率和人耳所聽到的聲音高低不成正比，而是與音調(人們為了描述聲音高低而定義的概念)成正比，聲音的頻率分布與臨界頻帶分布相一緻。梅爾頻率标度的機關是 Mel，它是為了描繪音調而被定義出來的，它更生動地反映出了頻率和音調的非線性關系。

MFCC是将人耳的聽覺感覺特性和語音産生機制相結合，是以目前大多數語音識别系統廣泛使用這種特征。對頻率軸不均勻劃分是MFCC特征差別于前面普通倒譜特征的最重要的特點，變換到Mel域後，Mel帶通濾波器組的中心頻率是按照Mel刻度均勻排列的。

語音的MFCC特征是基于人耳感覺實驗得到，将人耳當成特定的濾波器，隻考慮某些特定頻率成分。這些濾波器是在頻域上不均勻分布的。更多的濾波器聚集于低頻部分，高頻部分的濾波器較少。

梅爾頻率倒譜系數（Mel Frequency Cepstrum Coefficient, MFCC）考慮到了人類的聽覺特征，先将線性頻譜映射到基于聽覺感覺的Mel非線性頻譜中，然後轉換到倒譜上。

2. 梅爾倒譜系數的求解步驟概括

梅爾倒譜系數(MFCC)特征提取包含以下幾個步驟：

1、對語音信号進行預加重、分幀和加窗處理；

2、用周期圖(periodogram)法來進行功率譜(power spectrum)估計；（短時傅裡葉變換）

3、對功率譜用Mel濾波器組進行濾波，計算每個濾波器裡的能量；（梅爾頻譜）

4、對每個濾波器的能量取log；（log梅爾頻譜）

5、進行離散餘弦變換（DCT）變換；（梅爾倒譜）

6、保留DCT的第2-13個系數，去掉其它。（MFCC特征）

三、求梅爾倒譜系數詳細過程

1. 預加重、分幀和加窗處理

(1)讀入音頻

這裡我們選擇使用python的librosa子產品讀取音頻，注意将振幅乘于32767來抵消歸一化：

signal , sample_rate = librosa.load('chew.wav', sr=22050)#采樣率預設值是22050,注意讀取出來的資料，是做了32767的歸一化
signal = signal * 32767
print(len(signal))

axis_x = numpy.arange(0, signal.size, 1)
plt.plot(axis_x, signal, linewidth=5)
plt.title("Time domain plot")
plt.xlabel("Time", fontsize=14)
plt.ylabel("Amplitude", fontsize=14)
plt.tick_params(axis='both', labelsize=14)
plt.savefig('Time domain plot.png')
plt.show()
           

波形圖如下：

(2)預加重

對信号應用預加重濾波器，以放大高頻。預加重濾波器在幾種方面有用：

（1）平衡頻譜，因為高頻通常比低頻具有較小的幅度；

（2）避免在傅立葉變換操作期間出現數值問題；

（3）還可改善信号 噪聲比（SNR）。

可以使用以下公式中的一階濾波器将預加重濾波器應用于信号x：

y(t) = x(t) - α*x(t-1)，其中濾波器系數（α）的典型值為0.95或0.97

pre_emphasis = 0.97
emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])

axis_x = numpy.arange(0, emphasized_signal.size, 1)
plt.plot(axis_x, emphasized_signal, linewidth=5)
plt.title("Pre-Emphasis")
plt.xlabel("Time", fontsize=14)
plt.ylabel("Amplitude", fontsize=14)
plt.tick_params(axis='both', labelsize=14)
plt.savefig("Pre-Emphasis.png")
plt.show()
           

預加重後波形圖如下：

(3)分幀

經過預加重後，我們需要将信号分成短幀。 此步驟的基本原理是信号中的頻率會随時間變化，是以在大多數情況下，對整個信号進行傅立葉變換是沒有意義的，因為我們會随時間丢失信号的頻率輪廓。

為避免這種情況，我們可以假設信号的頻率在很短的時間内是固定的。 是以，通過在此短幀上進行傅立葉變換，可以通過串聯相鄰幀來獲得信号頻率輪廓較好的近似。

語音進行中的典型幀大小為20毫秒至40毫秒，連續幀之間有50％（+/- 10％）重疊。常見的設定是幀大小為25毫秒，frame_size = 0.025和10毫秒跨度（重疊15毫秒）

先計算一些長度，注意這裡的長度是指包含采樣點的數量。

frame_length是每一幀的長度，frame_step是幀的跨度，signal_length是音頻信号長度。

frame_stride = 0.01
frame_size = 0.025
frame_length, frame_step = frame_size * sample_rate, frame_stride * sample_rate  # Convert from seconds to samples
signal_length = len(emphasized_signal)
frame_length = int(round(frame_length))
frame_step = int(round(frame_step))
           

求幀的數量，這裡用到了numpy.ceil函數，即向上取整函數。

#numpy.ceil 向上取整函數
num_frames = int(numpy.ceil
				 (float(numpy.abs(signal_length - frame_length + frame_step)) / frame_step))  # Make sure that we have at least 1 frame
           

由于信号的長度可能不是幀長的整數倍，我們需要對原信号進行填充，将其填充至幀長的整數倍，下面代碼中的z即是要填充的0矩陣：

pad_signal_length = (num_frames - 1) * frame_step + frame_length
# 填充信号以確定所有幀具有相同數量的樣本，而不會截斷原始信号中的任何樣本
z = numpy.zeros((pad_signal_length - signal_length))
pad_signal = numpy.append(emphasized_signal, z)
           

填充前和填充後音頻信号的長度如下：

下面我們将音頻信号轉化為二維矩陣，每一行即是一個音頻幀的内容，frames即是分幀後的二維矩陣，在本次實驗中，幀矩陣的shape是325*551

#前一個 num_frames * frame_length 後一個 num_frames * frame_length,index是個二維矩陣，每一行是一個frame的值，比如第一行是0-550，第二行就是220-770
indices = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1)) + numpy.tile(
	numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T
print(indices.shape)
#将pad_signal轉成num_frames * frame_length的格式，indices中的值是pad_signal的index
frames = pad_signal[indices.astype(numpy.int32, copy=False)]
           

(4)加窗

将信号切成幀後，我們對每個幀應用諸如漢明窗之類的視窗函數。 Hamming視窗具有以下形式：w[n]=0.54-0.46cos(2pin/(N-1))

其中0<=n<=N-1, N是窗長，有很多原因需要将窗函數應用于這些幀，特别是要抵消FFT無限計算并減少頻譜洩漏。

使用numpy的窗函數對音頻信号加窗：

2. 用周期圖(periodogram)法來進行功率譜(power spectrum)估計；（短時傅裡葉變換）

現在，我們可以在每個幀上執行N點FFT來計算頻譜，這也稱為短時傅立葉變換（STFT），

其中N通常為256或512，NFFT = 512； 然後使用以下公式計算功率譜（周期圖）：P=|FFT(xi)|^2/N，其中，xi是信号x的第i幀。

NFFT = 512
mag_frames = numpy.absolute(numpy.fft.rfft(frames, NFFT))  # Magnitude of the FFT
#print(mag_frames.shape)
pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2))  # Power Spectrum,**2是平方
           

在本次實驗中，功率譜的shape是 325*257。

功率譜示例如下，它的橫坐标是幀下标，縱坐标是不同頻率，圖中顔色越深(比如紅色)，對應頻率的能量越大：

3.對功率譜用Mel濾波器組進行濾波，計算每個濾波器裡的能量，對每個濾波器的能量取log

聲譜圖往往是很大的一張圖，且依舊包含了大量無用的資訊，是以我們需要通過梅爾标度濾波器組（mel-scale filter banks）将其變為梅爾頻譜

(1)梅爾尺度（Mel Scale）是建立從人類的聽覺感覺的頻率——Pitch到聲音實際頻率直接的映射。頻率的機關是赫茲（Hz），人耳能聽到的頻率範圍是20-20000Hz，但人耳對Hz這種标度機關并不是線性感覺關系，例如，若把音調頻率從1000Hz提高到2000Hz，我們的耳朵隻能覺察到頻率似乎提高了一些而不是一倍。但是通過把頻率轉換成美爾尺度，我們的特征就能夠更好的比對人類的聽覺感覺效果。從頻率到梅爾頻率的轉換公式如下：

或者

我們可以觀察一下轉換後的映射圖，可以發現人耳對于低頻聲音的分辨率要高于高頻的聲音，因為赫茲到梅爾是log的關系，是以當頻率較小時，mel随Hz變化較快；當頻率很大時，mel的上升很緩慢，曲線的斜率很小。這說明了人耳對低頻音調的感覺較靈敏，在高頻時人耳是很遲鈍的，梅爾标度濾波器組啟發于此。

(2)梅爾濾波器

為了模拟人耳對聲音的感覺，人們發明的梅爾濾波器組。一組大約20-40(通常26)個三角濾波器組，它會對上一步得到的周期圖的功率譜估計進行濾波。而且區間的頻率越高，濾波器就越寬(但是如果把它變換到美爾尺度則是一樣寬的)。為了計算友善，我們通常把26個濾波器用一個矩陣來表示，這個矩陣有26行，列數就是傅裡葉變換的點數。

計算過程如下圖所示，最後我們會保留這26個濾波器的能量。圖(a)是26個濾波器；圖(b)是濾波後的信号；圖©是其中的第8個濾波器，它隻讓某一頻率範圍的信号通過；圖(d)通過它的信号的能量；圖(e)是第20個濾波器；圖(f)是通過它的信号的能量。

(3)濾波

選取最低頻率為0，最高頻率為采樣率的一半，設定26個區間，因為一個濾波器實際山要占據三個刻度，是以要劃分出28和區間。

使用公式進行頻率和梅爾刻度的轉換。

nfilt = 26
low_freq_mel = 0
high_freq_mel = (2595 * numpy.log10(1 + (sample_rate / 2) / 700))  # Convert Hz to Mel
mel_points = numpy.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # Equally spaced in Mel scale，劃分出nfilt+2個區間
#print(mel_points)
hz_points = (700 * (10 ** (mel_points / 2595) - 1))  # Convert Mel to Hz
#print(hz_points)
bin = numpy.floor((NFFT + 1) * hz_points / sample_rate)
           

bin儲存的是刻度對應的傅裡葉變換點數。

fbank的shape是26*257，用來存儲每個濾波器的值。

fbank = numpy.zeros((nfilt, int(numpy.floor(NFFT / 2 + 1))))
print(fbank.shape)
for m in range(1, nfilt + 1):
	f_m_minus = int(bin[m - 1])  # left
	f_m = int(bin[m])  # center
	f_m_plus = int(bin[m + 1])  # right

	for k in range(f_m_minus, f_m):
		fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1])
	for k in range(f_m, f_m_plus):
		fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m])
           

将功率譜與濾波器做點積，這時filter_bank的shape是325*26，達成了降維的目标。

再将filter_bank中的0值改為最小負數，防止運算出現問題，再對每個濾波器的能量取log即得到log梅爾頻譜。

filter_banks = numpy.dot(pow_frames, fbank.T)
filter_banks = numpy.where(filter_banks == 0, numpy.finfo(float).eps, filter_banks)  # Numerical Stability
filter_banks = 20 * numpy.log10(filter_banks)  # dB
           

3. 進行離散餘弦變換（DCT）變換，保留DCT的第2-13個系數

進行離散餘弦變換（DCT）變換，保留DCT的第2-13個系數

事實證明，在上一步中計算出的濾波器組系數是高度相關的，這在某些機器學習算法中可能會出現問題。是以，我們可以應用離散餘弦變換（DCT）去相關濾波器組系數，并産生濾波器組的壓縮表示。

通常，對于自動語音識别（ASR），結果倒譜系數2-13将保留，其餘的将被丢棄； num_ceps =12。丢棄其他系數的原因是它們代表濾波器組系數的快速變化，而這些細微的細節對自動語音識别（ASR）毫無幫助。

這裡在進行DCT變化時用到了scipy子產品的scipy.fftpack.dct函數。

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

"""
可以将正弦提升器1應用于MFCC，去加重過高的MFCCs，這被可以改善嘈雜信号中的語音識别。
"""
cep_lifter = 22
(nframes, ncoeff) = mfcc.shape
n = numpy.arange(ncoeff)
lift = 1 + (cep_lifter / 2) * numpy.sin(numpy.pi * n / cep_lifter)
mfcc *= lift  # *
           

四、代碼總和

import numpy
import librosa
from scipy.io import wavfile
from scipy.fftpack import dct
import matplotlib.pyplot as plt

signal , sample_rate = librosa.load('chew.wav', sr=22050)#采樣率預設值是22050,注意讀取出來的資料，是做了32767的歸一化
signal = signal * 32767
print(len(signal))

axis_x = numpy.arange(0, signal.size, 1)
plt.plot(axis_x, signal, linewidth=5)
plt.title("Time domain plot")
plt.xlabel("Time", fontsize=14)
plt.ylabel("Amplitude", fontsize=14)
plt.tick_params(axis='both', labelsize=14)
plt.savefig('Time domain plot.png')
plt.show()

"""
Pre-Emphasis 預加重
"""
pre_emphasis = 0.97
emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])

axis_x = numpy.arange(0, emphasized_signal.size, 1)
plt.plot(axis_x, emphasized_signal, linewidth=5)
plt.title("Pre-Emphasis")
plt.xlabel("Time", fontsize=14)
plt.ylabel("Amplitude", fontsize=14)
plt.tick_params(axis='both', labelsize=14)
plt.savefig("Pre-Emphasis.png")
plt.show()

frame_stride = 0.01
frame_size = 0.025
frame_length, frame_step = frame_size * sample_rate, frame_stride * sample_rate  # Convert from seconds to samples
signal_length = len(emphasized_signal)
frame_length = int(round(frame_length))
frame_step = int(round(frame_step))
'''num_frames = int(numpy.ceil
				 (float(numpy.abs(signal_length - frame_length)) / frame_step))  # Make sure that we have at least 1 frame
#numpy.ceil 向上取整函數
pad_signal_length = num_frames * frame_step + frame_length'''
#numpy.ceil 向上取整函數
num_frames = int(numpy.ceil
				 (float(numpy.abs(signal_length - frame_length + frame_step)) / frame_step))  # Make sure that we have at least 1 frame
pad_signal_length = (num_frames - 1) * frame_step + frame_length
# 填充信号以確定所有幀具有相同數量的樣本，而不會截斷原始信号中的任何樣本
z = numpy.zeros((pad_signal_length - signal_length))
pad_signal = numpy.append(emphasized_signal, z)
print(emphasized_signal.shape)
print(pad_signal.shape)
#前一個 num_frames * frame_length 後一個 num_frames * frame_length,index是個二維矩陣，每一行是一個frame的值，比如第一行是0-550，第二行就是220-770
indices = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1)) + numpy.tile(
	numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T
print(indices.shape)
#将pad_signal轉成num_frames * frame_length的格式，indices中的值是pad_signal的index
frames = pad_signal[indices.astype(numpy.int32, copy=False)]


frames *= numpy.hamming(frame_length)
# frames *= 0.54 - 0.46 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1))  # Explicit Implementation **

# 傅立葉變換和功率譜

NFFT = 512
mag_frames = numpy.absolute(numpy.fft.rfft(frames, NFFT))  # Magnitude of the FFT
#print(mag_frames.shape)
pow_frames = ((1.0 / NFFT) * ((mag_frames) ** 2))  # Power Spectrum,**2是平方
#print(pow_frames.shape)

# 濾波器組 Filter Banks

nfilt = 26
low_freq_mel = 0
high_freq_mel = (2595 * numpy.log10(1 + (sample_rate / 2) / 700))  # Convert Hz to Mel
mel_points = numpy.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # Equally spaced in Mel scale，劃分出nfilt+2個區間
#print(mel_points)
hz_points = (700 * (10 ** (mel_points / 2595) - 1))  # Convert Mel to Hz
#print(hz_points)
bin = numpy.floor((NFFT + 1) * hz_points / sample_rate)

fbank = numpy.zeros((nfilt, int(numpy.floor(NFFT / 2 + 1))))
print(fbank.shape)
for m in range(1, nfilt + 1):
	f_m_minus = int(bin[m - 1])  # left
	f_m = int(bin[m])  # center
	f_m_plus = int(bin[m + 1])  # right

	for k in range(f_m_minus, f_m):
		fbank[m - 1, k] = (k - bin[m - 1]) / (bin[m] - bin[m - 1])
	for k in range(f_m, f_m_plus):
		fbank[m - 1, k] = (bin[m + 1] - k) / (bin[m + 1] - bin[m])
filter_banks = numpy.dot(pow_frames, fbank.T)
filter_banks = numpy.where(filter_banks == 0, numpy.finfo(float).eps, filter_banks)  # Numerical Stability
filter_banks = 20 * numpy.log10(filter_banks)  # dB

plt.title("filter_banks")
plt.imshow(numpy.flipud(filter_banks.T), cmap=plt.cm.jet, aspect=0.1,
		   extent=[0, filter_banks.shape[1], 0, filter_banks.shape[0]])  # 畫熱力圖
plt.xlabel("Frames", fontsize=14)
plt.ylabel("Dimension", fontsize=14)
plt.tick_params(axis='both', labelsize=14)
plt.savefig('filter_banks.png')
plt.show()



# 梅爾倒譜Mel-frequency Cepstral Coefficients (MFCCs)

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


cep_lifter = 22
(nframes, ncoeff) = mfcc.shape
n = numpy.arange(ncoeff)
lift = 1 + (cep_lifter / 2) * numpy.sin(numpy.pi * n / cep_lifter)
mfcc *= lift  # *
#print(mfcc.shape)

plt.title("mfcc")
plt.imshow(numpy.flipud(mfcc.T), cmap=plt.cm.jet, aspect=0.05, extent=[0, mfcc.shape[1], 0, mfcc.shape[0]])  # 畫熱力圖
plt.xlabel("Frames", fontsize=14)
plt.ylabel("Dimension", fontsize=14)
plt.tick_params(axis='both', labelsize=14)
plt.savefig('mfcc.png')
plt.show()