当前位置：首页 > news >正文

STFT和梅尔频谱图

news 2025/8/31 19:35:49

一、STFT（短时傅里叶变换）的详细原理与公式

1.1 傅里叶变换的局限性与STFT的诞生

傅里叶变换（FT）能将时域信号 $x (t)$ 转换为频域表示 $X(ω)X(\omega)$ ，公式为：
$X(\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} dt \tag{1}$
其中 $ω=2πf\omega = 2\pi f$ 为角频率， $f$ 为频率（Hz）。

但FT的缺陷是完全丢失了时间信息——只能得到信号中“有哪些频率”，无法知道“这些频率在何时出现”。例如，一段包含钢琴声和鼓声的音频，FT只能告诉你总频谱中有哪些频率成分，却无法区分钢琴声和鼓声的时间先后。

为解决这一问题，STFT通过**“加窗分段”**的方式，在“时间局部性”和“频率分辨率”之间取得平衡：将长信号切割为多个短时片段（假设每个片段是“平稳的”），对每个片段单独做傅里叶变换，最终拼接出“时间-频率”平面的能量分布。

1.2 STFT的数学定义

连续时间下，STFT定义为：
$\text{STFT}(\tau, \omega) = \int_{-\infty}^{\infty} x(t) g^*(t - \tau) e^{-j\omega t} dt \tag{2}$
其中：

$τ\tau$ 为“时间中心”（窗函数的滑动位置，对应横轴时间）；
$g (t)$ 为窗函数（如汉宁窗、矩形窗）， $g^*(t)$ 为其共轭；
物理意义：用窗函数 $\tau)$ 截取 $\tau$ 附近的信号片段，再对该片段做傅里叶变换，得到该时刻的频率谱。

1.3 离散化实现（实际工程中）

实际处理数字音频时，信号是离散的（采样率为 $f_s$ ，单位Hz），需对STFT离散化：

分帧：将离散信号 $x [n]$ （ $n$ 为采样点索引）按窗长 $N$ 分成多个帧，第 $m$ 帧为 $xm[n]=x[n+m⋅S]⋅w[n]x_m[n] = x[n + m \cdot S] \cdot w[n]$ ，其中：
- $S$ 为帧移（相邻帧的采样点间隔，通常取 $S = N /2$ 以减少帧间信息丢失）；
- $w [n]$ 为离散窗函数（如汉宁窗： $\cos(2\pi n/N))$ ， $\leq n < N$ ）。
傅里叶变换：对每帧做离散傅里叶变换（DFT）：
$X_m[k] = \sum_{n=0}^{N-1} x_m[n] e^{-j2\pi kn/N} \tag{3}$
其中 $k$ 为频率索引，对应频率 $fk=k⋅fs/Nf_k = k \cdot f_s / N$ （Hz）。
STFT矩阵：最终得到的离散STFT是一个矩阵，行索引为帧号 $m$ （对应时间），列索引为频率索引 $k$ （对应频率），元素为 $X_m[k]$ 。

1.4 频谱图（Spectrogram）的定义

频谱图是STFT幅度的平方（能量），公式为：
$\text{Spectrogram}(m, k) = |X_m[k]|^2 \tag{4}$

横轴：时间（由帧号 $m$ 转换， $tm=m⋅S/fst_m = m \cdot S / f_s$ ）；
纵轴：频率（ $fk=k⋅fs/Nf_k = k \cdot f_s / N$ ）；
值：对应频率的能量（颜色越深，能量越高）。

1.5 关键特性：时间与频率分辨率的权衡

窗长 $N$ 决定了分辨率：

窗长越长（ $N$ 越大）：频率分辨率越高（ $Δf=fs/N\Delta f = f_s / N$ 越小），但时间分辨率越低（ $Δt=N/fs\Delta t = N / f_s$ 越大）；
窗长越短（ $N$ 越小）：时间分辨率越高，但频率分辨率越低。

这一权衡源于不确定性原理： $Δt⋅Δf≥1/4π\Delta t \cdot \Delta f \geq 1/4\pi$ ，即无法同时无限提高时间和频率分辨率。在音频中，窗长通常取20-50ms（如16kHz采样率下，20ms对应320个采样点）。

二、梅尔频谱图（Mel Spectrogram）的详细原理与公式

梅尔频谱图是在STFT基础上，结合人耳听觉特性（梅尔刻度，Mel Scale）优化的时频表示，核心是将“线性频率”转换为“符合人耳感知的非线性频率”。

