1. 文件结构

WaveletPeriod/
├── main_wavelet_period.m   % 一键运行
├── wavelet_power_spectrum.m % 小波功率谱 + 显著性
├── period_peak_detect.m     % 自动周期峰值
├── plot_wavelet_results.m   % 时频图 + 周期图
└── example/└── temp.csv             % 示例月平均温度

2. 核心函数

① 小波功率谱（wavelet_power_spectrum.m）

function [W, period, scale, sig95] = wavelet_power_spectrum(x, dt, mother)
% 输入：
%   x     : 1×N 时间序列
%   dt    : 采样间隔（月/日/年）
%   mother: 'Morlet'（默认）
% 输出：
%   W     : 小波功率谱（N×Nscale）
%   period: 周期向量（与 W 行对应）
%   sig95 : 95% 显著性水平if nargin<3, mother = 'Morlet'; end
% 1. 选择母小波
[wave, scale] = morlet_wavelet(length(x), dt);
% 2. 连续小波变换
W = cwt(x, scale, wave);
% 3. 显著性检验（红噪声背景）
sig95 = red_noise_significance(x, dt, scale);
end

② Morlet 母小波生成（morlet_wavelet.m）

function [wave, scale] = morlet_wavelet(N, dt)
s0   = 2*dt;           % 最小尺度
dj   = 0.125;          % 尺度步长
J    = floor(log2(N)*dj^-1);
scale = s0 * 2.^((0:J)*dj);
wave  = { @(s,t) pi^(-0.25) .* exp(1j*5*t./s) .* exp(-0.5*(t./s).^2) };
end

③ 显著性检验（红噪声背景）

function sig95 = red_noise_significance(x, dt, scale)
% Torrence & Compo 红噪声显著性
alpha = ar1_coeff(x);                % 估计 AR(1) 系数
varx  = var(x);
sig95 = varx * (1-alpha^2) ./ (1 - 2*alpha*cos(2*pi/(scale/dt)) + alpha^2);
sig95 = sig95 * 0.95;                % 95% 分位
end

3. 运行（main_wavelet_period.m）

clear; clc; close all; addpath('.');%% 1. 读入时间序列（月平均温度）
data = readmatrix('example/temp.csv');  % 两列：日期, 温度
t    = data(:,1);   x = data(:,2);
dt   = 1;           % 月间隔%% 2. 小波功率谱
[W, period, scale, sig95] = wavelet_power_spectrum(x, dt);%% 3. 周期峰值检测
[pkPer, pkVal] = period_peak_detect(W, period);%% 4. 可视化
plot_wavelet_results(t, x, W, period, sig95, pkPer, pkVal);