信号处理新手必看：如何用Python实现STFT替代FFT分析非平稳信号？

最新推荐文章于 2026-07-03 16:03:34 发布

原创最新推荐文章于 2026-07-03 16:03:34 发布 · 735 阅读

本内容遵循CC 4.0 BY-SA版权协议

信号处理新手必看：如何用Python实现STFT替代FFT分析非平稳信号？

如果你刚开始接触信号处理，很可能第一个学会的“大招”就是傅里叶变换（FFT）。它能将一个复杂的声音、一段脑电波，或者任何随时间变化的信号，分解成一个个不同频率的正弦波，让你一眼看清信号里“藏”着哪些频率成分。这感觉就像拿到了一把万能钥匙，很多问题迎刃而解。

但很快，你可能会遇到一个尴尬的局面：当你分析一段音乐、一段语音，或者机械振动信号时，FFT给出的频谱图看起来“一团糟”，或者根本无法解释信号在时间上的变化。比如，一段先有鸟鸣后有雷声的录音，FFT只能告诉你这段录音里包含了鸟鸣和雷声的所有频率，却无法告诉你鸟鸣在何时开始、雷声在何时响起。这就是FFT的“阿喀琉斯之踵”——它完美地揭示了频率，却彻底丢失了时间信息。

问题的根源在于信号的平稳性。FFT假设你分析的信号在整个时间段内，其统计特性（如频率成分）是不变的。然而，现实世界中绝大多数让我们感兴趣的信号，如语音、音乐、金融时间序列、生物医学信号（EEG/ECG），都是非平稳的，它们的频率内容会随着时间动态变化。这时，你需要一个能同时看清“何时”出现“何种”频率的工具。这就是短时傅里叶变换（STFT） 登场的时刻，它本质上是给FFT加了一个“滑动的窗口”，让我们得以一窥信号的时频全貌。本文将彻底抛弃MATLAB的视角，手把手带你用Python生态中的强大工具（librosa和scipy），从原理到实战，构建属于你自己的非平稳信号分析流程。

1. 核心理念：从全局频谱到时频局部化

理解STFT，关键在于理解“短时”和“加窗”这两个概念。我们放弃一次性分析整个漫长信号这种“上帝视角”，转而采取一种更务实的“局部观察”策略。

想象一下，你要分析一首长达五分钟的交响乐。FFT的做法是把整首曲子当作一个整体，一次性进行变换，结果得到的是这五分钟内所有乐器所有音符频率的总和，你无法区分小提琴独奏是在开头还是结尾。STFT则像是一个拿着放大镜沿着时间轴慢慢移动的乐评人。这个放大镜就是窗函数。他每次只观察一个很短的片段（比如0.1秒），对这个片段进行FFT，记录下这个瞬间有哪些频率。然后，他把放大镜稍微往前移动一点（可能会与之前观察的片段有重叠），再分析下一个瞬间。最终，他把所有瞬间的频谱结果按时间顺序排列起来，就得到了一张声谱图——一张以时间为横轴、频率为纵轴、颜色深浅表示能量强弱的图像。

这个过程中有几个关键参数决定了你观察的“清晰度”：

窗长：放大镜的宽度。窗越长，你对频率的分辨能力越强（能区分开两个非常接近的频率），但对时间变化的捕捉能力越弱（无法精确定位频率变化的瞬间）。
窗类型：放大镜的镜片材质。不同的窗函数（如汉明窗、汉宁窗、矩形窗）在时域和频域有不同的特性，主要影响频谱泄漏（能量扩散到邻近频率）的程度。
重叠长度：每次移动放大镜时，与上一次观察区域的重叠部分。更多的重叠能让时间轴上的变化更平滑，减少信息跳跃感，但也会增加计算量。

下表对比了FFT与STFT的核心差异，帮助你建立直观认识：

