语音信号处理一

前言

准备入门语音处理，记录下自己学习的情况，如有错误之处，欢迎各位大佬指正。

音频信号

声音是一种机械波，即声波，是一种信号，一般称为音频信号。原始的音频信号通常由于人类发声器官或者语音采集设备所带来的静音片段、混叠、噪声、高次谐波失真等因素，一定程度上会对语音信号质量产生影响。所以在正式使用声学模型进行语音识别之前，我们必须对音频信号进行预处理和特征提取。

最初始的预处理工作就是静音切除，也叫语音激活检测(Voice Activity Detection, VAD) 或者语音边界检测。目的是从音频信号流里识别和消除长时间的静音片段，在截取出来的有效片段上进行后续处理会很大程度上降低静音片段带来的干扰。除此之外，还有许多其他的音频预处理技术，大家可以找来信号处理相关的资料进行阅读，笔者这里不多说。

然后就是特征提取工作。音频信号中通常包含着非常丰富的特征参数，不同的特征向量表征着不同的声学意义，从音频信号中选择有效的音频表征的过程就是语音特征提取。常用的语音特征包括线性预测倒谱系数（LPCC）和梅尔频率倒谱系数（MFCC），其中 LPCC 特征是根据声管模型建立的特征参数，是对声道响应的特征表征。而 MFCC 特征是基于人的听觉特征提取出来的特征参数，是对人耳听觉的特征表征。所以，在对音频信号进行特征提取时通常使用 MFCC 特征。

MFCC 主要由预加重、分帧、加窗、快速傅里叶变换(FFT)、梅尔滤波器组、离散余弦变换几部分组成，其中FFT与梅尔滤波器组是 MFCC 最重要的部分。一个完整的 MFCC 算法包括如下几个步骤：

• 快速傅里叶变换（FFT）
• 梅尔频率尺度转换
• 配置三角形滤波器组并计算每一个三角形滤波器对信号幅度谱滤波后的输出
• 对所有滤波器输出作对数运算，再进一步做离散余弦变换(DTC)，即可得到MFCC。

傅里叶变换

一个在长度为P的定义域里的可积(integrable)实数函数 [公式] 可以展开成无限个三角函数的和，也就是傅里叶级数：
$$f(x)=\frac{a_0}{2}+\sum _{n=1}^N\left(a_n\cos \left(\frac{2\pi nx}{P}\right)+b_n\sin \left(\frac{2\pi nx}{P}\right)\right)$$
也可以变换为如下形式：
$$f(x)=\frac{A_0}{2}+\sum _{N=1}^N{A}_n\cos \left(\frac{2\pi nx}{P}-{\varphi }_n\right)$$
其中振幅序列$A_n=\sqrt{a_n^2+b_n^2}$，相位${\phi }_n=-\arctan \frac{b_n}{a_n}$，根据欧拉公式，可以将傅里叶级数写成复指数形式：
$$f(x)=\sum _{n=-\infty }^{\infty }{c}_n{e}^{\frac{2\pi inx}{T}},\text{ 其中 }{c}_k=\frac{1}{T}{\int }_{-\frac{T}{2}}^{\frac{T}{2}}f(t)e^{\frac{-2\pi ikt}{T}}dt$$
傅里叶级数可以逼近任意周期函数，从而得到函数在高维空间基的向量表示，然后根据其维度和向量的膜画出频谱图，而对于非周期函数，可以将周期看作$T=\infty$的周期函数，在上面的傅里叶级数的复指数表达式中，可以将$c_n$看做基坐标,$e^{\frac{2\pi inx}{T}}$看作正交基，周期为无穷时，可以得到：
$$\begin{array}{c}f(x)={\sum }_{n=-\infty }^{\infty }{c}_n\cdot e^{\frac{2\pi inx}{T}}\\ T=\infty \end{array}\}\Rightarrow f(x)={\int }_{-\infty }^{\infty }F(\omega )e^{i\omega x}d\omega$$
上面简化了一下，用$\omega$代表频率。
其中$F(\omega )$：
$$\begin{array}{c}{c}_n=\frac{1}{T}{\int }_{x_0}^{x_0+T}f(x)\cdot e^{-\frac{2\pi inx}{T}}dx\\ T=\infty \end{array}\}\Rightarrow F(\omega )=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }f(x)e^{-i\omega x}dx$$
$F(\omega )$就是傅立叶变换，得到的就是频域曲线。

