简介：面向气象与气候研究者的 MATLAB 小波功率谱分析工具包，专注于降水、气温、太阳黑子等时间序列的周期性检测。压缩包共 28 个文件，以 m 脚本为主，辅以 png 结果图、dat/txt 数据文件及 bak/asv 备份文件，总大小仅 184KB，轻量且便于按需调用。目前已有 2149 人学习下载。资料覆盖连续小波变换、显著性检验、功率谱计算与可视化等完整流程，包含 Morlet 小波等多种实现及配套示例数据，如气温序列、太阳黑子数、Nino 海温等。使用者可基于现有脚本直接修改参数，快速绘制小波功率谱图并提取显著周期，从而识别气候序列中的局地周期变化和非平稳特征。对于正在开展周期分析、需要可复现脚本作为起点的高校学生和科研人员，这份工具包提供了清晰的代码框架与真实数据演示，能够显著降低从方法理解到实际出图的门槛，也可作为进一步开展极端天气与气候变率研究的参考基础。

1. 小波功率谱：同一份降水数据，傅里叶只能“听完”，小波却能“拆开看”

一张某地几十年的逐月降水距平图，丢给傅里叶变换，你得到的是一个静态的周期列表：若干个峰值频率，每个峰值一个全局周期强度。可一旦把同一份数据画成小波功率谱，你会看到完全不同的信息：横轴是时间，纵轴是周期，颜色深浅代表振荡能量，某个周期可能在 1970—1985 年特别强，之后突然消失，又在 2000 年前后重新出现。这种“周期随时间变化”的现象，傅里叶谱根本表达不出来，因为它的基函数是无限延伸的正弦波，天然把时间信息全部平均掉了。小波功率谱之所以能成为气候周期分析里的主力工具，靠的是它同时保留“周期”与“时间位置”两个维度。这篇文章要拆的，正是这套机制的核心：连续小波变换为什么能同时定位时间和频率、功率谱的显著性检验到底在检验什么、以及 MATLAB 里那些参数和阴影区域背后的实现逻辑，把周期分析的黑匣子打开。

2. 连续小波变换：它如何做到“既看得到周期，又记得住时间”

2.1 内积视角：小波变换的本质是把信号投影到“局部振荡”上

连续小波变换的定义式并不复杂：

W(a,b) = (1/√a) ∫ f(t) ψ*((t - b)/a) dt

这里 ψ(t) 是小波母函数，a 是尺度，b 是平移量。和傅里叶变换的内积结构相比，关键差异全在这个 ψ((t-b)/a) 上：傅里叶基函数 e^(-iωt) 在整个时间轴上无限振荡，没有“位置”概念；小波则把自身压缩或拉伸成一个局部振荡团块，再用 b 把这个团块搬到时间轴的任意位置。a 大于 1 时窗口被拉宽，匹配低频振荡；a 小于 1 时窗口收缩，匹配高频细节。前面的 1/√a 是能量归一化因子，保证不同尺度下的小波具有相同能量，否则大尺度窗口天然占比大，功率谱在不同周期之间就没法比较。

这个“局部振荡”的特性，是小波能同时给出时间和频率的根本原因。如果信号里存在一个只在 1980—1995 年显著的 5 年周期，傅里叶变换把这个局部振荡摊到整个时间轴上，得到的峰值会被严重稀释；而小波只会让那一时段对应的尺度系数显著变大。你看到的时间-周期热力图，本质上是信号在每一个“局部片段”上与不同尺度小波的内积强度。

2.2 Morlet 小波：为什么周期分析几乎都选复数小波

气候周期分析里最常见的母函数是 Morlet 小波：

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

它是一条复正弦被高斯包络截断。选择复数小波有两个实际好处：一是可以同时提取振荡的瞬时幅度和相位；二是复数小波的功率谱在统计检验时天然对应两个自由度，可以和卡方分布对接，这点在第三章会展开。

ω₀ 是中心频率参数，默认值通常是 6。这个数字意味着小波窗口内大约能容纳三个完整振荡周期，窗口形状在时间定位和频率纯度的权衡中取了一个工程上的平衡点。ω₀ 越大，振荡次数越多，频率分辨率越好，但高斯包络衰减越慢，对周期随时间的变化越不敏感；ω₀ 越小，时间定位越准，但频谱泄露严重，周期峰值变得模糊。不确定性原理在这里是硬约束：Δt·Δω ≥ 1/2，你不可能同时把时间和频率分辨率无限提升。Morlet 只是把这条物理界限在“认出周期”和“定位时刻”之间做了实用的取舍。

2.3 尺度换周期：一个被 MATLAB“藏起来”的关键换算

直接调用 MATLAB 新版 cwt 时，尺度到周期的换算几乎是无感的，因为函数会直接返回周期数组：

