企业官网建设流程全解析

简介：压缩包提供了一套基于凯斯西储大学数据的轴承故障特征频率MATLAB计算源码，面向机械工程、设备维护及信号处理领域学习者与工程技术人员，旨在解决轴承故障诊断中特征频率计算与分析问题。整个压缩包共3个文件，含1个.m源代码和2张说明示意图（分别对应流程讲解与计算结果截图），大小仅91KB，轻量便携。源代码覆盖振动信号读取、数字滤波、傅里叶变换、频谱绘制及特征频率识别全流程，同时附有步骤分析基础讲解，详细说明基本旋转频率、滚动体通过频率和径向跳动频率等特征量如何由轴承结构参数推出，并区分内圈、外圈、滚动体故障在频谱上的倍数表现，为定位故障类型提供了清晰依据。已有4130人学习浏览，适合设备故障诊断课程设计、预防性维修研究及初学者快速上手，整体短小精悍，属于能直接运行并便于扩展的实用工具。

1. 项目概述

1.1 核心需求解析

做设备故障诊断的朋友，尤其是刚接触滚动轴承这一块的，大概率都绕不开一个经典数据集——凯斯西储大学（CWRU）轴承数据中心公开的振动信号数据。这套数据从上世纪90年代开始被广泛使用，几乎成了轴承故障诊断领域的“标准练习册”。但数据是公开了，真正用起来却有一道绕不过去的坎：拿到原始振动信号之后，怎么算出故障特征频率（BPFO、BPFI、BSF、FTF），并且把这些频率和频谱图对应起来，才是从“看数据”到“做诊断”的关键一步。

这个项目标题是“凯斯西储大学轴承故障特征频率MATLAB源代码.zip”，说白了就是一个帮你把CWRU轴承数据的特征频率计算从手算公式变成自动化代码的实用工具包。它解决的核心痛点非常明确：一是CWRU数据集的格式、采样率、转速等信息散落在各个文档里，整理起来费时；二是故障特征频率的计算必须严格依赖轴承几何参数，公式记错一个符号、单位搞错一个量级，后面的诊断结论就全歪了；三是即便算出了理论频率，如何跟实测频谱的峰值对齐，很多新手会卡在这一步。

这套源码适合谁来参考？我认为有三类人。第一类是刚入门的机械故障诊断方向的研究生，需要快速吃透CWRU数据并完成第一个实验；第二类是工业现场做设备状态监测的工程师，想用现成代码快速验证自己的诊断逻辑；第三类是上《机械故障诊断》《信号处理》课程的学生，拿来做课程设计或毕业设计的基础框架。这篇文章我会从轴承故障特征频率的物理原理讲起，再到CWRU数据集的实际参数梳理，最后完整拆解这套MATLAB源码的结构、计算逻辑和落地验证方法，把我在实际使用中踩过的坑一并交代清楚。

1.2 为什么选择MATLAB做这件事

之前有朋友问我，这种特征频率计算用Python写不行吗？当然可以，而且Python在机器学习后续处理上还有优势。但MATLAB在振动信号处理领域依然有不可替代的地位，至少有三点让我坚持用MATLAB来搞这套代码。

第一，MATLAB的信号处理工具箱非常成熟，pwelch、fft、envelope这些函数是现成的，而且对采样率、频率分辨率的处理逻辑封装得很好，出图质量高，用于论文和报告几乎不用二次加工。第二，MATLAB的交互式体验适合“探索性分析”。做故障诊断时经常要反复调整频带、观察谐波、对比不同故障类型的频谱特征，MATLAB的工作区变量管理和绘图窗口联动非常顺手，改参数立刻能看到结果。第三，CWRU数据集在学术圈的传播时间很早，大量经典论文的方法验证都是基于MATLAB完成的，很多对比实验的代码片段也是MATLAB写的，用MATLAB复现文献结果最省力。

当然，选择MATLAB也有代价——正版授权不便宜，而且代码对版本有一定依赖。但这套特征频率计算代码的核心逻辑只用到了基础语法和Signal Processing Toolbox，大家在R2016b及以上版本基本都能直接跑通，对版本要求并不苛刻。

2. 轴承故障特征频率的计算原理

2.1 四个特征频率的物理含义

