企业官网建设流程全解析

简介：面向机械工程、设备维护与信号处理学习者，帮助解决轴承故障诊断中特征频率计算与分析的实际问题，尤其适合需要把理论公式落实到代码层面的初学者。内容基于凯斯西储大学公开轴承数据，提供MATLAB实现，涵盖振动信号读取、数字滤波、傅里叶变换、频谱绘制与故障特征频率识别等关键环节，并附有“轴承故障特征频率步骤分析基础讲解”，对BPFI、BPFO等典型特征频率以及滚动体损伤、内外圈故障对应的频率倍数做了图示说明，便于对照理解物理含义。压缩包共3个文件，包括1个MATLAB源代码文件与2张说明截图，整体仅91KB，结构精简，可快速上手而不占用额外空间。已有4130人学习浏览，适合希望结合理论与实操入门轴承故障诊断的读者，运行代码即可观察不同损伤部位对应的频谱特征。 前阵子一个做设备状态监测的朋友拿了一组轴承振动数据来问，说频谱里那个尖峰到底该用哪个频率去对齐。我第一反应就是让他先把凯斯西储大学（CWRU）那套经典数据集的故障特征频率算明白，因为这套数据几乎是轴承故障诊断入门和算法验证绕不开的标准样本。可惜很多人在这一步就卡住了，要么是把公式套错，要么是算出的频率和官方给的对不上，后面做包络谱、训练分类模型全跟着翻车。这次我把自己整理的一套MATLAB源代码完整拆开讲一遍，从原理公式到代码实现，再配合CWRU真实数据做验证，把我踩过的坑也一并列出来，希望能帮你少走点弯路。

1. 先搞清楚：这套代码到底解决什么问题

1.1 CWRU数据集里有什么，为什么绕不开它

CWRU轴承数据中心提供的滚动轴承故障数据，可以说是故障诊断领域最常用的公开数据集之一。整套数据来自一个电机-轴承试验台，分别采集了驱动端和风扇端的振动加速度信号，故障轴承是SKF 6205-2RS JEM深沟球轴承，采用电火花加工的方式在内外圈、滚动体上制造了单点故障，故障直径有0.007英寸、0.014英寸、0.021英寸等几种规格，还设置了0到3马力四档负载工况。

这套数据集之所以被高频使用，是因为它结构完整、文件命名规律清晰，每个.mat文件里都包含时域振动信号和时间序列，采样率有12kHz和48kHz两档。对做算法验证的人来说，你不需要自己搭建试验台，就能拿到带明确故障标签的真实信号，用来测试特征提取、频域分析、机器学习分类的效果。但有个前提是，你得先把不同故障类型对应的特征频率算准，否则后续所有频域分析都是在乱猜。

1.2 计算特征频率的价值，不只是一个数学作业

轴承故障特征频率的本质，是故障点与滚动体、滚道接触时产生的周期性冲击在频谱上对应的位置。比如外圈固定不动，故障点每次经过承载区就产生一次冲击，这个冲击的重复频率就是外圈故障特征频率。有了这个频率，我们在频谱图里就能快速定位故障类型，而不是靠肉眼去看时域波形里的毛刺。

我实际处理振动数据时，特征频率还有一个更重要的作用：它是带通滤波和包络谱分析的参考中心。直接看原始振动加速度的频谱，往往会被宽频噪声和结构共振淹没，通常需要对目标频段做带通滤波再做包络解调，才能清晰地找到故障特征。而带宽怎么选、中心频率设在哪，都要以特征频率的理论值做参考。换句话说，这是整个故障诊断链条的地基。

2. 特征频率计算：公式和官方倍率到底该信哪个

2.1 四类故障特征频率的物理含义

滚动轴承的特征频率通常分四种，分别是外圈故障特征频率BPFO、内圈故障特征频率BPFI、滚动体故障特征频率BSF和保持架故障特征频率FTF。

外圈故障特征频率BPFO是故障点通过承载区时产生的冲击频率，外圈固定在轴承座里不动，因此这个频率只和滚动体公转速度有关。内圈故障特征频率BPFI则要复杂一些，因为内圈随轴一起旋转，故障点本身在转动，冲击频率是滚动体通过故障点的频率，所以它等于滚动体数乘以转频和保持架转频之差。滚动体故障特征频率BSF比较特殊，滚动体同时和内外圈接触，每自转一圈会产生两次冲击，这个频率和滚动体自转速度直接相关。保持架故障特征频率FTF其实是滚动体围绕轴心公转的速度，也就是保持架的转速，前三个频率都建立在它的基础之上。

这四种频率都正比于轴的转频fr，所以工程上习惯用“倍率”来表示，也就是特征频率等于多少倍的转频。转频根据实际转速算出来以后，乘以对应倍率就能得到具体的频率值。

