小波相干分析Matlab工具箱使用教程：从解压到跑通并读懂结果

简介：本资源是一套基于Matlab实现的小波相干性（Wavelet Coherence）分析完整代码包，面向本科及硕士阶段的信号处理、地球物理、气候时序分析等科研学习者，解决多变量非平稳时间序列间局部相关性与相位关系的可视化建模问题。压缩包共109个文件，包含35个核心.m函数脚本（含主分析流程、绘图封装与参数配置）、15份Markdown格式说明文档（涵盖算法原理、参数释义与典型用例）、12个文本配置文件及11个实测数据样本（如sst_nino3.dat），辅以HTML报告模板、PNG结果图与PDF参考文献，整体体积仅3.08MB，结构清晰、即装即用。已有140人下载学习，所有代码兼容Matlab 2014a/2019a，内附可直接运行的示例及对应结果图，特别适合初学者理解小波相干谱的计算逻辑与图形解读，并为后续扩展至神经网络预测、路径规划等跨领域时序建模提供底层工具支撑。 上周翻资料的时候，又看到那个躺在我硬盘里的wavelet-coherence matlab代码.zip。说实话，第一次拿到这个压缩包时我差点被劝退：没有说明文档，没有版本号，里面一堆.m文件和.dat文件，文件名还都是缩写。后来花了整整两天，把它从解压到跑通、再摸清参数，才发现这其实是目前分析两个时间序列耦合关系最实用的工具之一。这篇文章就把我从零开始折腾这套代码包的完整过程写下来，包括解压报错怎么处理、核心函数分别干什么、参数怎么调、结果图怎么读，以及我踩过的几个典型坑。如果你也刚好在找小波相干性的 Matlab 实现，或者下载了同名代码包还没跑通，这篇应该能帮你省下不少时间。

1. 为什么需要 wavelet-coherence：从相关到小波相干

很多做时间序列的人一开始都会有这个疑问：我明明会用 Pearson 相关，为什么还要折腾小波相干？这个问题的答案，恰恰是这套代码存在的前提。

1.1 相干性分析要解决什么问题

Pearson 相关算出来是一个数字，比如 0.6，它描述的是两个变量在整个时间段内的线性相关强度。但现实中的数据往往没这么老实。以气象数据为例，气温和降水之间可能存在“某几年相关很强、某几年又几乎无关”的现象，或者只在特定周期尺度上相关。比如两个站点之间的降水耦合，可能只对 2 到 4 年的周期显著，到了 8 年以上的尺度就完全脱钩。传统的相关分析会把所有这些信息压成一个平均值，等于把一个西瓜榨成汁，你喝到的是混合味，却分不清里面有哪些果肉。

小波相干（wavelet coherence）解决的就是这个问题：它在“时间—周期”二维平面上逐点计算两个序列的局部相关关系。横轴是时间，纵轴是周期（或者频率），颜色深浅代表相关强度。这样既能看出关系在什么时候出现、什么时候消失，也能分辨出这种关系主要体现在哪个周期带上。

这类分析最早在地球物理和气候领域用得最多，后来金融、脑电信号、机械故障诊断也都在用。只要你有两个等间隔采样的时间序列，想知道它们之间的耦合是否随时间变化、是否存在滞后，就可以用这套工具。

1.2 小波变换、交叉小波与小波相干的区别

很多新手会把连续小波变换（CWT）、交叉小波变换（XWT）和小波相干（WTC）混在一起。其实三层是递进的关系。

• 连续小波变换（CWT）：对单个序列做时频分解，得到的是“某个序列的能量在时间和频率上如何分布”。它回答的是“这个序列自身有哪些周期成分”。
• 交叉小波变换（XWT）：对两个序列的 CWT 结果做乘积，得到的是“两个序列共同的高能量区域”。它回答的是“两个序列在什么时候、什么尺度上同时有很大的功率”。
• 小波相干（WTC）：对交叉小波谱做归一化处理，本质上是在时频平面上算“局部相关系数”。它回答的是“就算能量不大，两个序列的波动模式是否也高度同步”。