要理解这套源码在算什么，先得搞清楚轴承故障特征频率到底是个什么东西。滚动轴承在工作时，内圈、外圈、滚动体和保持架之间相对运动，当某个部位出现局部损伤（比如剥落、裂纹、点蚀），滚动体经过损伤点的时候就会产生周期性的冲击。这个冲击的频率，由轴承的几何尺寸和转速唯一决定，就叫故障特征频率。

具体分四类：外圈故障频率（BPFO，Ball Pass Frequency of Outer Race），滚动体每经过外圈上的一个损伤点就产生一次冲击，这个频率主要跟滚动体数量、滚动体直径、节圆直径和接触角有关；内圈故障频率（BPFI，Ball Pass Frequency of Inner Race），类似的逻辑，只是损伤在内圈上，因为内圈随轴一起转，所以频率计算时会多一个转频的叠加项；滚动体故障频率（BSF，Ball Spin Frequency），滚动体上的损伤点在滚动体自转时反复接触内外圈滚道，频率是滚动体自转频率的两倍；保持架故障频率（FTF，Fundamental Train Frequency），保持架本身转速较慢，损伤后产生的冲击频率相对最低。

这四个频率的公式我在下面列出来。设D为节圆直径（滚动体中心所在圆的直径），d为滚动体直径，Z为滚动体数量，α为接触角，fr为转轴频率（Hz），则有：

• 外圈故障频率：BPFO = (Z × fr / 2) × (1 - d × cosα / D)
• 内圈故障频率：BPFI = (Z × fr / 2) × (1 + d × cosα / D)
• 滚动体故障频率：BSF = (D × fr / (2 × d)) × (1 - (d × cosα / D)²)
• 保持架故障频率：FTF = (fr / 2) × (1 - d × cosα / D)

注意，这里的fr不是转速的rpm，而是每秒多少转。如果给定转速是1797 rpm，那fr = 1797 / 60 = 29.95 Hz。这是新手最容易大意的地方，公式本身不复杂，但单位换算错了，频率直接差60倍，频谱上怎么找都找不到对应峰值。

2.2 CWRU数据集的轴承参数梳理

凯斯西储大学的这套公开数据用的测试轴承是SKF 6205-2RS JEM深沟球轴承（部分实验也用了NTN轴承做对比）。对于SKF 6205这个型号，业界通用的几何参数如下：

把接触角近似为0度处理，在工程上是合理的，因为深沟球轴承在纯径向载荷下接触角很小，对结果影响可以忽略。基于这些参数，在1797 rpm（约30 Hz）下计算出的理论特征频率大致是：BPFO约107.4 Hz，BPFI约162.2 Hz，BSF约70.6 Hz，FTF约11.9 Hz。

CWRU数据集的采样频率分两种：正常和故障数据大多是12 kHz采样，部分高速实验是48 kHz采样。驱动端轴承座采集的振动信号，是最常用的分析对象。还有一点必须提醒大家：数据文件名里面会标注故障直径（0.007、0.014、0.021英寸）、故障位置（内圈IR、外圈OR、滚动体B）、电机负载（0~3 hp）。做特征频率验证时，优先选择单一故障类型、单一负载工况的数据，不要一上来就混合分析，不然频谱上各种频率混叠在一起，新手很容易绕晕。

2.3 理论频率和实际频谱为什么有偏差

即便公式算得很准，实测频谱中故障特征频率的峰值也不会正好落在理论值上，通常会有1%到3%的偏移。原因有几个：一是轴承实际运行中温度升高，导致轴承游隙变化，滚动体和滚道的接触几何随之改变；二是轴承存在打滑现象，尤其是轻载工况下滚动体在滚道上的滑动更明显；三是转速本身有波动，CWRU数据集的转速是设定值，但采集过程中电机负载和电网波动会让实际转频有小幅漂移。

这就引出一个实操上的重要观念：算出来的特征频率是“指导线”，不是“精确刻度”。在看频谱时，应该在理论频率附近取一个窄带搜索（比如±3%），找到该范围内幅值最大的峰值作为实际故障特征频率，并观察它是否出现二倍频、三倍频等高次谐波。这套源码里也内置了这个思路，后面我会展开讲。

3. 源码整体设计与模块拆解

3.1 文件结构与运行流程

这个zip包解压之后，核心文件主要包含三类：主脚本、函数文件和参数配置文件。我建议不要把全部代码都堆在一个文件里，哪怕只是几行计算，分模块的好处是后续替换数据、调整参数都很灵活。整理后的目录结构大致如下：