2.2 理论几何公式与CWRU官方倍率存在差异的原因

理论上，特征频率可以由轴承几何参数直接计算，标准公式如下：

• FTF = (fr / 2) × (1 - (Bd / Dp) × cos α)
• BPFO = N × FTF
• BPFI = N × (fr - FTF)
• BSF = (Dp / (2 × Bd)) × fr × (1 - (Bd / Dp × cos α)²)

其中N是滚动体数量，Bd是滚动体直径，Dp是节圆直径，α是接触角。以CWRU常用的6205轴承为例，滚动体数量N=9，滚珠直径Bd约为0.3126英寸，节圆直径Dp约为1.537英寸，接触角按0°处理。把数值代入公式，可以得到理论倍率大约是：

• FTF ≈ 0.3983 × fr
• BPFO ≈ 3.5847 × fr
• BPFI ≈ 5.4150 × fr
• BSF ≈ 2.3565 × fr

但CWRU官网给出的官方倍率却是BPFO=3.048、BPFI=4.857、BSF=1.994、FTF=0.398。FTF对得上，但其他三个差距不小。我在第一次对比时也一度怀疑自己公式写错了，后来确认这就是正常现象。

差异主要来自几个方面：一是实际工况下滚动体存在滑动，不是纯滚动，导致公转转速和理论值产生偏差；二是轴承在承受载荷时接触角、弹性变形会改变，几何参数不能完全按静态值计算；三是电火花加工出的故障位置有限，冲击形式也不是理想周期脉冲，用不同解调方法测到的频率会有细微差别。因此在做CWRU数据验证时，最稳妥的方式是直接用官方倍率，因为后续研究者对比报告、论文里普遍使用的也基本都是这套倍率；而几何公式更适合在自行设计试验台、手里只有轴承图纸时做预判，两者用途不同。

3. MATLAB实现：从参数输入到频率输出

3.1 程序框架设计思路

为了实现一个既能在CWRU数据上直接验证、又能在工程现场灵活使用的工具，我的代码没有搞成死板的一言堂，而是把计算模式做成可选参数。一种是调用CWRU官方倍率，直接适用于这个经典数据集；另一种是输入实际轴承几何参数，用理论公式计算，方便移植到其他型号轴承上。

输出端用结构体统一返回BPFO、BPFI、BSF、FTF四个频率值，同时附上各自的转频倍率。结构体的好处是字段名一目了然，后续做频谱标注时不用记索引位置。代码量控制在几十行，没有依赖额外的工具箱，纯MATLAB基础函数就能跑，这对很多人来说比较重要，因为换一台没装专业工具箱的电脑也不至于报错。

3.2 核心函数代码与逐段解析

下面这个faultFreq函数是整套代码的核心，我注释写得很详细，建议直接复制到MATLAB里运行一遍。

function freq = faultFreq(mode, params, frHz) % 计算滚动轴承故障特征频率 % 输入： % mode - 'cwru'：使用凯斯西储大学官方倍率，适用于CWRU数据集 % 'geo' ：根据轴承几何参数理论公式计算 % params - struct，geo模式下需要包含字段： % Bd(滚珠直径)、Dp(节圆直径)、Z(滚珠数量)、alpha(接触角，单位度) % frHz - 轴的转频，单位Hz，frHz = 转速RPM / 60 % 输出： % freq - struct，包含BPFO、BPFI、BSF、FTF及各倍率 if strcmpi(mode, 'cwru') % 6205-2RS JEM SKF 在CWRU数据集中的官方倍率 k_BPFO = 3.048; k_BPFI = 4.857; k_BSF = 1.994; k_FTF = 0.398; freq.BPFO = k_BPFO * frHz; freq.BPFI = k_BPFI * frHz; freq.BSF = k_BSF * frHz; freq.FTF = k_FTF * frHz; % 把倍率也存一份，后面画图标注会用到 freq.k = [k_BPFO, k_BPFI, k_BSF, k_FTF]; elseif strcmpi(mode, 'geo') Bd = params.Bd; Dp = params.Dp; Z = params.Z; alpha = params.alpha * pi / 180; % 角度转弧度 % 中间量：滚珠直径与节圆直径的比值，再乘接触角余弦 ratio = Bd / Dp * cos(alpha); freq.FTF = 0.5 * (1 - ratio) * frHz; freq.BPFO = Z * freq.FTF; freq.BPFI = Z * (frHz - freq.FTF); freq.BSF = (Dp / (2 * Bd)) * (1 - ratio^2) * frHz; freq.k = [freq.BPFO / frHz, freq.BPFI / frHz, ... freq.BSF / frHz, freq.FTF / frHz]; else error('mode参数必须是cwru或geo'); end end

