Offer 驾到，掘友接招！我正在参与2022春招打卡活动，点击查看活动详情。

语音的声学处理通常称为特征提取或者信号分析，特征是表示语音信号的一个时间片的矢量。常见的特征类型有LPC（线性预测编码）特征和PLP（感知线性预测编码），该特征称为声谱特征，使用形成波形的不同频度的分布来表示这个波形，这种频度分布称为声谱。

1. 基本概念

（1）声波

说话者使用特定的方式通过声门由口腔或鼻腔流出，就造成了空气压力变化，我们将描画空气压力对于时间的变化情况的方法来表示声波。声波在空气中是一种纵波，它的振动方向和传播方向是一致的。

波有两个重要的特征：频度（frequency）和振幅（amplitude）

频度：是1s内波本身重复的次数，或称为‘周’（cycle）

振幅：空气压力变化的大小是振幅，振幅越高，压力越大

人类感知特性：：音高（pitch）是与频率有关，音强是和语音的响度有关。

（参考：语音识别：原理与应用）

（2）声谱

声谱特征是根据Fourier对于复杂波的深入研究提出来的。Fouries指出，每个复杂波都可以表示为很多频度不同的简单波的总和。声谱是一个波的这些不同频度成分的表示。声谱图也称为幅度谱或者幅值谱，横轴是时间，纵轴是频率，一般是指短时傅里叶变换后得到的不含相位信息的时频图，对于梅尔谱变换成为功率谱。

（3）频谱

频谱是频率密度的简称，就是频率的分布曲线。复杂振荡分解为振幅不同和频率不同的谐振荡，这些谐振荡的幅值按频率排列的图形叫频谱。

（4）周期图

周期图是一种信号功率谱密度估计方法。为得到功率谱固估值，先取信号序列的离散傅里叶变换，然后取其幅频特征的平方并除以序列长度N，由于序列的离散傅里叶具有周期性，因而这种功率谱也具有周期性，常称为周期图。

（5）常见指标

过零率（zero crossing rate） 是信号符号变化的比率，即，在每帧中，语音信号从正变为负或从负变为正的次数。过零率越大，频率近似越高。

谱质心（Spectral Centroid） 是描述音色属性的重要物理参数之一，是频率成分的重心，是在一定频率范围内通过能量加权平均的频率，其单位是Hz。它是声音信号的频率分布和能量分布的重要信息。谱质心描述了声音的明亮度，具有阴暗、低沉品质的声音倾向有较多低频内容，谱质心相对较低，具有明亮、欢快品质的多数集中在高频，谱质心相对较高。

声谱衰减是对声音信号形状（波形图）的一种衡量，表示低于总频谱能量的指定百分比的频率。

色度频率是音乐音频有趣且强大的表示，其中整个拼欧被投影到12个区间，代表音乐八度音的12个不同的半音（或色度）

2. 信号处理

输入的声波要进行数字化，就是进行模拟-数字转换，也称为声谱特征的抽取过程，抽取过程从声波本身开始，以特征矢量结束。模拟-数字转换的过程分为两步：采样（sampling）和量化（quantization）。

（1）采样

采样是指对模拟信号在时间域上的离散化过程，即把一个时间上连续、幅度上连续的模拟信号转换成时间上离散、幅度上连续的信号。

采样率是每秒采样点数，通常的采样率是8khz，16khz，采样的方法通常是奈奎斯特（Nyquist）定理，当采样率大于信号中最高频率的两倍时，采样之后的数字信号能够完整保留原始信号中的信息。

（2）量化

量化是将采样之后的实数值进行整数化的过程，量化之后还会有编码过程，常见的编码有PCM编码（脉冲编码调制），MP3编码等。

3.特征提取

语音识别需要进行特征提取，就是提取语音信号中有助于理解语音内容的部分而丢弃其他的东西，常见的语音特征提取方式如下图：

3.1 预处理

预处理过程是对语音信号进行预加重、分帧和加窗处理，目的是对信号进行合理分段。

（1）预加重

因为高频信号的能量通常较低，受到的抑制影响越大，因此需要增加高频部分的能量，即补偿语音信号高频部分的振幅。好处有：增加高频部分的能量使得能量分布更加均衡；防止傅里叶变换的数值计算不稳定；有可能增加信噪比，计算公式为

$y_{t}=x_{t}-\alpha x_{t-1}$

