MFCC语音特征值提取算法

博主简介

博主是一名大二学生，主攻人工智能研究。感谢让我们在CSDN相遇，博主致力于在这里分享关于人工智能，c++，Python，爬虫等方面知识的分享。 如果有需要的小伙伴可以关注博主，博主会继续更新的，如果有错误之处，大家可以指正。

专栏简介：   本专栏主要研究python在人工智能方面的应用，涉及算法，案例实践。包括一些常用的数据处理算法，也会介绍很多的Python第三方库。如果需要，点击这里   订阅专栏。

给大家分享一个我很喜欢的一句话：“每天多努力一点，不为别的，只为日后，能够多一些选择，选择舒心的日子，选择自己喜欢的人！”

目录

 背景引入

MFCC语音特征值提取算法简介

语音信号分帧

计算MFCC系数

对语音信号进行预加重

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

定义滤波器组

计算MFCC系数

 背景引入

特征值提取，在模式识别领域是很常见的一种算法和手段。特征值看起来好像很陌生，其实在我们日常生活中也很常见。我们使用的身份认证，ID，都可以视为不同系统下的特征值。

MFCC在语音识别领域就是一组特征向量，它通过对语音信号（频谱包络与细节）进行编码运算来得到。MFCC有39个系数，其中包括13个静态系数，13个一阶差分系数，以及13个二阶差分系数。差分系数用来描述动态特征，也就是声学特征在相邻帧间的变化情况。这些系数都是通过离散余弦变换（Discrete Cosine Transform，DCT）计算而来。

MFCC语音特征值提取算法简介

MFCC意为梅尔频率倒谱系数，顾名思义，MFCC语音特征提取包含两个关键步骤；将语音信号转化为梅尔频率，然后进行倒谱分析。梅尔频谱是一个可用来代表短期音频的频谱，梅尔刻度（Mel Scale）则是一种基于人耳对等距的音高变化的感官判断而确定的非线性频率刻度。梅尔频率和正常的频率f之间的关系：

mel(f)=25951g(1+f/700)

当梅尔刻度均匀分布，则对应的频率之间的距离会越来越大。梅尔刻度的滤波器组在低频部分的分辨率高，跟人耳的听觉特性比较相符，这也是梅尔刻度的物理意义。在梅尔频域内，人对音调的感知度为线性关系，如果两段语音的梅尔频率相差两倍，则人耳听起来两者的音调也相差两倍。

转化为梅尔频率时，首先对时域信号进行离散傅里叶变换，将信号转换到频域，然后再利用梅尔刻度的滤波器组对频域信号进行切分，使每个频率段对应一个数值。倒谱（Cepstrum）通过对一个时域信号进行傅里叶变换后取对数，并再次进行反傅里叶变换（Inverse Fast Fourier Transform，IFFT）得到。倒谱可分为复倒谱（Complex Cepstrum），实倒谱（Real Cepstrum）和功率倒谱（Power Cepstrum）。倒谱分析可用于信号分解，也就是将两个信号的卷积转化为两个信号的相加。

MFCC的物理含义，简而言之，可理解为语音信号的能量在不同频率范围的分布。

人的发声过程可以看作是肺里的气流通过声带这个线性系统。如果用e（t）表示输入声音的音高，h(t)表示声带的响应（也即我们需要获取的语音特征），那么听到的语音信号x(t)即为二者的卷积：

x(t)=e(t)*h(t)

x（t)为时域信号，对其进行离散傅里叶变换后可得到频域信号X（K），亦即频谱：

X(K)=DFT(x(t))

时域信号的卷积在频域内则可表示为二者的乘积：

X(K)=E(k)*H(k)

通常，在频域分析中我们只关注频谱的能量而忽略其相位信息，即：

||X(k)||=||E(k)||||H(k)||

对频谱进行对数运算：

log||X(k)||=log||E(k)||+log||H(k)||

然后进行反傅里叶计算：

c(n)=IFFT(log||X(k)||)=IFFT(log||E(k)||)+IFFT(log||H(k)||))

c(n)即为倒谱系数，已经和原始的时域信号x(t)不一样。并且时域信号的卷积关系已经转化为频域信号的线性相加关系。

