简介传统傅里叶变换难以刻画随时间变化的周期特征小波功率谱则能在时间与频率两个维度同时揭示信号变化。面向气象与气候研究者的MATLAB小波功率谱分析资源包专门用于降水、气温、太阳活动等时间序列的周期识别与非平稳信号分析。压缩包共28个文件以15个m脚本为核心另有png功率谱图、txt样本数据、dat数据文件、bak备份文件等整体仅184KB结构紧凑易检索。脚本覆盖数据读取、连续小波变换、功率谱计算、显著性检验和绘图全流程并自带小波基函数、显著性检验等常用子程序可在常见MATLAB环境中直接运行。包内提供太阳黑子数、Nino3区海温、气温等多个真实示例并配有标注95%置信水平的周期谱图便于对照验证方法效果。目前已有2149人下载学习通过研读源码和示例可掌握Morlet小波基选择、尺度参数设置、峰值提取与周期解释等关键操作进而迁移到自身研究数据中适合具备MATLAB基础的气候与环境科学研究者。1. 小波功率谱为什么 FFT 答不出周期从哪年开始一句话说FFT 把气象序列里的周期性「平均」掉了。两条观测序列一条前 30 年有显著 8 年周期、后 20 年变成随机波动另一条整个 50 年都有稳定 8 年周期它们的 FFT 幅值谱几乎长得一样只差在峰高上。可气象诊断里最常问的是「这个周期从哪一年开始、到哪一年结束」——降水、气温的驱动因子本身在演变周期不会老老实实站满整个序列。小波功率谱的厉害之处就在保留时间轴它把周期演化画成一张随时间铺开的谱图。MATLAB 里做这件事只需要几行代码但真正决定结果对错的是背后的尺度换算、边界效应和显著性检验这些恰恰是最容易被当成「默认设置」翻车的环节。2. 从傅里叶到连续小波功率谱的时间轴是怎么多出来的2.1 傅里叶变换的全局平均周期被「加」没了常规周期分析从 FFT 开始。对一个长度为 N 的序列 x(t)傅里叶系数写成X(f) Σ x(t) e^{-i2πft}求和符号把所有时段的样本叠加在一起输出的 X(f) 是这个频率分量在整个时间跨度上的「平均强度」。如果周期只在序列的前半段存在后半段是纯随机这个峰依然在只是变矮了。你从频谱图上能读出「有一个 8 年左右的准周期」却读不出它在什么年代活跃——时间索引 t 在求和过程中被消掉了。一个直观的例子构造 50 年的月降水序列前 30 年由「8 年周期 噪声」构成后 20 年只有噪声。FFT 幅值谱里 8 年峰依然清晰峰高大约只有全程稳定周期序列的六七成。你看这幅图无法判断峰是「全程均匀贡献」还是「前段集中爆发」。而小波功率谱图上会出现一条清晰的断裂带功率在 30 年位置戛然而止。这就是时频分析存在的理由。真实气候序列里这种「周期中断」非常常见。某个年代大气环流型调整原本稳定的周期模式会出现频率漂移、强度波动甚至暂时消失。FFT 做全场平均等于把一个随时间演化的过程压成一个静态统计值丢失的正是最关键的时间演化信息。把序列切段做滑动 FFT 能不能补救这就是短时傅里叶变换STFT的思路开一个固定长度的窗在时间轴上滑动每个窗内做一次傅里叶变换。问题在于窗口长度一旦固定时间分辨率和频率分辨率就互相锁死。窗短一个窗内装不下几个周期的振动频率刻度粗糙低频长周期根本分不出来窗长频率分辨率上去了窗内的时间细节又被平均掉周期何时变化依旧看不出来。气象里要分析的周期跨度往往从一个季度到几十年跨度几十倍一个固定窗不可能两头兼顾。2.2 CWT 的缩放平移机制尺度就是自适应窗宽连续小波变换的定义是W(a, b) (1/√a) ∫ x(t) ψ*((t-b)/a) dt其中 ψ 是母小波a 是尺度因子b 是平移因子。a 大小波基被拉宽对应低频、长周期a 小基函数被压缩对应高频、短周期。这个「窗宽随尺度自动伸缩」的性质正是 STFT 做不到的分析 8 年周期时小波基自动覆盖大约 8 年长度的有效窗口分析 3 年周期时窗口自动缩到 3 年。时间定位能力和频率分辨能力在每一个尺度上独立权衡而不是用同一个窗硬套所有频段。从计算角度看CWT 的每一个系数 W(a,b) 都是信号 x(t) 与尺度伸缩后的小波基做内积。内积大说明信号在这一小段上的局部形态「很像」这个小波基。由于小波基是衰减振荡型系数幅值就代表这个频段在这个时段的振幅模平方就是功率。MATLAB 工具箱底层的实现方式是滤波器组把母小波按尺度伸缩成一组带通滤波器中心频率随 1/a 变化带宽也随 1/a 变化品质因数近似恒定。这就是小波分析对非平稳信号比 STFT 更自然的原因它把序列里所有可能尺度的局部波动都做了一次匹配滤波。需要强调尺度 a 不等于周期。尺度 a 的量纲是时间数值等于母小波基在时域被拉伸的倍数它对应的傅里叶等效周期由母小波的频谱形状决定不同母函数换算不同这个细节放到第 3 章拆。2.3 功率归一化为什么 |W|² 不能拿来直接用小波功率谱的标准定义是 P(a,b)|W(a,b)|²表示时间 b 附近、尺度 a 所在频带内的信号能量强度。真正的坑在于如果直接把 |W|² 拿来画图低频端大尺度的功率会被系统性夸大。原因在 CWT 定义里的 1/√a 归一化因子。它保证不同尺度下小波基的 L2 范数一致却并没有保证功率谱的能量密度可比。直观理解大尺度基函数在时域铺得更宽对信号覆盖的时间范围更长叠加的样本更多模平方数值天然偏大。直接画图的结果就是低频端一片通红2 年周期的峰和 10 年周期的峰看似量级悬殊其实可能是归一化假象。解决方法是把功率谱除以尺度或等价地乘以频率power abs(coefs).^2 ./ scales;这相当于把「每尺度功率」换算成「每单位周期带宽的功率密度」。怎么验证这个归一化是必要的取一段白噪声序列跑完整条链路分别画 abs(coefs).^2 和它除以尺度后的谱前者在低频端明显上翘后者在整条尺度轴上才近似一条水平线。白噪声在各频带应有恒定的功率谱密度所以后者才对。到这一步你得到的是「归一化功率谱」但它的绝对值仍依赖序列方差。要想判断某个峰是不是真实周期还需要一个参照背景——这就是后面第 4 章的显著性检验也是整条分析链路里最容易被误用的环节。3. Morlet 母函数气象周期分析默认选它的三个底层原因3.1 波包结构一支复指数振荡乘一个高斯衰减窗气象周期分析提到小波默认母函数几乎总是 Morlet。它的定义是ψ₀(η) π^{-1/4} e^{iω₀η} e^{-η²/2}结构非常直白一支复指数振荡 e^{iω₀η} 外面乘一个高斯衰减窗 e^{-η²/2}。在时域它就是一个有有限能量的「波包」中心是主振荡向两侧按高斯曲线快速衰减。这个形状天然适合周期分析一支纯正弦信号只有在局部窗口内和这个波包形状相似时内积才会给出大系数高斯窗提供时域定位复指数核提供频率选择二者合一。前面的常数 π^{-1/4} 是让母小波的能量归一为 1。这个细节直接影响功率谱的绝对值取掉归一化系数所有功率值整体缩放显著性阈值也必须跟着变。所以不要徒手改写母小波的幅度系数直接用工具箱内置实现。Morlet 还满足小波的「容许条件」吗严格说纯复指数的傅里叶变换在零频处不衰减到零需要加一个修正项当 ω0 足够大通常 ω05时修正项可以忽略。气象应用里 ω0 取 6早已远离容许条件的边界这也是它敢被当作工程默认值的原因。3.2 尺度到周期的换算1.03 是推导出来的不是拍脑袋对 Morlet 小波尺度 a 对应的傅里叶周期 T 有解析关系。ω06 时傅里叶等效频率约等于 ω0/(2πa)于是T 1/f ≈ 2πa/ω0 ≈ 1.03a这个 1.03 不是经验系数是小波基傅里叶变换峰值的位置决定的。改变 ω0换算系数跟着变。教材里写「Morlet 的 Fourier period ≈ 1.03a」都默认 ω06。MATLAB 里不必手算。老版本小波工具箱用 scal2frq(scales, morl, dt) 返回尺度对应的等效频率周期取倒数即可新版本 cwt 直接返回带单位的时间-频率网格。要注意该函数返回的是「傅里叶等效频率」不是小波中心频率两者之间的带宽展宽关系不同。见过不少人把尺度轴直接标成周期画完后发现 10 年周期的峰出现在尺度 8 的位置跑来追问是不是程序写错了——其实是换算没做。3.3 中心频率与时频分辨率何以 ω06 是气象折中Morlet 的自由参数是中心角频率 ω0。它同时控制高斯窗的衰减速度ω0 越大振荡次数越多小波基在时域越长频率分辨率越好而时间分辨率越差ω0 越小基函数越短时间定位越准但频率峰展宽相邻周期分不开。ω0时间定位频率分辨率典型场景2好差周期峰严重展宽突发事件的时频定位6中足够分辨相邻周期气象/气候周期分析默认20差接近傅里叶线谱只需要高精度频率读数海森堡不确定性原理在这背后是硬约束小波基在时域和频域的有效支撑面积不小于常数两者乘积有下界不存在让时间和频率同时达到最优的参数。气象里选 6是因为降水、气温信号的周期往往不是严格单频而是有一定带宽的过程——比如准两年振荡、8 年左右的年代际波动——既需要时间定位能力去找断点又需要把相邻周期分开。这个折中在公开文献里被检验得最多选 6 也意味着你的结果和前人工作有统一的比较基准。3.4 复数小波与相位还能判断两个要素谁先变Morlet 是复值母函数W(a,b) 的幅角就是该时频点的相位。做降水与气温的联合分析时交叉小波功率谱可以给出两者在某个时间段的相位差进而判断「升温发生在降水减少之前还是之后」这类因果先后问题。实数小波如墨西哥帽擅长检测波峰位置但没有方向区分能力给不出相位。这是 Morlet 霸榜气象周期分析的另一个原因单序列分析只用模值双序列分析要复用同一套系数拿相位。4. MATLAB 实现链路从尺度轴、功率谱到红噪声显著性检验4.1 老版本 cwt 与新版 cwt接口分水岭先搞清楚MATLAB 的连续小波变换接口在 R2016b 前后有一次大改网上大量示例混用两套语法抄错接口是最常见的事故源。老版写法是coefs cwt(s, scales, morl);输入的是尺度向量输出是系数矩阵行对应尺度、列对应时间。新版推荐做法是[wt, f] cwt(s, amor, VoicesPerOctave, 8);直接传小波名返回系数矩阵 wt 和频率向量 f尺度由工具箱内部按倍频程自动生成。两个版本在振幅归一化、尺度轴方向上有系统性差异拿到结果先 whos 检查维度再确认 f 是递增还是递减再画图。下面按老版语法给一条完整链路新版的差异我在参数说明里标注。dt 1/12; % 采样间隔月数据单位是年 N length(s); % s 为预处理后的序列 % 尺度轴每倍频程 8 个尺度范围 2^1 到 2^6 scales 2.^(1:1/8:6); % 老版指定尺度和小波名返回系数矩阵 coefs cwt(s, scales, morl); % 功率谱必须除以尺度做能量归一化 power abs(coefs).^2 ./ scales; % 尺度转周期scal2frq 第三参数 dt 决定频率单位 freq scal2frq(scales, morl, dt); period 1 ./ freq; % 周期的单位与 dt 一致这里是年 % 边界效应e-folding 时间Morlet 近似为 sqrt(2)*scale coi sqrt(2) * scales * dt; % 距左右边界的最小可信时间 % 绘图横轴时间纵轴周期大周期在上 contourf(t, period, power, 40, LineStyle, none); set(gca, YDir, reverse);代码里每一行都在对应一个设计决策。scales 用了对数均匀分布而不是线性均匀分布周期轴跨度从 2 到 64线性取点会让短周期段分辨率过剩、长周期段分辨率不足对数轴让每个倍频程有相同数量的分析尺度。scal2frq 返回的是等效频率单位由 dt 决定dt1/12 时频率单位是「每 1/12 年」period1./freq 就是「年」。COI 那句是 Morlet 的 e-folding 时间近似表示距边界多远之内的小波系数受截断污染绘图时通常用阴影盖住 COI 之外区域。4.2 尺度范围怎么设不是越大越好边界会吃掉一半尺度下限取多少取决于你要看的最短周期。月降水做准 2 年周期分析尺度下限 2 足够硬要往下到 1 以下对应周期不足一年高频段几乎全是季节噪声谱图会被一片细碎高能纹理淹没。尺度上限受数据长度约束COI 随尺度线性增长假设序列长 50 年尺度 64 对应的 COI 约 90 年边界效应吞掉的不止中点整条时间轴都不可信。经验上限是让最大尺度对应的 COI 小于序列长度的一半即最大周期不超过数据长度的约三分之一到二分之一具体可以拿 coi 向量对半程检查。Voices per octave每倍频程尺度数在 8 到 16 之间常用。取太少相邻尺度的功率峰之间出现缝隙窄带周期可能漏检取太多相邻尺度又高度相关只是让图变平滑不提供新信息。8 到 10 足够覆盖气象周期峰的半宽跑起来也快。4.3 红噪声显著性把伪周期挡在报告外气象序列几乎都带正自相关相邻时段的观测不是独立的。AR(1) 红噪声模型把这种惯性考虑进来作为周期显著性检验的背景。完整公式是先估计序列滞后 1 自相关系数 alpha再按频率计算归一化红噪声谱然后用卡方分布给阈值。% 红噪声背景AR(1) 模型alpha 为滞后1自相关系数 alpha corr(s(1:end-1), s(2:end)); % 归一化红噪声谱freq 来自 scal2frq Pn (1 - alpha^2) ./ ... (1 - 2*alpha*cos(2*pi*freq*dt) alpha^2); % 功率谱先除以序列方差得到标准化谱 power_norm power ./ var(s); % 自由度复小波模平方为 295% 置信 dof 2; threshold Pn * chi2inv(0.95, dof) / dof; sig_mask power_norm threshold;为什么要除以序列方差标准化后白噪声在各尺度的期望功率为 1红噪声背景 Pn 也是无量纲谱阈值公式才成立。dof 取 2 是因为复值小波系数有实部、虚部两个自由度若用实数小波如墨西哥帽则 dof 取 1对时间轴做了平滑后自由度还要相应放大。画图时用 sig_mask 叠加等值线只圈出显著区域比整片填色更能说明问题。5. 避坑记录周期分析最容易翻车的四个环节5.1 边界上的伪峰被当成真周期现象功率谱图的低频段边缘出现一片连续的高功率区峰值比全图中部还亮看起来像发现了一个「漫长的气候趋势周期」。原因小波基在数据边界外没有数据内积计算等价于给信号乘了一个跳变窗信号在边界处突然中断会产生虚假振荡。这个污染的尺度范围随尺度 a 增大而增大COI 之外的功率本身就是截断伪迹。解决画图时把 COI 外侧区域用半透明色打阴影寻找峰值只在 COI 内部进行。比较多个序列的谱时只取几个 COI 的交集作为可信区。这不是可选项是必须项边界伪峰的功率常常高于真实信号。5.2 白噪声背景做检验显著峰多出一倍现象同一序列用白噪声背景检验出 5 个「显著周期」换成红噪声只剩 1 个。原因气温、降水这类连续气象变量普遍有正自相关AR(1) 模型下低频背景功率天然偏高。白噪声背景谱是一条水平线等于否认了序列的惯性低频段所有峰都被错误抬升。解决显著性阈值必须基于 AR(1) 红噪声背景估计。alpha 从数据本身估计滞后 1 自相关系数注意它也有抽样误差样本太短时给置信区间留有裕量不要把刚好压线的峰当成板上钉钉的周期。5.3 尺度轴直接当周期轴峰位对不上现象图中标注「2 年周期」的峰值位置在尺度 3 或尺度 1.5 处横坐标对不上换一种母小波后峰还跑到别的位置。原因尺度 a 与傅里叶周期 T 的换算随母函数而定。Morletω06的 T≈1.03a 还算接近墨西哥帽小波的换算系数接近 4直接拿尺度当周期会把峰位标偏到四分之一处。解决换母函数或换工具箱版本后先做一次合成正弦检验对已知周期信号跑 scal2frq确认峰值落在期望周期上再画真实数据。这个检验耗时五分钟能省掉所有坐标轴争论。5.4 滑动平均之后做小波低频峰是加工出来的现象原序列没有明显长周期做了 3 年滑动平均后低频功率谱末端冒出一个显著峰。原因滑动平均本质是低通滤波抹掉短期波动的同时把相邻低频分量叠加放大三次样条插值则会在缺测区制造平滑连通人为引入低频假象。小波分析就是对连续尺度的相似度打分这类「加工痕迹」会直接变成谱上的峰。解决小波分析用原始采样序列缺测段不能合理处理就截断别用插值硬填。滑动平均只用来画趋势参考线不送进小波。若必须去趋势优先用一阶差分多项式拟合会把趋势残差拆成多个低频伪周期。6. 时频分辨率的取舍为什么分析结果不能直接复用别人的参数不确定性原理决定了小波分析里没有「唯一正确」的参数组合。尺度轴密度、母小波中心频率、COI 的判定方式都影响最终图的纹理。两个研究组用不同版本的工具箱对同一条降水序列做分析主周期位置肉眼几乎一致但显著区域的形状、起始年份常有差异。这不是工具出 bug而是时频分析的物理本性你选择的观测窗口不同看到的信息天然一详一略。换个尺度轴密度峰的位置不挪峰的边界一定挪。所以我处理真实数据前有一道强制流程先构造合成序列做参数校准——前半段 8 年周期、中段纯噪声、后半段 3 年周期加白噪声用自己刚要跑真实数据的参数走完整条链路看小波功率谱是否重建出「前 8 后 3」的时频图案看边界伪峰是否落在 COI 外看显著性掩膜是否只圈住已知频段。通过这一步尺度轴、母函数、显著性置信度三个参数才算被「钉」住了。再往上走一层单序列周期分析只是小波在气象里的切口。要回答「降水减少和气温升高谁先出现」需要把两条序列的小波系数做交叉功率谱和相干分析相位差随时间的变化图就是因果先后判断的依据。这套扩展的核心机制仍然是 COI、红噪声背景、显著性检验三件套参数校准方式可以完全复用。说来好笑我早期做周期分析踩过第五章里几乎所有的坑拿着白噪声背景报过一串伪周期把 COI 外的边界强峰写进过报告还因为没做尺度换算把峰位标错过快一倍。从那以后我每次拿到一条气象序列都会先用合成信号走一遍校准确认峰值落在预期周期上、伪峰全部被阈值挡掉才敢上真实观测数据。机械地跑一遍不费事但它能拦下绝大多数让你事后想删报告的错误。希望帮到你。本文还有配套的精品资源点击获取