（超详细）语音信号处理之特征提取_语音特征提取

语音信号处理之特征提取

语音信号处理之特征提取要对语音信号进行分析，首先要分析并提取出可表示该语音本质的特征参数。有了特征参数才能利用这些特征参数进行有效的处理。

根据提取参数的方法不同，可将语音信号分析分为时域，频域，倒频域，和其他域的分析方法。根据分析方法的不同，可将语音信号分析分为模型分析方法和非模型分析方法。

本文主要以第一种分类方法。时域分析方法简单，计算量小，物理意义明确，但由于语音信号最最重要的感知特性在功率谱中，而人耳对相位变化并不敏感，所以频域分析更为重要。

通过最基本的特征，后面针对不同的任务演变出了各种特征。

语音数据

时域图就是语音波形，横轴时间，纵轴幅值，下是用Adobe Audition打开的音频的波形图，表示这段语音波形时采样率16000Hz，量化精度是16bit。

从图中可以得到各个音的起始位置，但很难看出更加有用的信息。但如果我们将其放大到100ms的场景下，可以得到下图所示。

可以看出在短时内波形是存在一定的周期的，不同的发音往往对应着不同的周期的变化，因此在短时域上我们可以将波形通过傅里叶变换转化为频域图，观察音频的周期特性，从而获取有用的音频特征。

额外补充（不感兴趣可跳过）
截取一段清音/s/：确实好像没有周期性规律

清音一般用白噪声建模：
白噪声（white noise）是指功率谱密度在整个频域内是常数的噪声。 所有频率具有相同能量密度的随机噪声称为白噪声。

再截取一段浊音/o/：短时内波形是存在一定的周期的

预处理

1. 分帧

为什么要分帧？
语音信号是非线性，非平稳，时变的随机信号。但短时内可看作是平稳时不变的，这样就有了分帧技术。

由于语音信号具有短时平稳性，我们通常对语音进行分帧加窗处理，截取短时音频片段，通常帧移的大小为5-10ms，窗长大小通常为帧移的2-3倍即20-30ms。

帧长设置有什么依据呢？
从宏观上看，它必须足够短来保证帧内信号是平稳时不变的。口型的变化是导致信号不平稳的原因，所以在一帧的期间内口型不能有明显变化，即一帧的长度应当小于一个音素的长度。正常语速下，音素的持续时间大约是 50~200 毫秒，所以帧长一般取为小于 50 毫秒。
从微观上来看，它又必须包括足够的振动周期，即包含若干个语音的基频，男声基频在 100 Hz左右，女声在 200Hz左右，换算成周期就是 10 ms和 5 ms，所以一般取至少 20 ms。

python代码：

def framesig(sig,frame_len,frame_step):
    """Frame a signal into overlapping frames.

    :param sig: the audio signal to frame.
    :param frame_len: length of each frame measured in samples.
    :param frame_step: number of samples after the start of the previous frame that the next frame should begin.
    :returns: an array of frames. Size is NUMFRAMES by frame_len.
    """
    slen = len(sig)
    frame_len = int(round_half_up(frame_len))
    frame_step = int(round_half_up(frame_step))
    if slen <= frame_len:
        numframes = 1
    else:
        numframes = 1 + int(math.ceil((1.0*slen - frame_len)/frame_step))

    padlen = int((numframes-1)*frame_step + frame_len)

    zeros = numpy.zeros((padlen - slen,))
    padsignal = numpy.concatenate((sig,zeros))

    indices = numpy.tile(numpy.arange(0,frame_len),(numframes,1)) + numpy.tile(numpy.arange(0,numframes*frame_step,frame_step),(frame_len,1)).T
    indices = numpy.array(indices,dtype=numpy.int32)
    frames = padsignal[indices]
    return frames

frames = framesig(sig=sig, frame_len=0.030 * sample_rate, # 取帧长为30ms
                                          frame_step=0.006 * sample_rate, # 取帧移为6ms
                                          )
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

2. 加窗