2.1 人耳听觉特性与梅尔刻度

人耳对频率的感知是非线性的：

对低频（如20-1000Hz）的分辨率高（能区分1Hz的差异）；
对高频（如1000-20000Hz）的分辨率低（可能需要几十Hz才能区分）。

梅尔刻度（Mel Scale）将线性频率（Hz）映射为“感知频率”（Mel），使“梅尔频率间隔相等”对应“人耳感知的音高间隔相等”。映射公式为：
$\text{Mel}(f) = 2595 \cdot \log_{10}\left(1 + \frac{f}{700}\right) \tag{5}$
逆映射（从Mel到Hz）为：
$f(\text{Mel}) = 700 \cdot \left(10^{\text{Mel}/2595} - 1\right) \tag{6}$

例如：

0Hz对应0Mel；
1000Hz对应 $2595 \cdot \log_{10}(1 + 1000/700) \approx 1000$ Mel（设计上使1000Hz≈1000Mel）；
2000Hz对应≈1700Mel（而非2000Mel，体现高频压缩）。

2.2 梅尔滤波器组（Mel Filter Banks）

为将STFT的线性频谱转换为梅尔频谱，需设计一组三角滤波器（梅尔滤波器组），每个滤波器覆盖一个梅尔频率区间，作用是“聚合”线性频率的能量到梅尔频率轴上。

步骤如下：

确定梅尔频率范围：设音频最高频率为 $fmaxf_{\text{max}}$ （如采样率的一半， $f_s/2$ ），则对应的梅尔频率范围为 $\text{Mel}(f_{\text{max}})]$ 。
生成梅尔中心频率：在梅尔频率范围内均匀取 $M + 2$ 个点（ $M$ 为滤波器数量，通常取20-128）：
$\text{mel}_0 = 0, \quad \text{mel}_1, \quad \dots, \quad \text{mel}_{M+1} = \text{Mel}(f_{\text{max}}) \tag{7}$
再通过逆映射公式（6）转换回线性频率，得到 $f0,f1,…,fM+1f_0, f_1, \dots, f_{M+1}$ （作为滤波器的边界频率）。
定义三角滤波器：第 $m$ 个滤波器 $H_m(k)$ （ $\leq m \leq M$ ）在频率索引 $k$ 处的响应为：
$H_m(k) = \begin{cases} 0 & \text{if } f_k < f_{m-1} \text{ or } f_k > f_{m+1} \\ \frac{f_k - f_{m-1}}{f_m - f_{m-1}} & \text{if } f_{m-1} \leq f_k \leq f_m \\ \frac{f_{m+1} - f_k}{f_{m+1} - f_m} & \text{if } f_m \leq f_k \leq f_{m+1} \end{cases} \tag{8}$
其中 $f_m$ 是第 $m$ 个滤波器的中心频率（线性频率）， $f_k$ 是STFT的频率轴（线性频率）。

2.3 梅尔频谱的计算

计算STFT功率谱：先通过STFT得到功率谱 $P(m, k) = |X_m[k]|^2$ （同公式4）。
梅尔滤波：将功率谱与每个梅尔滤波器卷积（点积），得到每个梅尔频段的能量：
$\sum_{k=0}^{N/2} P(t, k) \cdot H_m(k) \tag{9}$
其中 $t$ 为时间帧， $m$ 为梅尔滤波器索引（1≤m≤M）。
对数压缩：人耳对响度的感知接近对数关系，且对数能压缩能量范围（减少动态范围），因此通常取对数：
$\text{LogMel}(m, t) = \log\left(S(m, t) + \epsilon\right) \tag{10}$
其中 $ϵ\epsilon$ 是极小值（如1e-10），避免 $S (m, t) = 0$ 时的对数无穷大。实际中也常用 $10 \cdot \log_{10}(S(m, t) + \epsilon)$ （单位dB）。

2.4 梅尔频谱图的特性

横轴：时间（与STFT一致）；
纵轴：梅尔频率（非线性，由公式5映射）；
值：对数能量（符合人耳响度感知）。

相比STFT频谱图，梅尔频谱图更紧凑（高频被压缩）、更符合听觉特性，因此在语音识别、音频分类等任务中性能更优。

三、STFT与梅尔频谱图的核心区别（公式视角）

