简介:这份资源面向机械故障诊断方向的学生、科研人员与工程技术人员,提供一套基于MATLAB的轴承故障分析代码,重点针对滚动轴承内圈与外圈故障。资源围绕同一转动频率信号采用不同处理方法、以及同一方法对不同转动频率信号进行对比,涵盖时域无量纲参数分析、FFT分析、Hilbert包络谱分析、Haar小波分析与数学形态学分析等思路,帮助读者理解各方法在故障频率识别上的差异与适用边界。压缩包共30个文件,约1MB,包含5个m脚本文件、11个dat数据文件与14张jpg结果图,脚本对应不同分析流程,数据文件为实测振动信号,图片则记录频谱、包络图等中间结果,便于对照验证。目前已有12920人学习下载,适合希望快速复现轴承故障诊断实验、比较多种信号处理方法效果并积累排错经验的读者参考。
1. 轴承故障诊断的四种方法:从 FFT 到包络图再到 Haar 小波,一套 MATLAB 代码全跑通
振动信号里藏着轴承的健康状态,但原始时域波形往往看不出门道。这份 MATLAB 代码包围绕四种经典方法展开:直接 FFT 频谱分析、Hilbert 包络解调、Haar 小波变换、以及包络图与频谱的联合诊断。它解决的核心问题是:拿到一段轴承振动数据后,怎么一步步从时域波形走到能判断故障类型的特征频率。适合正在做旋转机械状态监测的工程师、需要快速验证算法的研究生,以及手头有实测数据但不确定该用哪种分析路径的从业者。代码不依赖额外工具箱,基础 MATLAB 环境即可运行,四种方法各自独立成脚本,方便对比和裁剪。
2. FFT 与包络解调:先看清频谱,再拆出冲击
2.1 为什么先做 FFT:从时域到频域的第一次筛选
轴承故障的振动信号本质上是周期性的冲击序列调制在系统共振频率上。时域波形里,早期故障的冲击幅值很小,淹没在噪声和转频振动里,肉眼几乎分辨不出。FFT 的作用是把信号从时间轴搬到频率轴,让周期成分以谱峰形式暴露出来。
但直接对原始信号做 FFT 有个前提:信号必须是稳态的。如果转速波动大或者负载变化剧烈,频谱会糊成一片。我一般会先截取一段转速稳定的数据,通常取 1 到 2 秒,采样率保证在 10 kHz 以上,这样高频共振区才不会被截断。
代码里 FFT 部分的核心逻辑是:去均值、加窗、做 FFT、取单边谱、转 dB。加窗这一步很多人会忽略,矩形窗在非整周期采样时泄漏严重,汉宁窗是更稳妥的选择。
% FFT 频谱分析:输入信号 x,采样率 fs x = x - mean(x); % 去直流分量 N = length(x); w = hanning(N); % 汉宁窗抑制频谱泄漏 xw = x .* w; X = fft(xw); P = abs(X(1:N/2+1)) / sum(w) * 2; % 单边幅值谱,窗增益归一化 f = (0:N/2) * fs / N; P_db = 20 * log10(P + eps); % 转 dB,eps 防止 log(0) plot(f, P_db); xlabel('频率 (Hz)'); ylabel('幅值 (dB)');这段代码里,sum(w)是窗函数的增益,不做归一化的话不同窗长之间幅值没法比。eps是 MATLAB 里最小的正浮点数,加它是为了防止零值取对数变成负无穷。取N/2+1个点是因为实信号的频谱对称,只看一半就够了。
FFT 能告诉你哪些频率有能量,但轴承故障的特征频率往往不在低频段,而是在高频共振区被调制。这时候光看 FFT 还不够,需要包络解调把调制信息拆出来。
2.2 Hilbert 包络解调:把冲击从载波上剥下来
包络解调的思路是:轴承故障产生的冲击会激发系统的高频共振,共振信号像一个载波,冲击的重复频率是调制信号。直接看频谱,你看到的是共振峰,看不到故障特征频率。Hilbert 变换可以构造解析信号,取模就得到包络。
具体步骤是:先带通滤波,把共振频带选出来;然后做 Hilbert 变换取包络;最后对包络再做一次 FFT,得到包络谱。包络谱里的峰值对应的就是故障特征频率,比如内圈故障频率 BPFI、外圈故障频率 BPFO、滚动体故障频率 BSF。
% Hilbert 包络解调 bp_low = 2000; bp_high = 5000; % 带通范围,需根据共振区调整 [b, a] = butter(4, [bp_low bp_high]/(fs/2), 'bandpass'); x_bp = filtfilt(b, a, x); % 零相位滤波,避免相位失真 env = abs(hilbert(x_bp)); % Hilbert 包络 env = env - mean(env); % 去直流 N_env = length(env); Env_spec = abs(fft(env .* hanning(N_env))); f_env = (0:N_env/2) * fs / N_env; plot(f_env, Env_spec(1:N_env/2+1)); xlabel('频率 (Hz)'); ylabel('包络幅值');filtfilt而不是filter,是因为filter会引入相位延迟,包络的时序位置会偏移,对后续定位冲击时刻不利。带通范围bp_low和bp_high是最关键的参数,选错了整个包络谱都是垃圾。常见做法是先看 FFT 频谱,找到共振峰所在频段,然后把带通设在这个峰的两侧。如果共振区不明显,可以用谱峭度图来辅助定位,但那需要额外代码,这份包里没包含。
包络谱出来之后,对照轴承的故障特征频率表,看峰值出现在哪里。内圈故障的包络谱峰在 BPFI 及其谐波处,外圈在 BPFO 处,滚动体在 BSF 处。如果峰值不明显,可能是带通选错了,或者故障还处于早期,冲击能量太弱。
注意:包络解调对带通参数非常敏感,同一段数据换一个带通范围,包络谱可能完全变样。建议先用 FFT 确认共振区,再动手调带通。
3. Haar 小波变换:用多分辨率分析抓冲击时刻
3.1 Haar 小波为什么适合冲击检测
Haar 小波是最简单的小波基,它的优势在于对突变信号极其敏感。轴承故障产生的冲击在时域上是一个陡峭的跳变,Haar 小波的阶跃特性正好匹配这种波形。相比 db 小波或 morlet 小波,Haar 的滤波器系数只有两个,计算量小,分解和重构都很快。
小波变换的核心思想是多分辨率分析:把信号分解到不同尺度上,粗尺度看趋势,细尺度看细节。轴承的冲击成分通常藏在细节系数里,尤其是第一层和第二层的细节系数。
代码里用的是单层离散小波变换,把信号拆成近似系数和细节系数。近似系数是低频趋势,细节系数是高频细节。对细节系数做阈值处理,把噪声压掉,再重构回去,冲击成分就凸显出来了。
% Haar 小波分解与重构 wname = 'haar'; [cA, cD] = dwt(x, wname); % 单层分解 % 对细节系数做软阈值处理 sigma = median(abs(cD)) / 0.6745; % 噪声标准差估计 thr = sigma * sqrt(2 * log(length(cD))); % 通用阈值 cD_thr = sign(cD) .* max(abs(cD) - thr, 0); % 软阈值 x_denoised = idwt(cA, cD_thr, wname); % 重构 % 对重构信号做包络分析 env_w = abs(hilbert(x_denoised)); env_w = env_w - mean(env_w);median(abs(cD)) / 0.6745是噪声标准差的鲁棒估计,0.6745 是正态分布中位数绝对偏差到标准差的转换系数。阈值公式sigma * sqrt(2 * log(n))是 Donoho 提出的通用阈值,n 是细节系数长度。软阈值比硬阈值更平滑,不会在重构信号里引入新的跳变。
分解层数可以根据需要调整,dwt只做一层,多层用wavedec。但层数不是越多越好,层数太多,细节系数里的冲击成分会被稀释,反而不容易检测。我一般做到三层就停,再往下意义不大。
3.2 小波重构后的包络谱对比
小波去噪之后再做包络解调,和直接做包络解调有什么区别?区别在于信噪比。如果原始信号噪声很大,直接做 Hilbert 包络,包络里混着大量噪声,包络谱的峰值可能被淹没。小波去噪先把高频噪声压掉,再取包络,包络谱的峰值会更干净。
但小波去噪也有风险:阈值设得太高,冲击成分会被一起压掉;阈值设得太低,噪声压不干净。软阈值虽然平滑,但会低估冲击幅值。实际用的时候,我一般会把原始包络谱和小波去噪后的包络谱叠在一起看,如果小波去噪后峰值更明显,说明去噪有效;如果峰值反而变小了,说明阈值过头了。
% 对比原始包络谱和小波去噪后的包络谱 env_raw = abs(hilbert(x)); env_raw = env_raw - mean(env_raw); spec_raw = abs(fft(env_raw .* hanning(length(env_raw)))); spec_raw = spec_raw(1:floor(length(env_raw)/2)+1); env_w = env_w(1:length(env_raw)); % 对齐长度 spec_w = abs(fft(env_w .* hanning(length(env_w)))); spec_w = spec_w(1:floor(length(env_w)/2)+1); f_axis = (0:length(spec_raw)-1) * fs / (2*length(spec_raw)-2); plot(f_axis, spec_raw, 'b', f_axis, spec_w, 'r'); legend('原始包络谱', '小波去噪后包络谱');这段代码里,floor(length(env_raw)/2)+1是取单边谱的点数,f_axis的计算要跟单边谱的长度匹配。两条曲线叠在一起,红色峰值如果比蓝色高,说明小波去噪起了正面作用。
提示:Haar 小波的细节系数对冲击敏感,但也对噪声敏感。如果信号里噪声占主导,先做小波去噪再取包络,比直接取包络更稳。
4. 四种方法的联合诊断流程:参数怎么设、结果怎么读
4.1 从原始数据到诊断结论的完整链路
四种方法不是孤立的,实际诊断时是一条链路:先看时域波形确认有没有冲击成分,再做 FFT 找共振频带,然后带通滤波加 Hilbert 包络解调得到包络谱,最后用 Haar 小波去噪提升信噪比。如果包络谱里能看到明显的故障特征频率及其谐波,基本可以下结论。
代码包里四个脚本可以按顺序跑,也可以单独跑。我一般会先跑 FFT 脚本,把频谱图存下来,肉眼找共振峰。找到之后把带通参数填进包络解调脚本,跑出包络谱。如果包络谱不干净,再跑小波去噪脚本,把去噪后的信号替换进去。
参数设置上,采样率fs必须和实际数据一致,这个搞错了后面全错。带通范围bp_low和bp_high根据 FFT 频谱里的共振峰位置来定,通常取峰两侧各 500 Hz 到 1 kHz。小波分解层数默认三层,阈值系数用通用阈值,不需要手动调。
| 方法 | 关键参数 | 适用场景 | 输出 |
|---|---|---|---|
| FFT | 窗函数、采样率 | 稳态信号、初步排查 | 频谱图 |
| Hilbert 包络 | 带通范围 | 冲击调制信号 | 包络谱 |
| Haar 小波 | 分解层数、阈值 | 噪声大、冲击弱 | 去噪信号 |
| 联合诊断 | 以上全部 | 实测数据、故障定位 | 特征频率 |
4.2 故障特征频率的计算与对照
包络谱出来之后,怎么判断是内圈故障还是外圈故障?需要计算轴承的故障特征频率。这些频率由轴承几何参数和转速决定,公式如下:
- 外圈故障频率 BPFO = (n/2) × fr × (1 - (d/D) × cosα)
- 内圈故障频率 BPFI = (n/2) × fr × (1 + (d/D) × cosα)
- 滚动体故障频率 BSF = (D/(2d)) × fr × (1 - ((d/D) × cosα)^2)
其中 n 是滚动体个数,fr 是转频,d 是滚动体直径,D 是节径,α 是接触角。这些参数从轴承型号手册里查,代码包里没有内置轴承参数表,需要自己填。
% 轴承故障特征频率计算 n = 8; % 滚动体个数 fr = 30; % 转频 (Hz) d = 8; % 滚动体直径 (mm) D = 40; % 节径 (mm) alpha = 0; % 接触角 (度) BPFO = (n/2) * fr * (1 - (d/D) * cosd(alpha)); BPFI = (n/2) * fr * (1 + (d/D) * cosd(alpha)); BSF = (D/(2*d)) * fr * (1 - ((d/D) * cosd(alpha))^2); fprintf('BPFO = %.2f Hz\n', BPFO); fprintf('BPFI = %.2f Hz\n', BPFI); fprintf('BSF = %.2f Hz\n', BSF);cosd是角度制余弦,不是弧度制。算出来的频率和包络谱里的峰值对照,哪个接近就是哪个故障。如果包络谱里同时出现多个特征频率,可能是复合故障,需要进一步分离。
注意:转频 fr 不是电机的额定转速,而是实际转速。如果转速波动大,特征频率会漂移,包络谱的峰值会变宽。这种情况下需要做阶次分析,但阶次分析不在这个代码包的范围内。
5. 避坑与排查:五个让包络谱变成废纸的常见操作
5.1 带通范围拍脑袋定,包络谱全是噪声
现象:包络谱里看不到任何明显的峰值,或者峰值出现在和故障特征频率无关的位置。
原因:带通范围没有对准共振频带,把噪声或者转频谐波选进来了。Hilbert 包络对带通极其敏感,带通选错,包络就是噪声的包络。
解决:先跑 FFT,把频谱图放大看高频段,找到能量最高的共振峰,把带通中心设在这个峰上,带宽取峰宽的 2 到 3 倍。如果不确定,可以多试几组带通参数,对比包络谱的干净程度。
5.2 忘了去直流,包络谱零频处一个大峰
现象:包络谱在 0 Hz 处有一个巨大的峰值,其他频率的峰值被压得看不见。
原因:Hilbert 包络的均值不为零,FFT 之后零频分量占了主导。包络信号本身是单边谱,均值不为零是正常的,但做 FFT 之前必须去均值。
解决:在hilbert之后加一行env = env - mean(env)。这行代码在包络解调脚本里有,但如果自己改代码容易漏掉。
5.3 小波阈值设太高,冲击被一起压掉
现象:小波去噪后的信号变得很平滑,包络谱里的峰值反而比去噪前更小。
原因:通用阈值sigma * sqrt(2*log(n))在噪声标准差估计偏大的时候会过高,软阈值把冲击成分也压缩了。
解决:把阈值系数从 1.0 降到 0.5 或者 0.3,重新跑一遍对比。或者改用硬阈值,硬阈值保留冲击幅值,但会引入振铃。我一般先用软阈值跑一遍,如果峰值被压得太狠,换硬阈值再试。
5.4 采样率填错,特征频率对不上
现象:包络谱里的峰值位置和计算出来的 BPFO、BPFI 对不上,差了一个固定的比例。
原因:采样率fs和实际数据不一致,导致频率轴整体缩放。这种情况在拿别人数据的时候特别常见,数据文件里没写采样率,自己估了一个错的。
解决:检查数据来源,确认采样率。如果实在不知道,可以用已知的转频反推:在 FFT 频谱里找到转频的峰值,用峰值位置反推采样率。但这个方法只在转频已知的情况下有效。
5.5 直接对原始信号做包络,没做带通滤波
现象:包络谱里低频段全是转频及其谐波,故障特征频率被淹没。
原因:原始信号里转频成分的能量远大于冲击成分,不做带通滤波直接取包络,包络主要反映的是转频的幅度调制,不是冲击的调制。
解决:包络解调之前必须做带通滤波,把转频和低频振动滤掉,只保留共振频带。带通滤波用filtfilt而不是filter,避免相位失真。
6. 进阶技巧:用包络谱的谐波结构区分故障类型
包络谱出来之后,光看基频峰值还不够,谐波结构才是区分故障类型的关键。内圈故障的包络谱里,BPFI 的谐波通常比较丰富,因为内圈旋转导致故障点进出承载区,冲击幅值被转频调制。外圈故障的包络谱里,BPFO 的基频峰值突出,谐波相对少,因为外圈固定,冲击幅值稳定。
我一般会看包络谱里前三个峰值的相对高度。如果 BPFI 的二次谐波比基频还高,内圈故障的可能性很大。如果 BPFO 的基频一枝独秀,外圈故障的概率高。滚动体故障的包络谱里,BSF 的峰值通常比较低,而且会出现保持架频率 FTF 的边带。
% 包络谱峰值提取与谐波对比 [pks, locs] = findpeaks(spec_w, 'MinPeakHeight', max(spec_w)*0.1); f_peaks = f_axis(locs); amp_peaks = pks; % 找 BPFI 及其二次谐波 [~, idx1] = min(abs(f_peaks - BPFI)); [~, idx2] = min(abs(f_peaks - 2*BPFI)); if ~isempty(idx1) && ~isempty(idx2) ratio = amp_peaks(idx2) / amp_peaks(idx1); fprintf('BPFI 二次谐波与基频幅值比: %.2f\n', ratio); if ratio > 1 fprintf('内圈故障可能性大\n'); end endfindpeaks的MinPeakHeight设成最大峰值的 10%,是为了过滤掉噪声引起的小峰。min(abs(f_peaks - BPFI))是找最接近 BPFI 的峰值位置,因为 FFT 的频率分辨率有限,峰值不一定正好落在 BPFI 上。
这套流程跑下来,从原始数据到故障类型判断,快的话十几分钟就能出结果。但前提是数据质量过得去,采样率、转速、轴承参数都齐全。如果数据本身有问题,再好的算法也救不回来。从那以后我每次拿到新数据,都强制先跑一遍 FFT 确认频谱干净,再往下走包络和小波。希望帮到你。
本文还有配套的精品资源,点击获取