最近在閱讀語音方向的論文，其中有個被提及很多的語音信號特征MFCC(Mel-Frequency Cepstral Coefficients)，找到了基于python的語音庫librosa(version=0.7.1)和python_speech_features(version=0.6)，下文對這兩個庫計算MFCC的流程細節稍作梳理。

LibROSA - librosa 0.7.1 documentation?librosa.github.iopython_speech_features?pypi.org

一、librosa

1.源語音信號，shape = wav.length

wav, sample_rate = librosa.load(path, sr=22050, mono=True, offset=0.0, duration=None,dtype=np.float32, res_type='kaiser_best') #加載語音文件得到原始數據

2.填充及分幀(無預加重處理)，分幀后所有幀的shape = n_ftt * n_frames

y = numpy.pad(array, pad_width, mode, **kwargs) # 原函數及參數 y = numpy.pad(wav, int(n_fft // 2), mode='reflect') # 默認，n_fft 為傅里葉變換維度 y = numpy.pad(wav, (0,pad_need), mode='constant') # 自定義

librosa調用numpy對源語音數據進行填充，默認模式是'reflect'進行鏡像填充，舉個例子：

對序列[1,2,3,4,5]進行左填充2個、右填充3個，左邊以1作對稱軸填充[3,2]，右邊以5作對稱軸填充[4,3,2]，最后結果為[3,2,1,2,3,4,5,4,3,2]。

自定義模式是按照python_speech_features的方式設置參數的，pad_need是經過計算得到的需要填充的數據數量，在下文說明。

y_frames = util.frame(y, frame_length=n_fft, hop_length=hop_length) # hop_length為幀移，librosa中默認取窗長的四分之一

一般來說 幀長 = 窗長 = 傅里葉變換維度，否則要進行填充或者截斷處理。通過閱讀源碼發現，librosa調用的分幀方式和python_speech_features一樣，這和librosa.feature.melspectrogram()參數里的center描述有所沖突，我把兩者對同一語音文件的分幀結果輸出到excel文件進行了對比，結果是一樣的，當然分幀之前的填充方式也修改成了一致。

3.對所有幀進行加窗，shape = n_ftt * n_frames。librosa中window.shape = n_ftt * 1

fft_window = librosa.filters.get_window(window, Nx, fftbins=True) # 原函數及參數 fft_window = get_window('hann', win_length, fftbins=True) # 窗長一般等于傅里葉變換維度，短則填充長則截斷

librosa加的窗函數調用的scipy，如scipy.signal.windows.hann。python_speech_features加的窗函數調用的numpy，如numpy.hanning，漢寧窗公式為：

給原信號加窗的實現方式為相乘：

frames *= 0.5 - 0.5 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1)) # 原信號乘以漢寧窗函數

