振动分析实验：滤波、功率谱、包络与时频方法

实验实践更新于 2026-10-02编辑本文

本页承接测量到诊断。所有数值实验均为合成数据，用于验证计算方法；没有硬件接入、现场故障标签或工程验收结论。

1. 选择方法之前先写出目标

目标	可采用的方法	必须说明的代价
抑制已知带外分量	低通、高通、带通、带阻	幅相变化、过渡带、边界瞬态
估计随机背景的频率分布	PSD、Welch 平均	频率分辨与估计方差的折中
观察冲击的重复规律	共振带选择与包络	频带选择、干扰、转速变化
观察频率随时间改变	STFT、连续小波	时间与频率分辨的折中
分析不同尺度	离散小波	小波基、分解级数、边界延拓
抵消可观测相关噪声	自适应滤波	独立参考通道、收敛与稳定条件

处理选择应对应假设，不能默认“滤波越多、曲线越平滑、诊断越准确”。

2. 滤波同时改变信号

LTI 滤波器满足 $Y(f)=H(f)X(f)$ ：噪声与有用成分只要位于同一频带，就都会受 $H(f)$ 影响。Butterworth 的平滑幅频响应并不意味着对所有波形都保真。

二阶节 SOS 表示适合数值实现。前后向滤波在边缘之外近似零相位，但其合成幅频响应为 $|H(f)|^2$ ，不能把单程的截止点增益直接当双程结果；它需要未来样本，属于离线处理。短记录的填充和首尾瞬态也必须检查。SciPy sosfiltfilt

带冲击的合成波形、低通结果与中值滤波结果

图中中值滤波削弱尖峰，低通抑制快速变化。如果尖峰来自轴承冲击，这也意味着诊断信息被削弱。不能仅因处理后更平滑便称其“正确去噪”。

移动平均也是滤波器。长度 $L$ 的等权移动平均在采样索引上通常带来 $(L-1)/2$ 的群延迟；采用 valid 卷积还会缩短输出，绘图时要对齐时间轴。推导见滤波与采集。

3. PSD 为什么能积分得到均方值

对实信号的双边功率谱密度 $S_{qq}(f)$ ，在一致的归一化条件下：

E[q^2]=\int_{-\infty}^{\infty}S_{qq}(f)\,df.

单边 PSD 把正负频率的贡献合到非负频率，因此带内均方值近似为：

q_{\mathrm{RMS,band}}^2 \approx\sum_{k\in\mathrm{band}}S^{(1)}_{qq}(f_k)\Delta f.

若输入是 mm/s，则 PSD 的单位是 $(\mathrm{mm/s})^2/\mathrm{Hz}$ ，乘 Hz 后成为均方速度，再开方才得到 mm/s。去均值会移除直流贡献；Welch 分段、加窗、平均也使估计值与有限记录的直接 RMS 不必完全相等。SciPy Welch

下面的两个整周期正弦不含噪声，选定记录长度使 Welch 各段也覆盖整周期，可检查缩放是否一致。

import numpy as np
from scipy import signal

fs = 4096.0
t = np.arange(8192) / fs
x = np.sin(2*np.pi*100*t) + 0.5*np.sin(2*np.pi*300*t)

f, psd = signal.welch(
    x, fs=fs, window="hann", nperseg=1024,
    noverlap=512, detrend=False, scaling="density"
)
df = f[1] - f[0]
mean_square_from_psd = np.sum(psd) * df
mean_square_direct = np.mean(x*x)
assert np.isclose(mean_square_from_psd, 0.625)
assert np.isclose(mean_square_direct, 0.625)
print("RMS:", np.sqrt(mean_square_from_psd))

对非整周期窄峰做带内积分时，需要考虑窗泄漏进入相邻频点的能量。PSD 峰高也随频率分辨和窗变化，不应直接当作正弦幅值。

4. 包络解调在测什么

对幅度缓慢变化的窄带模型：

x(t)=[1+m\cos(2\pi f_mt)]\cos(2\pi f_ct), \quad 0<m<1,\quad f_m\ll f_c,

展开乘积得到载频 $f_c$ 以及边频 $f_c\pm f_m$ 。在适当的频谱分离条件下，解析信号

z(t)=x(t)+j\mathcal H\{x(t)\}

的模近似为包络 $e(t)=|z(t)|=1+m\cos(2\pi f_mt)$ 。再对 $e(t)$ 去均值并求谱，可看到调制频率 $f_m$ 。SciPy Hilbert

轴承冲击激起的共振衰减不是严格的单一调幅余弦，但“高频振荡由低频重复事件调制”的思路相似。必须先选合适频带，再解释包络谱。公式边界见包络分析。

import numpy as np
from scipy import signal

fs = 4096.0
t = np.arange(8192) / fs
x = (1 + 0.4*np.cos(2*np.pi*30*t)) * np.cos(2*np.pi*1000*t)
sos = signal.butter(4, [700, 1300], btype="bandpass", fs=fs, output="sos")
filtered = signal.sosfiltfilt(sos, x)
envelope = np.abs(signal.hilbert(filtered))

