MFCC梅尔频率倒谱系数

于 2021-01-30 14:43:56 发布
阅读量2.1k 收藏 27
点赞数 3
文章标签： python 信号处理 算法 傅立叶分析
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/weixin_43545253/article/details/113309973
版权

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开始，直到到达语音文件的末尾，如果语音文件没有分成偶数个帧，则用零填充赝本集使帧个数成偶数个。
    帧长 = f s . 0.025 ( 采样点 ) 帧移 = f s . 0.01 ( 采样点 ) \text{帧长}=f_s.0.025\left( \text{采样点} \right) \\ \text{帧移}=f_s.0.01\left( \text{采样点} \right) 帧长=fs​.0.025(采样点)帧移=fs​.0.01(采样点)

  • 对于每个单个帧，为每个帧提取一组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 ) = ( 1 − a ) − a ⋅ cos ⁡ ( 2 ⋅ π ⋅ n N )    1 ⩾ n ⩾ N 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 W(n)=(1−a)−a⋅cos(2⋅π⋅Nn​)1⩾n⩾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个倒谱系数。

  • M F C C ( i , n ) = ∑ m = 1 M log ⁡ [ H ( i , m ) ] . cos ⁡ [ π . n . ( 2 m − 1 ) 2 M ] 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]} MFCC(i,n)=m=1∑M​log[H(i,m)].cos[2Mπ.n.(2m−1)​]

    其中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的备忘，读者可忽略或指正。

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

灰沿
关注
评论 2
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值