梅尔频率倒谱系数MFCC

梅尔倒谱系数（Mel-scale FrequencyCepstral Coefficients，简称MFCC）：依据人的听觉实验结果来分析语音的频谱
MFCC分析依据的听觉机理：

1. 第一梅尔刻度（Mel scale）：人耳感知的声音频率和声音的实际频率并不是线性的$\Rightarrow$从频率转换为梅尔刻度的公式为：$f_{mel}=2595\ast {\log }_{10}(1+\frac{f}{700})$；从梅尔回到频率：$f=700(10^{f_{mel}/2595}-1)$
   $f_{mel}$是以梅尔（Mel）为单位的感知频域（梅尔频域），$f$是以Hz为单位的实际语音频率
   
2. 第二临界带（Critical Band）：把进入人耳的声音频率用临界带进行划分，将语音在频域上划分成一系列的频率群，组成了滤波器组，即Mel滤波器组

人耳对不同频率的声波有不同的听觉敏感度$\Rightarrow$两个不同响度的声音作用于人耳时，响度较高的声音会影响到对响度较低的声音的感受
耳蜗的工作机制：耳蜗的基底膜对不同频率的声音振动有不同的响应

• 低频（低音）声音在基底膜上的行波传递距离较长，能在基底膜较远端引发振动
• 高频（高音）声音在基底膜上的传递距离较短，通常会在基底膜的靠近基部处引发振动

当低音和高音同时存在时，低频声波的传播会在更广泛的区域内产生振动，从而覆盖住高频部分，使得高频的声音被掩蔽$\Rightarrow$低频声音（如低音）具有较强的掩蔽能力，能掩盖掉频率相近的其他声音
临界带宽是指人耳在某个频率范围内能够区分的最小频率差，超过这个带宽外的声音可以被人耳区分，但在带宽内的声音容易被掩蔽掉

预处理

预处理包括预加重、分帧、加窗函数

import numpy
import scipy.io.wavfile
from scipy.fftpack import dctsample_rate, signal = scipy.io.wavfile.read('OSR_us_000_0010_8k.wav') 
signal = signal[0:int(3.5 * sample_rate)]  # 我们只取前3.5s

时域中的语音信号

预加重

预加重处理其实是将语音信号通过一个高通滤波器：$y(t)=x(t)-\alpha x(t-1)$，其中滤波器系数$\alpha$通常取值在 0.9 到 1 之间

emphasized_signal = numpy.append(signal[0], signal[1:] - pre_emphasis * signal[:-1])

预加重的作用：

1. 增强高频成分：语音信号的频谱特性通常是低频成分的能量较高，高频成分的能量较低。预加重通过增加高频成分的能量，使整个频谱更加平坦
2. 提高信噪比：高频成分往往更容易受到噪声的影响。预加重增强了这些高频分量，从而相对提高了它们的信噪比
3. 减少信号失真：低频成分可能会引起过度的信号能量集中，导致失真。预加重可以平衡频谱，减少这种失真

分帧

在大多数情况下，语音信号是非平稳的，对整个信号进行傅里叶变换是没有意义的
将信号帧化为20-40 ms帧，标准是25毫秒$\text{frame\_size}=0.025$，这意味着8kHz信号的帧长度为$0.025\ast 8000=200$个采样点。帧移通常为10毫秒$\text{frame\_stride}=0.01$（80个采样点）第一个语音帧从0开始，下一个语音帧从80开始，直到到达语音文件的末尾

• 帧长（frame length）：一帧的持续时间，如果帧长为 25 ms，采样率为 16,000 Hz，那么每一帧的采样点数为 400 个
• 帧移（frame step）：两帧之间的时间间隔，如果帧移为 10 ms，表示每一帧开始后过 10 ms 就会开始下一帧
• 重叠部分：帧长 25 ms，帧移 10 ms，那么每两帧之间有 15 ms 的重叠部分