特性维度	傅里叶变换 (FFT)	短时傅里叶变换 (STFT)
核心输出	全局频谱（频率 vs. 幅度/相位）	时频谱图（时间 vs. 频率 vs. 幅度）
时间信息	完全丢失	保留，可定位频率成分的发生时刻
适用信号	平稳信号（频率成分不随时间变化）	非平稳信号（频率成分随时间变化）
分析视角	全局、整体	局部、滑动
关键参数	无（或仅FFT点数）	窗函数、窗长、重叠点数、FFT点数
结果可视化	频谱图	声谱图/时频图

提示：STFT并非完美无缺。它受限于“海森堡测不准原理”在时频分析中的体现——你无法同时获得无限高的时间分辨率和频率分辨率。窗长的选择永远是一种权衡，这正是更高级方法（如小波变换）试图解决的问题。但对于绝大多数入门到中级应用，STFT已经足够强大和直观。

2. Python实战：用Librosa和Scipy轻松实现STFT

理论聊完，我们进入最激动人心的实战环节。Python在信号处理领域拥有媲美甚至超越MATLAB的生态，这里我们重点介绍两个最常用、最友好的库：Librosa（专注于音频/音乐分析）和Scipy（通用的科学计算工具包）。

2.1 环境搭建与示例信号生成

首先，确保你的环境中安装了必要的库。打开终端或Anaconda Prompt，执行以下命令：

pip install numpy scipy matplotlib librosa

我们将生成一个经典的线性扫频信号作为分析对象，它频率随时间线性增加，是非平稳信号的典型代表。

import numpy as np
import matplotlib.pyplot as plt

# 生成线性扫频信号
fs = 1000  # 采样率 1000 Hz
duration = 2  # 信号时长 2 秒
t = np.linspace(0, duration, int(fs * duration), endpoint=False)  # 时间向量

# 频率从 0 Hz 线性增加到 250 Hz
f0 = 0
f1 = 250
chirp_signal = np.sin(2 * np.pi * (f0 + (f1 - f0) * t / (2 * duration)) * t)

# 为了增加趣味性，我们再叠加一个在1秒时突然出现的100Hz单频信号
pure_tone = np.zeros_like(t)
pure_tone[int(fs * 1):] = 0.5 * np.sin(2 * np.pi * 100 * t[int(fs * 1):])
composite_signal = chirp_signal + pure_tone

# 绘制时域波形
fig, ax = plt.subplots(2, 1, figsize=(10, 6))
ax[0].plot(t, chirp_signal)
ax[0].set_title('线性扫频信号 (0Hz -> 250Hz)')
ax[0].set_xlabel('时间 [秒]')
ax[0].set_ylabel('幅度')
ax[0].grid(True)

ax[1].plot(t, composite_signal)
ax[1].set_title('复合信号 (扫频 + 1秒后出现的100Hz纯音)')
ax[1].set_xlabel('时间 [秒]')
ax[1].set_ylabel('幅度')
ax[1].grid(True)
plt.tight_layout()
plt.show()

运行这段代码，你会看到两个时域波形图。仅从波形上，我们很难直观看出第二个复合信号中包含了哪些频率成分以及它们何时出现。接下来，就是STFT大显身手的时候。

2.2 使用Scipy.signal.spectrogram进行基础分析

Scipy的signal模块提供了直接的spectrogram函数，其接口设计非常接近MATLAB，易于上手。

from scipy import signal

# 使用 scipy.signal.spectrogram
nperseg = 256  # 窗长度，对应时间分辨率
noverlap = 200  # 重叠点数，重叠率约为 200/256 ≈ 78%
nfft = 512     # FFT点数，通常 >= nperseg，用于提升频率插值精度

f_scipy, t_scipy, Sxx_scipy = signal.spectrogram(composite_signal, fs,
                                                  window='hamming',
                                                  nperseg=nperseg,
                                                  noverlap=noverlap,
                                                  nfft=nfft,
                                                  mode='magnitude')