# 排除本教学记录两端的滤波/解析变换边界；实际长度需按滤波器验证。
middle = envelope[512:-512]
middle = middle - middle.mean()
window = signal.windows.hann(middle.size, sym=False)
spectrum = np.abs(np.fft.rfft(middle * window))
f = np.fft.rfftfreq(middle.size, 1/fs)
band = (f > 5) & (f < 100)
peak_hz = f[band][np.argmax(spectrum[band])]
assert abs(peak_hz - 30) <= fs / middle.size
print("包络谱主峰约", peak_hz, "Hz")

合成重复冲击及其 Hilbert 包络，重复频率为 30 Hz

这张旧教学图展示冲击模型，而上面的代码使用便于验证的调幅模型。两者都不是某个实际轴承的故障测量。30 Hz 也只是演示参数，不能直接称为 BPFI 或 BPFO。

5. STFT 与小波的尺度

STFT 把信号分成局部加窗记录：

X(m,\omega)=\sum_n x[n]w[n-mR]e^{-j\omega n},

其中 $R$ 是帧移。长窗通常有利于区分相近频率，却会混合更长时间内的变化；短窗有利于定位事件时间，却使频率主瓣变宽。重叠增加时间取样密度，不消除这种折中。

离散小波用低通、高通滤波与抽取逐级分解。db4 的“五层分解”返回 $A_5,D_5,D_4,\ldots,D_1$ ，不是六路采样率相同的时域信号；系数索引也不能直接标成原始时间。边界延拓会影响长度和端点。PyWavelets DWT

db4 五层离散小波分解，横轴为各层系数索引

近似理想滤波器时，细节 $D_j$ 对应频带约为 $[f_s/2^{j+1},f_s/2^j]$ ；真实小波滤波器有过渡带和频谱重叠，不能把这当精确砖墙分频。

import numpy as np
import pywt

fs = 1000.0
t = np.arange(1024) / fs
x = np.sin(2*np.pi*10*t) + 0.5*np.sin(2*np.pi*100*t)
coeffs = pywt.wavedec(x, "db4", mode="symmetric", level=5)
rebuilt = pywt.waverec(coeffs, "db4", mode="symmetric")[:x.size]
assert np.allclose(rebuilt, x)
print("各层长度:", [len(c) for c in coeffs])

去噪还需要明确阈值估计、软／硬阈值、处理哪些层以及重构误差。分解成功不等于已成功去噪。连续小波与 STFT 采用不同的时频取舍，小波并非在所有频率都具有更高分辨率。Hilbert–Huang 方法还涉及经验模态分解、模态混叠和端点效应，应作为另一个有条件的方法研究，而非自动确诊步骤。

6. 自适应滤波需要参考信息

噪声抵消模型中，主通道 $d[n]=s[n]+v[n]$ ，另一个参考通道 $u[n]$ 应与噪声 $v[n]$ 相关、与目标 $s[n]$ 尽量不相关。令 $\mathbf u_n$ 为最近 $L$ 个参考样本：

y[n]=\mathbf w_n^\mathsf T\mathbf u_n,\qquad e[n]=d[n]-y[n].

最小化 $e[n]^2/2$ 的梯度为 $-e[n]\mathbf u_n$ ，所以 LMS 更新为 $\mathbf w_{n+1}=\mathbf w_n+\mu e[n]\mathbf u_n$ 。归一化 LMS 则把步长除以参考能量：

\mathbf w_{n+1}=\mathbf w_n+ \frac{\mu}{\varepsilon+\|\mathbf u_n\|^2}e[n]\mathbf u_n.

这里 $y$ 是估计噪声， $e$ 才是保留目标的候选输出。原稿把未知“纯信号”当现成输入，不能代表真实无参考的去噪。步长、延时、相关性与稳定性必须另行验证；参考混入目标时也会抵消真正的振动成分。

7. 高阶谱与模态各解决什么

零均值平稳信号的三阶累积量为 $C_3(\tau_1,\tau_2)=E[x(t)x(t+\tau_1)x(t+\tau_2)]$ 。其双变量傅里叶变换是双谱：

B(f_1,f_2)=\iint C_3(\tau_1,\tau_2) e^{-j2\pi(f_1\tau_1+f_2\tau_2)}\,d\tau_1d\tau_2.

它研究三频率间的相位耦合等统计特性，不是“把 FFT 幅值平方”；后者属于二阶功率谱的相关运算。有限样本估计需分段、平均、频率区域与显著性控制，不能只因双谱有峰便宣称某种故障。

模态分析关注结构的固有频率、阻尼和振型，常由输入力与响应估计 FRF 后拟合；与运行响应的 FFT、ODS 不同。继续阅读ODS 与模态分析。