MATLAB快速谱相干：从一维时间序列到旋转机械多通道分析

前几天我在一个设备诊断交流群里看到有人贴图同一条轴上的两路振动信号普通幅值谱看着都差不多在某个轴承故障特征频率附近却同时出现了一处明显的相干峰。下面跟了几条回复有人问“相干峰到底代表什么”有人说“这方法我试过但求出来全是一堆 1”还有人问“能不能直接套到多通道阵列上”。这个话题正好撞上标题里那件事——用 MATLAB 在一维时间序列上做快速谱相干并且把它从旋转机械诊断用到多元信号分析。我最早接触谱相干不是为了诊断而是处理两个麦克风采集的信号想找出它们在哪个频段真正“共用同一个振源”。后来做旋转机械测试时发现这东西比单纯看频谱更擅长识别周期性的故障冲击尤其在背景噪声很大的现场。写这篇东西是想把两条坐标轴串起来讲清楚谱相干算的是什么、MATLAB 里怎么写才能快、旋转机械场景怎么做、最后怎么扩展到多通道。适合刚接触振动分析信号处理的读者也适合已经会调用现成函数但说不清参数意义的人。1. 谱相干到底在算什么一个频域里的相关性1.1 从普通相干到幅值平方相干先回忆一下最基础的概念两个一维信号 x(t)、y(t) 的互相关度量它们在时间轴上“对齐”的程度而谱相干把这种相关性拆到每个频率分量上。工程里最常用的是幅值平方相干定义式是MSC(f) |Sxy(f)|² / (Sxx(f) · Syy(f))这里 Sxx(f)、Syy(f) 是两个信号的自功率谱Sxy(f) 是互功率谱取值的范围在 0 到 1 之间。0 表示这个频率点上两个信号完全不相关1 表示一个信号可以在该频率完全由另一个线性表示。可以把它理解成“频段里的 R²”回归里的决定系数到了频谱域。在实际计算里自谱和互谱都不是一次 FFT 就能得到的稳定结果而是需要做多次平均。这也是为什么谱相干天然带“统计”属性你把它看成一条曲线的时候背后其实是一组样本在支撑。样本越少曲线越飘越容易出现某个频点莫名其妙接近 1 的情况。1.2 旋转机械为什么要用谱相干旋转机械的振动信号通常由三部分构成确定性周期分量、循环平稳分量、随机噪声。确定性分量的频率位置固定比如轴转频率和它的谐波、齿轮啮合频率。循环平稳分量的典型例子是轴承局部损伤产生的冲击每次经过损伤点时产生一次冲击冲击的幅度还受轴承载荷调制。随机噪声则来自现场环境、电噪声和传感器自身的底噪。普通幅值谱把所有能量都画在同一张图上随机噪声虽然不集中但会抬高底噪把一些小特征淹没。谱相干的做法不一样它通过两个通道的“相位一致性”来判定这个频率上是否存在共同的周期性激励。如果两路信号里的随机噪声互不相关这部分在互谱里平均后会趋近于 0导致相干值接近 0而两路信号共同拾取的故障冲击其相位与幅值在多次平均中保持稳定相干值就高。换句话说谱相干对“非相干背景”天然免疫。带上它到现场最大的好处是能够把“大家都有的振动”和“真正共同相关的振动”区分开。这对大机组很关键因为相邻设备都会通过基础和管道传递振动只看幅值根本分不清来源。1.3 一维时间序列里“快”的含义标题里出现“快速”两个字重点不是算法理论多高深而是 MATLAB 实现方式的问题。我第一次写谱相干时老老实实用一个循环逐段处理四分钟的数据要跑将近半分钟。后来改成矩阵化分帧把整段数据切成若干行一次性做 FFT时间降到了两三秒。这里的关键不是某个函数而是“分段-加窗-FFT-平均”这个过程能否并行化。对一维时间序列多数现成函数比如内置的单通道相干估计函数也足够应付两个信号的场合。可一旦要扩展到八个通道、十六个通道或者要在每个频点上都做统计检验懂手写实现就非常重要你能控制中间量的存取方式能批量构建交叉谱矩阵还能在不算互谱的地方直接跳过计算。所谓“快”表面是代码效率本质是对计算流程的控制。2. MATLAB 快速实现抛开循环用矩阵化帧计算2.1 设计思路Welch 估计与手写函数谱相干的核心是谱估计而 Welch 方法就是“信号分帧、加窗、FFT、平均”。对比一下最直观的写法for k 1:nseg seg data(idx(k):idx(k)nfft-1); Xf fft(seg .* win); Sxx Sxx abs(Xf).^2; end这个写法没有任何理论问题问题出在约束环境下的效率。MATLAB 本身对向量化计算友好把“纵向时间”的循环改成“横向批次”的矩阵往往能获得十倍以上的加速。实现思路是先建立一个索引矩阵每段数据占一行然后把整块矩阵和窗函数逐行相乘最后沿着第一维做平均。我就是用这种方式写了自定义的谱相干函数只依赖基础功能不等同于工具箱内置的现成接口。这个函数可以接受两个一维序列返回频率轴和相干曲线。它的思想和厂商提供的现成函数完全一致但结构更透明方便在旋转机械场景里继续扩展。2.2 从分帧到结果的完整实现下面是我在实际项目里用过的一个版本。它分成三步先做信号对齐和去均值再通过索引矩阵把所有段一次性取出来最后用矩阵化 FFT 计算自谱和互谱。function [f, Cxy] fast_mscohere(x, y, Fs, nfft, win, ov) if nargin 6 || isempty(ov) ov 0.5; end if nargin 5 || isempty(win) win hann(nfft, periodic); end x x(:); y y(:); n min(length(x), length(y)); x x(1:n) - mean(x(1:n)); y y(1:n) - mean(y(1:n)); nadv round(nfft * (1 - ov)); noverlap nfft - nadv; nseg floor((n - noverlap) / nadv); if nseg 4 error(分段数太少请增大信号长度或减小nfft); end win win(:); idx (0:nseg-1) * nadv (1:nfft); XF fft(x(idx) .* repmat(win, nseg, 1), nfft, 2); YF fft(y(idx) .* repmat(win, nseg, 1), nfft, 2); Sxx mean(abs(XF).^2, 1); Syy mean(abs(YF).^2, 1); Sxy mean(conj(XF) .* YF, 1); Cxy abs(Sxy).^2 ./ (Sxx .* Syy eps); Cxy Cxy(1:nfft/21); f (0:nfft/2) * Fs / nfft; end有几个细节专门说一下。第一加窗后自谱和互谱都没有做能量归一化这是故意留的因为相干值的分子分母同时包含同一种平方项归一化因子会约掉不影响最终结果。第二互谱用的是 conj(XF) .* YF 而不是直接点乘这样得到的复数相位符合“以 x 为参考、看向 y”的约定。第三分母里加了一个 eps 防止除零同时不会影响正常物理量级。更稳妥的做法是把自己手写的结果和工具箱内置的相干函数做一个对比比如调用内置函数对等验证。两个结果在相同窗、相同重叠率下应该基本重合。我第一次对比时发现曲线几乎压在一起误差只有 10 的负 15 次方量级。如果你也想检验手写函数建议先跑一组随机数据再跑一组带有确定公共正弦分量的模拟信号。2.3 参数怎么定窗长、重叠与段的博弈这里的问题不是“哪个函数对”而是“参数为什么这么设”。我把常用参数列了一张表直接照抄也能用但更重要的是理解原因。参数常用取值范围影响窗型Hanning/Hann 周期性窗减小频谱泄漏区分故障特征频率附近的边带nfft512 到 16384 之间决定频率分辨率越长分辨率越细但单段包含的突发冲击变多重叠率50% 到 75%提高段落间的统计代表性因为加窗后边缘信息损失最少段数8 到 16 段以上段落太少相干值的统计波动会非常大去均值必须做直流分量会让零频附近相干指标失真频率分辨率是 Fs / nfft。以 12800 Hz 采样率为例nfft 取 8192分辨率约 1.5625 Hz。对于轴速 29.4 Hz、轴承故障特征频率约 95 Hz 的场景这个分辨率足够分辨故障频率和转频的边带。如果 nfft 取 256分辨率变成 50 Hz故障特征频率和边带会糊成一片相干峰也不明显。重叠率也不能只追求高。重叠过高会让相邻段落之间的样本相关性变强看上去段数很多实际有效独立样本并没有那么多。我自己的经验是普通诊断场景用 50% 到 75%如果信号本身比较平稳、只需要一条平滑的相干曲线可以选 75%如果希望保留冲击细节50% 更稳。3. 旋转机械场景从故障特征到事件同步处理3.1 一个典型的轴承故障分析流程把谱相干接到旋转机械诊断最常见的操作流程不是上来就算两路原始振动信号的相干而是先做带通滤波再做包络最后对包络信号或经过相位解调的对数包络结果做相干分析。原因很简单轴承早期故障的冲击能量通常分布在高频共振区原始信号里的低频轴频分量会占据相干度主导地位掩盖故障边带。完整步骤可以是这样的对两路加速度计信号做去趋势和去均值处理。根据轴承共振频带选择带通范围比如 3000 Hz 到 9000 Hz。用希尔伯特变换提取包络exp abs(hilbert(bandpass_filtered))。对两路包络信号计算谱相干重点观察轴转频率以及轴承故障特征频率处是否有尖峰。对比包络谱和相干谱如果包络谱有峰而相干谱没有说明这个峰可能来自单通道噪声或非共同激励源。我做过一个模拟测试假设转频 29.4 Hz外圈故障特征频率约 95.5 Hz设置一个在故障频率处周期性出现的衰减振荡再叠加随机噪声和工频干扰。单独看幅值谱时只有很微弱的故障频率峰位但计算两通道相干后在 95.5 Hz 位置出现明显的窄带尖峰而同位置的背景噪声相干值接近 0。这就是谱相干“从噪声里捞特征”的能力。3.2 变转速不要硬算转到域再相干有一个常见的坑是采集的原始数据转速并不稳定。转速一变故障特征频率在整个记录中就被“拉宽”了变成一条斜坡而不是一条线。此时直接在时间域上分段、平均、算相干结果会恶化因为同一频率处有时有特征、有时没有互谱的相位在不断变化平均后互相抵消。解决办法是把数据从“时间域”重采样到“转数域”。工程上叫阶次跟踪。如果你有键相脉冲比如每转一个脉冲就可以根据脉冲位置把每一转的信号插值成相同长度例如每转 1024 点然后把重采样后的转域信号输入谱相干函数。这时横轴不再是 Hz而是转速的倍数也就是阶次。轴承故障特征往往表现为固定的阶次转速变化不会影响其判读。没有键相脉冲时也可以用瞬时频率估计。先用短时傅里叶变化提取轴频曲线再用爬坡信号做转速参考最后插值重采样。这一套做下来谱相干的抗漂移能力会大幅提高。3.3 实操时我踩过的三个坑第一个坑是直接对原始振动信号算相干而没做任何滤波或包络。结果低频轴频分量在 1 到 30 Hz 区域相干值高达 0.98故障特征频段反而被淹没。后来我意识到谱相干分析的是“共同激励”如果轴频本身在两路信号中都足够强它当然会优先占据统计显著性。要做故障诊断需要把分析频带聚焦在故障激励所在的共振区。第二个坑是窗长选得太短。一次现场分析里我为了提升频率分辨率把 nfft 拉得很长结果窗口内包含了很多个轴的旋转周期局部冲击的位置不再固定于窗内样本点相位随机化进一步加剧相干峰反而消失了。后面改用“窗长至少覆盖两个故障周期且分段数尽量多”的原则才稳定下来。第三个坑是仿真验证时用完全相同的确定性信号做两条通道。比如 x sin(2pift)y sin(2pift)做出来的相干曲线处处等于 1。这个时候不要慌更不要怀疑算法而是要知道完全确定性的信号没有“不相关”的那部分相干当然接近 1。拿真实噪声或模拟随机故障源做验证才能看到有实际意义的谱相干形态。4. 多元信号扩展从两通道到传感阵列4.1 直接矩阵化计算交叉谱矩阵两通道谱相干的逻辑是“Sxy / sqrt(Sxx·Syy)”。扩展到多通道时自然要引入交叉谱矩阵。假设一次测量里有 ch 个传感器通道数据矩阵 X 的尺寸是 N×ch。在每个频率点 f 上我们可以构造一个 ch×ch 的厄米特矩阵 S(f)其中对角线元素是各通道的自谱非对角线元素是两两通道的互谱。这个矩阵直接用循环逐频点计算也可以但为了快最好还是走矩阵化分帧思路。下面是一个用于多通道交叉谱矩阵的函数它返回的是全部频点对应的交叉谱矩阵以及每个频点的特征值。function [f, Csm, eigVal] multichannel_csm(X, Fs, nfft, ov) [N, ch] size(X); if nargin 4 || isempty(ov) ov 0.5; end win hann(nfft, periodic); nadv round(nfft * (1 - ov)); noverlap nfft - nadv; nseg floor((N - noverlap) / nadv); nf floor(nfft / 2) 1; Csm zeros(nf, ch, ch); for k 1:nseg segIdx (k-1) * nadv (1:nfft); seg X(segIdx, :); seg seg - mean(seg, 1); F fft(seg .* repmat(win, 1, ch), nfft, 1); for m 1:nf z F(m, :).; Csm(m, :, :) squeeze(Csm(m, :, :)) (z * z) / nseg; end end f (0:nf-1) * Fs / nfft; eigVal zeros(nf, ch); for m 1:nf eigVal(m, :) sort(eig(squeeze(Csm(m, :, :))), descend); end end这个函数里没有做自谱的能量归一化因为没有必要特征值的相对大小和比值才是后续判断的核心。如果通道数较多建议把频点循环改成 parfor 并行或者只在关心的频带里计算避免在内存中堆一个巨大的三维矩阵。4.2 用特征值判断公共振源多通道交叉谱矩阵在数学上是半正定的厄米特矩阵它的特征值反映的是不同“振源子空间”的能量强度。假设传感器阵列收到一个公共周期源和几个独立的噪声源交叉谱矩阵会有一个较大的主导特征值对应公共源其余特征值接近于噪声功率的水平对应相对独立的背景噪声。这个性质非常实用。例如两组加速度计分别贴在轴承座的两侧理论上它们应该共享同一个轴承故障激励。如果算出来的交叉谱矩阵在某个频点有一个显著大的特征值说明该频点存在共同源如果特征值分布非常均匀像“平板”一样平坦则说明各通道基本独立这个频点大概率来自局部噪声。有一种工程做法是直接观察“最大特征值”与“第二大特征值”的比值。比值接近 1 说明通道间没有公共强源比值远大于 1 说明存在一个主要激励源。当然这只是一个快速判断手段严格的统计判断需要引入随机矩阵理论或显著性阈值。4.3 给多元脚本的几条工程建议第一不要把所有频点都堆成三维矩阵再去处理。如果通道数是 16频率点数是 4097一个三维矩阵就是 16×16×4097 个双精度数消耗约 8.6 MB单个也能放下但后面做特征循环时要反复挤压维度效率不高。更合理的做法是每算一个频点就立即处理特征值和主特征向量只保存少量结果。第二通道之间可能量级差异很大。某一路传感器灵敏度高另一路灵敏度低交叉谱矩阵会被高灵敏通道主导。计算前做标准化或者把每个通道除以自身标准差可以避免自谱差异掩盖真实的相关结构。第三多元场景下经常需要可视化但不要画十几条相干曲线挤在一张图里。常用做法是画“主特征值-频率”谱线再在显著频点标出对应的主特征向量。主特征向量各元素平方的分布能告诉你哪个传感器对这个振源贡献最大等于把阵列定位功能也一起做了。5. 排查与量化让谱相干结果可信5.1 相干值为什么普遍接近 1 或 0这是新手最常碰到的问题。如果整条谱线大量接近 1先按这个顺序查是否用了同一段信号的两个副本是否分段数太少是否两通道都是确定性仿真信号而没有加噪声是否忘记对信号去均值导致零频附近被直流主导每一种情况的原因和处理方式不同我整理成一个表现象典型原因处理方式几乎所有频点相干接近 1分段数太少或信号完全确定增加信号长度降低 nfft增加段数仿真时加独立噪声只在零频接近 1直流偏置去均值或高通滤波故障频段相干接近 0转速漂移、频带选择错误做转域重采样或带通滤波后再算整个曲线波动剧烈段数不足、重叠过低增加重叠率或分段数峰宽得离谱nfft 太小增大 nfft提高频率分辨率需要强调真实测量中谱相干等于 1 的频点很少。如果有往往意味着该频点只有单一确定性源且几乎没有随机误差项或者数据处理过程中犯了逻辑错误。不要轻易相信一个单频点的“完美相干值”最好看它周边频段的形态。5.2 显著阈值怎么算谱相干是一条统计曲线不能用频谱峰的高度直接判断“重要不重要”。工程上常用一个基于独立分段数的近似阈值。如果分段数为 K显著性水平为 alpha比如常见的 0.05那么经验阈值近似为threshold 1 - alpha^(1/(K-1))表面上看这个公式有点绕它其实是说在“两个通道完全不相关”的零假设下由于随机起伏某些频点的相干值也可能偶然超过某个水平。如果超过阈值就有理由认为该频点存在真实的耦合关系。举个例子分段数 K 20threshold 大约等于 0.14。如果某个频点的相干值为 0.35那它明显超过了偶然水平值得进一步分析。分段数 K 40threshold 会降到约 0.07。所以把 nseg 加到 40 以上不仅让曲线平滑也让显著性判断更严格。我在报告里一般会把阈值画成一条横线只有超过横线的峰才做进一步故障解释否则不予采信。5.3 一些可以复用的小技巧最后说我个人的三个习惯。第一个是每次写上手写函数之前先造两个完全无关的信号做“零对照”再制造一个只在某个窄频带共同相关、在其它频带无关的信号做“正对照”曲线形态对了才接真实数据。这个流程能省掉大量排查时间。第二个是保存中间量。在自写函数里我会把 Sxx、Syy、Sxy 都作为可选输出保留下来。表面上多占了内存但后面画对数功率谱、看互谱相位时可以直接复用这三个量不用重新算一遍分段 FFT。第三个小技巧是在旋转机械现场采集时尽量让两个传感器布置在相同振源传递路径上。比如轴承座的两个螺栓位置或者同一条轴系的两端。通道间传播路径差异过大的话共同振源的相位会严重扭曲相干峰即便存在也会被频响特性改变形状。这不是代码能弥补的必须从测点布置时就考虑进去。谱相干这个工具最吸引我的地方是它总能在“看起来都有噪声”的数据里告诉你哪些成分是共同且稳定的。我会在项目一开始就把它定位成“统计筛选器”先用它找频带、找激励源再针对找到的频带做深入量化诊断。每次准备做一套新的旋转机械测试或一组多通道阵列测量时我的固定习惯是先跑一遍谱相干再做其它精细分析相当于先给数据画一幅地图。