滤波器组FBanks特征 & 梅尔频率倒谱系数MFCC基于librosa, torchaudio

---

滤波器组(Filter Banks, FBanks)特征 & 梅尔频率倒谱系数(Mel Frequency Cepstral Coefficients, MFCC) 基于librosa, torchaudio

• 说明：FBanks & MFCC作为特征被广泛应用于语音识别领域。本文将使用librosa和torchaudio分别实现。计算流程如下图所示（此处暂不涉及PLP）。如有错误，欢迎指正。
• 目的：尽可能简单地理解与使用MFCC，以快速上手语音识别任务。
  

1. 准备工作：导入librosa，torchaudio等包及数据

• 本人使用版本说明：
  $\circ$python==3.8.9
  $\circ$librosa==0.9.2
  $\circ$torch==1.10.1+cu113
  $\circ$torchaudio==0.10.1+cu113
  $\circ$scipy==1.8.0
• 示例音频数据：
  $\circ$下载地址: Lab41-SRI-VOiCES-src-sp0307-ch127535-sg0042.wav
  $\circ$采样频率: 16k
• 参数说明：
  $\circ$sample_rate (sr): 采样频率，即每秒采样点数，单位为赫兹(Hz)。

import torch
import torchaudio
from torch import Tensor
import numpy as np
import librosa
from librosa.display import specshow
from matplotlib import pyplot as plt
from scipy.fftpack import dct

wave_pt, sample_rate_pt = torchaudio.load('Lab41-SRI-VOiCES-src-sp0307-ch127535-sg0042.wav')
# librosa的默认采样频率是22050，设None即可读取原始sr
wave_librosa, sample_rate_librosa = librosa.load('Lab41-SRI-VOiCES-src-sp0307-ch127535-sg0042.wav', sr=None)

print('torchaudio版: ', wave_pt.shape, '|', sample_rate_pt)
print('librosa版: ', wave_librosa.shape, '|', sample_rate_librosa)

torchaudio版:  torch.Size([1, 54400]) | 16000
librosa版:  (54400,) | 16000

time_axis = torch.arange(0, wave_pt.shape[1]) / sample_rate_pt
plt.figure(figsize=(15, 10))
plt.subplot(2, 1, 1)
plt.plot(time_axis, wave_pt.squeeze(), 'red')
plt.title('origin wave_pt')
plt.ylabel('Amplitude')
plt.xlabel('Time(s)')
plt.subplot(2, 1, 2)
plt.plot(time_axis, wave_librosa, 'blue')
plt.title('origin wave_librosa')
plt.ylabel('Amplitude')
plt.xlabel('Time(s)')
plt.show()

2. 预处理

2.1 分帧 (Framing)

• 说明：大多数情况下，对整个信号进行傅里叶变换是没有意义的，因为我们会随着时间的推移丢失信号的频率轮廓。假定信号在短时间($20ms\sim 40ms$)内是稳定的，并将这段时间的信号设为1帧，相邻的帧之间有50%±10%的重叠。通常帧长为$25ms$，帧移为$10ms$。

• 参数：
  $\circ$frame_size: 0.025
  $\circ$frame_stride: 0.01

2.2 加窗(Window)

• 说明：分帧的相当于对信号加了矩形窗，相应的频域会存在严重的频谱泄露，为了减少频谱泄露，对每帧信号加其他形式的窗。

$\circ$汉明窗(Hamming):
$\diamond$公式：
$${\mathrm{w}}_{\text{ham }}[n]=0.54-0.46\cos \left(\frac{2\pi n}{N}-1\right)$$
$\diamond$图形：

$\circ$汉宁窗(Hanning):
$\diamond$公式：
$${\mathrm{W}}_{han}[n]=0.5\left[1-\cos \left(\frac{2\pi n}{N}-1\right)\right]$$
$\diamond$图形：

3. 短时傅立叶变换(Short Time Fourier Transform, STFT)