% x 为逐月降水距平序列，dt = 1 表示采样间隔为 1 个月 % 'amor' 是 MATLAB 对复数 Morlet 小波的命名 [wt, period] = cwt(x, 'amor', dt);

但理解换算公式仍然必要。Morlet 小波的尺度 a 与周期 T 的关系为：

T = 4πa / (ω₀ + √(2 + ω₀²))

当 ω₀ = 6 时，系数约等于 1.03，也就是说用尺度 a 近似当作周期 T，误差只有 3% 左右。很多老脚本里直接把尺度当周期用，在这个参数下结论不会颠覆。但如果换了小波基，情况完全不同：墨西哥帽小波的换算系数接近 4，直接套用 Morlet 的近似值会得出完全错误的周期。这也是为什么在参数说明里永远要把“用的哪种小波、中心频率是多少、换算因子是多少”写清楚。

2.4 频域加速：把卷积变成逐尺度点乘

直接按定义式计算连续小波变换，每个尺度需要做一次与信号长度同阶的卷积，计算量很大。实际实现几乎都走频域路线：对信号做一次 FFT；对每个尺度 a，在频率轴上构造小波 ψ 的频域形式；二者点乘后做 IFFT。核心只有三行：

xhat = fft(x); % 信号变换到频域 f = (0:n/2) / (n * dt); % 单边频率轴 w(si, :) = ifft(xhat .* conj(psi_hat)); % 频域相乘后回到时域

整个变换的复杂度降到 O(S · N log N)，其中 S 是尺度个数。这也解释了为什么小波功率谱可以逐行在普通笔记本上算完几十年的月数据。理解这个流程之后，再看 MATLAB 的 cwt 输出就不会觉得它“玄学”，它只是在这个频域骨架上补充了归一化因子、边界处理和数值稳定性优化。

3. 功率谱与显著性检验：为什么小波图上的“红”不能直接当结论

3.1 小波功率谱与全局谱：能量如何在时频平面分布

小波系数 W(a,b) 是复数，取模平方 |W(a,b)|² 就得到小波功率谱。它的物理含义可以理解为：在某个时间位置 b 上，尺度 a 对应的振荡能量密度。逐点看这张图，颜色越深代表该时刻该周期的振荡越强，但它只是“功率大小”，不是“统计结论”。

把功率谱沿时间轴平均，得到的是全局小波功率谱。它和傅里叶功率谱的作用相似，但数值稳定性更好：傅里叶谱把每个频率贡献成一整条竖线，峰形尖锐；小波谱在频率方向上有高斯包络的展宽，峰形平滑。实际分析时，全局谱适合用来初筛候选周期，时频图用来确认周期存在的时段和连续性。

3.2 红噪声假设：气候序列的背景谱不该用白噪声

对功率谱做显著性检验，第一步是建立“如果没有真实周期，谱长什么样”的零模型。很多入门分析直接用白噪声背景，这是最常见的系统性错误。降水、气温这类气候序列存在很强的惯性，今天的值会影响明天的值，普遍可以用一阶自回归模型 AR(1) 描述：

x(t) = α·x(t-1) + ε(t)

其中 α 是滞后一阶自相关系数，用一行 MATLAB 就能估计：

alpha = corr(x(1:end-1), x(2:end));

AR(1) 过程的功率谱不是平的，而是低频强、高频弱的连续谱，这就是“红噪声”。α 越接近 1，低频能量占比越高。如果拿白噪声做背景检验，低频端的红噪声能量会被误判成显著周期，导致图上出现一片根本没有意义的“显著区”。

AR(1) 红噪声的理论功率谱写法是：

P(k) = σ² (1 - α²) / (1 + α² - 2α·cos(2πk/N))

其中 k 对应第 k 个频率分量，σ² 是序列方差。这个公式是显著性检验的背景谱基准。

3.3 卡方检验与置信阈值：从“颜色深浅”到“显著”

Morlet 是复数小波，系数实部和虚部分别近似服从正态分布，功率谱值 |W|²/σ² 因此近似服从自由度为 2 的卡方分布。所以 95% 置信水平的功率阈值是：

threshold = 红噪声背景谱 × χ²₂(0.95) / 2

MATLAB 里直接算临界乘子：

p = 0.95; th = chi2inv(p, 2) / 2; % 约等于 2.996

自由度为 2 的卡方分布在 95% 分位数约 5.991，除以 2 得约 2.996。也就是说，某点的功率只有超过红噪声背景谱的约 3 倍，才能说它在 95% 置信水平下“不等于红噪声”。99% 置信水平对应乘子约 4.605，要求更高。

