python np fft_语音MFCC提取：librosa && python_speech_feature(2019.12)

最近在阅读语音方向的论文，其中有个被提及很多的语音信号特征MFCC(Mel-Frequency Cepstral Coefficients)，找到了基于python的语音库librosa(version=0.7.1)和python_speech_features(version=0.6)，下文对这两个库计算MFCC的流程细节稍作梳理。

LibROSA - librosa 0.7.1 documentation​librosa.github.io python_speech_features​pypi.org

一、librosa

1.源语音信号，shape = wav.length

wav, sample_rate = librosa.load(path, sr=22050, mono=True, offset=0.0, duration=None,
         dtype=np.float32, res_type='kaiser_best')  #加载语音文件得到原始数据

2.填充及分帧(无预加重处理)，分帧后所有帧的shape = n_ftt * n_frames

y = numpy.pad(array, pad_width, mode, **kwargs)     # 原函数及参数
y = numpy.pad(wav, int(n_fft // 2), mode='reflect') # 默认，n_fft 为傅里叶变换维度
y = numpy.pad(wav, (0,pad_need), mode='constant')   # 自定义

librosa调用numpy对源语音数据进行填充，默认模式是'reflect'进行镜像填充，举个例子：

对序列[1,2,3,4,5]进行左填充2个、右填充3个，左边以1作对称轴填充[3,2]，右边以5作对称轴填充[4,3,2]，最后结果为[3,2,1,2,3,4,5,4,3,2]。

自定义模式是按照python_speech_features的方式设置参数的，pad_need是经过计算得到的需要填充的数据数量，在下文说明。

y_frames = util.frame(y, frame_length=n_fft, hop_length=hop_length) # hop_length为帧移，librosa中默认取窗长的四分之一

一般来说 帧长 = 窗长 = 傅里叶变换维度，否则要进行填充或者截断处理。通过阅读源码发现，librosa调用的分帧方式和python_speech_features一样，这和librosa.feature.melspectrogram()参数里的center描述有所冲突，我把两者对同一语音文件的分帧结果输出到excel文件进行了对比，结果是一样的，当然分帧之前的填充方式也修改成了一致。

3.对所有帧进行加窗，shape = n_ftt * n_frames。librosa中window.shape = n_ftt * 1

fft_window = librosa.filters.get_window(window, Nx, fftbins=True) # 原函数及参数
fft_window = get_window('hann', win_length, fftbins=True)         # 窗长一般等于傅里叶变换维度，短则填充长则截断

librosa加的窗函数调用的scipy，如scipy.signal.windows.hann。python_speech_features加的窗函数调用的numpy，如numpy.hanning，汉宁窗公式为：

给原信号加窗的实现方式为相乘：

frames *= 0.5 - 0.5 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1)) # 原信号乘以汉宁窗函数

