4、信号处理基础（二）：包络谱分析、小波变换入门、时频分析（短时傅里叶变换STFT）、Python实战
各位同行，欢迎来到信号处理的第二讲。
上一章我们聊了傅里叶变换，说白了就是把时域信号扔到频域里看。但干我们这行的都知道，风机这种旋转机械，信号里藏着大量非平稳成分。光看频谱，很多时候是看不全的。
举个例子。你测到一个振动信号，频谱上有个峰值。你以为是某个轴承的故障频率。结果拆开一看，是齿轮的啮合频率被调制了。为什么？因为频谱把时间信息给丢了。
所以这一章，我们来解决这个问题。我会带你搞懂包络谱分析、小波变换和短时傅里叶变换。这些都是我这些年做故障诊断时，最常用的几把刀。
4.1 包络谱分析：从调制信号里挖出故障特征
先说说包络谱。这玩意儿在轴承故障诊断里，几乎是必用的。
你想想看，一个正常的轴承，滚动体滚过外圈，振动是平稳的。但一旦外圈出现点蚀，滚动体每滚过一次，就会产生一个冲击。这个冲击会激起系统的固有频率，形成高频振荡。而冲击本身是低频的，比如几十赫兹。
这就形成了一个调制信号：高频载波（固有频率）被低频信号（故障频率）调制。
直接对原始信号做FFT，你看到的可能是高频段的一大片能量。低频的故障特征被淹没了。怎么办？
包络谱的思路很简单：先把高频载波“剥掉”，提取出低频的包络线，再对这个包络线做FFT。
核心步骤：
- 带通滤波：把信号中感兴趣的共振频带滤出来。我一般会先看频谱，找到能量集中的高频区域。
- 包络提取：对滤波后的信号做Hilbert变换，得到解析信号，取模就是包络。
- FFT分析：对包络信号做傅里叶变换，得到包络谱。
我在项目里遇到过一回。某台引风机，振动值不大，但声音有点“发闷”。频谱上看不出明显异常。我习惯性地做了个包络谱，结果在故障频率的1倍频、2倍频上看到了清晰的峰值。拆开一看，外圈已经剥落了一块。嗯，这就是包络谱的威力。
避坑指南：
我曾经犯过一个错：没做带通滤波，直接对整个信号做Hilbert变换。结果包络线里混入了大量低频干扰，包络谱一团糟。记住，滤波的频带选择很关键。选对了，故障特征一目了然；选错了，啥也看不到。
4.2 短时傅里叶变换（STFT）：给频谱加上时间轴
FFT有个硬伤：它假设信号是平稳的。但风机启动、停机、变工况时，信号是非平稳的。频率成分随时间变化。
STFT的思路很朴素：把信号切成一小段一小段，假设每一段内是平稳的，然后对每一段做FFT。最后把结果拼起来，就得到了一个“时间-频率”的二维图谱。
数学上，STFT的定义是：
STFT{x(t)}(τ, ω) = ∫ x(t) · w(t-τ) · e^(-jωt) dt其中 w(t) 是窗函数。窗函数的选择，直接决定了时频分辨率。
窗函数的选择：
- 窗越窄：时间分辨率高，能看清信号在何时变化。但频率分辨率低，频带变宽。
- 窗越宽：频率分辨率高，能看清细微的频率差异。但时间分辨率低，看不清变化时刻。
这就是海森堡测不准原理在信号处理里的体现。你不可能同时获得无限高的时间和频率分辨率。
我个人习惯，分析风机启动过程的振动信号时，会用较窄的窗（比如汉宁窗，长度128点），因为我想看清转速变化时，各频率成分出现的时刻。而分析稳态运行信号时，我会用较宽的窗（比如512点），追求频率精度。
4.3 小波变换入门：多分辨率分析的利器
STFT的窗宽是固定的。这带来一个问题：低频信号需要长窗口才能看清频率，高频信号需要短窗口才能定位时间。STFT做不到“自适应”。
小波变换就是来解决这个问题的。
它用一组可伸缩、可平移的“小波基”来代替STFT里的固定窗函数。低频时，小波自动变宽，获得高频率分辨率；高频时，小波自动变窄，获得高时间分辨率。
我画了一张图，帮你理解小波变换的核心逻辑：
你看，小波变换就像一把“变焦镜头”。看全局时用广角，看细节时用长焦。这在分析风机振动信号时特别有用。比如，一个冲击信号，在时域上可能只持续几毫秒，但在频域上覆盖很宽的频带。小波变换能同时捕捉到它的发生时刻和频率成分。
注意：小波基函数的选择没有标准答案。不同的小波基（如db4、db8、symlet、morlet）对同一信号的分析结果可能不同。我建议你多试几种，结合你的物理直觉来判断哪个结果更合理。
4.4 Python实战：动手实现包络谱与STFT
光说不练假把式。我们直接上代码。下面是一个完整的Python示例，模拟一个轴承外圈故障信号，并做包络谱分析和STFT时频图。
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.fft import fft, fftfreq
# 1. 生成模拟信号：轴承外圈故障
fs = 10000 # 采样率 10kHz
t = np.arange(0, 1, 1/fs)
# 故障频率 100Hz，载波频率 2000Hz
fault_freq = 100
carrier_freq = 2000
# 模拟冲击序列
impulses = np.zeros_like(t)
impulse_interval = int(fs / fault_freq)
impulses[::impulse_interval] = 1
# 生成载波信号
carrier = np.sin(2 * np.pi * carrier_freq * t)
# 调制：冲击激励载波
signal_raw = impulses * carrier
# 加一点噪声
signal_raw += 0.3 * np.random.randn(len(t))
# 2. 包络谱分析
# 带通滤波：围绕载波频率
b, a = signal.butter(4, [1800, 2200], btype='band', fs=fs)
filtered = signal.filtfilt(b, a, signal_raw)
# Hilbert变换求包络
analytic_signal = signal.hilbert(filtered)
envelope = np.abs(analytic_signal)
# 对包络做FFT
N = len(envelope)
yf = fft(envelope)
xf = fftfreq(N, 1/fs)[:N//2]
envelope_spectrum = 2.0/N * np.abs(yf[:N//2])
# 3. STFT时频分析
f_stft, t_stft, Zxx = signal.stft(signal_raw, fs, nperseg=256)
# 4. 绘图
plt.figure(figsize=(12, 8))
# 原始信号
plt.subplot(3, 1, 1)
plt.plot(t[:500], signal_raw[:500])
plt.title('原始信号 (时域)')
plt.xlabel('时间 (s)')
plt.ylabel('幅值')
# 包络谱
plt.subplot(3, 1, 2)
plt.plot(xf[:500], envelope_spectrum[:500])
plt.title('包络谱')
plt.xlabel('频率 (Hz)')
plt.ylabel('幅值')
plt.axvline(x=fault_freq, color='r', linestyle='--', label='故障频率 100Hz')
plt.legend()
# STFT时频图
plt.subplot(3, 1, 3)
plt.pcolormesh(t_stft, f_stft, np.abs(Zxx), shading='gouraud')
plt.title('STFT时频图')
plt.xlabel('时间 (s)')
plt.ylabel('频率 (Hz)')
plt.colorbar(label='幅值')
plt.tight_layout()
plt.show()
运行结果解读：
- 包络谱图上，你应该能在100Hz处看到一个明显的峰值。这就是我们模拟的故障频率。
- STFT时频图上，你会看到在2000Hz附近有一条水平的亮线，代表载波频率。同时，在100Hz的间隔上，有垂直的条纹，代表冲击的周期性发生。
这段代码我用了很多次。每次调试新算法，我都会先跑一遍这个模拟信号，确认流程没问题，再上真实数据。你想想看，如果连模拟信号都跑不出预期结果，那真实数据就更难说了。
4.5 三种方法的对比与选择
最后，我帮你梳理一下这三种方法的应用场景。别贪多，用对地方才是关键。
| 方法 | 适用场景 | 优点 | 缺点 |
|---|---|---|---|
| 包络谱分析 | 轴承、齿轮等调制故障诊断 | 能提取低频故障特征，抗干扰能力强 | 需要先知道共振频带，依赖带通滤波 |
| 短时傅里叶变换 | 变工况、启停机过程的时频分析 | 物理意义明确，实现简单 | 窗宽固定，时频分辨率无法兼顾 |
| 小波变换 | 非平稳信号、瞬态冲击、多尺度分析 | 多分辨率，自适应时频分析 | 小波基选择困难，计算量较大 |
我个人习惯是：先做包络谱，快速判断有没有调制故障。如果信号平稳，用STFT看全貌。如果信号复杂，或者想深挖细节，再上小波变换。别一上来就搞复杂的方法，简单方法往往最有效。
好了，这一章的内容就到这里。信号处理是故障诊断的基石，这些方法你多练几次，就会形成肌肉记忆。下次遇到真实的风机数据，你自然知道该用哪把刀。