企业官网建设流程全解析

简介：这份MATLAB代码资源面向轴承故障诊断与信号处理方向的工程师、研究生及设备维护人员，针对滚动轴承早期故障特征提取与识别问题，提供了基于Hilbert包络谱、Haar小波、数学形态学、时域无量纲参数及FFT五种分析思路的可执行脚本。资源压缩包共31个文件，大小仅1.04MB，其中包含5个.m主程序代码、11个.dat原始振动数据、14个.jpg结果图以及1个README说明文档，代码与数据一一对应，便于对照验证。这些脚本覆盖从时域指标计算到频域包络解调、从小波多分辨率分析到形态学滤波的完整链路，可直接用于外圈故障、内圈故障等典型工况的数据处理与特征可视化。目前已有188人学习下载，适合希望快速上手故障诊断算法、需要可运行示例参考的读者。

1. 轴承故障诊断 MATLAB 代码：从振动信号到故障判据的完整链路

采集系统报出「外圈疑似故障」的告警，现场却没人敢立刻停机，这种两难在设备维护里很常见。轴承故障诊断 MATLAB 代码解决的就是这个判断问题：振动信号里隐藏着可与几何参数对应起来的故障特征频率，把采集数据经过去噪、包络解调、频谱分析，就能落到一组可判读的峰值与指标上。这套代码链路覆盖外圈、内圈、滚动体、保持架四种故障频率的计算，并给出包络谱分析与时域特征分类的完整实现。适合设备状态监测工程师、故障诊断方向的研究生，以及想用手头振动数据快速验证诊断结论的 MATLAB 用户。

2. 轴承故障特征频率计算：从几何参数到可复用的 MATLAB 函数

2.1 四种特征频率的物理含义与公式

轴承故障诊断不是靠「听声音」或「看波形」，而是基于滚动轴承的运动学关系。滚珠滚过外圈或内圈的局部缺陷时会产生周期性冲击，冲击频率只取决于几何尺寸和转速，与传感器安装位置无关。四种特征频率的物理含义不同：外圈故障频率 BPFO 是滚珠依次通过外圈固定缺陷点的频率，外圈固定在轴承座上，谱图表现最稳定；内圈故障频率 BPFI 对应缺陷在内圈上与滚珠接触的频率，内圈随轴旋转，数值比 BPFO 高；滚动体故障频率 BSF 对应滚动体自转一周中缺陷与内外圈接触的频率；保持架故障频率 FTF 数值最低，一般小于 0.5 倍转频。

计算只需要四个几何参数加一个转速：滚动体个数 nb、滚动体直径 d、节圆直径 D、接触角 alpha，以及轴转频 fr：

BPFO = (nb*fr/2) * (1 - (d/D)*cos(alpha)) BPFI = (nb*fr/2) * (1 + (d/D)*cos(alpha)) BSF = (D*fr/(2*d)) * (1 - ((d/D)*cos(alpha))^2) FTF = (fr/2) * (1 - (d/D)*cos(alpha))

这套公式的前提是滚动体纯滚动、无滑移。实际运行中接触角随载荷变化，计算结果与真实故障频率之间通常有 1%～2% 的偏差，所以后续谱峰搜索必须预留容差带宽。

2.2 把公式封装成 MATLAB 函数

下面这个函数把四个公式封装在一起，是整条诊断代码链路的起点：

function [bpfo, bpfi, bsf, ftf] = bearingFaultFreq(nb, d, D, alpha, fr) % 计算滚动轴承四种故障特征频率 % 输入: % nb - 滚动体个数 % d - 滚动体直径 (m) % D - 节圆直径 (m) % alpha - 接触角 (rad) % fr - 转频 (Hz), 由转速 rpm/60 得到 % 输出: % bpfo, bpfi, bsf, ftf - 外圈、内圈、滚动体、保持架故障频率(Hz) bpfo = nb * fr / 2 * (1 - d/D*cos(alpha)); bpfi = nb * fr / 2 * (1 + d/D*cos(alpha)); bsf = D*fr / (2*d) * (1 - (d/D*cos(alpha))^2); ftf = fr / 2 * (1 - d/D*cos(alpha)); end

调用时注意单位统一：d 和 D 全部用毫米或全部用米，公式里只有比值，单位不影响结果；alpha 必须用弧度，直接传角度数是常见错误。转速换算同样容易出错：给定 rpm 时，fr = rpm / 60。

