#coding=utf-8
#計算每一幀的MFCC系數

import?numpy?as?np
from?preprocessing?import?audio2frame
from?preprocessing?import?pre_emphasis
from?preprocessing?import?spectrum_power
from?scipy.fftpack?import?dct
import?myplot
import?math

#適配3.x版本的xrange
try:
????xrange（1）
except:
????xrange=range

def?envelope（x）:
????n=len（x）
????env=np.zeros（n）
????for?i?in?range（n）:
????????j=max（0?i-5）
????????k=min（i+5?n-1）
????????env[i]=np.mean（x[j:k]）
????return?env

def?fbank（signal?samplerate=16000?win_length=0.025?win_step=0.01?filters_num=26?NFFT=512?\
????low_freq=0?high_freq=None?pre_emphasis_coeff=0.97）:
????‘‘‘
????計算音頻信號的經過梅爾三角濾波的能量譜
????參數說明：
????samplerate:采樣頻率
????win_length:窗長度，每一幀的時間間隔
????win_step:窗間隔
????filters_num:梅爾濾波器個數
????NFFT:FFT采樣點個數
????low_freq:最低頻率
????high_freq:最高頻率
????pre_emphasis_coeff:預加重系數
????‘‘‘
????high_freq=high_freq?or?samplerate/2?#計算音頻樣本的最大頻率
????signal_emp=pre_emphasis（signal?pre_emphasis_coeff）#預加重處理
????raw_spectrum=np.absolute（np.fft.fft（signal））
????emp_spectrum=np.absolute（np.fft.fft（signal_emp））
????#畫出預加重前后的頻域圖
????myplot.plot_spectrum（10*np.log10（raw_spectrum）?10*np.log10（emp_spectrum）?samplerate）

????#畫出原始和預加重之后的時域圖
????myplot.plot_signal1（signal?signal_emp?samplerate）

????frames=audio2frame（signal_emp?win_length*samplerate?win_step*samplerate）#得到幀數組
????#畫出分幀之后的每一幀加上hamming窗之后的時域圖，在此之前，在audio2frame中會先將原始的frame畫出來
????myplot.plot_frames（frames?np.shape（frames）[0]?np.shape（frames）[1]）
????#myplot.plot_hamming（frames[77:]）

????spec_power=spectrum_power（frames?NFFT）#得到每一幀的能量譜
????env=envelope（spec_power[1?:]）
????#畫出某一幀的頻譜包絡圖，這里選取第一幀
????myplot.plot_envelope（10*np.log10（spec_power[1?:]）?10*np.log10（env）?NFFT?samplerate）

????energy=np.sum（spec_power?1）#對每一幀的能量譜求和
????energy=np.where（energy==0?np.finfo（float）.eps?energy）#對能量為0的地方調整為eps，?方便進行對數處理
????#畫出每一幀沒有經過濾波器組的功率譜和熱圖
????myplot.plot_power_before_mel（10*np.log10（spec_power）?NFFT?np.shape（frames）[0]?samplerate）
????#myplot.plot_power_heatmap_before_mel（10*np.log10（spec_power））

????fb?freq_point=get_filter_banks（filters_num?NFFT?samplerate?low_freq?high_freq）#獲得每一個濾波器的頻率寬度
????#畫出濾波器組的形狀
????myplot.plot_fbank（fb?NFFT?filters_num?samplerate）

????feat=np.dot（spec_power?fb.T）#對濾波器和能量譜進行點乘
????#得到num_frames*num_filter大小的矩陣，也就是得到了每一幀用濾波器組計算（將能量用濾波器進行加重）之后的結果
????feat=np.where（feat==0?np.finfo（float）.eps?feat）#同樣不能出現0
????#畫出經過濾波器組濾波的每一幀功率譜，注意這個時候的頻率標度已經不是線性的了，是按照濾波器中心頻率來標的
????myplot.plot_power_after_mel（freq_point?np.arange（0?np.shape（frames）[0]?1）?10*np.log10（feat））
????return?feat?energy

def?get_filter_banks（filters_num=20?NFFT=512?samplerate=16000?low_freq=0?high_freq=None）:
????‘‘‘
????計算梅爾三角間距濾波器，該濾波器在第一個和第三個頻率處為0，在第二個頻率處為1
????參數說明：
????filters_num:?濾波器個數
????NFFT:FFT的總計的個數
????samplerate:采樣頻率
????low_freq:最低頻率
????high_freq:最高頻率
????‘‘‘
????#首先，將頻率hz轉化為梅爾頻率，因為人耳分辨聲音的大小與頻率并非線性正比，所以化為梅爾頻率再線性分割
????low_mel=hz2mel（low_freq）
????high_mel=hz2mel（high_freq）
????#需要在low_mel和high_mel之間等間距插入filters_num個點，一個filters_num+2個點

?屬性????????????大小?????日期????時間???名稱
-----------?---------??----------?-----??----

?????文件???????9378??2018-03-20?21:46??code\calcmfcc.py

?????文件???????4639??2018-03-20?22:20??code\myplot.py

?????文件???????5251??2018-03-20?10:04??code\preprocessing.py

?????文件????????313??2018-03-20?10:05??code\test.py

?????文件??????41964??2018-03-17?23:29??code\test.wav

?????目錄??????????0??2018-03-21?08:30??code

-----------?---------??----------?-----??----

????????????????61545????????????????????6