# 绘制声谱图
plt.figure(figsize=(10, 5))
# 使用 pcolormesh 绘制，坐标轴更准确
plt.pcolormesh(t_scipy, f_scipy, 10 * np.log10(Sxx_scipy + 1e-10), shading='gouraud')
# 注意：我们对幅度谱取对数并转换为分贝(dB)尺度，这是观察能量变化的常用方法，+1e-10是为了避免log10(0)
plt.colorbar(label='强度 [dB]')
plt.title('Scipy 生成的声谱图 (窗长=256, 重叠=200)')
plt.xlabel('时间 [秒]')
plt.ylabel('频率 [Hz]')
plt.ylim(0, 300)  # 限制频率显示范围，因为我们信号最高频率约250Hz
plt.tight_layout()
plt.show()

观察生成的声谱图，你应该能清晰地看到两条轨迹：一条从0Hz斜向上至250Hz的亮线（我们的扫频信号），另一条在1秒后出现的、位于100Hz的水平亮线（后加入的纯音）。STFT成功地将时间信息与频率信息关联了起来！

2.3 使用Librosa进行更高级的音频特性分析

如果你是处理音频信号（如.wav, .mp3文件），Librosa是更专业的选择。它封装了许多音频分析中常用的预处理步骤和更符合听觉特性的显示方式。

import librosa
import librosa.display

# 假设我们有一个音频文件 ‘example_audio.wav’，这里我们用之前生成的复合信号模拟
# librosa 通常期望采样率为 22050 Hz，这里我们重新采样以适配
composite_audio = composite_signal.astype(np.float32) # 确保数据类型

# 使用 librosa.stft
n_fft = 2048  # FFT窗口大小
hop_length = 512  # 帧移（窗长 - 重叠），这里重叠为 n_fft - hop_length
win_length = n_fft  # 窗长，通常等于n_fft

# 计算STFT矩阵 (复数)
D = librosa.stft(composite_audio, n_fft=n_fft, hop_length=hop_length, win_length=win_length, window='hann')

# 将复数谱转换为幅度谱
S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max)

# 使用librosa专属的显示函数绘制声谱图
plt.figure(figsize=(12, 6))
librosa.display.specshow(S_db, sr=fs, hop_length=hop_length, x_axis='time', y_axis='linear')
plt.colorbar(format='%+2.0f dB')
plt.title('Librosa 生成的声谱图 (更符合音频分析习惯)')
plt.xlabel('时间 [秒]')
plt.ylabel('频率 [Hz]')
plt.ylim(0, 300)
plt.tight_layout()
plt.show()

librosa.display.specshow 的显示效果通常更美观，并且其y_axis参数可以设置为‘log’来使用对数频率轴，这更接近人耳的听觉感知（Mel谱、CQT谱等是librosa的强项，常用于音乐信息检索和语音识别）。

3. 参数调优：窗函数、重叠与分辨率的艺术

生成一张声谱图很容易，但生成一张信息清晰、易于解读的声谱图，则需要理解参数背后的权衡。这是STFT分析中最具“手艺”的环节。

3.1 窗函数的选择：平衡主瓣与旁瓣

窗函数决定了你如何“截取”那一小段信号。不同的窗在时域乘上信号，在频域则相当于与信号的频谱进行卷积，这会导致频谱泄漏——某个单一频率的能量会“泄漏”到相邻的频率bin中。

矩形窗：相当于直接截断。主瓣最窄（频率分辨率潜在最佳），但旁瓣极高，泄漏严重，除非信号恰好是周期的整数倍，否则不推荐使用。
汉宁窗：旁瓣衰减很快，能有效抑制泄漏，但主瓣较宽（频率分辨率稍差）。是音频分析中的常用选择。
汉明窗：与汉宁窗类似，但第一个旁瓣更低，主瓣稍宽。是通用信号处理的默认选择之一。
布莱克曼窗：主瓣更宽，但旁瓣衰减极好，泄漏最小。