frame_length, frame_step = frame_size * sample_rate, frame_stride * sample_rate  # 从秒转换为采样点
signal_length = len(emphasized_signal)
frame_length = int(round(frame_length))
frame_step = int(round(frame_step))
# 确保至少有1帧
num_frames = int(numpy.ceil(float(numpy.abs(signal_length - frame_length)) / frame_step))  pad_signal_length = num_frames * frame_step + frame_length
z = numpy.zeros((pad_signal_length - signal_length))
# 填充信号，确保所有帧的采样数相等，而不从原始信号中截断任何采样
pad_signal = numpy.append(emphasized_signal, z) indices = numpy.tile(numpy.arange(0, frame_length), (num_frames, 1)) + numpy.tile(numpy.arange(0, num_frames * frame_step, frame_step), (frame_length, 1)).T
frames = pad_signal[indices.astype(numpy.int32, copy=False)]

加窗

在信号处理（尤其是音频处理）过程中，将信号分割成帧后，每帧的开头和结尾容易产生不连续性或边界效应，在频域中会产生不必要的伪影或噪声（称为“频谱泄漏”）
通过对每个帧乘以窗函数，尤其是渐变的窗函数，可以减少这种不连续性，保证帧的两端更“平滑”，从而减小频谱泄漏的影响
Hamming 窗口是常用的一种窗函数：$w(n)=(1-\alpha )-\alpha \cdot \cos \left(\frac{2\pi n}{N-1}\right)$，其中$n$是采样点的索引、$N$是窗口/帧的长度、$\alpha =0.46$
Hamming 窗口的特点：它的两端值接近 0，中间值接近 1，因此乘以 Hamming 窗口会使帧的两端平滑下降，而中间部分的信号保留得较完整

• 减小边界效应：将每一帧乘以 Hamming 窗口后，信号在帧的左右两端逐渐减小到接近 0，这使得每帧在拼接时不再产生突然的跳跃，信号的开头和结尾更加“平滑”，从而减少了不连续的边界问题$\Rightarrow$将信号划分成多个帧时，每一帧都被认为是一个独立的信号片段。如果帧的开始和结束处信号值存在较大的突变，就会产生不连续性
• 降低频谱泄漏：当对分帧后的信号进行傅里叶变换时，如果帧的边界不连续，会导致在频域中泄漏到不相关的频率分量。Hamming 窗口通过将帧的两端平滑过渡，减少了不连续性，进而降低了频谱泄漏的影响$\Rightarrow$傅里叶变换假设信号是周期性的，因此它会自动认为每一帧的结束点与下一帧的起始点是连续的。如果实际情况中，帧的两端不连续，傅里叶变换会试图将这个突变表示为多种不同频率的叠加，从而在频域中引入了不相关的频率分量。这些不相关的频率分量扩散到整个频谱中，形成所谓的“频谱泄漏”

帧的拼接：将处理后的每帧按帧移的间隔进行逐一加和。因为相邻的帧之间有重叠部分，需要将这些重叠的部分叠加起来，而不是简单地直接连接

• 第一帧 frame_1（400 点）直接放在结果信号的起始位置（即 0 到 399 采样点）
• 第二帧 frame_2 从 160 采样点处开始，它的前 240 个点和 frame_1 的后 240 个点会发生重叠，因此它们重叠部分的值需要相加，而不是替换
• 依次类推，直到所有帧都拼接完成

frames *= numpy.hamming(frame_length)
# frames *= 0.54 - 0.46 * numpy.cos((2 * numpy.pi * n) / (frame_length - 1))  # 内部实现

FFT（Fourier-Transform）

由于信号在时域上的变换通常很难看出信号的特性，通常对它做FFT变换转换为频域上的能量分布来观察，不同的能量分布，代表不同语音的特性
短时傅里叶变换（STFT）：对每一帧信号执行一个长度为$N$的傅里叶变换，最终得到一个包含$N$个频率分量的频谱，其中$N$通常为256或512：$S_i(k)={\sum }_{n=1}^N{s}_i(n)e^{-j\frac{2\pi kn}{N}},1\le k\le K$