这段代码里我特意把倍率数组freq.k也存了下来，因为在做频谱图标注时经常要把“几倍转频”直接打在图上，省得每次再除一次。CWRU模式下直接用了官方倍率，不管转速怎么变，只要转频输入正确，输出频率就是对的。

如果你手里有轴承图纸，想用理论公式算其他型号轴承的特征频率，就把mode改成geo，然后按照下面这样构造参数结构体：

params = struct('Bd', 0.3126, ... % 滚珠直径，英寸 'Dp', 1.537, ... % 节圆直径，英寸 'Z', 9, ... % 滚动体数量 'alpha', 0); % 接触角，度 frHz = 1797 / 60; freq = faultFreq('geo', params, frHz); disp(freq.BPFO);

这里有两个细节容易踩坑。第一，Bd和Dp的单位必须保持一致，你用英寸就全部用英寸，你用毫米就全部用毫米，因为公式里只出现两者的比值，单位不一致会导致结果完全错误。第二，接触角alpha的默认值在深沟球轴承里通常取0°，但如果遇到角接触球轴承，一定要查手册填实际接触角，否则计算结果会偏差很大。

3.3 示例运行：以1797 r/min工况为基准

CWRU数据集在0马力负载下，电机转速大约是1797转/分，转频frHz = 1797 / 60 ≈ 29.95 Hz。用CWRU官方倍率计算，可以得到：

这几个数值是后续所有频谱定位的基准。你拿到一个CWRU数据文件，做完FFT之后，先在频谱图上找这些频率位置是否出现明显峰值，就能快速判断和文件标签是否吻合。

如果用geo模式跑一遍，得到的外圈特征频率大约是107.37Hz，会比官方值高出不少。这个现象前面解释过了，不是代码问题，而是理论计算和实际工况的差异。CWRU数据做验证时，我更推荐直接用官方倍率，因为论文里通常也是这么对齐的。

4. 拿真实数据跑一遍：从.mat文件到频谱定位

4.1 数据文件的挑选和加载

CWRU数据集的.mat文件命名虽然直观，但不熟悉的人很容易拿错文件。驱动端轴承数据文件如105、106、107等，分别对应不同故障位置和尺寸，里面除了振动信号X105_DE_time，可能还有X105_FE_time、RPM等字段。读取方式很简单：

load('105.mat'); fs = 12000; % 驱动端采样率12kHz data = X105_DE_time; % 驱动端振动加速度信号

这里要注意采样率的选择。CWRU的驱动端数据采样率是12kHz，风扇端是48kHz，两者在文件里的信号命名也不同。如果你加载的是驱动端数据却误用48kHz采样率去计算频率轴，那频谱上的特征峰会全部偏到错误位置。我最初有一段时间就是在这里吃了亏，后来固定用文件对应的采样率，再没出过这种低级错误。

另外还要确认文件对应的转速。不同负载下转速不同，0、1、2、3马力工况对应的转速约为1797、1772、1750、1730转/分。文件标签里会写明负载条件，算转频时必须和工况匹配，不能所有文件都用1797算。有些数据文件的RPM字段可以直接读取，那就优先用RPM字段值。

4.2 FFT频谱定位的完整步骤

下面是一段可以直接跑通的流程，读取数据、计算特征频率、做FFT、在频谱上标注特征频率位置。

% 读取数据 load('130.mat'); % 外圈故障数据示例，0.007英寸，3点钟位置 fs = 12000; data = X130_DE_time; N = length(data); % 计算转频，动态读取文件中的RPM字段（如果存在） if exist('RPM', 'var') && ~isempty(RPM) rpm_value = RPM; else rpm_value = 1797; % 手动指定工况转速 end frHz = rpm_value / 60; % 计算特征频率 freq = faultFreq('cwru', [], frHz); % 做FFT并取单边频谱 L = min(N, 12000); % 取1秒数据足够 segment = data(1:L); win = hann(L); Y = fft(segment .* win); 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; plot(f, P1); xlim([0 500]); hold on; xline(freq.BPFO, 'r--', 'BPFO'); xline(freq.BPFI, 'g--', 'BPFI'); xline(freq.BSF, 'm--', 'BSF'); grid on; xlabel('频率 (Hz)'); ylabel('幅值'); legend('频谱', 'BPFO', 'BPFI', 'BSF');

这段代码有两个细节值得一提。一是加了hann窗，轴承振动信号里常常含有周期性冲击，加窗能减少频谱泄漏，特征峰更尖锐。二是取1秒数据做FFT而不是全部数据，在保证频率分辨率足够的前提下计算量更小，频率分辨率约1Hz，对特征频率定位已经够用。

如果你加载的是外圈故障文件，频谱图上应该能在91Hz附近看到一个明显的峰值，而145Hz和59Hz位置没有突出峰值，这说明文件标签和计算结果是吻合的。内圈故障文件则应该在145Hz附近出现明显峰值，而且这个峰值周围通常还会伴有一圈以转频为间隔的边带，这是内圈故障的典型特征。

