MFCC语音特征值提取算法(二)

简介: MFCC语音特征值提取算法(二)

计算MFCC系数


由于信号在时域上的变换很难看出特征,因此,我们通常将它转换为频域上的能量分布以便于观察。不同的能量分布,代表不同语音的特征。语音原信号在与窗函数(如汉明窗)相乘后,每帧还必须再经过快速傅里叶变换以得到频谱上的能量分布。对语音信号分帧加窗后的各帧的频谱,然后对频谱进行取模平方运算后即为语音信号的功率谱。


对信号幅度谱、功率谱以及对数功谱的计算实例代码如下:


import numpy
import logging
def msgspec(frames,NFFT):
    """
    计算帧中每个帧的幅度谱。如果帧为N*D,则输出N*(NFFT/2+1)
    """
    if numpy.shape(frames)[1]>NFFT:
        logging.warn('frame length (%d)is greater than FFT size(%d),frame will be truncated .Increase NFFT to avoid.',numpy.shape(frames)[1],NFFT)
        complex_spec=numpy.fft.rfft(frames,NFFT)
        return numpy.absolute(complex_spec)
def power_spectrum(frames,NFFT):
    return 1.0/NFFT*numpy.square(spectrum_magnitude(frames,NFFT))
def log_power_spectrum(frames,NFFT,norm=1):
    spec_power=power_spectrum(frames,NFFT)
    spec_power[spec_power<1e-30]
    log_spec_power=10*numpy.log10(spec_power)
    if norm:
        return log_spec_power-numpy.max(log_spec_power)
    else:
        return log_spec_power


此外,信号的每一帧的音量(即能量),也是语音的特征,而且非常容易计算。因此,通常会再加上一帧的能量,使得每一帧基本的语音特征增加一个维度,包括一个对数能量和倒谱参数。标准的倒谱参数MFCC,只反映了语音参数的静态特征,语音参数的动态特征可以用这些静态特征的差分普来描述。


MFCC的全部组成如下:N维MFCC系数(N/3 MFCC系数+N/3 一阶差分系数+N/3二阶差分系数)+帧能量。以语音识别中常用的39维MFCC为例,即为:13个静态系数+13个一阶差分系数(Delta系数)+13个二阶差分系数(Delta-Delta系数)。其中,差分系数用来描述动态特征,即声学特征在相邻帧间的变化情况。


在MFCC计算中还涉及频率与梅尔刻度之间的转换,其转换方式如下:


m=2595lg(1+f700)m=2595lg(1+f700)
def hz2mel(hz):
    return 2595*numpy.log10(1+hz/700.0)


同样,我们也可以推出下列公式:


f=700(10m/2595−1)f=700(10m/2595−1)


Delta系数的计算公式为:


dt=∑Nn=1n(ct+n−ct−n)2∑Nn=1n2dt=∑n=1Nn(ct+n−ct−n)2∑n=1Nn2


其中,dtdt为Delta系数,从帧t根据静态系数ct−N到ct+Nct−N到ct+N 计算而得。N一般取值为2。Delta-Delta(加速度)系数的计算方法相同,但他们是根据Delta而不是静态系数来进行计算得到的。计算Deltaa系数的示例代码如下:


def delta(feat,N):
    if N<1:
        raise ValueError('N must be an integer>=1')
    NUMFRAMES=len(feat)
    denominator=2*sum([i**2 for i in range(1,N+1)])
    delta_feat=numpy.pad(feat,((N,N),(0,0)),mode='edge')
    for t in range(NUMFRAMES):
        delta_feat[t]=numpy.dot(numpy.arange(-N,N+1),padded[t:t+2*N+1])/denominator
        return delta_feat


当然除了自己定义函数,也可以直接使用工具包中的API。


0535071f203c890b5be5ad29974466dc_b70fb94e68c44c91b8af0b80ee9a8c87.png


对语音信号进行预加重


import numpy as np
import matplotlib.pyplot as plt
from python_speech_features.sigproc import *
from python_speech_features import *
from scipy.fftpack import dct
import scipy.io.wavfile as wav
sample_rate,signal=wav.read('./test.wav')
#保留语音的前3.5秒
signal=signal[0:int(3.5*sample_rate)]
#信号预加重
emphasized_signal=preemphasis(signal,coeff=0.95)
#显示信号
plt.plot(signal)
plt.title("Original Signal")
plt.plot(emphasized_signal)
plt.title("Preemphasis Signal")
plt.show()