• $S_i(k)$：第$i$帧第$k$个频率分量
• $s_i(n)$：第$i$帧的第$n$个时域采样点的信号值，表示原始信号在时域的值
• $N$：：表示经过傅里叶变换后的频率分量的个数
• $k$：表示频率索引，取值范围为$1\le k\le K$，其中$K$通常是$N/2$，因为傅里叶变换的结果通常是对称的，前$K$项就包含所有独立的频率信息
• $e^{-j\frac{2\pi kn}{N}}$：傅里叶变换中的复数旋转因子
  $\frac{n}{N}$表示信号在一个周期内的相位变化，当$n$增加时，信号的相位也会变化，这种变化与频率$k$相关

功率谱

功率谱：对语音信号中的频谱取模平方（（取对数或者取平方）得到语音信号的谱线能量$P=\frac{|FFT(x_i){|}^2}{N}$，其中$x_i$是信号$X$的第$i$帧
在信号处理中，噪声往往集中在特定的频率范围内。通过查看信号的功率谱，可以识别出噪声频率并设计合适的滤波器去除这些噪声
信号在某一频率上的功率反映了该频率的强弱

滤波器组

计算Mel滤波器组：将功率谱通过一组Mel刻度(通常取40个滤波器，$nfilt=40$)的三角滤波器来提取频带(frequency bands)

• Mel刻度：人类对频率感知的非线性刻度$\Rightarrow$将实际的频率值转换为与人耳感知一致的数值，使得低频段的变化在Mel刻度上比高频段的变化更敏感$\Rightarrow$人耳对低频更敏感，能更好地区分相邻的频率，而对高频的分辨能力较差$\Rightarrow$$\mathrm{Mel}(f)=2595{\log }_{10}\left(1+\frac{f}{700}\right)$，其中$f$是频率，转换后的Mel值反映了人耳对不同频率的感知
• 临界带宽：在听觉系统中，能够分辨出的最小频率差$\Rightarrow$如果两个声音的频率差小于这个带宽，人耳就难以区分它们，它们会被感知为同一个声音
  依赖频率：临界带宽的大小与频率相关。在低频段，临界带宽较小，人耳可以分辨出较小的频率差异；而在高频段，临界带宽较大，意味着人耳难以分辨高频信号中的细微差异
  掩蔽效应：当低频声音和高频声音同时出现时，低频声音会“占据”较小的频率带宽，而高频声音的临界带宽较大，更容易被低频声音覆盖或掩蔽$\Rightarrow$在耳蜗中，低频声音的行波传播距离比高频声音更长，能够覆盖更广的区域，这意味着低频声音能更有效地激活耳蜗中较大的一部分，干扰甚至覆盖高频声音的感知；在声音的频谱中，低频声音的能量往往更大，且由于临界带宽较小，低频声音的能量集中，这种强度较高的低频声音更容易掩蔽掉能量相对较弱、临界带宽较大的高频声音$\iff$当低频和高频声音同时出现时，低频声音的强度和较小的临界带宽使得它更容易“占据”听觉频率空间，而高频声音由于较大的临界带宽和相对较弱的能量分布，更容易受到低频声音的干扰或掩蔽
• Mel滤波器组：就像人类的听觉感知系统（耳朵），人耳只关注某些特定的频率分量（人的听觉对频率是有选择性的）。它对不同频率信号的灵敏度是不同的，换言之，它只让某些频率的信号通过，而压根就直接无视它不想感知的某些频率信号。但是这些滤波器在频率坐标轴上却不是统一分布的，在低频区域有很多的滤波器，他们分布比较密集，但在高频区域，滤波器的数目就变得比较少，分布很稀疏。因此Mel刻度的目的是模拟人耳对声音的非线性感知，在较低的频率具有更强的辨别力
  滤波器：用于分离或增强某些频率的声音信号，常见的滤波器包括低通滤波器（只允许低频通过）、高通滤波器（只允许高频通过）以及带通滤波器（只允许一定范围的频率通过）
  Mel刻度：Mel刻度是基于人耳对频率感知的非线性频率刻度。它对低频敏感度更高，高频则较不敏感$\Rightarrow$Mel滤波器组（通常是三角滤波器）用于分割信号的频谱，每个滤波器的中心频率依据Mel刻度分布，低频滤波器更密集，高频滤波器较稀疏
  临界带宽：在滤波器设计中，临界带宽用来确定每个滤波器的带宽