4.3 边带和调制现象，没经验的容易忽略

很多人算对了特征频率，却在真实频谱里对不上号，问题往往出在调制边带上。内圈故障信号因为故障点随轴旋转，冲击幅度受载荷分布调制，在BPFI两侧会出现间隔为转频frHz的边带，看起来就像主峰周围的一梳子小峰。外圈故障虽然调制较轻，但有时也会出现间隔为保持架频率的边带。滚动体故障更复杂，故障点交替接触内外圈，调制频率可能是转频，也可能是保持架频率，而且由于滚动体位置滑动，特征峰还会出现加宽甚至有轻微漂移。

所以我实际识别时，不会只盯着特征频率那一根线，而是看它周围是否有规律的边带结构。如果BPFI附近边带不齐，就要怀疑是不是把内圈故障误判成了外圈故障。这类信号分析经验，只看理论公式是学不到的，必须亲手在频谱图上比对几组数据才能形成判断直觉。

5. 实际使用中的常见坑与排查

5.1 CWRU数据文件的那些坑

CWRU数据文件命名规则很规律，但真用起来还是有几个容易栽跟头的地方。比如0.007英寸和0.014英寸故障的文件号不同，对应关系依赖官方表格，记不住就容易对应错。另外部分文件是风扇端数据，采样率是48kHz，如果用12kHz去算频率轴，整个频谱就废了。还有一种常见情况：文件名里写的是正常轴承，但实际振动信号里依然有一些高频冲击成分，可能是电机本身或负载端引起的，这时候如果你没算特征频率就做频谱分析，很容易把无关峰值当成故障特征。

5.2 代码和计算层面的坑

代码层面的坑更多是细节问题。第一次用geo模式计算时，如果忘了把接触角从角度转弧度，算出来的结果可能完全不可用。我后来在代码里强制params.alpha按度为单位，内部再统一转弧度，至少保证外部输入时不用每次换算。

还有一个我经常遇到的坑是特征频率在频率分辨率不够的时候无法对齐。比如数据很短只有0.1秒，FFT频率分辨率为10Hz，91Hz和107Hz的差异就完全分辨不出来，很容易误判。解决办法是尽量截取超过1秒的数据，或者对原始信号做重叠分段、平均频谱，来获得更稳定的特征峰。

6. 代码还能怎么改造成更完整的小工具

6.1 自动判断故障类型

有了特征频率，再结合频谱峰值搜索，就可以让程序自动判断故障类型。思路是取频谱中几个特征频率附近的局部最大值，和该频段的高斯噪声底做对比，如果信噪比超过设定阈值，就认为该类型故障存在，可以直接输出到命令行：

[peak_bpfo, idx] = max(P1(idx_range_BPFO)); noise_floor = mean(P1(noise_idx)); if peak_bpfo / noise_floor > 5 fprintf('检测到外圈故障，特征频率 %.2f Hz，信噪比 %.2f dB\n', ... freq.BPFO, 10*log10(peak_bpfo/noise_floor)); end

这个逻辑再往后扩展，就能做成一个批量诊断脚本，一次性处理整个CWRU目录下几十个文件，把每段信号识别的结果和文件标注做比对，输出一个准确率报表，验证自己的诊断算法比手动看频谱高效得多。

6.2 包络谱模块，处理低频冲击被噪声淹没的情况

直接FFT看特征频率在故障早期往往不够灵敏，因为冲击能量分散在高频共振带上，低频特征峰会很小。这时候就需要包络谱处理：先对信号做带通滤波，选在结构共振明显的频率区间，比如2000到5000Hz，然后取希尔伯特变换的幅值，也就是包络信号，再对包络信号做FFT，特征频率就会在包络谱里变得非常清晰。

这个思路可以从CWRU数据扩展到现场采集数据。试验台数据相对干净，现场信号里混着齿轮啮合频率、电磁干扰、随机噪声，不做包络谱几乎没办法直接看特征频率。只要把这个模块加上，整套代码就从“算频率的小脚本”升级成了“能处理实际故障信号的小工具”。

另外，滚动体故障因为特征峰幅值小且容易波动，用包络谱的效果往往比直接频谱好很多，如果你手中正好有滚动体故障的数据文件，建议优先试一下这个流程。我在实际项目中，靠包络谱捕捉到过好几起早期外圈点蚀案例，都是直接FFT完全看不出来的。

这套代码我后续还在不断维护，比如加入滑移率修正参数、自动估计转频、批量导出报告等功能。核心的faultFreq函数基本稳定，你可以直接拿去用。算特征频率这件事本身不复杂，但把它和真实数据验证、包络分析串成一条完整的诊断链路，才是让它发挥价值的关键。

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