1. 为什么偏偏是小波交叉功率谱——当两路信号在哪个频段相关答不上来时做信号处理的朋友应该都遇到过这种尴尬手头有两路数据明显觉得它们之间有关系但传统方法一个值根本说不清楚。我去年处理两路水声信号时就撞上了这个问题。一路是目标辐射噪声的包络一路是拖曳阵某个通道的输出。按理说目标信号进来了两路之间应该有强相关但真正拿起互相关函数一看相关峰被各种干扰压得根本看不清再用FFT估计相干系数也只能得到一条平均意义上的曲线告诉我整个时间段里哪个频段相关性强。可问题是目标信号是时有时无的上一秒还很强的相关性下一秒就被环境噪声盖掉了。这种非平稳场景下全局相干性给出的结论几乎没法用。后来我换了思路把两路信号都做连续小波变换再在同一个时频网格上做交叉功率谱。效果立刻不一样了横轴是时间纵轴是频率颜色表示两路信号在这个时刻这个频带的共同能量箭头表示相位关系。哪一秒开始相关、持续多久、在哪个频段最强、相位滞后多少全部一目了然。这就是小波交叉功率谱分析Cross Wavelet TransformXWT的核心价值。这篇笔记我把整套可直接跑的MATLAB源代码拆开讲一遍从原理到代码实现再到我自己踩过的坑尽量让看完的人不只知道怎么按回车还知道每行代码在算什么、算出来的图该怎么读。适合正在做故障诊断、水声信号处理、脑电信号分析、气象海洋数据分析的朋友参考尤其是那些用传统相干方法算不出结果、正在找更精细时频分析手段的人。2. 交叉小波谱的数学骨架连续小波变换、交叉功率与相位箭头想用好交叉小波谱先得把三个概念理顺单个信号的连续小波变换、两个小波变换的交叉乘积、以及交叉谱的显著性检验。缺一个后面的图都容易读歪。2.1 Morlet小波为什么是首选连续小波变换本质上是用一族可伸缩平移的小波去匹配信号的局部特征。小波函数必须既在时域又频域有较好的局部性。常用的小波基里complex Morlet小波是交叉谱分析的绝对主力因为它是一个解析小波只有正频率分量相位信息完整这直接决定了我们能从交叉谱里提取两路信号的相位差。Morlet母小波的标准形式是ψ0(t) π^(-1/4) · e^(jω0 t) · e^(-t²/2)其中ω0是无量纲的中心频率通常取6。这个取值不是随便定的ω06时小波的时频面积比较平衡频率分辨率和时间分辨率都够用如果ω0取得更大波形更像正弦频率分辨率高但时间分辨率差边缘效应也更明显。交叉谱分析里我一般不轻易改这个参数默认6用到底。对信号x做连续小波变换得到的是二维复系数矩阵Wx(a,b)a是尺度参数对应频率b是平移参数对应时间。小波功率谱定义为|Wx|²。因为小波系数是复值所以它天然携带相位信息这就给后面的交叉相位分析留下了入口。2.2 交叉谱与相位箭头颜色只是能量箭头才是关系两路信号x和y各自做完连续小波变换得到复系数Wx和Wy后交叉小波谱的定义非常简洁Wxy Wx · conj(Wy)这里conj是取共轭。交叉小波功率就是|Wxy|。如果某个时频点上两路信号都有明显能量且相位一致|Wxy|就会很大如果两路信号在该处只有一个有能量或相位关系混乱交叉功率就小。而这个复数的幅角就是相位差Φxy arg(Wxy)图中用箭头表示箭头水平向右表示两路信号在该时频点同相垂直向上表示y信号超前x信号90°垂直向下表示滞后90°箭头方向随相位差连续旋转。这个信息对分析因果方向、传播延迟特别有用。但交叉谱有一个天然短板——它没有归一化的量纲输出数值大小同时受两路信号自身幅值影响。如果x一路信号幅值很大、y幅值很小交叉功率会虚高误导判断。这时候需要引入小波相干谱Wavelet CoherenceWTC把交叉功率用两路信号各自的平滑功率归一化得到0到1之间的相干系数。WTC本质上就是时频域的动态相关系数比XWT更适合比较不同幅值的信号。2.3 显著性检验别把噪声的偶然相关当结论交叉谱很容易犯一个错误两路独立的白噪声随机信号算出来的交叉功率也不是零而且图上看还挺热闹。如果直接上色标图这些随机起伏会被当成强相关来解读。所以标准的交叉谱分析必须做显著性检验。常见做法基于红噪声AR(1)过程背景假设先估计两路信号各自的lag-1自相关系数构造理论红噪声谱再依据卡方分布给出95%置信线。也可以更粗暴地用蒙特卡洛方法生成很多对不相关的随机信号统计交叉功率的95%分位数大于这个阈值的时频点才算显著。一句话总结显著性检验的意义是把偶然出现的相关和真实的共同振荡区分开。我第一次画交叉谱图时没做显著性检验图上大片黄色区域差点得出错误结论后来叠加上置信线真正显著的区域就缩成了几个清晰的斑块。这个环节千万别省。3. MATLAB源代码从零实现连续小波交叉谱如果只用现成工具箱很多人会忽略中间细节出了问题也不知道是算法问题还是参数问题。我下面的实现分两个版本主版本用MATLAB的Wavelet Toolboxcwt函数简洁、稳定适合日常使用备选版本给了一个手写的Morlet小波变换函数没有工具箱的时候可以直接替换。3.1 主程序完整可跑的MATLAB交叉谱脚本%% 小波交叉功率谱分析 - 主程序 % 功能计算两路信号的小波交叉谱、相位、并绘图 % 适用非平稳两通道信号观察何时、何频带存在共同振荡 clear; clc; close all; rng(2025); %% ---------- 1. 构造测试信号 ---------- Fs 256; % 采样率 256 Hz dt 1 / Fs; t 0:dt:4; % 4 秒数据 N numel(t); x zeros(1, N); y zeros(1, N); % 共同分量1~3秒内 8 Hz 正弦y 相比 x 滞后 45 度 idx t 1 t 3; x(idx) x(idx) sin(2*pi*8*t(idx)); y(idx) y(idx) 1.2*sin(2*pi*8*t(idx) - pi/4); % 各自独立分量x 在 0.5~2s 有 25Hz 成分y 有 3Hz 长时成分 seg t 0.5 t 2; x(seg) x(seg) 0.8*sin(2*pi*25*t(seg)); y y 0.6*sin(2*pi*3*t); % 加一点噪声模拟真实环境 x x 0.05*randn(1, N); y y 0.05*randn(1, N); %% ---------- 2. 连续小波变换 ---------- % 使用复Morlet小波amor返回复值系数矩阵 [wx, freqs] cwt(x, amor, Fs); [wy, ~] cwt(y, amor, Fs); % 如果需要更密的频率轴可以设置每倍频程小波数 % [wx, freqs] cwt(x, amor, Fs, VoicesPerOctave, 16); %% ---------- 3. 交叉小波功率与相位 ---------- wxy wx .* conj(wy); % 交叉小波谱复数矩阵 xwt_power abs(wxy).^2; % 交叉小波功率能量形式 phase angle(wxy); % 相位差单位 rad % 想看归一化相干谱时需要先对交叉谱/功率做平滑 % 这里直接输出XWT结果WTC在5.3节单独说明 %% ---------- 4. 可视化 ---------- figure(Color, w, Position, [200 200 980 560]); % pcolor 配合 log 频率轴比 imagesc 更准确 pcolor(t, freqs, xwt_power); shading interp; % 平滑着色 set(gca, YScale, log, YDir, normal); ylim([min(freqs) max(freqs)]); xlim([t(1) t(end)]); xlabel(时间 (s), FontSize, 12); ylabel(频率 (Hz), FontSize, 12); title(交叉小波功率谱 (XWT), FontSize, 14); colormap(parula); cb colorbar; ylabel(cb, 交叉功率); hold on; %% ---------- 5. 叠加相位箭头 ---------- % 箭头不能画太密否则全是黑线我一般取时间步长为总时长的1/20 % 频率步长为总行数的1/12 tStep round(N / 20); fStep round(size(wx, 1) / 12); [tg, fg] meshgrid(1:tStep:N, 1:fStep:size(wx, 1)); ph phase(fg, tg); % 注意行列索引 % 用 quiver 画箭头缩放系数0.5让箭头不互相压盖 quiver(t(tg), freqs(fg), 0.5*cos(ph), 0.5*sin(ph), 0, ... w, LineWidth, 1.1);这段跑完会得到一张时频交叉功率图图上叠加了相位箭头。cwt函数返回的freqs是自然对数间隔的频率向量所以纵轴用log刻度能明显改善低频段的分辨率。我强烈建议不要用imagesc硬画因为imagesc假定像素等间距在log频率轴上会把低频成分压得看不见。3.2 备选方案没有工具箱时手写Morlet小波变换如果你用的MATLAB没有Wavelet Toolbox可以把这个函数贴在脚本里替换内置cwtfunction wt my_morlet_cwt(x, Fs, freq, w0) % 手写Morlet连续小波变换频域实现 % 输入 % x : 输入信号行向量 % Fs : 采样频率 % freq : 要分析的中心频率向量 % w0 : Morlet小波带宽参数默认6 % 输出 % wt : 复小波系数矩阵size length(freq) x length(x) if nargin 4, w0 6; end x x(:).; n numel(x); N 2^(nextpow2(n) 1); % FFT点数 X fft(x, N); f (0:N-1) * Fs / N; wt zeros(numel(freq), n); for k 1:numel(freq) % 尺度与中心频率的关系 s w0 / (2*pi*freq(k)); % 小波函数的频域表示只保留正频率 psi exp(-0.5 * (s*2*pi*f - w0).^2); psi(f 0) 0; % 线性卷积由频域乘法完成sqrt(s)保证能量归一化 wt(k, :) ifft(X .* conj(psi) .* sqrt(s), symmetric); end wt wt(:, 1:n); end使用方式很简单[wx, freqs] my_morlet_cwt(x, Fs, 1:0.2:50)然后就进入了交叉谱计算流程。手写版的优势是频率范围完全自控缺点是速度比工具箱慢不少数据量大时要耐心。4. 仿真验证用已知信号检验你的交叉谱程序对不对代码写完之后别急着往真实数据上套先用一组已知答案的仿真信号验证。这部分工作花不了十分钟却能避免后面一整天对着错误图做无用功。4.1 测试信号为什么这么设计上面主程序里的信号有三个关键点两路信号在1~3秒的8Hz处有共同分量且y相对x相位滞后π/4用来验证交叉谱能否在正确时频位置上给出强能量和正确相位。x有独立的25Hz短时成分y有独立的3Hz长时成分用来测试交叉谱会不会把各自有能量误判成共同能量。正确结果是8Hz处交叉能量强25Hz和3Hz处应当弱。加了少量白噪声检验程序在噪声干扰下是否还能找到目标斑块。这组设计覆盖了交叉谱最重要的性能指标时频定位准确度、相位恢复能力、抗干扰能力。4.2 读图颜色、位置、箭头分别说明什么理论上跑出来的交叉功率谱应该呈现一个清晰的亮斑中心落在横轴1~3s、纵轴8Hz处。这就是两路信号存在共同振荡的核心证据。相位箭头在这个亮斑内部应当方向一致指向右上方——因为在我们的约定里arg(wx) - arg(wy) 0表示y滞后于x。确实构造信号时给y加了-π/4相位延迟所以箭头应该稳定落在右上方45度附近。如果把鼠标移到亮斑中心用MATLAB的datatip工具读出精确坐标时间应约为2s频率应非常接近8Hz。这说明小波交叉谱的时频定位能力是可靠的。4.3 三个自查指标每次写完交叉谱程序我习惯先跑仿真再判断算法是否正确强能量必须出现在设定的共同频率时间窗内。相位箭头的方向必须与设定延迟一致如果箭头方向乱得不成形多半是复数共轭的顺序反了记得wx .* conj(wy)不是conj(wx).* wy。独立成分处不能出现强交叉能量——如果25Hz处也亮成一片说明两路信号在该处能量太强而背景谱检验缺失这时就要加显著性检验。这是我屡试不爽的验收流程。真实数据就是没法给出标准答案的仿真这关过了才有底气让程序去处理那些既没有标准频谱也没有标准相位的数据。5. 实际项目里最容易翻车的四个细节仿真跑通了下面这些是从真实数据处理中沉淀出来的经验。顺序大致按画图之前→画图→读图来排。5.1 频率范围与频率分辨率怎么选cwt默认自动选择频率范围但自动档不一定适合你的问题。比如信号主要能量集中在5Hz以下默认范围却把高频段占了半张图低频细节反而看不清。建议显式设置频率范围[wx, freqs] cwt(x, amor, Fs, FrequencyLimits, [1 60]);频率下限不能低于1Hz数据长度4秒时再低就是伪信号了上限到奈奎斯特频率的一半即可。想要频率分辨率更细腻加参数VoicesPerOctave, 16。我常用的组合是频宽1~60Hz、每倍频程12~16个voice兼顾分辨率和计算速度。注意频率范围过大小波变换的边界效应也会扩大尤其在低频段几乎整段都会被影响锥COI覆盖结果没有意义。5.2 边界效应COI图边缘的结果不可信连续小波变换在时间轴两端会因为没有足够的数据支撑而产生虚假能量。这个区域叫影响锥Cone of InfluenceCOI。我在第一次画出交叉谱图时看到图左右边缘有一整条窄窄的亮带以为发现了强相关。后来一算COI亮带完全陷在边界区域里根本不具备显著性。真实结论只有图中间那块亮斑。处理办法很直接画图时把COI区域叠加上去读图时只关注COI之外的时频点。如果想标注边界可以近似用Morlet的e-folding时间计算边界随时间的变化再画两条对称的白色虚线把中间可信区域框出来。严谨的程序里这一步不能省。5.3 显著性检验的两种落地路径我在2.3节说过显著性检验的必要性这里给两条可操作的路径。路径一是解析法基于AR(1)红噪声谱假设用Torrence和Compo的经典公式计算单路小波功率谱阈值交叉谱的显著性需要同时考虑两路的背景谱公式稍复杂。这个方法的优点是计算快缺点是假设信号符合AR(1)模型对很多真实的非平稳信号并不严格成立。路径二是蒙特卡洛法直接在MATLAB里造若干对随机信号% 蒙特卡洛显著性检验示意代码 nsim 200; sig95 zeros(size(xwt_power)); for k 1:nsim nx 0.1*randn(1, N); % 随机噪声 ny 0.1*randn(1, N); [wx0, ~] cwt(nx, amor, Fs); [wy0, ~] cwt(ny, amor, Fs); cp abs(wx0 .* conj(wy0)).^2; sig95 max(sig95, quantile(cp, 0.95, 2)); % 取95%分位数 end % 在图上叠加95%置信线 hold on; contour(t, freqs, xwt_power ./ sig95, [1 1], k, LineWidth, 2);蒙特卡洛法的前提是构造不相关的替代信号最简单就是白噪声如果你想更贴近真实背景可以用带AR(1)特性的随机过程生成。这个方法的优点是直观且不依赖理论分布缺点是仿真次数越多计算越慢。实际项目中我常把两种方法都跑一遍结论一致时才能放心。5.4 相位箭头的两个常见问题箭头太密会糊成一片太稀疏又看不出相位变化趋势。我通常按数据长度的1/20取时间步长频率方向按行数的1/10到1/12取步长然后把最外圈的箭头COI区域里的剔除只保留可信区域内的箭头。这样图面干净信息量也够。另一个问题是相位在-π到π边界处的跳跃。两路信号相位差接近180度时箭头方向会在左右之间剧烈跳变看着就像噪声。这不是程序错了而是相位角度的固有环状特性。处理时可以用unwrap把相位沿频率方向解卷绕或者干脆接受这种跳变分析时避开临界频段。6. 从脚本到工具把交叉谱分析封装成能复用的函数上面的主程序是一锤子买卖一次只能处理一组数据。实际项目中我往往要对几十组数据批量做分析所以最后把这段逻辑整理成了一个独立函数约定好输入输出后面所有项目都直接调它。6.1 函数接口设计建议function result xwt_analysis(x, y, Fs, varargin) % 小波交叉功率谱分析函数 % 输入: % x, y - 等长的两路信号 % Fs - 采样率 % 可选参数: FreqRange, [fmin fmax]; VoicesPerOctave, n; % ShowPlot, true/false; SigTest, mc/none % 输出: % result - 结构体包含 freqs, time, power, phase, sig95返回的result结构体包含后续所有画图、报告所需的字段既可以把result保存成.mat也可以把核心字段导出成CSV给其他工具用。% 调用示例 r xwt_analysis(x, y, 256, ... FreqRange, [1 60], VoicesPerOctave, 16, ... ShowPlot, true, SigTest, mc);这个封装过程的收益是一次写清、长期复用。调试只做一次后面批量处理几十组数据时一个循环就完了。6.2 批量处理与结果导出批量处理时我通常把多通道数据放在矩阵里循环调用函数再把提取的特征存进表格。比如可以在循环里统计每个时频斑块的峰值频率、峰值时间、相位均值最后得到一张表通道对峰值频率峰值时间峰值相位差显著性ch1-ch28.1 Hz2.0 s-42度通过ch1-ch3无无无未通过这样一个循环下来几十组数据的分析就有了结构化的产出后续做统计对比就非常方便。顺便强调一下这个统计过程要和交叉谱图并行保留图是给人看的数据才是给报告用的。6.3 自写代码与现成工具包怎么选目前网上流传较广的是Grinsted等人在2004年发布的MATLAB小波相干工具包集成度高、显著性检验完整直接用确实省事。我的看法是如果你只需要一个最终结果图项目周期又紧直接用现成工具包没有任何问题。如果你要分析的数据比较特殊比如需要自定义小波基、需要非均匀时间轴、需要嵌入更大的处理流程或者你就是想把原理彻底弄透自己实现这一段比调黑盒更有价值。我之所以坚持保留自写版本就是因为有一次遇到非均匀重采样数据现成工具包直接罢工最后靠自写版本改了两行才解决问题。交叉谱分析说到底不是特别复杂的算法几个关键环节掌握之后灵活度远高于固定工具包。最后分享一个个人习惯每次跑交叉谱分析我都会同时保留三样东西——原始信号、交叉谱矩阵、显著性阈值矩阵。这三样数据在随时可以重画图、改配色、换频率范围而不需要重新跑一遍计算。图可以临时画中间数据丢了才真是欲哭无泪。