简介：面向气象、气候与环境科学领域的学习者与研究人员，这份MATLAB小波功率谱分析资料聚焦降水、气温等时间序列的局部周期性与瞬态变化识别，适合已具备一定信号处理基础、希望将小波方法落地到实际气象数据的中高级用户。资源包共28个文件，约184KB，以15个m脚本为核心，配套png结果图、txt与dat数据文件及若干备份文件，涵盖Morlet小波变换、功率谱计算、显著性检验与多尺度可视化等环节，可直接运行并对照图示理解不同尺度上的功率分布。已有2149人学习下载，说明其在同类气象周期分析资料中具有较好的参考价值。读者可从中获得一套可复用的分析脚本与示例数据，借助波基选择、系数平方求谱、峰值识别等实现思路，快速完成从数据预处理到周期解释的完整流程，并迁移至极端天气探测、气候变化诊断等场景。
1. 小波功率谱做周期分析：一份能直接跑的 MATLAB 资源包
手里有一组几十年的月降水或日均气温序列，想搞清楚它到底存在哪些周期、这些周期在哪些年份强、哪些年份弱，FFT 只能给你一条全局频谱，答不了「1998 年那个 8 年周期还在不在」这种问题。小波功率谱就是干这个的：它把一维时间序列摊到时频平面上，横轴是年份、纵轴是周期尺度，颜色深浅代表该周期在该时段的能量强弱，再叠一条红噪声显著性检验的等值线，就能判断哪些周期是「真信号」、哪些只是随机涨落的玄学。这份 MATLAB 资源包围绕降水、气温这类气象水文序列的周期分析展开，包含连续小波变换主函数、Morlet 母小波生成、红噪声背景谱与显著性检验、锥形影响区（COI）绘制和功率谱出图脚本，适合做气候诊断、水文序列周期识别的研究生和一线业务人员直接改数据复现。下面按「先立住原理、再动手跑通、最后避坑」的顺序拆开讲。
2. 小波功率谱的原理与母小波选型：为什么是 Morlet 而不是 Mexican Hat
2.1 从 FFT 的局限说到时频局部化
傅里叶变换把整条序列投影到一组全局正弦基上，得到的是「整个时段平均下来有哪些频率成分」，时间信息被完全积分掉了。对于降水、气温这种非平稳序列——周期振幅随年代变化、甚至周期本身在漂移——全局谱会给出误导性的结论：一个只在 1970—1985 年活跃的 6 年周期，在 FFT 谱上照样是一个尖峰，你无法区分它是持续存在还是阶段性出现。短时傅里叶变换（STFT）加窗能补回时间信息，但窗宽固定，低频要宽窗、高频要窄窗，二者不可兼得。
连续小波变换（CWT）的思路是用一组可伸缩、可平移的母小波去卷积原序列：尺度参数 a 控制小波被拉伸的程度（对应周期长短），平移参数 b 控制小波在时间轴上的位置。尺度大时小波被拉宽，频率分辨率高、时间分辨率低；尺度小时反过来。这正好契合气象序列「低频周期要看长时段、高频振荡要看短时段」的需求，也是这份资源包选择 CWT 而非 STFT 的根本理由。
离散形式下，对长度为 N、采样间隔为 dt 的序列 x_n，尺度 a 上的小波系数写成卷积：
W(a, b) = Σ x_n · ψ*[(n − b)·dt / a]
功率谱就是 |W(a,b)|²。资源包里的主函数做的就是这个卷积，只不过用 FFT 在频域实现以提速——直接时域卷积对 N=600、尺度 40 档的序列要跑几十秒，频域实现能压到一秒以内，这是实操中必须知道的性能点。
2.2 Morlet 母小波的复值与相位优势
母小波的选择直接决定谱的形态。常用候选有 Morlet、Mexican Hat、Paul 等，这份资源包默认用复 Morlet，原因有三条，都是踩过坑之后的选择：
第一，Morlet 是复小波，小波系数同时含模和相位。模的平方给出功率谱，相位可以用来做交叉小波、相干性分析，后续想扩展到「降水和气温的周期是否同步」时不用换工具。第二，Morlet 在时频域都是高斯包络，局部化性质好，谱图不会出现明显的边缘泄漏。第三，Morlet 的尺度与傅里叶周期之间有近似线性关系，尺度转周期只需乘一个常数，出图时纵轴标注直观。
复 Morlet 的标准形式是 ψ(t) = π^(−1/4) · exp(iω₀t) · exp(−t²/2)，其中 ω₀ 是无量纲中心频率，资源包里默认取 6。这个 6 不是随便定的：ω₀=6 时 Morlet 小波满足容许性条件近似成立（严格说 ω₀≥5 才够），同时时频分辨率折中较好。ω₀ 调大，频率分辨率升高但时间分辨率下降，谱图上的周期带会变窄变「锐」；调小则相反。做年际周期（2—8 年）分析时 6 是稳妥值，若专门盯 ENSO 相关的 3—7 年带，可以试 5 到 8 之间扫一遍看稳定性。
尺度与周期的换算关系是资源包里一个容易忽略但很关键的参数：
| 参数 | 含义 | 默认值 | 调整建议 |
|---|---|---|---|
| omega0 | Morlet 中心频率 | 6 | 年际周期分析保持 6，专盯窄带可试 5—8 |
| dt | 采样间隔（年） | 1 | 月数据填 1/12，务必与序列对齐 |
| dj | 尺度离散化步长 | 0.25 | 越小谱越平滑但计算量上升，0.125—0.25 常用 |
| s0 | 最小尺度 | 2·dt | 一般取 2 倍采样间隔，低于此受 Nyquist 限制 |
| sj | 最大尺度 | 由序列长度定 | 常取 N/2 或更小，超过 N/3 后 COI 覆盖严重 |
尺度序列按 s_j = s0 · 2^(j·dj) 生成，j 从 0 到 J，J = (1/dj)·log2(sj/s0)。这套几何级数采样保证每个倍频程内有 1/dj 个尺度，dj=0.25 即每倍频程 4 个尺度，谱图纵向不会出现台阶感。周期换算用 Fourier 周期 λ = 4π·a / (ω₀ + √(2+ω₀²))，ω₀=6 时约等于 1.03a，工程上直接当 a 用误差不大，但写论文时建议用精确式。
2.3 红噪声背景谱与显著性检验的逻辑
小波功率谱最容易翻车的地方不是算谱，是判显著。原始功率谱上到处都是「峰」，不做检验就宣称「存在 5 年周期」是典型误用。资源包的做法是：先估计序列的红噪声背景谱，再对每个尺度算 95% 置信水平的理论功率，最后在谱图上画等值线，只有落在等值线以内的区域才算显著。
红噪声假设序列是一阶自回归过程 AR(1)：x_t = α·x_{t−1} + ε_t，α 由 lag-1 自相关系数估计。气象水文序列普遍存在这种持续性，用白噪声做背景会高估显著性。资源包里背景谱的计算按 Torrence & Compo 的经典公式，对每个尺度 a 算理论红噪声功率，再乘卡方分布的 95% 分位系数。实操中要注意：α 的估计要用去趋势后的序列，若原序列有明显线性趋势不先去掉，α 会被高估，背景谱抬高，真周期反而被压到不显著——这是血泪经验里最常见的一条。
3. 把数据喂进小波函数：从原始序列到功率谱矩阵的四步
3.1 数据预处理与去趋势
拿到手的降水或气温序列，第一步不是直接进小波函数，而是补齐、去趋势、标准化。资源包假定输入是一维列向量，缺失值用线性插值补，然后减去最小二乘拟合的线性趋势。去趋势不是可选项：气候序列普遍有增温趋势或年代际漂移，不去趋势会让低频段功率虚高，谱图底部一片红，什么周期都看不出来。
% 读入序列并预处理：补缺、去趋势、标准化 raw = readmatrix('precip_monthly.csv'); % 单列，时间升序 x = raw(:); % 线性插值补缺（NaN 位置） idx = 1:numel(x); x = interp1(idx(~isnan(x)), x(~isnan(x)), idx, 'linear', 'extrap'); % 去线性趋势 t = (1:numel(x))'; p = polyfit(t, x, 1); x = x - polyval(p, t); % 标准化，便于不同量纲序列对比 x = (x - mean(x)) / std(x);这段代码的逻辑是：插值保证序列连续（小波卷积对 NaN 极其敏感，一个 NaN 会污染整段系数）；去趋势用一阶多项式拟合后相减，保留年际波动、去掉长期漂移；标准化让功率谱的量级与序列单位无关，方便多站点对比。参数上，interp1的'extrap'只在序列首尾缺测时起作用，若中间缺测超过连续 12 个点，建议整段剔除而不是硬插，插出来的假周期比缺数据更麻烦。
3.2 调用主函数生成小波系数
资源包的核心函数签名大致是[wave, period, scale, coi] = wavelet_ps(x, dt, dj, s0, sj, omega0)，返回小波系数矩阵 wave（复数，行对应尺度、列对应时间）、对应的 Fourier 周期向量 period、尺度向量 scale 和锥形影响区 coi。调用时最容易错的是 dt：月数据必须填 1/12，填 1 会让所有周期标签错位 12 倍，谱图形状对但纵轴全错，这种错误肉眼很难发现，因为图看起来「很正常」。
dt = 1/12; % 月数据，单位为年 dj = 0.25; % 每倍频程 4 个尺度 s0 = 2 * dt; % 最小尺度 sj = numel(x) * dt / 2; % 最大尺度取序列长度的一半（年） [wave, period, scale, coi] = wavelet_ps(x, dt, dj, s0, sj, 6); power = abs(wave).^2; % 功率谱这里sj取序列总时长的一半是经验值：尺度超过序列长度 1/3 后，卷积需要的数据超出序列边界，COI 会覆盖大部分区域，谱值不可信。若序列只有 30 年，最大可信周期约 10 年，想看 20 年以上的年代际周期，数据长度至少 60 年，这是硬约束，不是调参能解决的。
3.3 红噪声检验与显著性等值线
检验部分资源包单独封装，输入原始序列和尺度向量，输出每个尺度上的 95% 理论功率。核心是估 α、算背景谱、乘分位系数三步。α 用 lag-1 自相关估，注意用去趋势后的序列；背景谱按尺度逐点算；分位系数用卡方近似，自由度与尺度有关。
% 红噪声背景与 95% 显著性阈值 alpha = abs(atan(imag(...))); % 实际用 lag-1 自相关，此处示意 % 资源包内封装为： [signif, fft_theor] = rednoise_signif(x, dt, scale, 0.95); % 画等值线时用 signif 与 power 比较 contour(time, period, power ./ signif, [1 1], 'k', 'LineWidth', 1.2);逻辑说明：power ./ signif是功率与阈值的比值，等值线取 1 的位置，比值大于 1 的区域即显著。用比值而非直接比较，是为了让等值线在任意功率量级下都能画。参数上，0.95 是置信水平，做严格诊断可提到 0.99，但气象序列样本量有限，0.99 往往只剩零星显著点，实际报告里 0.95 更常用。注意rednoise_signif内部对 α 做了截断（限制在 0 到 0.99 之间），防止序列接近白噪声时 α 估计溢出。
3.4 出图：功率谱、COI 与显著性等值线三件套
一张合格的小波功率谱图必须同时有：填色功率谱、黑色显著性等值线、COI 边界线（通常再加一层半透明遮罩）。缺任何一样，图都不能直接进论文。资源包的绘图脚本用imagesc或contourf打底，contour叠等值线，plot画 COI 两条斜线，最后fill把 COI 以外区域遮灰。
figure; contourf(time, period, log2(power), 32, 'LineStyle', 'none'); hold on; contour(time, period, power ./ signif, [1 1], 'k', 'LineWidth', 1.2); plot(time, coi, 'k--', 'LineWidth', 1.5); % COI 边界 fill([time fliplr(time)], [coi fliplr(max(period)*ones(size(coi)))], ... [0.5 0.5 0.5], 'FaceAlpha', 0.4, 'EdgeColor', 'none'); set(gca, 'YDir', 'reverse', 'YScale', 'log'); ylabel('Period (year)'); xlabel('Time (year)');关键点在log2(power)和YScale log：功率谱动态范围大，线性色标下高频细节全被低频大值压掉，取对数后年际尺度的小结构才看得见；纵轴取对数是因为尺度本身按几何级数生成，线性纵轴会让低频段被拉得过宽。COI 遮罩用fill加FaceAlpha实现半透明，比直接裁掉数据更直观——读者能看到边界外的谱形，但知道那部分不可信。
4. 避坑与排查：小波功率谱五个高频翻车点
4.1 纵轴周期整体错位
现象：谱图形状合理，但标注的周期与实际对不上，比如明明 3—7 年带最亮，纵轴却标在 30—70 年。原因几乎都是 dt 单位填错，月数据填了 1 而不是 1/12，导致所有尺度换算成周期时放大 12 倍。解决：进函数前打印dt和period(1:5)核对，月数据的 period 首值应在 0.2 年左右，年数据在 2 年左右，对不上就是单位问题。
4.2 显著性等值线一条都画不出来
现象：contour报错或图上没有黑线。原因通常是power ./ signif的最大值都小于 1，即整条序列没有任何点超过 95% 阈值。这未必是代码错，可能是序列确实以噪声为主，也可能是 α 被高估导致背景谱过高。解决：先打印max(power(:)./signif(:))，若在 0.8 附近说明接近显著但没过线，可检查是否忘了去趋势；若远小于 1，考虑序列太短或该段本无稳定周期，换更长序列再试。
4.3 COI 遮罩把有效区域也盖住了
现象：图上大半区域被灰色遮住，只剩中间一条窄带。原因是最尺度 sj 取得太大，接近甚至超过序列长度，COI 在两端迅速抬升。解决：把 sj 压到序列时长的 1/3 以内，或换更长的数据。经验规则是：想可信地看 P 年周期，序列至少要有 3P 年长度，30 年数据看 10 年周期已是上限。
4.4 复 Morlet 的相位被当成功率用
现象：后续做交叉小波时结果全乱。原因是有时直接对wave取绝对值平方前忘了wave是复数矩阵，或者误把实部当功率。解决：功率一律用abs(wave).^2，相位用angle(wave)，两者不要混。资源包里主函数返回的是复系数，这一点在函数头注释里写明了，改代码时别顺手real()。
4.5 多站点批量跑时内存爆掉
现象：循环处理上百个站点时 MATLAB 卡死或报内存不足。原因是每个站点都保留了完整的 wave 复数矩阵（尺度数 × 时间长度），上百站点叠加很快到 GB 级。解决：循环内只保留power和signif的降采样版本或统计量（如各周期带的平均功率），wave用完即clear；或把尺度数 dj 调大到 0.5 减少矩阵行数。批量任务建议分块存盘，别指望一次全在内存里。
5. 进阶用法：交叉小波看降水与气温的周期同步性
单序列功率谱只能回答「我自己有哪些周期」，实际诊断里更常问的是「降水的 5 年周期和气温的 5 年周期是不是同一时段、同相位」。资源包的复 Morlet 系数正好支持交叉小波：对两条等长序列分别算 wave1、wave2，交叉小波谱为 wave1 .* conj(wave2)，其模平方是交叉功率，相位差是 angle(wave1 .* conj(wave2))。相位差接近 0 表示同相，接近 π 表示反相，这对判断「暖年降水偏少」这类关系在哪个周期带上成立很有用。
% 交叉小波：降水 x1 与气温 x2 [wave1, period, scale, coi] = wavelet_ps(x1, dt, dj, s0, sj, 6); [wave2, ~, ~, ~] = wavelet_ps(x2, dt, dj, s0, sj, 6); xwt = wave1 .* conj(wave2); % 交叉小波系数 xpow = abs(xwt); % 交叉功率 xphase = angle(xwt); % 相位差 % 交叉功率显著性用两序列各自红噪声谱的乘积做背景 signif_x = signif1 .* signif2; contourf(time, period, log2(xpow), 32, 'LineStyle', 'none'); hold on; contour(time, period, xpow ./ signif_x, [1 1], 'k', 'LineWidth', 1.2);参数上要注意两点：两条序列必须等长、同采样、同去趋势方式，否则相位差没有意义；交叉功率的显著性阈值是两条序列各自红噪声谱的乘积，不是简单取平均，资源包里rednoise_signif对单序列返回的 signif 直接相乘即可。相位箭头图（quiver）可选，但在周期带密集时容易糊成一团，我一般只在论文插图里对显著区域抽稀画箭头，正文分析用相位差的直方图更清楚。
验证方法上，我习惯做一次「已知周期注入」测试：往白噪声里塞一个振幅 0.5、周期 6 年的正弦，跑一遍完整流程，看谱峰是否落在 6 年、显著性等值线是否包住注入段。这一步能同时验证 dt、尺度换算、显著性检验三处是否一致，比盯着真实数据猜靠谱得多。从那以后我每次换数据集或改母小波参数，都强制先跑一遍这个注入测试再出正式图，省下的返工时间远超那几分钟。希望帮到你。
本文还有配套的精品资源，点击获取
