Python语音信号处理

    语言信息是多种信息的混合载体 ,其中包括内容信息、说话人信息和情感信息。 本文介绍了一些语音的基本知识,和使用Python进行处理。


时域特征

    使用wave模块读取wav音频文件,画图时域图像,代码如下。

import numpy as np
import matplotlib.pyplot as plt
import os
import wave

path='D://NLP//dataset//语音情感//test.wav'
f=wave.open(path,'rb')
params=f.getparams()
#通道数、采样字节数、采样率、采样帧数
nchannels,sampwidth,framerate,nframes=params[:4]
voiceStrData=f.readframes(nframes)
waveData = np.fromstring(voiceStrData,dtype=np.short)#将原始字符数据转换为整数
#音频数据归一化
waveData = waveData * 1.0/max(abs(waveData))
#将音频信号规整乘每行一路通道信号的格式,即该矩阵一行为一个通道的采样点,共nchannels行
waveData = np.reshape(waveData,[nframes,nchannels]).T # .T 表示转置
f.close()

time=np.arange(0,nframes)*(1.0/framerate)
plt.plot(time,waveData[0,:],c='b')
plt.xlabel('time')
plt.ylabel('am')
plt.show()

代码执行结果:

index.png


频域特征

    numpy模块自带了快速傅里叶变换的函数,对上面的音频数据进行傅里叶变换,代码如下:

fftdata=np.fft.fft(waveData[0,:])
fftdata=abs(fftdata)
hz_axis=np.arange(0,len(fftdata))
plt.figure()
plt.plot(hz_axis,fftdata,c='b')
plt.xlabel('hz')
plt.ylabel('am')
plt.show()

程序运行结果:

fft.png


语谱图

    语音的时域分析和频域分析是语音分析的两种重要方法,但是都存在着局限性。时域分析对语音信号的频率特性没有直观的了解,频域特性中又没有语音信号随时间的变化关系。语谱图的原理如下。

20130623205844093.jpg

    这段语音被分为很多帧,每帧语音都对应于一个频谱(通过短时FFT计算),频谱表示频率与能量的关系。在实际使用中,频谱图有三种,即线性振幅谱、对数振幅谱、自功率谱(对数振幅谱中各谱线的振幅都作了对数计算,所以其纵坐标的单位是dB(分贝)。这个变换的目的是使那些振幅较低的成分相对高振幅成分得以拉高,以便观察掩盖在低幅噪声中的周期信号)。

20130623205928484.jpg

    我们先将其中一帧语音的频谱通过坐标表示出来,如上图左。现在我们将左边的频谱旋转90度。得到中间的图。然后把这些幅度映射到一个灰度级表示,0表示黑,255表示白色。幅度值越大,相应的区域越黑。这样就得到了最右边的图。那为什么要这样呢?为的是增加时间这个维度,这样就可以显示一段语音而不是一帧语音的频谱,而且可以直观的看到静态和动态的信息。

    这样我们会得到一个随着时间变化的频谱图,这个就是描述语音信号的spectrogram语谱图或声谱图。

20130623210019187.jpg

    而语谱图综合了时域和频域的优点,明显的显示出了语音频谱随时间的变化情况、语谱图的横轴为时间,纵轴为频率,任意给定频率成分在给定时刻的强弱用颜色深浅来表示。颜色深的,频谱值大,颜色浅的,频谱值小。语谱图上不同的黑白程度形成不同的纹路,称之为声纹,不同讲话者的声纹是不一样的,可用作声纹识别。

    语谱图分为功率谱,幅度谱,相位谱。使用matplotlib可以直接获得语谱图,代码如下:

#帧长20~30ms
framelength = 0.025 
#每帧点数 N = t*fs,通常情况下值为256或512,要与NFFT相等
#而NFFT最好取2的整数次方,即framesize最好取的整数次方
framesize = framelength*framerate  
#找到与当前framesize最接近的2的正整数次方
nfftdict = {}
lists = [32,64,128,256,512,1024]
for i in lists:
    nfftdict[i] = abs(framesize - i)
sortlist = sorted(nfftdict.items(), key=lambda x: x[1])#按与当前framesize差值升序排列
framesize = int(sortlist[0][0])#取最接近当前framesize的那个2的正整数次方值为新的framesize
 
