作为一名在工业诊断领域摸爬滚打多年的工程师，我深知旋转机械就像工厂的“心脏”，它的每一次异常跳动都可能预示着巨大的隐患。而要从嘈杂的振动信号中，精准地捕捉到这些“心跳异常”，频域分析无疑是我们手中的“听诊器”。今天，我就和大家深入聊聊，如何运用MATLAB这把“瑞士军刀”，在频域世界里高效地揪出旋转机械的故障特征。

一、为何选择频域分析？从时域的“一团麻”到频域的“清晰谱”

想象一下，你拿到了一段从风机轴承座上采集到的振动加速度信号。在时域图里，它可能就是一串上下剧烈波动、看起来杂乱无章的曲线，就像一团纠缠不清的毛线。故障信息（比如微弱的冲击）被淹没在强大的背景噪声和工频振动中，肉眼几乎无法分辨。

这时，傅里叶变换（FFT）就派上用场了。它就像一台神奇的“信号分解机”，能把时域里那团“毛线”按照不同频率成分“梳理”开来，转换成一张“频谱图”。在这张图上，横轴是频率，纵轴是振幅（或功率）。于是，工频（转频）及其倍频（谐波）、轴承的故障特征频率、齿轮的啮合频率等，都会以一个个突起的“谱峰”形式清晰地呈现在我们面前。故障特征频率的出现及其幅值的增长，直接对应着机械部件的损伤，比如轴承滚珠剥落会产生特定的通过频率，齿轮断齿会引发啮合频率的调制边带。频域分析让我们从宏观的“混乱”走向微观的“秩序”，实现了故障特征的分离与凸显。

二、MATLAB频域分析核心技巧与完整示例

掌握了“为什么”，接下来就是“怎么做”。MATLAB的信号处理工具箱为我们提供了强大而便捷的工具链。下面，我将通过一个完整的示例，演示从数据加载到特征提取的全过程。

技术栈声明： 本文所有示例均基于 MATLAB 及其内置的 Signal Processing Toolbox。

假设我们有一段采样频率为 fs = 10000 Hz，时长2秒的振动信号 vibrationData，它来自一台转速约为3000 RPM（转频 fr = 50 Hz）的电机驱动泵，我们怀疑其驱动端轴承存在早期外圈故障。

% ========== 示例：轴承外圈故障特征提取 ==========
% 作者：资深故障诊断工程师
% 功能：演示使用MATLAB进行频域分析，提取轴承外圈故障特征频率
% 技术栈：MATLAB + Signal Processing Toolbox

%% 1. 模拟数据生成（在实际工作中，这部分由数据采集系统获得）
clear; clc; close all;
fs = 10000;             % 采样频率，10 kHz
t = 0:1/fs:2-1/fs;     % 时间向量，2秒时长
N = length(t);         % 信号点数

% 模拟正常振动：工频（50Hz）及其二次谐波，加随机噪声
fr = 50;               % 转频，50 Hz (3000 RPM / 60)
normal_vib = 1.5*sin(2*pi*fr*t) + 0.6*sin(2*pi*2*fr*t + pi/4);
noise = 0.8 * randn(1, N); % 高斯白噪声

% 模拟轴承外圈故障冲击：故障特征频率（BPFO）假设为 3.1倍转频
BPFO = 3.1 * fr;       % 外圈故障特征频率，155 Hz
% 生成周期性冲击序列（每1/BPFO秒一次冲击）
impulse_train = zeros(1, N);
impulse_interval = round(fs / BPFO);
impulse_train(1:impulse_interval:N) = 3.0; % 在冲击点置入幅值3
% 对冲击序列进行低通滤波，模拟实际冲击响应
[b, a] = butter(4, 500/(fs/2), 'low'); % 4阶低通，截止频率500Hz
fault_impulse = filter(b, a, impulse_train);

% 合成最终振动信号：正常振动 + 故障冲击 + 噪声
vibrationData = normal_vib + fault_impulse + noise;

%% 2. 时域初步观察
figure('Name', '时域信号观察', 'NumberTitle', 'off');
subplot(2,1,1);
plot(t, vibrationData);
xlabel('时间 (s)');
ylabel('振幅');
title('原始振动信号时域波形');
grid on;
% 局部放大观察冲击（查看0.1到0.2秒区间）
subplot(2,1,2);
zoom_start = find(t>=0.1, 1);
zoom_end = find(t>=0.2, 1);
plot(t(zoom_start:zoom_end), vibrationData(zoom_start:zoom_end));
xlabel('时间 (s)');
ylabel('振幅');
title('局部放大波形 (0.1s - 0.2s) - 试图观察冲击');
grid on;
% 评论：在时域中，微弱的周期性冲击被强大的工频成分和噪声淹没，难以直接识别。

%% 3. 基本的FFT频谱分析
L = N; % 使用整个数据长度
Y = fft(vibrationData, L); % 执行FFT
P2 = abs(Y/L);          % 双侧频谱
P1 = P2(1:L/2+1);       % 取单侧频谱（正频率部分）
P1(2:end-1) = 2*P1(2:end-1); % 幅值修正（除直流分量外）
f = fs*(0:(L/2))/L;     % 频率向量

