小波交叉功率谱与相位分析：MATLAB实现全解析

1. 为什么偏偏是"小波交叉功率谱"——当"两路信号在哪个频段相关"答不上来时

做信号处理的朋友应该都遇到过这种尴尬：手头有两路数据，明显觉得它们之间有关系，但传统方法一个值根本说不清楚。

我去年处理两路水声信号时就撞上了这个问题。一路是目标辐射噪声的包络，一路是拖曳阵某个通道的输出。按理说目标信号进来了，两路之间应该有强相关；但真正拿起互相关函数一看，相关峰被各种干扰压得根本看不清，再用FFT估计相干系数，也只能得到一条"平均意义"上的曲线，告诉我整个时间段里哪个频段相关性强。可问题是，目标信号是时有时无的，上一秒还很强的相关性，下一秒就被环境噪声盖掉了。这种非平稳场景下，全局相干性给出的结论几乎没法用。

后来我换了思路，把两路信号都做连续小波变换，再在同一个时频网格上做交叉功率谱。效果立刻不一样了：横轴是时间，纵轴是频率，颜色表示两路信号在这个时刻这个频带的共同能量，箭头表示相位关系。哪一秒开始相关、持续多久、在哪个频段最强、相位滞后多少，全部一目了然。这就是小波交叉功率谱分析（Cross Wavelet Transform，XWT）的核心价值。

这篇笔记我把整套可直接跑的MATLAB源代码拆开讲一遍，从原理到代码实现，再到我自己踩过的坑，尽量让看完的人不只知道怎么按回车，还知道每行代码在算什么、算出来的图该怎么读。适合正在做故障诊断、水声信号处理、脑电信号分析、气象海洋数据分析的朋友参考，尤其是那些用传统相干方法算不出结果、正在找更精细时频分析手段的人。

2. 交叉小波谱的数学骨架：连续小波变换、交叉功率与相位箭头

想用好交叉小波谱，先得把三个概念理顺：单个信号的连续小波变换、两个小波变换的交叉乘积、以及交叉谱的显著性检验。缺一个，后面的图都容易读歪。

2.1 Morlet小波为什么是首选

连续小波变换本质上是用一族可伸缩平移的"小波"去匹配信号的局部特征。小波函数必须既在时域又频域有较好的局部性。常用的小波基里，complex Morlet小波是交叉谱分析的绝对主力，因为它是一个解析小波，只有正频率分量，相位信息完整，这直接决定了我们能从交叉谱里提取两路信号的相位差。

Morlet母小波的标准形式是：

ψ0(t) = π^(-1/4) · e^(jω0 t) · e^(-t²/2)

其中ω0是无量纲的中心频率，通常取6。这个取值不是随便定的：ω0=6时，小波的时频面积比较平衡，频率分辨率和时间分辨率都够用；如果ω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 Coherence，WTC），把交叉功率用两路信号各自的平滑功率归一化，得到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，可以把这个函数贴在脚本里，替换内置cwt：

function 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 三个自查指标

每次写完交叉谱程序，我习惯先跑仿真再判断算法是否正确：

1. 强能量必须出现在设定的共同频率时间窗内。
2. 相位箭头的方向必须与设定延迟一致，如果箭头方向乱得不成形，多半是复数共轭的顺序反了（记得wx .* conj(wy)，不是conj(wx).* wy）。
3. 独立成分处不能出现强交叉能量——如果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 Influence，COI）。

我在第一次画出交叉谱图时，看到图左右边缘有一整条窄窄的亮带，以为发现了强相关。后来一算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 批量处理与结果导出

批量处理时，我通常把多通道数据放在矩阵里循环调用函数，再把提取的特征存进表格。比如可以在循环里统计每个时频斑块的峰值频率、峰值时间、相位均值，最后得到一张表：

这样一个循环下来，几十组数据的分析就有了结构化的产出，后续做统计对比就非常方便。顺便强调一下：这个统计过程要和交叉谱图并行保留，图是给人看的，数据才是给报告用的。

6.3 自写代码与现成工具包怎么选

目前网上流传较广的是Grinsted等人在2004年发布的MATLAB小波相干工具包，集成度高、显著性检验完整，直接用确实省事。我的看法是：

• 如果你只需要一个最终结果图，项目周期又紧，直接用现成工具包没有任何问题。
• 如果你要分析的数据比较特殊（比如需要自定义小波基、需要非均匀时间轴、需要嵌入更大的处理流程），或者你就是想把原理彻底弄透，自己实现这一段比调黑盒更有价值。

我之所以坚持保留自写版本，就是因为有一次遇到非均匀重采样数据，现成工具包直接罢工，最后靠自写版本改了两行才解决问题。交叉谱分析说到底不是特别复杂的算法，几个关键环节掌握之后，灵活度远高于固定工具包。

最后分享一个个人习惯：每次跑交叉谱分析，我都会同时保留三样东西——原始信号、交叉谱矩阵、显著性阈值矩阵。这三样数据在，随时可以重画图、改配色、换频率范围，而不需要重新跑一遍计算。图可以临时画，中间数据丢了才真是欲哭无泪。