简介本资源是一套基于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 文件。我把常见情况归纳成了一张排查表报错信息大概率原因处理方式file is not a zip file下载到的其实是个 HTML 错误页或者下载中断用文件管理器看扩展名用文本编辑器打开看到html开头就说明下错了需要重新下载invalid zip archive: could not find EOCDzip 文件被截断缺少末尾目录记录使用zip -FF 损坏文件.zip --out 修复文件.zip尝试修复修复失败就重新下载解压后某个.m文件打不开解压软件可能把文件名编码搞乱了换一个解压工具Windows 上我建议用 7-Zip 或 Bandizip解压路径有中文或空格某些老版本 Matlab 对中文路径支持很差把解压目录放在D:\tools这类纯英文路径下如果你是在 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越小时间定位越准但频率带宽会变宽。w06是一个折中的经典选择此时小波尺度与傅里叶周期近似相等解释起来很符合直觉。如果你发现某个周期性成分在图上被拉得很宽、无法区分可以尝试把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参数决定是否对数据末尾做零填充默认pad1也就是填充目的是让 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 时nsim200结果已经稳定。数据点 2000 以上nsim建议至少 500否则显著性区域边界会轻微抖动。调参阶段先用nsim20看趋势最后出图再用nsim1000。另外为了结果可复现调用前设置随机数种子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。但老代码的循环变量里可能包含随机数生成、数组索引需要先测试。我的建议是先跑一次nsim20验证结果再改成并行。如果你不希望改动原函数也可以降低默认密度先nsim100出一张初稿确认周期尺度、相位方向和显著性区域都符合预期后再用nsim1000出终稿。这样既能保证质量又能节省大量等待时间。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 里跑一遍自带的示例数据。跑通了再换自己数据换数据时先画单序列小波功率谱确认数据本身没有异常再去做相干分析。这套流程我是吃了不少亏才总结出来的照着走能省的不只是两天时间还有对着空白图发呆的无数个夜晚。本文还有配套的精品资源点击获取