“循环平稳信号”——这个词对很多刚接触信号处理的人来说有点绕口。但如果你的研究对象是旋转机械、轴承、齿轮箱那这个词迟早会出现在你面前。它们的故障振动信号天然就带有二阶循环平稳特性周期性的冲击激励起高频共振故障特征频率每隔固定间隔就把这个共振“调制”一次。放在频谱上看这些信息往往被噪声和低频大能量淹没根本看不出名堂但放在循环频率维度上它们会变成非常清晰、稳定的窄带峰。这也是为什么快速谱相关分析算法能在MATLAB环境下成为诊断早期微弱故障的利器。要分析这类信号快速谱相关是目前很实用的一条技术路线。它把常规频谱里的“一条线”扩展成“循环频率—谱频率”的二维平面故障特征和共振频带自动分离不需要像包络谱那样人工挑频带。这篇东西我会从循环平稳的原理讲起把快速谱相关算法的核心思路拆开给出完整的MATLAB实现和轴承故障仿真用例最后分享我调试这套方法时踩过的一些坑。想把循环平稳分析真正用在实测数据上的这篇应该能给你省不少时间。1. 先搞懂要分析的对象循环平稳信号1.1 平稳、非平稳、循环平稳区别在哪做信号处理的都知道平稳信号是“统计特性不随时间变化”的信号白噪声和稳定的正弦波就是典型。非平稳信号则相反统计特性随时在变比如语音、瞬态冲击。循环平稳刚好夹在中间它的统计特性不恒定但也不乱变而是有规律地周期变化。用数学化一点的语言说如果一个随机信号的自相关函数满足R_x(t, τ) E{x(tτ/2) x*(t−τ/2)}并且对任意时间点t都有 R_x(tT0, τ) R_x(t, τ)其中T0是某个固定周期那它就是二阶循环平稳信号。这个概念难在它不像频谱那样直观但它恰好能描述一大类工程信号的本质规律。1.2 轴承故障信号为什么会循环平稳以滚动轴承为例。外圈出现局部缺陷时滚动体每次经过缺陷位置都会撞一下产生一个瞬态冲击。这个冲击会激励起轴承座和传感器安装结构的固有频率形成一段高频衰减振荡。问题的关键在“每次”上冲击不是随机的而是以故障特征频率BPFO为周期反复出现。于是信号的均值虽然没有明显变化噪声背景大但信号的能量——或者说二阶矩——会随BPFO的周期而波动。这正是循环平稳的定义统计特性周期性变化。转频、齿轮啮合、不对中、松动等故障都会在振动信号中留下类似的循环平稳特征。换句话说循环平稳不是某个特定故障独有的性质而是旋转机械振动信号的一种普遍属性。1.3 传统频谱和包络谱为什么不够用我最早做轴承诊断时用的也是频谱和包络谱。频谱的好处是简单但问题在于冲击调制信息藏在高频共振带里频谱上只能看到一个宽宽的“鼓包”根本看不出120Hz这种调制频率。包络谱倒是能把调制频率翻出来可它有个前提——你得先人工指定共振频带。这个频带一旦选偏后面的结果就全完。谱相关分析的好处在于不依赖人工选带。它把信号在“谱频率f”和“循环频率α”两个维度上同时展开共振频带在f轴上自然呈现调制频率在α轴上自然呈现二维图上哪个位置有亮峰一目了然。相比包络谱这是一种更“全局”的分析方式。下面这张表可以把差异看得更清楚方法需要提前选频带抗噪能力输出维度对故障调制的揭示能力时域波形否弱时间轴几乎无功率谱否中频率轴差包络谱是中频率轴中谱相关/相干否强f-α二维平面强仅从这张表就能看出谱相关是更贴合循环平稳信号本质的工具。问题只在于传统的谱相关计算量很大直接用于工程信号不太现实所以才需要“快速”谱相关算法来做加速。2. 谱相关算法的原理以及“快”在哪里2.1 谱相关密度与循环谱相干谱相关密度的定义在不同教材里写法不太一样但意思类似。设信号x(t)的傅里叶变换为X(f)定义S_x(f, α) E{ X(f α/2) X*(f − α/2) }这里面f是谱频率α是循环频率。这个公式的含义可以这么理解在频谱上相距α的两个频率分量它们的相关性有多大。普通功率谱只考虑α0这一条线也就是同一频率分量自身的能量谱相关则把所有可能的频率间隔都翻了一遍。如果信号里存在周期为1/α的调制那么间隔为α的谱分量之间会有稳定的相位关系表现在S_x(f, α)上就形成一个峰。直接拿S_x的绝对幅值来做分析有个问题不同频带的能量本来就有差异强能量频带会把弱调制淹没。工程上更常用的指标是循环谱相干系数γ相当于把S_x归一化到0~1之间γ_x(f, α) S_x(f, α) / √[ S_x(fα/2, 0) · S_x(f−α/2, 0) ]分母实际上是两条“零循环频率线”上的功率谱乘积。这样一来不管频带本身能量高低只要调制关系够强γ都会接近1。我做故障诊断时基本都用γ因为它比谱密度幅值稳定得多也不受传感器灵敏度标定影响。2.2 传统谱相关计算为什么不实用理论上最直接的谱相关估计方法是把信号切成多段对每一段先做带通滤波再做互谱平均。这种办法有两个绕不开的瓶颈一是循环频率α的扫描数量很大二是每换一个α都要重新计算一组互谱总体计算复杂度接近O(N²)甚至更高。举个例子信号长度2秒采样率25.6kHz总共51200个采样点。如果要把0~500Hz的循环频率都扫一遍频率分辨率做到1Hz那就要计算500条谱相关线每条线上又是25600个频率点。直接算的话内存和时间都很难接受。这也是循环平稳分析在早期只停留在理论圈子的原因直到快速算法出现之后它才真正变成工程上可用的工具。2.3 快速谱相关的两条实现路线现在工程界常用的快速谱相关算法大致分两类。第一类基于短时傅里叶变换STFT思路是先对信号做一次STFT得到时频矩阵X(m, k)其中m是帧序号k是频率bin。然后对每个频率bin把相隔d个bin的两个STFT系数序列做互相关平均就能同时得到频带和循环频率信息。因为STFT本身可以用FFT加速整体复杂度是O(N log N)比传统方法低了一个量级。第二类基于滤波器组分解把信号分成多个窄带子信号再对每个窄带的包络做FFT得到循环频率特征。这两种路线各有优劣STFT实现简单、容易理解但频率分辨率和循环频率分辨率会被同一个窗长耦合在一起滤波器组方法实现稍复杂但分辨率可以独立控制结果也更干净。实际项目中我通常先用STFT方法快速摸底确认信号里确实有循环平稳成分之后再上滤波器组版本做精细分析。2.4 包络谱其实是一种退化的谱相关理解谱相关和包络谱关系有个很直观的角度包络谱本质上就是对某个选定的频带f0附近的STFT系数在时间方向做解调并求FFT相当于只在ff0处沿着α方向切了一刀。谱相关则是一刀都没切直接把整个f-α平面算出来。所以谱相关里天然包含了包络谱的信息还额外给出了“调制发生在哪个频带”这个信息。反过来看如果三维图中峰的α坐标指向某个故障特征频率f坐标指向某个共振频带那这个结果就同时回答了两个问题信号里调制周期是多少以及调制能量集中在哪里。这比单纯输出一条包络谱曲线信息量大得多。3. MATLAB环境下从零实现快速谱相关3.1 仿真轴承故障信号怎么构造先把测试信号造出来后续所有算法验证都在这个信号上进行。仿真对象为滚动轴承外圈单点故障参数尽量贴近实测fs 25600; % 采样率 25.6 kHz T 2; % 时长 2 秒 N round(fs * T); t (0:N-1) / fs; fr 30; % 轴转频 30 Hz bpfo 120; % 外圈故障特征频率 120 Hz f_res 4500; % 结构共振中心频率 4.5 kHz decay 800; % 冲击衰减系数 x zeros(1, N); impulse_len round(0.01 * fs); % 单次冲击约 10 ms impulse_t (0:impulse_len-1) / fs; impulse exp(-decay * impulse_t) .* sin(2*pi*f_res*impulse_t pi/6); % 按 BPFO 周期叠加冲击 for k 0 : ceil(T * bpfo) - 1 pos round(k / bpfo * fs) 1; if pos impulse_len - 1 N x(pos:posimpulse_len-1) x(pos:posimpulse_len-1) impulse; end end % 加入转频幅度调制模拟载荷变化 x x .* (1 0.25 * sin(2*pi*fr*t)); % 加入高斯白噪声 x x 0.6 * randn(size(t)); save(sim_bearing_fault.mat, x, fs, fr, bpfo);这段代码里有两个地方值得说明。第一冲击衰减系数800对应的时间常数为1/800秒约1.25ms5倍时间常数大概6.25ms10ms的冲击长度足够让振荡基本衰减完既能模拟真实冲击又不会让相邻冲击严重重叠。第二转频幅度调制是刻意加进去的它会在谱相关图的低循环频率区域产生30Hz、60Hz、90Hz等成分用来观察算法是否能区分故障调制和转频调制。3.2 基于STFT的快速谱相关核心函数下面给出完整的核心函数fsc_stft。它接受信号、采样率、窗长、重叠率和最大循环频率输出循环谱相干矩阵。function [gamma, f, alpha] fsc_stft(x, fs, nfft, overlap_ratio, alpha_max) % 基于STFT的快速谱相关分析 % 输入 % x - 信号列向量/行向量 % fs - 采样率单位 Hz % nfft - STFT窗长度同时也是FFT点数 % overlap_ratio - 帧重叠率建议0.7~0.85 % alpha_max - 需要分析的循环频率上限单位 Hz % 输出 % gamma - 循环谱相干系数矩阵尺寸 (nfft/21) x (dmax1) % f - 谱频率轴单位 Hz % alpha - 循环频率轴单位 Hz if isrow(x) x x(:); end win hann(nfft, periodic); noverlap round(nfft * overlap_ratio); [X, f] spectrogram(x, win, noverlap, nfft, fs); % 窗能量归一化避免窗函数影响功率幅度估计 win_norm sqrt(sum(win.^2)); X X / win_norm; [Nf, ~] size(X); P mean(abs(X).^2, 2); % 各频率bin的平均功率 da fs / nfft; % 循环频率网格间隔 dmax round(alpha_max / da); alpha (0:dmax) * da; gamma zeros(Nf, dmax1); for d 0:dmax for k 1:Nf-d % 间隔 d 个频率bin的STFT系数做时间方向共轭相乘并平均 sc_val mean(X(kd, :) .* conj(X(k, :)), 2); denom sqrt(P(k) * P(kd)); if denom eps gamma(k, d1) abs(sc_val) / denom; end end end % α0 对应普通功率谱置零以避免掩盖调制成分 gamma(:, 1) 0; end这个函数的核心逻辑非常朴素对每个频率偏移d把第k个和第kd个频率bin的STFT系数在时间方向做共轭乘积再平均除以两端功率谱的几何平均得到相干系数。d从0到dmax变化对应循环频率从0到alpha_max每一列就是一条沿频谱方向的“循环切片”。你可能会觉得这不算“快速”——毕竟内外两层循环还在。但要注意计算STFT的主循环只跑了一次FFT带来的复杂度是O(N log N)而原来传统算法需要对每个α单独做一次长傅里叶分析。实际运行起来这段代码在2秒、25.6kHz的信号上基本一两秒就能出结果这就是“快速”的含义。3.3 运行主流程并绘制谱相关图在MATLAB脚本中调用并把结果画出来clear; close all; clc; load(sim_bearing_fault.mat, x, fs, fr, bpfo); nfft 2048; [gamma, f, alpha] fsc_stft(x, fs, nfft, 0.75, 500); % 二维谱相关图 figure(Color, w, Position, [100 100 860 520]); imagesc(alpha, f/1000, gamma); axis xy; xlabel(循环频率 \alpha (Hz)); ylabel(谱频率 f (kHz)); title(快速谱相关分析循环谱相干); h colorbar; h.Label.String \gamma(f,\alpha); colormap(jet); ylim([2 8]); % 只显示共振频带附近区域 set(gca, FontSize, 12); % 循环频率等于BPFO处的切片 [~, idx] min(abs(alpha - bpfo)); fprintf(实际选取循环频率: %.1f Hz\n, alpha(idx)); figure(Color, w, Position, [100 100 700 420]); plot(f/1000, gamma(:, idx), LineWidth, 1.2); xlabel(谱频率 f (kHz)); ylabel(\gamma(f, \alpha BPFO)); title(BPFO对应的谱相关切片); grid on; xlim([0 13]);这里nfft取2048时循环频率网格间隔是fs/nfft12.5HzBPFO120Hz落在网格附近但不在网格正中间。实际运行时会发现选到的循环频率可能是125Hz而不是120Hz谱峰幅度会有轻微衰减。这个问题在算法层面能解决我放到后面“问题排查”那节细说。3.4 关键参数怎么选快速谱相关最关键的参数就是nfft、重叠率和alpha_max它们的判断逻辑和普通STFT不太一样我整理成了表格参数推荐范围说明nfft512 ~ 8192同时决定谱频率分辨率Δf和循环频率网格ΔαΔf需要权衡overlap_ratio0.7 ~ 0.85增加时间平均帧数降低谱估计方差代价是计算时间略增alpha_max0.5~2 kHz取决于故障特征频率的范围轴承外圈故障一般500Hz内够用信号长度尽量2秒以上总时长决定循环频率的理论分辨率1/Tnfft的选择要仔细想nfft取得越大频率分辨率越好但循环频率网格间隔也随之变大可能漏掉细微的调制特征nfft取得太小谱频率分辨率不够共振频带和调制成分会糊在一起。工程上我一般先用nfft1024或2048快速跑一遍看大致结构再根据想关注的故障频率尝试微调。alpha_max的选择比较简单只要略大于你关心的故障特征频率就行没必要把所有循环频率都算出来。算得越多矩阵越大计算时间成倍增加。我没有把alpha_max设到fs/2这种极端情况实际也没必要。4. 完整案例轴承外圈早期故障诊断4.1 案例流程与初筛用上一节的仿真信号做完整分析。流程分四步生成/读取信号先做时域和频谱观察。运行fsc_stft得到γ(f, α)。观察二维图定位亮斑坐标。取特征循环频率切片与包络谱对比。时域波形里能看到周期性冲击但冲击频率和噪声混在一起肉眼数出来的周期不完全等于BPFO。再看功率谱它显示的频谱在4~6kHz附近有一个明显的隆起这就是共振频带但功率谱里完全看不出120Hz的调制周期。4.2 二维谱相关的解读方法运行fsc_stft后二维图上会出现几个值得注意的元素。最明显的是位于f约4.5kHz、α约120Hz处的高亮条带这就是外圈故障的“调制—共振”特征能量集中在结构共振频带调制周期等于BPFO。这从颜色深浅就能直接看出来。与此同时在α30Hz、60Hz、90Hz等位置也有较弱的成分主要集中在f低于2kHz的区域。这些来自转频的幅度调制是载荷变化带来的正常现象。需要注意故障特征120Hz和转频调制30Hz在图上处在不同的频率区域这个信息在包络谱里很难分得这么清楚。画完γ(f, α)图后在α120Hz附近做一个切片就能看到γ沿频率f的分布。这个切片图上的峰所在位置正好对应共振频带中心4.5kHz峰幅值接近0.8。这等于谱相关方法自动帮你找到了“该解调哪个频带”——不需要人工选择。4.3 和人工选带包络谱放在一起比为了说明谱相关的优势我把传统的包络谱也跑一遍% 选对共振频带 [b1, a1] butter(4, [3600 5400]/fs*2, bandpass); x_env1 filtfilt(b1, a1, x); env1 abs(hilbert(x_env1)) - mean(abs(hilbert(x_env1))); env_spec1 abs(fft(env1 .* hann(N))); f_axis_env (0:N-1)/N*fs; % 选错频带 [b2, a2] butter(4, [200 1200]/fs*2, bandpass); x_env2 filtfilt(b2, a2, x); env2 abs(hilbert(x_env2)) - mean(abs(hilbert(x_env2))); env_spec2 abs(fft(env2 .* hann(N))); figure(Color, w, Position, [100 100 820 620]); subplot(2,1,1); plot(f_axis_env(1:N/2), env_spec1(1:N/2)); xlim([0 500]); grid on; title(包络谱共振频带选对 (3.6~5.4kHz)); xlabel(频率 (Hz)); ylabel(幅值); subplot(2,1,2); plot(f_axis_env(1:N/2), env_spec2(1:N/2)); xlim([0 500]); grid on; title(包络谱共振频带选错 (0.2~1.2kHz)); xlabel(频率 (Hz)); ylabel(幅值);选对频带时包络谱在120Hz处有明显的峰选错频带时包络谱上120Hz处几乎什么东西都没有。这说明包络谱完全依赖人工经验频带选错故障特征直接就丢了。而谱相关不需要这一步即使你提前不知道共振频带在哪只要跑一遍γ(f, α)图亮斑位置会同时告诉你“共振频带在哪”和“调制频率是多少”。这就是它的核心价值。4.4 对结果的工程判断拿到谱相关结果后工程上有两件事可以做。第一看γ(f, α)图的峰值是否稳定——如果多次测量、多段数据都能在同一(f, α)位置出现峰说明特征是可复现的可以作为故障判据如果只是随机噪声形成的孤立亮点多次测量一般位置会漂移。第二峰值对应的α频率是否和轴承参数计算的BPFO接近。轴承出厂参数能算出理论BPFO实测和理论偏差一般在百分之几以内偏差过大时反而要考虑是不是其他故障。我在实际项目中还习惯把γ(f, α)图转成灰度显示再用伪彩色调色板突出高强度区域。做故障预警时只提取特定循环频率处的γ值比如γ_120把它作为趋势量跟踪。如果γ值随运行时间缓慢爬升说明缺陷在扩展如果只是单点跳变可能是瞬态干扰要多看几段数据再下结论。5. 实际使用中踩过的坑和解决办法5.1 特征频率落在网格缝隙里这是我遇到最多的问题。STFT实现的循环频率网格间隔是fs/nfftnfft是2的整数幂时可选的间隔就是6.25Hz、12.5Hz、25Hz等。如果BPFO正好是120Hz用nfft2048时最近的网格是125Hz虽然差值不大但谱峰能量会被分散峰值读数和真实值有偏差。解决思路有几个第一换非2的整数幂窗长让da约等于目标频率的约数比如nfft2560时da10Hz120Hz恰好落在第12条循环频率线上第二对网格搜索结果做抛物线插值在峰值附近利用相邻点估计更精确的峰位第三改用滤波器组类型实现让循环频率分辨率由总时长决定与窗长无关。我一般先试第一种改几个nfft看结果是否收敛再用第三种的思路做最终版本。5.2 噪声太大二维图上看不到明显亮点谱相关看起来是二次统计量对噪声有一定抑制能力但信噪比特别低时γ图也会变得很花。这时候不要急着调小窗口先检查几个常规因素重叠率是不是太低、信号时长是不是太短、alpha_max是不是设得过大导致杂散循环频率太多。实际调试中我建议把γ阈值提高一些只看大于0.3的像素。另外可以做一个简单的二维中值滤波窗口取3×3或5×5去掉孤立亮点再找局部最大值。这个方法不严谨但很实用做故障初筛时能省很多时间。5.3 出现大量难以解释的循环频率成分实测数据里除了轴承故障还有电机转频、齿轮啮合、皮带轮转速等周期源。这些周期源都可能产生循环平稳成分导致γ图上出现很多峰。最典型的例子是50Hz电气噪声及其倍频它们会在α50Hz、100Hz、150Hz处形成很规则的竖线很容易被误读为故障特征。遇到这种情况一定要把已知的频率成分列出来在γ图上先标出“预期位置”再去看剩余的能量团。比如已知转频30Hz那就在α30Hz的倍数位置标记确认这些位置的峰值来自正常调制再关注α120Hz处的异常峰。做故障诊断永远先做排除法而不是看到峰就报警。5.4 计算时间太长与内存溢出当nfft选到8192、alpha_max选到2000Hz、信号长度又很长时gamma矩阵的尺寸和循环次数都会明显增加。这时候可以先降采样把采样率降到故障特征频率的5~10倍既不影响分析又能大幅减少数据量。其次把d的最大值限制在必要范围内比如轴承外圈故障特征频率不超过300Hz就设alpha_max500。第三用parfor加速外层的d循环如果并行池空闲速度能提升好几倍。内存方面gamma矩阵的元素都是double类型尺寸为(Nf)×(dmax1)。以nfft4096、alpha_max1000Hz为例da大约是6.25Hzdmax160Nf2049矩阵总元素约33万占用约2.6MB这个量级完全可控。真正占资源的是长信号的STFT结果X矩阵切片数量太大时建议改用信号分段处理每段算完就释放再合并结果。5.5 谱相关与包络谱结果不一致有时候谱相关图显示某个循环频率有能量但包络谱却没有对应峰值反过来也一样。这通常是因为包络谱只用了单一频带而谱相关在计算时会综合所有频带的贡献。选对了频带两者自然一致选错了频带只有谱相关能给出正确结果。反过来如果谱相关图出现峰而包络谱没有去检查峰值对应的f坐标把它作为窄带滤波的中心频率再做一次包络谱就能验证结果。这个交叉验证思路我在写故障诊断报告时经常用客户也更信服。最后再分享一个个人习惯。我在每次做完谱相关分析后都会顺手把γ(f, α)矩阵保存成MAT格式方便后面换不同阈值重新画。后面做数据分析时经常要反复调整显示范围、重点看某一段频率如果每次都重新算一遍很浪费时间直接读矩阵再切片效率完全不一样。还有第一次跑通算法后一定先用仿真信号验一遍确保参数和输出轴都对再上实测数据不然错误很可能藏在你看不见的地方排查起来非常痛苦。这套快速谱相关流程我建议你也在自己的项目里跑一遍试试。