为什么要加窗？
加窗的目的有两个：频谱泄露和栅栏效应。
即使语音短时内是准周期的，但分帧时很难保证每帧内就是完全周期的。非周期信号被截断后，经过傅里叶变换FFT，会出现原本信号并没有的频率，即多余的频率成分，这种现象叫做频谱泄露。通过加上合适的窗，可以抑制这些成分，即高次谐波，使得结果更圆滑一些。
对信号做FFT时，得到的是一系列离散的谱线，如果信号中的频率成分位于谱线之间，而不是谱线上，那么就会造成幅值和香味上的偏差，这就是栅栏效应。
理论上，这两种误差都无法消除，但可以选择适合的窗函数进行抑制。

为什么要重叠帧移？
加窗的代价是一帧信号两端的部分被削弱了，没有像中央的部分那样得到重视。弥补的办法是，帧不要背靠背地截取，而是相互重叠一部分。相邻两帧的起始位置的时间差叫做帧移。

加窗操作如下图所示：

一般有矩形窗、汉明（Hamming）窗、汉宁（Hanning）窗等。
定义矩形窗为w(m)：

那么对于语音信号 x(t)，其加窗分帧后第n帧语音信号 xn(m)为：

在该计算式中，n=0,T,2T,…,N为帧长，T为帧移长度。

加窗分帧总体：

加窗就是用一定的窗函数w(n)来乘s(n)，从而得到加窗语音信号
x(n)=s(n)∗w(n)。

python代码：

def framesig(sig,frame_len,frame_step,winfunc=lambda x:numpy.ones((x,))):
    """Frame a signal into overlapping frames.

    :param sig: the audio signal to frame.
    :param frame_len: length of each frame measured in samples.
    :param frame_step: number of samples after the start of the previous frame that the next frame should begin.
    :param winfunc: the analysis window to apply to each frame. By default no window is applied.
    :returns: an array of frames. Size is NUMFRAMES by frame_len.
    """
    slen = len(sig)
    frame_len = int(round_half_up(frame_len))
    frame_step = int(round_half_up(frame_step))
    if slen <= frame_len:
        numframes = 1
    else:
        numframes = 1 + int(math.ceil((1.0*slen - frame_len)/frame_step))

    padlen = int((numframes-1)*frame_step + frame_len)

    zeros = numpy.zeros((padlen - slen,))
    padsignal = numpy.concatenate((sig,zeros))

    indices = numpy.tile(numpy.arange(0,frame_len),(numframes,1)) + numpy.tile(numpy.arange(0,numframes*frame_step,frame_step),(frame_len,1)).T
    indices = numpy.array(indices,dtype=numpy.int32)
    frames = padsignal[indices]
    # 加窗操作
    win = numpy.tile(winfunc(frame_len),(numframes,1))
    return frames*win

frames = framesig(sig=sig, frame_len=0.030 * sample_rate, # 取帧长为30ms
                  frame_step=0.006 * sample_rate, # 取帧移为6ms
                  winfunc=np.hamming
                  )
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

3. 预加重

为什么要预加重？
语音经过说话人的口唇辐射发出，受到唇端辐射抑制，高频能量明显降低。一般来说，当语音信号的频率提高两倍时，其功率谱的幅度下降约6dB，即语音信号的高频部分受到的抑制影响较大。比如像元音等一些因素的发音包含了较多的高频信号的成分，高频信号的丢失，可能会导致音素的共振峰并不明显，使得声学模型对这些音素的建模能力不强。
具体做法：
预加重（pre-emphasis）是个一阶高通滤波器，可以提高信号高频部分的能量，给定时域输入信号x[n]，预加重之后信号为：

其中，a是预加重系数，一般取0.97或0.95。如下图4所示，元音音素 /aa/ 原始的频谱图（左）和经过预加重之后的频谱图（右）。

python代码：