语音信号分帧

语音信号属于准稳态信号，这也意味着，在一定的时间内，信号会保持稳定。这个时常对于我们人类来说，很短，一般只有10ms~30ms。在这一区间（即帧）内，可将语音信号看成稳态信号，只有稳态信号才能进行信号处理。

信号分帧一般会涉及到一个加窗的操作，即将原始信号与一个窗函数相乘。我们用计算机处理信号的时候，一般不会取无限长的信号，而是会取其中间的一段信号，这将会减少工作量，也会加快程序法分析的时间。

无限长的信号被截断后，其频谱会发生畸变，从而导致频谱能量泄露。 为了减少这种能量泄露，我们可采用不同的截取函数对信号进行截断。执行截断操作的函数称为窗函数，简称为窗。常用的窗函数有矩形窗，三角窗，汉明（Hamming）窗及汉宁窗等。

汉宁窗也叫升余弦窗，是很有用的窗函数。如果测试测试信号有多个频率分量，频谱表现非常复杂，测试目的更多在于关注频率点而非能量大小，则用汉宁窗。汉宁窗主瓣加宽并降低，旁瓣则显著减小，从减少泄漏的观点出发，汉宁窗明显优于矩形窗。但汉宁窗主瓣加宽，相当于分析带宽加宽，频率分辨率下降，他与矩形窗相比，泄露以及波动较小，选择性则相应较高。

汉明窗是用来加权余弦形成的锥形窗，也称之为改进的升余弦窗，只是加权系数不同，其旁瓣更小，但其旁瓣衰减速度比汉宁窗要慢。汉明窗是以著名的美国数学家理查德·卫斯理·汉明（Richard Wesley Hamming）的名字来命名：

w(n)=0.54-0.46cos(\frac{2 \pi n }{M-1}) (0\leq n\leqslant M-1)

下面的代码就是用python来生成汉明窗和汉宁窗：

import matplotlib.pyplot as plt
import scipy #信号处理工具包
plt.figure(figsize=(6,2))
plt.plot(scipy.hanning(512),"b--",label="Hanning") #绘制汉宁窗
plt.plot(scipy.hamming(512),"r--",label="Hamming") #绘制汉明窗
plt.title("Demo Hanning & Hamming Window")
plt.legend()
plt.show()

除了scipy模块可以实现汉宁窗和汉明窗，我们也可以用NumPy来实现汉宁窗和汉明窗。示例代码如下：


import numpy as np
import matplotlib.pyplot as plt
hanWing=np.hanning(512)#定义汉宁窗
hamWin=np.hamming(512) #定义汉明窗
plt.plot(hanWing,'y--',label="Hanning")
plt.plot(hamWin,'b--',label="Hamming")
plt.title("Hamming & Hanning window")
plt.ylabel("Amplitude")
plt.xlabel("Sample")
plt.legend()
plt.show()

信号加窗，从本质上而言，就是将原始信号与一个窗函数相乘。进行加窗操作之后，我们就可以对信号进行傅里叶展开。加窗的代价就是，一帧信号的两端部分将会被消弱。所以在进行信号分帧处理时，帧与帧之间需要有部分重叠。相邻两帧重叠后，其起始位置的时间差称之为帧移，即步长（Stride） 。

以下为简单的信号加窗操作示意图：


import numpy as np
import matplotlib.pyplot as plt
import scipy
x=np.linspace(0,10,1000)
originWav=np.sin(x**2) #示例原信号
win=scipy.hamming(1000) #定义一个窗函数，这里使用的汉明窗
winFrame=originWav*win
#结果可视化
plt.title("Signal Chunk with Hamming Windows")
plt.plot(originWav)
plt.plot(win)
plt.plot(winFrame)
plt.legend()
plt.show()

运行程序，其中蓝色波形为原信号，橙色波形为窗函数，绿色为加窗操作之后的信号。

假设x为语音信号，w为窗函数，则分帧信号为：

y(n)=\sum_{n=-(N/2)+1}^{N/2}x(m)w(n-m)