CWRU_Bearing_Fault_Frequency/ ├── main.m % 主脚本，一键运行 ├── config.m % 参数配置文件，轴承型号和工况参数 ├── load_cwru_data.m % 数据加载函数，自动读取.mat文件 ├── compute_fault_freq.m % 六频计算函数 ├── plot_spectrum.m % 频谱绘制函数 └── data/ % 存放CWRU原始数据文件

主脚本的运行流程很直接：先运行config设定参数，再加载数据，然后计算特征频率，最后绘制时域波形和频谱图，并在频谱图上把理论特征频率的位置用竖线标出来。整个流程不需要人工介入，跑完就能得到一张带频率标注的诊断图。

这种设计的好处在于“计算逻辑”和“数据处理”解耦。如果换了其他型号的轴承，只要改config.m里的几何参数就行；如果要处理自己的实测数据，只需要替换load_cwru_data.m的读取逻辑，后面的计算和绘图流程完全复用。很多初学MATLAB的朋友容易把代码写成一个大脚本，从数据读取到出图全揉在一起，改一处就要动全身。我之前也这么干过，后来维护起来真的头大，强烈建议从一开始就养成模块化的习惯。

3.2 参数配置模块：为什么单独拎出来

config.m看起来只是几个变量赋值，但我坚持把它单独成一个文件，原因是CWRU数据集的“坑”就在参数上。转速、负载、采样率、故障位置、故障直径、轴承型号这些参数如果硬编码在主脚本里，换一个数据文件就要去代码里翻找修改，效率低还容易漏改。

config.m里的核心参数至少包括：

%% 轴承几何参数（SKF 6205-2RS JEM） bearing.Z = 9; % 滚动体数量 bearing.d = 7.94; % 滚动体直径，单位mm bearing.D = 39.04; % 节圆直径，单位mm bearing.alpha = 0; % 接触角，单位度 %% 工况与数据参数 data_path = 'data/'; % 数据文件夹路径 file_name = 'IR007_0.mat'; % 实际文件名 fs = 12000; % 采样频率，单位Hz shaft_speed_rpm = 1797; % 电机转速，单位rpm %% 分析参数 freq_band = [0 600]; % 频谱分析频带，单位Hz line_width = 2; % 标注线宽度

特别注意CWRU数据文件的命名规则。比如IR007_0.mat，IR代表内圈故障（Inner Race），007代表故障直径0.007英寸，最后的0代表负载0 hp。如果是外圈故障，可能是OR007@6_0.mat，这个@6表示故障点在6点钟方向。CWRU外圈故障数据有3点钟、6点钟、12点钟三个加载方向，对特征频率本身没有影响，因为频带的分布主要是由损伤引起的冲击周期决定的，但不同方向的故障信号在幅值上会不同，分析时要注意。

这类命名规则如果第一次接触，很容易把文件编号当成转速或者别的参数，直接导致后面算错了频率还在频谱里使劲找峰。把数据文件的完整路径和命名规则写在config的注释里，是我在实际项目中增加的额外保险。

3.3 数据加载函数：处理CWRU原始.mat文件

CWRU数据文件虽然是.mat格式，但它的内部变量结构并不统一。有的版本存的是DE（驱动端加速度信号）、FE（风扇端加速度信号）、BA（基座加速度信号）这几个数组；有的版本还附带转速时间序列RPM。写加载函数的时候不能写死只取某一个变量，最好做一个自动查找的机制。

一个稳妥的加载函数写法参考：

function [data, fs, rpm] = load_cwru_data(file_path) % 加载CWRU轴承数据集.mat文件 % 返回：振动信号data，采样率fs，转速rpm S = load(file_path); fields = fieldnames(S); % 自动识别驱动端加速度信号 if isfield(S, 'DE') data = S.DE; elseif isfield(S, 'X100_DE_time') data = S.X100_DE_time; else % 取第一个数组类型字段 for i = 1:length(fields) if isnumeric(S.(fields{i})) && numel(S.(fields{i})) > 1000 data = S.(fields{i}); break; end end end % 采样率：CWRU公开数据中12k和48k两种 if length(data) > 100000 fs = 48000; else fs = 12000; end % 转速：可从文件名或RPM字段获取 if isfield(S, 'RPM') rpm = S.RPM(1); else rpm = []; end end

