简介:这份资源面向从事振动分析与旋转机械故障诊断的工程师、研究人员及高校学生,提供一套基于MATLAB的阶次分析完整实现代码。阶次分析可将随时间变化的振动信号转换为随旋转角度变化的阶次信号,从而有效区分机械故障特征频率与其他频率成分,适用于变转速工况下的故障诊断场景。压缩包共3个文件,包含1个m脚本、1个mat数据文件和1个jpg图示,整体约13.72MB,其中脚本承载阶次跟踪与频谱分析的核心算法,mat文件保存采集或预处理后的振动数据,图示则直观展示阶次分析的概念与流程。目前已有357人学习关注。读者可借助这套代码快速理解阶次转换、信号预处理与频谱识别的实现思路,对照数据文件复现分析流程,并在此基础上修改参数以适配自身采集的振动信号,为旋转机械故障诊断研究提供可直接运行的参考范例。
1. 从一段被噪声淹没的振动信号说起:阶次分析到底在解决什么
电机升降速试验台上采回来一组振动加速度信号,转速从 600 rpm 拉到 3000 rpm,采样率 25.6 kHz,数据看着挺干净,一做 FFT 就傻眼——谱线糊成一片,转频对应的峰值随着转速漂移,根本对不上号。这不是数据坏了,而是旋转机械的振动本质上是角度域的周期现象,用等时间间隔采样再去做频谱,转速一变,相位就错位,能量自然散开。阶次分析(Order Analysis)要干的事,就是把时间域的非平稳信号重采样到角度域,让「每转发生几次」这件事变成一根根稳定的阶次谱线,1 阶就是转频,2 阶可能是不对中,0.5 阶往往指向油膜涡动或气流激振。
这份 MATLAB 实现阶次分析完整代码文件,核心就是围绕「转速脉冲 → 角度重采样 → 阶次谱」这条链路展开的。它适合两类人:一类是做旋转机械状态监测、故障诊断的工程师,手里有键相脉冲和振动通道数据,想快速验证阶次谱;另一类是在校做振动分析课题的人,需要一套能跑通、能改参数、能对照结果的参考实现。代码不依赖任何工具箱之外的付费模块,信号处理工具箱是主要依赖,常见做法是 R2016b 以上版本都能直接跑。下面我按「先跑通、再拆解、最后避坑」的顺序,把这份资源里真正值得抄作业的部分讲清楚。
2. 阶次分析的理论骨架:从等时间采样到等角度重采样
2.1 为什么 FFT 在变速工况下会失效
等时间采样下,采样间隔 Δt 固定,转速 n(t) 变化时,每转对应的采样点数 N = 60·fs / n 是变的。假设 fs = 25.6 kHz,转速 600 rpm 时每转 2560 点,3000 rpm 时每转 512 点,差了 5 倍。对整段数据做 FFT,相当于把不同「每转点数」的片段强行拼在一起,转频峰值被抹平、展宽,阶次信息就丢了。更麻烦的是,齿轮啮合、轴承故障这些特征频率都是转频的整数或非整数倍,转速一漂,它们在时间谱上跟着跑,根本没法用固定频率窗去盯。
阶次分析的做法是引入一个参考通道——通常是键相(Keyphasor)或转速脉冲信号,每转产生一个或若干个脉冲。用这个脉冲去标定「角度」,再把振动信号按等角度间隔重新采样。重采样之后,每转固定采 N 个点,转速变化被归一化掉,做 FFT 得到的就是阶次谱,横轴单位是「阶」而不是 Hz。常见做法是每转采 256 或 512 点,兼顾分辨率和计算量。
2.2 角度域重采样的两种实现路径
工程上落地角度重采样,主流有两条路。一条是硬件阶次跟踪,用专门的阶次分析仪或采集卡,靠锁相环硬件直接输出等角度采样时钟,精度高但设备贵。另一条是软件阶次跟踪,也就是这份代码走的路:先采时间域信号和键相脉冲,再在 MATLAB 里做插值重采样。软件方案灵活、成本低,代价是对键相脉冲的质量和插值算法有要求。
软件重采样的关键步骤是:从键相脉冲提取每转的到达时刻 t_i,构造角度-时间映射 θ(t),再对振动信号 x(t) 做插值,得到等角度序列 x(θ)。插值方法常见有线性插值、三次样条插值、sinc 插值。线性插值快但高频误差大,三次样条在振动信号里用得最多,sinc 插值精度最高但计算慢。这份代码里默认给的是三次样条,参数可切换。
2.3 代码文件的结构与依赖
拿到这份资源,先别急着跑主脚本,花两分钟看清目录结构能省很多事。典型组织方式是:
| 文件/目录 | 作用 | 是否必须 |
|---|---|---|
| main_order_analysis.m | 主入口,串起整条链路 | 是 |
| gen_sim_signal.m | 生成仿真振动+键相信号 | 调试用 |
| extract_keyphasor.m | 键相脉冲时刻提取 | 是 |
| angle_resample.m | 角度域重采样核心 | 是 |
| order_spectrum.m | 阶次谱计算与绘图 | 是 |
| data/ | 存放实测或仿真数据 | 按需 |
| config.m | 采样率、每转点数等参数 | 是 |
依赖上,主要用到signal工具箱里的resample、spline、pwelch,以及基础绘图。如果你的 MATLAB 没装信号处理工具箱,resample会报错,替代方案是用interp1手写重采样,代码里留了注释开关。
提示:先跑
gen_sim_signal.m生成一段已知阶次的仿真数据,再跑主流程,能快速判断代码链路是否通。直接上实测数据,一旦键相提取出错,后面全错,排查成本高。
3. 把代码跑起来:键相提取、重采样与阶次谱三步落地
3.1 键相脉冲时刻提取:阈值法加去抖
键相通道通常是一个每转一个脉冲的方波或正弦,幅值随转速略有变化。提取时刻最稳的办法是找上升沿过阈值的时刻,再做去抖。代码里extract_keyphasor.m的核心逻辑如下:
function t_key = extract_keyphasor(kp, fs, thresh_ratio) % kp: 键相通道原始信号 % fs: 采样率 Hz % thresh_ratio: 阈值比例,默认 0.5,即幅值中点的 50% if nargin < 3 thresh_ratio = 0.5; end % 自适应阈值:取上下包络中点,抗幅值漂移 hi = prctile(kp, 95); lo = prctile(kp, 5); thresh = lo + (hi - lo) * thresh_ratio; % 找上升沿:前一点低于阈值,当前点高于阈值 above = kp > thresh; rise_idx = find(diff(above) == 1) + 1; % 去抖:相邻上升沿间隔小于最小转周期的一半则丢弃 min_gap = round(0.5 * fs / 50); % 假设最高 3000 rpm,转周期 20ms keep = [true; diff(rise_idx) > min_gap]; rise_idx = rise_idx(keep); % 线性插值到亚采样点精度 t_key = zeros(size(rise_idx)); for k = 1:length(rise_idx) i = rise_idx(k); if i > 1 && kp(i) ~= kp(i-1) frac = (thresh - kp(i-1)) / (kp(i) - kp(i-1)); t_key(k) = (i - 1 + frac) / fs; else t_key(k) = (i - 1) / fs; end end end逻辑说明:先用 5% 和 95% 分位数估上下包络,避免固定阈值在幅值漂移时漏检或误检;diff(above)==1抓上升沿;min_gap去抖,防止噪声在阈值附近反复穿越造成假脉冲;最后用线性插值把上升沿定位到亚采样点,提升时刻精度。参数上,thresh_ratio默认 0.5,如果键相是窄脉冲而非方波,可以调到 0.3 左右;min_gap里的 50 是假设最低转速 50 rpm,实际按你的工况改。
3.2 角度域重采样:构造角度-时间映射
拿到每转的键相时刻 t_key 后,角度-时间映射就确定了:第 i 个脉冲对应角度 2π(i-1)。等角度重采样的目标是生成一组均匀角度 θ_j = j·Δθ,反查对应时刻 t_j,再对振动信号插值。
function [x_ang, theta] = angle_resample(x, t, t_key, n_per_rev) % x: 振动信号 % t: 时间向量 % t_key: 每转键相时刻 % n_per_rev: 每转采样点数,如 256 n_rev = length(t_key) - 1; % 完整转数 theta = (0 : n_rev*n_per_rev - 1) * (2*pi / n_per_rev); % 构造角度-时间映射:每转内线性插值 t_ang = zeros(size(theta)); for r = 1:n_rev idx = (r-1)*n_per_rev + (1:n_per_rev); t_ang(idx) = linspace(t_key(r), t_key(r+1), n_per_rev + 1); t_ang(idx) = t_ang(idx(1:end-1)); end % 三次样条插值到等角度时刻 x_ang = interp1(t, x, t_ang, 'spline', 0); end逻辑说明:n_rev是完整转数,最后一段不完整的转丢弃,避免边界效应;theta是等角度序列;内层linspace假设每转内转速近似恒定,这是软件阶次跟踪的常见近似,转速变化剧烈时误差会变大,后面避坑章会讲怎么判断;interp1用spline做三次样条插值,外插值填 0 防止 NaN。参数n_per_rev决定阶次谱的上限,按采样定理,能分析的最高阶次是 n_per_rev/2,256 点对应 128 阶,一般够用。
3.3 阶次谱计算:加窗、FFT 与阶次轴对齐
重采样后的 x_ang 是等角度序列,直接做 FFT 就得到阶次谱。但直接 FFT 会有泄漏,需要加窗。代码里order_spectrum.m的做法:
function [order, amp] = order_spectrum(x_ang, n_per_rev, win_type) % x_ang: 等角度重采样信号 % n_per_rev: 每转点数 % win_type: 'hann' 或 'rect' if nargin < 3 win_type = 'hann'; end N = length(x_ang); switch win_type case 'hann' w = hann(N); otherwise w = ones(N, 1); end X = fft(x_ang(:) .* w); P = abs(X(1:floor(N/2)+1)) / sum(w) * 2; order = (0:floor(N/2))' / n_per_rev; amp = P; end逻辑说明:加 Hann 窗抑制泄漏,sum(w)归一化保证幅值可读;取单边谱乘 2 补回负频率能量;阶次轴用(0:N/2)/n_per_rev换算,第 k 根谱线对应 k/n_per_rev 阶。注意这里阶次分辨率是 1/n_per_rev,256 点就是 1/256 阶,想分辨密集边带可以把 n_per_rev 提到 512 或 1024,代价是计算量和内存。
注意:阶次谱的幅值只有在整周期采样且窗函数归一化正确时才可信。如果只关心阶次位置不关心幅值,窗类型影响不大;要做趋势监测,窗和归一化必须固定,否则今天 0.1 明天 0.12 你会以为是故障发展。
4. 避坑与排查:阶次分析里最容易翻车的五个地方
4.1 键相脉冲漏检,阶次谱整体偏移
现象:阶次谱上 1 阶不在 1.0,而是 0.98 或 1.03,所有阶次按比例缩放。原因:键相提取时漏掉了个别脉冲,或者多检了噪声毛刺,导致每转对应的角度基准错位。解决:把extract_keyphasor输出的 t_key 画出来,看相邻间隔是否均匀;用diff(t_key)的直方图判断,正常应该集中在转周期附近,出现双峰或长尾就是漏检或多检。调thresh_ratio和min_gap重新提取。
4.2 转速变化剧烈时,每转内恒速假设失效
现象:升降速很快的工况下,阶次谱上高阶次峰值展宽、幅值偏低。原因:angle_resample里每转内用linspace假设转速恒定,实际转速在转内也有变化,角度-时间映射不准。解决:把每转内的线性插值换成基于键相脉冲的样条插值,或者提高键相脉冲密度(每转多个脉冲)。常见做法是在键相盘上做多齿,每转 60 个脉冲,映射精度能提升一个量级。
4.3 重采样后没做抗混叠,高阶次出现假峰
现象:阶次谱在接近 n_per_rev/2 的位置出现不明峰值,换数据依然存在。原因:时间域信号里高于 n_per_rev/2 阶的频率成分,在重采样时混叠进来。解决:重采样前先对时间域信号做低通滤波,截止频率对应最高分析阶次。代码里可以在angle_resample前加一步lowpass或filtfilt,截止频率取n_per_rev/2 * min(转频)。
4.4 数据段太短,阶次谱分辨率不够
现象:想分辨 1 阶和 1.02 阶的边带,谱上却是一根粗线。原因:阶次分辨率 = 1/n_per_rev,但实际分辨率还受数据长度限制,总转数太少时,即使 n_per_rev 很大,谱线也是插值出来的假分辨率。解决:保证至少 20 转以上的数据,做趋势分析时用 50 转以上。如果数据不够,宁可降低 n_per_rev 换取更稳的平均。
4.5 仿真数据能跑,实测数据报错或全零
现象:gen_sim_signal跑得好好的,换实测数据后interp1报错或结果全零。原因:实测数据的键相通道可能是负脉冲、或者振动通道有直流偏置、或者时间向量不是从 0 开始。解决:先检查t_key是否单调递增且在t范围内;振动信号去均值x = x - mean(x);时间向量统一从 0 开始。这三步做完,九成实测数据能跑通。
5. 进阶技巧:用阶次谱做故障定位与趋势监测
跑通基础流程之后,真正有价值的是把阶次谱用起来。我一般会做两件事:一是阶次切片,二是阶次趋势。
阶次切片是盯住某几个关键阶次,看它们随转速或时间的变化。比如齿轮箱,啮合阶次是齿数,边带是啮合阶次 ± 转频阶次,边带幅值上升往往指向齿轮裂纹或磨损。代码里可以在order_spectrum外面套一层循环,对每段数据算阶次谱,再提取目标阶次附近的峰值:
target_orders = [1, 2, 0.5, 18, 19]; % 转频、2倍频、半阶、啮合及边带 order_trend = zeros(length(target_orders), n_seg); for s = 1:n_seg seg = x_ang((s-1)*seg_len + (1:seg_len)); [o, a] = order_spectrum(seg, n_per_rev, 'hann'); for k = 1:length(target_orders) [~, idx] = min(abs(o - target_orders(k))); order_trend(k, s) = max(a(max(1,idx-2):min(end,idx+2))); end end这段逻辑是分段算阶次谱,再在目标阶次附近取局部最大值,避免单根谱线受噪声影响。seg_len一般取 10 到 20 转对应的点数,太短谱不稳,太长时间分辨率差。target_orders按你的设备特征频率填,不确定就先看全阶次谱找峰值。
阶次趋势的判读有个血泪经验:先看 1 阶稳不稳。1 阶是转速基准,如果 1 阶幅值自己就在漂,说明要么转速测量有问题,要么重采样参数不对,这时候看高阶次趋势全是玄学。1 阶稳了,再看高阶次和边带,才有意义。
另一个技巧是阶次谱与时间谱对照。同一段数据,分别做时间域 FFT 和阶次谱,如果某个峰值在时间谱上随转速跑、在阶次谱上固定,那它就是阶次相关成分,大概率是旋转机械故障;如果它在两边都固定,那可能是电气干扰或结构共振,跟转速无关。这个对照能帮你快速排除一批假故障。
验证方法上,我习惯用仿真数据做闭环:gen_sim_signal里预设 1 阶、2 阶、0.5 阶和 18 阶成分,跑完阶次分析看这四个阶次是否准确出现在对应位置、幅值比例是否对得上。对得上,再上实测数据。从那以后我每次换新设备数据,都强制先跑一遍仿真闭环,确认代码链路没被参数改动带偏,再动真数据。希望帮到你。
本文还有配套的精品资源,点击获取