你可以通过一个简单的实验来感受差异：

windows = ['boxcar', 'hann', 'hamming', 'blackman']
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
axes = axes.ravel()

for idx, win in enumerate(windows):
    f, t, Sxx = signal.spectrogram(composite_signal, fs, window=win, nperseg=256, noverlap=200)
    ax = axes[idx]
    pcm = ax.pcolormesh(t, f, 10 * np.log10(Sxx + 1e-10), shading='gouraud')
    ax.set_title(f'窗函数: {win}')
    ax.set_xlabel('时间 [秒]')
    ax.set_ylabel('频率 [Hz]')
    ax.set_ylim(0, 300)
    fig.colorbar(pcm, ax=ax, format='%+2.0f dB')

plt.tight_layout()
plt.show()

观察四张图，你会发现使用boxcar（矩形窗）时，频率轨迹可能显得更“细”但背景噪声（泄漏）更明显；而blackman窗的轨迹可能略“粗”，但背景更干净。对于包含强单频成分的信号，汉明或汉宁窗通常是安全的起点。

3.2 窗长、重叠与时空分辨率权衡

这是STFT最核心的权衡，可以用一个简单的比喻：窗长是望远镜的焦距，重叠是拍摄视频的帧率。

长窗：高频率分辨率，低时间分辨率。适合分析变化缓慢的频率成分。
短窗：高时间分辨率，低频率分辨率。适合捕捉快速的瞬时事件（如敲击声、爆破音）。
高重叠：时间轴上的变化更平滑，减少“块状”效应，计算量增大。
低重叠：可能导致时间信息丢失，特别是对于短窗情况。

下面的代码演示了不同窗长对同一信号的影响：

window_lengths = [64, 128, 256, 512]  # 分别对应 64ms, 128ms, 256ms, 512ms (在fs=1000时)
fig, axes = plt.subplots(2, 2, figsize=(13, 9))
axes = axes.ravel()

for idx, nseg in enumerate(window_lengths):
    f, t, Sxx = signal.spectrogram(composite_signal, fs, window='hamming',
                                     nperseg=nseg, noverlap=int(nseg*0.75), nfft=max(512, nseg))
    ax = axes[idx]
    pcm = ax.pcolormesh(t, f, 10 * np.log10(Sxx + 1e-10), shading='gouraud')
    ax.set_title(f'窗长 = {nseg} 点 ({1000*nseg/fs:.0f} ms)')
    ax.set_xlabel('时间 [秒]')
    ax.set_ylabel('频率 [Hz]')
    ax.set_ylim(0, 300)
    fig.colorbar(pcm, ax=ax, format='%+2.0f dB')

plt.tight_layout()
plt.show()

仔细观察：

窗长64：100Hz的水平线（纯音）在时间上定位非常精准（在1秒处突然开始），但扫频信号的斜线变得很“粗”，频率分辨率差，你几乎看不清它从0到250Hz的连续变化。
窗长512：扫频信号的斜线变得非常清晰、锐利，频率分辨率极高。但是，100Hz纯音的开始时间变得模糊，它从1秒处“扩散”开，时间分辨率下降。

注意：在实际项目中，没有“最佳”参数，只有“最适合当前分析目标”的参数。如果你关心心跳信号中R波的位置（瞬时事件），就用短窗；如果你关心脑电信号中Alpha波（8-13Hz）的强度变化（相对缓慢），就用长窗。通常需要通过反复试验来确定。

4. 超越基础：实用技巧与常见陷阱

掌握了基本操作和参数调整后，一些进阶技巧和避坑指南能让你的分析工作更加稳健和专业。

4.1 从STFT结果中重构信号

STFT通常是可逆的（在满足一定条件下），这意味着你可以修改时频谱图，然后通过逆变换得到修改后的时域信号。这在音频处理（如降噪、时频滤波）中非常有用。

# 使用 librosa 进行重构示例
# 1. 计算STFT
D = librosa.stft(composite_audio, n_fft=2048, hop_length=512, win_length=2048, window='hann')