2.3 参数取值表与工程边界

以深沟球轴承 6205 为例，一组公开数据集中常见的参数如下：

对应转频 fr ≈ 29.95 Hz，代入函数得到：BPFO ≈ 107.4 Hz，BPFI ≈ 162.2 Hz，BSF ≈ 70.6 Hz，FTF ≈ 11.9 Hz。注意：也有资料把滚动体故障频率写成 141.2 Hz，那是按「缺陷同时接触内外圈」把冲击频率计成两倍的约定，两种定义都能用，关键是整个项目保持一致，否则谱峰搜索会整体偏掉。

边界情况再提两点。接触角为零时公式退化为简化形式，这是深沟球轴承的默认假设，角接触球轴承必须给实际接触角；如果四个特征频率之间存在近似整数倍关系，多半是几何参数有误，回轴承型号手册核对。实际诊断的搜索带宽建议取 ±1.5%，低频段因频率分辨率限制可放宽到 ±3%。

提示：把特征频率计算独立成函数而不是写进主脚本，批量诊断不同型号轴承时只需要改一行调用。

3. 振动信号预处理：采样率校验、去直流与带通滤波

3.1 数据读取与采样率校验

现场采集和公开数据集最常见的格式是 .mat、CSV 和纯文本。.mat 文件一行代码就能读：

data = load('bearing_outer_fault.mat'); fs = data.fs; % 采样率, 单位 Hz acc = data.vibration; % 加速度信号, 单位 g 或 m/s^2 t = (0:length(acc)-1) / fs; % 时间轴, 用于时长核对与绘图

CSV 用 readtable 读取，注意第一行是否带表头，时间列是否被自动识别成 datetime，识别错误就按数值矩阵读。数据进来之后先校验采样率再开算：把 fs 与传感器配置值比对，同时用 length(acc) / fs 核对记录时长。低于 0.5 秒的数据做包络谱分析往往只有两三次冲击周期，谱线糊成一团，先补数据再做诊断。

3.2 去直流与带通滤波

原始加速度信号至少有两个必须去掉的成分：直流分量和高频噪声。直流来自传感器偏置，会让 FFT 零频出现巨大谱线，掩盖低频特征；高频噪声抬高包络谱底噪，降低谱峰信噪比。处理顺序固定为三步：去均值、带通滤波、取包络：

% 去直流 x = acc - mean(acc); % 带通滤波: 4 阶 Butterworth, 通带 2000 - 10000 Hz % 采样率必须高于 2*fc_high, 即 fs > 20000 Hz; 数据自带 fs 则优先用数据值 fc_low = 2000; fc_high = 10000; [b, a] = butter(4, [fc_low fc_high] / (fs / 2), 'bandpass'); x_filtered = filtfilt(b, a, x);

butter 把上下边频除以奈奎斯特频率 fs/2 归一化，归一化结果必须在 0 到 1 之间，所以 fc_high 超过 fs/2 时 MATLAB 直接报错。四阶是工程最常见的折中，阶数太高非线性相位畸变大，太低阻带衰减不足。filtfilt 相对 filter 的好处是零相位，正反各滤一次后谱峰位置不被相移拉偏，代价是数据两端有一点边缘效应，处理短数据时首尾各舍掉几十个点。

3.3 滤波参数设置与常见失败表现

带通上下边频没有万能值，应该贴近轴承固有共振频带。经验规则如下：

不知道共振频带时，我一般会先对原始信号做一次全频段 FFT，找到能量集中的高频峰群，把带通设置在该峰群的半功率带宽附近。这一步是包络分析里最容易被跳过却又最关键的地方，直接决定后续谱峰信噪比。不滤波直接做 Hilbert 变换是最典型的错误，得到的包络谱低频段被轴频及其谐波占满，故障频率即使存在也被淹没。

4. 包络谱分析：Hilbert 变换、谱峰搜索与故障判据

4.1 基于 Hilbert 变换的包络谱计算

故障冲击激起轴承座固有共振，表现为高频衰减振荡，故障频率本身以低频调制形式藏在幅值里。直接对原始信号做 FFT，高频共振能量把故障频率完全压住，谱图上看不到有效峰值；对包络做 FFT，故障频率作为幅值变化频率就凸显出来。Hilbert 变换把实信号变成解析信号，取模得到包络：