有一个很关键的点：交叉小波的能量高，不代表两个序列一定强相关；小波相干值接近 1，也不代表两个序列的绝对振幅一致。XWT 关注的是共同能量，WTC 关注的是波形一致性。实际分析时，我会把这两张图并排看，才能得到完整判断。

1.3 这套代码包的定位

从函数命名和代码风格来看，这个 zip 里装的应该是经典的 Grinsted 等人开发的小波相干工具箱，核心算法基于 Torrence 和 Compo 的连续小波变换实现。这个工具箱在地球科学领域几乎算事实标准，后来也被移植到 Python 等语言里。它的优势很直观：内置了基于 AR(1) 红噪声假设的显著性检验，可以直接输出带置信区间的图，甚至还能画出相位箭头，告诉我们两个序列之间的滞后方向。

所以，如果你只是想快速得到一张“能写在论文里”的小波相干图，这个代码包比从零写一个要靠谱得多。后面我讲的步骤，都是针对这个经典版本。

2. 代码包内部结构与核心函数角色

解压之后先别急着运行，花十分钟把文件结构看清楚，后面能少很多莫名其妙的报错。这个 zip 里的文件通常不长，但每个文件都有自己的分工。

2.1 解压后你会看到什么

一个典型的小波相干 Matlab 工具箱解压后会有以下这些文件：

• wt.m：连续小波变换函数，返回单个序列的小波功率谱、尺度、周期和 COI。
• xwt.m：交叉小波变换函数，返回两个序列的交叉小波谱和相位角。
• wtc.m：小波相干函数，返回相干系数 Rsq、周期、尺度、COI 和相位角，也是我们最常用的入口。
• ar1.m：用最大似然或最小二乘估计序列的 AR(1) 滞后自相关系数。
• rednoise.m：生成红噪声替代序列，供蒙特卡洛显著性检验使用。
• wavelet.m：底层小波母函数实现，包含 Morlet、Paul、DOG 等选项。
• bight.m：绘图中用来生成黑线等高线的辅助函数，缺了它画不出置信区间边界。
• 示例数据文件：通常是.mat或.txt，用来让你快速跑通整个流程。
• 可能还有README.txt或license.txt，有必要先扫一眼。

这些文件之间存在依赖关系。你如果只拷贝了wtc.m而没有ar1.m和wavelet.m，一运行就会报Undefined function。所以最稳妥的做法是保持压缩包内的目录结构，整个文件夹都放进 Matlab 路径。

2.2 核心函数 wt、xwt、wtc 分别承担什么任务

我自己在使用时，最核心的依赖关系是这样的：

• 如果你只想看单个序列的周期成分，调用wt(x)，它会先对序列做连续小波变换，然后返回小波功率谱power、尺度数组scale、傅里叶周期period和影响锥coi。

[power, period, scale, coi] = wt(x, 'dt', dt);

• 如果你想知道两个序列的共同高能区域，调用xwt(x, y)，它会先分别对 x 和 y 做 CWT，再计算交叉谱Wxy和相位差phase。

[Wxy, period, scale, coi, phase] = xwt(x, y, 'dt', dt);

• 如果你想看局部相关性，也就是绝大多数论文里那个红蓝相间的图，调用wtc(x, y)，它内部会做交叉谱的平滑处理，然后计算相干系数。

[Rsq, period, scale, coi, phase] = wtc(x, y, 'dt', dt, 'nsim', 200);

这里我要特别强调一下wtc.m的默认行为：如果不输入输出参数直接运行wtc(x, y)，它会在计算完成后画出一张完整的图；如果指定了输出参数，则只返回数据，不自动画图。这是很多人第一次使用时容易懵的地方——直接运行有图，一加上分号保存数据反而没图了。

2.3 数据输入格式与约定

这个代码包的输入约定说简单也简单，说严格也严格。函数接受的x和y必须是列向量，也就是size(x)应该是[N, 1]。如果你导入的数据是行向量，比如从表格里复制出来是[1, N]，不转换直接跑，大概率会报矩阵维度不匹配的错误，或者得到一张完全空白的图。