4.STFT处理得到spectrum(频谱，实际是多帧的)，shape = (n_ftt // 2 +1) * n_frames

fft = librosa.core.fft.get_fftlib()
stft_matrix = fft.rfft(fft_window * frames)

librosa和python_speech_features都是调用的numpy进行傅里叶变换，numpy.fft.rfft(frames,NFFT)。

5.取绝对值得到magnitude spectrum/spectrogram(声谱，包含时间维度，即多帧),shape = (n_ftt // 2 +1) * n_frames

magnitude_spectrum = numpy.abs(stft_matrix)    # 承接上一步的STFT
magnitude_spectrum = numpy.abs(librosa.core.spectrum.stft(wav, n_fft=n_fft, hop_length=hop_length,
                        win_length=win_length, center=center,
                        window=window, pad_mode=pad_mode))**power  # librosa封装的计算magnitude spectrum函数，power参数默认为1.0

6.取平方得到power spectrum/spectrogram(声谱，包含时间维度，即多帧),shape = (n_ftt // 2 +1) * n_frames

power_spectrum = numpy.square(magnitude_spectrum)  # 承接上一步的magnitude_spectrum
power_spectrum = numpy.abs(librosa.core.spectrum.stft(wav, n_fft=n_fft, hop_length=hop_length,
                        win_length=win_length, center=center,
                        window=window, pad_mode=pad_mode))**power  # librosa封装的计算power spectrum函数，power参数设置为2.0

7.构造梅尔滤波器组，shape = n_mels * (n_ftt // 2 +1)

mel_basis = librosa.filters.mel(sr, n_fft, n_mels=128, fmin=0.0, fmax=None, htk=False,
        norm=1, dtype=np.float32)     # 原函数及参数

librosa和python_speech_features构造mel_filtersbank的方式不同，librosa构造细节尚未掌握。

8.矩阵乘法得到mel_spectrogram，shape = n_mels * n_frames

mel_spectrogram = numpy.dot(mel_basis, power_spectrum)

[ n_mels ，(n_ftt // 2 +1) ] * [ (n_ftt // 2 +1) ，n_frames ] = [ n_mels，n_frames]

9.对mel_spectrogram进行log变换，shape = n_mels * n_frames

log_mel_spectrogram = librosa.core.spectrum.power_to_db(mel_spectrogram, ref=1.0, amin=1e-10, top_db=80.0)

S_db ~= 10 * log10(S) - 10 * log10(ref)

librosa采用以上方式计算log映射，python_speech_features定义mfcc( )函数时直接取自然常数e为底的对数。

10.IFFT变换，实际采用DCT得到MFCC，shape = n_mels * n_frames

mfcc = scipy.fftpack.dct(log_mel_spectrogram, type=2, n=None, axis=0, norm=None, overwrite_x=False) 
# n表示计算维度，需与log_mel_spectrogram.shape[axis]相同,否则作填充或者截断处理。axis=0表示沿着自上而下的方向，分别选取每一行所在同一列的元素进行运算。

与python_speech_features相同，librosa也是调用scipy对log_mel_spectrogram进行离散余弦变换：scipy.fftpack.dct()。

11.取MFCC矩阵的低维(低频)部分，shape = n_mfcc * n_frames

mfcc = mfcc[ :n_mfcc]  # 取低频维度上的部分值输出，语音能量大多集中在低频域，数值一般取13。

二、python_speech_features

1.源语音信号，shape = wav.length

sample_rate,signal = scipy.io.wavfile.read(filename, mmap=False) # scipy加载语音文件
signal, sample_rate = librosa.load(path, sr=22050, mono=True, offset=0.0, duration=None,
         dtype=np.float32, res_type='kaiser_best')            # librosa加载语音文件

2.预加重，shape = wav.length

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

语音信号的预加重，目的是为了对语音的高频部分进行加重，去除口唇辐射的影响，增加语音的高频部分的分辨率和信噪比，计算方式为：y(n)）=x(n)-a*x(n-1)，系数a一般取0.97。

3.分帧和加窗，shape = n_frames * n_ftt

frames = python_speech_features.sigproc.framesig(signal,frame_len,frame_step,winfunc=lambda x:numpy.ones((x,))) # 分帧及加窗函数

python_speech_features中默认帧长取25ms、帧移10ms，对应实际的采样点为：

帧长点数 = 帧长时间 * 采样率 ，帧移点数 = 帧移时间 * 采样率

numframes = 1 + int(math.ceil((1.0*slen - frame_len)/frame_step))  # 计算总分帧数

总分帧数量 = 1 + 向上取整( ( 原信号长度 - 帧长点数 ) / 帧移点数 )

padlen = int((numframes-1)*frame_step + frame_len)     # 计算待填充点的数量
zeros = numpy.zeros((padlen - slen,))                  # 在原数据末尾进行零填充

待填充点数 = ( 总分帧数量 - 1 ) * 帧移点数 + 帧长点数

这种填充方式和上面librosa的自定义的那种填充方式( y = numpy.pad(wav, (0,pad_need), mode='constant') )效果是一样的，可以看出librosa库的功能还是更加多样化的。

加窗的话，python_speech_features默认不加窗，但提供了调用numpy中窗函数的参数接口，经测试numpy.hanning窗函数和scipy.signal.windows.hann窗函数的数值是一致的，只不过前者为矩阵形式(元素相同的多个向量构成)，后者为向量形式。

mfcc = python_speech_features.base.mfcc(signal,samplerate=16000,winlen=0.025,winstep=0.01,numcep=13,
         nfilt=26,nfft=512,lowfreq=0,highfreq=None,preemph=0.97,ceplifter=22,appendEnergy=True,
         winfunc=lambda x:numpy.ones((x,))) 
# winfunc参数即为窗函数接口，例如 winfunc=numpy.hanning

4.STFT处理得到spectrum(频谱，实际是多帧的)，shape = n_frames * (n_ftt // 2 +1)

complex_spec = numpy.fft.rfft(frames,n_fft)    # 离散傅里叶变换得到频谱图

和librosa一致，python_speech_features也是调用numpy下的函数做离散傅里叶变换。

5.取绝对值得到magnitude spectrum/spectrogram(声谱，包含时间维度，即多帧),shape = n_frames * (n_ftt // 2 +1)

mag_spec = numpy.absolute(complex_spec)       # 承接上一步取绝对值得到幅度谱

6.取平方得到power spectrum/spectrogram(声谱，包含时间维度，即多帧),shape = n_frames * (n_ftt // 2 +1)

power_spec = 1.0 / n_fft * numpy.square(mag_spec) # 承接上一步取平方得到功率谱

python_speech_features计算功率谱的方式和libraosa不一致，这里额外除以了傅里叶变换的维度n_fft。

7.构造梅尔滤波器组，shape = n_mels * (n_ftt // 2 +1)

fb = python_speech_features.base.get_filterbanks(nfilt=20,nfft=512,samplerate=16000,lowfreq=0,highfreq=None)

关于梅尔滤波器的构造细节，参考了这篇文章的部分介绍。具体细节如下：

赫兹转梅尔：mel = 2595 * log10( 1+ hz / 700. )

梅尔转赫兹：hz = 700 * ( 10** ( mel / 2595.0 ) - 1 )

(1)假设构造n_mels=10个梅尔滤波器，建立坐标系纵坐标表示mel刻度、横坐标表示Hz刻度

(2)横坐标(频率)设置最低值300hz、最高值8000hz(采样率的二分之一)，对应的纵坐标的最低和最高值分别为401.97mel和2840.02mel。

(3)把纵坐标从低到高划分为12(10+2)个离散的点melpoints = [ 401.97， 623.61，845.25，1066.89，1288.54，1510.18， 1731.82，1953.46， 2175.10，2396.74，2618.38，2840.02] ，这些点映射回横坐标上分别为[ 300，517.34，781.91，1103.98，1496.06，
1973.34，2554.36，3261.65，4122.66，5170.80，6446.75， 8000. ]。

(4)上面横坐标这12(0-11)个点正好对应10个三角滤波器的三个点：例如0、1、2分别对应第一个滤波器的左底点、顶点、右底点，1、2、3对应第二个滤波器，以此类推第十个滤波器对应9、10、11。这十个三角滤波器的函数解析式如下：

(5)对信号的滤波处理本质上是乘法运算，由于待滤波的spectrogram实际是离散的数据，所以要对连续的滤波器函数进行离散的采样，从源代码来看python_speech_features和librosa的采样方式是不同的，经测试发现两者所构造的梅尔滤波器数据确实不一样。

(6)拿第一个滤波器举例来说，其横坐标的采样点选择是根据经傅里叶变换后的spectrogram维度确定的：

bin = numpy.floor( ( n_fft + 1 ) * mel2hz( melpoints ) / samplerate ) # 例子中n_fft取值512,samplerate取值16000
# 构造梅尔采样容器，实际上就是把原横坐标"samplerate/2"上的12个离散点映射到新横坐标"(n_ftt+1)/2"上面，这样做保证了后面的乘法运算维度一致。
# bin = [  9.  16.  25.  35.  47.  63.  81. 104. 132. 165. 206. 256.]，例子中n_fft取值512

经过以上映射后，第一个滤波器的函数表达式就确定了：在[9,16)区间为y=(x-9)/(16-9)，在[16,25]区间为y=(25-x)/(25-16)，其余点y取值为0。在这个滤波器的257(0-256)个离散点中，沿着三角形左边y值逐渐递增，在顶点处达到最大值"1"，然后顺着右边y值逐渐递减到0，其余部分点对应的纵坐标取值都为0。

(7)例子中(n_mel=10,n_fft=512,samplerate=16000)构造的梅尔滤波器矩阵具体数值如下(作了转置处理)：

8.矩阵乘法得到mel_spectrogram，shape = n_frames * n_mels

mel_spectrogram = numpy.dot(power_spec,fb.T)

[ n_frames ，(n_ftt // 2 +1) ] * [ (n_ftt // 2 +1) ，n_mels ] = [ n_frames，n_mels]

mel_spectrogram矩阵中的每个元素其实就是每一帧spectrogram向量和每个滤波器向量的内积。mel_spectrogram(0,1)即为第1帧的spectrogram同第2个三角滤波器内积运算的结果。

9.对mel_spectrogram进行log变换，shape = n_frames * n_mels

log_mel_spectrogram = numpy.log(mel_spectrogram) 
log_mel_spectrogram = python_speech_features.sigproc.logpowspec(frames,n_ftt,norm=1)

python_speech_features.base.mfcc( )函数默认取以自然常数e为底的对数，作者还提供了另一种取对数域的接口python_speech_features.sigproc.logpowspec( )，计算方式与librosa略有差异：log_S = 10 * log10(S)

10.IFFT变换，实际采用DCT得到MFCC，shape = n_frames * n_mels

mfcc = scipy.fftpack.dct(log_mel_spectrogram, type=2, n=None, axis=0, norm=None, overwrite_x=False) 
# n表示计算维度，需与log_mel_spectrogram.shape[axis]相同,否则作填充或者截断处理。axis=0表示自上而下方向分别选取每一行所在同一列的元素进行计算。

11.取MFCC矩阵的低维(低频)部分，shape = n_frames * n_mfcc

mfcc = mfcc[:,:n_mfcc]  # 取低频维度上的部分值输出，语音能量大多集中在低频域，数值一般取13。

python_speech_features.base.mfcc(appendEnergy=True )函数中，appendEnergy参数控制是否把MFCC的第一个倒谱系数替换为每一帧总能量的对数，每一帧总能量的计算方式为：

energy = numpy.sum(power_spec,axis=1) # axis=1表示沿着从左到右的方向,分别选取每一列所在同一行的元素进行运算。

12.倒谱提升，shape = n_frames * n_mfcc

mfcc = python_speech_features.base.lifter(cepstra=mfcc, L=22)

倒谱提升系数默认设置为22，具体实现方式为：

nframes,ncoeff = numpy.shape(cepstra)              # ncoeff=n_mfcc
n = numpy.arange(ncoeff)                           # n=0,1,2...n_mfcc-1
lift = 1 + (L / 2.) * numpy.sin(numpy.pi * n / L)  # lift.shape = n_mfcc * 1
return lift * cepstra                              # shape = n_frames * n_mfcc

13.微分：mfcc的动态特征提取，shape = n_frames * n_mfcc

mfcc_delta_1 = python_speech_features.base.delta(feat=mfcc, N=1)  # 计算mfcc的一阶微分
mfcc_delta_2 = python_speech_features.base.delta(feat=mfcc, N=2)  # 计算mfcc的二阶微分

有时会把MFCC的基础特征同其一阶、二阶微分数据结合起来使用，以做到特征层面的动静结合。