指标	STFT频谱图	梅尔频谱图
频率轴定义	线性频率 $fk=k⋅fs/Nf_k = k \cdot f_s / N$	梅尔频率 $Mel(fk)\text{Mel}(f_k)$ （公式5）
能量计算	$X_m[k])^2$ （公式4）	$log⁡(∑P(t,k)⋅Hm(k))\log\left(\sum P(t,k) \cdot H_m(k)\right)$ （公式9+10）
设计依据	数学时频分析（无感知优化）	人耳听觉特性（梅尔刻度+对数响度）
维度压缩	无（保留全部线性频率）	有（压缩为M个梅尔频段，M≪N/2）

总结：STFT是“忠实于物理频率”的时频分析工具，而梅尔频谱图是“忠实于人耳感知”的优化版本，两者的核心差异体现在频率轴的定义和能量聚合方式上。

使用python进行提取

在Python中，提取STFT和梅尔频谱图的常用库包括 Librosa（音频处理专用）、SciPy（科学计算）、Torchaudio（PyTorch生态）等。以下是各库的实现代码及说明：

1. Librosa（推荐，专为音频处理设计）

Librosa是最常用的音频特征提取库，封装了STFT和梅尔频谱图的便捷接口，且默认参数符合音频处理最佳实践。

安装

pip install librosa matplotlib  # matplotlib用于可视化

代码示例

import librosa
import librosa.display
import matplotlib.pyplot as plt
import numpy as np# 1. 加载音频文件（支持wav、mp3等格式）
# 这里使用Librosa自带的示例音频，也可替换为本地文件路径
audio_path = librosa.example('libri1')  # 示例音频：一段语音
y, sr = librosa.load(audio_path, sr=None)  # y: 音频信号数组，sr: 采样率# 2. 计算STFT并生成频谱图
n_fft = 2048  # FFT窗口大小（决定频率分辨率）
hop_length = 512  # 帧移（决定时间分辨率，通常为n_fft的1/4）
stft = librosa.stft(y, n_fft=n_fft, hop_length=hop_length)  # 复数矩阵：(1 + n_fft/2, 帧数)
stft_mag = np.abs(stft)  # 取幅度（能量）
stft_db = librosa.amplitude_to_db(stft_mag, ref=np.max)  # 转换为分贝值（增强动态范围）# 3. 计算梅尔频谱图
n_mels = 128  # 梅尔滤波器数量（梅尔频段数）
mel_spectrogram = librosa.feature.melspectrogram(y=y, sr=sr, n_fft=n_fft, hop_length=hop_length, n_mels=n_mels
)  # 梅尔频谱：(n_mels, 帧数)
mel_db = librosa.amplitude_to_db(mel_spectrogram, ref=np.max)  # 转换为分贝值# 4. 可视化
plt.figure(figsize=(15, 10))# 绘制STFT频谱图
plt.subplot(2, 1, 1)
librosa.display.specshow(stft_db,sr=sr,hop_length=hop_length,x_axis='time',  # 横轴：时间y_axis='hz'     # 纵轴：线性频率（Hz）
)
plt.colorbar(format='%+2.0f dB')
plt.title('STFT 频谱图')# 绘制梅尔频谱图
plt.subplot(2, 1, 2)
librosa.display.specshow(mel_db,sr=sr,hop_length=hop_length,x_axis='time',  # 横轴：时间y_axis='mel'    # 纵轴：梅尔频率
)
plt.colorbar(format='%+2.0f dB')
plt.title('梅尔频谱图')plt.tight_layout()
plt.show()

关键参数说明

n_fft：FFT窗口大小（如2048），决定频率分辨率（值越大，频率分辨率越高）。
hop_length：帧移（如512），决定时间分辨率（值越小，时间分辨率越高）。
n_mels：梅尔滤波器数量（如128），决定梅尔频谱的频率维度。
amplitude_to_db：将线性能量转换为分贝（dB），符合人耳对响度的对数感知。

2. SciPy（科学计算库，适合底层实现）

SciPy的signal模块提供STFT的底层实现，灵活性高，但需手动处理梅尔频谱图的转换（通常结合Librosa的梅尔滤波器组）。

安装

pip install scipy librosa matplotlib

代码示例

