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

一、核心代码框架

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

%% 1. 数据加载与预处理% 加载两个时间序列数据（示例：温度与降水数据）load('climate_data.mat');% 假设数据包含变量temp和precipdata1=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/contentcsq/55082.html

五、扩展功能代码

1. 相位时延计算

   phase=angle(cross_wave);% 提取相位谱[m,n]=size(phase);time_delay=zeros(m,1);fori=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);