另一个硬性要求是等间隔采样。小波变换本身不关心时间轴上的具体日期，它只依赖等间隔假设和采样时间间隔dt。如果你的数据缺了几个月，别直接把缺的那行删掉就完事。删除会让后续数据提前，等效于改变了时间轴，低频周期会算错。正确的做法是先插值补齐，再进入分析流程。

此外，NaN 是这个工具箱的“天敌”。wt.m内部做 Fourier 变换时遇到 NaN 会直接把整条谱变成 NaN。所以拿到数据后，第一件事就是检查缺失情况，用isnan定位并处理，不要指望它能自动跳过缺失值。

3. 从下载到跑通第一张图的完整步骤

很多人下载这个 zip 之后，第一步就卡住了：压缩包解不开。这不是玩笑，我见过大量求助帖都是同一个标题：file is not a zip file或者invalid zip archive: could not find EOCD。我把从解压到跑通第一张图的完整路径写在这里。

3.1 正确解压 zip 文件的姿势与报错排查

先说你最可能遇到的解压问题。如果你在 Windows 上双击压缩包时提示“文件已损坏”或“压缩包无效”，多半不是压缩包真的坏了，而是它压根不是一个完整的 zip 文件。

我把常见情况归纳成了一张排查表：

如果你是在 Linux 下操作，命令也很常规：

unzip wavelet-coherence-matlab.zip

如果unzip提示End-of-central-directory signature not found，先运行：

file wavelet-coherence-matlab.zip

看看它到底是 zip 还是 HTML 或纯文本。看输出再决定是重新下载还是修文件，不要盲目重复试。

3.2 添加路径与依赖关系

解压完成后，打开 Matlab，把当前目录切换到解压后的文件夹，然后把它加入搜索路径：

cd('D:\tools\wavelet-coherence-matlab') addpath(genpath(pwd));

这里用genpath(pwd)而不是简单的addpath(pwd)，是因为工具箱内部可能还有子目录，递归添加能避免找不到依赖函数的问题。

添加完路径，建议先确认关键函数能被找到：

which wtc which wavelet which rednoise

只要这三个都不报错，说明基础依赖齐全。如果提示某个函数未找到，回到压缩包检查是不是完整解压，或者从其它发布渠道补一个同名函数回来。

3.3 用内置示例跑通第一张图

接下来我想强烈建议：先用工具箱自带的示例数据跑通，再替换成自己的数据。很多人习惯下载完直接用自己数据测，一旦出错，根本分不清是数据格式问题还是代码本身问题。

示例数据可能是.txt也可能直接是.mat。如果是文本文件，用load或readmatrix读进来：

data = readmatrix('example_data.txt'); x = data(:, 1); y = data(:, 2); dt = 1; % 根据数据采样频率调整，月数据写 1/12，年数据写 1

然后直接调用：

wtc(x, y, 'dt', dt);

如果一切正常，Matlab 会弹出一张小波相干图：横轴是时间索引，纵轴是周期（对数刻度），颜色表示相干系数，还有一个半透明的锥形区域。看到这张图，你的代码路径就彻底跑通了。

3.4 跑通过程中常见的运行报错

我把自己在跑通流程里遇到的高频报错整理一下：

• Undefined function 'bight'：这是绘图辅助函数缺失，常见于从某个网站单独下载的代码包被精简过。解决办法是找完整版本，补上bight.m。
• Error using * inner matrix dimensions must agree：数据方向不对。用size(x)检查，如果是[1, N]，就改成x = x(:)。
• Error using movmean或类似函数不存在：你的 Matlab 版本太老，部分平滑函数是老版本没有的。可以手动实现滑动平均，或者升级到新一点版本。
• 运行很久不出图：蒙特卡洛模拟默认次数是 200，数据点几千个时确实需要时间。先用短序列测试，或者把nsim临时改小到 20。

4. 参数选择与自定义：决定结果的是这些细节

代码跑通只是第一步，真正影响论文结论的是参数设置。同样是两个序列，不同母小波、不同尺度范围、不同显著性检验设置，得出来的图可能有肉眼可见的差别。下面几个参数是我每次使用都必须检查的。

4.1 母小波选择：为什么默认用 Morlet

