简介这份资源面向气象学、气候学及环境科学领域的学习者与研究者提供基于MATLAB的小波功率谱分析工具用于揭示降水、气温等时间序列中的局部周期性与瞬态变化特征。包内共28个文件以m脚本为主辅以png结果图、txt数据说明、dat与powerspec数据文件及bak备份压缩包约184KB涵盖小波变换、功率谱计算、显著性检验与可视化等完整流程。资源包含Morlet小波、连续与离散变换、功率谱密度计算及周期峰值识别等核心环节并配有太阳黑子、Nino3海温等示例数据与多组绘图脚本便于读者直接运行、对照结果并迁移到自身研究数据中。已有2147人学习下载适合希望掌握非平稳信号周期分析、提升气象数据处理能力的中高级用户参考。1. 小波功率谱做周期分析一份能直接跑的 MATLAB 资源做气象或水文时间序列的人大概都遇到过这种场面手里攥着某站点几十年的月降水或日均气温画出来一条上下翻飞的曲线肉眼能看出有年循环但要说清楚“除了年周期还有没有 3 到 7 年的振荡、这个振荡在哪几年特别强”光靠傅里叶变换是不够的——傅里叶只告诉你整个时段平均有哪些频率成分不告诉你这些成分什么时候出现、什么时候消失。小波功率谱就是干这个的它把一维时间序列摊到“时间-尺度”二维平面上让你看到不同周期信号的强度随时间怎么变。这份 MATLAB 资源围绕降水、气温这类气象序列的周期分析展开核心是小波变换加功率谱计算附带显著性检验和边界效应处理。适合手头有长序列观测数据、想快速出图出结论的从业者也适合刚接触小波分析、需要一个能跑通模板的新手。下面按“原理选型 → 代码落地 → 避坑排查 → 进阶技巧”的顺序拆开讲。2. 小波功率谱的选型逻辑为什么不用傅里叶为什么选 Morlet2.1 傅里叶变换在非平稳序列上的硬伤气象序列几乎都是非平稳的。年降水不会每年一模一样厄尔尼诺相关的几年振荡也不是全程恒定。傅里叶变换把整个时间轴积分掉得到的是全局频谱等于把“什么时候发生”这个信息彻底抹平。举个具体场景某站 1960—2020 年的降水序列傅里叶谱在 35 年处有个峰但你没法判断这个峰是 1970 年代贡献的还是 2000 年代贡献的。如果要做年代际对比傅里叶直接出局。短时傅里叶变换STFT加窗能部分解决这个问题但窗长固定窗短了频率分辨率差窗长了时间分辨率差二者不可兼得。小波变换用可伸缩的母小波替代固定窗低频用宽窗、高频用窄窗天然适配“低频周期长、高频周期短”的气象信号特征。这是选小波而不是 STFT 的根本理由。2.2 母小波选 Morlet 的工程理由母小波有很多种Morlet、Mexican hat、Paul、DOG 等。做周期分析绝大多数气象文献用 Morlet原因有三条。第一Morlet 是复小波能同时给出振幅和相位振幅平方就是功率谱相位可以用来判断两个序列的滞后关系。第二Morlet 在频率域近似高斯型频率局部化好适合识别“某个周期附近的能量”。第三Morlet 的尺度与傅里叶周期之间有近似线性关系换算方便。Morlet 小波的定义里有个中心频率参数 ω₀常见取值 6。ω₀ 越大频率分辨率越高、时间分辨率越低ω₀ 越小则相反。ω₀6 是气象领域的经验平衡点满足容许性条件的同时兼顾两类分辨率。如果你做的是高频突变信号比如逐小时降水可以降到 5 甚至更低做年代际振荡可以升到 8。这个参数在代码里通常写作w0或omega0。2.3 功率谱、显著性检验与边界效应三件套小波功率谱本身只是|W|²但一张能发文章、能写报告的图必须叠加两样东西显著性检验线和影响锥COI。显著性检验的零假设通常是“红噪声背景谱”因为气象序列普遍存在红噪声特征低频能量天然偏高。检验方法是用蒙特卡洛或卡方分布构造 95% 置信线功率超过这条线的区域才算“显著周期”。影响锥是边界效应区小波在序列两端做卷积时数据不足边缘区域的功率不可信COI 曲线以下的区域要打折扣看。这两样缺一个结论就可能翻车。常见做法是先对序列做距平或标准化再算小波系数然后逐尺度计算红噪声背景谱和 95% 阈值最后画填色功率谱 黑色等值线 COI 曲线 显著性线。下面进入代码落地。3. 从原始序列到功率谱图MATLAB 代码逐段拆解3.1 数据准备与预处理拿到一份降水或气温的月值序列第一步不是直接扔进小波函数而是检查缺失值、去趋势、标准化。缺失值用线性插值补齐去趋势用最小二乘拟合线性项后减去标准化用(x-mean)/std。这三步不做功率谱会被趋势项污染低频区出现虚假高能。% 假设原始数据为列向量 x时间单位为月 % 1. 缺失值线性插值 nanIdx isnan(x); x(nanIdx) interp1(find(~nanIdx), x(~nanIdx), find(nanIdx), linear, extrap); % 2. 去线性趋势 t (1:length(x)); p polyfit(t, x, 1); x_detrend x - polyval(p, t); % 3. 标准化 x_norm (x_detrend - mean(x_detrend)) / std(x_detrend);逻辑说明interp1的extrap选项处理序列首尾缺失避免插值后仍有 NaN。polyfit一次多项式拟合趋势polyval还原后相减。标准化让不同站点的功率谱可以横向比较。参数上如果序列本身已经是距平值第 2 步可跳过如果做的是气温去趋势尤其重要因为增温趋势会在低频区造出一个假峰。3.2 小波变换核心计算MATLAB 没有内置的连续小波变换函数直接输出 Morlet 功率谱cwt在新版本里改了接口所以常见做法是自己写卷积或者用 FFT 实现。下面用 FFT 方式效率高且尺度可自定义。% 参数设置 dt 1; % 采样间隔月值序列取1 s0 2*dt; % 最小尺度 dj 0.25; % 尺度间隔 J floor(log2(length(x_norm)*dt/s0)/dj); % 最大尺度层数 scales s0 * 2.^(dj*(0:J)); % 尺度序列 w0 6; % Morlet中心频率 % FFT准备 n length(x_norm); x_fft fft(x_norm); omega 2*pi*[0:floor(n/2), -ceil(n/2)1:-1]/n; % 角频率 % 逐尺度计算小波系数 wavelet_coef zeros(length(scales), n); for a 1:length(scales) s scales(a); % Morlet小波在频率域的表达式 psi_hat pi^(-1/4) * (omega - w0) .* exp(-(s*omega - w0).^2 / 2) .* (omega 0); % 卷积等价于频域相乘 wavelet_coef(a,:) ifft(x_fft .* psi_hat); end % 功率谱 power abs(wavelet_coef).^2;逻辑说明s0是最小分析尺度月值序列取 2 个月dj控制尺度离散程度0.25 是常用值越小越平滑但计算量越大J决定最大尺度受序列长度限制一般不超过log2(N)。psi_hat是 Morlet 在频率域的解析形式(omega 0)保证只取正频率。ifft后得到复数小波系数取模平方得功率。参数上w0改大则频率分辨率提高dj改小则尺度更密。3.3 显著性检验与影响锥这一步是区分“能看的图”和“能写结论的图”的关键。红噪声背景谱用 AR(1) 模型拟合序列的自相关系数然后按卡方分布算 95% 阈值。% 拟合AR(1)系数 r1 corrcoef(x_norm(1:end-1), x_norm(2:end)); alpha r1(1,2); % 逐尺度计算红噪声背景谱和95%阈值 variance 1; % 标准化后方差为1 background zeros(size(scales)); threshold zeros(size(scales)); for a 1:length(scales) s scales(a); % 红噪声背景谱 background(a) variance * (1 - alpha^2) / (1 - 2*alpha*cos(2*pi*dt/s) alpha^2); % 95%卡方阈值自由度为2 threshold(a) background(a) * chi2inv(0.95, 2) / 2; end % 影响锥边界处尺度对应的可信边界 coi 1 ./ scales * sqrt(2) * dt; % 简化形式实际按e折时间计算逻辑说明alpha是 lag-1 自相关系数反映红噪声强度。background公式是 AR(1) 的理论谱。chi2inv(0.95,2)是自由度为 2 的卡方分布 95% 分位数除以 2 是因为功率谱的实部和虚部各贡献一个自由度。coi这里给的是简化计算严谨做法是按sqrt(2)*s*dt逐尺度算 e 折时间。参数上alpha越接近 1红噪声背景越高显著区域越少这是符合物理的——强持续性序列本来就难识别出独立周期。3.4 出图与周期换算最后把功率谱画成填色图叠加显著性线和 COI再把尺度轴换算成周期轴。% 尺度转周期 period 4*pi*scales / (w0 sqrt(2 w0^2)); % 绘图 figure; contourf(1:n, period, power, 20, LineColor, none); hold on; plot(1:n, coi, k--, LineWidth, 1.5); % COI曲线 contour(1:n, period, power, [1 1]*mean(threshold), k, LineWidth, 1); % 显著性线 set(gca, YScale, log); ylabel(周期 (月)); xlabel(时间 (月)); colorbar;逻辑说明period换算公式来自 Morlet 小波的尺度-周期关系w06时约等于scales*1.03。contourf画填色contour画显著性等值线coi用虚线标出。Y 轴取对数让周期分布更均匀。参数上填色层数 20 是经验值太少粗糙、太多拖慢渲染显著性线取mean(threshold)是为了画一条近似水平线严谨做法是逐尺度画曲线。4. 避坑与排查小波功率谱最容易翻车的五个地方4.1 现象低频区一大片高能但全是假的原因序列没去趋势或者去趋势不彻底。气温序列的增温趋势、降水序列的长期变化都会在最大尺度附近堆出高功率看起来像“几十年周期”实际是趋势项。解决做小波前强制去线性趋势必要时去二次趋势出图后检查最大尺度处的功率是否沿时间轴均匀分布均匀分布基本就是趋势。4.2 现象边缘区域出现奇怪的竖条纹原因影响锥没画或者没遮。小波在序列两端数据不足卷积结果不可信但功率值照样算出来画出来就是边缘高能假象。解决出图必须叠加 COI 曲线COI 以下区域要么打阴影要么直接不解读。常见做法是用patch函数把 COI 以下区域填成半透明灰色。4.3 现象显著性检验线画出来显著区域大得离谱原因红噪声背景谱的alpha估计错了。如果序列本身已经标准化且去趋势alpha应该重新算不能沿用原始序列的。另外自由度取 2 是复数小波的标准做法如果误用自由度 1阈值会偏低显著区域虚大。解决检查alpha是否在 0 到 1 之间检查chi2inv的自由度参数是否为 2。4.4 现象换一段序列周期位置对不上原因尺度序列scales的dj或s0改了但周期换算公式没同步。period依赖scales和w0任何一个变了都要重算。解决把s0、dj、w0、period换算写成一个函数改参数时统一调用不要手抄公式。4.5 现象代码跑得动但功率谱全是 NaN原因序列里有 Inf 或者复数。降水数据一般不会但如果是比值序列或者做了对数变换可能出现 Inf。解决预处理阶段加一句x(~isfinite(x)) NaN再走插值流程。另外检查fft输入是否全为实数复数输入会导致ifft结果异常。5. 进阶技巧多序列交叉小波与尺度平均谱5.1 交叉小波功率谱看两个序列的共变周期单序列小波功率谱只能看“自己有什么周期”如果想知道“降水和气温在 35 年周期上是不是同步”需要交叉小波。做法是把两个序列各自做小波变换然后计算W_xy W_x * conj(W_y)取模平方得交叉功率谱相位角angle(W_xy)反映滞后关系。MATLAB 里把第 3.2 节的循环跑两遍存两个wavelet_coef再逐点相乘即可。相位角解读时注意正值表示 x 超前 y负值相反单位是弧度换算成时间要除以对应尺度的角频率。% 假设已有 wavelet_coef_x 和 wavelet_coef_y cross_power abs(wavelet_coef_x .* conj(wavelet_coef_y)).^2; phase angle(wavelet_coef_x .* conj(wavelet_coef_y)); % 相位转滞后时间以月为单位 lag_time phase ./ (2*pi ./ period(:));参数说明conj取共轭angle返回 -π 到 π 的相位。lag_time逐尺度逐时间点算正值表示 x 超前。注意交叉功率谱的显著性检验更复杂常用蒙特卡洛模拟代码量比单序列大建议先跑通单序列再上交叉。5.2 尺度平均谱把二维图压成一维曲线二维功率谱信息量大但不好做定量对比。常见做法是对时间轴取平均得到“全局小波谱”横轴周期、纵轴平均功率形状类似傅里叶谱但保留了尺度局部化优势。再进一步可以对某个周期带比如 37 年做时间平均得到该带的平均能量随时间的变化用来做年代际对比。% 全局小波谱 global_spectrum mean(power, 2); % 3-7年周期带月值36-84月的平均能量随时间变化 band_idx period 36 period 84; band_energy mean(power(band_idx, :), 1); plot(1:n, band_energy);逻辑说明mean(power,2)对时间维平均得到每个尺度的平均功率。band_idx选出目标周期带mean(...,1)对尺度维平均得到该带能量随时间的曲线。参数上周期带的选取要结合研究区域气候背景37 年是 ENSO 典型带10 年以上是年代际带不要随意切。5.3 一个我常犯的错早些年做小波分析图一出来看到显著区域就兴奋直接写“存在 5 年显著周期”结果审稿人一句“COI 以下区域怎么解释”就卡住了。后来养成习惯出图后先看显著区域是否落在 COI 以上再看它是否连续跨越多个时间段最后才下结论。还有一个血泪经验是dj设得太小尺度密到 200 层图是平滑了但计算量翻倍而且相邻尺度高度相关显著性检验会偏乐观。现在固定dj0.25够用且不拖。每次改参数我都强制把s0、dj、w0、period四者一起检查一遍确认换算没断链再出图。希望帮到你。本文还有配套的精品资源点击获取