MFCC语音特征值提取算法

news2024/11/24 12:51:49

博主简介

博主是一名大二学生，主攻人工智能研究。感谢让我们在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热力图

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

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

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

Java学习之路003——集合

Java学习之路003——集合

1、 代码演示 【1】新增一个类&#xff0c;用来测试集合。先创建一组数组&#xff0c;数组可以存放不同的数据类型。对于Object类型的数组元素&#xff0c;可以通过.getClass方法获取到具体类型。 【2】如果数组指定类型为int的时候&#xff0c;使用.getClass()就会提示错误。 …
阅读更多...
官方不推荐@Autowired

官方不推荐@Autowired

1用lombok注解 2 构造器
阅读更多...
云原生是什么？核心概念和应用方法解析

云原生是什么？核心概念和应用方法解析

什么是云原生&#xff1f; 云原生是一种基于容器、微服务和自动化运维的软件开发和部署方法。它可以使应用程序更加高效、可靠和可扩展&#xff0c;适用于各种不同的云平台。 如果要更直接通俗的来解释下上面的概念。云原生更准确来说就是一种文化&#xff0c;是一种潮流&…
阅读更多...
使用邮箱验证实现登录功能（发送邮件，redis）

使用邮箱验证实现登录功能（发送邮件，redis）

目录 概述 前端搭建 后端搭建 生成验证码-存入redis&#xff08;主要过程代码&#xff09; 发送邮件&#xff08;主要过程代码&#xff09; 登录验证-取出redis中数据进行验证&#xff08;主要代码&#xff09; 完整代码一-LoginController 完整代码二-LoginService 完…
阅读更多...
域基础和基本环境搭建

域基础和基本环境搭建

1.1 名词解释 域和工作组的区别&#xff1a; 工作组中所有的计算机都是对等的&#xff0c;也就是没有服务器和客户机的之分&#xff0c;所以工作组并不存在真正的集中管理作用&#xff1b;域是一个有安全边界的计算机集合&#xff0c;安全边界指的是一个域中的用户无法访问到另…
阅读更多...
六、并发集合

六、并发集合

文章目录并发集合ConcurrentHashMap存储结构存储操作put方法putVal方法-散列算法putVal方法-添加数据到数组&初始化数组putVal方法-添加数据到链表扩容操作treeifyBin方法触发扩容tryPreSize方法-针对putAll的初始化操作tryPreSize方法-计算扩容戳并且查看BUGtryPreSize方法…
阅读更多...
[git可视化软件]gitkraken平替:GitAhead

[git可视化软件]gitkraken平替:GitAhead

日期2023-02-28 gitkraken6.5.1已经不能登陆使用了!! 6.5.1免费版已经无法使用!!! 现在是2023-02-28 这款工具已经废除了6.5.1版本的使用功能了,我直接卡在登陆界面进不去项目了. 要想继续管理私有项目,只能升级最新版的软件,并且开通会员.会员费用高的一批,一年要59.4美元.约…
阅读更多...
前端进阶JS运行原理

前端进阶JS运行原理

JS运行原理 深入了解V8引擎原理 浏览器内核是由两部分组成的&#xff0c;以webkit为例&#xff1a; WebCore&#xff1a;负责HTML解析、布局、渲染等等相关的工作&#xff1b;JavaScriptCore&#xff1a;解析、执行JavaScript代码&#xff1b; 官方对V8引擎的定义&#xff1…
阅读更多...
docker部署zabbix6.2.7+grafana

docker部署zabbix6.2.7+grafana

目录 1、下载docker 2、下载相关镜像文件 3、创建一个供zabbix系统使用的网络环境 4、创建一个供mysql数据库存放文件的目录 5、启动mysql容器 6、为zabbix-server创建一个持久卷 7、启动zabbix-server容器 8、创建语言存放目录 9、启动zabbix-web容器 10、启动zabbix…
阅读更多...
探析集团企业 1+N 模式，重新定义集团型CRM

探析集团企业 1+N 模式，重新定义集团型CRM

目录 一、客户经营、运营监控 二、流程驱动、业务成长 三、规则规范 业务治理 什么是集团型CRM【1N】&#xff1f;本文中我们可以把集团看作为“1”&#xff0c;其他分公司或组织看作为“N”。本篇我们主要分析集团CRM业务定位。 我们从企业集团总部的职能定位确定集团CRM…
阅读更多...
攻不下dfs不参加比赛（七）

攻不下dfs不参加比赛（七）