import numpy as np
import scipy.signal as signal
import matplotlib.pyplot as plt
import librosa  # 用于加载音频和梅尔滤波器组# 1. 加载音频
y, sr = librosa.load(librosa.example('libri1'), sr=None)# 2. 计算STFT（SciPy实现）
n_fft = 2048
hop_length = 512
f, t, stft = signal.stft(y,fs=sr,nperseg=n_fft,  # 窗口大小（同n_fft）noverlap=n_fft - hop_length,  # 重叠部分（n_fft - hop_length）window='hann'  # 汉宁窗（常用音频窗函数）
)
stft_mag = np.abs(stft)  # 幅度谱
stft_db = 20 * np.log10(stft_mag + 1e-10)  # 转换为分贝（手动计算）# 3. 可视化STFT
plt.figure(figsize=(10, 5))
plt.pcolormesh(t, f, stft_db, shading='gouraud')
plt.colorbar(label='dB')
plt.xlabel('时间 (s)')
plt.ylabel('频率 (Hz)')
plt.title('SciPy STFT 频谱图')
plt.tight_layout()
plt.show()
# 用SciPy+Librosa结合计算梅尔频谱图（SciPy无内置梅尔转换）# 1. 用SciPy计算STFT幅度谱
stft_mag = np.abs(stft)  # 来自上面的STFT结果# 2. 用Librosa生成梅尔滤波器组
mel_filters = librosa.filters.mel(sr=sr,n_fft=n_fft,n_mels=128  # 梅尔频段数
)  # 滤波器组形状：(n_mels, 1 + n_fft/2)# 3. 应用梅尔滤波（将线性频谱转换为梅尔频谱）
mel_spectrogram = np.dot(mel_filters, stft_mag)  # 矩阵乘法：(n_mels, 帧数)
mel_db = 20 * np.log10(mel_spectrogram + 1e-10)  # 转换为分贝# 4. 可视化梅尔频谱图
plt.figure(figsize=(10, 5))
librosa.display.specshow(mel_db,sr=sr,hop_length=hop_length,x_axis='time',y_axis='mel'
)
plt.colorbar(label='dB')
plt.title('SciPy+Librosa 梅尔频谱图')
plt.tight_layout()
plt.show()

3. Torchaudio（PyTorch生态，适合深度学习）

Torchaudio是PyTorch的音频处理库，支持GPU加速，输出为张量（可直接输入神经网络）。

安装

pip install torchaudio torch matplotlib

代码示例

import torch
import torchaudio
import torchaudio.transforms as T
import matplotlib.pyplot as plt# 1. 加载音频（返回张量格式）
waveform, sr = torchaudio.load(torchaudio.utils.download_asset("tutorial-assets/Lab41-SRI-VOiCES-src-sp0307-ch127535-sg0042.wav"))
# waveform形状：(1, 采样点数)，1表示单通道# 2. 计算STFT
n_fft = 2048
hop_length = 512
stft_transform = T.Spectrogram(n_fft=n_fft,hop_length=hop_length,power=2.0  # 1.0为幅度谱，2.0为功率谱
)
stft = stft_transform(waveform)  # 形状：(1, 1 + n_fft/2, 帧数)
stft_db = T.AmplitudeToDB()(stft)  # 转换为分贝# 3. 计算梅尔频谱图
mel_transform = T.MelSpectrogram(sample_rate=sr,n_fft=n_fft,hop_length=hop_length,n_mels=128
)
mel_spectrogram = mel_transform(waveform)  # 形状：(1, n_mels, 帧数)
mel_db = T.AmplitudeToDB()(mel_spectrogram)# 4. 可视化（需转换为numpy格式）
plt.figure(figsize=(15, 10))# STFT频谱图
plt.subplot(2, 1, 1)
plt.imshow(stft_db[0].numpy(),  # 取第0个通道origin='lower',aspect='auto',extent=[0, stft_db.shape[2] * hop_length / sr, 0, sr/2]
)
plt.colorbar(label='dB')
plt.xlabel('时间 (s)')
plt.ylabel('频率 (Hz)')
plt.title('Torchaudio STFT 频谱图')# 梅尔频谱图
plt.subplot(2, 1, 2)
plt.imshow(mel_db[0].numpy(),origin='lower',aspect='auto',extent=[0, mel_db.shape[2] * hop_length / sr, 0, 1]
)
plt.colorbar(label='dB')
plt.xlabel('时间 (s)')
plt.ylabel('梅尔频率')
plt.title('Torchaudio 梅尔频谱图')plt.tight_layout()
plt.show()