def preemphasis(y, coef=0.97, zi=None, return_zf=False):
    """Pre-emphasize an audio signal with a first-order auto-regressive filter:

        y[n] -> y[n] - coef * y[n-1]


    Parameters
    ----------
    y : np.ndarray
        Audio signal

    coef : positive number
        Pre-emphasis coefficient.  Typical values of ``coef`` are between 0 and 1.

        At the limit ``coef=0``, the signal is unchanged.

        At ``coef=1``, the result is the first-order difference of the signal.

        The default (0.97) matches the pre-emphasis filter used in the HTK
        implementation of MFCCs [#]_.

        .. [#] http://htk.eng.cam.ac.uk/

    zi : number
        Initial filter state.  When making successive calls to non-overlapping
        frames, this can be set to the ``zf`` returned from the previous call.
        (See example below.)

        By default ``zi`` is initialized as ``2*y[0] - y[1]``.

    return_zf : boolean
        If ``True``, return the final filter state.
        If ``False``, only return the pre-emphasized signal.

    Returns
    -------
    y_out : np.ndarray
        pre-emphasized signal

    zf : number
        if ``return_zf=True``, the final filter state is also returned
    """

    b = np.asarray([1.0, -coef], dtype=y.dtype)
    a = np.asarray([1.0], dtype=y.dtype)

    if zi is None:
        # Initialize the filter to implement linear extrapolation
        zi = 2 * y[..., 0] - y[..., 1]

    zi = np.atleast_1d(zi)

    y_out, z_f = scipy.signal.lfilter(b, a, y, zi=np.asarray(zi, dtype=y.dtype))

    if return_zf:
        return y_out, z_f

    return y_out


wave_data, self.sr = librosa.load(input_file, sr=sr)  # 音频全部采样点的归一化数组形式数据
wave_data = preemphasis(wave_data, coef=preemph)  # 预加重，系数0.97
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62

时域特征

1.短时能量、短时平均幅度

短时能量通常指的是，时域上每帧能量的幅度的平方。
第n 帧语音信号 xn(m)的短时能量En为：

使用幅值平方将对高幅值信号具有较大的敏感度，为了降低敏感度，定义短时平均幅度函数Mn为：

短时能量En和短时平均幅度函数Mn的用途：

• 区分浊音和清音。浊音相比较于清音的En具有较大的数值，因而可用于区分浊音和清音。

• 区分有声段和无声段，也可对声母和韵母分界，对无间隙的连字分界。如在简单的语音活动检测（Voice Activity Detection，VAD）中，直接利用能量特征：能量大的音频片段是语音，能量小的音频片段是非语音（包括噪音、静音段等）。这种VAD的局限性比较大，正确率也不高，对噪音非常敏感。
  例：语音“蓝天白云”
  

• 在语音识别任务中作为特征，表示能量特征和超音频信息。

python代码：

def __init__(self, input_file, sr=None, frame_len=512, n_fft=None, win_step=2 / 3, window="hamming"):
        """
        初始化
        :param input_file: 输入音频文件
        :param sr: 所输入音频文件的采样率，默认为None
        :param frame_len: 帧长，默认512个采样点(32ms,16kHz),与窗长相同
        :param n_fft: FFT窗口的长度，默认与窗长相同
        :param win_step: 窗移，默认移动2/3，512*2/3=341个采样点(21ms,16kHz)
        :param window: 窗类型，默认汉明窗
        """
        self.input_file = input_file
        self.frame_len = frame_len  # 帧长，单位采样点数
        self.wave_data, self.sr = librosa.load(self.input_file, sr=sr)
        self.window_len = frame_len  # 窗长512
        if n_fft is None:
            self.fft_num = self.window_len  # 设置NFFT点数与窗长相等
        else:
            self.fft_num = n_fft
        self.win_step = win_step
        self.hop_length = round(self.window_len * win_step)  # 重叠部分采样点数设置为窗长的1/3（1/3~1/2）,即帧移(窗移)2/3
        self.window = window

def energy(self):
        """
        每帧内所有采样点的幅值平方和作为能量值
        :return: 每帧能量值，np.ndarray[shape=(1，n_frames), dtype=float64]
        """
        mag_spec = np.abs(librosa.stft(self.wave_data, n_fft=self.fft_num, hop_length=self.hop_length,
                                       win_length=self.frame_len, window=self.window))
        pow_spec = np.square(mag_spec) # [frequency, time (n_frames)]
        energy = np.sum(pow_spec, axis=0) # [n_frames]
        energy = np.where(energy == 0, np.finfo(np.float64).eps, energy)  # 避免能量值为0，防止后续取log出错(eps是取非负的最小值), 即np.finfo(np.float64).eps = 2.220446049250313e-16
        return energy