• 说明：为了分离不同频率的信号，傅里叶变换的意义：任何连续测量的时序或信号，都可以表示为不同频率的正弦波信号的无限叠加。
  

• 计算公式：

$\circ$欧拉公式：$${\mathrm{e}}^{\mathrm{j}\omega \mathrm{n}}=\cos (\omega \mathrm{n})+\mathrm{j}\sin (\omega \mathrm{n})$$

$\circ$离散傅里叶变换(Discrete Fourier Transform, DFT)第$k$个点的计算方式：
$$\mathrm{X}[k]=\sum _{n=0}^{N-1}x[n]{\mathrm{e}}^{-\frac{j2\pi kn}{K}}=\sum _{n=0}^{N-1}x[n]{\mathrm{e}}^{-jwn}$$
$\diamond$$\omega =\frac{2\pi k}{K}\in [0,2\pi ]$为角频率；
$\diamond$$x[n]$是时域波形第$n$个采样点；
$\diamond$$X[k]$是傅里叶频谱第$k$个点；
$\diamond$$N$是采样的点数，$K$是DFT的大小，其中$K\ge N$。

$\circ$利用欧拉公式：$${\mathrm{e}}^{-\frac{j2\pi kn}{K}}=\cos \left(\frac{2\pi kn}{K}\right)-\mathrm{j}\sin \left(\frac{2\pi kn}{K}\right)$$ \qquad\qquad\qquad 则：$$\mathrm{X}[k]={\mathrm{X}}_{\mathrm{r}\mathrm{e}\mathrm{a}\mathrm{l}}[k]-\mathrm{j}{\mathrm{X}}_{\mathrm{i}\mathrm{m}\mathrm{a}\mathrm{g}}[k]$$ $${\mathrm{X}}_{\text{real }}[k]=\sum _{n=0}^{N-1}x[n]\cos \left(\frac{2\pi kn}{K}\right)$$ $${\mathrm{X}}_{\text{imag }}[k]=\sum _{n=0}^{N-1}x[n]\sin \left(\frac{2\pi kn}{K}\right)$$

$\circ$振幅频谱(Magnitude Spectrum)：$${\mathrm{X}}_{\text{magnitude }}[k]=\sqrt{\left({\mathrm{X}}_{\text{real }}[k{]}^2+{\mathrm{X}}_{\text{imag }}[k{]}^2\right)}$$
$\circ$能量频谱(Power Spectrum)：$${\mathrm{X}}_{\text{power }}[k]={\mathrm{X}}_{\text{real }}[k{]}^2+{\mathrm{X}}_{\text{imag }}[k{]}^2$$
$\circ$分贝(Decibel)：$$N_{\mathrm{d}\mathrm{B}}[k]=10\lg \frac{{\mathrm{X}}_{\text{power }}[k]}{ref}$$

• 代码：torch.stft及librosa.stft

$\circ$参数说明：
$\diamond$n_fft: FFT窗口长度，则对应的分帧时间为$\frac{n_{fft}}{sample\_rate}$秒，默认为2048。
$\diamond$hop_length: 帧移，如果未指定，则默认$wi{n}_{length}//4$。
$\diamond$win_length: 每一帧音频都由窗函数加窗，窗长$wi{n}_{length}$，且$wi{n}_{length}<=n_{fft}$，用零填充以匹配$n_{fft}$，默认$wi{n}_{length}=n_{fft}$。
$\diamond$window: 窗函数，类型为string, tuple, number, function, shape=(n_fft,)的向量，默认为'hann'。
$\diamond$center: 如果为True，则假定输出复数矩阵D具有居中的帧，如果False，则假定D具有左对齐的帧，默认False