figure('Name', '基本FFT幅值谱', 'NumberTitle', 'off');
plot(f, P1);
xlabel('频率 (Hz)');
ylabel('幅值');
title('振动信号的单边幅值谱');
xlim([0, 500]); % 聚焦在0-500Hz低频段
grid on;
hold on;
% 标记理论转频和故障频率位置
plot([fr, fr], [0, max(P1)], 'r--', 'LineWidth', 1.2, 'DisplayName', sprintf('转频 %.1f Hz', fr));
plot([BPFO, BPFO], [0, max(P1)], 'g--', 'LineWidth', 1.2, 'DisplayName', sprintf('BPFO %.1f Hz', BPFO));
legend('Location', 'best');
hold off;
% 评论：在基本谱图中，50Hz和100Hz的工频及其谐波非常突出。155Hz处的BPFO峰虽然存在，但相对于工频峰和背景噪声，并不十分显著，容易被忽略。

%% 4. 功率谱密度（PSD）估计 - 使用pwelch方法，效果更佳
% 关联技术详解：pwelch采用Welch平均周期图法，通过分窗、重叠、平均，能有效平滑随机噪声，得到更稳定的频谱估计，是工程实践中的首选。
figure('Name', 'Welch功率谱密度估计', 'NumberTitle', 'off');
window = hann(round(fs/2)); % 汉宁窗，窗长0.5秒对应数据点
noverlap = round(length(window)*0.5); % 50%重叠
nfft = max(256, 2^nextpow2(length(window))); % FFT点数

[pxx, f_welch] = pwelch(vibrationData, window, noverlap, nfft, fs);
plot(f_welch, 10*log10(pxx)); % 纵坐标转换为分贝(dB)刻度，便于观察小信号
xlabel('频率 (Hz)');
ylabel('功率/频率 (dB/Hz)');
title('Welch方法估计的功率谱密度 (PSD)');
xlim([0, 500]);
grid on;
hold on;
plot([fr, fr], ylim, 'r--', 'LineWidth', 1.2);
plot([BPFO, BPFO], ylim, 'g--', 'LineWidth', 1.2);
legend('PSD', sprintf('转频 %.1f Hz', fr), sprintf('BPFO %.1f Hz', BPFO), 'Location', 'best');
hold off;
% 评论：PSD图在分贝刻度下，动态范围更广。可以观察到155Hz处有一个清晰的谱峰，明显高于周围的噪声基底，故障特征的可辨识度大大增强。

%% 5. 包络谱分析 - 提取冲击特征的神器
% 关联技术详解：轴承故障的冲击会调制在高频共振带上。包络分析（也称解调分析）先通过带通滤波隔离共振频带，然后提取信号的包络（即幅值变化），再对包络做频谱分析。这个包络谱能直接清晰地展现故障冲击的频率（如BPFO）。
figure('Name', '包络谱分析', 'NumberTitle', 'off');
% 步骤1：带通滤波，选取一个共振频带（例如2000-3000Hz）。这里我们由于模拟数据限制，对原信号进行希尔伯特变换直接求包络。
analytic_signal = hilbert(vibrationData); % 希尔伯特变换得到解析信号
envelope = abs(analytic_signal);         % 解析信号的模即为包络线

% 步骤2：对包络信号进行频谱分析（同样使用pwelch）
[pxx_env, f_env] = pwelch(envelope-mean(envelope), window, noverlap, nfft, fs); % 减去均值去除直流
subplot(2,1,1);
plot(t, envelope);
xlabel('时间 (s)');
ylabel('包络幅值');
title('振动信号的包络线');
grid on;
xlim([0, 0.5]); % 显示前0.5秒

subplot(2,1,2);
plot(f_env, 10*log10(pxx_env));
xlabel('频率 (Hz)');
ylabel('功率/频率 (dB/Hz)');
title('包络信号的功率谱（包络谱）');
xlim([0, 500]);
grid on;
hold on;
% 在包络谱中，故障特征频率及其倍频应该非常突出
plot([BPFO, BPFO], ylim, 'g--', 'LineWidth', 2, 'DisplayName', sprintf('BPFO %.1f Hz', BPFO));
plot([2*BPFO, 2*BPFO], ylim, 'm--', 'LineWidth', 1.5, 'DisplayName', sprintf('2*BPFO %.1f Hz', 2*BPFO));
plot([fr, fr], ylim, 'r:', 'LineWidth', 1, 'DisplayName', sprintf('转频 %.1f Hz', fr));
legend('Location', 'best');
hold off;
% 评论：包络谱是诊断轴承、齿轮局部损伤类故障的利器。图中，在155Hz (BPFO) 和310Hz (2*BPFO) 处出现了明显的谱峰，而工频（50Hz）成分很弱。这强有力地指示了轴承外圈存在周期性冲击，即故障特征。

