简介面向机械工程、设备维护与信号处理学习者帮助解决轴承故障诊断中特征频率计算与分析的实际问题尤其适合需要把理论公式落实到代码层面的初学者。内容基于凯斯西储大学公开轴承数据提供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 × FTFBPFI N × (fr - FTF)BSF (Dp / (2 × Bd)) × fr × (1 - (Bd / Dp × cos α)²)其中N是滚动体数量Bd是滚动体直径Dp是节圆直径α是接触角。以CWRU常用的6205轴承为例滚动体数量N9滚珠直径Bd约为0.3126英寸节圆直径Dp约为1.537英寸接触角按0°处理。把数值代入公式可以得到理论倍率大约是FTF ≈ 0.3983 × frBPFO ≈ 3.5847 × frBPFI ≈ 5.4150 × frBSF ≈ 2.3565 × fr但CWRU官网给出的官方倍率却是BPFO3.048、BPFI4.857、BSF1.994、FTF0.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 - structgeo模式下需要包含字段 % Bd(滚珠直径)、Dp(节圆直径)、Z(滚珠数量)、alpha(接触角单位度) % frHz - 轴的转频单位HzfrHz 转速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官方倍率计算可以得到故障类型倍率特征频率Hz外圈 BPFO3.04891.28内圈 BPFI4.857145.47滚动体 BSF1.99459.72保持架 FTF0.39811.92这几个数值是后续所有频谱定位的基准。你拿到一个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/21); 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去算频率轴整个频谱就废了。还有一种常见情况文件名里写的是正常轴承但实际振动信号里依然有一些高频冲击成分可能是电机本身或负载端引起的这时候如果你没算特征频率就做频谱分析很容易把无关峰值当成故障特征。易错点具体表现解决办法文件标注与故障类型不匹配频谱特征峰与理论频率对不上先核对官方文件索引表确认故障位置和故障尺寸采样率用错频率轴整体偏移特征峰位置不对驱动端用12kHz风扇端用48kHz固定写死并加注释转速RPM与实际值不符转频偏差导致所有特征频率偏差优先读取.mat文件里的RPM字段或根据负载工况手动匹配加载后信号字段选错用了风扇端信号做驱动端分析检查字段名DE_TIME是驱动端FE_TIME是风扇端5.2 代码和计算层面的坑代码层面的坑更多是细节问题。第一次用geo模式计算时如果忘了把接触角从角度转弧度算出来的结果可能完全不可用。我后来在代码里强制params.alpha按度为单位内部再统一转弧度至少保证外部输入时不用每次换算。还有一个我经常遇到的坑是特征频率在频率分辨率不够的时候无法对齐。比如数据很短只有0.1秒FFT频率分辨率为10Hz91Hz和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函数基本稳定你可以直接拿去用。算特征频率这件事本身不复杂但把它和真实数据验证、包络分析串成一条完整的诊断链路才是让它发挥价值的关键。本文还有配套的精品资源点击获取