df309e78d4330ae3325f27901bde0b91_4959d18714ce4d43aab644430332af67.png


上述示例代码,对信号进行预加重处理的是preemphasis(signal,coeff)函数,除了这个函数,也可以使用以下代码实现:


264b526e632240a5064fecae38ddd063_b8e874baf90b41db812eb18801abab9a.png


pre_emphasis=0.95

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


源代码:


import numpy as np
import matplotlib.pyplot as plt
from python_speech_features.sigproc import *
from python_speech_features import *
from scipy.fftpack import dct
import scipy.io.wavfile as wav
sample_rate,signal=wav.read('./test.wav')
pre_emphasis=0.95
emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1])
#保留语音的前3.5秒
#signal=signal[0:int(3.5*sample_rate)]
#信号预加重
#emphasized_signal=preemphasis(signal,coeff=0.95)
#显示信号
plt.plot(signal)
plt.title("Original Signal")
plt.plot(emphasized_signal)
plt.title("Preemphasis Signal")
plt.show()


通过上面的程序可知,两种函数都可以进行预加重处理,可以自行选择合适的方法。


对语音信号进行短时傅里叶变换

在对语音信号进行处理之前,我们需要对不稳定的语音信号进行短时分帧以获取傅里叶变换必需的稳定信号。语音处理范围内的典型帧大小范围为20ms~40ms,连续帧之间重叠50%左右。因此一般将帧长度设置为25ms。短时傅里叶变换(Short-Time Fourier Transform,SIFT)在MFCC计算过程中主要用于短时分帧处理后,通过对信号进行时域到频域的转换来获取语音信号的频谱。


#对信号进行短时分帧处理
frame_size=0.025 #设置帧长
#计算帧对应采样数(frame_length)以及步长对应采样数(frame_step)
frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate
signal_length=len(emphasized_signal) #信号总采样数
frame_length=int(round(frame_length)) #帧采样数
frame_step=int(round(frame_step))
#num_frames为总帧数,确保我们至少有一个帧
num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step))
pad_signal_length=num_frames*frame_step+frame_length
z=np.zeros((pad_signal_length-signal_length))
#填充信号以后确保所有的帧的采样数相等
pad_signal=np.append(emphasized_signal,z)
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[indices.astype(np.int32,copy=False)]


信号经过短时分帧之后,可通过短时傅里叶变换得到各种频谱


NFFT=512
mag_frames=np.absolute(np.fft.rfft(frames,NFFT))
pow_frames=((1.0/NFFT)*((mag_frames)**2))
log_pow_frames=logpowspec(pow_frames,NFFT,norm=1)
#保留语音的前3.5秒
#signal=signal[0:int(3.5*sample_rate)]
#信号预加重
#emphasized_signal=preemphasis(signal,coeff=0.95)
#显示信号
plt.plot(mag_frames)
plt.title("Mag_Spectrum")
plt.plot(emphasized_signal)
plt.show()
plt.plot(pow_frames)
plt.title("Power_Spectrum")
plt.show()
plt.plot(pow_frames)
plt.title("Log_Power_Spectrum")
plt.show()


运行上面的程序,就可以得到处理结果,下面展示原有的所有代码:


import numpy as np
import matplotlib.pyplot as plt
from python_speech_features.sigproc import *
from python_speech_features import *
from scipy.fftpack import dct
import scipy.io.wavfile as wav
sample_rate,signal=wav.read('./test.wav')
pre_emphasis=0.95
emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1])
#对信号进行短时分帧处理
frame_size=0.025 #设置帧长
frame_stride=0.1
#计算帧对应采样数(frame_length)以及步长对应采样数(frame_step)
frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate
signal_length=len(emphasized_signal) #信号总采样数
frame_length=int(round(frame_length)) #帧采样数
frame_step=int(round(frame_step))
#num_frames为总帧数,确保我们至少有一个帧
num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step))
pad_signal_length=num_frames*frame_step+frame_length
z=np.zeros((pad_signal_length-signal_length))
#填充信号以后确保所有的帧的采样数相等
pad_signal=np.append(emphasized_signal,z)
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[indices.astype(np.int32,copy=False)]
NFFT=512
mag_frames=np.absolute(np.fft.rfft(frames,NFFT))
pow_frames=((1.0/NFFT)*((mag_frames)**2))
log_pow_frames=logpowspec(pow_frames,NFFT,norm=1)
#保留语音的前3.5秒
#signal=signal[0:int(3.5*sample_rate)]
#信号预加重
#emphasized_signal=preemphasis(signal,coeff=0.95)
#显示信号
plt.plot(mag_frames)
plt.title("Mag_Spectrum")
plt.plot(emphasized_signal)
plt.show()
plt.plot(pow_frames)
plt.title("Power_Spectrum")
plt.show()
plt.plot(log_pow_frames)
plt.title("Log_Power_Spectrum")
plt.show()