这套检验依赖高斯假设，如果序列有明显的季节循环、极端值或趋势，AR(1) 参数会被污染，检验结论就不可信。这也是为什么预处理和显著性检验必须绑定出现，而不是分开的两个步骤。

3.4 显著周期判读：连续带状区域才算数

显著性检验解决的是“单点是否显著”，但画出来的图里偶尔会有孤立的显著点穿过阈值线。这些点往往来自单一极端值或边界失真，不能解释为周期信号。判读时要坚持两条规则：第一，显著区域必须是连续的带状，沿时间轴跨越至少两到三个完整振荡周期；第二，显著带出现的时间区间本身要有规律，比如持续十年以上的 4—6 年周期带，才值得被写进结论。

4. MATLAB 实操：从 cwt 到自建算法，关键参数背后的实现逻辑

4.1 新版 cwt 的用法与返回值语义

新版 MATLAB 的 cwt 接口把采样间隔和复 Morlet 小波封装成了两个参数：

% dt 单位是月，返回的 period 单位与 dt 一致 [wt, period, coi] = cwt(x, 'amor', dt);

第三个输出 coi（cone of influence）是边界影响锥，绘制时 MATLAB 会自动把它渲染成图上的阴影区。需要注意：老版本接口 cwt(x, scales, 'morl') 返回的是尺度序列而不是周期序列，混用新旧接口是脚本报错和数据解读错误的常见来源。绘图时用对数周期轴通常是更好的选择，因为低频端相邻尺度在小波频域上非常接近，线性轴会把低频显著区挤成一团。

4.2 频域自建 Morlet CWT 的骨架

对于想确认“小波变换内部到底做了什么”的开发者，自建一个教学级实现比文档更有说服力。尺度序列的对数等距设计与 COI 上限是其中的关键：

function [w, period] = morlet_cwt_demo(x, dt) n = numel(x); x = x(:); % 最小尺度取 2 倍采样间隔，小于它窗口内装不下一个振荡 s0 = 2 * dt; % 最大尺度取序列长度的 1/4，更长时边界效应主导整段数据 smax = floor(n / 4) * dt; % 对数等距尺度序列：每个倍频程 16 个尺度点 scales = 2 .^ (log2(s0):1/16:log2(smax)); omega0 = 6; f = (0:n/2) / (n * dt); % 单边频率轴 xhat = fft(x); w = zeros(numel(scales), n); for si = 1:numel(scales) a = scales(si); % Morlet 小波频域形式：只在中心频率附近有能量 psi_hat = exp(-0.5 * (a * 2 * pi * f - omega0).^2); w(si, :) = ifft(xhat .* conj(psi_hat)); end % omega0=6 时尺度到周期的近似换算，误差约 3% period = scales * 1.03; end

这段代码说明了两件事：第一，频域实现中每个尺度独立做一次点乘和 IFFT，累加顺序不影响结果，便于并行；第二，ψ_hat 的形状决定了小波等效于一个带宽随尺度变化的带通滤波器，尺度 a 越大，通带中心频率越低、带宽越窄，这正是低频周期分辨率的来源。真实研究建议仍用 MATLAB 自带的 cwt，它在归一化因子、零频修正和边界处理上更严格，这段骨架用于理解机制和调试参数。

4.3 COI 边界锥：藏在图上但决定可信度的阴影

COI 的本质是小波能量在边界处发生泄漏的范围。Morlet 小波的高斯包络衰减到 e^(-1) 的时间约为 √2a，因此边界效应影响的时间宽度随尺度增大而增大。a = 12 个月时，边界影响区间约为 17 个月；a = 60 个月时，边界影响膨胀到约 7 年。这意味着低频周期即使在图上有大范围深色区域，只要它接近数据两端，就有可能是信号截断造成的伪影。判读时不看阴影区内的功率，是周期分析的底线操作。

4.4 预处理三件套：季节性、趋势、方差

降水气温序列直接用原始值做小波分析，低频端几乎必然出现虚假强信号。季节循环本身就是一个确定性周期，它会把 AR(1) 的 α 抬高，让红噪声背景谱失真。标准流程是：先做逐月距平，减去每个月的历史平均值，消除季节循环；再对距平序列做标准化，让不同季节方差可比；最后检查是否有长期线性趋势，有则去趋势。这里给出一个常见方案对照：

5. 周期分析常见问题：五条踩坑记录与排查思路

5.1 边界处的伪高功率

现象：小波功率谱图在数据两端出现巨大的红色三角区，低频尤其明显，甚至掩盖了内部真实信号。

原因：信号在截断处的不连续产生了高频泄漏，而低频小波窗口很长，边界附近的系数有一半窗口悬在数据外，被零填充强行压平，制造出虚假的低频强能量。

