基于Matlab的多维小波相干（MWC)、偏小波（PWC)与全局相干性

码农世界 2024-05-19 前端 239 次浏览 0个评论

文章目录

前言
偏小波相干（PWC）
多小波相干（MWC）
全局小波相干（Global Coherence）

部分代码
出图效果

参考文献

前言

在前面我们提到了两个变量的交叉小波分析，一个自变量和一个因变量，但在实际中，某一变量往往受多个其它变量的影响。我们往往关注在不同尺度与时间/空间位置上，因变量如何受多个变量的影响。虽然一些多变量方法例如多变量谱一致性（multiple spectral coherence）、多变量经验模态分解（multivariate empirical mode decomposition）可以分析多变量之间在不同尺度上的相关性，但是这些方法都假设相关性特征不随时间/空间位置变化。另外，与传统相关性方法一样，当其它变量与自变量相关时，其它变量的存在可能会导致两变量小波相干性结果具有误导性，从而不能揭示内在过程的真正机理。

在信号处理和频域分析中，小波变换是一种有效的工具，能够揭示时间序列数据的频率特征和变化规律。小波相干（Wavelet Coherence，WTC）及其衍生的偏小波相干（Partial Wavelet Coherence，PWC）、多小波相干（Multiple Wavelet Coherence，MWC）以及全局小波相干（Global Coherence）等概念，提供了更深入的频域相关性分析，有助于理解信号之间的关联程度及其时间-频率特性。本文将对这些概念进行系统解析，并介绍在MATLAB中如何实现这些分析。

偏小波相干（PWC）

偏小波相干是对小波相干的一种改进，它在分析两个信号之间的相干性时，考虑了第三个信号（或多个信号）的影响。PWC能够剔除第三个信号对两个信号之间关系的干扰，更加准确地反映出两个信号之间的真实相关性。

在MATLAB中，计算PWC的方法与WTC类似，只是需要在计算小波相干度时，加入第三个信号的信息。以下是一个简单的示例代码：

多小波相干（MWC）

多小波相干是对多个信号之间相干性的分析方法，它考虑了多个信号之间的交互关系，能够揭示复杂系统中信号之间的频域耦合情况。MWC的计算与PWC类似，只是需要将多个信号的信息纳入到小波变换和相干度计算中。

全局小波相干（Global Coherence）

全局小波相干是对整个信号集合的相干性进行综合分析的方法。它不仅考虑了信号之间的相互关系，还考虑了整体信号集合的频域特征和共同振荡模式。全局小波相干性分析能够帮助我们理解复杂系统中信号之间的整体协调性和耦合程度。

部分代码

function varargout=pwc(X,varargin)
%---------parse function arguments------------------------------------------
 
[row,col]=size(X);
[y1,dt1]=formatts(X(:,1));
mm=y1(1,1);
nn=y1(end,1);
 
[y2,dt2]=formatts(X(:,2));
mm2=y2(1,1);
nn2=y2(end,1);
 
for i=3:col;
[x,dtx]=formatts(X(:,i));
 
if (dt1~=dtx |dt2~=dtx)
    error('timestep must be equal between time series');
end
 
mm1=x(1,1);
nn1=x(end,1);
 
 
mm=max([mm,mm1,mm2]);
nn=min([nn,nn1,nn2]);
 
x1(:,(i-2))=x(:,1);
x2(:,(i-2))=x(:,2);
 
end
 
t=(mm:dt1:nn)';
 
 
%common time period
if length(t)<4
    error('The three time series must overlap.');
end 
  
n=length(t);
 
%----------default arguments for the wavelet transform-----------
Args=struct('Pad',1,...      % pad the time series with zeroes (recommended)
            'Dj',1/12, ...    % this will do 12 sub-octaves per octave
            'S0',1*dt1,...    % this says start at a scale of 2 years
            'J1',[],...
            'Mother','Morlet', ...
            'MaxScale',[],...   %a more simple way to specify J1
            'MakeFigure',(nargout==0),...
            'MonteCarloCount',300,...
            'BlackandWhite',0,...
            'AR1','auto',...
            'ArrowDensity',[30 30],...
            'ArrowSize',1,...
            'ArrowHeadSize',1);
 
Args=parseArgs(varargin,Args,{'BlackandWhite'});
 