% 对滤波信号取包络, 再计算包络谱 env = abs(hilbert(x_filtered)); % 包络信号, 反映冲击幅值变化 N = length(env); win = hann(N); % 汉宁窗抑制频谱泄漏 spec = abs(fft(env .* win)); % 对包络做 FFT, 得到包络谱 freq = (0:N-1) / N * fs; % 频率轴 % 实数频谱对称, 只保留单边谱 spec = spec(1:floor(N/2)+1); freq = freq(1:floor(N/2)+1);

hann 窗是包络谱默认选择，旁瓣衰减快，适合特征频率与邻近噪声峰需要区分的场景。不加窗且数据长度不是故障频率周期整数倍时，谱峰能量泄漏到相邻频点，幅值偏低容易漏判。频率分辨率 Δf = fs / N，等价于记录时长的倒数，所以 1 秒数据给 1 Hz 分辨率，0.5 秒数据只有 2 Hz。BPFO 的 1% 容差约 1 Hz，采集前先按这个公式确认需要的记录长度。

4.2 特征频率处的峰值搜索与故障判据

人工看谱图可以，批量诊断必须把「看谱」变成「搜峰」。搜索逻辑：在特征频率附近开容差窗口，找窗口内最大幅值，与全谱底噪均值对比：

target_freq = bpfo; tol = 0.015 * target_freq; % ±1.5% 容差 idx = find(freq >= target_freq - tol & ... freq <= target_freq + tol); if ~isempty(idx) [peak_amp, peak_pos] = max(spec(idx)); peak_freq = freq(idx(peak_pos)); % 窗口内峰值对应的实际频率 mask = true(size(freq)); mask(idx) = false; % 从底噪统计中剔除目标窗口 noise_floor = mean(spec(mask)); fprintf('峰值 %.2f Hz, 幅值 %.3f, 信噪比 %.1f dB\n', ... peak_freq, peak_amp, 20*log10(peak_amp / noise_floor)); end

信噪比阈值经验值是 6 dB，即峰值幅值达到底噪两倍以上才判有效，低于这个值只能记疑似，需要结合时域峭度再判断。单个峰值不够，真实故障冲击是周期窄脉冲，谱上同时出现二倍频、三倍频。判定逻辑写成基频信噪比超 6 dB 且二倍频或三倍频存在才判故障，只出现基频峰很可能是随机冲击。

4.3 内外圈故障在谱图上的区分

外圈故障冲击幅值基本恒定，包络谱在 BPFO 整数倍处出现等间距谱峰，谱线干净。内圈故障缺陷随轴旋转，经过承载区时冲击被转频调制，BPFI 谐波两侧出现间隔为 fr 的边带。区分技巧是检查 BPFI 谐波两侧边带：在 n*BPFI ± fr 范围内搜索，边带幅值达到主峰三分之一以上基本可确认内圈故障。滚动体故障叠加自转影响，谱峰出现 BSF ± FTF 边带，且滚动体滑移让谱峰比外圈故障宽得多，呈「馒头状」而非尖峰状。

提示：包络谱只能证明存在周期性冲击，故障类型必须结合谐波阶次和边带结构判断，只盯基频最容易误判。

5. 时域特征提取与批量诊断：从单条谱线到设备健康度判断

5.1 峭度、均方根值与峰值因子计算

包络谱覆盖的是故障已发展到可辨识冲击的阶段，更早的微弱磨损要靠统计特征筛查。三个最常用特征：RMS 反映能量水平，峭度对早期冲击敏感，峰值因子把冲击高度与能量水平归一化：

function [kurt, rms, cf] = timeDomainFeatures(x) % 输入 x 为加速度信号(可滤波后使用) x = x - mean(x); rms = sqrt(mean(x.^2)); % 有效值, 单位与信号一致 sd = std(x); kurt = mean(x.^4) / sd^4; % 峭度, 健康轴承约等于 3 cf = max(abs(x)) / rms; % 峰值因子 end

峭度的物理意义是信号概率分布的尾部厚度。健康振动近似高斯分布，峭度约 3；故障初期冲击让波形尾部变厚，峭度迅速升到 5 以上。这个特征不依赖轴承几何参数，缺乏型号数据时也能做第一层筛查。峰值因子对随机噪声敏感，建议滤波后计算，否则高频噪声撑大 RMS、降低区分度。

