简介：这份资源面向气象学、气候学及环境科学领域的学习者与研究者，提供基于MATLAB的小波功率谱分析工具，用于揭示降水、气温等时间序列中的局部周期性与瞬态变化特征。包内共28个文件，以m脚本为主，辅以png结果图、txt数据说明、dat与powerspec数据文件及bak备份，压缩包约184KB，涵盖小波变换、功率谱计算、显著性检验与可视化等完整流程。资源包含Morlet小波、连续与离散变换、功率谱密度计算及周期峰值识别等核心环节，并配有太阳黑子、Nino3海温等示例数据与多组绘图脚本，便于读者直接运行、对照结果并迁移到自身研究数据中。已有2147人学习下载，适合希望掌握非平稳信号周期分析、提升气象数据处理能力的中高级用户参考。
1. 小波功率谱做周期分析：一份能直接跑的 MATLAB 资源
做气象或水文时间序列的人，大概都遇到过这种场面：手里攥着某站点几十年的月降水或日均气温，画出来一条上下翻飞的曲线，肉眼能看出有年循环，但要说清楚“除了年周期，还有没有 3 到 7 年的振荡、这个振荡在哪几年特别强”，光靠傅里叶变换是不够的——傅里叶只告诉你整个时段平均有哪些频率成分，不告诉你这些成分什么时候出现、什么时候消失。小波功率谱就是干这个的：它把一维时间序列摊到“时间-尺度”二维平面上，让你看到不同周期信号的强度随时间怎么变。这份 MATLAB 资源围绕降水、气温这类气象序列的周期分析展开，核心是小波变换加功率谱计算，附带显著性检验和边界效应处理。适合手头有长序列观测数据、想快速出图出结论的从业者，也适合刚接触小波分析、需要一个能跑通模板的新手。下面按“原理选型 → 代码落地 → 避坑排查 → 进阶技巧”的顺序拆开讲。
2. 小波功率谱的选型逻辑：为什么不用傅里叶，为什么选 Morlet
2.1 傅里叶变换在非平稳序列上的硬伤
气象序列几乎都是非平稳的。年降水不会每年一模一样，厄尔尼诺相关的几年振荡也不是全程恒定。傅里叶变换把整个时间轴积分掉，得到的是全局频谱，等于把“什么时候发生”这个信息彻底抹平。举个具体场景：某站 1960—2020 年的降水序列，傅里叶谱在 3～5 年处有个峰，但你没法判断这个峰是 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 小波的尺度-周期关系，w0=6时约等于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 交叉小波功率谱看两个序列的共变周期
单序列小波功率谱只能看“自己有什么周期”，如果想知道“降水和气温在 3～5 年周期上是不是同步”，需要交叉小波。做法是把两个序列各自做小波变换，然后计算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 尺度平均谱：把二维图压成一维曲线
二维功率谱信息量大但不好做定量对比。常见做法是对时间轴取平均，得到“全局小波谱”，横轴周期、纵轴平均功率，形状类似傅里叶谱但保留了尺度局部化优势。再进一步，可以对某个周期带（比如 3～7 年）做时间平均，得到该带的平均能量随时间的变化，用来做年代际对比。
% 全局小波谱 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)对尺度维平均，得到该带能量随时间的曲线。参数上，周期带的选取要结合研究区域气候背景，3～7 年是 ENSO 典型带，10 年以上是年代际带，不要随意切。
5.3 一个我常犯的错
早些年做小波分析，图一出来看到显著区域就兴奋，直接写“存在 5 年显著周期”，结果审稿人一句“COI 以下区域怎么解释”就卡住了。后来养成习惯：出图后先看显著区域是否落在 COI 以上，再看它是否连续跨越多个时间段，最后才下结论。还有一个血泪经验是dj设得太小，尺度密到 200 层，图是平滑了，但计算量翻倍，而且相邻尺度高度相关，显著性检验会偏乐观。现在固定dj=0.25，够用且不拖。每次改参数，我都强制把s0、dj、w0、period四者一起检查一遍，确认换算没断链，再出图。希望帮到你。
本文还有配套的精品资源，点击获取