离散傅里叶变换

从时域的采样数据转变为频域的算法，称为离散傅里叶变换（DFT）。DFT将采样信号的时域和频域联系起来，DFT广泛应用于谱分析、应用力学、光学、医学图像、数据分析、仪器及远程通信等方面。假设有N个时域采样信号，对这N个样本进行DFT变换，结果仍将为N个样本，但它却是频域表示法。时域的N个样本与频域的N个样本之间的关系如下：

假设信号采样率为$f_s$，采样间隔为$\Delta t$，有$\Delta t=1/f_s$，采样信号表示为$x_n,0⩽n⩽N-1$(即有$N$个样本)，对这N个样本进行傅里叶变化，公式如下：
$$X_k=\sum _{n=0}^{N-1}{x}_n{e}^{-j2\pi kn/N}k=0,1,2,\dots ,N-1$$

其结果$X_k$为$x_n$对应的频域表示。注意时域和频域中均有$N$个样本，同时域中的时间间隔对应的频域间隔$\Delta f$为：
$$\Delta f=\frac{f_s}{N}=\frac{1}{N\Delta t}$$

$\Delta f$也被称为频率分辨率，增加采样次数$N$或减小采样频率$f_s$均能减小$\Delta f$（提高分辨率)。

快速傅里叶变换

快速傅里叶变换 (fast Fourier transform), 即利用计算机计算离散傅里叶变换（DFT)的高效、快速计算方法的统称，简称FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT算法计算量的节省就越显著。对于长度为N的输入矢量，FFT是$\mathrm{O}(\mathrm{N}\log \mathrm{N})$级的，而普通DFT算法是$\mathrm{O}({\mathrm{N}}^2)$级的。FFT是怎么快速计算的呢？答案就在于它利用了对称性。从DFT计算公式可看出，不管输入信号$x_n$是实数还是复数，$X_k$总是复数(也可能虚数部分为零),它包含两部分:幅值和相位。可以证明，对于实信号$x_n$，其DFT具有对称的性质：
$$\begin{array}{c}\left|X_k\right|=\left|X_{N-k}\right|\\ \mathrm{phase}\left(X_k\right)=-\mathrm{phase}\left(X_{N-k}\right)\end{array}$$
即幅值偶对称，相位奇对称（下标为k和N-k的两个复数共轭）。
下面的例子演示了这一个规律，先以rand随机产生有8个元素的实数数组x，然后用fft对其运算之后，观察其结果为8个复数：

可以看出第(1、7)，(2、6)，(3、5)个复数互为共轭复数。DFT变换后的复数数组的另一规律是：下标为0和$N/2$的两个复数的虚数部分为0。下面让我们来看看FFT变换之后的那些复数都代表什么意思。

• 首先下标为0的实数（数组中第一个元素）表示了时域信号中的直流成分的多少
• 下标为k的复数a+b*j表示时域信号中周期为N/k个取样值的正弦波和余弦波的成分的多少，其中a表示cos波形的成分，b表示sin波形的成分

我们再来看一下$X_{N+k}$的值，根据DFT公式，可以推出：
$$\begin{array}{rl}{X}_{N+k}& =\sum _{n=0}^{N-1}{x}_n\cdot e^{-i2\pi (N+k)n/N}\\ & =\sum _{n=0}^{N-1}{x}_n\cdot e^{-i2\pi n}\cdot e^{-i2\pi kn/N}\\ & =\sum _{n=0}^{N-1}{x}_n\cdot e^{-i2\pi kn/N}\end{array}$$