语音频率和Mel频率之间转换：

• 从频率转换为梅尔刻度的公式为：$\mathrm{Mel}(f)=2595{\log }_{10}\left(1+\frac{f}{700}\right)$
• 从梅尔回到频率：$f=700(10^{f_{mel}/2595}-1)$

定义一个有M个三角滤波器的滤波器组（滤波器的个数和临界带的个数相近），$M$通常取22-40，26是标准。滤波器组中的每个滤波器都是三角形的，中心频率为$f(m)$，中心频率处的响应为1，并向0线性减小，直到达到两个相邻滤波器的中心频率，其中响应为0，各$f(m)$之间的间隔随着m值的增大而增宽$\Rightarrow$每个滤波器的响应形状是三角形，这意味着在中心频率$f(m)$处，滤波器的响应最强（即响应值为1），而随着频率远离中心，滤波器的响应会线性下降，直到在相邻滤波器的中心频率处降为0。这样的设计使得相邻滤波器在频谱上有所重叠，确保信号中的频率信息平滑过渡，不会突然丢失

• 中心频率$f(m)$：每个滤波器都有一个中心频率，它是该滤波器对输入信号最敏感的频率
• 线性减小：滤波器对中心频率的两侧频率的响应逐渐减弱。这种线性下降的设计让每个滤波器只对特定的频率区间有效
• 三角滤波器的频率响应定义为：$H_m(k)=\{\begin{array}{cc}0& k<f(m-1)\\ & \\ \frac{k-f(m-1)}{f(m)-f(m-1)}& f(m-1)\le k<f(m)\\ & \\ 1& k=f(m)\\ & \\ \frac{f(m+1)-k}{f(m+1)-f(m)}& f(m)<k\le f(m+1)\\ & \\ 0& k>f(m+1)\end{array}$

奈奎斯特频率是离散信号系统采样频率的一半

# --*-- coding:utf-8 --*--
# @Author : 一只楚楚猫
# @File : 01滤波器组.py
# @Software : PyCharmimport numpy as np# 示例参数
sample_rate = 16000  # 采样率 16kHz
NFFT = 512  # FFT点数
nfilt = 40  # 滤波器组的滤波器数量
low_freq_mel = 0  # 最低的Mel频率
high_freq_mel = (2595 * np.log10(1 + (sample_rate / 2) / 700))  # 最高Mel频率，对应Nyquist频率# 1. 生成Mel刻度上的点
mel_points = np.linspace(low_freq_mel, high_freq_mel, nfilt + 2)  # 在Mel刻度上均匀生成42个点
# hz_points = [[   0.           44.37407701   91.56109503  141.73937073  195.09852453,  251.84019719  312.17881177  376.34238398  444.57338374  517.12965156,  594.28537283  676.33211398  763.57992429  856.35850754  955.01846792, 1059.93263499 1171.49747253 1290.13457677]
# hz_points.shape: (nfilt + 2, )
hz_points = 700 * (10 ** (mel_points / 2595) - 1)  # 将Mel刻度转换回Hz刻度# 2. 计算每个频率点在FFT频率bin中的索引
# bins = [  0   1   2   4   6   8  10  12  14  16  19  21  24  27  30  33  37  41,  45  49  54  59  64  69  75  81  88  95 103 110 119 128 137 148 158 170, 182 195 209 224 239 256]
# bins.shape: (nfilt + 2, )
bins = np.floor((NFFT + 1) * hz_points / sample_rate).astype(int)  # 每个Hz点对应的频率bin# 3. 创建滤波器组的矩阵
fbank = np.zeros((nfilt, int(np.floor(NFFT / 2 + 1))))for m in range(1, nfilt + 1):f_m_minus = bins[m - 1]  # 左边界f_m = bins[m]  # 中间频率f_m_plus = bins[m + 1]  # 右边界# 左侧线性上升for k in range(f_m_minus, f_m):fbank[m - 1, k] = (k - f_m_minus) / (f_m - f_m_minus)# 右侧线性下降for k in range(f_m, f_m_plus):fbank[m - 1, k] = (f_m_plus - k) / (f_m_plus - f_m)# 4. 假设我们有信号的功率谱 pow_frames (形状: num_frames x NFFT/2+1)
# 假设 pow_frames 是任意生成的矩阵
num_frames = 100  # 假设有100帧信号
# pow_frames.shape: (num_frames, NFFT/2+1)
pow_frames = np.random.rand(num_frames, int(NFFT / 2 + 1))  # 随机生成的功率谱矩阵# 5. 将滤波器组应用到功率谱上
# fbank.shape: (nfilt, int(np.floor(NFFT / 2 + 1)))
filter_banks = np.dot(pow_frames, fbank.T)# 6. 防止数值稳定性问题，将0替换为非常小的数值
# filter_banks.shape: (num_frames, nfilt)
filter_banks = np.where(filter_banks == 0, np.finfo(float).eps, filter_banks)# 7. 转换为分贝(dB)
filter_banks = 20 * np.log10(filter_banks)# 打印滤波器组的响应
print(f"Mel滤波器组的响应形状: {filter_banks.shape}")

requency bins：频域中样本之间的间隔

• 例如，如果采样率为100 Hz，FFT大小为100，则在[0 100）Hz之间有100个点。因此，将整个100 Hz范围划分为100个间隔，如0-1 Hz、1-2 Hz等。每个这样的小间隔，例如0-1 Hz，都是一个frequency bin

梅尔频率倒谱系数（MFCC）

可以把 MFCC 想象成将声音“简化成一张声音的名片”，这张名片包含了声音的主要特征，比如音调、音色等，而忽略掉了一些不太重要的细节
$C(n)={\sum }_{m=0}^{N-1}\text{log}(x_{mel})\cos (\frac{\pi n(m-0.5)}{M}),n=1,2,\dots ,L$，其中$C_n$表示第$n$个梅尔频率倒谱系数、$N$是窗口/帧的长度、$M$是三角滤波器个数、$L$阶指MFCC系数阶数$\Rightarrow$MFCC 特征向量由不同阶的倒谱系数组成。每个系数捕捉了语音信号在频率空间中的不同信息，通常只保留低阶的 MFCC 系数（如前 12 或 13 阶），因为这些系数包含了主要的音频特征，高阶系数则多是噪声或细节，不利于识别任务


时间分辨率：时域图横坐标上最小的时间间隔，即模拟信号采样后得到的离散信号中的密集程度，提高采样率可以提高时域内的时间分辨率

均值归一化

为了平衡频谱并改善信噪比（SNR），可以简单地从所有帧中减去每个系数的平均值

filter_banks -= (numpy.mean(filter_banks, axis=0) + 1e-8)

对于MFCC：

mfcc -= (numpy.mean(mfcc, axis=0) + 1e-8)

在处理音频信号时，背景噪声常常表现为某些频率成分上的恒定噪声，这些噪声在整个音频中几乎保持一致。通过从所有帧中减去每个频率上的平均值，可以消除或减弱这种静态噪声，使得音频信号中的有用成分更加突出

总结

计算Mel刻度滤波器组Filter Banks由语音信号的性质和人类对这些信号的感知所驱动的；计算MFCC由一些机器学习算法的限制所驱动的
Filter Banks和MFCC对比：

• 计算量：MFCC是在FBank的基础上进行的
• 特征区分度：FBank特征相关性较高（相邻滤波器组有重叠），MFCC具有更好的判别度
• 信息量：Filter Banks包含比MFCC更多的信息

如果机器学习算法不易受高度相关输入的影响，则使用Mel刻度滤波器组。如果机器学习算法易受相关输入的影响，则使用MFCC

参考文献

1、梅尔频率倒谱系数(MFCC)的原理讲解及python实现
2、时间分辨率、频率分辨率