其中，w(n-m)为窗口序列，当n去不同的值时，窗口w(n-m)沿x(m)。因此，w(n-m)是一个“滑动的”窗口。y(n)为短时傅里叶变换（SIFT）。由于窗口是有限长度的，满足绝对可和条件，所以这个变幻的前提条件是存在的，这也是信号分帧的理论依据。

以下示例代码从指定文件夹读取一个音频文件，然后将该音频文件分帧并显示其中一个分帧信号的波形：


#读取指定音频文件
import matplotlib.pyplot as plt
import numpy as np
import wave #导入波形处理工具包
import os
import soundfile
def audioSignalFrame(signal,nw,inc):
    '''
    signal:原始音频信号
    nw:每一帧的长度
    inc:相邻帧的间隔
    '''
    #信号总长度
    signal_length=len(signal)
    #若信号长度小于一个帧的长度，则帧数定义为1
    if signal_length<=nw:
        nf=1
    else:
        nf=int(np.ceil((1.0*signal_length-nw+inc)/inc))
        #所有帧加起来的总的铺平的长度
    pad_length=int((nf-1)*inc+nw)
    #长度不够时，使用0填补，类似于FFT中的扩充数组长度
    zeros=np.zeros((pad_length-signal_length,))
    #填补后的信号
    pad_signal=np.concatenate((signal,zeros))
    #相当于对所有帧的时间点进行抽取，得到nf*nw的长度的矩阵
    indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T
    #将indices转化为矩阵
    indices=np.array(indices,dtype=np.int32)
    #得到帧信号
    frames=pad_signal[indices]
    #窗函数，这里默认取1
    return frames
def readSignalWave(filename):
    f=wave.open(filename,'rb')
    params=f.getparams()
    nchannels,sampwidth,frammerate,nframes=params[:4]
    #读取音频，字符串格式
    strData=f.readframes(nframes)
    #将字符串转化为int
    waveData=np.fromstring(strData,dtype=np.int16)
    f.close()
    #信号幅值归一化
    waveData=waveData*1.0/(max(abs(waveData)))
    waveData=np.reshape(waveData,[nframes,nchannels]).T
    return waveData
if __name__=='__main__':
    filepath="./test.wav"
    #dirname=os.listdir(filepath)
    #filename=filepath+dirname[3]
    data=readSignalWave(filepath)
    #初始化每帧长度及帧间隔
    nw=512
    inc=128
    Frame=audioSignalFrame(data[0],nw,inc)
    #显示原始信号
    plt.plot(data[0])
    plt.title("Original Signal")
    plt.show()

    #显示第一帧信号
    plt.plot(Frame[0])
    plt.title("First Frame")
    plt.show()

上面的代码中，没有对信号进行加窗处理，若要执行信号加床操作，只需将分帧函数稍作修改，


def audioSignalFrame(signal,nw,inc,winfunc):
    '''
    signal:原始音频信号
    nw:每一帧的长度
    inc:相邻帧的间隔
    '''
    #信号总长度
    signal_length=len(signal)
    #若信号长度小于一个帧的长度，则帧数定义为1
    if signal_length<=nw:
        nf=1
    else:
        nf=int(np.ceil((1.0*signal_length-nw+inc)/inc))
        #所有帧加起来的总的铺平的长度
    pad_length=int((nf-1)*inc+nw)
    #长度不够时，使用0填补，类似于FFT中的扩充数组长度
    zeros=np.zeros((pad_length-signal_length,))
    #填补后的信号
    pad_signal=np.concatenate((signal,zeros))
    #相当于对所有帧的时间点进行抽取，得到nf*nw的长度的矩阵
    indices=np.tile(np.arange(0,nw),(nf,1))+np.tile(np.arange(0,nf*inc,inc),(nw,1)).T
    #将indices转化为矩阵
    indices=np.array(indices,dtype=np.int32)
    #得到帧信号
    frames=pad_signal[indices]
    #窗函数，这里默认取1
    win=np.tile(winfunc,(nf,1))

    return frames*win

当然，随着函数的改变，主函数中对函数的调用也需要改变，只需要改变参数即可。除了调用工具包中的汉明窗函数，也可以使用公式来定义。