根据欧拉公式，对于所有的整数$n$，$e^{-i2\pi n}=1$，则有:$X_{N+k}=X_k$或者$X_{k+i\times N}=X_k$，这体现了变换系数(矩阵)的周期性。利用对称性和周期性，DFT运算中有些项便可以合并，并能将长序列的DFT分解为短序列的DFT进行运算。而DFT的运算量是与$N^2$成正比的，所以$N$越小越有利。
DFT的计算可以分为两部分。设序列长度为N，并且满足N为2的整数次幂。按照n的奇偶性把$x_n$分解为两个$\frac{N}{2}$点的子序列:
$$\{\begin{array}{c}x(2m)=x_1(m)\\ x(2m+1)=x_2(m)m=0,1,2,\dots ,\frac{N}{2}-1\end{array}$$
则从DFT的定义可得：
$$\begin{array}{rl}{X}_k& =\sum _{n=0}^{N-1}{x}_n\cdot e^{-i2\pi kn/N}\\ & =\sum _{m=0}^{N/2-1}{x}_{2m}\cdot e^{-i2\pi k(2m)/N}+\sum _{m=0}^{N/2-1}{x}_{2m+1}\cdot e^{-i2\pi k(2m+1)/N}\\ & =\sum _{m=0}^{N/2-1}{x}_{2m}\cdot e^{-i2\pi km/(N/2)}+e^{-i2\pi k/N}\sum _{m=0}^{N/2-1}{x}_{2m+1}\cdot e^{-i2\pi km/(N/2)}\\ & =\sum _{m=0}^{N/2-1}{x}_1(m)\cdot e^{-i2\pi km/(N/2)}+e^{-i2\pi k/N}\sum _{m=0}^{N/2-1}{x}_2(m)\cdot e^{-i2\pi km/(N/2)}\end{array}$$

从而$N$个点序列的离散傅里叶变换分解为$N/2$个点序列的离散傅里叶变换来实现。用序列$X_1(k)、X_2(k)$分别表示$x_1(m)$和$x_2(m)$的$N/2$点DFT，即：
$$\begin{array}{rl}& X_1(k)=DFT\left[x_1(m)\right]=\sum _{m=0}^{N/2-1}{x}_1(m)e^{-i2\pi km(N/2)},k=0,1,2,\dots ,N/2-1\\ & X_2(k)=DFT\left[x_2(m)\right]=\sum _{m-0}^{N/2-1}{x}_2(m)e^{-i2\pi km(N/2)},k=0,1,2,\dots ,N/2-1\end{array}$$
所以
$$X(k)=X_1(k)+e^{-i2\pi k/N}\cdot X_2(k),k=0,1,2,\dots ,N/2-1(a)$$
利用$e^{-i2\pi (k+N/2)/N}=-e^{-i2\pi k/N}$和序列$X_1(k),X_2(k)$隐含的周期性(以$N/2$为周期)可以得到：
$$X\left(k+\frac{N}{2}\right)=X_1(k)-e^{-i2\pi k/N}\cdot X_2(k),k=0,1,2,\dots ,N/2-1(b)$$
这样将$N$点的序列DFT变换分解为计算两个$N/2$点的DFT变换$X_1(k),X_2(k)$代入式(a)，可以求得$X_{(}k)$前一半$(k=0,1,2,...,N/2-1)$项数的结果，再计算(b)式得到$X_{(}k)$的后一半$(k=N/2,...,N-1)$项数的结果。由于$N$为2的整数次幂，因而$N/2$仍是偶数，可以进一步把每个$N/2$点的子序列再按奇偶性分组为两个$N/4$点的子序列，这样就大大的加快了DFT的计算速度。

总结

本文主要介绍了音频信号预处理中通过傅里叶变换将时域信号转为频域信号，后面将继续学习并更新梅尔频谱的原理。

参考

1. 语音识别——一份简短的技术综述
2. 傅里叶级数和傅里叶变换是什么关系？
3. 快速傅里叶变换(FFT)