标题 为什么练dfs题目总结重点为什么练dfs 相信学过数据结构的朋友都知道dfs（深度优先搜索）是里面相当重要的一种搜索算法，可能直接说大家感受不到有条件的大家可以去看看一些算法比赛。这些比赛中每一届或多或少都会牵扯到dfs，可能提到dfs大家都知道但是我们为了避免眼高手…
阅读更多...
数据结构(Java)---队列

数据结构(Java)---队列

1.简介 队列是一种操作受限的线性表&#xff0c;和栈类似&#xff0c;队列的限制主要是运行在表的一端进行插入操作&#xff0c;而在表的另一端进行删除的操作&#xff08;栈的操作受限是只能在一端进行插入和删除&#xff09;。队列中进行插入的一端称为队尾 &#xff0c;删除…
阅读更多...
【JavaWeb】Servlet基础

【JavaWeb】Servlet基础

文章目录1.Tomcat服务器安装注意事项2.编写WebApp3.BS系统角色和协议4.模拟Servlet4.1模拟sun公司4.2模拟Tomcat服务器4.3模拟WebApp开发者5.开发一个带有Servlet的WebApp5.1创建一个名为crm的项目5.2 在项目中创建一个名为WEB-INF的文件&#xff08;必须&#xff09;5.3在WEB-…
阅读更多...
JAVA开发（java类加载过程）

JAVA开发（java类加载过程）

1、java语言的平台无关性。 因为java语言可以跑在java虚拟机上&#xff0c;所以只要能装java虚拟机的地方就能跑java程序。java语言以后缀名 .java为文件扩展名。通过java编译器javac编译成字节码文件.class 。java字节码文件通过java虚拟机解析运行。所以java语言可以说是编译…
阅读更多...
解决windows安装wxPython安装失败、速度过慢及PyCharm上wx包爆红问题

解决windows安装wxPython安装失败、速度过慢及PyCharm上wx包爆红问题

网上关于wxPython安装失败&#xff0c;安装速度过慢&#xff0c;以及安装成功后PyCharm中import wx仍然爆红的文章有很多&#xff0c;也特别杂&#xff0c;解决起来特别困难&#xff0c;今天在这里对问题的处理进行一个整合&#xff0c;希望能帮助到大家。 安装wxPython这里运用…
阅读更多...
模拟微信聊天-课后程序(JAVA基础案例教程-黑马程序员编著-第八章-课后作业)

模拟微信聊天-课后程序(JAVA基础案例教程-黑马程序员编著-第八章-课后作业)

【案例9-1】 模拟微信聊天 【案例介绍】 1.案例描述 在如今&#xff0c;微信聊天已经人们生活中必不可少的重要组成部分&#xff0c;人们的交流很多都是通过微信来进行的。本案例要求&#xff1a;将多线程与UDP通信相关知识结合&#xff0c;模拟实现微信聊天小程序。通过监…
阅读更多...
华为OD机试题【翻转单词顺序】用 C++ 进行编码 （2023.Q1）

华为OD机试题【翻转单词顺序】用 C++ 进行编码 （2023.Q1）

最近更新的博客 华为od 2023 | 什么是华为od，od 薪资待遇，od机试题清单） 文章目录 最近更新的博客使用说明翻转单词顺序题目输入输出示例一输入输出示例二输入输出说明示例三输入输出说明示例四
阅读更多...
解决跑pytorch代码报错AttributeError: module ‘distutils‘ has no attribute ‘version‘

解决跑pytorch代码报错AttributeError: module ‘distutils‘ has no attribute ‘version‘

跑pytorch代码报错AttributeError: module ‘distutils’ has no attribute ‘version’ Traceback (most recent call last): File “D:/pycharm_envir/gaozhiyuan/Segmentation/pytorch_segmentation/deeplabv3-plus-pytorch-main/train.py”, line 16, in from utils.callb…
阅读更多...
【C++修炼之路】22.哈希

【C++修炼之路】22.哈希

每一个不曾起舞的日子都是对生命的辜负 哈希一.哈希概念及性质1.1 哈希概念1.2 哈希冲突1.3 哈希函数二.哈希冲突解决2.1 闭散列/开放定址法2.2 开散列/哈希桶三.开放定址法代码3.1 插入Insert3.2 查找Find3.3 删除Erase3.4 映射的改良&完整代码四.开散列代码4.1 插入Inser…
阅读更多...
Linux之init.d、rc.d文件夹说明

Linux之init.d、rc.d文件夹说明

备注&#xff1a;Ubuntu没有rc.d文件夹&#xff0c;原因看问题四 Linux的几个重要文件 rc.d&#xff0c;init.d文件夹的说明 今天在研究mysql的安装的时候&#xff0c;最后一步要创建一个软连接&#xff0c;使得mysql服务可以自启动&#xff0c;代码如下&#xff1a; ln -s…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023