def short_time_energy(self):
        """
        计算语音短时能量：每一帧中所有语音信号的平方和
        :return: 语音短时能量列表(值范围0-每帧归一化后能量平方和，这里帧长512，则最大值为512)，
        np.ndarray[shape=(1，无加窗，帧移为0的n_frames), dtype=float64]
        """
        energy = []  # 语音短时能量列表
        energy_sum_per_frame = 0  # 每一帧短时能量累加和
        for i in range(len(self.wave_data)):  # 遍历每一个采样点数据
            energy_sum_per_frame += self.wave_data[i] ** 2  # 求语音信号能量的平方和
            if (i + 1) % self.frame_len == 0:  # 一帧所有采样点遍历结束
                energy.append(energy_sum_per_frame)  # 加入短时能量列表
                energy_sum_per_frame = 0  # 清空和
            elif i == len(self.wave_data) - 1:  # 不满一帧，最后一个采样点
                energy.append(energy_sum_per_frame)  # 将最后一帧短时能量加入列表
        energy = np.array(energy)
        energy = np.where(energy == 0, np.finfo(np.float64).eps, energy)  # 避免能量值为0，防止后续取log出错(eps是取非负的最小值)
        return energy
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

2. 短时过零率

短时过零率表示一帧语音中波形信号穿过零值的次数。对于连续信号，过零意味着波形通过时间轴，而对于离散信号，过零意味着相邻采样点的符号改变。

首先定义符号函数sgn[·]为：

则第n帧语音信号 xn(m)的短时过零率Zn为：

由于短时过零率容易受到低频干扰，可设置相关门限T，将过零修改为穿过正负门限的次数，即：

门限的存在使得短时过零率Zn具有一定的扛干扰能力，避免随机噪声导致的虚假过零。

短时过零率的用途：

• 浊音能量集中于3kHz内的低频率段，清音能量集中于高频率段，而短时过零率可以一定程度反映频率高低，因而浊音段相对于清音段，其短时过零率减低。

• 将短时过零率和短时能量结合实现端点检查。短时能量适用于背景噪声较小的情况，而短时过零率适用于背景噪声较大的情况。实际中，通常结合两个参数实现语音起点和终点的判断。

语音“我回来了”

python代码：

def zero_crossing_rate(self):
        """
        计算语音短时过零率：单位时间(每帧)穿过横轴（过零）的次数
        :return: 每帧过零率次数列表，np.ndarray[shape=(1，无加窗，帧移为0的n_frames), dtype=uint32]
        """
        zcr = []  # 语音短时过零率列表
        counting_sum_per_frame = 0  # 每一帧过零次数累加和，即过零率
        for i in range(len(self.wave_data)):  # 遍历每一个采样点数据
            if i % self.frame_len == 0:  # 开头采样点无过零，因此每一帧的第一个采样点跳过
                continue
            if self.wave_data[i] * self.wave_data[i - 1] < 0:  # 相邻两个采样点乘积小于0，则说明穿过横轴
                counting_sum_per_frame += 1  # 过零次数加一
            if (i + 1) % self.frame_len == 0:  # 一帧所有采样点遍历结束
                zcr.append(counting_sum_per_frame)  # 加入短时过零率列表
                counting_sum_per_frame = 0  # 清空和
            elif i == len(self.wave_data) - 1:  # 不满一帧，最后一个采样点
                zcr.append(counting_sum_per_frame)  # 将最后一帧短时过零率加入列表
        return np.array(zcr, dtype=np.uint32)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

3.短时自相关函数

语音信号xn(m)的短时自相关函数Rn(k)为：

其中，若信号xn(m)具有周期性，则短时自相关函数Rn(k)也具有周期性，且两者周期相同；Rn(k)为偶函数，当k=0时，自相关函数具有最大值。