工具箱支持多种母小波，最常见的是 Morlet，实际代码里的默认设置通常也是'mother', 'Morlet'。Morlet 由一个复正弦波乘以高斯窗组成，它的好处是能同时提取到幅度和相位信息，所以小波相干图中的相位箭头只有选择复数小波才能画出来。

Morlet 有一个关键参数是波数w0，官方代码里默认为 6。这个值决定了小波的振荡次数：w0越大，频率分辨率越高，时间分辨率越低；w0越小，时间定位越准，但频率带宽会变宽。w0=6是一个折中的经典选择，此时小波尺度与傅里叶周期近似相等，解释起来很符合直觉。

如果你发现某个周期性成分在图上被拉得很宽、无法区分，可以尝试把w0调大到 8 或 10；反过来，如果你更关心某个突变事件的时间位置，可以把w0调小。注意w0不是代码里的直接参数，需要在小波函数定义里改，或者查wtc.m是否支持'mother'参数的扩展字段。

4.2 尺度范围与 pad 设置

在小波变换中，“尺度”约等于频率的倒数。代码默认会从最小尺度s0开始，按指数增长生成一系列尺度，直到覆盖数据允许的最大尺度。默认设置通常能自动适应数据长度，但有时自动生成的尺度范围太窄，导致低频区域被过多 COI 覆盖。

可以手动控制尺度范围：

[Rsq, period, scale, coi, phase] = wtc(x, y, 'dt', dt, 's0', 2*dt, 'dj', 0.125, 'J', 40);

这里的s0是最小尺度，dj是相邻尺度间隔，J是总尺度数。dj越小，尺度网格越密，图越平滑，但计算量越大。通常0.125已经够用。

pad参数决定是否对数据末尾做零填充，默认pad=1，也就是填充，目的是让 FFT 长度变为 2 的幂次，加快计算。绝大多数情况下不需要改它。但要注意，零填充只是加快卷积计算，并不会改变边界效应，所以 COI 依然存在。

4.3 去趋势与 AR(1) 红噪声背景

这套显著性检验的默认零假设是：两个序列都是一阶自回归红噪声，也就是 AR(1) 过程。这个假设对很多气候、经济序列是合理的，因为它符合“上一期的状态会影响下一期”的特点。

不过，如果你的序列带有明显线性趋势，或者存在强烈的季节性周期，直接丢进wtc里会让低频段出现大面积“伪显著”。比如你有 30 年的月均温数据，整体升温趋势会让 8 年以上周期的能量巨大，但它们其实是趋势的产物，不是两个序列的真实耦合关系。

我的习惯是：先对序列做标准化和去趋势，再进入小波相干分析。比如用内置detrend去掉线性项，再减去均值除以标准差：

x = detrend(x); x = (x - mean(x)) / std(x); y = detrend(y); y = (y - mean(y)) / std(y);

ar1.m函数会自动估计序列的滞后一阶自相关系数，并用于生成红噪声模拟序列，你不用手动传入太多东西。但如果你的数据包含强周期成分，建议先做季节差分或用滤波手段去掉周期，否则显著性检验会被“已知周期”污染。

4.4 显著性检验参数：蒙特卡洛模拟次数到底设多少

wtc.m的显著性检验是蒙特卡洛方法：生成很多对红噪声序列，计算它们的小波相干分布，然后用实际观测值跟分布比较。默认的模拟次数在经典版本里是 200 或者由nsim指定。

nsim越大，显著性边界越稳定，但耗时线性增长。我的实际经验是：

• 数据点少于 500 时，nsim=200结果已经稳定。
• 数据点 2000 以上，nsim建议至少 500，否则显著性区域边界会轻微抖动。
• 调参阶段先用nsim=20看趋势，最后出图再用nsim=1000。

另外，为了结果可复现，调用前设置随机数种子：

rng(1);

这样别人运行你的代码能得到完全相同的显著性区域，而不是每次重新抽样都有一点点不同。这一点在学术交流中很重要。

5. 读懂输出图：从颜色、箭头到置信区域

跑出了图，别急着放论文，先搞清楚图里的每个元素是什么意思。小波相干图看起来花花绿绿，但真正要读的信息就那么几块。