5.2 特征-健康状态映射与阈值

三个特征组合覆盖不同健康阶段：

阈值必须随工况调整。转速波动大、载荷变动频繁的设备，RMS 本身就在波动，单看绝对值没有意义，要和同转速同载荷的历史数据对比。润滑状态与温度同样影响振动幅值，现场项目的阈值是标定出来的，不是抄教科书的。批量部署第一步永远是收集正常工况数据，先算出基线特征分布，再定报警线。

5.3 批量诊断脚本与阈值联动

现场诊断通常是几十个文件一起处理。批量循环如下：

files = dir(fullfile('data', '*.mat')); % 同一型号轴承的特征频率只需算一次 [bpfo, ~, ~, ~] = bearingFaultFreq(9, 7.94e-3, 39.04e-3, 0, 1797/60); results = table(); for i = 1:length(files) d = load(fullfile(files(i).folder, files(i).name)); x = d.vibration - mean(d.vibration); % 1. 时域特征 [kurt, rms, cf] = timeDomainFeatures(x); % 2. 带通滤波 + Hilbert 包络谱 [b, a] = butter(4, [2000 10000] / (d.fs / 2), 'bandpass'); env = abs(hilbert(filtfilt(b, a, x))); spec = abs(fft(env .* hann(length(env)))); freq = (0:floor(length(spec)/2)) / length(spec) * d.fs; spec = spec(1:floor(length(spec)/2)+1); % 3. BPFO 窗口内峰值与底噪对比 tol = 0.015 * bpfo; idx = find(freq >= bpfo - tol & freq <= bpfo + tol); peak = max(spec(idx)); snr = peak / mean([spec(1:idx(1)-1), spec(idx(end)+1:end)]); % 4. 峭度与谱峰联合判定 if kurt > 8 || (kurt > 5 && snr > 3) status = '严重故障'; elseif kurt > 4 && snr > 1.5 status = '疑似故障'; else status = '正常'; end results = [results; table({files(i).name}, kurt, rms, cf, peak, {status})]; end

滤波系数和窗函数在循环内重复计算，文件上百个时可以挪到循环外。判定顺序先看峭度再看包络谱信噪比，两者都在临界区时排除外部敲击这类瞬时冲击：瞬时冲击让峭度单次跳变，但包络谱特征频率处没有稳定峰值。两条线索合起来判断，误报率明显低于只看任一指标。

6. 用仿真信号验证 MATLAB 代码链路的正确性

没有现成故障数据时，仿真验证是排查代码逻辑错误的第一步。构造一个外圈故障仿真信号：以 5 kHz 共振频率的衰减振荡模拟轴承冲击响应，按 BPFO 周期触发，叠加白噪声：

fs = 48000; t = (0:fs*2-1)' / fs; % 2 秒数据 bpfo = 107.4; % 由 bearingFaultFreq 计算 zeta = 500; fn = 5000; T = 1 / bpfo; x = zeros(size(t)); for k = 0:floor(max(t) / T) % 每个冲击周期构造衰减振荡 idx = (t >= k*T) & (t < k*T + 0.02); % 冲击后 20 ms 响应窗 tt = t(idx) - k*T; x(idx) = x(idx) + exp(-zeta*tt) .* sin(2*pi*fn*tt); end x = x + 0.05 * randn(size(t)); % 白噪声模拟测量噪声

把 x 当作实测数据跑第 3、4 章的预处理与包络谱分析，验证标准有三条：107.4 Hz 处出现显著峰值；峰值与设定值偏差小于频率分辨率；二倍频位置存在次峰。仿真信号都检不出故障时，问题几乎一定出在滤波器通带或窗函数长度上。再把 bpfo 改小三分之一重跑，峰值跟随移动才能排除假峰。真实项目里用同样方法生成内圈、滚动体故障数据，写批量脚本前先验证阈值逻辑，避免拿现场数据试错。

提示：仿真验证的是代码链路没有逻辑错误，不代表算法对真实故障同样鲁棒。实际数据里的滑移、变载荷和温度漂移，必须靠现场标定覆盖。

本文还有配套的精品资源，点击获取