★ FEATURED ARTICLE 简介旋转机械振动分析与故障诊断中阶次分析能解决变转速信号频率混叠问题将时域信号转换为角域阶次谱从而分离不同转速下的故障特征频率帮助识别齿轮、轴承等部件异常。这份MATLAB实现的阶次分析代码包面向振动监测、故障诊断和相关方向初学者可直接运行参考。压缩包共3个文件整体大小约13.72MB包含一个m脚本、一个mat数据文件和一张示意图片。其中m脚本是核心算法实现覆盖信号去噪、阶次跟踪与角域重采样等步骤mat文件可载入验证算法示意图便于理解流程。已有357人学习下载该资源具备一定实用价值。使用者可获得一套可运行、可修改的代码模板从原始信号到阶次谱的完整链路均可跑通也能替换自身数据调整参数用于教学、实验或工程故障特征提取为后续诊断分析提供基础。1. 阶次分析要解决的变转速工况下FFT为什么会失灵处理旋转机械振动时稳态工况下的FFT频谱非常直观峰值频率和故障特征一一对应。可一旦遇到升降速、启停机这类变转速工况同一台设备、同一个故障频谱里的峰值会跟着转速来回“漂移”用固定频率轴去解释越看越糊涂。阶次分析解决的就是这个问题把振动信号从时域重采样到旋转角度域用“转速的倍数”也就是阶次来代替绝对频率。齿轮啮合阶次、轴承故障特征阶次在变转速下仍然稳定这才是诊断能进行下去的前提。这套MATLAB实现从转速提取、等角度重采样到阶比谱和瀑布图输出覆盖了完整链路适合旋转机械升降速工况分析与故障诊断场景。2. 转速提取与等角度重采样整套代码的地基就这两块2.1 阶次分析的基本逻辑为什么坐标轴要换成旋转角度先回到原理层面想清楚一件事旋转机械的振动激励源大多和转角强相关。齿轮每转一圈啮合冲击发生的次数是固定的滚动轴承某个缺陷每转一圈也只在特定角度区间内被激励一次。也就是说振动信号在“角度域”里具有周期性而在“时间域”里只有当转速恒定时才呈现周期性。一旦转速变化同一故障对应的时间间隔忽长忽短FFT自然把能量摊开峰值变得又宽又矮。阶次分析的基本做法是建立瞬时转角与时间的映射关系转角为时间的积分θ(t) ∫ω(t)dt其中 ω(t) 是瞬时角速度在每个相等的转角间隔内重新采样振动信号得到角度域信号对角度域信号做FFT横轴单位就变成“每转的循环数”也就是阶次一个常被忽略的约束是阶次分析的精度上限取决于每转采样点数。若每转重采样1024点能分析的极限阶次大约是512阶奈奎斯特条件。工程中常见的高频齿轮啮合阶次、滚动轴承高频特征基本都被这个上限罩住先估算再选采样参数比事后返工省事得多。在这个MATLAB代码包里我习惯把实现拆成几个独立函数extract_rpm.m负责从转速脉冲信号提取转速曲线order_resample.m负责等角度重采样order_spectrum.m负责阶比谱计算order_waterfall.m负责瀑布图绘制。主脚本只需要按顺序调用调试时也能单独检查每一段输出的中间变量不会整个项目变成黑匣子。2.2 从测速脉冲到转速曲线一个干净利落的提取函数多数测试台架会给一路转速计脉冲信号tacho比如每转一个方波脉冲。转速提取的核心就是检测这些脉冲的到达时刻相邻两个脉冲的时间间隔对应一圈转速直接由 60/周期 算出来。MATLAB里最顺手的是用findpeaks定位脉冲位置再差分得到周期序列。function rpm_curve extract_rpm(tacho, fs, t) % 从转速计脉冲信号提取转速曲线 % 输入: % tacho - 转速计脉冲信号单位V % fs - 采样率单位Hz % t - 时间轴向量单位s % 输出: % rpm_curve - 两列矩阵第1列为时刻第2列为转速(rpm) % 检测脉冲上升沿位置MinPeakHeight按实际幅值一半左右设置 [~, locs] findpeaks(tacho, MinPeakHeight, 2.5, ... MinPeakDistance, round(fs/5000)); % 相邻脉冲间隔对应一转 pulse_t t(locs); rot_period diff(pulse_t); rpm_instant 60 ./ rot_period; % 中值滤波去掉单点毛刺 rpm_med medfilt1(rpm_instant, 5); % 转速值对应两个脉冲的中点时刻 t_mid (pulse_t(1:end-1) pulse_t(2:end)) ./ 2; rpm_curve [t_mid(:), rpm_med(:)]; end逻辑说明findpeaks返回的locs是脉冲峰值在原始信号中的索引转成时间后做一阶差分。rot_period里的每个元素是“这一圈用了多少秒”取倒数再乘60就是这一圈的平均转速。medfilt1的窗宽设为5对单点毛刺有很好的抑制效果又不会把真实的转速变化抹平。参数说明MinPeakHeight建议设为脉冲幅值的50%以上——如果脉冲是5V方波设在2.5V比较稳妥MinPeakDistance按最高转速下的脉冲间隔来估算比如最高转速6000rpm时脉冲间隔是0.01s采样率25600Hz下对应256个采样点设成round(fs/5000)取了约5个点这里有个细节fs/5000在25600Hz下是5.12取整是5明显太小。正确估算应该是fs/(max_rpm/60)即最高转速对应的采样点数。代码里这个参数如果设错会把同一个脉冲的多次振荡当成多个脉冲。实际使用时干脆设成fs*60/(max_rpm*0.9)留出冗余更安全。2.3 等角度重采样的完整实现插值边界必须处理干净有了转速曲线下一步是把时间域的振动信号映射到等角度网格上。常见做法是通过积分得到累计转角再在累计转角-时间曲线上插值找出等角度网格对应的时间点最后在这些时间点上对原始振动信号插值。这里每一步都要留意数据边界。function y_angle order_resample(data, fs, rpm_curve, spp) % 等角度重采样 % 输入: % data - 原始振动信号列向量 % fs - 采样率单位Hz % rpm_curve - n×2矩阵第1列为时间(s)第2列为转速(rpm) % spp - 每转采样点数建议1024或2048 % 输出: % y_angle - 角度域重采样后的信号 if nargin 4, spp 1024; end t_rpm rpm_curve(:,1); omega rpm_curve(:,2) / 60 * 2 * pi; % 瞬时角速度 rad/s theta cumtrapz(t_rpm, omega); % 累计转角 rad % 构造等角度网格 delta_theta 2 * pi / spp; theta_max floor(theta(end) / delta_theta) * delta_theta; theta_grid (0:delta_theta:theta_max); % 阶段1插值由累计转角反查时间点 t_interp interp1(theta, t_rpm, theta_grid, pchip); % 阶段2插值由时间点取值振动信号 y_angle interp1((0:length(data)-1)/fs, data, t_interp, linear); end逻辑说明cumtrapz做的是梯形积分转速曲线本身就是离散点用梯形积分比简单矩形积分更稳。theta_max取整是为了让末端不产生超出数据范围的插值点否则interp1会返回NaN后面FFT就全乱套了。插值分两步走第一步把角度网格转换成时间网格第二步按时间网格从原始信号取值。两段插值用不同方法是因为转速曲线平滑用pchip能避免spline在端点处的过冲振动信号本身是宽带信号用linear反而是最稳妥的避免插值自身引入虚假振荡。参数说明spp直接决定角度域采样密度一般取2的幂便于FFT。1024点对应理论最高512阶已覆盖绝大多数旋转机械特征如果关注的是高频轴承故障阶次比如300阶以上建议直接上2048。要注意的是interp1在第二步里若t_interp接近信号末端线性插值可能拿不到足够平滑的值实际工程中我会把末尾0.5%的数据直接丢弃避免边界效应污染阶比谱。3. 阶比谱、瀑布图与特征阶次提取让结果变成能用的诊断结论3.1 阶比谱计算对角度域信号做FFT坐标轴换成阶次角度域信号拿到后阶比谱的计算本质上就是一次普通FFT。区别只在横轴单位角度域信号的采样率是“每圈多少点”那么FFT后的频率轴单位就是“每圈多少次循环”即阶次。这部分最容易被半路出家的代码坑掉——直接用fft之后拿原始时间频率Hz当横轴那跟没做阶次分析没有区别。function [order_axis, spec_amp] order_spectrum(y_angle, spp, n_window) % 计算阶比谱 % 输入: % y_angle - 等角度重采样后的信号 % spp - 每转采样点数 % n_window - FFT点数可选默认取信号长度 % 输出: % order_axis - 阶次轴 % spec_amp - 单边幅值谱 if nargin 3, n_window length(y_angle); end n_window min(n_window, length(y_angle)); % 加Hanning窗抑制频谱泄漏 win hanning(n_window, periodic); y_seg y_angle(1:n_window) .* win; Y fft(y_seg); Y Y(1:floor(n_window/2)1); spec_amp 2 * abs(Y) / n_window; spec_amp(1) spec_amp(1) / 2; % 直流分量不乘2 % 阶次分辨率 1 / 参与FFT的总转数 total_rev n_window / spp; order_axis (0:floor(n_window/2)) / total_rev; end逻辑说明汉宁窗是阶比谱的默认配置因为升降速工况下信号本身就有较大的非平稳成分矩形窗的频谱泄漏会把相邻阶次的能量混成一团。spec_amp乘2是单边谱惯例直流分量单独处理否则0阶处会出现一个虚高尖峰。阶次轴的计算是关键FFT点数除以每转采样点数得到的total_rev是参与FFT的旋转圈数频率轴每个bin对应1/total_rev阶。比如信号总长是2048个角度点spp1024那这段信号只有2圈阶次分辨率就是0.5阶——这就是为什么想分辨0.1阶的细微差别必须让参与FFT的转数足够多。参数说明n_window默认取信号全长但如果总圈数太多比如持续升速几分钟直接整段FFT会丢掉转速变化中的细节信息。工程上广泛采用分段FFT的做法每段覆盖固定的转数比如8圈或16圈分段之间有50%重叠最后得到一系列阶比谱切片。这样做既保证阶次分辨率又能保留时间维度上的变化。3.2 阶比瀑布图把升降速过程按时间切片铺开阶比瀑布图是变转速分析中最直观的输出——横轴是阶次纵轴是时间或转速颜色深浅代表幅值。实现思路不复杂把重采样后的角度域信号按固定转数分段每段计算阶比谱拼成一个矩阵再用waterfall或pcolor绘图。function order_waterfall(y_angle, spp, rev_per_seg, overlap) % 绘制阶比瀑布图 % 输入: % y_angle - 等角度重采样信号 % spp - 每转采样点数 % rev_per_seg - 每段覆盖的转数推荐8~16 % overlap - 段重叠率0~0.75 n_seg floor(length(y_angle) / (spp * rev_per_seg * (1-overlap))) - 1; seg_len spp * rev_per_seg; spec_matrix zeros(floor(seg_len/2)1, n_seg); order_axis (0:floor(seg_len/2)) / rev_per_seg; for i 1:n_seg idx_start round((i-1) * seg_len * (1-overlap)) 1; idx_end idx_start seg_len - 1; if idx_end length(y_angle), break; end [~, spec_amp] order_spectrum(y_angle(idx_start:idx_end), spp, seg_len); spec_matrix(:, i) spec_amp; end % 用pcolor显示避免waterfall的线条遮挡 figure; pcolor(1:n_seg, order_axis, spec_matrix); shading interp; xlabel(时间切片序号); ylabel(阶次); colorbar; end逻辑说明分段步长等于段长乘以(1-overlap)保证相邻段之间有重叠谱图随时间连续变化。每段独立做窗函数和FFT幅值归一化方式一致拼出来的矩阵可以直接插值着色。用pcolor而不是waterfall是为了避免三维线条把峰值区域遮挡住——实际看谱时颜色深浅比线条高低更容易读。参数说明rev_per_seg是瀑布图最重要的参数。段覆盖转数越少时间分辨率越高但阶次分辨率越低二者由1/rev_per_seg直接绑定。比如每段8转阶次分辨率是0.125阶足够区分常见的齿轮啮合阶次边带如果每段只有2转分辨率只有0.5阶两条接近的阶次就分不开了。我的习惯是先用16转粗看全局再针对异常转速区间用8转甚至4转局部加密分析。3.3 特征阶次自动提取用findpeaks代替肉眼找峰瀑布图看完了最终要落实到“哪个阶次有异常能量”。手动缩放视图找峰效率太低尤其是数据段很多的时候。一个常用的自动提取方案是在每个阶比谱切片上用findpeaks定位局部极大值再按阶次位置做累计。function [peak_orders, peak_amps] extract_order_peaks(order_axis, spec_amp, ... order_range, min_prom) % 在指定阶次范围内搜索峰值 % 输入: % order_axis - 阶次轴 % spec_amp - 阶比幅值谱 % order_range - 搜索范围如[0 100] % min_prom - 最小突出度幅值低于该值的峰忽略 % 输出: % peak_orders - 峰值阶次位置 % peak_amps - 峰值幅值 mask order_axis order_range(1) order_axis order_range(2); [peak_amps, locs] findpeaks(spec_amp(mask), ... MinPeakProminence, min_prom, ... MinPeakDistance, 3); peak_orders order_axis(mask); peak_orders peak_orders(locs); end逻辑说明MinPeakProminence比MinPeakHeight更实用——它衡量的是峰值相对于两侧谷底的突出程度而不只是绝对大小。在不同转速段幅值差异很大时固定高度阈值会误判突出度阈值相对稳定。MinPeakDistance设为3避免同一个真实阶次的多个邻近点被重复识别。参数说明min_prom需要根据数据幅值定没有通用标准。我一般先画出一两个典型切片的阶比谱观察噪声基底幅值把阈值设在噪声基底的2~3倍。比如噪声基底是0.02阈值取0.05比较合理。这个值是“玄学”成分最多的参数需要针对每批数据调整但比固定高度阈值可靠得多。4. 避坑阶次分析里见过最多的五个参数与数据陷阱4.1 现象1转速脉冲毛刺导致重采样后信号呈锯齿状现象做完等角度重采样画出时域波形发现一段段的锯齿阶比谱在1阶、2阶附近出现很宽的裙边峰值高度明显下降。原因转速脉冲信号上有毛刺findpeaks把毛刺也识别成了脉冲导致某一段的脉冲间隔被错误缩短算出的转速瞬间飙高累计转角曲线出现局部抖动重采样时间点也跟着乱掉。解决先对tacho信号做一次低通滤波再用滞后比较施密特触发恢复方波最后才进findpeaks。我后来习惯在extract_rpm里直接加一道中值滤波窗宽7~9再用diff检查脉冲间隔把明显小于正常间隔90%的脉冲对丢弃效果很直接。如果毛刺密集建议回到原始信号看波形确认是不是测速传感器接触不良。4.2 现象2插值方法用错阶比谱出现虚假峰值现象角度重采样后的信号看起来平滑但阶比谱在高阶区域出现一组等间隔的虚假峰像是某种调制边带回头检查原始振动信号却没有这类特征。原因interp1用了spline插值。转速变化剧烈时三次样条在陡峭区域会产生过冲振荡这些振荡被当作真实振动信号进入FFT形成周期性假峰。spline的特点是连续且光滑但保形性差振动信号本身是宽带信号任何插值引入的人为振荡都会污染全频段。解决第一步插值角变时用pchip它能保持单调性不产生过冲第二步插值时域取振动值用linear线性插值不引入新频率分量。如果计算资源允许也可以改用resample函数配合去噪但核心原则不变插值方法选保守的别选花哨的。4.3 现象3转速曲线平滑窗口设太大阶次细节被抹掉现象阶比谱整体形状看着对但某些高阶级次比如50阶以上峰值比预期低很多边带结构也几乎看不见。原因转速曲线经过过度平滑瞬时转速的微小波动被平均掉。重采样时把这些波动当成“噪声”滤掉但真实机械的瞬时转速本身就包含扭振信息这些信息耦合在振动相位里平滑过头等于把信号源切断。解决平滑窗口的长度不要超过整个升速段数据点数的1%比如转速曲线有2万个点窗口设200点左右顶天。实际更推荐用movmedian而不是movmean中值对局部异常值稳健又能保留真实转速脉动。判断窗口是否合适的快速方法是把原始瞬时转速和平滑后的转速画在一张图上如果平滑曲线完全掩盖了脉动细节调小窗口。4.4 现象4记录时长不足阶次分辨率达不到需求现象分析目标是要区分35.2阶和35.5阶两个相邻特征但阶比谱里只有一个宽峰怎么调窗函数都分不开。原因阶次分辨率等于参与FFT的总转数的倒数。35.2和35.5只差0.3阶要分辨至少需要3.3转以上的数据参与FFT而如果整个数据段只有2转物理上限就摆在那里什么窗函数都救不了。解决这是原理性约束只能从数据端下手。一是延长采样时间保证目标转速区间覆盖足够多的转数二是把分析段集中在转速变化较慢的区间同等时间内转数更多三是用Zoom-FFT思路先对目标阶次附近窄带做带通滤波再降采样放大细节虽然分辨率仍受总转数限制但能改善局部谱的平整度。事前估算公式很简单所需总转数 1 / 阶次间隔设计测试时长时先按这个算一遍。4.5 现象5跳过抗混叠滤波高阶次区域出现镜像峰现象阶比谱在靠近理论最高阶次的区域出现成对的镜像峰值既有真实频率的成分又有往上折返后落到低阶区域的影子。原因等角度重采样的“采样率”是每转采样点数奈奎斯特上限是spp/2阶。原始时域信号里如果含有远超这个上限的高频成分重采样后会发生频谱混叠高阶次能量被折叠到其他位置。很多人只关注时域采样率是否够忽略了角度域同样有混叠风险。解决在重采样前对原始时域信号做低通滤波。截止频率按最高关注阶次乘以最高转速对应的旋转频率来算例如关注到200阶、最高转速6000rpm100Hz截止频率设为20kHz左右如果原始采样率只有25.6kHz那要么降低关注阶次上限要么提高每转采样点数。我在实际代码里会把这一步放在extract_rpm之后、order_resample之前先滤波再重采样顺序不能反。5. 计算阶次跟踪COT没有转速脉冲时的进阶替代方案5.1 为什么需要COT现场装不上转速计的场景太常见了不是所有设备都能方便地安装键相器或光电转速传感器。电机驱动端可能有现成的编码器但齿轮箱中间轴、负载端往往没有测速条件。还有一种情况是tacho信号在采集过程中丢失或严重污染重新采集一次代价太高。计算阶次跟踪COT的思路是用振动信号本身来估计参考轴的瞬时转速再走同样的等角度重采样流程。MATLAB实现COT并不是调用一个神奇函数而是自己完成“瞬时频率估计重采样”这两件事。5.2 从振动信号估计瞬时频率希尔伯特变换的实用做法最常见的做法是选择一个与转频强相关的窄带成分用希尔伯特变换提取瞬时相位再对相位差分得到瞬时频率。具体来说如果振动信号里有一个明显的1阶分量转频振动将它带通滤波后解析信号的相位就是转角的线性函数求导就能恢复转速。function rpm_est estimate_rpm_from_vib(data, fs, t) % 从振动信号估计瞬时转速用于无tacho场景 % 输入: % data - 振动信号 % fs - 采样率 % t - 时间轴 % 输出: % rpm_est - 两列矩阵第1列为时间第2列为转速 % 粗看频谱确定基频位置 N length(data); f_axis (0:floor(N/2)) * fs / N; Y abs(fft(data)); [~, idx] max(Y(2:floor(N/2)1)); f_peak f_axis(idx 1); % 带通滤波中心频率为峰值频率带宽取峰值的±15% bw f_peak * 0.15; data_f bandpass(data, [f_peak-bw f_peakbw], fs); % 希尔伯特变换提取瞬时相位 analytic hilbert(data_f); inst_phase unwrap(angle(analytic)); % 相位差分得到瞬时频率 inst_freq diff(inst_phase) / (2*pi) * fs; rpm_est [t(1:end-1), inst_freq * 60]; end逻辑说明先用整段频谱粗定位1阶基频峰再用窄带带通滤波把其他阶次成分滤掉避免它们干扰相位解调。hilbert得到解析信号后angle给出瞬时相位unwrap消除跳变差分就能恢复每时刻的瞬时频率。这个方法成立的前提是1阶振动分量存在且幅值明显高于噪声和其他阶次。如果设备在某个转速段1阶振动很弱估计出来的转速曲线会出现跳变需要额外做中值滤波。参数说明带宽系数0.15是经验值带宽太窄会丢失真实转速脉动太宽会引入相邻阶次干扰。如果1阶分量不突出也可以选某个齿轮啮合阶次比如25阶作为参考成分——阶次越高频率对转速变化越敏感但带通滤波的难度也越大。我在现场处理时一般先做一次短时傅里叶变换看哪个频率成分在整个升降速过程中连续且能量稳定再决定用哪个做COT源不盲目用1阶。5.3 COT的边界条件什么场景下别硬上COT看着方便但它有三个绕不开的短板选型前得想清楚。第一它依赖信噪比。如果振动信号里背景噪声和随机冲击较大带通滤波后仍无法获得干净的相位转速估计会频繁跳变。第二带通参数的设定直接影响结果同一个数据带宽取5%和取20%得到的转速曲线差异可能很大这是COT被叫作“玄学”的原因。第三COT估计的是平均意义上的转速无法捕捉转轴瞬间扭转振动带来的细节对于需要分析扭振的场合物理转速计仍然不可替代。工程上的判断标准很简单能用转速计的绝不用COT装不上转速计时优先用齿轮啮合的高阶次分量做参考而不是1阶因为高阶次信噪比通常更好如果数据采集规划阶段就知道要跑COT尽量提高采样率并延长记录时间给带通滤波留足空间。这些边界条件我会写在代码注释的头部提醒后面接手的人先做可行性判断再跑流程。6. 验证方法仿真合成信号先过一遍再上现场数据6.1 构造已知特征阶次的仿真信号阶次分析代码最容易出现的问题是“看起来在跑但结果谁也说不清对不对”。我的习惯是先造一组包含已知阶次成分的仿真信号把完整流程跑一遍确认输出的峰值位置和幅值和设定一致再换现场数据。fs 25600; t (0:fs*2-1) / fs; rpm 1500 (3000 - 1500) / 2 * t; % 线性升速 theta 2 * pi * cumtrapz(t, rpm / 60); % 累计转角 % 1阶、2阶、27.5阶三个特征分量加噪声 y 1.0 * sin(theta) ... 0.5 * sin(2 * theta) ... 0.2 * sin(27.5 * theta) ... 0.05 * randn(size(t));6.2 用输出结果反推参数合理性跑完阶次分析检查三个指标27.5阶峰值的实际位置偏差应小于阶次分辨率的1/21阶和2阶的幅值比应接近2:1因为原始幅值就是1和0.5瀑布图里27.5阶的谱线应当随转速变化保持水平而不是上下漂移。如果这三项都满足说明重采样链路没问题。之后换现场数据时如果结果异常回到第4章的五个排查点逐项核对而不是盲目调参数。从那以后我每次接到变转速振动分析任务都强制先跑一遍仿真合成信号验证模块确认代码链路本身没有引入额外噪声再换现场数据进同一流程。这套MATLAB实现不是只能处理理想数据——只要转速提取和重采样这两道关口把住了大多数工程数据都能给出稳定的阶比谱。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