5.1 小波功率谱图怎么看

虽然我们要的是相干图，但很多时候第一步应该先看每个序列的单序列小波功率谱。功率谱用颜色表示强度，通常暖色代表能量高，冷色代表能量低。横轴是时间，纵轴是周期，周期轴一般用对数刻度，所以图下方是高分辨率、图上方是低分辨率。

图上有两个关键标注：一是黑色粗线圈出来的是 5% 显著性水平区域，代表在这个时空区域里，序列的能量不是红噪声随机波动能解释的；二是底部弧线围起来的半透明区域，叫影响锥 COI，表示受边界影响不可靠的区域。

对于小波功率谱，我通常会先找“哪个周期带最稳定、能量最强”。比如某站点降水在 2 到 4 年周期带上持续有显著能量，说明该地降水存在明显的年际变率。这个结论直接决定后续相干分析该关注哪些周期段。

5.2 交叉小波与相干性图的区别

交叉小波图 XWT 和相干图 WTC 经常被放一起，但它们的含义完全不同。XWT 显示共同高能量区域，适合找“两个序列同时强振荡”的时段；WTC 显示局部相关强度，适合找“两个序列同步变化”的时段，哪怕振荡幅度很小。

我举一个实际例子：A 序列振幅很大，B 序列振幅很小，但两者形态非常一致，只是比例不同。这种情况下 XWT 可能在 A 能量高的区域有一点信号，但整体上共同能量不高；WTC 则可能非常高，因为它比的是变化模式，不是振幅。

反过来，如果两个序列在某个时段都有巨大能量，但波形来自不同周期，XWT 会显示高能，WTC 反而很低。所以我在报告结论时一定同时给 XWT 和 WTC，避免只看单一指标得出偏颇判断。

5.3 相位箭头与滞后关系

小波相干图里的小箭头是相位差信息，也是这套代码包最迷人的地方。箭头方向由交叉谱的相位角决定。常见的约定是：

• 箭头指向右方：两个序列在该时间尺度上同相，呈正相关。
• 箭头指向左方：两个序列反相，呈负相关。
• 箭头指向上方：第一个序列（x）领先第二个序列（y）约四分之一个周期。
• 箭头指向下方：第二个序列（y）领先第一个序列（x）约四分之一个周期。

我这里用了“常见约定”，因为不同代码实现的相位角可能基于atan2(imag(Wxy), real(Wxy))，也可能互换参数顺序。严谨的做法是在使用前打开xwt.m或wtc.m，看一眼相位角是怎么算的。绝大多数 Grinsted 版本遵循上述约定，但你替换数据后最好用一组已知相位差的数据验证一下。

如果图上的箭头方向很乱，不要强行解释。我一般只在显著性区域内统计平均相位方向，再结合领域的物理含义判断谁驱动谁。箭头只有落在黑线圈住的范围里才有讨论价值，黑圈外的箭头可能是随机噪声造成的。

5.4 锥形影响区 COI 的意义

COI 是“cone of influence”的缩写，翻译成影响锥。小波变换需要对有限数据两端做截断处理，所以靠近首尾的谱值会混入大量边界假信号。COI 就是这样一个边界分界线：COI 内部的能量和相干大概率受到边界污染，不能作为严格结论的证据。

识别方法很简单：看图中曲线阴影锥形区域。周期越大，边界影响越深入数据内部。比如你用 10 年逐月数据，周期为 8 年的信号几乎整条都被 COI 覆盖，因此不能声称你发现了 8 年的显著耦合。

处理 COI 的策略通常有三个：

• 数据尽量取长一些，让感兴趣的周期段完整落在 COI 外部。
• 分析时把结论集中在 COI 以外的时空区域。
• 如果必须讨论低频，建议用更长时段的数据作为补充证据。

个人经验：不要试图通过调小pad来消除 COI，因为边界效应来自小波本身，不是零填充导致的。COI 是诚实的警告灯，硬把它去掉等于自欺欺人。

6. 避坑记录：我在使用这个代码包时踩过的坑