其中 $\alpha$ 是预加重的系数，默认为0.97

（2）分帧

语音信号是一个非平稳信号，但是考虑在发浊音时声带有规律地振动，即基音频率在短时间范围内是相对固定的，语音信号具有短时平稳特性，认为10ms~30ms的语音信号是一个准稳态过程。短时分析主要采用分帧方式，相邻两帧之间的基因可能发生变化，为保证声学特征参数的平滑性，一般采用重叠取帧的方式。

将语音信号进行切割，且保证一定的重复率。核心参数有帧长（每帧的长度）和帧移（移动的距离）。对信号进行分帧，范围是20-40ms，通常采用25ms，考虑到信号的重叠性，每次移动10ms，通常一帧（信号的采样频率是16kHZ）有1600025/1000=400个样本点。移动10ms，则1600010/1000=160个样本点，则每次移动160个样本。

（换算详情：1s=1000ms，1khz=1000hz，赫兹是频率单位），如果最后一帧不够400个样本点，可以考虑用0补齐。

分帧方式相当于对语音信号进行了加矩形窗的处理，易出现频谱泄露， 因此需要加窗函数。

（3）加窗

考虑到语音信号的短时平稳性，对每帧信号进行加窗处理，常见的窗函数有：汉明窗、汉宁窗、布莱克曼窗。公式为

$x_{l}[n]=w[n]x[n+lL]$

其中，w[n]是窗函数，N是窗长，l是帧索引，L是帧移。

3.2 各类特征（更多细节实现看参考资料）

（1）前期资料

美尔尺度（Mel Scale）

美尔尺度是建立从人类的听觉感知的频率到声音实体频率直接的映射。通过把频率转换为美尔尺度，特征能够更好的匹配人类的听觉感知效果。从频率到梅尔频率的转换公式，

从梅尔频率到频率的转换公式为：

梅尔滤波器组

梅尔滤波器的长度与功率谱相等，每个滤波器只有对于需要采集的频率范围是非零，其余都是0。

计算梅尔滤波器组的参数

转换，使用梅尔尺度公式将最小和最大频率转换为美尔尺度的频率，假设最值是300Hz，8kHz，梅尔频率是401.25Mel，2834.99Mel
划分，若我们有10个滤波器，则需要12个点（加上最大、最小频率），并在最值范围内平均分配

将梅尔频率转换为频率，

创建滤波器组

第一个滤波器开始于第一个点（300Hz），峰值在第二个点（517Hz）结束于第三个点（781Hz），第一个滤波器开始于第二个点（517Hz），峰值在第三个点（781Hz）结束于第四个点（1103Hz），公式为

Bark频率

临时频带是在某一段频率范围内纯音和噪声功率相等，则该纯音处于刚好能被听到的临界状态，一个临界频带的宽带被称为一个Bark，Bark频率 z(f) 和线性频率****的对应关系如下：

（2）FBank特征

FBank特征提取的步骤：

对语音信号进行分帧处理
用周期图法进行功率谱估计 （短时傅里叶）
对功率谱用Mel滤波器组进行滤波，计算每个滤波器里的能量（fBank特征）
对每个滤波器的能量取log

（3）MFCC特征

MFCC特征提取的步骤：

对语音信号进行分帧处理
用周期图法进行功率谱估计 （短时傅里叶）
对功率谱用Mel滤波器组进行滤波，计算每个滤波器里的能量（fBank特征）
对每个滤波器的能量取log
进行离散余弦变换(DCT)
保留DCT的第2-13个系数，去掉其他

（3）PLP特征

对语音信号进行分帧处理
用周期图法进行功率谱估计 （短时傅里叶）
临界频带分析（Bark滤波器）
等响度预加重
强度响度转换
逆傅里叶变换
线性预测,通过使线性预测到的采样点值在最小均方误差约束下逼近实际语音采样点值，可以求取一组唯一的预测系数，简称为线性预测编码（Linear prediction coding，LPC）

（4）区别与联系

滤波器：MFCC、FBank是用梅尔滤波器，PLP是Bark滤波器。

STFT：三者都采用短时傅里叶变换。

优点：FBank保留更多原始特征，适用于DNN建模；MFCC去相关性较好，适合GMM建模；PLP抗燥性更强。

4.librosa中的实现

4.1 流程图

4.2 相关函数的参数

（1）MFCC参数

4.3 实现细节