这里用isfield和fieldnames做自动识别，是处理这种“版本混杂”数据集的实用技巧。我自己在整理CWRU数据的过程中发现，不同渠道下载的数据集变量命名可能不一样，有的叫X100_DE_time，有的叫X200_DE_time，前缀对应的是不同的实验工况编号。所以不要完全依赖变量名，必要时直接遍历所有字段，找出那个“长度明显是振动信号”的数值数组。

另一个重点：CWRU的12k数据和48k数据的文件大小差别很大，判断采样率可以用信号长度做近似，但这不严谨。如果文件名或者文档中有明确标注，优先用文件信息确定采样率。更稳妥的方式是准备一个文件名-采样率对照表，把它维护在config.m里，一劳永逸。

3.4 特征频率计算函数：公式的代码实现

计算函数是这套源码的核心，代码本身不长，但每一步都要对应到公式，注释一定要写清楚，方便后续检查。

function [freq, labels] = compute_fault_freq(bearing, fr) % 计算滚动轴承四个故障特征频率 % 输入：bearing结构体（Z, d, D, alpha），fr转轴频率（Hz） % 输出：freq为四个频率值，labels为对应的故障类型名称 Z = bearing.Z; d = bearing.d; D = bearing.D; alpha = bearing.alpha * pi / 180; % 角度转弧度 cos_a = cos(alpha); % 外圈故障频率 BPFO BPFO = Z * fr / 2 * (1 - d * cos_a / D); % 内圈故障频率 BPFI BPFI = Z * fr / 2 * (1 + d * cos_a / D); % 滚动体故障频率 BSF BSF = D * fr / (2 * d) * (1 - (d * cos_a / D)^2); % 保持架故障频率 FTF FTF = fr / 2 * (1 - d * cos_a / D); freq = [BPFO, BPFI, BSF, FTF]; labels = {'BPFO', 'BPFI', 'BSF', 'FTF'}; end

有两点值得注意。第一，接触角我做了度转弧度的处理，因为config里人的输入习惯是角度制，但MATLAB的cos函数要求弧度制。这种“输入友好、内部严苛”的做法能减少低级错误。第二，BSF和FTF在实际信号中往往没有BPFO、BPFI那么明显，因为滚动体和保持架的故障冲击信号在传递路径上衰减更多，频谱上经常只出现幅值很低的峰。

如果你用的是NTN轴承或者别的型号，轴承的Z、d、D参数要重新查手册，不要沿用6205的参数。因为CWRU数据集里虽然主力是SKF 6205，但也有用NTN轴承做的对比实验，两种轴承的几何参数是不同的，直接用6205的参数算NTN数据的特征频率，结果就会对不上。这是我在对比实验时踩过的坑，当时频谱上找不到理论频率，一度怀疑数据有问题，后来查了文档才发现是轴承型号没有切换。

3.5 频谱可视化：把理论频率叠加到实测频谱上

计算出了四个特征频率只是第一步，真正让结论“立得住”的是把它画到频谱图上。plot_spectrum.m函数负责这部分工作：

function plot_spectrum(data, fs, bearing_freqs, labels, freq_band) % 绘制振动信号频谱，并标注理论故障特征频率 N = length(data); if mod(N, 2) ~= 0 data = data(1:end-1); end % 去除直流分量 data = data - mean(data); % 计算单边幅值谱 Y = fft(data); P2 = abs(Y / N); P1 = P2(1:N/2+1); P1(2:end-1) = 2 * P1(2:end-1); f = fs * (0:(N/2)) / N; % 限定频带 idx = f >= freq_band(1) & f <= freq_band(2); f_band = f(idx); P1_band = P1(idx); figure('Color', 'white', 'Position', [100 100 900 500]); plot(f_band, P1_band); xlabel('Frequency (Hz)'); ylabel('Amplitude (m/s^2)'); title('Vibration Spectrum with Fault Characteristic Frequencies'); grid on; hold on; % 标注理论特征频率 colors = {'r', 'm', 'g', 'b'}; for i = 1:length(bearing_freqs) if bearing_freqs(i) >= freq_band(1) && bearing_freqs(i) <= freq_band(2) xline(bearing_freqs(i), '--', colors{mod(i-1, length(colors))+1}, ... 'LineWidth', 1.5, 'Label', labels{i}, 'FontSize', 10); end end hold off; end