最后这部分是纯经验，全是实际操作中撞出来的。如果你能把前面的步骤都跑通，大概率不会再犯我犯过的错，但还是有几个坑值得单独强调。

6.1 数据长度和缺失值处理

我第一次用这个代码包时，拿着 15 年的逐月数据（180 个点）去算，结果低频周期 5 年以上全在 COI 里，能看的只有 1 到 5 年周期段。这不是代码的问题，而是数据长度在物理上就无法支撑低频段的分析。小波变换的频率分辨率有限，数据越短，低频可分辨性越差。

对于缺失值，我的建议是不要偷懒用均值填充。均值填充会在填充点附近制造平坦段，破坏高频成分，导致小波功率谱出现假的低频能量。更好的做法是：

x = fillmissing(x, 'linear');

如果缺失段太长，比如连续缺了 20% 以上，直接用插值已经不可靠，最好换成不需要完整序列的方法，或者截取连续完整的子段分析。

6.2 版本兼容问题：老代码遇上新 Matlab

这套代码最早的版本是 2004 年前后写出来的，Matlab 这些年更新了很多绘图和数值特性。我在 R2022b 上运行，遇到过两个现象：一是默认配色从jet变成了parula，导致小波图颜色风格和其他论文对不上；二是个别老函数调用方式会在命令行刷 warning，但不影响结果。

简单解决方式：

colormap(jet); % 恢复经典配色 warning off; % 出图前关掉警告，眼不见心不烦

如果遇到Subscript indices must either be real positive integers or logicals，别急着改代码，先检查数据里有没有 NaN 或 0。这类错误通常是数据传进出错，不是函数本身坏了。另外，在 Linux 虚拟机里跑 Matlab 会比较慢，尤其是蒙特卡洛模拟，必要时减小nsim或者用原生 Linux 版 Matlab，而不是虚拟机。

6.3 计算时间与并行优化

当你把nsim调大，数据点又长，运行时间可能从几十秒变成十几分钟。如果只是生成一张图，等几分钟还能接受；但如果要批量分析多个站点组合，就必须优化。

最快的优化方式是把wtc.m里的蒙特卡洛随机模拟循环改成parfor。但老代码的循环变量里可能包含随机数生成、数组索引，需要先测试。我的建议是先跑一次nsim=20验证结果，再改成并行。

如果你不希望改动原函数，也可以降低默认密度：先nsim=100出一张初稿，确认周期尺度、相位方向和显著性区域都符合预期后，再用nsim=1000出终稿。这样既能保证质量，又能节省大量等待时间。

6.4 结果批量导出与重绘

默认生成的图其实已经不错，但如果你想统一风格，或者想给多个子图拼到一张大图里，最好使用返回值自己画。常用片段：

[Rsq, period, scale, coi, phase] = wtc(x, y, 'dt', dt, 'nsim', 200); figure('Color', 'w'); imagesc(t, log2(period), Rsq); set(gca, 'YDir', 'reverse', 'YTick', log2(period(1:4:end)), 'YTickLabel', period(1:4:end)); colormap(jet); colorbar;

需要注意，自己画图时要把周期轴设成对数刻度，并注意YDir方向。默认imagesc的纵轴是从上往下递增，而小波图的周期轴一般是低周期在下、高周期在上，所以要反向。

批量导出时，我一般用exportgraphics而不是saveas，因为前者能保证 300dpi 分辨率以及字体嵌入：

exportgraphics(gcf, 'wtc_result.png', 'Resolution', 300);

我个人的习惯是，所有绘图参数都写在一个独立的脚本里，只改数据文件名，就能批量出图。这样可以避免每次手动调整坐标轴范围。

最后再分享一个可能帮你少走弯路的小技巧：拿到这种历史悠久的 Matlab 代码包，第一步不是看代码，而是先在 Matlab 里跑一遍自带的示例数据。跑通了再换自己数据，换数据时先画单序列小波功率谱，确认数据本身没有异常，再去做相干分析。这套流程我是吃了不少亏才总结出来的，照着走，能省的不只是两天时间，还有对着空白图发呆的无数个夜晚。

本文还有配套的精品资源，点击获取