#!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile
from scipy.fftpack import dct
# 原始数据
sample_rate,signal=scipy.io.wavfile.read('test.wav')
#截取
original_signal=signal[0:int(1.5*sample_rate)]
#预加重
pre_emphasis=0.97
emphasized_signal=np.append(original_signal[0],original_signal[1:]-pre_emphasis*original_signal[:-1])
#分帧
frame_size=0.025
frame_stride=0.01
frame_length=int(round(frame_size*sample_rate))
frame_step=int(round(frame_stride*sample_rate))
signal_length=len(emphasized_signal)
num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step))
pad_signal_length=num_frames*frame_step+frame_length
pad_signal=np.append(emphasized_signal,np.zeros((pad_signal_length-signal_length)))
indices=np.tile(np.arange(0,frame_length),(num_frames,1))+np.tile(np.arange(0,num_frames*frame_step,frame_step),(frame_length,1)).T
frames=pad_signal[np.mat(indices).astype(np.int32,copy=False)]
#加窗
frames*=np.hamming(frame_length)
#快速傅里叶变换FFT和功率谱
NFFT=512
mag_frames=np.absolute(np.fft.rfft(frames,NFFT)) # Magnitude of the FFT
pow_frames=(1.0/NFFT)*(mag_frames**2)
#三角带通滤波器（Mel滤波），将频率转换为Mel频率，filterbank到此为止
low_freq_mel=0
nfilt=40
high_freq_mel=(2595*np.log10(1+(sample_rate/2)/700))
mel_points=np.linspace(low_freq_mel,high_freq_mel,nfilt+2) # Equally spaced in Mel scale
hz_points=(700*(10**(mel_points/2595)-1)) # Convert Mel to Hz
bin = np.floor((NFFT+1)*hz_points/sample_rate)
fbank=np.zeros((nfilt,int(np.floor(NFFT/2+1))))
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=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
# 经离散余弦变换（DCT）得到MFCC系数
num_ceps=12
mfcc=dct(filter_banks,type=2,axis=1,norm='ortho')[:,1:(num_ceps)+1]
(nframes,ncoeff)=mfcc.shape
# 将正弦升降1应用于MFCC以降低已被声称在噪声信号中改善语音识别的较高MFCC
n=np.arange(ncoeff)
cep_lifter=22
lift=1+(cep_lifter/2)*np.sin(np.pi*n/cep_lifter)
mfcc*=lift
# 平均归一化滤波器组
filter_banks-=(np.mean(filter_banks,axis=0)+1e-8)
# 平均归一化MFCC
mfcc-=(np.mean(mfcc,axis=0)+1e-8)
signal_num=np.arange(len(signal))
sample_num=np.arange(len(original_signal))
emphasized_signal_num=np.arange(len(emphasized_signal))
pad_signal_num=np.arange(len(pad_signal))
plt.figure(figsize=(11,7))
plt.subplot(811)
plt.plot(signal_num/sample_rate,signal,color='black')
plt.plot(sample_num/sample_rate,original_signal,color='blue')
plt.ylabel('Amplitude')
plt.title('signal of Voice')
plt.subplot(812)
plt.plot(sample_num/sample_rate,original_signal,color='blue')
plt.xlabel('Time(sec)')
plt.ylabel('Amplitude')
plt.title('1.5s signal of Voice')
plt.subplot(813)
plt.plot(emphasized_signal_num/sample_rate,emphasized_signal,color='blue')
plt.xlabel('Time(sec)',fontsize=14)
plt.ylabel('Amplitude',fontsize=14)
plt.title('emphasized signal of Voice',fontsize=14)
plt.subplot(814)
plt.plot(pad_signal_num/sample_rate,pad_signal,color='blue')
plt.xlabel('Time(sec)',fontsize=14)
plt.ylabel('Amplitude',fontsize=14)
plt.title('pad signal of Voice',fontsize=14)
plt.subplot(815)
plt.imshow(np.flipud(filter_banks.T),cmap=plt.cm.jet,aspect=0.2,extent=[0,filter_banks.shape[1],0,filter_banks.shape[0]]) # 画热力图
plt.title('MFCC')
plt.subplot(816)
plt.imshow(np.flipud(mfcc.T),cmap=plt.cm.jet,aspect=0.2,extent=[0,mfcc.shape[0],0,mfcc.shape[1]])# 热力图
plt.title('MFCC')
plt.show()
参考链接：语音识别第4讲：语音特征参数MFCC
本文创建于2022.10.25/21.50，修改于2022.10.25/21.50