画频谱的两个细节影响很大。第一是FFT前一定要去直流分量，否则频谱在0 Hz处会有一个很高的峰值，整个纵轴被压扁，有效频段的幅值特征都看不清了。第二是用xline而不是手工写plot([f f], [ylim_min ylim_max])，因为xline的标注文字会自动避开坐标轴，并且支持标签显示，代码也简洁很多。

绘制频谱时，数据段的长度直接影响频率分辨率。频率分辨率Δf = fs / N，采样率12 kHz，如果取48000个点（即4秒数据），Δf = 0.25 Hz，这个分辨率对识别BPFO、BPFI足够了。如果数据点太少（比如1024个点），Δf = 11.7 Hz，这会把BPFO和BPFI之间的间隔（约55 Hz）给糊掉，频谱上完全看不出两个峰。所以建议做FFT时至少取2万个点以上，保证频率分辨率在1 Hz以内。

4. 实操过程与核心环节实现

4.1 零基础上手：环境准备与运行

在开始跑代码之前，先把环境准备好。我用的环境是MATLAB R2021a，操作系统是Windows 10。理论上R2016b以上版本都能运行，因为用到的xline函数是R2016a才引入的，fieldnames、fft这些基础的不能再基础了。如果你的版本更老，可以用plot函数手动代替xline。

下载CWRU数据集的时候注意找对下载入口。凯斯西储大学轴承数据中心官网提供了所有实验数据，按文件名索引下载。有些第三方平台也做了镜像整理，但一定要注意数据文件是否完整、采样率标注是否正确。我更推荐自己写一个下载脚本按需拉取，避免一次性下十几个GB的数据。没有脚本的话，用浏览器逐条下载也行，只是效率低一些。

解压zip包之后，把CWRU的.mat数据文件放进data/文件夹，然后在MATLAB里打开main.m，选择当前工作目录为项目根目录，直接点击“运行”即可。主脚本在执行完计算和绘图之后，会在命令行输出一个汇总表格，格式大致如下：

故障类型 理论频率(Hz) 频谱实测峰值(Hz) 误差(%) BPFO 107.36 106.80 0.52 BPFI 162.18 162.50 0.20 BSF 70.58 69.20 1.95 FTF 11.93 N/A N/A

输出的实测峰值是通过在理论频率附近搜索局部最大值得到的，搜索带宽默认设为2%。这个表格能快速评估计算结果的可靠性。

4.2 频带选择与窗函数：频谱质量的细节

在实际处理中，频谱质量的好坏跟两个因素直接相关：参与FFT的数据段选择和窗函数。

CWRU数据是平稳旋转机械的振动信号，理论上用矩形窗直接截取就能得到不错的频谱，但这会引入频谱泄漏。你可以这样理解：矩形窗相当于强制把一个非整数周期的信号段当作一个周期来处理，截断处的不连续在频域上表现为旁瓣泄漏，把原本干净的频峰会“糊”出一圈毛刺。更实际的做法是用汉宁窗（Hann）对数据段进行加权，主瓣会变宽一些，但旁瓣大幅降低，峰值定位更稳定。

代码里加窗的操作非常容易：

window = hann(N, 'periodic'); data_windowed = data .* window; Y = fft(data_windowed);

关于频带选择，CWRU在12 kHz采样下，分析到600 Hz已经能覆盖前几阶故障特征频率和谐波。但如果要做包络分析（envelope analysis），需要把分析频带扩到更高范围，因为故障冲击会激起轴承结构的高频共振，解调之后才会在低频段看到特征频率。这个包络分析是另一个话题了，这里不展开。总之，初版代码用0到600 Hz这个频带，对0.007英寸故障这类早期故障已经足够看到明显的BPFO峰值。

4.3 完整代码示例：主脚本串联全流程

把所有模块串起来的主脚本main.m完整参考如下：