假设语音信号xn(m)的周期为T，那么短时自相关函数Rn(k)将在k=T,2T…取值时出现峰值。若要出现第一个峰值（即k=T)，根据下式：

需要取到信号中x(m=2T的样本点，即语音帧宽至少应大于两个周期，否则第一个峰值将无法较好的显示。例语音最小基频为80Hz，最大周期为12.5ms，两倍周期为25ms，因此10kHz的采样信号的帧宽至少为250个采样点。

另一方面，考虑到语音信号的短时性，应设置较低的帧长，因此可使用修正短时自相关函数，其定义为：

注意：

相比于短时自相关函数Rn(k)，在修正短时自相关函数中，第一项xn(m)与Rn(k)中的xn(m)相同，而第二项x’n(m)与Rn(k)中的xn(m)相比，差异在于额外向后包括了k个样本点。

在严格定义中，修正短时自相关函数是一个互相关函数，其不满足自相关函数的性质（偶函数性），但其仍在周期整数倍上具有峰值。

短时自相关函数的用途：

• 浊音的自相关函数具有周期性，而清音的自相关函数类似于高频白噪声，没有周期性。
• 根据自相关函数的第一个峰值的位置，估算浊音的基音频率。

4.短时平均幅度差函数

短时自相关函数使用大量乘法运算，计算时间较长，短时平均幅度差Fn(k)使用减法代替了乘法，大大减少了运算量，大量运用于实时语音处理方案上，其定义为：

对于周期为T的语音信号，短时平均幅度差Fn(k)在k=T,2T…等取值上具有周期性的极小值。类似的，修正短时平均幅度差为：

短时平均幅度差Fn(k)和Rn(k)具有数值关系：

其中，β(k)对不同的语音段，其数值在0.6-1.0之间变化。

短时平均幅度差的用途：

• 基音周期的检测，该方法比短时自相关方法的计算更为简单。

频域特征

1.短时傅里叶变换

短时傅里叶变换（STFT）是最经典的时频域分析方法。所谓短时傅里叶变换，顾名思义，是对短时的信号做傅里叶变化。由于语音波形只有在短时域上才呈现一定周期性，因此使用的短时傅里叶变换可以更为准确的观察语音在频域上的变化。傅里叶变化实现的示意图如下：

一般实现都用快速傅里叶变化FFT。
这里音频经过快速傅里叶变换返回的是复数，其中实部表示的频率的振幅，虚部表示的是频率的相位。
包含FFT函数的库有很多，简单列举几个：

python代码：

import librosa
import torch
import scipy

x_stft = librosa.stft(wav, n_fft=fft_size, hop_length=hop_size,win_length=win_length)
x_stft = torch.stft(wav, n_fft=fft_size, hop_length=hop_size, win_length=win_size)
x_stft = scipy.fftpack.fft(wav)
1
2
3
4
5
6
7

2.语谱图

也叫声谱图，就是声音的可视化图。先将语音信号作FFT，然后以横轴为时间，纵轴为频率，用颜色表示幅值即可绘制出语谱图。在一幅图中表示信号的频率、幅度随时间的变化，故也称“时频图”。

详细原理看我博客：声谱图原理

3. 短时功率谱密度

信号经过FFT得到了频谱，而要反映信号的功率就要求信号的功率谱密度PSD。
功率谱定义：

可见对于有限的信号，功率谱之所以可以估计，是基于两点假设：1）信号平稳; 2)随机信号具有遍历性。
相关函数与功率谱密度是互为傅里叶变换。

基频

基频（fundamental frequency）是声带振动周期的倒数。比如语音中的元音和浊辅音，发音时喉咙是振动的。
那音高（pitch）和基频什么关系呢？它们指的是一回事，但音高是从听觉感知层面出发，主要出现在语音学、音系学等领域；基频是从物理层面出发，主要出现在信号处理、计算机等领域。

