MFCC.rar_MFCC_P6W_feature_python mfcc_python mfcc图

import wave as we import numpy as np from scipy import * import matplotlib.pyplot as plt import math from scipy.fftpack import dct #读取wav文件 def wavred(): file=we.open('a0001.wav','r') s=file.getparams() framesrate,framesnum=s[2],s[3] datawav=file.readframes(s[3]) file.close() datause=np.fromstring(datawav,dtype=np.short) datatime=np.arange(0,framesnum)*(1.0/framesrate) return datause,datatime,framesrate #信号预加重 def pre_emphsis(signal,coefficient=0.95): return np.append(signal[0],signal[1:]-coefficient*signal[:-1]) #分帧及加窗 def data2frame(signal,frame_length,frame_step): signal_length=len(signal) #信号总长度 frame_length=int(round(frame_length)) #帧长 frame_step=int(round(frame_step)) #帧移 #若信号长度小于一个帧的长度，帧数为1，否则计算帧的总个数 if signal_length<=frame_length: frame_num=1 else: frame_num=1+int(math.ceil((1.0*signal_length-frame_length)/frame_step)) #所有帧加起来的总的长度 pad_length=int((frame_num-1)*frame_step+frame_length) #多出的长度使用0填补 zeros=np.zeros((pad_length-signal_length,)) #将数组signal和zeros进行拼接，扩充了数组 pad_signal=np.concatenate((signal,zeros)) #对所有帧的时间点进行抽取，得到frames_num*frame_length长度的矩阵 indices=np.tile(np.arange(0,frame_length),(frame_num,1))+np.tile(np.arange(0,frame_num*frame_step,frame_step),(frame_length,1)).T #将indices转化为矩阵 indices=np.array(indices,dtype=np.int32) #帧矩阵 frames=pad_signal[indices] #加窗 win = np.hamming(frame_length) for i in range(0, frame_num - 1): frames[i] = frames[i] * win return frames #计算每一帧傅里叶变换后的频谱的幅度 #NFFT：FFT的大小 def spectrum_magnitude(frames,NFFT): complex_spectrum=np.fft.rfft(frames,NFFT) return np.absolute(complex_spectrum) #计算每一帧傅里叶变换以后的功率谱 def spectrum_power(frames,NFFT): return 1.0/NFFT*np.square(spectrum_magnitude(frames,NFFT)) #计算每一帧的功率谱的对数形式 def log_spectrum_power(frames,NFFT,norm=1): spec_power=spectrum_power(frames,NFFT) spec_power[spec_power<1e-30]=1e-30 log_spec_power=10*np.log10(spec_power) if norm: return log_spec_power-np.max(log_spec_power) else: return log_spec_power #把频率转化为梅尔频率 def hz2mel(hz): return 2595*np.log10(1+hz/700.0) #把梅尔频率转化为频率 def mel2hz(mel): return 700*(10**(mel/2595.0)-1) #构建梅尔三角间距滤波器 def get_filter_banks(filters_num,NFFT,samplerate,low_freq,high_freq): #将最低和最高频率转化为梅尔频率 low_mel=hz2mel(low_freq) high_mel=hz2mel(high_freq) #在low_mel和high_mel之间等间距插入filters_num个点,一共filters_num+2个点 mel_points=np.linspace(low_mel,high_mel,filters_num+2) #将梅尔频率转化为hz频率,并找到对应的hz位置 hz_points=mel2hz(mel_points) #将hz_points对应倒fft中的位置 bin=np.floor((NFFT+1)*hz_points/samplerate) #建立滤波器的表达式,每个滤波器在第1点和第3点处均为0,其余为三角形形状 fbank=np.zeros([filters_num,int(NFFT/2)+1]) for j in range(0,filters_num): for i in range(int(bin[j]),int(bin[j+1])): fbank[j, i] = (i - bin[j]) / (bin[j + 1] - bin[j]) for i in range(int(bin[j + 1]), int(bin[j + 2])): fbank[j, i] = (bin[j + 2] - i) / (bin[j + 2] - bin[j + 1]) return fbank #升倒谱函数,cepstra:MFCC系数,L:升系数,默认为22 def lifter(cepstra,L=22): if L>0: nframes,ncoeff=np.shape(cepstra) n=np.arange(ncoeff) lift=1+(L/2)*np.sin(np.pi*n/L) return lift*cepstra else: return cepstra #计算一阶系数或者加速系数的一般变换公式 def derivate(feat,big_theta=2,cep_num=13): result=np.zeros(feat.shape) denominator=0 for theta in np.linspace(1,big_theta,big_theta): denominator=denominator+theta**2 denominator=denominator*2 for row in np.linspace(0,feat.shape[0]-1,feat.shape[0]): #feat.shape[0]=890 tmp=np.zeros((cep_num,)) numerator=np.zeros((cep_num,)) for t in np.linspace(1,cep_num,cep_num): a=0 b=0 s=0 for theta in np.linspace(1,big_theta,big_theta): if (t+theta)>cep_num: a=0 else: a=feat[int(row)][int(t)+int(theta)-1] if (t-theta)<1: b=0 else: b=feat[int(row)][int(t)-int(theta)-1] s+=theta*(a-b) numerator[int(t)-1] =s tmp=numerator*1.0/denominator result[int(row)]=tmp return result #读取wav文件，wavdata是71332个采样点的频率值,wavtime是采样时间 wavdata, wavtime,fs= wavred() fft_size=256 #傅里叶变换的采样个数,即一帧的长度 filters_num=26 #梅尔滤波器个数 high_freq=fs/2 #声音信号频率的最大值 low_freq=0 #声音信号频率的最小值 cep_num=13 #倒谱系数的个数 cep_lifter=26 appendEnergy=True #预加重 sig1=pre_emphsis(wavdata) #分帧和加窗 sig2=data2frame(sig1,0.128*fs,0.04*fs) #傅里叶变换 spec_power=spectrum_power(sig2,fft_size) energe=np.sum(spec_power,1) energe=np.where(energe==0,np.finfo(float).eps,energe) #对能量为0的地方调整为esp,便于进行对数处理 #计算13个MFCC #获得一个梅尔滤波器组 fb=get_filter_banks(filters_num,fft_size,fs,low_freq,high_freq) #使用梅尔刻度滤波器组过滤 feat=np.dot(spec_power,fb.T) feat=np.where(feat==0,np.finfo(float).eps,feat) #取对数 feat=np.log(feat) #进行离散余弦变换,只取前13个系数 feat=dct(feat,type=2,axis=1,norm='ortho')[:,:cep_num] #升倒谱 feat=lifter(feat,cep_lifter) #只取2-13个系数,第一个用能量的对数来代替 if appendEnergy: feat[:,0]=np.log(energe) #计算差分 result1=derivate(feat) result2=derivate(result1) result3=np.concatenate((feat,result1),axis=1) result=np.concatenate((result3,result2),axis=1)