MATLAB小波交叉功率谱分析源代码实现

一、核心代码框架

基于MATLAB小波工具箱（Wavelet Toolbox），实现两个时间序列的小波交叉功率谱分析，包含数据预处理、参数设置、交叉谱计算及可视化。

%% 1. 数据加载与预处理
% 加载两个时间序列数据（示例：温度与降水数据）
load('climate_data.mat'); % 假设数据包含变量temp和precip
data1 = temp; % 时间序列1（如温度）
data2 = precip; % 时间序列2（如降水）

% 数据归一化（消除量纲影响）
data1 = (data1 - mean(data1)) / std(data1);
data2 = (data2 - mean(data2)) / std(data2);

% 确保数据长度一致
n = min(length(data1), length(data2));
data1 = data1(1:n);
data2 = data2(1:n);

%% 2. 设置小波变换参数
dt = 1.0; % 时间采样间隔（年）
pad = 1; % 填充零（推荐）
dj = 0.1; % 尺度分辨率（0.1表示更精细的尺度划分）
s0 = 2*dt; % 最小尺度（对应最低频率）
j1 = 7/dj; % 最大尺度数（根据数据长度调整）
mother = 'morlet'; % 小波基函数（Morlet适用于周期分析）

%% 3. 执行连续小波变换
[c1, period, scale, coi1] = wavelet(data1, dt, pad, dj, s0, j1, mother);
[c2, ~, ~, coi2] = wavelet(data2, dt, pad, dj, s0, j1, mother);

% 计算交叉功率谱（复数形式）
cross_wave = c1 .* conj(c2); % 交叉小波系数
cross_power = abs(cross_wave).^2; % 交叉功率谱密度

%% 4. 计算显著性水平
% 自功率谱显著性水平（用于交叉谱显著性参考）
[signif1, fft_theor1] = wave_signif(1.0, dt, scale, 0, 0.72, -1, -1, mother);
[signif2, fft_theor2] = wave_signif(1.0, dt, scale, 0, 0.72, -1, -1, mother);

% 交叉谱显著性（假设独立噪声，需根据实际调整）
signif_cross = sqrt(signif1' .* signif2'); % 假设噪声独立
sig95_cross = (signif_cross') * (ones(1, n)); % 扩展至全矩阵

%% 5. 可视化结果
figure;

% 子图1：时间序列原始数据
subplot(3,1,1);
plot(data1, 'r', 'LineWidth', 1.5); hold on;
plot(data2, 'b', 'LineWidth', 1.5);
xlabel('Time (year)');
ylabel('Normalized Value');
legend('Temperature', 'Precipitation');
title('Original Time Series');

% 子图2：交叉小波功率谱
subplot(3,1,2);
contourf(period, log2(scale), log2(cross_power), 20);
hold on;
contour(period, log2(scale), sig95_cross, [1,1], 'k'); % 显著性轮廓
plot(period, log2(coi1), 'k', 'LineWidth', 1.5); % 影响锥
title('Cross-Wavelet Power Spectrum');
xlabel('Time (year)');
ylabel('Scale (year)');
colorbar;
set(gca, 'YDir', 'reverse');

% 子图3：全局交叉谱密度
subplot(3,1,3);
plot(global_ws_cross, log2(scale));
hold on;
plot(global_signif_cross, log2(scale), '--r');
xlabel('Cross-Wavelet Power');
ylabel('Scale (year)');
title('Global Cross-Wavelet Spectrum');
colorbar;
set(gca, 'YDir', 'reverse');

%% 6. 保存结果
save('cross_wavelet_results.mat', 'cross_power', 'period', 'scale', 'coi1', 'global_ws_cross');

---

二、关键参数说明

1. 小波基选择

   • mother = 'morlet'：Morlet小波（时频局部化最佳，推荐用于周期分析）。

   • 其他选项：'dgauss'（高斯导数）、'paul'（Paul小波）。

2. 尺度参数

   • dj：尺度分辨率，值越小频率分辨率越高（建议0.1-0.25）。

   • s0：最小尺度需满足 s0 ≥ 2*dt，避免频谱泄漏。

3. 显著性计算

• wave_signif函数基于蒙特卡洛模拟生成理论谱，需确保数据平稳性。

• 交叉谱显著性可结合自谱显著性计算（如几何平均或独立噪声假设）。

---

三、应用场景与优化建议

1. 气候分析

   • 案例：分析气温与降水的跨尺度相关性，识别ENSO事件中的同步性。

   • 优化：调整dj和j1以平衡频率分辨率与计算效率。

2. 地球物理信号处理

   • 案例：地震波与地磁扰动的小波交叉谱分析，提取传播时延。

   • 优化：使用phase = angle(cross_wave)提取相位差，结合线性拟合计算时间偏移。

3. 金融时间序列

• 案例：股票指数与大宗商品价格的关联性分析。

• 优化：采用'paul'小波增强高频成分捕捉短期波动。

---

四、常见问题与解决方案

1. 数据长度不足

   • 问题：小波变换要求数据长度为2的幂次。

   • 解决：使用pad = 1自动填充零，或截断数据至最近幂次。

2. 边界效应

   • 问题：小波变换边缘区域受噪声干扰。

   • 解决：通过coi（Cone of Influence）标记影响区域，忽略边界数据。

3. 频率混叠

• 问题：高频信号被错误解析为低频。

• 解决：增加采样率或使用'dgauss'小波抑制混叠。

参考代码 matlab小波交叉功率谱分析源代码 www.youwenfan.com/contentcnq/55082.html

五、扩展功能代码

1. 相位时延计算

   phase = angle(cross_wave); % 提取相位谱
   [m, n] = size(phase);
   time_delay = zeros(m, 1);
   for i = 1:m
       % 线性拟合相位随时间变化（频率固定时）
       p = polyfit(time, phase(:,i), 1);
       time_delay(i) = -p(1) * (2*pi)/scale(i); % 时延公式
   end
   

2. 动态显著性阈值

   % 基于局部噪声估计的显著性（替代全局阈值）
   noise_level = median(abs(cross_wave(:,2:end-1)))/(0.6745);
   sig95_dynamic = 1.96*repmat(noise_level, 1, n);