# 2. 假设我们做一个简单的操作：将100Hz以上的频率成分衰减（低通滤波）
# 获取频率轴
freqs = librosa.fft_frequencies(sr=fs, n_fft=2048)
# 找到100Hz对应的索引
cutoff_bin = np.argmax(freqs > 100)
# 创建滤波器掩码
mask = np.ones_like(D)
mask[cutoff_bin:, :] = 0.1  # 将100Hz以上的成分衰减到原来的10%

D_filtered = D * mask

# 3. 逆STFT重构时域信号
y_reconstructed = librosa.istft(D_filtered, hop_length=512, win_length=2048, window='hann')

# 对比原始信号与重构信号
plt.figure(figsize=(12, 4))
plt.plot(t[:len(y_reconstructed)], composite_audio[:len(y_reconstructed)], alpha=0.7, label='原始信号')
plt.plot(t[:len(y_reconstructed)], y_reconstructed, alpha=0.7, label='滤波后重构信号')
plt.title('STFT滤波与信号重构对比')
plt.xlabel('时间 [秒]')
plt.ylabel('幅度')
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

你会发现，重构后的信号中，高频的扫频成分大幅减弱，而1秒后的100Hz纯音基本保留。这演示了时频滤波的基本思路。

4.2 避免常见陷阱

混叠效应：确保你的信号采样率fs至少是信号最高频率的两倍（奈奎斯特采样定理）。STFT无法解决原始信号已混叠的问题。
边界效应：在信号的开始和结束处，窗函数可能覆盖到信号之外的区域（补零）。这会导致这些时间点的频谱估计不可靠。通常可以忽略边缘部分，或使用librosa的center=False参数。
分贝尺度的使用：直接绘制幅度谱abs(S)时，弱信号可能被强信号淹没。使用分贝尺度20*log10(abs(S))或librosa.amplitude_to_db可以极大地增强动态范围，让弱成分可见。但要注意，分贝值依赖于参考值ref。
计算效率：对于超长信号，直接计算整个STFT矩阵可能内存不足。可以考虑分块处理，或者使用scipy.signal.spectrogram的return_onesided=False参数并结合nperseg合理控制矩阵大小。

4.3 扩展到真实世界应用

STFT是许多高级应用的基石：

语音识别：梅尔频率倒谱系数就是从STFT谱图演变而来的。
音乐信息检索：和弦检测、节拍跟踪都依赖于时频分析。
故障诊断：通过分析机械振动信号的时频谱，可以检测轴承故障等瞬时冲击事件。
生物医学工程：分析EEG信号中特定频段（如Delta, Theta, Alpha波）的能量随时间的变化，用于睡眠分期或脑机接口。

例如，在分析一段包含敲击声的音频时，你可以设置一个非常短的窗（如10ms）和高重叠，来精准定位每一次敲击发生的时刻，即使它们的频率内容很宽。而在分析一段悠长的提琴音时，则使用长窗来精确测量其基频和泛音。

走到这里，你已经不再是那个只会用FFT看全局频谱的新手了。STFT赋予了你一双能够透视信号时频结构的“眼睛”。记住，所有的参数——窗函数、窗长、重叠——都是你调节这双眼睛“焦距”和“敏感度”的旋钮。真正的掌握，源于在具体项目中的反复尝试和思考：我到底想从信号中看到什么？是瞬态事件的发生时刻，还是缓慢振荡的精确频率？当你开始带着这个问题去调整nperseg和noverlap时，你就已经踏上了信号分析工程师的道路。我最初在分析野外鸟类录音时，为了分离重叠的鸟鸣声，花了整整一个下午调整窗长，最终发现对于那种短促的啁啾声，一个大约30ms的汉明窗配合70%的重叠效果最好。这种经验，是任何教程都无法直接给你的，它来自于你对自己数据和目标的深刻理解。

标签