轴承故障特征频率计算与MATLAB频谱诊断实战

简介这份资源面向机械工程、设备维护与信号处理方向的学习者和研究人员聚焦轴承故障诊断中的特征频率计算问题提供可直接运行的MATLAB实现方案。压缩包共3个文件包含1个m脚本文件与2张png图片整体约91KB脚本用于特征频率计算与频谱仿真图片则展示计算截图及JEMSKF深沟球轴承结构便于对照理解。内容围绕凯斯西储大学轴承数据展开涉及基本旋转频率、滚动体通过频率、径向跳动频率等关键指标并结合傅里叶变换、数字滤波与调制解调等信号处理手段帮助读者从振动信号中提取故障特征。配套的步骤分析基础讲解覆盖数据预处理、特征提取与故障判断流程适合初学者按图索骥。目前已有4135人学习下载可作为轴承故障分析入门与课题研究的实用参考工具。1. 轴承故障特征频率这件事为什么值得你花一个下午跑通如果你手头有一台电机、一个加速度传感器和一段振动数据却不确定频谱里那几根尖峰到底代表内圈、外圈还是滚动体故障那这套轴承故障特征频率的 MATLAB 代码就是为你准备的。它解决的不是能不能采集到振动的问题而是采到之后怎么把故障类型和特征频率对上号的问题。凯斯西储的轴承数据集在旋转机械故障诊断圈子里几乎是入门标配采样率、转速、故障尺寸都标得清清楚楚配合特征频率公式一算频谱上该出现峰的位置就能提前预测出来。适合谁做设备状态监测的、写诊断算法的、带学生做课题的只要你能装 MATLAB 并会改几个参数就能在本地复现整套流程。我见过太多人卡在数据下载了但不知道怎么用这一步这篇就把这条路铺平。2. 特征频率到底怎么算从几何参数到频谱尖峰2.1 四个公式背后的物理含义轴承故障特征频率不是拍脑袋定的它由轴承几何和转速共同决定。核心参数有四个滚动体直径 d、节径 D、接触角 α、滚动体个数 n。外圈固定、内圈随轴旋转时各部位的通过频率如下外圈故障频率 BPFO (n/2) × fr × (1 - (d/D)cosα)内圈故障频率 BPFI (n/2) × fr × (1 (d/D)cosα)滚动体故障频率 BSF (D/(2d)) × fr × (1 - ((d/D)cosα)²)保持架故障频率 FTF (fr/2) × (1 - (d/D)cosα)其中 fr 是轴的转频单位 Hz等于转速 rpm 除以 60。注意 BSF 是滚动体自转频率实际频谱里常看到的是它的两倍因为滚动体每转一圈会进出承载区各一次。这个细节很多人第一次做会漏掉导致对不上峰。凯斯西储数据里驱动端轴承是 6205-2RS 深沟球轴承风扇端是 6203-2RS。以 6205 为例d0.3126 英寸D1.537 英寸n9α≈0。这些几何参数在代码里通常写成常量你换成自己的轴承时只改这几个数就行。2.2 用 MATLAB 把公式变成可调用的函数不要每次手算写成一个函数输入转速和几何参数输出四个频率。下面是我常用的写法function [BPFO, BPFI, BSF, FTF] bearing_freq(rpm, d, D, n, alpha) % bearing_freq 计算轴承四个故障特征频率 % rpm : 轴转速单位转/分 % d : 滚动体直径单位与D一致 % D : 节径单位与d一致 % n : 滚动体个数 % alpha : 接触角单位弧度深沟球轴承近似为0 fr rpm / 60; % 转频 Hz ratio (d / D) * cos(alpha); % 直径比项 BPFO (n / 2) * fr * (1 - ratio); BPFI (n / 2) * fr * (1 ratio); BSF (D / (2 * d)) * fr * (1 - ratio^2); FTF (fr / 2) * (1 - ratio); end逻辑说明先算转频再把直径比单独提出来四个公式共用。参数说明alpha 用弧度如果你手头是角度先乘 pi/180。这个函数不依赖任何工具箱纯基础运算放到任何 MATLAB 版本都能跑。调用示例1730 rpm 下 6205 轴承rpm 1730; d 0.3126; D 1.537; n 9; alpha 0; [BPFO, BPFI, BSF, FTF] bearing_freq(rpm, d, D, n, alpha); fprintf(BPFO%.2f BPFI%.2f BSF%.2f FTF%.2f\n, BPFO, BPFI, BSF, FTF);跑出来 BPFO 约 103.4 HzBPFI 约 156.1 HzBSF 约 135.4 HzFTF 约 11.5 Hz。这几个数就是你后面在频谱上要找的目标。2.3 把理论频率标到实测频谱上光有数字没用得画出来对比。假设你已经加载了一段振动信号 x 和采样率 fs做 FFT 并叠加理论频率线N length(x); X abs(fft(x - mean(x))) / N * 2; % 去直流并归一化 f (0:N-1) * (fs / N); % 频率轴 half 1:floor(N/2); % 只看前半段 plot(f(half), X(half)); hold on; xline(BPFO, r--, BPFO); xline(BPFI, g--, BPFI); xline(BSF, b--, BSF); xlabel(Frequency (Hz)); ylabel(Amplitude); legend(Spectrum,BPFO,BPFI,BSF);逻辑说明去均值是为了压掉 0 Hz 的大直流分量否则频谱第一根柱会盖住一切。归一化除以 N 再乘 2 是单边谱的标准做法。参数说明fs 必须和采集时一致凯斯西储驱动端数据常见 12 kHz 和 48 kHz 两档用错 fs 频率轴整体缩放峰就对不上了。xline 需要 R2018b 及以上老版本用 plot([BPFO BPFO], ylim, r--) 代替。3. 从原始数据到诊断结论一套可复现的处理链3.1 数据加载与通道选择凯斯西储的数据文件通常是 .mat 格式里面存着驱动端和风扇端的振动信号还有转速。加载时先 whos 看一眼变量名别硬猜。load(normal_0.mat); % 正常基线数据 whos; % 查看变量名和维度常见变量名是 X123_DE_time 这种格式DE 表示驱动端FE 表示风扇端。选驱动端数据做分析因为故障轴承通常装在驱动端。如果文件里同时有 DE 和 FE用 DE。加载后先看长度和采样率长度除以 fs 就是时长正常应该有几秒到十几秒。3.2 分段、加窗与平均谱直接对整段数据做一次 FFT噪声大、谱线毛刺多。工程上常用分段平均把信号切成若干段每段加汉宁窗各自 FFT 后取平均。这样随机噪声被压下去周期性故障特征更突出。fs 12000; % 采样率 segLen 2048; % 每段点数 overlap segLen / 2; % 50% 重叠 win hann(segLen); % 汉宁窗 nSeg floor((length(x) - segLen) / (segLen - overlap)) 1; P zeros(segLen, 1); for k 1:nSeg idx (k-1)*(segLen-overlap) (1:segLen); seg x(idx) .* win; P P abs(fft(seg - mean(seg))).^2; end P P / nSeg; f (0:segLen-1) * (fs / segLen);逻辑说明加窗是为了减少截断引起的频谱泄漏汉宁窗是通用选择。重叠 50% 能让分段之间过渡更平滑。取平方得到功率谱平均后噪声方差降低。参数说明segLen 决定频率分辨率2048 点在 12 kHz 下分辨率约 5.86 Hz够看低频特征如果想分辨密集边带加到 4096 或 8192。overlap 不要超过 segLen 的 75%否则段间相关性太强平均效果打折。3.3 包络谱把高频载波里的故障信息挖出来轴承早期故障的冲击很弱直接看原始频谱可能被其他振动盖住。常用做法是选一个高频共振带做带通滤波然后取包络再做 FFT这样故障特征频率及其谐波会以低频形式清晰出现。band [2000 4000]; % 共振带按实际调整 [b, a] butter(4, band / (fs/2), bandpass); xf filtfilt(b, a, x); % 零相位滤波 env abs(hilbert(xf)); % 希尔伯特包络 env env - mean(env); N length(env); E abs(fft(env)) / N * 2; fEnv (0:N-1) * (fs / N); half 1:floor(N/2); plot(fEnv(half), E(half)); xline(BPFO, r--); xline(BPFI, g--);逻辑说明filtfilt 做双向滤波避免相位失真。hilbert 得到解析信号取模就是包络。包络谱上 BPFO 或 BPFI 的基频及 2 倍、3 倍频会成串出现这是判据。参数说明band 选在系统共振区通常看原始频谱里 2-5 kHz 哪段能量高就选哪段。butter 的阶数 4 够用太高容易数值不稳定。如果包络谱上基频不明显但谐波清楚说明故障已经比较严重。3.4 用表格对照不同故障类型的频谱表现跑完几组数据后把结果整理成表方便横向对比。下面是我做 6205 轴承、1730 rpm 时的典型观察故障类型理论基频 (Hz)包络谱突出位置边带特征外圈103.4103、207、310边带少峰较干净内圈156.1156、312、468围绕转频 fr 有边带滚动体135.4270、5402×BSF边带间隔为 FTF正常—无明显周期峰以转频及谐波为主这张表不是死的转速变了频率跟着变但规律一致。内圈故障因为随轴旋转载荷方向变化频谱上会出现以转频为间隔的边带这是区分内圈和外圈的关键。滚动体故障看 2×BSF别只盯 BSF。4. 避坑与排查那些让频谱对不上号的常见问题4.1 采样率用错频率轴整体偏移现象理论 BPFO 是 103 Hz但频谱上峰出现在 51 Hz 或 206 Hz 附近。原因fs 设成了实际的一半或两倍。凯斯西储数据有 12 kHz 和 48 kHz 两种文件名里通常有标注但有人直接默认 12 kHz 去处理 48 kHz 的数据。解决加载后先看数据说明或文件命名确认 fs 再算频率轴。如果不确定用已知转频反推正常数据里转频 fr 的峰应该出现在 rpm/60 附近对不上就是 fs 错了。4.2 转速没换算成转频现象所有特征频率都偏大 60 倍。原因公式里 fr 是 Hz有人直接把 rpm 代进去了。解决fr rpm/60这一步别省。1730 rpm 对应 28.83 Hz不是 1730 Hz。4.3 包络谱里找不到基频只有谐波现象包络谱上 207 Hz 有峰但 103 Hz 几乎没有。原因故障冲击太强基频被调制或淹没或者带通频带选得不对。解决先确认带通频带覆盖了共振区可以画原始频谱看哪段能量集中。如果谐波清楚直接用谐波间隔反推基频间隔就是 BPFO。另外检查滤波器是否把低频特征也滤掉了带通下限不要设太高。4.4 分段长度选得太短频率分辨率不够现象BPFO 和 BPFI 靠得近时频谱上两个峰糊在一起分不开。原因segLen 太小分辨率 fs/segLen 太大。12 kHz 下 1024 点分辨率约 11.7 Hz而 BPFO 和 BPFI 可能只差 50 Hz勉强能分但边带看不清。解决segLen 至少 4096分辨率到 2.9 Hz边带间隔 FTF 约 11.5 Hz 也能分辨。代价是每段时长增加数据短的话段数变少平均效果下降自己权衡。4.5 忽略传感器安装方向现象同一故障水平方向测得到垂直方向测不到。原因轴承故障冲击的方向性传感器没装在承载方向。解决驱动端轴承通常径向承载传感器应装在径向最好是水平或垂直其中一个能收到最强信号的方向。如果只有一个传感器先试径向。这个坑在实测里很常见但用公开数据集时因为数据已经采好反而容易忽略等到自己上手采集就翻车。5. 进阶技巧用谱峭度自动选带通频带前面带通频带靠人眼看频谱选换一组数据就得重选。谱峭度能自动找到冲击成分最强的频带省去手动调参。核心思路是计算不同频带的峭度值峭度最大的频带就是共振带。function [fc, bw] select_band(x, fs) % 简易谱峭度选带在若干候选带宽里找峭度最大者 bwList [500 1000 2000]; % 候选带宽 fcList 1000:500:(fs/2 - 1000); % 候选中心频率 bestK -inf; fc 0; bw 0; for b bwList for c fcList lo max(c - b/2, 1); hi min(c b/2, fs/2 - 1); [bb, aa] butter(4, [lo hi]/(fs/2), bandpass); y filtfilt(bb, aa, x); k kurtosis(y); if k bestK bestK k; fc c; bw b; end end end end逻辑说明峭度对冲击成分敏感正常信号峭度接近 3含冲击时升高。遍历候选中心和带宽取峭度最大的组合。参数说明bwList 和 fcList 按你的 fs 调整fs 高的时候候选可以密一些。这个简易版够用如果要做完整谱峭度图可以用 kurtogram 函数需要 Wavelet Toolbox。选完带通后接前面的包络谱流程自动化程度就上来了。我自己的习惯是拿到一组新数据先跑一遍谱峭度选带再画包络谱最后对照理论频率表。这样即使转速或轴承型号变了流程不用大改。还有一点理论频率算出来是给你指方向的实际频谱上峰可能偏几赫兹因为转速有波动、轴承有滑移别死磕小数点后一位看趋势和倍频关系更靠谱。希望帮到你。本文还有配套的精品资源点击获取