简介：这份资源是原创经典全息谱Matlab源程序，面向从事旋转机械故障诊断、信号处理与状态监测方向的高校师生及工程技术人员，用于实现全息诊断（holospectrum）方法，该方法由西安交大已故屈院士发明，适合作为课程设计、毕业设计或科研复现的参考代码。压缩包共39个文件，约58KB，以m脚本文件为核心，配合dv、dh、dv1、dh1、dv2、dh2等数据文件及h、hlv、hlh等辅助文件，另含htm使用说明与url链接，覆盖二维、三维全息谱绘制及窗函数处理等环节。资源已有868人学习下载，热度稳定。读者可从中获得完整的全息谱计算与可视化脚本、配套示例数据、窗函数实现及使用帮助文档，便于快速理解全息诊断原理、复现经典算法并迁移到自身振动信号分析任务中。
1. 全息谱诊断到底在诊断什么：从一根轴的振动说起
一台旋转机械的振动信号，本质上是一组随时间变化的波形。传统做法是把它做 FFT，看各个频率上的幅值大小，然后对照特征频率判断故障类型。这套流程用了几十年，够用，但有个盲区：它把相位信息丢掉了。两个频率成分幅值相同、相位差不同，在普通频谱上看起来一模一样，可转子系统里，相位差恰恰是区分不平衡、不对中、松动这些故障的关键线索。
全息谱（Holospectrum）要解决的就是这个问题。它把转子同一截面两个正交方向（比如水平和垂直）的振动信号同时采集，做 FFT 之后，把同一频率下两个方向的幅值和相位合成一个椭圆——这个椭圆就是该频率下的全息谱图。椭圆的形状、偏心率、主轴方向，直接对应转子在那一频率下的运动轨迹形态。不平衡表现为接近正圆的椭圆，不对中表现为偏心率很大的扁椭圆，油膜涡动则会出现特殊的进动方向。
这套方法适合谁？做旋转机械状态监测和故障诊断的工程师，手里有双通道以上的振动采集设备，能拿到同一截面两个正交方向的同步信号。如果你只有单通道数据，全息谱的核心优势发挥不出来，得先解决采集端的问题。Matlab 是实现全息谱诊断最顺手的工具，矩阵运算天然适合做 FFT 和椭圆拟合，画图也方便。下面从数据准备开始，一步步把全息谱诊断的完整链路走通。
2. 从双通道振动信号到全息谱椭圆：Matlab 实现全链路
2.1 数据准备：同一截面两个正交方向的同步采样
全息谱的输入数据有两个硬性要求：第一，两个通道必须是同一截面的正交方向，通常标为 X 和 Y；第二，两个通道必须严格同步采样，不能有时间偏移。如果采集设备支持同步采样，直接导出两列数据即可。如果只有分别采集的数据，需要先做时域对齐，这一步后面避坑章节会展开。
假设你已经拿到了两列数据，存成 CSV 或者 MAT 文件。我一般习惯先把数据整理成一个结构体，后面调用方便：
% 读取双通道振动数据 % 假设数据文件为两列：第一列X方向，第二列Y方向 % 采样率 fs = 20000 Hz，采样时长 2 秒 data = readmatrix('vibration_data.csv'); x = data(:, 1); % X方向振动信号 y = data(:, 2); % Y方向振动信号 fs = 20000; % 采样频率，单位Hz N = length(x); % 采样点数 t = (0:N-1) / fs; % 时间向量 % 去均值，消除直流分量 x = x - mean(x); y = y - mean(y); % 检查两通道数据长度是否一致 assert(length(x) == length(y), '两通道数据长度不一致，检查采集设置');这段代码做了三件事：读数据、去均值、校验长度。去均值这一步容易被忽略，但直流分量会在 FFT 的 0 Hz 处产生一个大峰值，虽然不影响其他频率的分析，但会让后续画图时的坐标范围变得很难看。assert那行是血泪经验——两通道长度不一致通常意味着采集时有一路丢了触发或者提前停止，这种数据直接扔掉重采，不要试图截齐。
2.2 FFT 与相位提取：为什么不能只用 abs(fft(x))
很多人做频谱分析习惯直接abs(fft(x))，只看幅值。全息谱必须保留相位，所以要用完整的复数结果。Matlab 的fft返回复数，实部是余弦分量，虚部是正弦分量，angle函数取出相位角。
% 对两通道分别做FFT X_fft = fft(x) / N; % 归一化，除以N使幅值反映真实物理量 Y_fft = fft(y) / N; % 取单边谱（正频率部分） half_N = floor(N/2) + 1; X_spec = X_fft(1:half_N); Y_spec = Y_fft(1:half_N); % 频率轴 f = (0:half_N-1) * fs / N; % 提取幅值和相位 X_amp = 2 * abs(X_spec); % 乘2恢复单边幅值 X_phase = angle(X_spec); % 相位，弧度 Y_amp = 2 * abs(Y_spec); Y_phase = angle(Y_spec); % 直流分量不乘2 X_amp(1) = X_amp(1) / 2; Y_amp(1) = Y_amp(1) / 2;这里有几个参数需要说清楚。除以N是归一化，让幅值反映信号的实际大小，不然 FFT 结果的幅值和采样点数相关，换个采样长度数值就变了。取单边谱是因为实信号的负频率部分是冗余的，只取正频率。乘 2 是为了补偿丢掉的那一半能量，但直流分量（索引 1）不能乘 2，因为它没有对称的负频率分量。
相位提取用的是angle，返回的是[-π, π]范围内的弧度值。全息谱合成椭圆时，两个方向的相位差才是关键，绝对相位值本身受触发时刻影响，但相位差是稳定的。
2.3 全息谱椭圆合成：从两个复数到一个椭圆
全息谱的核心操作在这里：对每一个感兴趣的频率点，取 X 和 Y 的幅值和相位，合成一个椭圆。椭圆的参数方程可以写成：
x(t) = X_amp * cos(2πft + X_phase) y(t) = Y_amp * cos(2πft + Y_phase)消去时间 t，得到椭圆的一般方程。但实际画图时不需要显式求椭圆方程，直接在一个周期内采样画出来就行：
% 选择要分析的目标频率（比如转频） target_freq = 50; % Hz，假设转频为50Hz [~, idx] = min(abs(f - target_freq)); % 找到最接近的频率索引 % 提取该频率下的幅值和相位 xa = X_amp(idx); xp = X_phase(idx); ya = Y_amp(idx); yp = Y_phase(idx); % 合成椭圆：在一个周期内采样200个点 theta = linspace(0, 2*pi, 200); x_ellipse = xa * cos(theta + xp); y_ellipse = ya * cos(theta + yp); % 绘制全息谱椭圆 figure; plot(x_ellipse, y_ellipse, 'b-', 'LineWidth', 1.5); hold on; plot(0, 0, 'k+', 'MarkerSize', 10); % 标出原点 axis equal; % 等比例坐标轴，椭圆形状才真实 xlabel('X方向位移'); ylabel('Y方向位移'); title(sprintf('全息谱椭圆 - %d Hz', target_freq)); grid on;axis equal这行必须加。不加的话 Matlab 会自动缩放坐标轴，一个正圆可能被拉成椭圆，一个扁椭圆可能看起来像正圆，形状判断完全失真。这是新手最容易翻车的地方之一。
2.4 批量生成全息谱图：从单频到全频段
实际诊断中不会只看一个频率，而是要把整个频段的全息谱都画出来。常见做法是选若干个幅值较大的频率点，每个画一个椭圆，按频率排列。也可以画成三维形式，Z 轴是频率：
% 选取幅值超过阈值的前若干个频率点 threshold = 0.05 * max(X_amp); % 阈值设为最大幅值的5% peak_idx = find(X_amp > threshold); % 限制显示数量，避免图太乱 max_peaks = 10; if length(peak_idx) > max_peaks [~, sort_idx] = sort(X_amp(peak_idx), 'descend'); peak_idx = peak_idx(sort_idx(1:max_peaks)); end % 逐个绘制全息谱椭圆 figure; for k = 1:length(peak_idx) i = peak_idx(k); xa = X_amp(i); xp = X_phase(i); ya = Y_amp(i); yp = Y_phase(i); theta = linspace(0, 2*pi, 200); xe = xa * cos(theta + xp); ye = ya * cos(theta + yp); subplot(2, 5, k); plot(xe, ye, 'b-', 'LineWidth', 1.2); hold on; plot(0, 0, 'k+'); axis equal; title(sprintf('%.1f Hz', f(i))); grid on; end阈值设为最大幅值的 5% 是个经验值。太低会把噪声频率也画进来，图上一堆杂乱的小椭圆；太高会漏掉一些幅值不大但故障特征明显的频率成分。实际使用时可以调整这个比例，比如 3% 到 10% 之间。
提示：如果某个频率下 X 或 Y 的幅值接近零，椭圆会退化成一条直线。这本身也是一种诊断信息——说明该频率下转子只在一个方向振动，可能是支撑刚度不对称导致的。
3. 全息谱诊断的判据：椭圆形状怎么读
3.1 正进动与反进动：椭圆旋转方向的含义
全息谱椭圆不只是形状有含义，旋转方向也有。转子在某一频率下的进动方向，可以通过椭圆参数方程中 X 和 Y 的相位差来判断。当 Y 方向相位落后 X 方向 90 度时，转子做正进动（与旋转方向相同）；落后 -90 度时，做反进动。
在 Matlab 里判断进动方向，可以计算两个方向相位差的正弦值：
% 判断进动方向 phase_diff = yp - xp; % 相位差 if sin(phase_diff) > 0 direction = '正进动'; elseif sin(phase_diff) < 0 direction = '反进动'; else direction = '退化为直线'; end fprintf('频率 %.1f Hz 处的进动方向: %s\n', f(idx), direction);正进动通常是正常的，反进动往往意味着存在异常，比如动静碰摩或者某些类型的油膜振荡。但这不是绝对的，有些故障也会表现为正进动但椭圆形状异常，所以进动方向要和椭圆形状结合起来看。
3.2 椭圆偏心率与主轴方向：不平衡和不对中的区分
椭圆的偏心率定义为长轴和短轴的比值。偏心率接近 1 说明是正圆，转子在该频率下各方向振动均匀；偏心率远大于 1 说明是扁椭圆，某个方向振动占主导。
计算偏心率需要先求椭圆的长轴和短轴。对于参数方程形式的椭圆，可以通过矩阵方法求解：
% 计算椭圆偏心率 % 构造椭圆参数矩阵 A = [xa*cos(xp), xa*sin(xp); ya*cos(yp), ya*sin(yp)]; % 奇异值分解求半轴长度 [~, S, ~] = svd(A); semi_major = S(1,1); % 长半轴 semi_minor = S(2,2); % 短半轴 eccentricity = semi_major / semi_minor; fprintf('椭圆偏心率: %.2f\n', eccentricity);偏心率在 1.0 到 1.2 之间，通常认为接近正圆，不平衡的可能性较大。偏心率大于 2.0，说明椭圆很扁，不对中或者支撑刚度差异大的可能性较高。主轴方向则反映了振动能量集中的方向，和传感器的安装角度、支撑结构有关。
3.3 全息谱瀑布图：变转速工况下的诊断利器
固定转速下的全息谱只能看到当前状态。如果设备在变转速工况下运行，比如启停机过程，把不同转速下的全息谱按转速排列，就得到全息谱瀑布图。这张图能清晰展示各频率成分随转速的变化趋势，对识别临界转速、共振区域特别有效。
% 假设有多个转速下的数据，存在cell数组中 % data_cell{k} 包含第k个转速下的 x, y, fs % speeds 存储对应的转速值 figure; hold on; for k = 1:length(data_cell) x = data_cell{k}.x; y = data_cell{k}.y; fs = data_cell{k}.fs; N = length(x); X_fft = fft(x - mean(x)) / N; Y_fft = fft(y - mean(y)) / N; half_N = floor(N/2) + 1; X_spec = X_fft(1:half_N); Y_spec = Y_fft(1:half_N); f = (0:half_N-1) * fs / N; % 找转频对应的峰值 [~, idx] = min(abs(f - speeds(k)/60)); xa = 2*abs(X_spec(idx)); xp = angle(X_spec(idx)); ya = 2*abs(Y_spec(idx)); yp = angle(Y_spec(idx)); theta = linspace(0, 2*pi, 100); xe = xa * cos(theta + xp); ye = ya * cos(theta + yp); % 在三维空间中按转速排列 plot3(xe, ye, speeds(k)*ones(size(xe)), 'b-', 'LineWidth', 1); end xlabel('X方向'); ylabel('Y方向'); zlabel('转速 (rpm)'); title('全息谱瀑布图'); grid on; view(45, 30);这段代码把每个转速下的转频椭圆画在对应的 Z 高度上，形成三维瀑布图。观察椭圆随转速的演变，如果某个转速附近椭圆突然变大或者形状突变，那个转速就值得重点关注。
4. 避坑与排查：全息谱诊断中最容易翻车的五个地方
4.1 两通道相位不同步：椭圆形状完全失真
现象：画出来的全息谱椭圆形状怪异，和预期完全不符，换个时间段采的数据椭圆形状又变了。
原因：两个通道不是同步采样，存在固定的时间延迟。时间延迟会导致相位差偏移，椭圆形状随之改变。这种偏移在普通频谱上看不出来，因为幅值不受影响，但全息谱对相位极其敏感。
解决：检查采集设备是否支持同步采样。如果设备只有分别采样模式，需要在时域做互相关对齐。用xcorr计算两通道的延迟量，然后对其中一路做时移补偿：
% 互相关对齐两通道 [cc, lag] = xcorr(x, y); [~, max_idx] = max(abs(cc)); delay = lag(max_idx); % 对y通道做时移补偿 if delay > 0 y = [y(delay+1:end); zeros(delay, 1)]; elseif delay < 0 y = [zeros(-delay, 1); y(1:end+delay)]; end对齐之后重新做 FFT 和椭圆合成。如果延迟量超过一个采样周期，说明采集系统本身有问题，建议换设备。
4.2 采样率不够：高频成分混叠成假椭圆
现象：全息谱图上出现一些莫名其妙的频率成分，椭圆形状无法解释。
原因：采样率低于信号最高频率的两倍，发生混叠。高频成分被折叠到低频区域，在全息谱上表现为虚假的频率成分。
解决：采集前确认采样率满足奈奎斯特条件。对于旋转机械振动信号，一般采样率设为最高关注频率的 2.5 到 5 倍。如果已经采了数据，检查频谱上是否有异常的高频成分突然截断，那通常是混叠的迹象。补救办法是加抗混叠滤波器重采，已经混叠的数据无法通过后期处理恢复。
4.3 只看转频椭圆：漏掉真正的故障特征
现象：转频处的椭圆看起来正常，但设备实际有故障。
原因：很多故障的特征频率不是转频，而是转频的分数倍或倍数。比如油膜涡动特征频率约为转频的 0.42 到 0.48 倍，不对中的特征频率是 2 倍转频。只看转频会漏掉这些关键信息。
解决：全息谱分析要覆盖整个关注频段，至少包括 0.1 倍到 10 倍转频的范围。用前面 2.4 节的批量绘图方法，把所有幅值显著的频率点都画出来。重点关注 0.5X、1X、2X、3X 这些位置。
4.4 相位基准漂移：不同时间采的数据无法对比
现象：同一台设备，今天采的数据和昨天采的数据，全息谱椭圆形状差异很大，但设备状态并没有明显变化。
原因：相位是相对于触发时刻的。如果每次采集的触发时刻不固定，绝对相位值就会漂移。全息谱椭圆虽然主要看相位差，但相位差也会受到触发时刻不一致的影响。
解决：使用键相传感器提供固定的相位基准。每次采集都以键相信号为触发，保证相位参考一致。如果没有键相传感器，至少保证每次采集的触发条件相同，比如都用同一个振动通道的过零触发。
4.5 椭圆拟合时频率选错：相邻谱线干扰
现象：椭圆形状不稳定，稍微改变分析频率，椭圆就完全变样。
原因：FFT 的频率分辨率有限，如果两个频率成分靠得很近，选错谱线索引会取到相邻频率的幅值和相位，导致椭圆形状跳变。
解决：提高频率分辨率，增加采样点数或者降低采样率。同时用峰值搜索而不是固定频率索引来定位目标频率：
% 在目标频率附近搜索局部峰值 search_range = 5; % 搜索范围，单位Hz freq_mask = abs(f - target_freq) < search_range; local_idx = find(freq_mask); [~, peak_offset] = max(X_amp(local_idx)); idx = local_idx(peak_offset);这样即使目标频率和 FFT 谱线不完全对齐，也能找到最近的峰值位置，避免取到相邻谱线的数据。
5. 进阶技巧：用全息谱做趋势监测与自动报警
前面讲的都是单次分析。实际工程中更有价值的是连续监测——每天或每周采一次数据，看全息谱椭圆参数的变化趋势，在故障早期就发出预警。
我一般会提取几个关键特征量做趋势跟踪：转频椭圆的偏心率、长轴长度、进动方向。这三个量对不平衡、不对中、碰摩等常见故障比较敏感。下面是一个简单的趋势监测脚本框架：
% 趋势监测：假设有连续多天的数据 % data_log 是一个结构体数组，每天一条记录 % 每条记录包含 x, y, fs, timestamp num_days = length(data_log); ecc_trend = zeros(num_days, 1); major_trend = zeros(num_days, 1); direction_trend = zeros(num_days, 1); for d = 1:num_days x = data_log(d).x; y = data_log(d).y; fs = data_log(d).fs; N = length(x); X_fft = fft(x - mean(x)) / N; Y_fft = fft(y - mean(y)) / N; half_N = floor(N/2) + 1; X_spec = X_fft(1:half_N); Y_spec = Y_fft(1:half_N); f = (0:half_N-1) * fs / N; % 假设转频已知为50Hz [~, idx] = min(abs(f - 50)); xa = 2*abs(X_spec(idx)); xp = angle(X_spec(idx)); ya = 2*abs(Y_spec(idx)); yp = angle(Y_spec(idx)); % 计算椭圆参数 A = [xa*cos(xp), xa*sin(xp); ya*cos(yp), ya*sin(yp)]; [~, S, ~] = svd(A); ecc_trend(d) = S(1,1) / max(S(2,2), 1e-6); major_trend(d) = S(1,1); direction_trend(d) = sign(sin(yp - xp)); end % 绘制趋势图 figure; subplot(3,1,1); plot(ecc_trend, 'b-o'); ylabel('偏心率'); title('转频椭圆偏心率趋势'); grid on; subplot(3,1,2); plot(major_trend, 'r-o'); ylabel('长半轴'); title('转频椭圆长半轴趋势'); grid on; subplot(3,1,3); plot(direction_trend, 'k-o'); ylabel('进动方向'); title('进动方向趋势 (1=正进动, -1=反进动)'); grid on;这段代码每天提取三个特征量，画成趋势曲线。报警阈值可以这样设：偏心率连续三天超过基线值的 20%，或者长半轴连续一周单调递增，或者进动方向发生翻转。这些阈值需要根据具体设备的历史数据来调整，没有万能值。
注意：趋势监测的前提是每次采集的工况一致。如果转速、负载、温度变化很大，趋势曲线会有很大波动，需要先做工况归一化，或者只在相同工况下对比。
还有一个实用技巧：把全息谱椭圆和普通频谱放在一起看。频谱告诉你哪些频率有能量，全息谱告诉你这些频率下的振动形态。两者结合，诊断准确率会明显提高。我自己的习惯是，先扫一遍频谱找异常频率，再对这些频率画全息谱椭圆，最后对照趋势曲线判断故障发展阶段。
这套方法我用了几年，最深的体会是：全息谱的价值不在于单张图多好看，而在于它把相位这个被忽视的维度利用起来了。很多在频谱上看不出来的问题，全息谱上一眼就能看出不对劲。但前提是数据质量要过关，同步采样、相位基准、采样率这三件事没做好，后面画出来的椭圆都是玄学。希望帮到你。
本文还有配套的精品资源，点击获取