if isempty(Args.J1)
    if isempty(Args.MaxScale)
        Args.MaxScale=(n*.17)*2*dt1; %auto maxscale;
    end
    Args.J1=round(log2(Args.MaxScale/Args.S0)/Args.Dj);
end
 
ad=mean(Args.ArrowDensity);
Args.ArrowSize=Args.ArrowSize*30*.03/ad;
%Args.ArrowHeadSize=Args.ArrowHeadSize*Args.ArrowSize*220;
Args.ArrowHeadSize=Args.ArrowHeadSize*120/ad;
 
if ~strcmpi(Args.Mother,'morlet')
    warning('MWC:InappropriateSmoothingOperator','Smoothing operator is designed for morlet wavelet.');
end
 
if strcmpi(Args.AR1,'auto')
      for i=1:col
arc(i)= ar1nv(X(:,i));
       end
Args.AR1=arc
    if any(isnan(Args.AR1))
        error('Automatic AR1 estimation failed. Specify it manually (use arcov or arburg).');
    end
end
 
%-----------:::::::::::::--------- ANALYZE ----------::::::::::::------------
 
%Calculate and smooth wavelet spectrum Y1, Y2, and X
  
[Y1,period,scale,coiy1] = wavelet(y1(:,2),dt1,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);
sinv=1./(scale');
smY1=smoothwavelet(sinv(:,ones(1,n)).*(abs(Y1).^2),dt1,period,Args.Dj,scale);
dte=dt1*.01; 
idx=find((y1(:,1)>=(t(1)-dte))&(y1(:,1)<=(t(end)+dte)));
Y1=Y1(:,idx);
smY1=smY1(:,idx);
coiy1=coiy1(idx);
coi1=coiy1;
 
 
[Y2,period,scale,coiy2] = wavelet(y2(:,2),dt1,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);
sinv=1./(scale');
smY2=smoothwavelet(sinv(:,ones(1,n)).*(abs(Y2).^2),dt1,period,Args.Dj,scale);
dte=dt2*.01; 
idx=find((y2(:,1)>=(t(1)-dte))&(y2(:,1)<=(t(end)+dte)));
Y2=Y2(:,idx);
smY2=smY2(:,idx);
coiy2=coiy2(idx);
coi2=coiy2;
 
 
 
for  i=3:col
 [XS,period,scale,coix] = wavelet(x2(:,(i-2)),dt1,Args.Pad,Args.Dj,Args.S0,Args.J1,Args.Mother);
 
idx=find((x1(:,(i-2))>=(t(1))-dte)&(x1(:,(i-2))<=(t(end)+dte)));
XS=XS(:,idx);
coix=coix(idx);
 
XS1(:,:,(i-2))=XS;
coi0=min(coi1,coi2);
coi=min(coi0,coix);
end
 
% -------- Calculate Cross Wavelet Spectra----------------------------
 
% -- between dependent variable (or independent variables) and excluding
% factors------
 
for i=1:(col-2)
Wy1x=Y1.*conj(XS1(:,:,i));
sWy1x=smoothwavelet(sinv(:,ones(1,n)).*Wy1x,dt1,period,Args.Dj,scale);
sWy1x1(:,:,i)=sWy1x;
 
Wy2x=Y2.*conj(XS1(:,:,i));
sWy2x=smoothwavelet(sinv(:,ones(1,n)).*Wy2x,dt1,period,Args.Dj,scale);
sWy2x1(:,:,i)=sWy2x;
end
 
% ---- between dependent variable and independent variables------
Wy1y2=Y1.*conj(Y2);
sWy1y2=smoothwavelet(sinv(:,ones(1,n)).*Wy1y2,dt1,period,Args.Dj,scale);
 
% ----between excluding variables and excluding variables------
for i=1:(col-2);
for j=1:(col-2);
Wxx=XS1(:,:,i).*conj(XS1(:,:,j));
sWxx=smoothwavelet(sinv(:,ones(1,n)).*Wxx,dt1,period,Args.Dj,scale);
sWxx1(:,:,i,j)=sWxx;
end
end
 
% --------------- Partial wavelet coherence ---------------------------------
% calculate the partial wavelet coherence 
m=length(scale);
Rsq=zeros(m,n);
Wy1y2x=zeros(m,n);
nofe=col-2; %number of excluding factors
 
if nofe==1;    
% ------Partial wavelet coherence for one excluding variable ------------------
%%%%%%%%% |1-R^2 y,x穁 (s,t)|^2   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ii=1:m;
    for jj=1:n;
        sWxx1_i(ii,jj)=inv(sWxx1(ii,jj)); % inversed sWxx1
    end
end
 
sWy1x1_sWxx1=bsxfun(@times,sWy1x1,sWxx1_i); % sWy1x1*sWxx1_i
sWy2x1_c=conj(sWy2x1); %conjugate of smoothed cross-wavelet power spectra
% between independent variable and excluding variables
sWy1x1_sWxx1_sWy2x1=bsxfun(@times,sWy1x1_sWxx1,sWy2x1_c); 
% sWy1x1*sWxx1_i*sWy2x1_c
R2y1y2x=sWy1x1_sWxx1_sWy2x1./sWy1y2; %(Ry,x.Z)^2
abs_one_minus_R2y1y2x=(abs(1-R2y1y2x)).^2; %squared absolute of 1-(Ry,x.Z)^2
%%%%%%%%%%%%%%%% R^2 y,x (s,t)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
R2y1y2=(sWy1y2.*conj(sWy1y2))./(smY1.*smY2); %(Ry,x)^2
%%%%%%%%%%%%%%%%%%% 1-R^2 y,Z(s,t)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sWy1x1_c=conj(sWy1x1);
sWy1x1_sWxx1_sWy1x1=bsxfun(@times,sWy1x1_sWxx1,sWy1x1_c);
R2y1x=sWy1x1_sWxx1_sWy1x1./smY1; %(Ry,Z)^2
one_minus_R2y1x=(1-R2y1x);
 
%%%%%%%%%%%%%%%%%%%%%  1-R^2 x,Z(s,t) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sWy2x1_sWxx1=bsxfun(@times,sWy2x1,sWxx1_i);
sWy2x1_sWxx1_sWy2x1=bsxfun(@times,sWy2x1_sWxx1,sWy2x1_c);
R2y2x=sWy2x1_sWxx1_sWy2x1./smY2; %(Rx,Z)^2
one_minus_R2y2x=(1-R2y2x);
%%%%%%%%%%%%%%%%%%  R^2 y,x穁 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Rsq_num=roundn(bsxfun(@times,abs_one_minus_R2y1y2x,R2y1y2),-16); 
%numerator part of the equation for squared pwc, if the value is less than 
%10^(-16), 0 is assigned
Rsq_den=real(bsxfun(@times,one_minus_R2y1x,one_minus_R2y2x));  
%denominator part of the equation for squared pwc
Rsq= bsxfun(@rdivide,Rsq_num,Rsq_den); %squared partial wavelet coherence
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

出图效果

土壤湿度的小波相干性和三因子组合被认为是解释每个研究点土壤湿度变化的最具代表性的。上排子图为表土层(0-10 cm)的分析，第二行对应第二层土层(10-40 cm)，第三行对应第三层土层(40-100 cm)，下行为第四层土层(100-200 cm)。淡色区域是影响锥(COI)，其中边缘效应可能会扭曲结果。粗轮廓对红色噪声的显著性水平为5%，颜色表示一致性的强度。

欢迎扫描关注我的公众号，博士期间日常分享日常，获取更多科研代码和前沿论文资讯等相关内容

若需要以上小波分析的代码，请关注公众号回复“小波分析”关键词

参考文献

Hu, W., and B.C. Si (2021), Technical Note: Improved partial wavelet coherency for understanding scale-specific and localized bivariate relationships in geosciences, Hydrol. Earth Syst. Sci.,25,321-331.

Gu, X., Jamshidi, S., Sun, H., & Niyogi, D. (2021). Identifying multivariate controls of soil moisture variations using multiple wavelet coherence in the U.S. Midwest. Journal of Hydrology, 602, 126755

Zhao, R., Biswas, A., Zhou, Y., Zhou, Y., Shi, Z., & Li, H. (2018). Identifying localized and scale-specific multivariate controls of soil organic matter variations using multiple wavelet coherence. Science of the total environment, 643, 548-558

转载请注明来自码农世界，本文标题：《基于Matlab的多维小波相干（MWC)、偏小波（PWC)与全局相干性》