$\circ$输出：
$\diamond$D: 复数矩阵，shape$=(n_{channels},1+\frac{n_{fft}}{2},n_{frame})$, $n_{frame}=int(1+(y.shape[-1]-wi{n}_{length})//ho{p}_{length})$

d_pt = torch.stft(
    wave_pt,
    n_fft=2048,
    hop_length=None,
    win_length=None,
    window=torch.hann_window(window_length=2048),
    center=False,
    pad_mode="constant",
    return_complex=True
)

d_librosa = librosa.stft(
    wave_librosa,
    n_fft=2048,
    hop_length=None,
    win_length=None,
    window="hann",
    center=False,
    dtype=None,
    pad_mode="constant"
)
print('torch版复数矩阵: ', d_pt.shape)
print('librosa版复数矩阵: ', d_librosa.shape)

torch版复数矩阵:  torch.Size([1, 1025, 103])
librosa版复数矩阵:  (1025, 103)

def power_to_db(power_spec, ref=1.0, amin=1e-10, top_db=80.0) -> Tensor:
    """定义torch版的功率谱转到分贝"""
    log_spec = 10.0 * torch.log10(torch.clip(power_spec, min=amin))
    log_spec -= 10.0 * torch.log10(torch.tensor(np.maximum(amin, ref), device=power_spec.device))
    log_spec = torch.maximum(log_spec, log_spec.max() - top_db)
    return log_spec


magnitude_pt = torch.abs(d_pt)
power_pt = torch.pow(magnitude_pt, 2)
db_pt = power_to_db(power_pt, 1.0)

magnitude_librosa = np.abs(d_librosa)
power_librosa = np.power(magnitude_librosa, 2)
db_librosa = librosa.power_to_db(power_librosa, ref=1)

print('torch版能量谱|分贝: ', power_pt.shape, '|', db_pt.shape)
print('librosa版能量谱|分贝: ', power_librosa.shape, '|', db_librosa.shape)

torch版能量谱|分贝:  torch.Size([1, 1025, 103]) | torch.Size([1, 1025, 103])
librosa版能量谱|分贝:  (1025, 103) | (1025, 103)

def plot_db_spectrum(spectrum_db, title=None, y_label="freq_bin"):
    """绘制db谱图"""
    fig, axs = plt.subplots(1, 1)
    axs.set_title(title or "Spectrum (db)")
    axs.set_ylabel(y_label)
    axs.set_xlabel("frame")
    im = axs.imshow(spectrum_db, origin="lower", aspect="auto")
    fig.colorbar(im, ax=axs, format='%+2.0f dB')
    plt.show(block=False)


print('能量谱: ')
plot_db_spectrum(db_pt[0], title='spectrum_db_pt')
plot_db_spectrum(db_librosa, title='spectrum_db_librosa')
print('librosa自带画图：')
librosa.display.specshow(db_librosa, sr=sample_rate_librosa,
                         y_axis='linear', x_axis='time')

能量谱:

librosa自带画图：

4 梅尔滤波器组 (Mel Filter Banks)

4.1 Mel频率

• 定义：1Mel为1kHz 音调感知程度的 1/1000。
• 目的：模拟人耳对不同频率语音的感知。人类对不同频率语音有不同的感知能力：1kHz 以下，与频率成线性关系；1kHz 以上，与频率成对数关系。
• 公式：$$\mathrm{Mel}(f)=2595\lg (1+\frac{f}{700})$$

4.2 Mel Filter Banks

• 定义：定义一个有$M$个三角滤波器的滤波器组（滤波器的个数和临界带的个数相近），$M$通常取22-40，26是标准，示例图取$n_{mels}=40$。滤波器组中的每个滤波器都是三角形的，中心频率为$f(m)$ ，中心频率处的响应为1，并向0线性减小，直到达到两个相邻滤波器的中心频率，其中响应为0。第1个滤波器组将从第1个点开始，在第2个点达到峰值，然后在第3个点返回零；第2个滤波器组将从第2个点开始，在第3个点达到最大值，然后在第4个点为零，以此类推。各$f(m)$之间的间隔随着$m$值的增大而增宽，如图所示：
  

• 公式：$$H_m(k)=\{\begin{array}{cl}0& k<f(m-1)\\ \frac{k-f(m-1)}{f(m)-f(m-1)}& f(m-1)\le k\le f(m)\\ \frac{f(m+1)-k}{f(m+1)-f(m)}& f(m)\le k\le f(m+1)\\ 0& k>f(m+1)\end{array}$$
  $\circ$其中，$m$是滤波器组数；
  $\circ$$f(\cdot )$是$m+2$个经过$Mel$尺度变换的$Mel$频率;

• 代码参数(torchaudio.functional.melscale_fbanks & librosa.filters.mel)：
  $\circ$n_fft(librosa)：FFT窗口长度；
  $\circ$n_freqs(pt)：要滤波的频率数量；
  $\circ$n_mels：Mel滤波器组个数；
  $\circ$f_min：最小截止频率(Hz)；
  $\circ$f_max：最大截止频率(Hz)；
  $\circ$norm：如果是slaney，将三角形的mel权重除以mel带的宽度(面积归一化)。默认为None。

4.3 FilterBanks特征

• 公式：$$FilterBanks(k)={\mathrm{X}}_{\text{power }}(k)\cdot H_m(k)$$
• 直接计算代码参数(torchaudio.transforms.MelSpectrogram&librosa.feature.melspectrogram)：
  $\circ$sample_rate, n_fftn_mels, hop_length, win_length, window等同上。

n_fft = 2048
n_mels = 40
sample_rate = 16000

mel_filters_pt = torchaudio.functional.melscale_fbanks(
    int(n_fft // 2 + 1),
    n_mels=n_mels,
    f_min=0.0,
    f_max=sample_rate / 2.0,
    sample_rate=sample_rate,
    norm="slaney",
)

mel_filters_librosa = librosa.filters.mel(
    sr=sample_rate,
    n_fft=n_fft,
    n_mels=n_mels,
    fmin=0.0,
    fmax=sample_rate / 2.0,
    norm="slaney",
    htk=True,
).T

print('torch版滤波器组: ', mel_filters_pt.shape)
print('librosa版滤波器组: ', mel_filters_librosa.shape)

torch版滤波器组:  torch.Size([1025, 40])
librosa版滤波器组:  (1025, 40)

print('Mel滤波器组：')
plot_db_spectrum(power_to_db(mel_filters_pt), title='spectrum_db_pt')
plot_db_spectrum(librosa.power_to_db(mel_filters_librosa), title='spectrum_db_librosa')

Mel滤波器组：

filter_banks_pt = torch.mm(power_pt[0].T, mel_filters_pt).T
filter_banks_librosa = np.matmul(power_librosa.T, mel_filters_librosa).T
print('torch版FilterBanks特征: ', filter_banks_pt.shape)
print('librosa版FilterBanks特征: ', filter_banks_librosa.shape)

torch版FilterBanks特征:  torch.Size([40, 103])
librosa版FilterBanks特征:  (40, 103)

print('FilterBanks特征: ')
plot_db_spectrum(power_to_db(filter_banks_pt), title='spectrum_db_pt')
plot_db_spectrum(librosa.power_to_db(filter_banks_librosa), title='spectrum_db_librosa')

FilterBanks特征:

MelSpectrogram = torchaudio.transforms.MelSpectrogram(
    sample_rate=sample_rate,
    n_fft=2048,
    n_mels=40,
    hop_length=512,
    win_length=None,
    window_fn=torch.hann_window
)
mel_spec_pt = MelSpectrogram(wave_pt)

mel_spec_librosa = librosa.feature.melspectrogram(
    y=wave_librosa,
    sr=sample_rate,
    n_fft=2048,
    n_mels=40,
    hop_length=512,
    win_length=None,
    window="hann",
)
print('torch版直接API计算的FilterBanks特征: ', mel_spec_pt.shape)
print('librosa版直接API计算的FilterBanks特征: ', mel_spec_librosa.shape)

torch版直接API计算的FilterBanks特征:  torch.Size([1, 40, 107])
librosa版直接API计算的FilterBanks特征:  (40, 107)

print('直接API计算FilterBanks特征: ')
plot_db_spectrum(power_to_db(mel_spec_pt[0]), title='spectrum_db_pt')
plot_db_spectrum(librosa.power_to_db(mel_spec_librosa), title='spectrum_db_librosa')

直接API计算FilterBanks特征:

5. 梅尔频率倒谱系数(Mel-frequency Cepstral Coefficients, MFCC)

• 说明：FilterBanks系数之间是高度相关的，对某些任务会出现问题，利用离散余弦变换(Discrete Cosine Transform, DCT)可以将FilterBanks系数去相关并且可以压缩信息。对于自动语音识别(ASR)任务，一般保留$2\sim 13$个系数，其余丢掉的原因是由于FilterBanks系数的快速变化，表征了声音细微的细节，对ASR任务意义不大。
• DCT公式：$$C(n)=\sum _{m=0}^{N-1}s(m)\cos (\frac{\pi n(m-0.5)}{M}),n=1,2,...,L$$
  $\circ$$L$阶指MFCC系数阶数，通常取$2\sim 13$；
  $\circ$$M$是Mel滤波器个数。
• 代码参数(torchaudio.functional.create_dct&librosa使用的是scipy下的dct)：
  $\circ$n_mfcc：MFCC系数阶数。
• 直接计算MFCC代码参数(torchaudio.transforms.MFCC&librosa.feature.mfcc)：
  $\circ$n_mfcc：MFCC系数阶数；其他参数同FilterBanks参数。

n_fft = 2048
win_length = None
hop_length = 512
n_mels = 40
n_mfcc = 12

dct_pt = torchaudio.functional.create_dct(n_mfcc=n_mfcc, n_mels=n_mels, norm=None)
mfcc_pt = torch.mm(power_to_db(filter_banks_pt.T), dct_pt).T

mfcc_librosa = dct(x=librosa.power_to_db(filter_banks_librosa), axis=-2)[..., :n_mfcc, :]

print('torch版的MFCC特征: ', mfcc_pt.shape)
print('librosa版的MFCC特征: ', mfcc_librosa.shape)

torch版的MFCC特征:  torch.Size([12, 103])
librosa版的MFCC特征:  (12, 103)

print('MFCC特征：')
plot_db_spectrum(power_to_db(mfcc_pt), title='spectrum_db_pt')
plot_db_spectrum(librosa.power_to_db(mfcc_librosa), title='spectrum_db_librosa')

MFCC特征：

mfcc_transform = torchaudio.transforms.MFCC(
    sample_rate=sample_rate,
    n_mfcc=n_mfcc,
    melkwargs={
        "n_fft": n_fft,
        "n_mels": n_mels,
        "hop_length": hop_length,
        "mel_scale": "htk",
    },
)
mfcc_api_pt = mfcc_transform(wave_pt)

mfcc_api_librosa = librosa.feature.mfcc(
    y=wave_librosa,
    sr=sample_rate,
    n_mfcc=n_mfcc,
    dct_type=2,
    norm="ortho",
    # --- mel kwargs --- #
    n_fft=n_fft,
    n_mels=n_mels,
    hop_length=hop_length,
    htk=True,
)

print('torch版直接API计算的MFCC特征: ', mfcc_api_pt.shape)
print('librosa版直接API计算的MFCC特征: ', mfcc_api_librosa.shape)

torch版直接API计算的MFCC特征:  torch.Size([1, 12, 107])
librosa版直接API计算的MFCC特征:  (12, 107)

print('直接API计算MFCC特征: ')
plot_db_spectrum(power_to_db(mfcc_api_pt[0]), title='mfcc_db_pt')
plot_db_spectrum(librosa.power_to_db(mfcc_api_librosa), title='mfcc_db_librosa')

直接API计算MFCC特征:

6. 参考

1. Speech Processing for Machine Learning: Filter banks, Mel-Frequency Cepstral Coefficients (MFCCs) and What’s In-Between
2. Mel Frequency Cepstral Coefficient (MFCC) tutorial
   
3. Torchaudio Documentation
4. librosa api
5. librosa语音信号处理