解决：绘图时始终显示 COI 阴影区，并明确声明结论只讨论 COI 范围内的谱；如果站点序列较短，优先降低最大尺度或缩短结论周期范围。看到两端红色区域的第一反应不应该是“好显著的周期”，而是“这是边界效应”。

5.2 小波谱与傅里叶谱的峰值周期对不上

现象：同一份数据，傅里叶谱主峰在 11 年，小波全局谱主峰却落在 13 年附近。

原因：尺度到周期的换算不一致。新版 cwt 返回周期数组但内部用滤波组实现，老代码用尺度直当周期，或者中心频率 ω₀ 不是 6，都会让峰值偏移。

解决：确认脚本用的是哪一代接口；统一用 cwt 返回的 period 数组作为横轴；手动算尺度时用 T = 4πa / (ω₀ + √(2+ω₀²))，不要拍脑袋近似。偏差在 10% 以内时，先怀疑尺度换算而不是数据。

5.3 显著性检验后整张图全不显著

现象：功率谱颜色很饱满，但套用红噪声检验后几乎没有区域越过 95% 阈值。

原因：预处理偷工减料。季节循环和趋势残差让 AR(1) 的 α 被高估，红噪声背景谱被抬高，真实周期信号的能量相对背景就不够突出了。

解决：把预处理做到位——逐月距平、去趋势、标准化，再重新估计 α。还有一个高频翻车点：α 的估计用了全序列而序列本身含缺失值，缺失插值不当会把自相关拉高。数据质量差时用插补后的序列做检验，不如直接在检验时报告缺失时段。

5.4 同一组数据，两张谱图看起来完全不同

现象：两个脚本跑同一份数据，一个周期带清晰，另一个模糊成一片。

原因：Morlet 小波的 ω₀ 不同或尺度采样密度不同。ω₀ 从 6 调到 3，频率分辨率明显下降，周期带变宽；尺度网格太疏则低频端漏掉峰值。

解决：固定参数并在报告里写明“复 Morlet，ω₀=6，每倍频程 16 尺度”。比较两张谱时先归一化尺度轴，再看峰值位置是否一致，而不是盯着色带粗细。小波谱本来就是“加上参数假设的谱估计”，参数不同结果不同不是 bug，是设计。

5.5 把归一化功率谱拿去和文献数字直接对比

现象：某开发者想把本站点的小波功率谱数值与论文里的图做定量对比，发现量级差了几十倍。

原因：MATLAB 新版 cwt 默认输出的是归一化功率谱，单位是“相对于总方差的能量比例”，老版本或其他语言实现可能是原始能量谱，两者差一个与数据方差和尺度网格有关的常数因子。

解决：对比时只比较“显著周期的位置和时段”，比绝对值没有意义。如果一定要比较，先在同一标准化流程下重建功率谱并做全局谱归一化，再看形状。

6. 从显著周期到可解释信号：小波基取舍与周期信号重建

6.1 参数取舍：时间分辨率还是频率分辨率

小波功率谱图最终呈现的效果，很大程度上在选小波基时就已经定型。ω₀ 取大，周期定位更锐利，但时间轴上周期迁移的细节被模糊；ω₀ 取小，时变特征更清楚，代价是不同周期的功率相互涂抹。没有一套参数能同时做到“周期精确到 0.1 年”和“时间边界精确到 1 个月”。更可取的做法是把小波分析当作假设生成器：先用 ω₀=6 的复 Morlet 找到候选显著周期带，再用更高 ω₀ 的小波确认周期数值，必要时变换到显著周期对应的窄带做带通滤波，验证这个周期不是小波基函数“造”出来的。

6.2 从谱图到波形：把显著周期带还原成可分析的序列

功率谱证明了“某时段存在显著周期”之后，还差一步：把这段周期信号本身提取出来，才能讨论它的振幅变化和相位迁移。一个实用做法是针对显著周期带做带通滤波，再将滤波结果和小波显著带对比。二者时间上的包络一致，才能确认提取的不是数学伪影。整个流程走完，结论才从“图上有点红”变成“该站 1975 年后存在显著的 4—6 年周期振荡，振幅在 1980 年前后达到峰值”。

多年前我在某站降水数据上第一次做出漂亮的小波功率谱图时，直接按全局谱峰值写了“主周期约 28 个月”的结论，后来复查发现主峰几乎全部落在 COI 阴影区内，补做红噪声显著性检验后，显著区域只剩短短一截。从那以后，我每次做周期分析都强制走完整流程：逐月距平、去趋势、AR(1) 估计、功率谱显著性检验、COI 内判读、滤波重建验证，全部通过才写结论。小波功率谱的价值在于给你一个诚实看待周期的方式，但也只在流程完整时才可靠。希望这套拆解能帮你在下一次打开 cwt 时，心里多一层对参数、边界和统计阈值的把握。

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