简介：面向气象、海洋与气候研究者的MATLAB小波功率谱分析工具包，聚焦降水、气温等时间序列的周期识别与突变检测。小波功率谱结合小波变换的多分辨率优势与功率谱密度，可有效刻画非平稳信号中的局部周期、年际变化及瞬态事件，特别适用于季节振荡、ENSO信号等非固定周期现象的诊断。压缩包共28个文件，以15个.m功能脚本为主，包含连续小波变换、显著性检验、功率谱计算等核心流程，并配有txt/dat格式的示例观测数据、png结果图、bak/asv备份及rar数据存档，整包仅184KB，轻巧便于调试。已有2149人学习下载。除通用小波基函数wave_bases.m与核心分析例程wavelet.m外，还提供太阳黑子数、Nino3海温、TSI总太阳辐照度、气温与降水等多类真实数据驱动案例，可直接运行出含95%置信检验的小波功率谱图，便于对照理解尺度—周期对应关系，也可替换数据快速用于自选站点分析。
1. 小波功率谱做周期分析：为什么我劝你别再用滑动平均盯降水序列
做气候诊断或者水文分析时，周期提取是个绕不开的环节。很多人习惯先用滑动平均把序列磨平，再用傅里叶变换找显著周期，结果常常是主周期没找到，倒是在低频端凭空多出几个“假周期”。我自己拿到某流域近50年降水数据时，也走过这段弯路。后来换成连续小波变换（CWT）的功率谱来做周期分析，同一组数据跑下来，主周期信号清晰得多，而且能直接看到周期随时间的变化——这是傅里叶变换给不了的。这套基于MATLAB的小波功率谱分析资源，核心就是把降水、气温这类气候序列的周期结构拆开，用Morlet小波基计算功率谱，再叠加显著性检验，把“哪些周期是真的、哪些只是随机扰动”这件事说清楚。适合正在做气候诊断、水文趋势分析、生态遥感时序处理的人，也适合被导师催着交图的研究生。
2. 连续小波变换的原理与选型：Morlet为什么是默认选项
2.1 小波变换与傅里叶变换的本质差别
傅里叶变换把时间序列拆成无限长的正弦波，得到的是全局频谱——它假设信号在整段时间内统计特性不变。可降水序列偏偏是典型的非平稳信号，旱季雨季交替、年代际波动叠加，同一个周期在1970年代可能显著，到了2000年代就消失了。傅里叶变换对此毫无办法，它会把这些时变特性平均掉，结果就是频谱上出现一堆莫名其妙的宽峰。
小波变换的核心思路是用有限长度的波包（小波基）去匹配信号。这个小波基可以沿时间轴平移，也可以伸缩尺度。平移对应时间定位，伸缩对应频率分析。所以小波变换输出的是一个二维矩阵：横轴是时间，纵轴是尺度（可换算成周期），颜色深浅代表该时段的波动能量强弱。这就是所谓的小波功率谱，也叫尺度图（scalogram）。
如果信号里某个周期成分在某段时间增强，功率谱上就会在对应的“时间-周期”坐标处出现高亮斑块；如果某个周期信号随时间消失，斑块也随之褪色。这个特性直接回答了“某周期何时显著”的问题，而不是只告诉你“整体上有哪些周期”。
2.2 小波基函数的选型标准与参数说明
MATLAB的Wavelet Toolbox里提供的小波基类型很多，但做周期分析时，Morlet小波几乎是事实上的标准选择。原因是Morlet是一个复值小波，由高斯包络调制复指数波构成，它的实部和虚部之间有固定的90度相位差，可以同时提取振幅和相位信息。对于气候序列来说，相位信息常用来判断两个要素之间的滞后关系，比如降水领先气温还是气温领先降水，这只有复小波能做。
另一个理由是Morlet在时间域和频率域都有较好的局域性。它的波形接近真实的振荡信号，边缘效应相对可控，计算效率也高。相比之下，Mexican Hat小波虽然形状简单，但它没有复数形式，只能看能量，拿不到相位；DB系列小波是离散正交基，结构上适合信号重构和压缩，用在连续功率谱分析上反而不自然。
实际代码里最常用的是cmor小波，也就是复数Morlet。MATLAB中构造它的方式是：
% 构造复数Morlet小波基，带宽参数1.5，中心频率1.0 [psi, x] = cmorwavf(-8, 8, 1000, 1.5, 1.0);这段代码生成了一个离散化的Morlet小波波形。cmorwavf的五个参数分别是：支撑区间的左端点、右端点、采样点数、带宽参数和中心频率。带宽参数越大，小波在频率域的解析度越差，但在时间域的定位能力越强；中心频率决定小波振荡的快慢。对我处理过的年降水序列来说，带宽1.5、中心频率1.0这套组合是最常用的起点配置，它能在周期分辨率上得到一个相对均衡的观察效果。
2.3 连续小波变换系数计算的MATLAB实现
拿到小波基之后，实际计算用的是cwt函数。MATLAB从R2016b开始支持带频率参数的新版cwt语法，它不再返回传统意义上的尺度向量，而是直接返回周期或频率向量，这对做气候分析的人来说直观得多。
% 读取降雨距平序列，时间步长为1个月 load precipitation_anomaly.mat; Fs = 1; % 采样频率，单位：次/月 [wt, f] = cwt(precip_anomaly, Fs, 'voicesperoctave', 12);参数说明：precip_anomaly是降水距平序列，建议先做标准化处理，否则不同量纲的要素之间无法比较功率谱的绝对数值；Fs设为1，含义是每月一个采样点；voicesperoctave表示每个倍频程内划分的尺度数，我一般取12或者16，数值越大尺度越密，图像纵轴过渡越平滑，但计算量也会线性增加。返回的f是伪频率向量，单位与Fs一致。
这里有个关键点：用新版cwt算出来的功率谱是归一化后的相对值，不能直接拿来和文献里的绝对功率对比。如果需要绝对功率，应该用旧版语法加scal2frq做尺度-频率转换，或者手动实现完整的CWT积分。大多数周期分析的应用场景里，我们更关心的是相对能量分布和显著性，所以归一化版本反而更方便作图。
3. 小波功率谱的显著性检验：别把白噪声的随机起伏当成周期
3.1 为什么必须做显著性检验
小波功率谱是一个二维连续表面，功率值在时间和周期两个维度上都是平滑变化的。这就带来一个棘手的问题：即使输入信号完全是随机白噪声，功率谱表面也会出现若干高值斑块。如果不做检验，一眼看上去所有斑块都像“周期信号”，实际上其中相当一部分只是随机游走产生的统计涨落。
气候序列通常还带有红噪声特性——即低频段能量天然偏高，因为气候系统内部存在长期记忆效应。降水序列更是如此，干旱年的信号往往在之后若干年仍有余波，这种自相关性会让功率谱的低频端虚假抬升。Torrence和Compo在1998年提出的方法就是用红噪声作为背景谱，对每个尺度逐一检验功率值是否超出了95%置信线。
3.2 红噪声背景谱的构造
红噪声谱的强度由滞后1自相关系数alpha决定。计算alpha的公式是：alpha = (sum(x_t - mean)(x_{t+1} - mean)) / (sum(x_t - mean)^2)。在MATLAB中可以直接用autocorr函数取lag=1处的值。
% 计算序列滞后1自相关系数 [acf, lags] = autocorr(precip_anomaly, 'NumLags', 1); alpha = acf(2);这对降水距平序列来说基本都会是正值，体现在功率谱上就是低频（长周期）区域背景功率偏高。如果直接用白噪声做背景谱，会在低频端出现大量假阳性。
有了alpha之后，理论红噪声功率谱的表达式为：
P_k = (1 - alpha^2) / (1 + alpha^2 - 2alphacos(2pik/N))
其中k是频率序号，N是序列长度。这个公式给出的是一条平滑曲线，作为每个尺度上功率比较的基线。实际检验时，将每个时间点、每个尺度上的功率值除以对应尺度的背景功率，然后与卡方分布的95%分位数比较。
3.3 显著性检验的完整代码实现
这里给出一个已封装好的检验函数，直接输入cwt的功率谱矩阵和alpha值，返回显著性掩膜矩阵：
function sig_mask = wavelet_significance(wt_power, alpha, scale, dof) % 输入： % wt_power : cwt输出的功率谱矩阵（时间 x 频率） % alpha : 序列滞后1自相关系数 % scale : 与wt_power行数对应的尺度向量 % dof : 自由度，通常取2（复小波） % 输出： % sig_mask : 逻辑矩阵，1表示通过95%显著性检验 nScale = size(wt_power, 1); nTime = size(wt_power, 2); sig_mask = zeros(nScale, nTime); % 计算每个尺度对应的红噪声背景功率 P_background = (1 - alpha^2) ./ ... (1 + alpha^2 - 2*alpha*cos(2*pi*(1:nScale)/nScale)') ... .* ones(nScale, nTime); % 卡方检验：功率/背景功率 服从卡方分布 chi2_val = chi2inv(0.95, dof) / dof; % 逐点判断 ratio = wt_power ./ P_background; sig_mask(ratio > chi2_val) = 1; end逻辑说明：第一步把红噪声背景功率扩展到与功率谱相同的维度；第二步根据自由度查卡方分布的95%分位数，复小波的实部和虚部各贡献一个自由度，所以dof=2；第三步逐个像素点比较，功率超过背景功率的卡方阈值就判定为显著。
这里容易踩的一个坑是dof的取值。Torrence的经典方法在平滑后的功率谱上做检验时，dof取2是近似值，实际应该根据小波函数的尺度平滑特性修正。对cmor1.5-1.0来说，dof可以取2.0到3.5之间的值，取2是保守做法，会稍微高估显著性。如果发现检验结果显著区域过大，可以适当调高dof值重新计算。
4. 完整周期分析流程与绘图：从原始数据到显著性斑图
4.1 数据预处理的关键步骤
周期分析里最容易被忽略的就是预处理。很多人把原始降水序列直接丢进cwt函数，算出来的功率谱低频端被年际趋势严重污染，主周期信号反而被压得看不清。拿我当时处理某站点的月降水数据来说，序列长度约600个月，如果不做距平和标准化，功率谱上就会有一条巨大的能量带横跨所有周期，那是均值不为零导致的直流分量泄漏。
标准的做法分三步：第一步减去多年月均值，得到距平序列；第二步用序列标准差做归一化；第三步对距平序列做去趋势处理，一般减去线性趋势即可。对于有明显季节周期的要素，建议先按月标准化——即先算每个月的气候态均值，再用各月距平除以该月的标准差。这一步能消除季节方差差异对功率谱的干扰。
% 月降水数据预处理 data = readmatrix('precip_monthly.csv'); time = data(:, 1); precip = data(:, 2); % 重构为年-月矩阵，计算每月气候态 precip_mat = reshape(precip, 12, []); clim_mean = mean(precip_mat, 2); clim_std = std(precip_mat, 0, 2); % 逐月标准化：消除季节方差差异 precip_norm = (precip_mat - clim_mean) ./ clim_std; precip_anomaly = precip_norm(:)'; % 去线性趋势 precip_detrend = detrend(precip_anomaly, 'linear');这段代码的要点在于reshape成12行的矩阵后，每一行对应一个月份的多年序列，clim_mean和clim_std就是每个月份的气候均值和标准差。逐月标准化之后，功率谱上的季节信号不会被过度压制，也不会被异常放大。detrend这一步对趋势明显的序列必须做，否则最长的周期成分会始终占据主导。
4.2 功率谱计算与周期脊线提取
预处理完成之后，进入核心计算环节。这里我不建议直接调用MATLAB的cwt函数完事，而是配合自己写的小波基扫描函数，好处是尺度范围可以精确控制，低频端不会因为小波支撑不足而被截断。
% 自定义尺度范围，周期从2个月到序列长度的一半 dt = 1; % 月 s0 = 2*dt; % 最小尺度 sj = s0 * 2.^(0:0.125:log2(length(precip_detrend)*0.5/dt)); % 计算小波功率谱（基于FFT的快速算法） [wt_coef, f] = cwt(precip_detrend, 'morse', ... 'FrequencyLimits', [1/(length(precip_detrend)*0.8), 0.5]);这里我把尺度范围设成指数递增序列，周期覆盖从2个月到序列长度的一半左右。最小周期设2个月是因为对于月数据来说，周期短于2个月的信号本身不可分辨；最大周期设到序列长度的一半是为了避开边缘效应极严重的区域，那个区域的功率值基本不可信。
关于cwt返回的f向量和尺度向量之间，存在f = scal2frq(scale, wavelet, dt)的换算关系，其中wavelet是选取的小波名。如果后续要做显著性检验，建议自己用这个公式把尺度换算出来，再按尺度逐个做检验，免得频率和尺度顺序搞混。
4.3 显著性斑图与周期-时间剖面图绘制
绘图是周期分析里检验结果的最直观方式。一张标准的小波功率谱图由三个部分组成：功率谱等值线填色图、显著性区域黑色等高线、锥形影响区域（COI）阴影。这三个缺一不可，否则读图的人无从判断哪些区域的结果可信。
% 绘制小波功率谱 figure('Position', [100, 100, 1200, 400]); % 子图1：功率谱与显著性叠加 subplot(1, 2, 1); contourf(time, periods, wt_power, 40, 'LineStyle', 'none'); set(gca, 'YScale', 'log', 'YDir', 'reverse'); ylim([2, 64]); hold on; % 叠加显著性区域 contour(time, periods, sig_mask, [1, 1], 'LineColor', 'k', 'LineWidth', 1.5); % 绘制COI边界 patch(coi_x, coi_y, 'w', 'FaceAlpha', 0.6, 'EdgeColor', 'none'); colormap(jet); colorbar; xlabel('时间（年）'); ylabel('周期（月）'); title('小波功率谱（黑色等高线内为95%显著区域）'); % 子图2：全局平均功率谱 subplot(1, 2, 2); global_power = mean(wt_power, 2); plot(global_power, periods, 'b-', 'LineWidth', 1.2); hold on; % 叠加红噪声95%置信线 chi2_line = chi2inv(0.95, 2) / 2 * mean(P_background, 2); plot(chi2_line, periods, 'r--', 'LineWidth', 1); set(gca, 'YScale', 'log', 'YDir', 'reverse'); ylim([2, 64]); xlabel('全局平均功率'); ylabel('周期（月）'); title('全局小波功率谱与95%置信线');绘图参数里最值得留意的是ylim的选择。如果周期范围设置过大，边缘效应区域几乎占满整张图，中间的可信区域反而被压缩成一条窄带。我通常把显示范围控制在[2, N/4]的周期区间，宁可牺牲长周期段的展示，也要保证主周期附近的图形分辨率。COI的绘制用patch函数填充半透明白色，表示该区域受边界截断影响，功率值不可信，不必强行解读。
5. 小波分析避坑指南：六个血泪经验帮你保住结果的科学性
5.1 边界效应导致低频端整片伪显著
现象：功率谱图低频端出现一大片贯穿整个时间轴的显著区域，周期越长颜色越深，看起来像是存在一个超长周期信号。原因：小波变换在序列两端需要做延拓，MATLAB默认使用对称延拓。降水序列两端恰好处于趋势上升段时，对称延拓会人为制造出虚假的周期成分。解决：处理前先做去趋势，并检查序列两端是否为趋势的极值点。如果去趋势之后低频端伪显著区域仍然存在，需要对比不同延拓方式的结果，建议手动缩短分析时间范围，只保留序列中间的80%作为可信区间。
5.2 距离平没做标准化导致季节周期掩盖其他信号
现象：功率谱上1年周期处出现宽度异常的能量带，其他周期的信号全部被遮蔽。原因：降水量级在季节间差异悬殊，比如某地夏季月降水300毫米，冬季月降水不到20毫米，这种方差差异使冬季的弱周期信号在功率谱上完全不可见。解决：必须做逐月标准化，把每个月的序列按照该月自身的均值和标准差归一化。标准化做完之后，原来被淹没的年代际周期会浮出来。
5.3 显著性检验用错背景谱导致假阳性
现象：所有短周期全部通过95%检验，功率谱图几乎全黑。原因：背景谱误用了白噪声模型。气候序列普遍存在自相关性，白噪声背景谱在低频端远低于实际背景功率，导致短周期段的功率值几乎全部超过阈值。解决：用滞后1自相关系数构造红噪声背景谱后再检验。在代码里检查alpha值是否落在0到0.8之间，如果alpha为负值，说明序列经过过度差分，需要退回原始数据重新处理。
5.4 连续小波用了离散小波滤波器组
现象：功率谱呈块状分布，尺度轴不连续，纵轴看起来像阶梯。原因：MATLAB的cwt函数在高版本中存在两套实现——默认的滤波器组方法和经典积分小波变换方法。滤波器组方法适合做信号重构和滤波，但输出的功率谱尺度分辨率低。解决：在调用cwt时显式指定小波名称，例如cwt(x, 'morse')或cwt(x, 'cmor1.5-1')，这样会走积分小波变换路径，得到连续光滑的谱图。
5.5 小波基带宽参数选错导致周期脊线偏移
现象：图中主周期能量带的位置与已知主周期偏差20%以上。原因：cmor小波的带宽参数越大，频率分辨率越差，主峰会被展宽。带宽取3.0时，本应是12个月的周期信号可能被映射到10到15个月的宽带中。解决：对已知的ENSO周期信号做测试，分别用带宽1.0、1.5、2.0计算同一段序列，观察主峰位置的变化。我用带宽1.5得到的脊线与理论周期吻合度最好。
5.6 对标准化后的序列做累积距平导致周期失真
现象：功率谱显示出一个异常显著的长周期成分，对应周期接近序列全长。原因：对标准化序列再做累积距平，等效于对序列做一次积分滤波，这个操作会把低频端功率放大若干倍，制造出虚假的年代际周期。解决：周期分析输入序列必须是距平或标准化后的原始值，而不是累积值。累积序列只适合用于突变检验，不能送入cwt。
6. 多要素联合周期分析：交叉小波谱与相位滞后判断
6.1 两要素滞后关系的可视化方法
做了单要素的功率谱分析之后，下一步很自然地会问：降水和气温的周期信号之间是什么关系？在年尺度上谁领先谁？这时单要素功率谱已经没有能力回答，需要引入交叉小波变换（XWT）和小波相干（WTC）。
交叉小波谱的实质是两个序列的小波系数共轭相乘：W_xy = W_x * conj(W_y)。它的功率值反映两个序列在某个时段的共同振荡能量，相位角则反映两者的滞后关系。如果相位角集中在零附近，说明两者同步；如果相位角稳定在正区间，说明序列x领先序列y约四分之一个波长。
% 交叉小波谱计算 [wt_x, f] = cwt(precip_detrend, 'morse', 'VoicesPerOctave', 12); [wt_y, f] = cwt(temp_detrend, 'morse', 'VoicesPerOctave', 12); % 交叉谱：共轭相乘 xwt = wt_x .* conj(wt_y); xwt_power = abs(xwt) .^ 2; % 相位角计算 phase_angle = angle(xwt);逻辑说明：xwt_power里的每个像素点数值表示该时段该周期下两个序列的联合振荡强度；phase_angle取的是复数辐角，范围在[-pi, pi]之间。画联合功率谱时用xwt_power填充色块，再在色块上用箭头场叠加相位信息。箭头方向向右表示两序列同相，向上表示第二个序列领先90度，箭头向左表示反相。
绘制时需要注意，箭头应该在显著性区域内才绘制，否则满屏箭头完全没法读。显著性的判断方法与单要素类似，但背景谱要改用两个红噪声序列的联合置信水平。MATLAB中可以用bootstrapping方法生成随机相位扰动序列，反复计算交叉谱得到经验分布，然后取95%分位数作为阈值。
% 基于bootstrap的交叉谱显著性 nBoot = 200; xwt_sig = zeros(size(xwt_power)); parfor b = 1:nBoot % 随机打乱其中一个序列的相位 phase_rand = exp(1i * 2 * pi * rand(size(precip_detrend))); x_rand = real(ifft(fft(precip_detrend) .* phase_rand)); [wt_xr, ~] = cwt(x_rand, 'morse', 'VoicesPerOctave', 12); xwt_rand = abs(wt_xr .* conj(wt_y)) .^ 2; xwt_sig = xwt_sig + (xwt_power > xwt_rand); end xwt_sig = xwt_sig / nBoot;bootstrap参数说明：nBoot取200次，每次随机打乱序列x的傅里叶相位后再与y做交叉谱，得到的随机交叉谱集合构成了零假设分布。如果真实交叉谱功率高于95%的随机结果，判定为显著。parfor是并行循环，如果有并行计算工具箱可以加速，没有的话改成普通for循环即可，只是多等一会儿。
6.2 验证与结果解读的工程习惯
整套流程做完之后，我不会直接拿去写报告。我会先把结果跟已知的物理规律做一次对照，比如某区域的降水主周期如果是准两年振荡，功率谱上应该在24-32个月区间有一个清晰且连续的显著带，相位分析也应该能对得上大气环流的已知滞后关系。对照不上的部分，优先怀疑数据预处理环节是不是某一步出了问题。
从那以后我处理每一组气候序列，都强制走一遍“预处理-距平标准化-红色噪声检验-显著性斑图-COI检查”的完整流程，分量级的参数记录也单独存一份。这套流程帮我筛掉过至少三次因为边界效应导致的伪周期误判。
对于你手里的这套MATLAB小波功率谱分析资源，建议拿到手后先用自带示例数据跑通一遍，输出功率谱图对比一下形态是否和预期一致，再把代码里的尺度范围、显著性阈值、绘图参数替换成你自己的数据配置。前40分钟看起来像在调试，实际上是在给你的数据寻找一个可靠的周期分析基线。希望帮到你。
本文还有配套的精品资源，点击获取