def hamming(n):

    return 0.54-0.46*cos(2*pi/n*(arange(n)+0.5))

 语音信号在进行分帧之前，一般需要进行一个与加重操作。语音信号的预加重，是为了对语音的高频部分进行加重，使信号变得平坦，保持在地频到高频的整个频带中能用同样的信噪比求频谱。同时也为了消除发声过程中声带和口唇辐射效应，补偿语音信号受到发音系统所抑制的高频部分，增加语音的高频分辨率。

我们一般通过一阶有限长单位冲激响应（Finite Impulse Response，FIR）高通数字滤波器来实现预加重。FIR滤波器H(z)=1-az^{-1}作为传递函数，其中a为预加重系数，0.9<a<1.0。假设t时刻的语音采样值为x(t)，经过预加重处理后的结果为y(t)=x(t)-ax(t-1).a一般默认取0.95.

信号的预加重处理示例代码：

def preemphasis(signal,coeff=0.95):
    '''
    signal:要滤波的输入信号
    coeff:预加重系数。0表示无过滤，默认为0.95
    返回值：滤波信号
    '''
    return numpy.append(signal[0],signal[1:]-coeff*signal[:-1])

计算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+\frac{f}{700})

def hz2mel(hz):

    return 2595*numpy.log10(1+hz/700.0)

同样，我们也可以推出下列公式：

 f=700(10^{m/2595}-1)

 Delta系数的计算公式为：

d_{t}=\frac{\sum_{n=1}^{N}n(c_{t+n}-c_{t-n})}{2\sum_{n=1}^{N}n^{2}}

其中，d_{t}为Delta系数，从帧t根据静态系数c_{t-N}到c_{t+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。

对语音信号进行预加重


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()

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

 

 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()

                                                        (a)幅度谱

                                                             (b)功率谱

 

                                                             （c)功率对数谱

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

 

定义滤波器组

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

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