（1）MFCC函数

def mfcc(
    y=None, sr=22050, S=None, n_mfcc=20, dct_type=2, norm="ortho", lifter=0, **kwargs
):

    if S is None:
        S = power_to_db(melspectrogram(y=y, sr=sr, **kwargs))

    M = scipy.fftpack.dct(S, axis=0, type=dct_type, norm=norm)[:n_mfcc]

    if lifter > 0:
        M *= (
            1
            + (lifter / 2)
            * np.sin(np.pi * np.arange(1, 1 + n_mfcc, dtype=M.dtype) / lifter)[
                :, np.newaxis
            ]
        )
        return M
    elif lifter == 0:
        return M
    else:
        raise ParameterError(
            "MFCC lifter={} must be a non-negative number".format(lifter)
        )

（2）功率谱计算（melspectrogram）+滤波

该函数里边包含了功率谱计算和Mel滤波器进行滤波的过程，其中_spectrogram是进行功率谱计算，filters.mel则建立滤波器。

def melspectrogram(y=None,sr=22050,S=None,n_fft=2048,hop_length=512,win_length=None,window="hann",center=True,pad_mode="reflect",power=2.0,**kwargs):
    S, n_fft = _spectrogram(
        y=y,
        S=S,
        n_fft=n_fft,
        hop_length=hop_length,
        power=power,
        win_length=win_length,
        window=window,
        center=center,
        pad_mode=pad_mode,
    )

    # Build a Mel filter
    mel_basis = filters.mel(sr, n_fft, **kwargs)

    return np.dot(mel_basis, S)

功率谱计算主要是用短时傅里叶变换stft计算：

def _spectrogram(y=None,S=None,n_fft=2048,hop_length=512,power=1,win_length=None,window="hann",center=True,pad_mode="reflect"):

    if S is not None:
        # Infer n_fft from spectrogram shape
        n_fft = 2 * (S.shape[0] - 1)
    else:
        # Otherwise, compute a magnitude spectrogram from input
        S = (
            np.abs(
                stft(
                    y,
                    n_fft=n_fft,
                    hop_length=hop_length,
                    win_length=win_length,
                    center=center,
                    window=window,
                    pad_mode=pad_mode,
                )
            )
            ** power
        )

    return S, n_fft

mel滤波器是通过filter.mel进行滤波器建立：

def mel(sr,n_fft,n_mels=128,fmin=0.0,fmax=None,htk=False,norm="slaney",dtype=np.float32):

    if fmax is None:
        fmax = float(sr) / 2

    # Initialize the weights
    n_mels = int(n_mels)
    weights = np.zeros((n_mels, int(1 + n_fft // 2)), dtype=dtype)

    # Center freqs of each FFT bin
    fftfreqs = fft_frequencies(sr=sr, n_fft=n_fft)

    # 'Center freqs' of mel bands - uniformly spaced between limits
    mel_f = mel_frequencies(n_mels + 2, fmin=fmin, fmax=fmax, htk=htk)

    fdiff = np.diff(mel_f)
    ramps = np.subtract.outer(mel_f, fftfreqs)

    for i in range(n_mels):
        # lower and upper slopes for all bins
        lower = -ramps[i] / fdiff[i]
        upper = ramps[i + 2] / fdiff[i + 1]

        # .. then intersect them with each other and zero
        weights[i] = np.maximum(0, np.minimum(lower, upper))

    if norm == "slaney":
        # Slaney-style mel is scaled to be approx constant energy per channel
        enorm = 2.0 / (mel_f[2 : n_mels + 2] - mel_f[:n_mels])
        weights *= enorm[:, np.newaxis]
    else:
        weights = util.normalize(weights, norm=norm, axis=-1)

    # Only check weights if f_mel[0] is positive
    if not np.all((mel_f[:-2] == 0) | (weights.max(axis=1) > 0)):
        # This means we have an empty channel somewhere
        warnings.warn(
            "Empty filters detected in mel frequency basis. "
            "Some channels will produce empty responses. "
            "Try increasing your sampling rate (and fmax) or "
            "reducing n_mels."
        )

    return weights

5. 实例

import librosa
import librosa.display as display

v = librosa.load('./english.wav')

mfcc = librosa.feature.mfcc(v[0])

import matplotlib.pyplot as plt

display.specshow(mfcc)
plt.colorbar()
# display.waveshow(v[0])