精华内容
下载资源
问答
  • 主要为大家详细介绍了梅尔倒谱系数(MFCC)实现,具有一定的参考价值,感兴趣的小伙伴们可以参考一下
  • MFCC梅尔频率倒谱系数

    2021-01-30 14:43:56
    MFCC梅尔频率倒谱系数 文章目录MFCC梅尔频率倒谱系数MFCC梅尔频率倒谱系数的原理计算MFCC的步骤1、音频信号文件读取及预加重2、信号分帧加窗3、计算帧功率谱的周期图估计4、将梅尔滤波器组应用于功率谱,将每个...

    MFCC梅尔频率倒谱系数原理及python实现

    MFCC梅尔频率倒谱系数的原理

    MFCC梅尔频率倒谱系数的应用:MFCC通过提取语音信号特征(识别音频信号,丢弃背景噪声、情绪信息)准确表示短时功率谱图,进一步表现语音声道形状,有助于我们准确表示声道形状发音产生的音素,被广泛的应用于以线性预测系数(LPC)、线性预测倒谱系数(LPCC)(倒谱和LPCC教程)以及隐马尔可夫模型HMM分类器主要特征类型的自动语音识别(ASR)和说话人识别。

    计算MFCC的步骤

    1. 音频信号文件读取及预加重。

    2. 信号分帧加窗。

    3. 计算帧功率谱的周期图估计

    4. 将梅尔滤波器组应用于功率谱,将每个滤波器的能量求和。

    5. 取所有滤波器组能量的对数。

    6. 取对数滤波器组能量的离散余弦变换DCT。

    7. 保持DCT系数2-13,其余部分丢弃。

      梅尔音阶

      Volkmann和Newmann提出了一个音高单位,以使相等的音高距离听起来与听众相等, 这称为梅尔音阶。 对频率执行数学运算,以将其转换为mel标度。

      梅尔音阶将感觉到的纯音的频率或音高与其实际测量的频率相关联。与低频时相比,人类在分辨低频时音调的细微变化方面要好得多。合并此比例使我们的功能与人类听到的声音更加匹配。从频率转换为梅尔标度的公式为:
      在这里插入图片描述
      从梅尔回到频率:
      在这里插入图片描述
      对语音信号进行16kHz采样。
      1125取值有时也取2595或者其他值是两个不同的度量标准。

    1、音频信号文件读取及预加重

    • 依据所选取的音频文件,采样频率一般为8000Hz和16000Hz,需大于真实信号最大频率的2倍。
    • 预加重处理其实是将语音信号通过一个高通滤波器:
      在这里插入图片描述
      预加重的目的是提升高频部分,使信号的频谱变得平坦,保持在低频到高频的整个频带中,能用同样的信噪比求频谱,避免在后边的FFT操作中出现数值问题。同时,也是为了消除发生过程中声带和嘴唇的效应,来补偿语音信号受到发音系统所抑制的高频部分,也为了突出高频的共振峰。其中其中α一般取值为0.97、0.95。

    使用librosa.get_duration读取信号的持续时间,采样librosa.load返回语音信号采样率与音频采样数组,采样率需要大于真实信号的频率的2倍才可无失真还原。

    在预加重避免在FFT操作中出现数值问题,需要对高频信息进行加强,因为一般高频能量比低频小,取系数不同加强的程度也不同。
    在这里插入图片描述

    #导入音频处理的python工具包
    import numpy as np
    import librosa
    import time
    import math
    from scipy.fftpack import fft, ifft
    #导入音频及绘图显示包
    import librosa.display
    #导入绘图工作的函数集合
    import matplotlib.pyplot as plt
    #读取音频文件,Librosa默认的采样率是22050Hz
    #如果需要读取原始采样率44100Hz,需要设定参数sr=None
    #如果需要重采样,只需要将采样率参数sr设定为你需要的值
    times= librosa.get_duration(filename='./DrinkWater.wav')#获取音频时长单位为秒
    y, sr = librosa.load('DrinkWater.wav', sr=16000,offset=0.0,duration=None)#返回音频采样数组及采样率
    print("sr=",sr)
    print("times=",times)
    #绘图
    #plt.figure()
    #显示波形图
    #librosa.display.waveplot(y, sr, x_axis='time')#画出波形的振幅线
    #plt.title("DrinkWater wave")
    PointNumbers=int(times*sr)+1#在此由于时间精确度较高,总取样点数取整后加一,保证不丢失信息量,代价是多处理一个采样点
    #在此抛出一个时间精确度的问题
    print("PointNumbers=",PointNumbers)
    x1= np.arange(0,PointNumbers,1)#采样点刻度
    x2=np.arange(0,times,1/sr)#时间刻度
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('DrinkWater.wav', fontsize=12, color='black')
    plt.plot(x1,y)
    plt.show()
    plt.figure()
    plt.xlabel("times")
    plt.ylabel("wave")
    plt.title('DrinkWater.wav', fontsize=12, color='black')
    plt.plot(x2,y)
    plt.show()
    

    在此可根据语音处理的需求设置音频读取的时间精度来确定所取的采样点的个数。
    求得:

    sr= 16000#采样频率
    times= 26.992426303854874#音频时间
    PointNumbers= 431879#采样点数
    

    注意到返回的y为长度为431879的array数组。

    array([0.60138094, 0.9744118 , 0.8517753 , ..., 0.8879393 , 0.8933856 , 0.        ], dtype=float32)#
    

    音频DrinkWater.wav的采样图(采样点刻度):
    在这里插入图片描述
    音频DrinkWater.wav的采样图(采样时间刻度):
    在这里插入图片描述
    在此由于时间精确度较高,以sr=16000求出的总采样点数PointNumbers不是整数数,可以对求出总采样点数直接取整或是取整后后加1,二者不同之处在于,取整加一后保证不丢失信息量,代价是多处理一个采样点。在图中末尾的拖尾0点就是就是这个采样点带来的效果。

    从上面的分析,在此可在采样时间的精度上进行优化,由采样率为16000可知,依据采样频率可取时间精度为千分之一来保证采样点个数为整数,而不用增加或舍去采样点。

    #预加重x=y(i)-0.97y(i-1)
    PreEmphasis=y
    PointNumbers=int(PointNumbers)#转化为整型
    for i in range(1,PointNumbers,1):#range(PointNumbers),PointNumbers需为整型
        PreEmphasis[i]=PreEmphasis[i]-0.97*PreEmphasis[i-1]
    #plt.figure()
    #显示预加重波形图
    #print("音频文件drink_water.wav预加重波形图:")
    #librosa.display.waveplot(PreEmphasis, sr)
    #plt.title("PreEmphasis RecordVower1 wave")
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('PreEmphasis DrinkWater.wav', fontsize=12, color='black')
    plt.plot(x1,PreEmphasis)
    plt.show()
    plt.figure()
    plt.xlabel("time")
    plt.ylabel("wave")
    plt.title('PreEmphasis DrinkWater.wav', fontsize=12, color='black')
    plt.plot(x2,PreEmphasis)
    plt.show()
    

    音频DrinkWater.wav的预加重采样图(采样点刻度):
    在这里插入图片描述
    音频DrinkWater.wav的预加重采样图(时间刻度):
    在这里插入图片描述
    由图可知预加重后波形的幅值减小了。

    2、信号分帧加窗

    音频信号是非平稳信号,在短时范围内(20-40 ms)可假设语音信号是统计平稳的,语音信号分帧时间范围也为(20-40 ms),选取帧长太短时,我们将没有足够的样本来获得可靠的频谱估计;选取帧长太长时,则信号在整个帧中变化太大。

    • 25 ms是标准配置,故将语音信号分成25 ms每帧,例如:16kHz信号的帧长为0.025 * 16000 = 400个样本,帧步长通常约为10毫秒(160个样本)(帧会有一些重叠)前400个样本帧从样本0开始,下一个400样本帧从样本从160开始,直到到达语音文件的末尾,如果语音文件没有分成偶数个帧,则用零填充赝本集使帧个数成偶数个。
      帧长=fs.0.025(采样点)帧移=fs.0.01(采样点) \text{帧长}=f_s.0.025\left( \text{采样点} \right) \\ \text{帧移}=f_s.0.01\left( \text{采样点} \right)

    • 对于每个单个帧,为每个帧提取一组12个MFCC系数。简而言之,就是符号:S(n)。装帧后,Si(n)可以确定n的范围是1-400(如果我们的帧是400个样本),i范围是帧数。当我们计算复数DFT时,我们得到Si(k)其中i表示与时域帧相对应的帧号。Pi(k)是i的功率谱。

    • 要获取帧的离散傅立叶变换:
      在这里插入图片描述
      其中h(n)是N样本长序列的分析窗口(例如汉明窗口)。


    #分帧
    LenFrame=int(sr*0.025)#每帧长(持续0.025秒)
    PatFrame=int(sr*0.001)#帧移(持续0.001秒)
    Discover=int(LenFrame-PatFrame)#每帧之间的非重叠部分
    NumFrames=(PointNumbers-LenFrame)/Discover
    #分成整数帧
    if (NumFrames-int(NumFrames))==0 :
        NumFrames=int(NumFrames)
    else:
        NumFrames=int(NumFrames)+1#帧号数
    #if (NumFrames%2)==0:
    #    NumFrames=int(NumFrames)
    #else:
    #   NumFrames=int(NumFrames)+1
    NumFillZero=(NumFrames*Discover+LenFrame)-PointNumbers#填充0的数目
    NumFrames=int(NumFrames)+1#在此表示帧数
    print("sr=",sr)
    print("times=",times)
    print("LenFrame=",LenFrame)
    print("PatFrame=",PatFrame)
    print("Discover=",Discover)
    print("NumFrames=",NumFrames)
    print("PointNumbers=",PointNumbers)
    print("NumFillZero=",NumFillZero)
    print("sumPointNumbers=",PointNumbers+NumFillZero)
    FillZeros=np.zeros(NumFillZero)#生成填充序列
    PreEmphasis1=np.concatenate((PreEmphasis,FillZeros)) #填补后的信号记为PreEmphasis1,np.concatenate()连接两个维度相同的矩阵
    #plt.figure()
    #显示预加重填充波形图
    #print("音频文件drink_water.wav预加重填充波形图:")
    #librosa.display.waveplot(PreEmphasis1, sr)
    #plt.title("PreEmphasis1 RecordVower1 wave")
    x3= np.arange(0,PointNumbers+NumFillZero,1)#填充后采样点刻度
    x4= np.arange(0,(PointNumbers+NumFillZero)/sr,1/sr)#填充后时间刻度
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('PreEmphasis1 DrinkWater.wav', fontsize=12, color='black')
    plt.plot(x3,PreEmphasis1)
    plt.show()
    plt.figure()
    plt.xlabel("times")
    plt.ylabel("wave")
    plt.title('PreEmphasis1 DrinkWater.wav', fontsize=12, color='black')
    plt.plot(x4,PreEmphasis1)
    plt.show()
    

    结果得:

    sr= 16000
    times= 26.992426303854874
    LenFrame= 400
    PatFrame= 16
    Discover= 384
    NumFrames= 1125
    PointNumbers= 431879
    NumFillZero= 137
    sumPointNumbers= 432016
    

    可知采样数据点为PointNumbers有431879个,填充零为NumFillZero有137个,此时的总数据点数为sumPointNumbers=PointNumbers+NumFillZero有432016个。

    音频DrinkWater.wav的预加重填充采样图(采样点刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样图(时间刻度):

    在这里插入图片描述

    #开始分帧
    indices=np.tile(np.arange(0,LenFrame),(NumFrames,1))+np.tile(np.arange(0,NumFrames*Discover,Discover),(LenFrame,1)).T  #相当于对所有帧的时间点进行抽取,得到NumFrames*LenFrame长度的矩阵
    indices=np.array(indices,dtype=np.int32) #将indices转化为矩阵
    Frames=PreEmphasis1[indices] #得到帧信号
    x5= np.arange(PointNumbers+NumFillZero-LenFrame,PointNumbers+NumFillZero,1)#最后一帧采样点刻度
    x6= np.arange((PointNumbers+NumFillZero-LenFrame)/sr,(PointNumbers+NumFillZero)/sr-1/sr,1/sr)#最后一帧时间刻度
    #x9=np.arange((times+NumFillZero/sr)-(LenFrame/sr),times+NumFillZero/sr,1/sr)#某帧时间刻度
    plt.figure()
    #显示最后一帧波形图
    #print("音频文件drink_water.wav预加重填充最后一帧波形图:")
    #librosa.display.waveplot(Frames[NumFrames-1], sr)
    #plt.title("PreEmphasis1 RecordVower1 LastFrames  wave")
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('PreEmphasis1 DrinkWater LastFrames  wave', fontsize=12, color='black')
    plt.plot(x5,Frames[NumFrames-1])
    plt.show()
    plt.figure()
    plt.xlabel("times")
    plt.ylabel("wave")
    plt.title('PreEmphasis1 DrinkWater LastFrames  wave', fontsize=12, color='black')
    plt.plot(x6,Frames[NumFrames-1])
    plt.show()
    x7= np.arange(0,LenFrame,1)#第一帧采样点刻度
    x8= np.arange(0,(LenFrame)/sr,1/sr)#第一帧时间刻度
    #显示第一帧波形图
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('PreEmphasis1 DrinkWater FirstFrames  wave', fontsize=12, color='black')
    plt.plot(x7,Frames[0])
    plt.show()
    plt.figure()
    plt.xlabel("times")
    plt.ylabel("wave")
    plt.title('PreEmphasis1 DrinkWater FirstFrames  wave', fontsize=12, color='black')
    plt.plot(x8,Frames[0])
    plt.show()
    #显示第二帧波形图
    x9= np.arange(LenFrame,2*LenFrame,1)#第二帧采样点刻度
    x10= np.arange(LenFrame/sr,2*LenFrame/sr,1/sr)#第二帧时间刻度
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('PreEmphasis1 DrinkWater SecondFrames  wave', fontsize=12, color='black')
    plt.plot(x9,Frames[1])
    plt.show()
    plt.figure()
    plt.xlabel("times")
    plt.ylabel("wave")
    plt.title('PreEmphasis1 DrinkWater SecondFrames  wave', fontsize=12, color='black')
    plt.plot(x10,Frames[1])
    plt.show()
    

    分帧后的Frames为一个1125*400的多维数组(1125维,每维代表一帧)。

    结果:

    音频DrinkWater.wav的预加重填充采样数据分帧后的最后一帧图(采样点刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样数据分帧后的最后一帧图(时间刻度):

    在这里插入图片描述

    可在图中26.990到26.995秒之间看到一个负数的点,这个点是第一节抛出的问题引起的,可改变下面语句参数:

    times= librosa.get_duration(filename='./DrinkWater.wav')#获取音频时长单位为秒
    y, sr = librosa.load('DrinkWater.wav', sr=16000,offset=0.0,duration=None)#返回音频采样数组及采样率
    

    将第一句中的times的精度对应采样频率保留三位小数,将第二句中duration=None改为duration=times(此times为三位小数精度的时间)使填充前最后一个点的取值不为0。

    音频DrinkWater.wav的预加重填充采样数据分帧后的第一帧图(采样点刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样数据分帧后的第一帧图(时间刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样数据分帧后的第二帧图(采样点刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样数据分帧后的第二帧图(时间刻度):

    在这里插入图片描述

    #加窗
    HanningWindow=np.hanning(LenFrame)  #调用汉明窗参数为帧长
    #a=0.46
    #myhanming=np.array([0.]*LenFrame)  #创建一个1行LenFrame列的一维零数组
    #for i in range(0,LenFrame,1):    
    #    myhanming[i]=(1-a)-a*(math.cos(2*(math.pi)*i/LenFrame))
    plt.figure()
    #显示汉明窗
    #print("显示汉明窗:")
    #librosa.display.waveplot(HanningWindow, sr)
    #plt.title("HanningWindow(librosa.display.waveplot)  wave")
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('HanningWindow  wave)', fontsize=12, color='black')
    plt.plot(x7,HanningWindow)
    plt.figure()
    plt.xlabel("time")
    plt.ylabel("wave")
    plt.title('HanningWindow  wave', fontsize=12, color='black')
    plt.plot(x8,HanningWindow)
    #print(HanningWindow)
    HanningWindows=np.tile(HanningWindow,(NumFrames,1))#将窗函数平铺成NumFrames行
    #print(HanningWindows)
    #分别对每一帧加窗
    Frames=np.multiply((np.mat(Frames)),(np.mat(HanningWindows)))
    Frames=np.array(Frames)
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('SecondFrames with HanningWindow wave', fontsize=12, color='black')
    plt.plot(x9,Frames[1])
    plt.show()
    plt.figure()
    plt.xlabel("times")
    plt.ylabel("wave")
    plt.title('SecondFrames with HanningWindow wave', fontsize=12, color='black')
    plt.plot(x10,Frames[1])
    plt.show()
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("wave")
    plt.title('LastFrames with HanningWindow wave', fontsize=12, color='black')
    plt.plot(x5,Frames[NumFrames-1])
    plt.show()
    plt.show()
    plt.figure()
    plt.xlabel("times")
    plt.ylabel("wave")
    plt.title('LastFrames with HanningWindow wave', fontsize=12, color='black')
    plt.plot(x6,Frames[NumFrames-1])
    plt.show()
    

    汉明窗公式为:
    W(n)=(1a)acos(2πnN)  1nN W\left( n \right) =\left( 1-a \right) -a\cdot \cos \left( 2\cdot \pi \cdot \frac{n}{N} \right) \,\, 1\geqslant n\geqslant N

    a的取值不同对应的窗口也不同,一般情况下a = 0.46 。

    结果:

    汉明窗图(采样点刻度):

    在这里插入图片描述

    汉明窗图(时间刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样数据分帧加窗后的第二帧图(时间刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样数据分帧加窗后的第二帧图(时间刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样数据分帧加窗后的最后一帧图(采样点刻度):

    在这里插入图片描述

    音频DrinkWater.wav的预加重填充采样数据分帧加窗后的最后一帧图(时间刻度):

    在这里插入图片描述

    由加窗后的图像可知,经过加窗处理后,音频部分被抑制,第二帧加窗后体现不出信号的幅值变化,抛出一个问题。

    另选一个信号验证本程序加窗前后波形如下:

    在这里插入图片描述

    在这里插入图片描述

    得出出现上面抛出的问题的原因可能是加窗后整体帧信号的局部变化幅值较小造成的。

    3、计算帧功率谱的周期图估计

    计算每个帧的功率谱,周期图估计值根据声音的频率在不同的位置振动情况(对应人体耳蜗中振动的位置-摆动小毛发),确定帧中存在哪些频率。

    • K是DFT的长度。语音帧的基于周期图的功率谱估计Pi(K)由下式给出:
      在这里插入图片描述
      这称为功率谱的周期图估计。N为窗口长度,我们取复数傅里叶变换的绝对值,并将结果平方。

    • 计算滤波器组的梅尔间隔,取20-40(26个标准)三角形滤波器应用于周期图功率谱估计,滤波器组以26个长度为N的向量的形式出现,每个向量大多为零,但在频谱的特定部分非零。


    #帧功率谱的周期图
    #求傅里叶变换
    #x1 = np.linspace(0, 1, LenFrame)#创建0-1的等差数组 LenFrame个
    x2 = np.arange(0,LenFrame,1)#创建0-LenFrame的等差数组,步长为1
    FftFrames = np.fft.fft(Frames)#fft计算
    AbsFrames = np.abs(FftFrames)#求模
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("amplitude")
    plt.title('SecondFrames fft', fontsize=12, color='black')
    plt.plot(x9,AbsFrames[1])
    plt.show()
    #功率谱的周期图估计
    Power=(np.square(AbsFrames))/LenFrame
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("amplitude")
    plt.title('Power', fontsize=12, color='black')
    plt.plot(Power)
    plt.show()
    

    对于加窗后的矩阵Frames,它是一个1125**400的矩阵,也就是说,它有1125帧数据,且每一帧数据都有400个采样点,那么我们接下来就要对这1125帧的每一帧都要进行快速傅里叶变换,得到一个大小为1125400大小的矩阵FftFrames,其帧数还是1125帧,对每一帧的400个数据点分别取模再取平方,然后除以LenFrame为400;便得到能量谱密度Power,其大小为1125x400,对每一帧得到的能量进行相加,代表每一帧能量的总和。

    结果:

    Power的维数为(1125*400)

    显示第二帧的傅里叶变换:

    在这里插入图片描述

    显示所有帧的功率谱的周期图估计:

    在这里插入图片描述

    4、将梅尔滤波器组应用于功率谱,将每个滤波器的能量求和

    周期图频谱估计包含许多自动语音识别(ASR)不需要的信息,表现为耳蜗不能辨别两个间隔较小的信号频率之间的差异,随着频率的增加,这种影响变得更加明显,此时需要将成簇的周期图进行汇总,以了解各个频率区域中存在多少能量。

    梅尔滤波器组执行周期图能量汇总:第一个滤波器非常窄,可以指示0赫兹附近存在多少能量。随着频率的升高,我们对滤波器的关注也越来越小,滤波器也变得越来越宽。对每个位置产生多少能量,梅尔刻度可以准确地告诉我们,间隔滤光片组以及滤光片组的宽度下面给出了有关如何计算间距的信息。
    在这里插入图片描述
    梅尔滤波器组和开窗功率谱图

    • 为了计算滤波器组的能量,将每个滤波器组乘以功率谱,然后将系数相加。完成此操作后,我们剩下26个数字,这些数字可以指示每个滤波器组中的能量。

    • 实际上常常使用26-40个滤波器组。

    • 为了获得上图(a)所示的滤波器组,选择一个较低和较高的频率,较低的频率为300Hz,较高的频率为8000Hz。如果语音以8000Hz采样,则上限频率将限制为4000Hz。

    • 使用梅尔音阶公式将高频和低频转换为梅尔兹。在此,300Hz为401.25 Mels,而8000Hz为2834.99 Mels。

    • 在此示例中,做10个滤波器组,为此需要12个点。需要在401.25和2834.99之间线性间隔的10个额外点。

      结果是:

    m(i)= 401.25、622.50、843.75、1065.00、1286.25、1507.50、1728.74, 
           1949.99、217.24、239.49、2613.74、2834.99
    
    • 使用梅尔音阶公式将其转换回赫兹:
    h(i)= 300,517.33,781.90,1103.97,1496.04,1973.32,2554.33, 
           3261.62,4122.63,5170.76,6446.70,8000
    

    注意:频率的起点和终点都在我们想要的频率上。

    • 我们没有将滤波器放置在上面计算出的精确点上所需的频率分辨率,因此我们需要将这些频率四舍五入到最接近的FFT档。此过程不会影响功能的准确性。要将频率转换为FFT bin(离散线谱间隔),FFT大小和采样率如下:
    f(i)=下限((nfft + 1)* h(i)/采样率)
    

    结果如下:

    f(i)= 9、16、25、35、47、63、81、104、132、165、206、256
    

    可以看到最终的滤波器组在bin 256处结束,这对应于512点FFT大小的8kHz。

    • 创建滤波器组:第一个滤波器组将从第一个点开始,在第二个点达到峰值,然后在第三个点返回零。第二个滤波器组将从第二个点开始,在第三个点达到最大值,然后在第四个点为零,依此类推。计算公式如下:
      在这里插入图片描述
      其中M是滤波器数量,f()是M +2梅尔间隔频率列表。

    • 彼此重叠的所有10个过滤器的最终图是:
      在这里插入图片描述

    上图包含了10个过滤器的梅尔过滤器组,该滤波器组以0Hz开始,以8000Hz结束。这仅是一个指南,实际工作示例始于300Hz。

    • 滤波器的第二种方法:创建Hm矩阵,这个矩阵得定义如下所示:
      在这里插入图片描述
      发现该梅尔滤波器,在频率很小的时候,滤波器宽度很窄,随着其频率的增大,滤波器的宽度也会逐渐增大,但其幅值也会逐渐减小。
      在这里插入图片描述

    #梅尔滤波器组应用于功率谱
    fl = 0
    fh =int(sr)
    bl = 1125*np.log(1+fl/700) # 把 Hz 变成 Mel
    bh = 1125*np.log(1+fh/700)
    #梅尔频率转换
    x11 = np.arange(fl,fh,1)#创建0-LenFrame的等差数组,步长为1
    MelFrequency= 1125*np.log(1+x11/700)
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("mel-frequency")
    plt.title('mel-frequency transfornation', fontsize=12, color='black')
    plt.plot(x11,MelFrequency)
    plt.show()
    

    其在低频范围内增长速度很快,但在高频范围内,梅尔值的增长速度很慢。每一个频率值都对应着一个梅尔值。

    结果:

    梅尔频率图:

    在这里插入图片描述

    #梅尔滤波器组一
    p = 26#滤波器个数
    B = bh-bl
    y = np.linspace(0,B,p+2)# 将梅尔刻度等间隔
    #print(y)
    Fb = 700*(np.exp(y/1125)-1)# 把 Mel 变成 Hz
    #print(Fb)
    W2 = int(LenFrame )#LenFrame内对应的FFT点数
    df = sr/LenFrame
    freq = []#采样频率值
    for n in range(0,W2):
        freqs = int(n*df)
        freq.append(freqs)#采样频率值
    bank = np.zeros((p,W2))#生成一个p行W2列的全零数组
    for k in range(1,p+1):
        f1 = Fb[k-1]
        f2 = Fb[k+1]
        f0 = Fb[k]
        n1=np.floor(f1/df)#f(m-1)在频域中的谱线索引号
        n2=np.floor(f2/df)#f(m+1)在频域中的谱线索引号
        n0=np.floor(f0/df)#f(m)在频域中的谱线索引号。f(m)是从0开始,而在MATLAB中数组的索引是从1开始,所以要加1,否则会出现index=0的错误
        for i in range(1,W2):
            if i>=n1 and i<=n0:
                bank[k-1,i]=(i-n1)/(n0-n1)
            elif i>n0 and i<=n2:
                bank[k-1,i]=(n2-i)/(n2-n0)
        # print(k)
        # print(bank[k-1,:])
        plt.plot(freq,bank[k-1,:])
    plt.show()
    #将梅尔滤波器组和功率谱相乘
    #for k in range(1,p+1):
    #H=Power*(bank.T)
    #H=mat(Power)*mat((bank.T))#H(NumFrames*p)
    H=np.mat(Power)*np.mat((bank.T))#H(NumFrames*p)
    #print(Power)#Power(NumFrames*LenFrame)
    #print(bank.T)#bank(p*LenFrame)
    #print(H)
    #梅尔滤波器组二
    p2 = 26#滤波器个数
    B2 = bh-bl
    y2 = np.linspace(0,B2,p2+2)# 将梅尔刻度等间隔
    #print(y2)
    Fb2 = 700*(np.exp(y2/1125)-1)# 把 Mel 变成 Hz
    #print(Fb2)
    W22 = int(LenFrame )#LenFrame内对应的FFT点数
    df2 = sr/LenFrame
    freq2 = []#采样频率值
    for n2 in range(0,W22):
        freqs2 = int(n2*df2)
        freq2.append(freqs2)#采样频率值
    bank2 = np.zeros((p2,W22))#生成一个p行W2列的全零数组
    for k2 in range(1,p2+1):
        f12 = Fb2[k2-1]
        f22 = Fb2[k2+1]
        f02 = Fb2[k2]
        n12=np.floor(f12/df2)#f(m-1)在频域中的谱线索引号
        n2=np.floor(f22/df2)#f(m+1)在频域中的谱线索引号
        n02=np.floor(f02/df2)#f(m)在频域中的谱线索引号。f(m)是从0开始,而在MATLAB中数组的索引是从1开始,所以要加1,否则会出现index=0的错误
        for i2 in range(1,W22):
            if i2>=n12 and i2<=n02:
                bank2[k2-1,i2]=(2*(i2-n12))/((n2-n12)*(n02-n12))
            elif i2>n02 and i2<=n2:
                bank2[k2-1,i2]=(2*(n2-i2))/((n2-n12)*(n02-n12))
        # print(k2)
        # print(bank2[k2-1,:])
        plt.plot(freq2,bank2[k2-1,:])
    plt.show()
    

    结果:

    梅尔滤波器组一:

    在这里插入图片描述

    梅尔滤波器组二:

    在这里插入图片描述

    抛出问题,不同的梅尔滤波器组,不同滤波器参数及作用。

    梅尔滤波器组的作用是放大低频信息,抑制高频频信息,尽可能保留低频部分信息,低频信息的跨度小于高频部分,幅度更体现低频的放大与高频的缩小 。

    5、取所有滤波器组能量的对数

    对数运算允许我们使用倒谱均值减法,这是一种通道归一化技术,人类听力听不到线性范围的响度,有了滤波器组能量,我们就可以取它们的对数,要将声音的感知音量加倍,我们需要将成(>8)倍的能量投入其中,这意味着,如果声音很大,开始时能量的大变化听起来可能并没有那么大的不同,这种(压缩)操作使我们的对数处理功能与人类实际听到的声音更加接近。

    • 取指示滤波器能量的26个数字的对数,得到26个log滤波器组能量。

    #将梅尔滤波器组和功率谱相乘
    #将梅尔滤波器组和功率谱相乘
    H2=np.mat(Power)*np.mat((bank2.T))#H2(NumFrames*p2)
    #取所有滤波器组能量的对数
    H3=np.log(H2)#第二个滤波器
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("mfcc")
    plt.title('mfcc-filter*powe', fontsize=12, color='black')
    plt.plot(H2)
    plt.show()
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("mfcc")
    plt.title('mfcc-log(filter*powe)', fontsize=12, color='black')
    plt.plot(H3)
    plt.show()
    

    结果:

    H2、H3的维数为:1125*26

    梅尔滤波器组和功率谱相乘:

    在这里插入图片描述

    所有滤波器组能量的对数图:

    在这里插入图片描述

    6、取对数滤波器组能量的离散余弦变换DCT

    离散余弦变换具有信号谱分量丰富,能量集中,不需要对语音进行相位估算等优点,在较低的计算复杂度下获得较好的语音增强的效果

    滤波器组都是重叠的,故滤波器组的能量彼此相关性极强,DCT对能量进行去相关,这意味着可以使用对角协方差矩阵对HMM分类器中的特征进行建模。

    • .对26个对数滤波器组能量进行离散余弦变换(DCT),得到26个倒谱系数。

    • MFCC(i,n)=m=1Mlog[H(i,m)].cos[π.n.(2m1)2M] MFCC\left( i,n \right) =\sum_{m=1}^M{\log \left[ H\left( i,m \right) \right] .\cos \left[ \frac{\pi .n.\left( 2m-1 \right)}{2M} \right]}

      其中H为我们上边得到的矩阵H,M代表梅尔滤波器个数,i代表第几帧数据(取值为1-NumFrames),n代表第i帧的第n列(n取值范围为1-26)。


    #print(mfcc)#完成离散余弦变换
    #在Numpy库中,函数的输入x不仅可以为单独一个数,还可以是一个列表,一个Numpy数组
    #其结果为Numpy数组。也就是说Numpy中的函数不需要循环就可以实现每个元素的批量处理。
    mfcc=np.zeros(shape=(NumFrames,p2))#创建0数组
    for i3 in range(0,NumFrames-1):
        for i4 in range(0,p2-1):
            sum=0#求每帧能量总和
            for i5 in range(0,p2-1):#做离散余弦变换
                sum=sum+H3[i3,i5]*math.cos(((math.pi)*i4)*(2*i5-1)/(2*p2))
            mfcc[i3,i4]=((2/p2)**0.5)*sum
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("mfcc")
    plt.title('mfcc-DCT', fontsize=12, color='black')
    plt.plot(mfcc)
    plt.show()
    

    结果:

    得到1125*26的数组

    对数滤波器组能量的离散余弦变换DCT图:

    在这里插入图片描述

    7、保持DCT系数2-13,其余部分丢弃

    保留26个DCT系数中的12个,这是因为较高的DCT系数表示滤波器组能量的快速变化,并且事实证明,这些快速变化实际上会降低ASR性能,因此通过降低DCT系数可以得到很小的改进。

    • 仅保留26个系数中较低的12-13,产生的特征(每帧12个数字)称为“梅尔频率倒谱系数”。

    #保持DCT系数2-13,其余部分丢弃
    #mfcc=mfcc[[0,2],:] #取某几行
    #mfcc=mfcc[:,[end-12,end]] #取后13列
    mfcc=mfcc[:,[0,1,2,3,4,5,6,7,8,9,10,11,12]] #取前13列
    #mfcc3=mfcc[:,[0,1,2,3]] #取前13列
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("mfcc")
    plt.title('mfcc', fontsize=12, color='black')
    plt.plot(mfcc)
    plt.show()
    

    结果:

    得到1125*13的数组

    取前13列图:

    在这里插入图片描述

    得到1125*2的数组

    取2列图:

    在这里插入图片描述

    8、求升倒谱系数及mfcc的第一组、第二组及第三组

    因为大部分的信号数据一般集中在变换后的低频区,所以对每一帧只取前13个数据就好了。其公式表达如下:
    在这里插入图片描述

    其中L为升倒谱系数,一般取值为22,我们将其带入即可求得矩阵K,这是一个1x13大小的矩阵,这一部分的升倒谱的其实现代码如下:


    #求升倒谱取升倒谱系数为22
    L=22
    K=[0,0,0,0,0,0,0,0,0,0,0,0,0]
    for i in range(0,12,1):
        K[i]=1+(L/2)*math.sin((math.pi)*(i+1)/L)
    #print(K)
    #得到二维数组feat,这是mfcc的第一组数据,默认为三组
    feat=mfcc1
    for i3 in range(0,NumFrames-1):
        for i in range(0,12):
            feat[i3,i]= feat[i3,i]*K[i]
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("feat")
    plt.title('feat', fontsize=12, color='black')
    plt.plot(feat)
    plt.show()
    
    

    在这里插入图片描述

    包括第一组参数:
    在这里插入图片描述

    接下来就是求取MFCC的第二组,第三组参数了。第二组参数其实就是在已有的基础参数下作一阶微分操作,第三组参数在第二组参数下作一阶微分操作:

    在这里插入图片描述

    第二组参数:

    #接下来求取第二组(一阶差分系数),这也是mfcc参数的第二组参数
    dtfeat=np.array([[0]*13]*(NumFrames))  #创建一个NumFrames-1行13列的二维零数组
    #求差分
    for i3 in range(0+3,NumFrames-1-2):
        for i in range(0,12):
            dtfeat[i3,i]=feat[i3+1,i]-feat[i3-1,i]+2*feat[i3+2,i]-2*feat[i3-2,i]
    dtfeat=dtfeat/10
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("dtfeat")
    plt.title('dtfeat', fontsize=12, color='black')
    plt.plot(dtfeat)
    plt.show()
    

    在这里插入图片描述

    
    #求取二阶差分系数,mfcc参数的第三组参数
    #二阶差分系数就是对前面产生的一阶差分系数dtfeat再次进行操作
    dttfeat=np.array([[0]*13]*(NumFrames))  #创建一个NumFrames-1行13列的二维零数组
    #求差分
    for i3 in range(0+3,NumFrames-1-2):
        for i in range(0,12):
            dttfeat[i3,i]=dtfeat[i3+1,i]-dtfeat[i3-1,i]+2*dtfeat[i3+2,i]-2*dtfeat[i3-2,i]
    dttfeat=dttfeat/10#这里的10是根据数据确定的
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("dttfeat")
    plt.title('dttfeat', fontsize=12, color='black')
    plt.plot(dttfeat)
    plt.show()
    

    在这里插入图片描述

    9、将三组数据合并得到特征参数

    由于一阶求导和二阶求导的前两帧和后两帧数据都为0,于是我们就要对在mfcc_final中去除这四帧数据:

    #将三组参数拼接
    mfcc_final=(np.concatenate((feat.T,dtfeat.T,dttfeat.T))).T#1125*39 
    #填补后的信号记为PreEmphasis1,np.concatenate()连接几个个维度相同的矩阵
    plt.figure()
    plt.xlabel("frequency")
    plt.ylabel("mfcc_final")
    plt.title('mfcc_final', fontsize=12, color='black')
    plt.plot(mfcc_final)
    plt.show()
    #由于一阶求导和二阶求导的前两帧和后两帧数据都为0,于是我们就要对在mfcc_final中去除这四帧数据
    mfcc24=mfcc_final[3:NumFrames-1-2,:]
    #选取每一帧的mfcc系数的第一个数得到MFCC0
    mfcc0=mfcc24[:,1]#选取mfcc系数的第一个数,组成新的特征参数mfcc0
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("mfcc0")
    plt.title('mfcc0', fontsize=12, color='black')
    plt.plot(mfcc0)
    plt.show()
    mfcc00=(mfcc0-80)/2
    plt.figure()
    plt.xlabel("point")
    plt.ylabel("mfcc00")
    plt.title('mfcc00', fontsize=12, color='black')
    plt.plot(mfcc00)
    plt.show()
    
    

    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    总的有1125帧数据,剪掉四帧后 ,得到1121点的数据。
    可以看到,原始信号之前是431879个采样点的数据,而现在的M F C C 0 0参数图形大致与原始信号一致,并且其点数只有1121个点,这也就说明通过此方法M F C C 0 0,我们可以提取出语音信号的特点以及走向趋势,也就是说在某个程度上我们可以用这1121个点来代替431879点的数据。

    10、验证

    matlab验证:
    Matlab调用h=melbankm(p,n,fs,fl,fh,w),输入音频为DrinkWater.wav,数据每帧长n为400,滤波器组数p为26,采样频率fs为波器最16000、最低频为0,最高频fl取0.5,窗口采用汉明窗m(三角窗取t海宁窗取n)。fs与fl为采用fs归一化的频率。(问题:在书本中fl一般取0.5,对应例题的采样频率 为8000)f(i)=上限((nfft + 1)* h(i)/采样率)。
    在这里插入图片描述

    另一段语音验证:

    采用语音RecordVower1.wav
    在这里插入图片描述

    mfcc(matlab)
    在这里插入图片描述

    mfcc(mypython)
    在这里插入图片描述


    总结

    本次任务,也驱动于初学python,注释掉的代码部分是行不通的试错代码或不是最优需改进的代码,在此未作删除留作学习python的备忘,读者可忽略或指正。

    注意:对每个小节分割线上为理论资料,分割线下为具体实现,请注意资料对应的举例说明参数仅便于理解方法,具体音频实例参数对照注释及第二小节第一个图以及代码的语义注释。

    展开全文
  • 语音信号处理之(四)梅尔频率倒谱系数,语音信号处理之(四)梅尔频率倒谱系数语音信号处理之(四)梅尔频率倒谱系数语音信号处理之(四)梅尔频率倒谱系数语音信号处理之(四)梅尔频率倒谱系数语音信号处理之(四...
  • 从periodogram estimate of the power spectrum计算得到的倒谱系数,可以用于基音追踪(pitch tracking),然而,从AR power spectral estimate计算得到的倒谱系数可以用于语音识别(现在已经被MFCCs所替代)。...

    倒谱是表示一帧语音数据特征的一个序列。从periodogram estimate of the power spectrum计算得到的倒谱系数,可以用于基音追踪(pitch tracking),然而,从AR power spectral estimate计算得到的倒谱系数可以用于语音识别(现在已经被MFCCs所替代)。 One of the benefits of cepstrum and LPCCs over e.g. LPCs is that you can do cepstral mean subtraction (CMS) on cepstral coefficients to remove channel effects。

    倒谱(Cepstrum)是什么?

    这一节将会描述如何从功率谱的periodogram estimate来计算倒谱。首先,我们介绍自相关系数(autocorrelattion),然后,介绍如何使用类似的方法来计算倒谱。最后,我们介绍LPCCs的相关知识。

    倒谱可以看做和自相关序列类似的东西。如果,我们获得了信号的功率谱,我们可以使用Wiener-Khinchin theorem来计算信号的自相关序列。用数学公式表示为:$x(n)$是信号的时域表示,$X(k)$是信号的复数谱,$P(k)$是信号的功率谱,$A(n)$是信号的自相关序列。

    通过对信号$x(n)$取离散傅里叶变换,可以得到信号的复数谱: $$DFT(x(n)) {\rightarrow} X(k) $$

    对复数谱做逆离散傅里叶变换,可以得到信号的时域表示: $$IDFT(X(k)) \rightarrow x(n)$$

    通过对信号$x(n)$DFT的绝对值取平方,可以得到信号的功率谱: $$|DFT(x(n))|^2 \rightarrow P(k)$$

    如果,我们对信号的功率谱做$IDFT$,可以得到信号的自相关序列: $$IDFT(P(k)) \rightarrow A(n)$$

    如果,我们在对信号的功率谱做$IDFT$之前,对功率谱取对数,则可以获得信号的倒谱: $$IDFT(log(P(k))) \rightarrow C(n)$$

    因此,可以将倒谱理解为自相关序列的对数压缩,因为其携带有和自相关序列类似的信息(如“信号的周期性),但是,信号的倒谱是由其对数功率谱而非标准功率谱计算得到的。(注:有资料显示,对功率谱取对数,可以将乘性信号转变为加性信号。)

    13e35089a20331c1cab5241c85bbdb61.png

    在上述图片中,展示了原始语音信号的一帧数据、对数功率谱、自相关序列和倒谱。倒谱当中的峰值对应自相关序列当中的一个峰值,但是更为清晰。这个峰值的位置在58sample位置处,对应的基音频率为(16000/58)=275Hz (信号采样率为16kHz)。这是一个相当高的基音频率,原始语音信号是从一位女性说话人得到的。因为倒谱的强峰值,所以经常被用于基音检测。

    在MATLAB当中,我们可以使用下面代码得到信号的倒谱:

    PowerSpectrum = abs(fft(SpeechFrame,1024)).^2;

    AutoCorrelation = ifft(PowerSpectrum,1024);

    Cepstrum = ifft(log(PowerSpectrum),1024);

    什么是LPCCs(线性预测倒谱系数)?

    在前面的小节,我们介绍了标准倒谱(standard cepstrum),线性预测倒谱系数(Linear Prediction Cepstral Coefficients)的计算与其类似,除了LPCCs是由光滑自回归功率谱(smoothed Auto-Regressive power spectrum)计算得到的,而非功率谱的peridogram estimate。例:10阶AR谱估计,Levinson Curbin算法应用前10个自相关系数来计算10个线性预测系数。在MATLAB中,使用下面代码进行实现:

    [lp,g] = lpc(frame,10)

    前一节所展示信号的线性预测系数为:

    1.00 -2.22 1.68 0.05 -1.28 1.32 -0.30 -0.76 1.35 -1.19 0.44

    从AR谱估计来计算倒谱:

    [lp,g] = lpc(SpeechFrame,10);

    ARPowerSpectrum = g ./ abs(fft(lp,1024)).^2;

    Cepstrum = ifft(log(ARPowerSpectrum),1024);

    **注意到,上述计算步骤和前一节的计算过程完全相同,除了计算的功率谱不同。**前十个倒谱系数为:

    -11.78 2.22 0.79 -0.12 0.38 0.03 -0.20 0.04 -0.42 -0.11

    以这种方法计算的系数,第一个系数是被忽略的,仅仅决定于$g$。这是一帧数据的线性预测倒谱系数。如果想要计算信号的LPCCs,不要使用这种方法进行计算,因为下面有更加有效的方法。

    这一节的内容旨在展示倒谱系数和线性预测倒谱系数之间的区别和联系。

    从LPCs计算LPCCs

    这是一个简单的从LPCs计算LPCCs的递归公式,不用计算DFT。在下式当中,$a_k$表示线性预测系数,在前面小节中表示为$\text{lp}$。 $$ c(n)= \begin{cases} 0 &n<0\ ln(G) &n=0\ a_n+\sum{{k=1}^{n-1}(\frac{k}{n})c(k)a{n-k}} &0p\ \end{cases} $$

    从一个有限LPC系数,可以得到无限长度的线性预测倒谱系数。研究显示,12-20个倒谱系数对语音识别任务已经足够。

    参考文献

    展开全文
  • mel倒谱系数的提取

    2012-12-04 15:38:13
    用matlab编写 音频中mel倒谱系数的提取
  • 信号 基本Mel频率倒谱系数计算器
  • 主要为大家详细介绍了语音识别之梅尔频率倒谱系数及Python实现,具有一定的参考价值,感兴趣的小伙伴们可以参考一下
  • 加密域梅尔频率倒谱系数和脆弱的音频水印
  • 此文件夹应放在包含所有语音处理练习的同一文件夹中。
  • Mel倒谱系数

    2013-11-24 20:34:00
    Mel倒谱系数:MFCC Mel频率倒谱系数(Mel Frequency Cepstrum Coefficient)的缩写是MFCC,Mel频率是基于人耳听觉特性提出来的,它与Hz频率成非线性对应关系。Mel频率倒谱系数(MFCC)则是利用它们之间的这种关系,...

    Mel倒谱系数:MFCC

     

    Mel频率倒谱系数(Mel Frequency Cepstrum Coefficient)的缩写是MFCC,Mel频率是基于人耳听觉特性提出来的,它与Hz频率成非线性对应关系。Mel频率倒谱系数(MFCC)则是利用它们之间的这种关系,计算得到的Hz频谱特征。

     

    用录音设备录制一段模拟语音信号后,经由自定的取样频率(如8000 Hz、16000 Hz等)采样后转换(A/D)为数字语音信号。由于在时域(time domain)上语音信号的波形变化相当快速、不易观察,因此一般都会在频域(frequency domain)上来观察,其频谱是随着时间而缓慢变化的,因此通常可以假设在一较短时间中,其语音信号的特性是稳定的,通常我们定义这个较短时间为一帧(frame),根据人的语音的音调周期值的变化,一般取10~20ms。

     

    Mel-frequency cepstrum coefficient

    作用:和线性预测倒谱系数LPCC一起用于描述语音特征的参数:能量,基音频率,共振峰值等。

    详解几个概念:

    1.Mel频率:

    是模拟人耳对不同频率语音的感知。

    人类对不同频率语音有不同的感知能力:对1kHz以下,与频率成线性关系,对1kHz以上,与频率成对数关系。频率越高,感知能力就越差了。因此,在应用中常常只使用低频MFCC,而丢弃中高频MFCC。

    在Mel频域内,人对音调的感知能力为线性关系,如果两段语音的Mel频率差两倍,则人在感知上也差两倍。 转换公式:B(f)=1125ln(1+f/700) 其中f为频率,B为Mel-频率。

    2.倒谱:

    同态处理的结果,分为复数和实数倒谱,常用实数倒谱,是语音识别中的重要系数。

    具体过程:傅里叶变换----->对数运算----->傅里叶反变换。

    语音的产生用源、滤波器模型来表示,即把声带振动看作激励源e(n),把声道看成一个滤波器h(n),两者在时域进行卷积,得到语音信号s(n)。为了更好地处理语音,则需要分析s(n)以分别得到e(n)和h(n),这个过程称为解卷过程。abs(DFT)、Log、IDFT三步为一个卷积同态信号处理系统,经过这三步之后,倒谱域上
    s\'(n)=e\'(n)+h\'(n)
    激励源和信道已经变成了相加的关系,这时候通过一个倒滤波器,即图中的Cepstral Liftering,就可以把激励源和信道分开了。

    3.提取流程 (摘):

    MFCC参数的提取包括以下几个步骤:

    1. 预滤波:CODEC(所谓Codec,就是编码-解码器“Coder-Decoder”的缩写。说得通俗一点,对于音频就是A/D和D/A转换。)前端带宽为300-3400Hz(语音能量主要集中在250~4500Hz)。的抗混叠滤波器。

    工程测量中采样频率不可能无限高也不需要无限高,因为一般只关心一定频率范围内的信号成份。为解决频率混叠,在对模拟信号进行离散化采集前,采用低通滤波器滤除高于1/2采样频率的频率成份。实际仪器设计中,这个低通滤波器的截止频率(fc) 为:

     

      截止频率(fc)= 采样频率(fs) / 2.56

    2. A/D变换:8kHz的采样频率,12bit的线性量化精度。

    3. 预加重:通过一个一阶有限激励响应高通滤波器,使信号的频谱变得平坦,不易受到有限字长效应的影响。

    许多实际的消息信号,例如语言、音乐等,它们的功率谱随频率的增加而减小,其大部分能量集中在低频范围内。这就造成消息信号高频端的信噪比可能降到不能容许的程度。但是由于消息信号中较高频率分量的能量小,很少有足以产生最大频偏的幅度,因此产生最大频偏的信号幅度多数是由信号的低频分量引起。平均来说,幅度较小的高频分量产生的频偏小得多。所以调频信号并没有充分占用给予它的带宽。因为调频系统的传输带宽是由需要传送的消息信号(调制信号)的最高有效频率和最大频偏决定的。然而,接收端输入的噪声频谱却占据了整个调频带宽。这就是说,在鉴频器输出端噪声功率谱在较高频率上已被加重了。 为了抵消这种不希望有的现象,在调频系统中人们普遍采用了一种叫做预加重和去加重措施,其中心思想是利用信号特性和噪声特性的差别来有效地对信号进行处理。即在噪声引入之前采用适当的网络(预加重网络),人为地加重(提升)发射机输入调制信号的高频分量。然后在接收机鉴频器的输出端,再进行相反的处理,即采用去加重网络把高频分量去加重,恢复原来的信号功率分布。在去加重过程中,同时也减小了噪声的高频分量,但是预加重对噪声并没有影响,因此有效地提高了输出信噪比。

    4. 分帧:根据语音的短时平稳特性,语音可以以帧为单位进行处理,实验中选取的语音帧长为32ms,帧叠为16ms。

    5. 加窗:采用哈明窗对一帧语音加窗,以减小吉布斯效应的影响。

    吉布斯现象Gibbs phenomenon(又叫吉布斯效应): 将具有不连续点的周期函数(如矩形脉冲)进行傅立叶级数展开后,选取有限项进行合成。当选取的项数越多,在所合成的波形中出现的峰起越靠近原信号的不连续点。当选取的项数很大时,该峰起值趋于一个常数,大约等于总跳变值的9%。这种现象称为吉布斯现象。

    6. 快速傅立叶变换(Fast Fourier Transformation, FFT):将时域信号变换成为信号的功率谱。

     

    7. 三角窗滤波:用一组Mel频标上线性分布的三角窗滤波器(共24个三角窗滤波器),对信号的功率谱滤波,每一个三角窗滤波器覆盖的范围都近似于人耳的一个临界带宽,以此来模拟人耳的掩蔽效应。

     

    8. 求对数:三角窗滤波器组的输出求取对数,可以得到近似于同态变换的结果。

    9. 离散余弦变换(Discrete Cosine Transformation, DCT):去除各维信号之间的相关性,将信号映射到低维空间。

    10. 谱加权:由于倒谱的低阶参数易受说话人特性、信道特性等的影响,而高阶参数的分辨能力比较低,所以需要进行谱加权,抑制其低阶和高阶参数。

    11. 倒谱均值减(Cepstrum Mean Subtraction, CMS):CMS可以有效地减小语音输入信道对特征参数的影响。 
    12. 差分参数:大量实验表明,在语音特征中加入表征语音动态特性的差分参数,能够提高系统的识别性能。在本系统中,我们也用到了MFCC参数的一阶差分参数和二阶差分参数。 
    13. 短时能量:语音的短时能量也是重要的特征参数,本系统中我们采用了语音的短时归一化对数能量及其一阶差分、二阶差分参数。

     总结如下:

    Mel频率倒谱系数(MFCC)参数的提取步骤

    (1) 预加重(pre-emphasis)

    将经采样后的数字语音信号s(n)通过一个高通滤波器(high pass filter): H(z)= 1 – a×z -1 , 0.9  a  1.0 (一般取0.95左右)。经过预加重后的信号为: (n)= s(n)– a×s(n-1)。 因为发声过程中声带和嘴唇的效应,使得高频共振峰的振幅低于低频共振峰的振幅,进行预加重的目的就是为了消除声带和嘴唇的效应,来补偿语音信号的高频部分。

    (2) 分帧(frame blocking)

    一般取10-20ms为一帧,为了避免窗边界对信号的遗漏,因此对帧做偏移时候,要有帧迭(帧与帧之间需要重叠一部分)。一般取帧长的一半作为帧移,也就是每次位移一帧的二分之一后再取下一帧,这样可以避免帧与帧之间的特性变化太大。

    (3) 计算短时能量(energy)

    短时能量代表着音量的高低,亦即声音振幅的大小,可以根据此能量的值来过滤掉语音信号中的一些细微噪声。当一帧的能量值低于我们定的门槛值(threshold)时,则将此帧作为静音段(silence)。

    (4) 加窗(hamming window)

    语音在长范围内是不停变动的,没有固定的特性无法做处理,所以将每一帧代入窗函数,窗外的值设定为0,其目的是消除各个帧两端可能会造成的信号不连续性。常用的窗函数有方窗、汉明窗和汉宁窗等,根据窗函数的频域特性,常采用汉明窗。公式是在加窗范围内,w(n)=0.54-0.46*cos(2*pi*n/(n-1))。

    (5) 快速傅立叶变换(FFT transform)

    由于语音信号在时域上的变化快速而不稳定,所以通常都将它转换到频域上来观察,此时它的频谱会随着时间作缓慢的变化。所以通常将加窗后的帧经过FFT (Fast Fourier Transform)求出每帧的频谱参数。

    (6) 三角形带通滤波器(triangular band-pass filter)

    将每帧的频谱参数通过一组N个三角形带通滤波器(N一般为20~30个)所组成的梅尔刻度滤波器,将每个频带的输出取对数,求出每一个输出的对数能量(log energy),k = 1,2,… N。 再将此N个参数进行余弦变换(cosine transform)求出L阶的Mel-scale cepstrum参数。

    转载于:https://www.cnblogs.com/gogly/p/3440441.html

    展开全文
  • 分数梅尔倒谱系数在欺骗性语音检测中的应用
  • 梅尔倒谱系数

    2019-01-24 16:05:00
    MFCC梅尔倒谱系数是说话人识别、语音识别中最为常用的特征。我曾经对这个特征困惑了很久,包括为什么步骤中要取对数,为什么要最后一步要做DCT等等,以下将把我的理解记录下来,我找到的参考文献中最有价值的要数【1...

    MFCC梅尔倒谱系数是说话人识别、语音识别中最为常用的特征。我曾经对这个特征困惑了很久,包括为什么步骤中要取对数,为什么要最后一步要做DCT等等,以下将把我的理解记录下来,我找到的参考文献中最有价值的要数【1】了。是CUM一个教授做的PPT。

     

    整个流程如下:

    时域的波形图如下

     

    图1. 时域波形图

     

    第一步

    获得语谱图,语谱图是一个非常有力的工具,因为人耳就是进行的频率分析。

     

    图2. 语谱图

     

    第二步

    经过梅尔滤波器组。为什么要经过梅尔滤波器组?答:上面的图需要降维。根据生理学的发现,上面的语谱图实际上可以用经过一系列的梅尔滤波器组来进行降维。

     

    图3. 梅尔滤波器组

    滤波后的图像如下,假如一共有24个滤波器组,那么在下图在纵向上就降成了24维。

     

    图4. 经过梅尔滤波器组后的频谱图

     

    第三步

    取对数。为什么要取对数?解答如下。

    人类的发声系统发出的信号是由基音信息与声道信息卷积而成。记作"s卷积v"

    经过语谱图FFT变换后,卷积变成了乘法。即"FFT(s)*FFT(v)"。

    取对数后,乘法变成了加法。即"Log(FFT(s))+Log(FFT(v))"

    把卷积信号转换成加性信号,这就是取FFT和对数的原因。

     

    图5. 取对数后

     

    第四步

    DCT(离散余弦变换)

    在上一步中,我们成功地把基音信息与声道信息变成了加性的。那么如何分离呢?它们有如下性质:

    频谱图中(注意是一帧FFT变换内)

    (1)基音信息在频域是快速变化的。

    (2)声道信息在频域是缓慢变化的。

    因此再做一次DCT可以将其分离。我们称之为"倒谱域"。因此倒谱域的低频部分刻画了声道信息,高频部分刻画了基音信息。为什么是DCT而不是FFT?因为DCT变换之后的值仍为实数,因此更方便。

     

    图6. DCT变换后

     

    第五步

    对DCT变换后的谱图进行降维。

    (1)去掉第0维,因为第0维只是图5的均值,并不包含任何信息。

    (2)去掉13-24维,因为DCT本身就是用来去相关的,而图5没有太高频的成分,因此可以去掉。

     

    图7. 降维后的MFCC谱图

     

    图7就是最终的MFCC特征了!

     

    小结

    1. MFCC特征适用于说话人分类、语音识别,并且已经有了较好的识别结果。

    2. 虽然MFCC是个不错的特征,但是同时也丢掉了很多细节(图2至图4的过程),因此并不是非常完美。

     

    参考资料

    【1】http://download.csdn.net/detail/richard2357/6664585 (CMU的PPT,写的非常详细,我看过这个之后才真正理解)
    ---------------------
    作者:richard2357
    来源:CSDN
    原文:https://blog.csdn.net/richard2357/article/details/17147249
    版权声明:本文为博主原创文章,转载请附上博文链接!

    转载于:https://www.cnblogs.com/fpzs/p/10315044.html

    展开全文
  • MFCC(Mel 倒谱系数)

    千次阅读 2018-06-20 11:43:07
    Mel倒谱系数Mel倒谱系数:MFCC Mel频率倒谱系数(Mel Frequency Cepstrum Coefficient)的缩写是MFCC,Mel频率是基于人耳听觉特性提出来的,它与Hz频率成非线性对应关系。Mel频率倒谱系数(MFCC)则是利用它们之间的...
  • 基于深信度网络和小波包倒谱系数的语音情感识别
  • 基于梅尔频率倒谱系数的非负矩阵分解的音频哈希函数
  • Desicion-Maker2:转换为倒谱系数声音的定义
  • MFCC梅尔倒谱系数

    2018-06-08 14:18:58
    MFCC梅尔倒谱系数阅读数:7386MFCC梅尔倒谱系数是说话人识别、语音识别中最为常用的特征。我曾经对这个特征困惑了很久,包括为什么步骤中要取对数,为什么要最后一步要做DCT等等,以下将把我的理解记录下来,我找到...
  • MFCC倒谱系数

    千次阅读 2016-09-25 20:09:44
    MFCC是Mel频率倒谱系数(melfrequency cepstrum,MFCC)的缩写,Mel频率是基于人耳听觉特性提出来的,它与Hz频率成非线性对应关系。MFCC则是利用它们之间的这种关系计算得到频率特征,MFCC已经广泛应用在语音识别领域...
  • 本代码实现读入语音信号,提取该信号的梅尔倒谱系数,为后面的声音模板匹配打基础
  • 从periodogram estimate of the power spectrum计算得到的倒谱系数,可以用于基音追踪(pitch tracking),然而,从AR power spectral estimate计算得到的倒谱系数可以用于语音识别(现在已经被MFCCs所替代)。...
  • 针对LPCC只反应语音静态特征且不能突出其低频局部特征问题,提出一种以HHT倒谱系数为特征的说话人识别算法,HHT的经验模态分解使语音的低频局部特征得到更好的描述,Hilbert变换能够刻画语音动态特性,改进了LPCC的...
  • Mel频率倒谱系数

    千次阅读 2015-09-22 09:16:02
    MFCC:Mel频率倒谱系数(Mel Frequency Cepstrum Coefficient,MFCC)的缩写。Mel(美尔)是主观音高的单位,而Hz(赫兹)则是客观音高的单位。Mel频率是基于人耳听觉特性提出来的,它与Hz频率成非线性对应关系。Mel...
  • 参考文章:梅尔频率倒谱系数(MFCC) 学习笔记、声学特征(二) MFCC特征原理。   MFCC概述 在语音识别(SpeechRecognition)和话者识别(SpeakerRecognition)方面,最常用到的语音特征就是梅尔倒谱系数(Mel-...
  • 正弦信号的matlab代码梅尔频率倒谱系数 梅尔频率倒谱系数 该代码按照与Matlab中相同的步骤(功能:mfcc)来计算梅尔频率倒谱系数。 该代码使用默认的40频段滤波器组,其范围大约为133 Hz至6864 Hz,如Matlab中所述。...
  • Mel频率倒谱系数(Mel Frequency Cepstrum Coefficient)的缩写是MFCC,Mel频率是基于人耳听觉特性提出来的,它与Hz频率成非线性对应关系。Mel频率倒谱系数(MFCC)则是利用它们之间的这种关系,计算得到的Hz频谱特征...
  • 针对现有的煤矸界面识别技术采用的γ射线法不适用于顶板不含放射性元素或者放射性元素含量较低的工作面,而雷达探测法探测范围小、信号衰减严重的问题,提出了一种基于Mel频率倒谱系数和遗传算法的煤矸界面识别方法。...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 1,804
精华内容 721
关键字:

倒谱系数