H_m(k)=\left\{\begin{matrix} \frac{2(k-f(m-1))}{(f(m+1)-f(m-1)(f(m)-f(m-1)))} &f(m-1)\leq k\leq f(m+1) \\ \frac{2(f(m+1)-k)}{(f(m+1)-f(m-1)(f(m+1)-f(m-1)))}& f(m)\leq k\leq f(m+1) \end{matrix}\right.

对于其他的情况，例如，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()

三角形滤波器有两个主要功能，其一，对频谱进行平滑并消除谐波的作用，突显原先语音的共振峰；其二，用以降低运算量。如图所示的滤波器组中的每个滤波器在中心频率处响应为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()

 

 

对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()

                                 归一化的MFCC热力图

好了，本篇文章介绍就到此结束了，拜了个拜！ 

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：/a/5435.html

如若内容造成侵权/违法违规/事实不符，请联系我们进行投诉反馈qq邮箱809451989@qq.com，一经查实，立即删除！

相关文章

银行数字化转型导师坚鹏：银行产品经理技能快速提升之道

银行数字化转型导师坚鹏：银行产品经理技能快速提升之道

银行产品经理技能快速提升之道 ——以推动银行战略目标实现为核心&#xff0c;实现知行果合一课程背景&#xff1a; 很多银行都在开展产品经理技能提升工作&#xff0c;目前存在以下问题急需解决&#xff1a; 不知道银行产品经理掌握哪些关键知识&#xff1f; 不清楚如何有效…
阅读更多...
C++修炼之练气期第二层——缺省参数

C++修炼之练气期第二层——缺省参数

目录 1.缺省参数的概念 2.缺省参数的分类 全缺省参数 半缺省参数 实用场景示例 1.缺省参数的概念 缺省参数是声明或定义函数时为函数的参数指定一个缺省值。 在调用该函数时&#xff0c;如果没有指定实参则采用该形参的缺省值&#xff0c;否则使用指定的实参。 #inclu…
阅读更多...
四大排序算法之归并排序

四大排序算法之归并排序

说明 为了自己学习方便&#xff0c;我这里总结了四大排序算法涵盖了七种排序算法 分类算法名称时间复杂度 空间复杂度稳定性插入排序 直接插入排序 希尔排序 O(n^2) O(1) O(n^2/3) O(1) 稳定 不稳定 选择排序 选择排序 堆排序 O(n^2) O(1) O(nlogn) O(1) 不稳…
阅读更多...
linux查看进程、端口

linux查看进程、端口

1、先查看进程pidps -ef | grep 进程名如果已知pid&#xff0c;想看详情&#xff0c;则用 ps -ef pid2、通过pid查看占用端口(mac)netstat -na | grep 端口netstat -nap tcp | grep 进程pidnetstat -nap udp | grep 进程pid不加tcp或者udp的话mac上会报错&#xff1a;netstat常…
阅读更多...
基于ASP的反垃圾邮件管理系统的设计与实现

基于ASP的反垃圾邮件管理系统的设计与实现

随着Internet的迅速普及&#xff0c;电子邮件以其快捷、方便、低成本的特点逐渐成为人们进行信息交流的主要媒介之一&#xff0c;但是随之而来的垃圾邮件也越来越泛滥。垃圾邮件占用了有限的存储、计算和网络资源&#xff0c;耗费了用户大量的处理时间&#xff0c;影响和干扰了…
阅读更多...
程序员OKR学习法

程序员OKR学习法

版权声明 本文原创作者&#xff1a;谷哥的小弟作者博客地址&#xff1a;http://blog.csdn.net/lfdfhl OKR管理法 OKR&#xff08;Objectives and Key Results&#xff09;管理法是一种目标管理方法&#xff0c;旨在通过制定明确的目标和可量化的关键结果来帮助组织、团队和个人…
阅读更多...
RocketMQ的架构图

RocketMQ的架构图

文章目录RocketMQ 技术架构中有四大角色 NameServer 、Broker 、Producer 、Consumer 。我来向大家分别解释一下这四个角色是干啥的。 Broker&#xff1a; 主要负责消息的存储、投递和查询以及服务高可用保证。说白了就是消息队列服务器嘛&#xff0c;生产者生产消息到 Broker…
阅读更多...
Hive实战 --- 电子商务消费行为分析

Hive实战 --- 电子商务消费行为分析

目录 数据结构 Customer表 Transaction表 Store表 Review表 上传数据 创建目录用于存放数据 把本地文件上传到HDFS上 创建外部表 创建数据库 创建表 数据清洗 对transaction_details中的重复数据生成新ID 过滤掉store_review中没有评分的数据 找出PII (personal …
阅读更多...
【web前端初级课程】第八章 什么是事件？

【web前端初级课程】第八章 什么是事件？

目录 一、事件情况汇总 二、标签绑定 三、使用DOM0事件模型 四、使用DOM2事件模型 五、相关练习&#xff1a;图片切换 一、事件情况汇总 事件分为三部分&#xff1a;事件源&#xff1a;绑定事件的标签、事件对象&#xff1a;就是事件产生的相关数据、事件处理函数 二、标…
阅读更多...
Java使用功能方法交换a，b的值，通过构造方法输出姓名、年龄、家庭地址

Java使用功能方法交换a，b的值，通过构造方法输出姓名、年龄、家庭地址

目录 前言 一、使用功能方法交换a&#xff0c;b的值 1.1运行流程&#xff08;思想&#xff09; 1.2代码段 1.3运行截图 二、通过构造方法输出姓名、年龄、家庭地址 1.1运行流程&#xff08;思想&#xff09; 1.2代码段 1.3运行截图 前言 1.因多重原因&#xff0c;所以我…
阅读更多...
愚人节，聊聊那些正在坑人的“新型AI”

愚人节，聊聊那些正在坑人的“新型AI”

几年前的一个愚人节&#xff0c;我们和大家聊过AI技术被作为诈骗工具的情况。很不幸&#xff0c;当时讨论的一些苗头&#xff0c;现在都成了电诈犯罪中屡见不鲜的手段。更可气的是&#xff0c;随着AI技术与应用本身的发展&#xff0c;犯罪分子的AI手段不减反增。一些“新型AI”…
阅读更多...
（七）Tomcat源码阅读：Host组件分析

（七）Tomcat源码阅读：Host组件分析

一、概述 Host类中比较重要的类就是HostConfig其它类实现的功能和之前的组件差不多&#xff0c;这里就不多介绍了。 二、阅读源码 1、HostConfig &#xff08;1&#xff09;重要方法 lifecycleEvent&#xff1a; 根据对应的方法设置对应的属性&#xff0c;并调用对应的方…
阅读更多...
自己写gpt的软件教程-国内最好的chatgpt软件

自己写gpt的软件教程-国内最好的chatgpt软件

GPT-3是一种非常强大的自然语言处理技术&#xff0c;可以为用户生成高质量的文本内容。虽然GPT-3最初是为英文而设计的&#xff0c;但是近年来&#xff0c;GPT-3在中文领域也变得越来越流行。在本篇教程中&#xff0c;我们将详细介绍如何在GPT-3中生成中文内容。 一、准备工作 …
阅读更多...
第二天并发篇

第二天并发篇

一、线程状态 1.新建&#xff08;New&#xff09;&#xff1a;创建线程对象时 2.就绪&#xff08;Runnable&#xff09;&#xff1a;线程调用start方法&#xff0c;有执行资格没有执行权 3.运行&#xff1a;当就绪状态时抢到cpu的执行权之后&#xff0c;进入运行状态 4.阻塞&am…
阅读更多...
过程控制系统中的模块技术MTP

过程控制系统中的模块技术MTP

在过程自动化行业中&#xff0c;模块化设备概念近年来越来越受欢迎。其中最热门的是MTP。MTP称为模块类型封装&#xff0c;它是过程工业自动化技术用户协会&#xff08;NAMUR&#xff09;提出的过程自动化行业的模块化标准&#xff0c;通过这种模型&#xff0c;开发工作的重点从…
阅读更多...
C++(Qt)软件调试---linux下生成/调试Core文件（3）

C++(Qt)软件调试---linux下生成/调试Core文件（3）

#软件调试 C(Qt)软件调试—linux下生成/调试Core文件&#xff08;3&#xff09; 文章目录C(Qt)软件调试---linux下生成/调试Core文件&#xff08;3&#xff09;前言1、C生成Core和使用GDB调试1、环境2、C生成Core文件3、使用gdb工具调试core可定位段错误位置&#xff1b;4、修…
阅读更多...
【创作赢红包】你是真的“C”——C语言中文件操作函数使用的详细讲解【上篇】

【创作赢红包】你是真的“C”——C语言中文件操作函数使用的详细讲解【上篇】

你是真的“c”——C语言中文件操作函数使用的详细讲解~&#x1f60e;前言&#x1f64c;一、 为什么使用文件&#xff1a;&#x1f64c;二、 什么是文件&#xff1a;&#x1f64c;2.1 程序文件2.2 数据文件2.3 文件名3. 文件的打开和关闭3.1 文件指针3.2 文件的打开和关闭4. 文件…
阅读更多...
【ansible】实施任务控制

【ansible】实施任务控制

目录 实施任务控制 一&#xff0c;循环&#xff08;迭代&#xff09;--- loop 1&#xff0c;利用loop----item循环迭代任务 2&#xff0c;item---loop循环案例 1&#xff0c;定义item循环列表 2&#xff0c;通过变量应用列表格式 3&#xff0c;字典列表&#xff08;迭代嵌套子…
阅读更多...
一个ESP32小东西

一个ESP32小东西

之前发了ESP8266&#xff0c;有人评论说玩下ESP32然后就买了几个回来&#xff0c;当然&#xff0c;也想着和大家一起玩介绍下这个开发板开发板Github项目链接https://github.com/Xinyuan-LilyGO/T-QT把仓库的代码下载到本地我们可以用ESP-IDF和Arduino两个SDK来开发ESP32S3ESP-…
阅读更多...
回溯算法思想、回溯算法解题模板与回溯算法题目索引（不断更新）

回溯算法思想、回溯算法解题模板与回溯算法题目索引（不断更新）

回溯算法 回溯算法是一种试探性的搜索算法&#xff0c;它在解决某些组合问题、优化问题、路径问题等&#xff0c;非常有效。回溯算法的核心思想是通过递归和深度优先搜索&#xff08;DFS&#xff09;来搜索问题的解空间。 细说一下这些问题&#xff1a; 组合问题&#xff1a;N…
阅读更多...
最新文章

Copyright @ 2022~2023