4. 基于FFT的轴承故障特征频率提取

任务说明

对轴承振动信号做FFT，提取幅值谱，定位外圈故障特征频率（BPFO）及其谐波，并对比正常/故障两种状态的频谱差异。核心步骤：生成/读取振动信号 → 去均值 → 加窗 → FFT → 幅值谱 → 峰值检测 → 与理论BPFO比对。

依赖与环境

pip install numpy scipy matplotlib

实现代码

import numpy as np
from scipy.fft import rfft, rfftfreq
from scipy.signal import find_peaks, windows

FS = 12000          # 采样率 Hz
DUR = 1.0           # 信号时长 s
N = int(FS * DUR)

# 轴承参数（外圈故障特征频率 BPFO）
N_BALL, D_PITCH, D_BALL, CONTACT = 9, 39.04, 7.94, 0.0  # mm
def bpfo(fr):
    return (N_BALL / 2) * fr * (1 - (D_BALL / D_PITCH) * np.cos(np.radians(CONTACT)))

FR = 25.0           # 转频 Hz
F_BPFO = bpfo(FR)   # ≈ 107.4 Hz

def make_signal(fault=False, seed=0):
    rng = np.random.default_rng(seed)
    t = np.arange(N) / FS
    sig = 0.5 * np.sin(2 * np.pi * FR * t)          # 转频分量
    sig += 0.05 * rng.standard_normal(N)            # 背景噪声
    if fault:
        # 外圈故障：BPFO 及其 2、3 次谐波，带衰减冲击调制
        for k, amp in zip((1, 2, 3), (1.0, 0.5, 0.25)):
            sig += amp * np.sin(2 * np.pi * k * F_BPFO * t)
    return t, sig

def amplitude_spectrum(x):
    x = x - x.mean()
    xw = x * windows.hann(len(x))
    X = rfft(xw)
    f = rfftfreq(len(x), 1 / FS)
    mag = 2.0 * np.abs(X) / np.sum(windows.hann(len(x)))
    return f, mag

def detect_peaks(f, mag, fmin=20, fmax=1000, top=5):
    band = (f >= fmin) & (f <= fmax)
    fb, mb = f[band], mag[band]
    idx, _ = find_peaks(mb, height=mb.max() * 0.2)
    order = np.argsort(mb[idx])[::-1][:top]
    return fb[idx[order]], mb[idx[order]]

def report(name, sig):
    f, mag = amplitude_spectrum(sig)
    pf, pm = detect_peaks(f, mag)
    print(f"[{name}] 理论 BPFO = {F_BPFO:.2f} Hz")
    for fi, mi in zip(pf, pm):
        tag = " <== BPFO" if abs(fi - F_BPFO) < 5 else ""
        print(f"  峰值 {fi:8.2f} Hz  幅值 {mi:.4f}{tag}")
    return f, mag

if __name__ == "__main__":
    _, s_norm = make_signal(fault=False, seed=1)
    _, s_fault = make_signal(fault=True, seed=1)
    f, m_norm = report("正常", s_norm)
    _, m_fault = report("故障", s_fault)

    # 判定：故障信号在 BPFO 附近幅值应显著高于正常信号
    band = (f > F_BPFO - 5) & (f < F_BPFO + 5)
    ratio = m_fault[band].max() / (m_norm[band].max() + 1e-12)
    print(f"BPFO 幅值比(故障/正常) = {ratio:.2f}  ->  {'通过' if ratio > 5 else '未通过'}")

运行与验证

python fft_bearing.py

观察输出：故障信号在约 107 Hz（BPFO）及其 214 Hz、321 Hz 谐波处出现明显峰值，正常信号在该频段仅有噪声底。判定条件：BPFO 频带内故障/正常幅值比 > 5 即视为特征频率成功提取。

处理流程

关键点

• 加 Hann 窗抑制频谱泄漏，幅值按窗函数和归一化，保证峰值幅值可比。
• BPFO 由轴承几何参数与转频计算，是判定故障的参考基准。
• 故障信号在 BPFO 及谐波处能量集中，正常信号仅含转频与噪声底，二者频谱差异即为识别依据。