★ FEATURED ARTICLE 简介这份原创经典全息谱Matlab源程序面向从事旋转机械故障诊断、信号处理与状态监测方向的高校师生及工程技术人员核心解决全息谱holospectrum分析方法的编程实现问题。该方法由西安交大已故屈院士提出是转子振动全息诊断领域的经典理论适合具备一定Matlab基础、希望复现或二次开发全息诊断算法的读者。压缩包共39个文件约58KB以dv、dh系列数据文件为主配套7个m脚本文件另含htm使用说明、hlv与hlh等辅助文件数据与源码配套齐全便于直接运行验证。资源已有868人学习下载热度稳定。读者可从中获得全息谱二维与三维绘图的完整实现思路、窗函数处理与数据读取脚本以及配套的实测振动数据能够快速搭建全息诊断实验环境对照经典方法理解幅值谱、相位谱与全息谱之间的转换关系为论文复现或工程排故提供可运行的参考代码。1. 全息谱诊断的工程落地从一份 Matlab 源程序说起旋转机械的故障诊断里频谱分析是最常用的手段但它有个绕不开的短板——只处理单一通道的振动信号把转子系统各方向的振动割裂开来看。实际工况中转子在水平方向和垂直方向的振动是耦合的单纯看某一方向的频谱很容易漏掉一些只在两个通道之间才显现出来的特征。全息谱的思路就是把这层关系补上它把互相垂直的两个振动信号做融合处理在极坐标平面上还原出轴心轨迹的频域形态让不同频率成分的椭圆轨迹一目了然。这套方法在失衡、不对中、油膜涡动这类故障的识别上比传统单通道频谱多了一个维度的判断依据。这份 Matlab 源程序打包了一套全息谱诊断的完整实现从信号预处理、频谱变换到全息谱绘制和特征提取都有对应脚本。适合做旋转机械状态监测的工程师、带实验课的高校教师以及正在做故障诊断方向课题的研究生。如果你手头有双通道振动采集数据想快速验证全息谱方法能不能用起来这份代码能省掉从公式推导到编程实现的大段弯路。下面按实际拆包和跑通的顺序把关键环节和容易翻车的地方逐一讲清楚。2. 全息谱的数学骨架与代码结构拆解2.1 全息谱到底在算什么从两路信号到椭圆轨迹全息谱的核心操作对象是两个正交方向的振动信号通常记为 x(t) 和 y(t)。对这两个信号分别做 FFT得到各自的幅值和相位然后在每一个频率点上把 x 方向的幅值相位和 y 方向的幅值相位合成一个椭圆。这个椭圆的形状、倾角、长短轴比例就携带了该频率成分下的转子运动特征。具体来说对于某个频率 f设 x 方向在该频率的幅值为 Ax、相位为 φxy 方向为 Ay、相位为 φy则椭圆的长半轴、短半轴和倾角可以通过以下关系计算长半轴 R_major sqrt(Ax² Ay² 2·Ax·Ay·sin(Δφ)) 的某种组合形式短半轴 R_minor 对应 Ax² Ay² - 2·Ax·Ay·sin(Δφ)倾角 θ 由 Ax、Ay 和相位差 Δφ φy - φx 共同决定不同实现里公式的等价形式有差异但物理含义一致椭圆的扁率反映两个方向振动的相关程度正圆说明两个方向能量相当且相位差接近 90°退化成一条直线则说明振动集中在一个方向。工程上判断不对中时重点看 2X 频率处椭圆的形状是否变得扁平或出现明显的方向偏转。这份源程序里全息谱的合成逻辑集中在几个核心函数中输入是两路时域信号和采样率输出是各频率点对应的椭圆参数矩阵。理解了这个矩阵的维度含义后面调参和排错就有了抓手。2.2 源程序的文件组织与调用链拆开压缩包后代码大致分为三层数据层、算法层和展示层。数据层负责读取振动信号文件常见格式是文本或 mat 文件每列对应一个通道算法层包含 FFT 变换、全息谱合成、特征量计算等函数展示层负责绘制全息谱图、轴心轨迹图和瀑布图。典型的调用链是这样的% 主脚本全息谱诊断入口 clear; clc; close all; % 1. 加载双通道振动数据 % data 为 N×2 矩阵第一列 x 方向第二列 y 方向 % fs 为采样频率单位 Hz load(vibration_data.mat); % 假设文件内含 data 和 fs fs 2000; % 根据实际采集设置修改 % 2. 预处理去直流分量 data detrend(data, constant); % 3. 计算全息谱 % holospec 返回结构体包含频率向量、椭圆长轴、短轴、倾角 holospec compute_holospectrum(data(:,1), data(:,2), fs); % 4. 绘制全息谱图 plot_holospectrum(holospec); % 5. 提取特征频率处的椭圆参数 freq_of_interest 50; % 关注 1X 频率假设转频为 50Hz idx find_closest_freq(holospec.freq, freq_of_interest); fprintf(1X 处长轴: %.4f, 短轴: %.4f, 倾角: %.2f 度\n, ... holospec.major_axis(idx), holospec.minor_axis(idx), ... holospec.tilt_angle(idx));这段主脚本的逻辑很直白先加载数据去直流然后调用全息谱计算函数最后画图和提取特征。参数方面fs必须和采集时的实际采样率一致否则频率轴会整体偏移freq_of_interest根据你的转子转频来设比如转频 3000 rpm 对应 50 Hz那么 1X 就是 50 Hz2X 是 100 Hz。compute_holospectrum内部一般会做加窗处理常见的是 Hanning 窗窗长和重叠率会影响谱线的平滑程度。如果代码里没有自动加窗建议在 FFT 之前手动加窗否则泄漏会比较明显。2.3 关键参数采样率、窗函数与频率分辨率全息谱的质量很大程度上取决于三个参数采样率、窗函数类型和 FFT 点数。采样率决定了你能看到的最高频率。根据奈奎斯特采样定理可分析的最高频率是采样率的一半。如果转子系统的高频故障特征在 500 Hz 以上采样率至少要到 2000 Hz 才稳妥。这份代码里默认的采样率可能是一个常见值但实际使用时一定要按你的采集设备设置来改。窗函数的选择直接影响谱线的幅值精度和频率分辨率。Hanning 窗适合大多数旋转机械信号旁瓣衰减快泄漏小如果信号里包含明显的冲击成分Flat Top 窗在幅值精度上更有优势但频率分辨率会下降。代码里如果提供了窗函数选项建议先用 Hanning 窗跑一遍观察谱线是否干净。FFT 点数决定了频率分辨率即相邻两个频率点的间隔。分辨率 采样率 / FFT 点数。比如采样率 2000 HzFFT 点数 2048分辨率约为 0.98 Hz。如果两个故障特征频率靠得很近比如转频和它的边带就需要提高 FFT 点数来分辨。但点数太高会导致计算量增加而且如果数据长度不够补零虽然能加密频率轴但不会真正提高分辨率。提示改 FFT 点数时确保数据长度至少等于 FFT 点数否则补零后的谱线只是插值不能提升实际分辨能力。3. 跑通第一张全息谱图数据准备与参数调试3.1 数据格式对齐两路信号的时间同步与幅值标定全息谱对数据的要求比单通道频谱更严格因为两路信号必须严格时间同步。如果 x 和 y 方向的采集不是同一时刻触发的相位差就会引入额外误差椭圆的形状会失真。常见做法是用同一块采集卡的两个通道同时采样或者用同步触发信号确保两路数据起点一致。数据加载后先检查两路信号的均值和量纲。如果一路是电压信号V另一路是加速度信号g直接做全息谱没有物理意义。需要先统一到同一物理量比如都转成位移μm或都转成加速度m/s²。代码里如果有标定系数输入务必填对。% 检查两路信号的均值和标准差 fprintf(x 通道均值: %.6f, 标准差: %.6f\n, mean(data(:,1)), std(data(:,1))); fprintf(y 通道均值: %.6f, 标准差: %.6f\n, mean(data(:,2)), std(data(:,2))); % 如果量纲不一致手动标定 % 假设 x 通道灵敏度为 100 mV/gy 通道为 10 mV/μm % 需要根据实际传感器参数换算到同一物理量 sensitivity_x 100; % mV/g sensitivity_y 10; % mV/μm % 换算示例将电压转成加速度和位移再统一这段检查代码看起来简单但实际项目中经常因为忽略量纲一致性导致全息谱图完全不可读。我一般会在加载数据后先画两路信号的时域波形肉眼确认没有明显的幅值量级差异。3.2 第一次运行从默认参数到可读谱图第一次跑代码时建议先用默认参数把图出出来再说。如果代码里没有示例数据可以自己构造一段模拟信号来验证流程是否跑通% 构造模拟双通道振动信号 fs 2000; % 采样率 2000 Hz T 1; % 信号时长 1 秒 t (0:1/fs:T-1/fs); % 时间向量 f1 50; % 转频 50 Hz f2 100; % 2X 频率 100 Hz % x 方向1X 幅值 12X 幅值 0.3 x 1.0*sin(2*pi*f1*t) 0.3*sin(2*pi*f2*t pi/6); % y 方向1X 幅值 0.8相位差 45 度2X 幅值 0.5 y 0.8*sin(2*pi*f1*t pi/4) 0.5*sin(2*pi*f2*t pi/3); data [x, y];用这段模拟数据跑全息谱预期在 50 Hz 处看到一个椭圆在 100 Hz 处看到另一个椭圆。如果出来的图里椭圆退化成直线或者形状明显不对先检查相位差的计算是否正确。模拟信号里 x 和 y 的 1X 相位差是 45°全息谱图上椭圆的倾角应该和这个相位差有对应关系。3.3 谱图解读椭圆形状与故障类型的对应关系全息谱图出来后怎么读是关键。工程上常见的对应关系有这么几条1X 处椭圆接近正圆且长轴方向稳定通常说明转子存在不平衡且两个方向的振动能量相当。1X 处椭圆扁平长轴方向偏向某一侧可能是不对中或支撑刚度不对称。2X 处椭圆明显增大不对中的典型特征尤其是联轴器不对中。出现 0.4X0.5X 的低频椭圆且椭圆形状随转速变化油膜涡动或油膜振荡的嫌疑很大。椭圆的长短轴比例接近 1 但倾角频繁变化可能是松动或摩擦引起的非线性振动。这些对应关系不是绝对的需要结合转速、负载和历史数据综合判断。代码里如果有特征提取函数可以把这些椭圆参数批量算出来做成趋势图比单张谱图更有说服力。注意全息谱的椭圆参数对相位误差很敏感如果两路信号的相位差有系统性偏差椭圆形状会整体旋转或变形。建议先用已知相位差的标准信号校准一次。4. 避坑与排查全息谱代码常见的五个翻车点4.1 现象全息谱图上全是噪声看不到明显椭圆原因最常见的是两路信号没有做去直流处理或者去直流不彻底。直流分量在 FFT 后会变成 0 Hz 处的一个大峰值如果代码里没有把它剔除全息谱的低频区域会被这个峰值淹没。另一个原因是信号幅值太小量化噪声占了主导。解决在 FFT 之前强制去均值用detrend或直接减去mean。如果信号幅值确实小先做归一化再分析。代码里如果有remove_dc选项确保它是打开的。4.2 现象椭圆的长短轴比例和理论值对不上原因相位差的计算可能用了错误的参考方向。全息谱里 x 和 y 的相位差定义是 φy - φx如果代码里写反了椭圆的旋转方向会反过来长短轴比例也可能变化。另一个可能是 FFT 后的相位没有做 unwrap导致跳变。解决检查相位差的计算顺序确认和你的坐标系定义一致。对于相位跳变在计算相位差之前先做unwrap。如果代码里没有这个步骤手动加上。4.3 现象改采样率后频率轴整体偏移原因代码里可能硬编码了采样率或者频率轴的计算用了错误的公式。常见的是频率向量写成(0:N-1)*fs/N但如果fs没有同步更新频率轴就会错位。解决把采样率作为参数传入所有相关函数避免全局变量。频率向量的计算统一用(0:N/2)*fs/N取正半轴。改完参数后用已知频率的模拟信号验证一次。4.4 现象全息谱图只显示一半负频率部分丢失原因FFT 的结果是对称的正负频率各一半。全息谱通常只关心正频率部分但如果代码里直接用了完整的 FFT 结果而没有取正半轴绘图时可能出现镜像或重叠。解决在计算全息谱之前先取 FFT 结果的前 N/21 个点对应的频率向量也取正半轴。如果代码里已经做了这个处理检查索引是否从 1 开始Matlab 索引从 1 开始容易和 0 混淆。4.5 现象换一组数据后代码报错提示维度不匹配原因不同来源的数据文件格式可能不一样有的是一列一个通道有的是行向量有的包含表头。代码里如果写死了data(:,1)和data(:,2)遇到转置的数据就会出错。解决加载数据后先判断维度如果是行向量就转置。用size检查行列数确保是 N×2 的矩阵。如果数据文件有表头用readmatrix或importdata时跳过表头行。% 数据维度自适应处理 if size(data, 1) size(data, 2) data data; % 转置成 N×2 end if size(data, 2) 2 error(数据必须包含两个通道); end5. 进阶用法批量处理与特征趋势跟踪5.1 批量跑多组数据循环调用与结果汇总实际项目中往往有几十组甚至上百组振动数据需要分析。一组一组手动跑效率太低把主脚本改成批量处理模式更实际。思路是遍历数据文件夹对每个文件调用全息谱计算函数把关键频率处的椭圆参数存到一个结果矩阵里。% 批量处理文件夹下的所有 mat 文件 data_dir vibration_data/; file_list dir(fullfile(data_dir, *.mat)); results []; for i 1:length(file_list) file_path fullfile(data_dir, file_list(i).name); loaded load(file_path); data loaded.data; fs loaded.fs; % 去直流 data detrend(data, constant); % 计算全息谱 hs compute_holospectrum(data(:,1), data(:,2), fs); % 提取 1X 和 2X 处的椭圆参数 f1 50; % 转频根据实际工况修改 f2 100; % 2X idx1 find_closest_freq(hs.freq, f1); idx2 find_closest_freq(hs.freq, f2); results [results; ... f1, hs.major_axis(idx1), hs.minor_axis(idx1), hs.tilt_angle(idx1), ... f2, hs.major_axis(idx2), hs.minor_axis(idx2), hs.tilt_angle(idx2)]; end % 保存结果到 CSV csvwrite(holospectrum_results.csv, results);这段代码的关键在于find_closest_freq函数它负责在频率向量里找到最接近目标频率的索引。如果频率分辨率不够高找到的索引可能偏离实际频率导致椭圆参数不准。解决办法是提高 FFT 点数或者在目标频率附近做插值。5.2 特征趋势跟踪从单张谱图到状态监测单张全息谱图只能反映某一时刻的状态要做状态监测需要把多组数据的特征参数按时间顺序排列观察趋势变化。比如 1X 处长轴幅值随时间逐渐增大通常是不平衡发展的信号2X 处短轴比例突然下降可能是不对中加剧。把批量处理的结果矩阵按时间排序用plot画出各特征参数的趋势曲线。如果代码里有趋势分析函数直接调用没有的话几行代码就能补上% 假设 results 按时间顺序排列 time 1:size(results, 1); figure; subplot(2,1,1); plot(time, results(:,2), b-o, time, results(:,6), r-s); legend(1X 长轴, 2X 长轴); xlabel(数据组序号); ylabel(长轴幅值); title(全息谱特征趋势); subplot(2,1,2); plot(time, results(:,4), b-o, time, results(:,8), r-s); legend(1X 倾角, 2X 倾角); xlabel(数据组序号); ylabel(倾角 (度));趋势图比单张谱图更容易发现渐变型故障。我一般会在趋势图上加一条报警线当特征值超过阈值时自动标记出来。阈值怎么定需要结合历史数据和设备厂家的建议没有统一标准。5.3 一个容易忽略的细节转速波动对全息谱的影响全息谱假设分析窗口内转速是稳定的。如果转速波动较大1X 频率会在一个范围内漂移FFT 后的谱线会变宽椭圆的形状也会模糊。这种情况下直接做全息谱效果不好需要先做阶次分析或转速跟踪。常见做法是用转速脉冲信号做键相把时域信号重采样成等角度间隔的信号再做全息谱。这份代码里如果没有阶次分析模块可以先用转速稳定的数据段做分析避开启停阶段。如果必须分析变速工况建议先做角域重采样再调用全息谱函数。提示转速波动超过 ±2% 时全息谱的椭圆参数可信度会明显下降。分析前先看一眼转速信号的稳定性。从那以后我每次拿到一组新的振动数据都会先画时域波形和转速趋势确认数据段稳定、量纲一致、两路同步再跑全息谱。这个习惯帮我省掉了很多次返工。希望这份源程序和你手头的数据能顺利对上帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