%% CWRU轴承故障特征频率分析主脚本 % 作者：公众号/博客同名 % 功能：一键计算并可视化CWRU数据集的轴承故障特征频率 % 适用：MATLAB R2016b及以上版本 clear; clc; close all; %% 1. 加载参数配置 run('config.m'); %% 2. 自动加载数据 full_path = fullfile(data_path, file_name); [data, fs, rpm_from_data] = load_cwru_data(full_path); % 如果数据文件本身有转速信息，优先使用数据中的转速 if ~isempty(rpm_from_data) shaft_speed_rpm = rpm_from_data; end fr = shaft_speed_rpm / 60; % 转轴频率，单位Hz fprintf('数据文件: %s\n', file_name); fprintf('采样频率: %d Hz\n', fs); fprintf('转轴频率: %.2f Hz\n', fr); %% 3. 计算故障特征频率 [freqs, labels] = compute_fault_freq(bearing, fr); for i = 1:length(freqs) fprintf('%s: %.2f Hz\n', labels{i}, freqs(i)); end %% 4. 绘制频谱并标注特征频率 plot_spectrum(data, fs, freqs, labels, freq_band); %% 5. 搜索理论频率附近的实测峰值 search_band_ratio = 0.02; % 搜索带宽：理论频率的±2% N = length(data); if mod(N, 2) ~= 0 data = data(1:end-1); end data = data - mean(data); Y = abs(fft(data) / N); P1 = Y(1:N/2+1); P1(2:end-1) = 2 * P1(2:end-1); f_axis = fs * (0:(N/2)) / N; fprintf('\n理论频率附近的实测峰值搜索：\n'); fprintf('%-8s %-12s %-14s %-8s\n', '故障类型', '理论频率(Hz)', '实测峰值(Hz)', '误差(%)'); for i = 1:length(freqs) target = freqs(i); lo = target * (1 - search_band_ratio); hi = target * (1 + search_band_ratio); idx_range = f_axis >= lo & f_axis <= hi; if sum(idx_range) > 0 [peak_val, peak_idx_local] = max(P1(idx_range)); idx_full = find(idx_range); peak_freq = f_axis(idx_full(peak_idx_local)); err = abs(peak_freq - target) / target * 100; fprintf('%-8s %-12.2f %-14.2f %-8.2f\n', labels{i}, target, peak_freq, err); else fprintf('%-8s %-12.2f %-14s %-8s\n', labels{i}, target, 'N/A', 'N/A'); end end

这个主脚本把从配置、加载、计算到验证的完整流程走了一遍，运行完成后你会在当前文件夹看到多张频谱图。如果是第一次跑，建议先用0.007英寸内圈故障的数据（比如IR007_0.mat）做验证，因为这类故障的特征频率峰值非常明显，最容易确认代码是否正确。

4.4 参数调整实战：换数据文件、换轴承型号

我拿一个实际操作场景来演示。假设我想分析外圈故障0.014英寸、负载2 hp的数据（文件名可能是OR014@6_2.mat），只需要修改config.m里的两项：

file_name = 'OR014@6_2.mat'; shaft_speed_rpm = 1750; % 2 hp负载对应的转速

负载和转速的对应关系，CWRU数据集文档中有明确说明：0 hp对应1797 rpm，1 hp对应1772 rpm，2 hp对应1750 rpm，3 hp对应1730 rpm。这个对应关系在不同故障类型下略有差异，但大致如此。如果没注意到这个规律，会发现明明用同一个数据集，但是频谱上的特征频率整体偏移了一段，因为转速变了，所有特征频率都跟着变了。

这个问题的背后是量化分析方法。我在项目里会额外增加一个可视化：把转频fr做了归一化处理，四个特征频率都除以fr，得到一组无量纲的“频率倍数”。因为BPFO/BFI/BSF/FTF都正比于fr，归一化之后理论值就固定了，不随转速变化。这在分析变转速工况时非常有用。源码里也可以加这样一个功能，方便在多个转速下快速验证。

5. 常见问题与排查技巧实录

5.1 频谱上找不到理论频率的峰值

这是最多人遇到的问题，也是特征频率分析里最让人沮丧的情景。我先按照自己的排查经验给出一张速查表：

基于实测经验，最常犯的错误是轴承参数没有跟着数据文件切换。比如分析NTN轴承数据时还在用SKF 6205的参数。还有一种情况是误用了风扇端（FE）的数据，风扇端轴承型号和驱动端不一样，特征频率也不同。如果CWRU文件名里没有明确标注，可以从文档确认当前文件对应哪个测试点。

5.2 外圈故障的“@6”方向标记是什么

