简介这份资源是原创经典全息谱Matlab源程序面向从事旋转机械故障诊断、信号处理与状态监测方向的高校师生及工程技术人员用于实现全息诊断holospectrum方法该方法由西安交大已故屈院士发明适合作为课程设计、毕业设计或科研复现的参考代码。压缩包共39个文件约58KB以dv、dh系列数据文件为主配合7个m脚本文件及dv1、dh1、dv2、dh2、h、hlv、hlh、v等辅助数据另含htm使用说明与url链接数据与脚本配套便于直接运行与验证。目前已有868人学习下载。资源包含holo2D.m、holo3D.m、dire.m、sorbit.m、gr.m、rg.m、hanning_window.m等核心脚本覆盖二维与三维全息谱绘制、方向判别、轴心轨迹合成及加窗处理等环节并附带多组实测数据可帮助读者理解全息谱的构造流程与诊断逻辑快速搭建实验环境并对照结果排查问题。1. 全息谱诊断从一段被丢进回收站的代码说起很多人第一次听到“全息谱”这三个字是在旋转机械故障诊断的论文里而不是在工程现场。它听起来像某种玄学实际上是一套把转子振动信号从“单通道频谱”升级成“多通道全息图”的分析方法。传统 FFT 只看一个测点、一个方向的幅值谱丢掉了相位而全息谱把同一轴承截面两个正交方向水平 X、垂直 Y的振动信号合起来用椭圆轨迹的方式表达每一阶频率分量幅值、相位、进动方向全在里面。这套方法对不平衡、不对中、油膜涡动、碰摩这些典型故障的区分度比单通道谱高出一截。标题里的“全息诊断”加“Matlab 源程序”指向的是一份能直接跑起来的工程实现读振动数据、做频谱校正、合成全息谱、画椭圆、输出诊断结论。适合两类人——做旋转机械状态监测的现场工程师以及需要复现全息谱算法做二次开发的研究生。下面我按“先立住原理、再动手复现、最后讲坑”的顺序把这份方案拆开讲清楚。2. 全息谱的数学骨架为什么必须做频谱校正2.1 从单通道 FFT 到正交合成单通道 FFT 给出的是某个方向上的幅值谱和相位谱但相位是相对于各自通道触发时刻的两个通道之间没有可比性。全息谱的核心操作是对 X、Y 两个通道分别做 FFT提取同一阶频率的幅值和相位然后合成一个椭圆。椭圆的半长轴、半短轴、倾角、进动方向就是这阶分量的全息特征。数学上设 X 通道在频率 f 处的复幅值为 A_x·e^(jφ_x)Y 通道为 A_y·e^(jφ_y)则合成信号可以写成x(t) A_x·cos(2πft φ_x) y(t) A_y·cos(2πft φ_y)消去时间 t得到椭圆方程。椭圆的形状参数直接对应故障类型正进动椭圆偏圆说明不平衡反进动分量大说明碰摩或油膜失稳椭圆很扁说明某一方向刚度不对称。2.2 为什么不做校正全息谱就是错的这是最容易翻车的地方。FFT 的离散频率分辨率是 fs/N如果转子工频不是恰好落在谱线上就会发生泄漏。泄漏导致幅值和相位同时失真而全息谱对相位极其敏感——相位差几度椭圆倾角和进动方向就可能反过来。常见做法是加 Hanning 窗做幅值修正再用比值法或能量重心法做频率校正。我一般用比值法找到峰值谱线 k 及其左右邻线 k-1、k1用它们的幅值比反推真实频率偏移量。校正后再重新计算该频率处的复幅值而不是直接取峰值谱线的值。% 比值法频率校正输入峰值谱线索引 k幅值谱 X_mag function [f_corr, A_corr, phi_corr] ratio_correct(X_mag, k, fs, N) % 取左右邻线幅值 y1 X_mag(k-1); y2 X_mag(k); y3 X_mag(k1); % 比值法求偏移量 delta范围 [-0.5, 0.5] delta (y3 - y1) / (y1 2*y2 y3) * 2; % 校正后的频率 f_corr (k - 1 delta) * fs / N; % 幅值修正用窗函数主瓣能量归一化 A_corr (y1 2*y2 y3) / 2 * (pi * delta / sin(pi * delta)); % 相位修正补偿偏移带来的相位旋转 phi_corr angle(X_mag(k)) - pi * delta; end这段代码里delta是真实频率相对于峰值谱线的偏移比例A_corr用三项幅值之和做能量归一化phi_corr补偿了偏移导致的相位旋转。参数fs是采样率N是 FFT 点数。注意delta的符号决定频率往哪边偏如果y1 y3说明真实频率在峰值左侧delta为负。2.3 全息谱的合成与椭圆参数提取校正完两个通道的复幅值后合成椭圆。设 X 通道复幅值为 XcY 通道为 Yc则椭圆的长半轴、短半轴、倾角可以用以下公式计算% 全息谱椭圆参数提取 function [a, b, theta, direction] holo_ellipse(Xc, Yc) % Xc, Yc 为复数分别代表 X、Y 通道在目标频率的复幅值 % 构造椭圆参数 A abs(Xc); B abs(Yc); phi_x angle(Xc); phi_y angle(Yc); % 相位差 dphi phi_y - phi_x; % 长半轴和短半轴 a sqrt((A^2 B^2 sqrt((A^2 - B^2)^2 4*A^2*B^2*cos(dphi)^2)) / 2); b sqrt((A^2 B^2 - sqrt((A^2 - B^2)^2 4*A^2*B^2*cos(dphi)^2)) / 2); % 倾角 theta 0.5 * atan2(2*A*B*cos(dphi), A^2 - B^2); % 进动方向dphi 在 (0, pi) 为正进动(-pi, 0) 为反进动 if sin(dphi) 0 direction 正进动; else direction 反进动; end enda和b分别是长半轴和短半轴theta是椭圆倾角direction由相位差的正弦符号决定。这里的关键是dphi必须用校正后的相位否则进动方向会判反。我见过不止一次现场数据没做校正全息谱画出来是反进动拆开一看是正进动白折腾一天。3. 用 Matlab 跑通全息谱诊断的最小闭环3.1 数据准备与预处理假设你手头有两通道振动数据采样率 2000 Hz采样点数 4096工频大约 50 Hz。数据格式可能是 CSV、MAT 或者二进制。我一般先做去直流和趋势项消除避免低频漂移污染谱线。% 读取两通道振动数据并预处理 load(vib_data.mat); % 假设变量名为 x_signal, y_signal fs 2000; % 采样率 N 4096; % FFT 点数 t (0:N-1) / fs; % 时间向量 % 去直流 x_signal x_signal - mean(x_signal); y_signal y_signal - mean(y_signal); % 去趋势项线性拟合后减去 x_signal detrend(x_signal); y_signal detrend(y_signal); % 加 Hanning 窗 win hanning(N); x_win x_signal(1:N) .* win; y_win y_signal(1:N) .* win;detrend去掉线性趋势hanning窗减少泄漏。注意窗函数会改变幅值后面校正时要补偿。如果数据长度不够 N要么补零要么降低 N但补零会改变谱线间隔校正时要注意。3.2 双通道 FFT 与频谱校正对两个通道分别做 FFT找到工频附近的峰值谱线然后调用前面的校正函数。% 双通道 FFT X_fft fft(x_win); Y_fft fft(y_win); X_mag abs(X_fft(1:N/21)); Y_mag abs(Y_fft(1:N/21)); f_axis (0:N/2) * fs / N; % 找工频附近峰值假设 45-55 Hz 范围 search_range find(f_axis 45 f_axis 55); [~, idx] max(X_mag(search_range)); k search_range(idx); % 峰值谱线索引 % 对 X 通道校正 [f_x, A_x, phi_x] ratio_correct(X_mag, k, fs, N); % 对 Y 通道校正同一频率处 [f_y, A_y, phi_y] ratio_correct(Y_mag, k, fs, N); % 合成复幅值 Xc A_x * exp(1j * phi_x); Yc A_y * exp(1j * phi_y);这里k是峰值谱线索引两个通道用同一个k保证频率一致。如果两个通道峰值不在同一根谱线说明系统有各向异性需要分别校正后再统一到同一频率。ratio_correct返回的f_x和f_y应该很接近如果差太多检查数据质量。3.3 全息谱绘制与诊断输出拿到Xc和Yc后调用椭圆参数提取函数然后画图。% 提取椭圆参数 [a, b, theta, direction] holo_ellipse(Xc, Yc); % 绘制全息谱椭圆 figure; t_ell linspace(0, 2*pi, 200); % 椭圆参数方程考虑倾角 x_ell a * cos(t_ell) * cos(theta) - b * sin(t_ell) * sin(theta); y_ell a * cos(t_ell) * sin(theta) b * sin(t_ell) * cos(theta); plot(x_ell, y_ell, b-, LineWidth, 1.5); axis equal; grid on; xlabel(X 方向位移); ylabel(Y 方向位移); title(sprintf(全息谱椭圆%s长半轴 %.3f短半轴 %.3f, direction, a, b)); % 诊断输出 fprintf(工频 %.2f Hz\n, f_x); fprintf(长半轴 %.4f短半轴 %.4f倾角 %.2f 度\n, a, b, rad2deg(theta)); fprintf(进动方向%s\n, direction); if b / a 0.3 strcmp(direction, 正进动) fprintf(诊断建议疑似不平衡建议做动平衡。\n); elseif b / a 0.7 fprintf(诊断建议椭圆接近圆形检查是否存在松动或碰摩。\n); elseif strcmp(direction, 反进动) fprintf(诊断建议存在反进动分量重点排查油膜涡动或碰摩。\n); end这段代码把椭圆画出来同时根据长短轴比和进动方向给出初步诊断。b/a 0.3且正进动典型不平衡b/a 0.7接近圆形可能是松动或碰摩反进动则指向油膜失稳。这些阈值不是绝对的现场要结合历史数据对比。4. 避坑指南全息谱诊断里最容易翻车的五个点4.1 相位没校正进动方向判反现象全息谱显示反进动但现场振动特征明明是不平衡。原因两个通道的相位没有做频谱校正或者校正时用了不同的频率偏移量。解决确保两个通道用同一个峰值谱线索引校正函数返回的相位要补偿偏移量。如果两个通道峰值不在同一谱线先统一频率再校正。4.2 窗函数用错幅值偏小一半现象全息谱椭圆比预期小很多长半轴只有理论值的一半。原因加了 Hanning 窗但没做幅值恢复。Hanning 窗的相干增益是 0.5幅值要乘 2。解决在ratio_correct里做能量归一化或者直接用A_corr A_corr * 2补偿。注意不同窗函数增益不同矩形窗是 1Hanning 是 0.5Hamming 是 0.54。4.3 采样率不够工频混叠现象全息谱在工频处没有明显峰值反而在低频出现异常椭圆。原因采样率低于工频的两倍或者抗混叠滤波器失效。解决检查采样率是否满足 Nyquist 定理旋转机械工频通常几十赫兹采样率至少 1000 Hz。如果数据已经混叠没法恢复只能重新采集。4.4 只取峰值谱线不做插值现象全息谱椭圆参数随数据长度变化很大重复性差。原因直接取峰值谱线的幅值和相位没有做频率插值。解决用比值法或能量重心法做校正校正后的参数对数据长度不敏感。我一般用比值法计算量小精度够用。4.5 忽略转速波动全息谱糊成一团现象全息谱椭圆边界模糊像一团乱麻。原因转子转速在采集期间波动工频不是恒定值。解决用键相脉冲做整周期采样或者用短时傅里叶变换做时频分析提取每一时刻的工频。如果转速波动超过 1%全息谱方法本身就不适用考虑阶次分析。5. 进阶技巧用全息谱做多阶分量叠加与趋势跟踪全息谱不止能看工频。把 2X、3X、0.5X 等分量的椭圆画在同一张图上可以直观看到各阶分量的相对大小和进动方向。我一般会做一个“全息谱瀑布图”横轴是频率纵轴是时间每个时间点画一个椭圆颜色表示进动方向。这样能跟踪故障的演化过程。% 多阶全息谱叠加提取 0.5X、1X、2X、3X 分量 harmonics [0.5, 1, 2, 3]; figure; hold on; colors lines(length(harmonics)); for i 1:length(harmonics) f_target harmonics(i) * f_x; k_target round(f_target * N / fs) 1; if k_target 1 k_target N/2 [~, A_x, phi_x] ratio_correct(X_mag, k_target, fs, N); [~, A_y, phi_y] ratio_correct(Y_mag, k_target, fs, N); Xc A_x * exp(1j * phi_x); Yc A_y * exp(1j * phi_y); [a, b, theta, ~] holo_ellipse(Xc, Yc); t_ell linspace(0, 2*pi, 100); x_ell a * cos(t_ell) * cos(theta) - b * sin(t_ell) * sin(theta); y_ell a * cos(t_ell) * sin(theta) b * sin(t_ell) * cos(theta); plot(x_ell, y_ell, Color, colors(i,:), LineWidth, 1.2); end end legend(0.5X, 1X, 2X, 3X); axis equal; grid on; title(多阶全息谱叠加);这段代码把 0.5X、1X、2X、3X 的椭圆画在一起。0.5X 分量大通常指向油膜涡动2X 分量大指向不对中3X 及以上指向松动或碰摩。注意k_target要取整如果目标频率不在谱线上校正函数会给出偏移量但偏移太大时精度下降。趋势跟踪的做法是把采集数据按时间分段每段做一次全息谱提取长半轴、短半轴、进动方向画成趋势曲线。长半轴突然增大说明该阶分量在恶化进动方向从正变反说明故障类型在转变。我一般会设一个报警阈值长半轴超过基线 2 倍或者进动方向翻转就触发预警。最后说一个血泪教训全息谱诊断的准确性八成取决于数据质量两成取决于算法。传感器安装角度偏了、键相脉冲抖了、采集期间转速飘了后面算法再精细也救不回来。我现在的习惯是拿到数据先看时域波形和转速脉冲确认没问题再进全息谱流程。希望帮到你。本文还有配套的精品资源点击获取