aa330eeef5e027a43952fe5f2ffb1adb_b50dedfc08ca498f881df43c352285f0.png


                                                       (a)幅度谱



9ee14eabbf43e7a36a640bf7619a4c98_5e39d86817a646b4a9a9167dc131da53.png

                                                            (b)功率谱


e031a38d2fed08de9ecb5654a1409efb_d47d02c4dbe745d3afc680f6db1e2f9b.png

                                                            (c)功率对数谱


音频文件使用不同,最终结果也会不同,大家自己使用自己的音频,注意音频格式为“.wav”



定义滤波器组


将信号通过一组梅尔刻度的三角形滤波器组,采用的滤波器为三角形滤波器,中心频率为f(m),m=1,2,3,```````,M,M通常取22~26. 各f(m)之间的间隔随着m值的减少而减少。随着m值的增大而增大。如图:


387a879ee0f9064ebb34724295656d64_56f588c01d7c4b6e8a40c5ce45549153.png


三角形滤波器的频率响应定义公式:4


Hm(k)=⎧⎩⎨⎪⎪2(k−f(m−1))(f(m+1)−f(m−1)(f(m)−f(m−1)))2(f(m+1)−k)(f(m+1)−f(m−1)(f(m+1)−f(m−1)))f(m−1)≤k≤f(m+1)f(m)≤k≤f(m+1)Hm(k)={2(k−f(m−1))(f(m+1)−f(m−1)(f(m)−f(m−1)))f(m−1)≤k≤f(m+1)2(f(m+1)−k)(f(m+1)−f(m−1)(f(m+1)−f(m−1)))f(m)≤k≤f(m+1)


对于其他的情况,例如,k<f(m-1)和k>=f(m+1)则为0,当k=f(m)时为1.


定义梅尔刻度的三角形滤波器组的示例代码为:


low_freq_MEL=0 #将频率转换为梅尔刻度
nfilt=40 #窗的数目
#计算m=2595*log10(1+f/700)
high_freq_mel=(2595*np.log10(1+(sample_rate/2)/700))
mel_points=np.linspace(low_freq_MEL,high_freq_mel,nfilt+2) #梅尔刻度的均匀分布
#计算f=700(10**(m/2595)-1)
hz_points=(700*(10**(mel_points/2595)-1))
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]) #三角形滤波器左边频率f(m-1)
    f_m=int(bin[m]) #三角形滤波器中间频率fm
    f_m_plus=int(bin[m+1]) #三角形滤波器右边频率f(m-1)
    for k in range(f_m_minus,f_m):
        fbank[m-1,k]=(k-bin[m-1])/(bin[m+1]-bin[m])
plt.plot(fbank.T)
plt.show()


f3afcc27d9682a6785fd0c9289ae90e8_5014668d42d1453fb4fc98ed6cdd85ac.png


三角形滤波器有两个主要功能,其一,对频谱进行平滑并消除谐波的作用,突显原先语音的共振峰;其二,用以降低运算量。如图所示的滤波器组中的每个滤波器在中心频率处响应为1,并朝着0线性减少,直至达到响应为0的两个相邻滤波器的中心频率。


计算MFCC系数


如果计算出的滤波器组系数高度相关,则在某些机器学习算法中可能会存在问题。我们可用离散余弦变换对滤波器组系数进行去相关,并产生滤波器组的压缩表示。滤波器组输出的对数能量经离散余弦变换后,即可得到MFCC系数。示例代码如下:


import numpy as np
import matplotlib.pyplot as plt
from python_speech_features.sigproc import *
from python_speech_features import *
from scipy.fftpack import dct
import scipy.io.wavfile as wav
sample_rate,signal=wav.read('./test.wav')
pre_emphasis=0.95
emphasized_signal=numpy.append(signal[0],signal[1:]-pre_emphasis*signal[:-1])
#对信号进行短时分帧处理
frame_size=0.025 #设置帧长
frame_stride=0.1
#计算帧对应采样数(frame_length)以及步长对应采样数(frame_step)
frame_length,frame_step=frame_size*sample_rate,frame_stride*sample_rate
signal_length=len(emphasized_signal) #信号总采样数
frame_length=int(round(frame_length)) #帧采样数
frame_step=int(round(frame_step))
#num_frames为总帧数,确保我们至少有一个帧
num_frames=int(np.ceil(float(np.abs(signal_length-frame_length))/frame_step))
pad_signal_length=num_frames*frame_step+frame_length
z=np.zeros((pad_signal_length-signal_length))
#填充信号以后确保所有的帧的采样数相等
pad_signal=np.append(emphasized_signal,z)
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[indices.astype(np.int32,copy=False)]
NFFT=512
mag_frames=np.absolute(np.fft.rfft(frames,NFFT))
pow_frames=((1.0/NFFT)*((mag_frames)**2))
log_pow_frames=logpowspec(pow_frames,NFFT,norm=1)
#保留语音的前3.5秒
#signal=signal[0:int(3.5*sample_rate)]
#信号预加重
#emphasized_signal=preemphasis(signal,coeff=0.95)
#显示信号
'''
plt.plot(mag_frames)
plt.title("Mag_Spectrum")
plt.plot(emphasized_signal)
plt.show()
plt.plot(pow_frames)
plt.title("Power_Spectrum")
plt.show()
plt.plot(log_pow_frames)
plt.title("Log_Power_Spectrum")
plt.show()
'''
low_freq_MEL=0 #将频率转换为梅尔刻度
nfilt=40 #窗的数目
#计算m=2595*log10(1+f/700)
high_freq_mel=(2595*np.log10(1+(sample_rate/2)/700))
mel_points=np.linspace(low_freq_MEL,high_freq_mel,nfilt+2) #梅尔刻度的均匀分布
#计算f=700(10**(m/2595)-1)
hz_points=(700*(10**(mel_points/2595)-1))
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]) #三角形滤波器左边频率f(m-1)
    f_m=int(bin[m]) #三角形滤波器中间频率fm
    f_m_plus=int(bin[m+1]) #三角形滤波器右边频率f(m-1)
    for k in range(f_m_minus,f_m):
        fbank[m-1,k]=(k-bin[m-1])/(bin[m+1]-bin[m])
plt.plot(fbank.T)
plt.show()
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)
num_ceps=12 #取12个系数
#通过DCT计算MFCC系数
mfcc=dct(filter_banks,type=2,axis=1,norm='ortho')[:,1:(num_ceps+1)]
#对MFCC进行倒谱提升可以改善噪声信号中的语音识别
(nframes,ncoeff)=mfcc.shape
n=np.arange(ncoeff)
cep_lifter=22 #倒谱滤波系数,定义倒谱所用到的滤波器组内滤波器个数
lift=1+(cep_lifter/2)*np.sin(np.pi*n/cep_lifter)
mfcc*=lift
mfcc-=(np.mean(mfcc,axis=0)+1e-8)
plt.imshow(np.flipud(mfcc.T),cmap=plt.cm.jet,aspect=0.2,extent=[0,mfcc.shape[0],0,mfcc.shape[1]]) #绘制MFCC热力图
plt.show()


52f10299a78b594c644b224a8f24f575_aca983c4e4b74dbcb3b370f1d5d6796a.png


对MFCC进行如下的归一化操作,运行操作,其相应的热力图如下:


filter_banks-=(np.mean(filter_banks,axis=0)+1e-8)
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.show()


4763d55458d74bba963e811be78131f0_2f370267378a41bb8d73c4b0ff31a2d5.png

                                归一化的MFCC热力图

目录
打赏
0
0
0
0
12
分享
相关文章
【MATLAB】语音信号识别与处理:移动中位数滤波算法去噪及谱相减算法呈现频谱
【MATLAB】语音信号识别与处理:移动中位数滤波算法去噪及谱相减算法呈现频谱
222 2
【MATLAB】语音信号识别与处理:一维信号NLM非局部均值滤波算法去噪及谱相减算法呈现频谱
【MATLAB】语音信号识别与处理:一维信号NLM非局部均值滤波算法去噪及谱相减算法呈现频谱
265 1
基于Logistic-Map混沌序列的数字信息加解密算法matlab仿真,支持对文字,灰度图,彩色图,语音进行加解密
本项目实现了一种基于Logistic Map混沌序列的数字信息加解密算法,使用MATLAB2022A开发并包含GUI操作界面。支持对文字、灰度图像、彩色图像和语音信号进行加密与解密处理。核心程序通过调整Logistic Map的参数生成伪随机密钥序列,确保加密的安全性。混沌系统的不可预测性和对初值的敏感依赖性是该算法的核心优势。示例展示了彩色图像、灰度图像、语音信号及文字信息的加解密效果,运行结果清晰准确,且完整程序输出无水印。
基于Logistic-Map混沌序列的数字信息加解密算法matlab仿真,支持对文字,灰度图,彩色图,语音进行加解密
|
12月前
|
语音去噪 python 传统算法
7月更文挑战第14天
217 1
基于语音信号MFCC特征提取和GRNN神经网络的人员身份检测算法matlab仿真
**语音识别算法概览** MATLAB2022a中实现,结合MFCC与GRNN技术进行说话人身份检测。MFCC利用人耳感知特性提取语音频谱特征,GRNN作为非线性映射工具,擅长序列学习,确保高效识别。预加重、分帧、加窗、FFT、滤波器组、IDCT构成MFCC步骤,GRNN以其快速学习与鲁棒性处理不稳定数据。适用于多种领域。
基于WOA鲸鱼优化的BiLSTM双向长短期记忆网络序列预测算法matlab仿真,对比BiLSTM和LSTM
本项目基于MATLAB 2022a/2024b实现,采用WOA优化的BiLSTM算法进行序列预测。核心代码包含完整中文注释与操作视频,展示从参数优化到模型训练、预测的全流程。BiLSTM通过前向与后向LSTM结合,有效捕捉序列前后文信息,解决传统RNN梯度消失问题。WOA优化超参数(如学习率、隐藏层神经元数),提升模型性能,避免局部最优解。附有运行效果图预览,最终输出预测值与实际值对比,RMSE评估精度。适合研究时序数据分析与深度学习优化的开发者参考。
基于GA遗传优化的BiLSTM双向长短期记忆网络序列预测算法matlab仿真,对比BiLSTM和LSTM
本内容包含基于BiLSTM与遗传算法(GA)的算法介绍及实现。算法通过MATLAB2022a/2024b运行,核心为优化BiLSTM超参数(如学习率、神经元数量),提升预测性能。LSTM解决传统RNN梯度问题,捕捉长期依赖;BiLSTM双向处理序列,融合前文后文信息,适合全局信息任务。附完整代码(含注释)、操作视频及无水印运行效果预览,适用于股票预测等场景,精度优于单向LSTM。
基于PSO粒子群优化的BiLSTM双向长短期记忆网络序列预测算法matlab仿真,对比BiLSTM和LSTM
本项目基于MATLAB2022a/2024b开发,结合粒子群优化(PSO)算法与双向长短期记忆网络(BiLSTM),用于优化序列预测任务中的模型参数。核心代码包含详细中文注释及操作视频,涵盖遗传算法优化过程、BiLSTM网络构建、训练及预测分析。通过PSO优化BiLSTM的超参数(如学习率、隐藏层神经元数等),显著提升模型捕捉长期依赖关系和上下文信息的能力,适用于气象、交通流量等场景。附有运行效果图预览,展示适应度值、RMSE变化及预测结果对比,验证方法有效性。
基于遗传算法的256QAM星座图的最优概率整形matlab仿真,对比优化前后整形星座图和误码率
本内容展示了基于GA(遗传算法)优化的256QAM概率星座整形(PCS)技术的研究与实现。通过Matlab仿真,分析了优化前后星座图和误码率(BER)的变化。256QAM采用非均匀概率分布(Maxwell-Boltzman分布)降低外圈星座点出现频率,减小平均功率并增加最小欧氏距离,从而提升传输性能。GA算法以BER为适应度函数,搜索最优整形参数v,显著降低误码率。核心程序实现了GA优化过程,包括种群初始化、选择、交叉、变异等步骤,并绘制了优化曲线。此研究有助于提高频谱效率和传输灵活性,适用于不同信道环境。
45 10
基于遗传优化ELM网络的时间序列预测算法matlab仿真
本项目实现了一种基于遗传算法优化的极限学习机(GA-ELM)网络时间序列预测方法。通过对比传统ELM与GA-ELM,验证了参数优化对非线性时间序列预测精度的提升效果。核心程序利用MATLAB 2022A完成,采用遗传算法全局搜索最优权重与偏置,结合ELM快速训练特性,显著提高模型稳定性与准确性。实验结果展示了GA-ELM在复杂数据中的优越表现,误差明显降低。此方法适用于金融、气象等领域的时间序列预测任务。

热门文章

最新文章

AI助理

你好,我是AI助理

可以解答问题、推荐解决方案等