CWRU外圈故障数据中，文件名里会带@3、@6、@12的后缀，代表外部损伤的加载方向。这个方向对特征频率的理论计算没有影响，但它会影响实测信号的幅值分布，因为故障点相对传感器和载荷区的位置不同，冲击信号的传递路径衰减不一样。

我遇到过这样一个情况：用OR007@6_0.mat分析时，BPFO的理论频率在频谱上很明显，但再用OR007@3_0.mat时，同一频率处峰值大幅降低。这不是算法错了，而是故障点的位置在载荷区外，冲击能量没有充分传递到传感器。这种情况下建议不要只看单一方向的频谱，可以综合分析多个方向的数据，或者改用包络谱来放大故障冲击特征。把这一点写进代码注释里，能省去很多不必要的排查时间。

5.3 频率分辨率不足导致双峰无法区分

BPFO和BPFI在很多时候只差几十赫兹，理论上好区分，但如果FFT点数太少，频率分辨率不够，两个峰就糊在一起了。

举个例子，若只取4096个点，12 kHz采样下Δf = 2.93 Hz。虽然BPFO和BPFI相差约55 Hz，理论上还是能分开，但如果你分析的是滚动体故障（BSF）和保持架故障（FTF）这种低频且幅值较弱的特征，由于频谱泄漏和噪声干扰，2.93 Hz的分辨率会让峰值定位误差较大。所以我建议默认取2万到5万个点，如果想做谐波分析或者更精细的诊断，甚至可以取10万以上个点。数据量多不是问题，CWRU每个文件都有几十万个点，完全够用。此外，如果两个频率峰距离非常近，可以尝试用零填充（zero padding）提高频谱的插值精度，但要注意，零填充不会提高真实频率分辨率，它只是让峰形更平滑。

5.4 命令行报错“未定义函数或变量”

这类报错通常有三个原因。第一个是当前工作目录不对，MATLAB找不到函数文件。解决办法是在运行main.m之前用cd命令切到项目根目录，或者在MATLAB编辑器中右键main.m选择“运行文件”，会自动将文件所在目录加入搜索路径。第二个是函数名拼写错误，比如compute_fault_freq写成了compute_fault_frequencies。第三个是config.m中的变量没有加载到工作区，因为你用的是run('config.m')，这种方式加载的变量会覆盖工作区同名变量，但如果config里故意写了clear语句，其他变量会被清掉。建议config.m里不要写clear，保持纯参数赋值。

我能给的最大建议是：任何报错，先看错误信息里提到的行号，跳过去对照代码检查，不要一上来就怀疑数据或算法。80%的错误都是路径、变量名、参数单位这些“低级”问题，但排查起来最耗时间，模块化代码配合清晰注释能帮你大幅减少这类问题。

6. 扩展思路：从特征频率计算到故障智能诊断

6.1 包络分析：让早期微弱故障“现形”

直接对原始振动信号做FFT，早期故障的特征频率峰值往往被淹没在噪声和振动能量中。原因是故障冲击本身能量不大，且传递路径上被结构衰减。包络分析（Envelope Analysis，也叫解调分析）的思路是：先对信号做带通滤波，滤出包含故障冲击共振的高频带，然后取包络（希尔伯特变换或检波），再对包络信号做FFT，这个时候低频段的故障特征频率就非常明显了。

MATLAB里做包络分析非常方便。可以用bandpass函数先滤出共振频带，再用hilbert提取包络：

% 带通滤波：中心频率和带宽需要根据频谱特征选择 data_filt = bandpass(data, [2000 5000], fs); % 包络提取 envelope_sig = abs(hilbert(data_filt)); % 包络谱 N = length(envelope_sig); Y = abs(fft(envelope_sig - mean(envelope_sig)) / N);

在上面这段代码里，带通滤波的频带选择是关键。选低了可能滤掉共振频段的能量，包络信号平坦无特征；选高了可能引入其他噪声。CWRU数据中，0.007英寸早期故障的共振频带大致在2 kHz到5 kHz，但这只是一个经验区间，实际应结合信号频谱的“能量聚集区”来确定。还有一个实用技巧：可以先看原始信号的频谱，找到振幅明显抬高的频带，那就是共振区，带通滤波就选那里。

6.2 用特征频率自动化判定故障类型

特征频率计算出来之后，下一步的自然延伸是“自动判断轴承有没有故障、是什么故障”。一个接地气的做法是：在频谱中搜索四个特征频率附近的峰值幅值，如果某个特征频率处的幅值明显高于周围底噪，并且它有整数倍的谐波，就判定该位置存在故障。

