梅尔倒谱系数(MFCC)实现
本文实例为大家分享了梅尔倒谱系数实现代码,供大家参考,具体内容如下
""" @author:zoutai @file:mymfcc.py @time:2018/03/26 @description: """ frommatplotlib.colorsimportBoundaryNorm importlibrosa importlibrosa.display importnumpy importscipy.io.wavfile fromscipy.fftpackimportdct importmatplotlib.pyplotasplt importnumpyasnp #第一步-读取音频,画出时域图(采样率-幅度) sample_rate,signal=scipy.io.wavfile.read('OSR_us_000_0010_8k.wav')#Fileassumedtobeinthesamedirectory signal=signal[0:int(3.5*sample_rate)] #plotthewave time=np.arange(0,len(signal))*(1.0/sample_rate) #plt.plot(time,signal) plt.xlabel("Time(s)") plt.ylabel("Amplitude") plt.title("SignalintheTimeDomain") plt.grid('on')#标尺,on:有,off:无。 #第二步-预加重 #消除高频信号。因为高频信号往往都是相似的, #通过前后时间相减,就可以近乎抹去高频信号,留下低频信号。 #原理:y(t)=x(t)−αx(t−1) pre_emphasis=0.97 emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1]) time=np.arange(0,len(emphasized_signal))*(1.0/sample_rate) #plt.plot(time,emphasized_signal) #plt.xlabel("Time(s)") #plt.ylabel("Amplitude") #plt.title("SignalintheTimeDomainafterPre-Emphasis") #plt.grid('on')#标尺,on:有,off:无。 #第三步、取帧,用帧表示 frame_size=0.025#帧长 frame_stride=0.01#步长 #frame_length-一帧对应的采样数,frame_step-一个步长对应的采样数 frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate#Convertfromsecondstosamples 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))#Makesurethatwehaveatleast1frame pad_signal_length=num_frames*frame_step+frame_length z=numpy.zeros((pad_signal_length-signal_length)) pad_signal=numpy.append(emphasized_signal,z)#PadSignaltomakesurethatallframeshaveequalnumberofsampleswithouttruncatinganysamplesfromtheoriginalsignal #ConstructanarraybyrepeatingA(200)thenumberoftimesgivenbyreps(348). #这个写法太妙了。目的:用矩阵来表示帧的次数,348*200,348-总的帧数,200-每一帧的采样数 #第一帧采样为0、1、2...200;第二帧为80、81、81...280..依次类推 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 frames=pad_signal[indices.astype(numpy.int32,copy=False)]#Copyofthearrayindices #frame:348*200,横坐标348为帧数,即时间;纵坐标200为一帧的200毫秒时间,内部数值代表信号幅度 #plt.matshow(frames,cmap='hot') #plt.colorbar() #plt.figure() #plt.pcolormesh(frames) #第四步、加汉明窗 #傅里叶变换默认操作的时间段内前后端点是连续的,即整个时间段刚好是一个周期, #但是,显示却不是这样的。所以,当这种情况出现时,仍然采用FFT操作时, #就会将单一频率周期信号认作成多个不同的频率信号的叠加,而不是原始频率,这样就差生了频谱泄漏问题 frames*=numpy.hamming(frame_length)#相乘,和卷积类似 ##frames*=0.54-0.46*numpy.cos((2*numpy.pi*n)/(frame_length-1))#ExplicitImplementation** #plt.pcolormesh(frames) #第五步-傅里叶变换频谱和能量谱 #_raw_fft扫窗重叠,将348*200,扩展成348*512 NFFT=512 mag_frames=numpy.absolute(numpy.fft.rfft(frames,NFFT))#MagnitudeoftheFFT pow_frames=((1.0/NFFT)*((mag_frames)**2))#PowerSpectrum #plt.pcolormesh(mag_frames) # #plt.pcolormesh(pow_frames) #第六步,FilterBanks滤波器组 #公式:m=2595*log10(1+f/700);f=700(10^(m/2595)−1) nfilt=40#窗的数目 low_freq_mel=0 high_freq_mel=(2595*numpy.log10(1+(sample_rate/2)/700))#ConvertHztoMel mel_points=numpy.linspace(low_freq_mel,high_freq_mel,nfilt+2)#EquallyspacedinMelscale hz_points=(700*(10**(mel_points/2595)-1))#ConvertMeltoHz bin=numpy.floor((NFFT+1)*hz_points/sample_rate) fbank=numpy.zeros((nfilt,int(numpy.floor(NFFT/2+1)))) forminrange(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 forkinrange(f_m_minus,f_m): fbank[m-1,k]=(k-bin[m-1])/(bin[m]-bin[m-1]) forkinrange(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)#NumericalStability filter_banks=20*numpy.log10(filter_banks)#dB;348*26 #plt.subplot(111) #plt.pcolormesh(filter_banks.T) #plt.grid('on') #plt.ylabel('Frequency[Hz]') #plt.xlabel('Time[sec]') #plt.show() # #第七步,梅尔频谱倒谱系数-MFCCs num_ceps=12#取12个系数 cep_lifter=22#倒谱的升个数?? mfcc=dct(filter_banks,type=2,axis=1,norm='ortho')[:,1:(num_ceps+1)]#Keep2-13 (nframes,ncoeff)=mfcc.shape n=numpy.arange(ncoeff) lift=1+(cep_lifter/2)*numpy.sin(numpy.pi*n/cep_lifter) mfcc*=lift#* #plt.pcolormesh(mfcc.T) #plt.ylabel('Frequency[Hz]') #plt.xlabel('Time[sec]') #第八步,均值化优化 #tobalancethespectrumandimprovetheSignal-to-Noise(SNR),wecansimplysubtractthemeanofeachcoefficientfromallframes. filter_banks-=(numpy.mean(filter_banks,axis=0)+1e-8) mfcc-=(numpy.mean(mfcc,axis=0)+1e-8) #plt.subplot(111) #plt.pcolormesh(mfcc.T) #plt.ylabel('Frequency[Hz]') #plt.xlabel('Time[sec]') #plt.show() #直接频谱分析 #plotthewave #plt.specgram(signal,Fs=sample_rate,scale_by_freq=True,sides='default') #plt.ylabel('Frequency(Hz)') #plt.xlabel('Time(s)') #plt.show() plt.figure(figsize=(10,4)) mfccs=librosa.feature.melspectrogram(signal,sr=8000,n_fft=512,n_mels=40) librosa.display.specshow(mfccs,x_axis='time') plt.colorbar() plt.title('MFCC') plt.tight_layout() plt.show()
以上就是本文的全部内容,希望对大家的学习有所帮助,也希望大家多多支持毛票票。