三、高级技巧与实战注意事项

掌握了基本流程后，一些高级技巧和细节能让你如虎添翼。

1. 分辨率与量程的权衡： FFT的频率分辨率 df = fs / N。要分辨两个靠近的频率（比如转频的边带），需要足够长的数据（N大）。但N太大会增加计算量，且可能包含非平稳过程。通常通过pwelch中的窗长来控制分辨率。频率量程（能看到的最大频率）由奈奎斯特定律决定，为 fs/2。确保 fs > 2 * f_max，f_max是你关心的最高频率。

2. 窗函数的选择： FFT默认假设信号是周期性的。如果截取的不是整数个周期，就会发生“频谱泄漏”，导致能量分散到旁频。加窗（如汉宁窗、汉明窗）可以抑制泄漏，但会降低频率分辨率并轻微改变幅值。pwelch函数内部已经处理了加窗的影响。

3. 阶次分析： 对于变速运行的设备（如启动/停机过程），频率是变化的，固定频率的谱线会模糊。此时需要将时域信号重采样为与轴旋转角度等间隔的“角域信号”，再进行FFT，得到“阶次谱”。这需要转速脉冲信号（键相信号）。MATLAB的 orderAnalysis 函数族可以完成此工作。

4. 自动化特征提取： 在大规模监测中，需要自动识别谱峰。可以使用 findpeaks 函数。

% 在包络谱中自动寻找故障频率峰值示例
[peak_amps, peak_locs] = findpeaks(pxx_env, f_env, 'MinPeakHeight', max(pxx_env)/20, 'MinPeakDistance', BPFO*0.8);
% 找到最接近理论BPFO的峰值
[~, idx] = min(abs(peak_locs - BPFO));
if abs(peak_locs(idx) - BPFO) / BPFO < 0.05 % 频率误差在5%以内
    fprintf('警报：检测到显著的轴承外圈故障特征频率在 %.2f Hz，幅值为 %.2e。\n', peak_locs(idx), peak_amps(idx));
end

四、应用场景、优缺点与总结

应用场景： 频域分析是旋转机械状态监测与故障诊断（PdM）的基石。它广泛应用于：

  • 风力发电机： 监测齿轮箱、主轴承、发电机的齿轮啮合故障、轴承损伤、电气不平衡。
  • 大型离心压缩机/泵： 诊断转子不平衡、不对中、喘振、叶片通过频率异常。
  • 机床主轴： 检测刀具磨损、主轴轴承损坏、动平衡失效。
  • 轨道交通： 分析列车轮对轴承、齿轮箱的健康状态。

技术优缺点：

  • 优点：
    1. 直观清晰： 将复杂的时域波形分解为独立的频率成分，故障特征一目了然。
    2. 抗干扰能力强： 通过PSD、包络谱等方法，能有效抑制随机噪声，凸显周期性故障特征。
    3. 成熟可靠： 理论基础坚实，在工业界有数十年的成功应用历史，是标准化的诊断手段。
    4. 便于量化： 特征频率、幅值、边带间距等均可精确测量，用于设定报警阈值和趋势分析。
  • 缺点与局限：
    1. 对非平稳信号乏力： 传统的FFT假设信号是平稳的。对于剧烈变转速、瞬态冲击过程，需要结合阶次分析、时频分析（如小波变换、短时傅里叶变换）。
    2. 对早期微弱故障不敏感： 当故障特征幅值极低时，可能被噪声完全掩盖，需要更先进的降噪和特征增强技术。
    3. 需要先验知识： 为了识别特征频率，需要知道设备的结构参数（如轴承几何尺寸、齿数）和实时转速。

注意事项：

  1. 数据质量是根本： 确保传感器安装正确、接线可靠、采样参数（fs， 量程）设置合理。垃圾数据进，垃圾结果出。
  2. 理解设备机理： 必须清楚你要诊断的对象的理论特征频率计算公式，以及不同故障在频谱上的典型表现模式（如调制边带、谐波簇）。
  3. 建立基线频谱： 在设备健康状态下采集频谱作为“基线”或“指纹”。后续的频谱与之对比，变化量往往比绝对值更有意义。
  4. 综合判断： 不要孤立地看待一个谱峰。结合时域指标（如峰值、峭度）、多测点信息、设备运行工况（负载、温度）进行综合诊断。

文章总结： 总而言之，MATLAB为旋转机械的频域故障诊断提供了一个极其强大且灵活的环境。从基础的FFT到更稳健的Welch PSD估计，再到针对冲击性故障的“杀手锏”——包络谱分析，我们拥有清晰的路径去穿透噪声，直达故障的本质。关键在于理解每种方法背后的物理和数学意义，并将其与具体的机械故障机理相结合。通过本文提供的完整示例和实战技巧，希望你能快速上手，将MATLAB的频域分析能力转化为保障设备安全稳定运行的可靠工具。记住，诊断既是科学，也是艺术，需要不断的实践和经验的积累。祝你在这个有趣的领域不断发现，精准排故！