4.STFT處理得到spectrum(頻譜，實際是多幀的)，shape = (n_ftt // 2 +1) * n_frames

fft = librosa.core.fft.get_fftlib() stft_matrix = fft.rfft(fft_window * frames)

librosa和python_speech_features都是調用的numpy進行傅里葉變換，numpy.fft.rfft(frames,NFFT)。

5.取絕對值得到magnitude spectrum/spectrogram(聲譜，包含時間維度，即多幀),shape = (n_ftt // 2 +1) * n_frames

magnitude_spectrum = numpy.abs(stft_matrix) # 承接上一步的STFT magnitude_spectrum = numpy.abs(librosa.core.spectrum.stft(wav, n_fft=n_fft, hop_length=hop_length,win_length=win_length, center=center,window=window, pad_mode=pad_mode))**power # librosa封裝的計算magnitude spectrum函數，power參數默認為1.0

6.取平方得到power spectrum/spectrogram(聲譜，包含時間維度，即多幀),shape = (n_ftt // 2 +1) * n_frames

power_spectrum = numpy.square(magnitude_spectrum) # 承接上一步的magnitude_spectrum power_spectrum = numpy.abs(librosa.core.spectrum.stft(wav, n_fft=n_fft, hop_length=hop_length,win_length=win_length, center=center,window=window, pad_mode=pad_mode))**power # librosa封裝的計算power spectrum函數，power參數設置為2.0

7.構造梅爾濾波器組，shape = n_mels * (n_ftt // 2 +1)

mel_basis = librosa.filters.mel(sr, n_fft, n_mels=128, fmin=0.0, fmax=None, htk=False,norm=1, dtype=np.float32) # 原函數及參數

librosa和python_speech_features構造mel_filtersbank的方式不同，librosa構造細節尚未掌握。

8.矩陣乘法得到mel_spectrogram，shape = n_mels * n_frames

mel_spectrogram = numpy.dot(mel_basis, power_spectrum)

[ n_mels ，(n_ftt // 2 +1) ] * [ (n_ftt // 2 +1) ，n_frames ] = [ n_mels，n_frames]

9.對mel_spectrogram進行log變換，shape = n_mels * n_frames

log_mel_spectrogram = librosa.core.spectrum.power_to_db(mel_spectrogram, ref=1.0, amin=1e-10, top_db=80.0)

S_db ~= 10 * log10(S) - 10 * log10(ref)

librosa采用以上方式計算log映射，python_speech_features定義mfcc( )函數時直接取自然常數e為底的對數。

10.IFFT變換，實際采用DCT得到MFCC，shape = n_mels * n_frames

mfcc = scipy.fftpack.dct(log_mel_spectrogram, type=2, n=None, axis=0, norm=None, overwrite_x=False) # n表示計算維度，需與log_mel_spectrogram.shape[axis]相同,否則作填充或者截斷處理。axis=0表示沿著自上而下的方向，分別選取每一行所在同一列的元素進行運算。

與python_speech_features相同，librosa也是調用scipy對log_mel_spectrogram進行離散余弦變換：scipy.fftpack.dct()。

11.取MFCC矩陣的低維(低頻)部分，shape = n_mfcc * n_frames

mfcc = mfcc[ :n_mfcc] # 取低頻維度上的部分值輸出，語音能量大多集中在低頻域，數值一般取13。

二、python_speech_features

1.源語音信號，shape = wav.length

sample_rate,signal = scipy.io.wavfile.read(filename, mmap=False) # scipy加載語音文件 signal, sample_rate = librosa.load(path, sr=22050, mono=True, offset=0.0, duration=None,dtype=np.float32, res_type='kaiser_best') # librosa加載語音文件

2.預加重，shape = wav.length

signal = numpy.append(signal[0],signal[1:] - coeff * signal[:-1])

語音信號的預加重，目的是為了對語音的高頻部分進行加重，去除口唇輻射的影響，增加語音的高頻部分的分辨率和信噪比，計算方式為：y(n)）=x(n)-a*x(n-1)，系數a一般取0.97。

3.分幀和加窗，shape = n_frames * n_ftt

frames = python_speech_features.sigproc.framesig(signal,frame_len,frame_step,winfunc=lambda x:numpy.ones((x,))) # 分幀及加窗函數

python_speech_features中默認幀長取25ms、幀移10ms，對應實際的采樣點為：

幀長點數 = 幀長時間 * 采樣率 ，幀移點數 = 幀移時間 * 采樣率

numframes = 1 + int(math.ceil((1.0*slen - frame_len)/frame_step)) # 計算總分幀數

總分幀數量 = 1 + 向上取整( ( 原信號長度 - 幀長點數 ) / 幀移點數 )

padlen = int((numframes-1)*frame_step + frame_len) # 計算待填充點的數量 zeros = numpy.zeros((padlen - slen,)) # 在原數據末尾進行零填充

待填充點數 = ( 總分幀數量 - 1 ) * 幀移點數 + 幀長點數

這種填充方式和上面librosa的自定義的那種填充方式( y = numpy.pad(wav, (0,pad_need), mode='constant') )效果是一樣的，可以看出librosa庫的功能還是更加多樣化的。

加窗的話，python_speech_features默認不加窗，但提供了調用numpy中窗函數的參數接口，經測試numpy.hanning窗函數和scipy.signal.windows.hann窗函數的數值是一致的，只不過前者為矩陣形式(元素相同的多個向量構成)，后者為向量形式。

mfcc = python_speech_features.base.mfcc(signal,samplerate=16000,winlen=0.025,winstep=0.01,numcep=13,nfilt=26,nfft=512,lowfreq=0,highfreq=None,preemph=0.97,ceplifter=22,appendEnergy=True,winfunc=lambda x:numpy.ones((x,))) # winfunc參數即為窗函數接口，例如 winfunc=numpy.hanning

4.STFT處理得到spectrum(頻譜，實際是多幀的)，shape = n_frames * (n_ftt // 2 +1)

complex_spec = numpy.fft.rfft(frames,n_fft) # 離散傅里葉變換得到頻譜圖

和librosa一致，python_speech_features也是調用numpy下的函數做離散傅里葉變換。

5.取絕對值得到magnitude spectrum/spectrogram(聲譜，包含時間維度，即多幀),shape = n_frames * (n_ftt // 2 +1)

mag_spec = numpy.absolute(complex_spec) # 承接上一步取絕對值得到幅度譜

6.取平方得到power spectrum/spectrogram(聲譜，包含時間維度，即多幀),shape = n_frames * (n_ftt // 2 +1)

power_spec = 1.0 / n_fft * numpy.square(mag_spec) # 承接上一步取平方得到功率譜

python_speech_features計算功率譜的方式和libraosa不一致，這里額外除以了傅里葉變換的維度n_fft。

7.構造梅爾濾波器組，shape = n_mels * (n_ftt // 2 +1)

fb = python_speech_features.base.get_filterbanks(nfilt=20,nfft=512,samplerate=16000,lowfreq=0,highfreq=None)

關于梅爾濾波器的構造細節，參考了這篇文章的部分介紹。具體細節如下：

赫茲轉梅爾：mel = 2595 * log10( 1+ hz / 700. )

梅爾轉赫茲：hz = 700 * ( 10** ( mel / 2595.0 ) - 1 )

(1)假設構造n_mels=10個梅爾濾波器，建立坐標系縱坐標表示mel刻度、橫坐標表示Hz刻度

(2)橫坐標(頻率)設置最低值300hz、最高值8000hz(采樣率的二分之一)，對應的縱坐標的最低和最高值分別為401.97mel和2840.02mel。

(3)把縱坐標從低到高劃分為12(10+2)個離散的點melpoints = [ 401.97， 623.61，845.25，1066.89，1288.54，1510.18， 1731.82，1953.46， 2175.10，2396.74，2618.38，2840.02] ，這些點映射回橫坐標上分別為[ 300，517.34，781.91，1103.98，1496.06，
1973.34，2554.36，3261.65，4122.66，5170.80，6446.75， 8000. ]。

(4)上面橫坐標這12(0-11)個點正好對應10個三角濾波器的三個點：例如0、1、2分別對應第一個濾波器的左底點、頂點、右底點，1、2、3對應第二個濾波器，以此類推第十個濾波器對應9、10、11。這十個三角濾波器的函數解析式如下：

(5)對信號的濾波處理本質上是乘法運算，由于待濾波的spectrogram實際是離散的數據，所以要對連續的濾波器函數進行離散的采樣，從源代碼來看python_speech_features和librosa的采樣方式是不同的，經測試發現兩者所構造的梅爾濾波器數據確實不一樣。

(6)拿第一個濾波器舉例來說，其橫坐標的采樣點選擇是根據經傅里葉變換后的spectrogram維度確定的：

bin = numpy.floor( ( n_fft + 1 ) * mel2hz( melpoints ) / samplerate ) # 例子中n_fft取值512,samplerate取值16000 # 構造梅爾采樣容器，實際上就是把原橫坐標"samplerate/2"上的12個離散點映射到新橫坐標"(n_ftt+1)/2"上面，這樣做保證了后面的乘法運算維度一致。 # bin = [ 9. 16. 25. 35. 47. 63. 81. 104. 132. 165. 206. 256.]，例子中n_fft取值512

經過以上映射后，第一個濾波器的函數表達式就確定了：在[9,16)區間為y=(x-9)/(16-9)，在[16,25]區間為y=(25-x)/(25-16)，其余點y取值為0。在這個濾波器的257(0-256)個離散點中，沿著三角形左邊y值逐漸遞增，在頂點處達到最大值"1"，然后順著右邊y值逐漸遞減到0，其余部分點對應的縱坐標取值都為0。

(7)例子中(n_mel=10,n_fft=512,samplerate=16000)構造的梅爾濾波器矩陣具體數值如下(作了轉置處理)：

8.矩陣乘法得到mel_spectrogram，shape = n_frames * n_mels

mel_spectrogram = numpy.dot(power_spec,fb.T)

[ n_frames ，(n_ftt // 2 +1) ] * [ (n_ftt // 2 +1) ，n_mels ] = [ n_frames，n_mels]

mel_spectrogram矩陣中的每個元素其實就是每一幀spectrogram向量和每個濾波器向量的內積。mel_spectrogram(0,1)即為第1幀的spectrogram同第2個三角濾波器內積運算的結果。

9.對mel_spectrogram進行log變換，shape = n_frames * n_mels

log_mel_spectrogram = numpy.log(mel_spectrogram) log_mel_spectrogram = python_speech_features.sigproc.logpowspec(frames,n_ftt,norm=1)

python_speech_features.base.mfcc( )函數默認取以自然常數e為底的對數，作者還提供了另一種取對數域的接口python_speech_features.sigproc.logpowspec( )，計算方式與librosa略有差異：log_S = 10 * log10(S)

10.IFFT變換，實際采用DCT得到MFCC，shape = n_frames * n_mels

mfcc = scipy.fftpack.dct(log_mel_spectrogram, type=2, n=None, axis=0, norm=None, overwrite_x=False) # n表示計算維度，需與log_mel_spectrogram.shape[axis]相同,否則作填充或者截斷處理。axis=0表示自上而下方向分別選取每一行所在同一列的元素進行計算。

11.取MFCC矩陣的低維(低頻)部分，shape = n_frames * n_mfcc

mfcc = mfcc[:,:n_mfcc] # 取低頻維度上的部分值輸出，語音能量大多集中在低頻域，數值一般取13。

python_speech_features.base.mfcc(appendEnergy=True )函數中，appendEnergy參數控制是否把MFCC的第一個倒譜系數替換為每一幀總能量的對數，每一幀總能量的計算方式為：

energy = numpy.sum(power_spec,axis=1) # axis=1表示沿著從左到右的方向,分別選取每一列所在同一行的元素進行運算。

12.倒譜提升，shape = n_frames * n_mfcc

mfcc = python_speech_features.base.lifter(cepstra=mfcc, L=22)

倒譜提升系數默認設置為22，具體實現方式為：

nframes,ncoeff = numpy.shape(cepstra) # ncoeff=n_mfcc n = numpy.arange(ncoeff) # n=0,1,2...n_mfcc-1 lift = 1 + (L / 2.) * numpy.sin(numpy.pi * n / L) # lift.shape = n_mfcc * 1 return lift * cepstra # shape = n_frames * n_mfcc

13.微分：mfcc的動態特征提取，shape = n_frames * n_mfcc

mfcc_delta_1 = python_speech_features.base.delta(feat=mfcc, N=1) # 計算mfcc的一階微分 mfcc_delta_2 = python_speech_features.base.delta(feat=mfcc, N=2) # 計算mfcc的二階微分

有時會把MFCC的基礎特征同其一階、二階微分數據結合起來使用，以做到特征層面的動靜結合。

總結

以上是生活随笔為你收集整理的librosa能量_语音MFCC提取：librosa amp;amp; python_speech_feature(2019.12)的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。