简介面向信号处理与数据分析的小波分析MATLAB工具包集成交叉小波、小波相干、交叉谱与相关分析模块适用于研究两个信号在不同时间尺度上的同步性、相位关系与相关性尤其适合地震学、金融时序、医学影像等非平稳信号场景。压缩包共3个m文件分别承担小波变换、交叉小波变换以及综合相关分析功能整体仅9KB代码精简且结构清晰易于阅读和二次修改。目前已有764人学习下载。资源以wtc.m为分析主线完整演示从小波基选取、信号分解、小波系数相关性计算到生成小波图谱与交叉小波图谱的处理流程支持针对不同数据特点调整小波参数便于观察多尺度上的统计关联与相位差异。对于需要快速搭建交叉小波分析框架的研究人员或工程师这套脚本能显著缩短算法实现与调试验证的时间是把小波理论落地到实际数据的实用工具。1. 交叉小波分析到底解决什么问题把两组非线性时序的关联拆到每个尺度做时间序列相关分析久了你一定碰过这种翻车两组数据整体相关系数算出来 0.7画在一起却明显感觉走势并不同步Pearson 相关只给一个数值把高频抖动和低频趋势全混在一起根本说不清“在哪个时间尺度上谁领先谁”。交叉小波分析cross wavelet transform配合小波相干 WTC要解决的正是这个黑匣子问题它把相关分析拆到“时间-尺度”二维平面上横轴是时间纵轴是周期或尺度颜色代表两列信号在该时段该周期上的共同能量和相关强度箭头代表相位超前滞后。水文、气象、金融、地震信号分析都在用这一套做信号关联诊断。这篇笔记给要立刻跑数据的从业者从公式到实现到坑位一次讲明白。2. 先把三个概念拆开交叉谱、交叉小波、小波相干不是一回事很多人把“交叉谱”“交叉小波”“小波相干”混着叫实际它们的用途差别很大。这章先把数学关系理清后面调参数和读图才不踩坑。2.1 交叉谱把相关性变成频率的函数但丢掉局部时间信息传统交叉谱本质是频域里的互相关对两列时序 x(t)、y(t) 做傅里叶变换得到 X(f)、Y(f)交叉谱密度定义为S_xy(f) E[ X(f) · Y*(f) ]它的模叫交叉幅值谱反映两个信号在频率 f 上共同振荡的能量大小辐角就是相位谱反映该频率上 y 相对 x 的滞后。经典的互相关函数和交叉谱只能回答“整个时间窗内哪个频率上有关系、滞后多少”回答不了“这段关系是否只在某一时段成立”。真实数据里这恰恰是最常见的问题比如降雨-径流关系在枯季和汛期完全不同脑电信号在清醒和睡眠不同频段的耦合也完全两样。傅里叶基函数在时间上无限延伸一旦信号里存在突变、断点或非平稳趋势交叉谱就会把这些局部变化平均掉低频和高频混在一张谱图上。结果看着平滑实际信息已经错位。这就是为什么小波分析流行之前做交叉谱的工程师经常对着一条光溜溜的谱线挠头谱峰是真实关系还是多段信号拼出来的假象完全不知道。2.2 交叉小波功率谱加上时间定位却保留频率定位小波变换给每个尺度 a 和时间位置 τ 定义一个基函数。最常用的是复值 Morlet 小波ψ(η) π^(-1/4) e^(i ω₀ η) e^(-η²/2)其中 ω₀ 通常取 6这个值让 Morlet 在时间分辨率与频率分辨率之间保持比较好的平衡ω₀ 越大频率越细但时间定位越粗。对 x(t) 做连续小波变换W_x(s, τ) ∫ x(t) (1/√s) ψ*((t-τ)/s) dt注意 s 是尺度不是频率。对 Morlet 且 ω₀6尺度与等效傅里叶周期近似满足 T ≈ 1.03 s画图时纵轴通常直接用周期而不是尺度本身。分别对两列信号做 CWT 之后交叉小波谱定义为W_xy(s, τ) W_x(s, τ) · W_y*(s, τ)其中 * 是复共轭。|W_xy|² 就是交叉小波功率表示两个序列在某个时间 τ、某个尺度 s 上共同波动的强度arg(W_xy) 则是相位差反映该时频点二者之间的领先滞后关系。傅里叶交叉谱只有一个频率轴交叉小波谱多出一个时间轴这是它与交叉谱最本质的差异。2.3 小波相干 WTC消除信号强弱影响得到逐尺度的“相关系数”交叉小波功率有个麻烦它受各自信号幅值影响。如果 x 的幅值是 y 的十倍交叉谱会被大幅值信号主导小幅值信号里本来很强的相关性被淹没。小波相干就是在交叉谱基础上用两列信号各自的功率做归一化R²(s, τ) |S(W_xy(s, τ))|² / [ S(|W_x(s, τ)|²) · S(|W_y(s, τ)|²) ]这里的 S 是平滑算子一般做法是尺度方向做高斯平滑、时间方向做移动平均或高斯平滑。R² 的值域是 0 到 1解释起来和相关系数类似但它是每一个时间-尺度格点上的“局部相关系数”。也正因为是局部估计小样本区域很容易出现虚高所以不能只看颜色深浅必须搭配显著性检验看。实际操作中交叉小波功率和 WTC 要分开解读交叉小波功率突出“共同能量大的区域”WTC 突出“关系稳定且显著的区域”。有的报告把两张图叠在一起讲信息重复读者看着费劲。我一般建议诊断阶段两张都画汇报阶段只保留 WTC 加显著边界线。2.4 显著性检验白噪声还是红噪声直接决定显著区面积小波相干值天然处处非零哪怕两组独立随机序列也能算出一堆 R²。判断哪些区域真正有意义通常用蒙特卡洛方法对 x、y 分别拟合 AR(1) 红噪声模型生成几百组代理序列对计算代理序列的 WTC 分布再把真实 WTC 值与 95% 或 99% 分位数比较。用白噪声代替 AR(1) 是新手最常见的错误。自然时间序列普遍存在惯性低频段能量天然比高频高白噪声零假设给不出这种谱形会把低频区域大量普通相关误判为显著。后面第 4 章会专门讲这个坑。3. 用 Python 复现 wtc-r16 流程从数据准备到画出时频相关图这章给出一套最小可跑的流程。我用 Python 加 PyWavelets 和 SciPy 实现代码结构按“预处理 → 小波变换 → 交叉小波与相干 → 绘图”四步展开。3.1 数据预处理先做差分、去趋势再谈相关性交叉小波对非平稳分量非常敏感。如果两组数据都带线性趋势低频区域会产生一大片虚假的相关条带看起来高度相关其实只是同向漂移。我先对序列做两步检查一是画时序图粗看趋势和突变二是做 ADF 平稳性检验。ADF 检验 p 值大于 0.05 时通常先做一阶差分或线性去趋势。差分会让低频信息丢失去趋势则会保留低频波动但去掉长期漂移怎么选取决于分析目的。研究年际到年代际关系时保留低频更重要我倾向用局域回归去趋势而不是差分研究高频耦合关系时可以直接差分。下面代码以去均值和线性去趋势为例。import numpy as np from scipy.signal import detrend # 构造两列合成信号便于验证代码正确性 dt 0.02 # 采样间隔单位秒 N 1024 t np.arange(N) * dt # x5 Hz 和 15 Hz 两个主频成分持续整段 x (0.8 * np.sin(2 * np.pi * 5 * t) 0.6 * np.sin(2 * np.pi * 15 * t) 0.05 * np.random.randn(N)) # y5 Hz 与 x 同频但滞后 0.3 秒15 Hz 在 t 8 秒后失耦 y (0.7 * np.sin(2 * np.pi * 5 * (t - 0.3)) np.where(t 8, 1, 0) * 0.6 * np.sin(2 * np.pi * 15 * t) 0.05 * np.random.randn(N)) # 去均值和线性趋势让序列近似平稳 x detrend(x) - np.mean(x) y detrend(y) - np.mean(y)这段代码先构造了已知答案的合成数据5 Hz 分量有固定 0.3 秒滞后15 Hz 分量在 8 秒后消失。后面算出来的图若能明显看到这两段结构说明流程可信。detrend默认去除线性成分np.mean再把均值归零对强非平稳序列建议先做 ADF 检验或一阶差分不要跳过这一步。3.2 Morlet 小波变换与 r16 参数每个倍频程 16 个子尺度“wtc-r16”里的 r16常见做法指每个倍频程细分成 16 个子尺度也就是尺度间隔参数 dj 1/16。J 是总尺度步数一般由数据长度和最小尺度决定目标是最大尺度能覆盖到约序列长度的一半。import pywt def wtc_r16_scales(dt, N, dj1/16): 按 r16 方案生成尺度序列 s0: 最小尺度取 2 倍采样间隔 J: 总尺度步数按最大尺度约等于 N*dt/2 反推 s0 2 * dt J int(np.floor((np.log2(N * dt / s0)) / dj)) scales s0 * np.power(2.0, np.arange(0, J 1) * dj) return scales scales wtc_r16_scales(dt, N) # 对 x、y 做连续小波变换morl 即 Morlet 小波 Wx, freqs pywt.cwt(x, scales, waveletmorl, sampling_perioddt) Wy, _ pywt.cwt(y, scales, waveletmorl, sampling_perioddt)pywt 的cwt接受自定义尺度序列和sampling_period做卷积时已经处理了尺度与信号的采样关系。返回的freqs单位是采样间隔的倒数实际频率要用freqs / dt换算画周期图时也直接用1 / freqs对应周期。dj1/16是 r16 的关键参数它控制尺度轴的分辨率想更细可以改成 1/32但计算量成倍上升图上的显著区域也不会有本质变化。3.3 交叉小波功率、相干、相位一次性算出有了Wx和Wy后面几行就能算完核心量。交叉小波谱是 Wx 乘 Wy 的共轭交叉功率是它的模方相位是复角。小波相干要额外做平滑平滑宽度直接影响结果这里单独写一个函数说明。from scipy.ndimage import gaussian_filter1d # 交叉小波谱 Wxy Wx * np.conj(Wy) power np.abs(Wxy) ** 2 # 交叉小波功率 phase np.angle(Wxy) # 相位差值域 [-pi, pi] def smooth_by_scale(a, scales, dt, time_win0.5, scale_win0.5, dj1/16): 时间维按各尺度等效周期的固定比例做高斯平滑 尺度维按倍频程固定宽度做高斯平滑 sm np.empty_like(a) for i in range(len(scales)): # Morlet(w06) 的尺度到傅里叶周期近似为 1.03 倍尺度单位先转成采样点 period_pts 1.03 * scales[i] / dt sm[i] gaussian_filter1d(a[i], sigmatime_win * period_pts, modeconstant) # 尺度方向再平滑sigma 用倍频程宽度除以 dj 换算成“格点数” sm gaussian_filter1d(sm, sigmascale_win / dj, axis0, modeconstant) return sm # 分别平滑三个量再算相干 sm_xy smooth_by_scale(Wxy, scales, dt) sm_xx smooth_by_scale(np.abs(Wx) ** 2, scales, dt) sm_yy smooth_by_scale(np.abs(Wy) ** 2, scales, dt) coherence np.abs(sm_xy) ** 2 / (sm_xx * sm_yy 1e-10)平滑宽度是 WTC 结果“敏感度”的来源。time_win0.5表示时间平滑窗等于半个等效周期scale_win0.5表示尺度方向平滑半个倍频程。time_win 调大相干图更平滑但时间分辨率变差调小则噪点增多。两个参数没有绝对最优值我习惯先用 0.5 跑一遍再对比 0.7 的结果确认主要结构没有崩塌。import matplotlib.pyplot as plt period 1.0 / freqs fig, ax plt.subplots(figsize(10, 5)) cf ax.contourf(t, period, coherence, levelsnp.linspace(0, 1, 31), cmapjet, extendboth) ax.set_yscale(log) ax.set_ylim(period[-1], period[0]) ax.set_ylabel(period (s)) ax.set_xlabel(time (s)) cb fig.colorbar(cf, axax) cb.set_label(WTC coherence)这段图配上 log 纵轴便于同时观察高频和低频的周期结构。合成数据里5 Hz 对应周期 0.2 秒图上应在 0.2 秒附近出现一条贯穿全时段的红色高相干带15 Hz 对应周期约 0.067 秒图上应在 8 秒前有显著条带8 秒后消失。看到这两段结构说明参数设置正确可以把这套流程换到自己的数据上。4. 交叉小波排坑记录COI、趋势项、显著性假阳性和相位平均这章写给正在跑数据的人。以下四个问题是我在实际项目里反复遇到、也帮同事排查过多次的典型坑每条按“现象 → 原因 → 解决”拆开讲。4.1 边界陷阱图两侧的“高相干”大多是假的现象WTC 图左右两侧靠近数据起止位置出现竖条状高相干区域看着像是数据在端点附近突然变得高度相关但业务上完全解释不通。原因小波变换在信号两端会截断小波窗口超出数据边界部分被强制补零。补零导致边界处小波功率被人为压低而交叉小波功率是两个信号功率的乘积两边都在边界衰减比值相干反而可能虚高。这种效应集中在“锥形影响区”COI 内COI 之外的数据基本不可信。解决绘图时务必画出 COI 边界并习惯性地忽略边界外区域。对 Morlet 小波COI 的时间范围可以用 e-folding 时间近似取 τ √2 · s。实际操作时我用类似下面几行做标记# 计算 COI 边界每个尺度对应的时间影响范围 time_axis t coi 1.03 * np.sqrt(2) * scales[:, None] / dt # 影响半径单位采样点 # 填充边界外的区域让它视觉上与主图区分开 ax.fill_between(time_axis, period[-1], period, wherenp.arange(N) coi.min(), colorwhite, alpha0.7) ax.fill_between(time_axis, period[-1], period, wherenp.arange(N) N - coi.min(), colorwhite, alpha0.7)注意 COI 在不同尺度宽度不同高频短周期COI 窄低频长周期COI 宽。填色时如果只画一条固定竖线会漏掉低频段的边界问题。正确做法是让白色遮罩在低频段向中间扩展。4.2 低频大片显著区趋势项没处理干净的典型症状现象两组完全不相关的序列WTC 低频区域出现大面积 R² 大于 0.8且显著区从图底一直蔓延到图中间颜色像一条“地毯”。原因两条序列都包含低频趋势或缓慢漂移。小波相干在小尺度高频能区分噪声但在大尺度低频上趋势成分占据了主导能量平滑后局部相关天然接近 1。这不是“发现了真实关系”而是共趋势造成的虚假相关。解决在进入交叉小波分析前检查每个序列的自回归系数和 ADF 检验结果。带趋势就做线性去趋势或一阶差分带季节周期就做季节差分或用 STL 分解去掉周期成分。处理后再跑一遍如果低频显著区明显缩小说明之前看到的确实是趋势假象。还有一个笨办法但很有效把 x 随机打乱顺序和 y 重算 WTC如果低频仍显著那基本可以断定是趋势或边界效应所致。4.3 显著性检验用错谱模型白噪声零假设给的阈值太低现象用某开源包默认参数跑 WTC显著区面积大得夸张几乎整个图都是 95% 显著。拿同一数据换用另一个工具显著区却缩水一大半。原因WTC 显著性检验的核心是零假设。若零假设是白噪声其谱是平的低频阈值很低自然序列低频本来就有能量轻松越过阈值。若零假设是 AR(1) 红噪声其谱低频能量更高阈值上移普通相关就不再显著。自然时间序列绝大多数符合红噪声特性所以必须用 AR(1) 代理检验。解决常见做法是对每组序列分别拟合 AR(1) 模型生成 200 到 1000 组代理序列对计算代理 WTC 的 95% 分位数绘制显著性等高线from statsmodels.tsa.ar_model import AutoReg # 对 x 拟合 AR(1)生成一个代理序列的简单示例 def ar1_surrogate(x): order 1 model AutoReg(x, lagsorder).fit() ar_coef model.params[1] resid model.resid innov np.random.choice(resid, sizelen(x), replaceTrue) surrogate np.zeros_like(x) surrogate[0] x[0] for i in range(1, len(x)): surrogate[i] ar_coef * surrogate[i-1] innov[i] return surrogate实际显著性检验需要生成大量代理对逐个算 WTC再取每个格点的 95% 分位数。这个流程计算量不小200 组代理在 1024 点数据上大约几十秒到几分钟可接受。代理数量太少时显著性边界抖动很大我一般至少跑 500 组。4.4 相位箭头乱转对相位求算术平均等于找死现象WTC 图上相位箭头方向在相邻格点间跳动有时从 179° 到 -179°画出来的箭头“满天乱飞”没法解释谁领先谁。原因相位是圆周量存在 ±π 周期性。一个格点为 3.0 rad另一个格点为 -3.0 rad算术平均算出来是 0实际两者只差 0.28 rad应该接近 3.0。直接求平均会把正确的滞后关系抹掉产生很多伪“相位突变”。解决相位计算全程使用圆统计。先算每个格点的相位再沿时间轴做圆均值def circular_mean(angles): 一组角度的圆平均返回平均值和一致性长度 r r_vec np.mean(np.exp(1j * np.array(angles))) return np.angle(r_vec), np.abs(r_vec)r 越接近 1说明该区域相位越集中、领先滞后关系越稳定r 接近 0说明相位散乱即使 WTC 显著也不该解释成稳定滞后。我在报告里一般把 r 大于 0.5 的区域才标上箭头小于 0.5 的置灰。这个习惯避免了大量“看着像有滞后、实际说不清方向”的结论。5. 用合成测试验证完整流程再上真实数据最后说一个我的固定习惯任何交叉小波分析先不做真实数据而是先跑合成测试。选定一个已知的时变滞后结构比如 5 Hz 滞后 0.3 秒、15 Hz 在 8 秒后失耦按第 3 章的代码生成数据并跑通全流程。检验点有三个一是 WTC 高频带是否出现在预期周期位置二是滞后换算是否等于 0.3 秒三是 8 秒后失耦频段是否在真实数据中“消失”。三个点都对再换真实数据。真实数据上我还会额外输出一张“区域平均相位差表”。做法是选定一个显著区域用圆平均算出该时段该尺度的平均相位 φ换算成时间滞后 Δt φ · T / (2π)其中 T 是该尺度对应的周期。表格给出区域起止时间、周期范围、圆平均相位、时间滞后和一致性长度 r这样业务人员不需要看复杂的时频图也能直接读到“第 3 到第 5 小时周期 8 到 16 分钟的波动y 比 x 平均领先约 1.2 分钟一致性 0.8”这类结论。输出这类定量结论时我始终保留两条底线一是所有显著结论必须带 AR(1) 蒙特卡洛显著性边界二是所有被解释的箭头必须满足圆一致性 r 大于 0.5。不满足这两条的区域一律不下因果判断。交叉小波分析最怕的就是把相关图讲成因果故事一个边界假象或趋势残留就可能让整个结论翻车。希望这套流程和踩坑笔记能帮你少走这圈弯路。本文还有配套的精品资源点击获取