NFFT = framesize #NFFT必须与时域的点数framsize相等,即不补零的FFT
overlapSize = 1.0/3 * framesize #重叠部分采样点数overlapSize约为每帧点数的1/3~1/2
overlapSize = int(round(overlapSize))#取整
spectrum,freqs,ts,fig = plt.specgram(waveData[0],NFFT = NFFT,Fs =framerate,window=np.hanning(M = framesize),noverlap=overlapSize,mode='default',scale_by_freq=True,sides='default',scale='dB',xextent=None)#绘制频谱图         
plt.ylabel('Frequency')
plt.xlabel('Time(s)')
plt.title('Spectrogram')

    plt.specgram()函数原型如下:

matplotlib.pyplot.specgram(x,NFFT=None,Fs=None,Fc=None,detrend=None,window=None,noverlap=None,cmap=None,xextent=None,
pad_to=None,sides=None,scale_by_freq=None,mode=None,scale=None,vmin=None,vmax=None,*,data=None,**kwargs)
参数如下:
x : 1-D阵列或序列包含数据的数组或序列。
Fs : 采样频率(每个时间单位的采样)。它用于计算每个时间单位的周期的傅里叶频率,频率。默认值为2。
window : callable或ndarray长度为NFFT的函数或向量。要创建窗口向量看到 window_hanning,window_none,numpy.blackman,numpy.hamming,numpy.bartlett,scipy.signal,scipy.signal.get_window,等默认为window_hanning。如果函数作为参数传递,则它必须将数据段作为参数并返回该段的窗口版本。
sides: {'default','onesided','twosided'}指定要返回的频谱的哪一侧。default为默认行为,对于实际数据和复杂数据都返回单侧行为。'onesided'迫使单侧频谱返回,而'twosided'强制双侧。
pad_to : 执行FFT时数据段填充的点数。这可以与NFFT不同,后者指定使用的数据点的数量。虽然没有增加频谱的实际分辨率(可分辨峰值之间的最小距离),但这可以在图中给出更多的点,从而允许更多细节。这对应于fft()调用中的n参数。默认值是无,这将pad_to 等于NFFT。
NFFT : 每个块中用于FFT的数据点数。2的指数值是最有效的。默认值为256.这不应该用于获得零填充,否则结果的缩放将不正确。请改用pad_to。
detrend : {'none','mean','linear'}或callable,默认'none',该函数在fft之前应用于每个段,旨在消除均值或线性趋势。与MATLAB不同,其中 detrend参数是向量,在Matplotlib中它是一个函数。该mlab模块定义detrend_none,detrend_mean和detrend_linear,但您也可以使用自定义函数。您还可以使用字符串选择其中一个函数:'none'调用detrend_none。'mean'的调用detrend_mean。'linear'调用detrend_linear。
scale_by_freq : bool,可选指定是否应通过缩放频率缩放得到的密度值,缩放频率以Hz ^ -1为单位给出密度。这允许在返回的频率值上进行积分。对于MATLAB兼容性,默认值为True。
mode : {'default','psd','magnitude','angle','phase'}使用什么样的频谱。默认为'psd',它取功率谱密度。“magnitude”返回幅度谱。'angle'返回相位谱而不展开。'phase'通过展开返回相位谱。
noverlap : 块之间的重叠点数。默认值为128。
scale : {'default','linear','dB'}规范中值的缩放。'linear'没有缩放。'dB'以dB刻度返回值。当模式为“psd”时,这是dB功率(10 * log10)。否则这是dB幅度(20 * log10)。如果mode为'psd'或'magnitude',则'default'为'dB',否则为'linear'。如果模式为“angle”或“phase”,则必须为“linear”。
Fc : x的中心频率(默认为0),它偏移图的x范围,以反映采集信号时使用的频率范围,然后进行滤波并下采样到基带。
CMAP:一个matplotlib.colors.Colormap例子; 如果为None,则使用由rc确定的默认值
xextent : None或(xmin,xmax,沿x轴的图像范围。默认设置将xmin设置为第一个bin(频谱列)的左边界,将xmax设置为最后一个bin的右边界。请注意,对于noverlap> 0,bin的宽度小于segment的宽度。
** kwargs其他关键字参数传递给imshow,这将生成specgram图像。

程序运行结果:

sp.png

    语谱图分为窄带和宽带语谱图。“窄带”,顾名思义,带宽小,则时宽大,则短时窗长,窄带语谱图就是长窗条件下画出的语谱图。“宽带”,正好相反。至于“横竖条纹”,窄带语谱图的带宽窄,那么在频率上就“分得开”,即能将语音各次谐波“看得很清楚”,即表现为“横线”。“横”就体现出了频率分辨率高。分辨率可以直观的看做“分开能力”。“频率分辨率”高就是在频率上将各次谐波分开的能力高,表现为能分辨出各次谐波的能力高,频率分辨率越高,越容易分辨各次谐波。类似的,宽带语谱图的时宽窄,那么在时间上就“分得开”,即能将语音在时间上重复的部分“看得很清楚”,即表现为“竖线”。“竖”就体现出了时间分辨率高。时间分辨率越高,谱图上的竖线看得越清楚。图1和图2分别示出了一条语音句子的窄带语谱图和宽带语谱图。短时窗长度分别是20ms和2ms。

    下图1是一条语音及其窄带语谱图。上是语音时间波形,下是上的窄带语谱图:

20171223221201396.jpg

    下图2是同一条语音及其宽带语谱图,上是语音时间波形,下是上的宽带语谱图:

A20171223221347005.jpg


梅尔语谱图

    人耳对声音频谱的响应是非线性的,经验表明:如果我们能够设计一种前端处理算法,以类似于人耳的方式对音频进行处理,可以提高语音识别的性能。FilterBank分析就是这样的一种算法。FBank特征提取要在预处理之后进行,这时语音已经分帧,我们需要逐帧提取FBank特征,FBank特征也叫梅尔语谱图,与前面的线性语谱图相区别。

    梅尔频率是基于人耳听觉特性提出来的, 它与Hz频率成非线性对应关系。梅尔语谱图的提取过程和线性语谱图差不多,流程如下:

0.预加重。

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

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

1、对波形图分帧。

    我们需要将不定长的音频切分成固定长度的小段,这一步称为分帧。一般取10-30ms为一帧,为了避免窗边界对信号的遗漏,因此对帧做偏移时候,要有帧迭(帧与帧之间需要重叠一部分)。 一般取帧长的一半作为帧移,也就是每次位移一帧的二分之一后再取下一帧,这样可以避免帧与帧之间的特性变化太大。通常的选择是25ms每帧,帧迭为10ms。接下来的操作是对单帧进行的。

    要分帧是因为语音信号是快速变化的,而傅里叶变换适用于分析平稳的信号。在语音识别中,一般把帧长取为10~30ms,这样一帧内既有足够多的周期,又不会变化太剧烈。每帧信号通常要与一个平滑的窗函数相乘,让帧两端平滑地衰减到零,这样可以降低傅里叶变换后旁瓣的强度,取得更高质量的频谱。帧和帧之间的时间差常常取为10ms,这样帧与帧之间会有重叠,否则,由于帧与帧连接处的信号会因为加窗而被弱化,这部分的信息就丢失了。傅里叶变换是逐帧进行的,为的是取得每一帧的频谱。一般只保留幅度谱,丢弃相位谱。

    例如我们取40ms位一帧的宽度,对于一个44.1kHz采样的信号,一帧就包含0.040*44100=1764个采样点,帧移通常去帧宽的二分之一,也就是20ms,这样就允许没两帧之间有一半的overlap。这样一来,第一帧就是从第一个采样点到第1764个采样点,第二帧就是从第882个采样点到第2646个采样点...直到最后一个采样点,如果音频长度不能被帧数整除,在最后补0 。对于一个30s的音频文件,可以得到44100*30 / 882 = 1500帧。

2、对每一帧进行加窗。加窗的目的是平滑信号,使用汉明窗加以平滑的话,相比于矩形窗函数,会减弱FFT以后旁瓣大小以及频谱泄露。例如使用汉明窗(hamming window)对信号进行加窗处理:

1.png

2.png

    从上面的例子可以看出来,如果不进行加窗,那么某一帧的结束值和下一帧的开始值会有一个gap(因为有overlap),在频谱图上来看,峰值会变得比较“宽”,而加过hamming window的帧在频谱图上的峰值就会变得更加sharp也更容易辨认。设加窗函数为h(n)。

3、对每一帧进行离散傅里叶变化(DFT)。

    我们分帧之后得到的仍然是时域信号,为了提取fbank特征,首先需要将时域信号转换为频域信号。傅里叶变换可以将信号从时域转到频域。傅里叶变换可以分为连续傅里叶变换和离散傅里叶变换,因为我们用的是数字音频(而非模拟音频),所以我们用到的是离散傅里叶变换。

    信号频率 : 组合产生复杂信号的简单信号的频率,通常简单信号平率范围很广

    采样频率 : 模拟到数字的转换过程中,需要对模拟信号进行采样,每秒内的采样点数量就是采样频率

    Nyquist定理(奈奎斯特): 如果想要从数字信号无损转到模拟信号,我们需要以最高信号频率的2倍的采样频率进行采样。通常人的声音的频率大概在3kHz~4kHz ,因此语音识别通常使用8k或者16k的wav提取特征。16kHz采样率的音频,傅里叶变换之后的频率范围为0-8KHz。

    离散傅里叶变换公式如下:

3.png

其中s(n)为波形信号,S(n)为幅度谱。因为有欧拉公式如下:

4.png

5.png

得DFT公式的另一种形式:

6.png

所以其实DFT变换就是两个“相关(correlation)”操作,一个是与频率为k的cos序列相关,一个是与频率为k的sin序列相关,然后两者叠加就是与频率k的正弦波相关的结果,如果得到的值很大,就表明信号包含频率为k的能量很大。

4、计算功率谱。

    例如,我们从一个1764个点的FFT计算得到功率谱以后,只保留前1764/2+1=883个系数。

5、计算Mel-spaced filterbank。

频率和mel频率之间的转化公式为:

    mel滤波器组是一组非线性分布的滤波器组,它在低频部分分布密集,高频部分分布稀疏,这样的分布是为了更好得满足人耳听觉特性。

7.png

将这样一组三角滤波器(例如128个)作用到一帧上,就将一个883维的向量转化为128维的向量。

8.png


6、对上述128维的mel功率谱取log,得到128维的 log-mel filer bank energies。这样做的原因是由于人耳对声音的感知并不是线性的,用log这种非线性关系更好描述,另外,取完log以后才可以进行倒谱分析。

    经过上面这几步就得到了梅尔语谱图。


梅尔频率倒谱系数

    MFCCs(Mel-frequency cepstral coefficients)即梅尔频率倒谱系数。mfcc特征在梅尔语谱图的基础上增加了dct倒谱环节。FBank特征已经很贴近人耳的响应特性,但是仍有一些不足:FBank特征相邻的特征高度相关(相邻滤波器组有重叠),因此当我们用HMM对音素建模的时候,几乎总需要首先进行倒谱转换,通过这样得到MFCC特征。

    离散余弦变换DCT。DCT和DFT类似,但是只使用实数,不涉及复数运算。

其中

引入ai是为了使ci正交化。将ci表示为矩阵形式:

    这样得到的C(n)矩阵中,较大的值都集中再靠近左上角的低能量部分,其余部分会产生大量的0或者接近0的数,这样,我们可以进行进一步数据压缩。这也表明DCT有很好的能量聚集效应。相比于傅里叶变换(FFT),离散余弦变换的结果没有虚部,更好计算(DCT也可以理解为没有虚部的FFT)。利于对于ASR任务,通常取低13维的系数。

    这样我们就得到了13banks 的MFCC。


    通过python提取MFCCs有两种方式,通过librosa模块或者python_speech_features模块,代码如下。

    1.通过python_speech_features提取mfcc

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import wavfile
from python_speech_features import mfcc, logfbank

# 读取输入音频文件
sampling_freq, audio = wavfile.read(path)
# 提取MFCC和滤波器组特征
mfcc_features = mfcc(audio, sampling_freq)
filterbank_features = logfbank(audio, sampling_freq)
print('\nMFCC:\n窗口数 =', mfcc_features.shape[0])
print('每个特征的长度 =', mfcc_features.shape[1])
print('\nFilter bank:\n窗口数 =', filterbank_features.shape[0])
print('每个特征的长度 =', filterbank_features.shape[1])
# 画出特征图,将MFCC可视化。转置矩阵,使得时域是水平的
mfcc_features = mfcc_features.T
plt.matshow(mfcc_features)
plt.title('MFCC')
# 将滤波器组特征可视化。转置矩阵,使得时域是水平的
filterbank_features = filterbank_features.T
plt.matshow(filterbank_features)
plt.title('Filter bank')
plt.show()

运行结果:

mfcc1.png


    2.通过librosa提取mfcc

    需要说明的是,librosa.load()函数是会改变声音的采样频率的。如果 sr 缺省,librosa.load()会默认以22050的采样率读取音频文件,高于该采样率的音频文件会被下采样,低于该采样率的文件会被上采样。因此,如果希望以原始采样率读取音频文件,sr 应当设为 None。该函数返回的参数y是经过归一化的声音数据

import librosa 

y,sr = librosa.load(path,sr=None)
'''
librosa.feature.mfcc(y=None, sr=22050, S=None, n_mfcc=20, dct_type=2, norm='ortho', **kwargs)
y:声音信号的时域序列
sr:采样频率(默认22050)
S:对数能量梅尔谱(默认为空)
n_mfcc:梅尔倒谱系数的数量(默认取20)
dct_type:离散余弦变换(DCT)的类型(默认为类型2)
norm:如果DCT的类型为是2或者3,参数设置为"ortho",使用正交归一化DCT基。归一化并不支持DCT类型为1
kwargs:如果处理时间序列输入,参照melspectrogram
返回:
M:MFCC序列
'''
mfcc_data = librosa.feature.mfcc( y,sr,n_mfcc=13)

plt.matshow(mfcc_data)
plt.title('MFCC')

运行结果:

mfcc2.png

    从上面的代码可以看到,这两个库提取出的mfcc是不一样的。哪里不一样呢?

    librosa.feature.mfcc这个函数内部是这样的:

# -- Mel spectrogram and MFCCs -- #
def mfcc(y=None, sr=22050, S=None, n_mfcc=20, **kwargs):
    if S is None:
        S = logamplitude(melspectrogram(y=y, sr=sr, **kwargs))
    return np.dot(filters.dct(n_mfcc, S.shape[0]), S)

    可以看出,如果只给定原始的时域信号(即S参数为None),librosa会先通过melspectrogram()函数先提取时域信号y的梅尔语谱图,存放到S中,再通过filters.dct()函数做dct变换得到y的梅尔倒谱系数。melspectrogram()函数是这样的,librosa.feature.melspectrogram(y=Nonesr=22050S=Nonen_fft=2048hop_length=512win_length=Nonewindow='hann'center=Truepad_mode='reflect'power=2.0**kwargs).

参数说明:

9.png

    可以看到其默认的窗长是2048点,这和python_speech_features库是不同的。


差分

    由于语音信号是时域连续的,分帧提取的特征信息只反应了本帧语音的特性,为了使特征更能体现时域连续性,可以在特征维度增加前后帧信息的维度。MFCC很好的表达了语音的特征,但只是静态的特征。提取动态特征,一般都采用一阶二阶差分。差分公式如下:

    用python的librosa库实现一次差分:

a = [1 2 3 4 5 1 3 5 7] 
b = liborsa.feature.delta(a, width=3) 
b= [ 0.5 1. 1. 1. -1.5 -1. 2. 2. 1. ]

    例如对于a中的第一个5,b中对应的值是前后两个值的差值除以(宽度-1): (1-4)/2

    二阶差分则是对delta再做一次差分运算。


    整个梅尔语音特征的提取过程如下:

20130623210522390.jpg



参考文献

[1]fengzhonghen.浅谈MFCC.https://blog.csdn.net/fengzhonghen/article/details/51722555. 2016年06月20日

[2]开飞机的猪猪.MFCC(Mel Frequency Cepstral Coefficient)提取过程详解.https://www.jianshu.com/p/73ac4365b1d3.  2017.09.26

[3]Tingwei_chen.python生成语谱图.https://www.cnblogs.com/tingweichen/p/9863092.html2018-10-27

[4]wxysunshy.FBank与MFCC.https://blog.csdn.net/sun___shy/article/details/82994455.2018年10月10日

[5]图图噜.语谱图 基频 共振峰 .https://blog.csdn.net/lzrtutu/article/details/78882715. 2017年12月23日

[6] zkl99999.梅尔频率倒谱系数(MFCC) 学习笔记,很全面很详细的一篇入门.https://blog.csdn.net/zkl99999/article/details/80723755. 2018年06月18日


首页 所有文章 机器人 计算机视觉 自然语言处理 机器学习 编程随笔 关于