这种阈值判断思路简单粗暴，但很有效。比如设定底噪水平为频带内振幅的均值加上3倍标准差，超过这个阈值的候选峰值才视为故障特征。如果BPFO处有峰且有BPFO×2、BPFO×3等高次谐波，就判定为外圈故障。这个逻辑用MATLAB实现起来是一个for循环加条件判断的事，而且可以直接接到本套源码后面。对于课程设计或毕业设计，这种“规则判断+可视化”的方式比直接套深度学习模型更能展示对底层原理的理解。

6.3 数据增强：从单点故障到多故障耦合

CWRU数据集本身是单点故障为主，但真实工业场景中轴承经常出现叠加故障（内圈外圈同时损伤），或者故障特征被齿轮箱等其他部件干扰。如果想把项目做得更深入，可以在现有代码基础上模拟多故障信号的叠加：把两个不同故障位置的理论脉冲序列叠加在一起，再添加噪声和转频谐波，就能构造出多故障耦合的仿真信号。这本质上是“信号级数据增强”，可以用来验证多故障诊断算法在CWRU原始数据之外的泛化能力。

但这里要提醒一点：仿真信号的故障机理相对理想，不能完全替代真实数据，所以实际项目里应该以CWRU原始数据验证为主，仿真数据作为辅助。我自己做项目时，用仿真信号验证算法逻辑的完备性，再用CWRU和实测数据去验证算法在真实噪声下的鲁棒性，二者互补，效果最好。

6.4 这套代码的局限性

说实话，这套特征频率计算代码解决的是“入门验证”和“常规分析”场景，它也有一些明显的边界。第一，它不支持转速变化工况，如果转速是波动的，直接FFT频谱会模糊，需要做阶比跟踪（order tracking）。第二，它没有内嵌自动特征提取和机器学习分类模块，故障类型的判定仍然需要人工确认。第三，CWRU数据本身是高信噪比的实验室数据，与工业现场的复杂环境有差距，用这套代码跑出来的结论不能直接等同于现场诊断结论。

如果要把这套代码往工程化方向推进，建议补上三个模块：一是基于包络谱的自动特征提取，二是多工况数据的批量处理和结果汇总，三是与数据库或报表系统对接，让每次分析的结果能自动归档。这几个方向在后续文章中我可以逐个展开讲。

7. 写在最后的工程建议

7.1 想跑通和跑好是两码事

我有一次帮朋友调代码，他拿到CWRU数据之后直接用默认参数跑，出了一张频谱图，图上四个理论频率标注线整整齐齐，但数据是内圈故障，频谱上最大的峰值却不在BPFI附近。排查了很久才发现，他把config里的轴承型号参数填的是另一个型号的，压根不是SKF 6205。从那之后，我自己用这套代码的习惯就固定成三步：先打印轴承参数和转速，确认“我用的参数是什么”；再打印理论频率值，确认“理论上该在哪里看到峰”；最后才跑到频谱图阶段。这三步各花几秒钟，但能省掉后面几十分钟的无效排查。

7.2 学会“读”频谱，而不是只看标注线

代码能把四个特征频率的竖线画在图上，但这不代表你就能直接下结论。看频谱要有层次感：先看整体能量分布，确认哪些频带能量高（是转频谐波还是共振频带）；再看故障特征频率的峰值幅值是否突出，并且观察它的高次谐波；最后对比不同工况、不同故障类型的数据，这样判断才靠谱。也就是说，这个代码给你的是“线索”，不是“答案”。

7.3 后续扩展方向简述

如果这篇入门代码你已经跑通了，后面有三条路可以选：第一条是从CWRU数据切换到真实工业信号，挑战会急剧上升，但收获也最大；第二条是结合时域指标（均方根值、峭度、峰值因子）做多域特征融合，这是传统诊断转向智能诊断的必由之路；第三条是把这套MATLAB脚本包装成一个简单的GUI工具或者函数库，方便团队内部其他人直接使用。我个人的实际体会是，哪怕只是把这几百行代码整理清楚、注释补全，对你的信号处理基本功都是很好的训练。轴承故障诊断这条路上没有捷径，但把特征频率这块地基打牢了，后面学包络分析、时频分析、深度学习诊断都会顺手很多。

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