下图是语音：“巴拔把爸”的波形（上）和语谱图（下）
基频就是语谱图中最下面的一条线（红色箭头所指），其他的线叫作谐波，或谐频。谐波是基频对应的整数次频率成分，人们对音高的感知取决于第一谐波数值减去第二谐波数值，奇数谐波的能量下降，iu会影响到音高的感知。

基频也是蓝色的线，注意：语谱图的坐标在左边是0-4000Hz，基频的左边在右边0-800Hz。

基频的计算方法看我其他博客。
这里列举一个librosa库的方法。

def pitch(self, ts_mag=0.25):
        """
        获取每帧音高，即基频，这里应该包括基频和各次谐波，最小的为基频（一次谐波），其他的依次为二次、三次...谐波
        各次谐波等于基频的对应倍数，因此基频也等于各次谐波除以对应的次数，精确些等于所有谐波之和除以谐波次数之和
        :param ts_mag: 幅值倍乘因子阈值，>0，大于np.average(np.nonzero(magnitudes)) * ts_mag则认为对应的音高有效,默认0.25
        :return: 每帧基频及其对应峰的幅值(>0)，
                 np.ndarray[shape=(1 + n_fft/2，n_frames), dtype=float32]，（257，全部采样点数/(512*2/3)+1）
        """
        mag_spec = np.abs(librosa.stft(self.wave_data, n_fft=self.fft_num, hop_length=self.hop_length,
                                       win_length=self.frame_len, window=self.window))
        pitches, magnitudes = librosa.piptrack(S=mag_spec, sr=self.sr, threshold=1.0, ref=np.mean,
                                               fmin=50, fmax=500)  # 人类正常说话基频最大可能范围50-500Hz
        ts = np.average(magnitudes[np.nonzero(magnitudes)]) * ts_mag
        pit_likely = pitches
        mag_likely = magnitudes
        pit_likely[magnitudes < ts] = 0
        mag_likely[magnitudes < ts] = 0
        return pit_likely, mag_likely
pitches, mags = self.pitch()  # 获取每帧基频
f0_likely = []  # 可能的基频F0
for i in range(pitches.shape[1]):  # 按列遍历非0最小值，作为每帧可能的F0
    try:
        f0_likely.append(np.min(pitches[np.nonzero(pitches[:, i]), i]))
    except ValueError:
        f0_likely.append(np.nan)  # 当一列，即一帧全为0时，赋值最小值为nan
f0_all = np.array(f0_likely)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

共振峰

共振峰（formant）：声门处的准周期激励进入声道时会引起共振特性，产生一组共振频率，这一组共振频率称为共振峰频率或简称共振峰。
共振峰包含在语音的频谱包络中，频谱包络的极大值就是共振峰。第一个峰就是第一共振峰，元音一般取前3-5个共振峰表示。

下图是语音：“巴拔把爸”的波形（上）和语谱图（下）
共振峰就是语谱图中的红色点线，因为坐标在0-4000Hz,大概只包含了前三个共振峰。

声强和声强级（声压和声压级）

单位时间内通过垂直于声波传播方向的单位面积的平均声能，称作声强，声强用P表示，单位为“瓦/平米”。实验研究表明，人对声音的强弱感觉并不是与声强成正比，而是与其对数成正比，所以一般声强用声强级来表示：

python代码：

def intensity(self):
        """
        计算声音强度，用声压级表示：每帧语音在空气中的声压级Sound Pressure Level(SPL)，单位dB
        公式：20*lg(P/Pref)，P为声压（Pa），Pref为参考压力(听力阈值压力)，一般为1.0*10-6 Pa
        这里P认定为声音的幅值：求得每帧所有幅值平方和均值，除以Pref平方，再取10倍lg
        :return: 每帧声压级，dB，np.ndarray[shape=(1，无加窗，帧移为0的n_frames), dtype=float64]
        """
        p0 = 1.0e-6  # 听觉阈限压力auditory threshold pressure: 2.0*10-5 Pa
        e = self.short_time_energy()
        spl = 10 * np.log10(1 / (np.power(p0, 2) * self.frame_len) * e)
        return spl
1
2
3
4
5
6
7
8
9
10
11