www.gusucode.com > wavelet工具箱matlab源码程序 > wavelet/wavelet/modwtxcorr.m

    function [xcseq,xcseqci,lags] = modwtxcorr(w1,w2,varargin)
%MODWTXCORR Wavelet cross-correlation sequence estimates using the MODWT.
%   XCSEQ = MODWTXCORR(W1,W2) returns the wavelet cross-correlation
%   sequence estimates for the MODWT transforms in W1 and W2. W1 and W2
%   must be the same size and obtained with the same wavelet. By default,
%   MODWTXCORR assumes the 'sym4' wavelet was used. XCSEQ is a cell array
%   of (2NJ-1)-by-1 vectors where NJ is the number of nonboundary
%   coefficients by level. The length of XCSEQ is equal to the number of
%   levels in W1 and W2 with nonboundary coefficients. This level is the
%   minimum of size(W1,1) and floor(log2(N/(L-1)+1)) where N is the length
%   of the data and L is the filter length. The elements in each cell of
%   XCSEQ correspond to cross-correlation sequence estimates at lags from
%   -(NJ-1) to (NJ-1). The 0-th lag element is located at the index
%   floor((2*NJ-1)/2)+1. If there are sufficient nonboundary coefficients
%   at the final level, MODWTXCORR returns the scaling cross-correlation
%   sequence estimate in the final cell of XCSEQ.
%
%   XCSEQ = MODWTXCORR(W1,W2,WAV) uses the wavelet WAV to determine the
%   number of boundary coefficients by level. WAV can be a string
%   corresponding to a valid wavelet or a positive even integer indicating
%   the length of the wavelet filter. If WAV is unspecified or specified 
%   as empty, the default 'sym4' wavelet is used. The wavelet filter length
%   must match the length used in the MODWT of the inputs.
%
%   [XCSEQ,XCSEQCI] = MODWTXCORR(...) returns 95% confidence intervals for
%   the cross-correlation sequence estimates in XCSEQ. XCSEQCI is a cell
%   array of (2NJ-1)-by-2 matrices. The first column of the k-th element of
%   XCSEQCI contains the lower 95% confidence bound for the
%   cross-correlation coefficient at each lag. The second column contains
%   the upper 95% confidence bound. Confidence intervals are computed using
%   Fisher's Z-transformation. The N in the standard error of Fisher's Z
%   statistics, sqrt(N-3), is determined by the equivalent number of
%   coefficients in the critically sampled DWT, floor(size(w1,2)/2^LEV),
%   where LEV is the level of the wavelet transform. MODWTXCORR returns
%   NaNs for the confidence bounds when N-3 is less than or equal to zero.
%
%   [...] = MODWTXCORR(W1,W2,WAV,ConfLevel) uses ConfLevel for the
%   coverage probability of the confidence interval. ConfLevel is a real
%   scalar strictly greater than 0 and less than 1. If ConfLevel is
%   unspecified or specified as empty, the coverage probability defaults 
%   to 0.95. 
%
%   [XCSEQ,XCSEQCI,LAGS] = MODWTXCORR(...) returns the lags for the wavelet
%   cross-correlation sequence estimates. LAGS is a cell array of column
%   vectors the same length as XCSEQ. The elements of LAGS range from
%   -(NJ-1) to (NJ-1) where NJ is the number of nonboundary coefficients at
%   level J.
%
%   [...] = MODWTXCORR(...,'reflection') reduces the number of wavelet and
%   scaling coefficients at each scale by 1/2. Use this option when the
%   MODWT of W1 and W2 were obtained using the 'reflection' boundary
%   condition. You must enter the entire string 'reflection'. If you added
%   a wavelet named 'reflection' using the wavelet manager, you must rename
%   that wavelet prior to using this option. 'reflection' may be placed in
%   any position in the input argument list after the input transforms W1
%   and W2. Specifying the 'reflection' option in MODWTXCORR is equivalent
%   to first obtaining the MODWT of W1 and W2 with 'periodic' boundary
%   handling and then computing the wavelet cross-correlation sequence
%   estimates.
%
%   %Example 1:   
%   %   Plot the cross-correlation sequence estimate for the Southern
%   %   Oscillation Index and Truk Island pressure data for scale 2^5.
%   load soi
%   load truk
%   wsoi = modwt(soi);
%   wtruk = modwt(truk);
%   [xcseq,xcseqci,lags] = modwtxcorr(wsoi,wtruk);
%   plot(lags{5},xcseq{5},'linewidth',2)
%   hold on
%   plot(lags{5},xcseqci{5},'r--')
%   set(gca,'xlim',[-120 120]);
%   xlabel('Lag (Days)'); 
%   grid 
%   title({'Cross-Correlation Sequence Level 5'; 'Scale: 32-64 Days'});
%   hold off
%   
%   %Example 2:
%   %   Plot wavelet correlation for two 5-Hz sine wave signals with
%   %   additive noise.
%   dt = 0.01;
%   t = 0:dt:6;
%   x = cos(2*pi*5*t)+1.5*randn(size(t));
%   y = cos(2*pi*5*t-pi)+2*randn(size(t));
%   wx = modwt(x,'fk14',5);
%   wy = modwt(y,'fk14',5);
%   modwtcorr(wx,wy,'fk14')
%
%   % For the significant correlation at scale 4, plot the wavelet
%   % cross-correlation sequence
%   J = 4;
%   [xcseq,xcseqci,lags] = modwtxcorr(wx,wy,'fk14');
%   zerolag = floor(numel(lags{J})/2)+1;
%   lagidx = zerolag-30:zerolag+30;
%   plot(lags{J}(lagidx).*dt,xcseq{J}(lagidx));
%   hold on;
%   grid 
%   plot(lags{J}(lagidx).*dt,xcseqci{J}(lagidx,:),'r--');
%   xlabel('Lag (Seconds)');
%   title('Scale: 0.32-0.16 Seconds');
%   
%   See also MODWTCORR, MODWTVAR, MODWTMRA, MODWT, IMODWT


% MODWTXCORR accepts between 2 and 5 inputs
narginchk(2,5);

% Ensure that the input has at least two rows
if (isrow(w1) || iscolumn(w1))
    error(message('Wavelet:modwt:InvalidCFSSize'));
end

% Ensure that the input matrices are the same size
% Ensure that the input matrices are the same size
if (numel(w1) ~= numel(w2))
    error(message('Wavelet:modwt:CFSMatrixSize'));
end

%Validate that the inputs are double-precision, real-valued with no
% NaNs or Infs
validateattributes(w1,{'double'},{'real','nonnan','finite'});
validateattributes(w2,{'double'},{'real','nonnan','finite'});

% Parse inputs
params = parseinputs(varargin{:});
ConfLevel = params.ConfLevel;
filtlen = params.L;
boundary = params.boundary;

% Level of the MODWT
level = size(w1,1)-1;

% Extract scaling coefficients
V1 = w1(end,:);
V2 = w2(end,:);
scalingcorr = false;

% Keep just the wavelet coefficients
w1 = w1(1:end-1,:);
w2 = w2(1:end-1,:);

% If the boundary is specified as 'reflection', remove the last N/2
% coefficients
if strcmpi(boundary,'reflection')
    if mod(size(w1,2),2)
        error(message('Wavelet:modwt:EvenLengthInput'));
    end
    N = size(w1,2)/2;    
else
    N = size(w1,2);
end

% For an unbiased estimate, keep only scales with nonboundary coefficients
Jmax = floor(log2((N-1)/(filtlen-1)+1));
if (Jmax<1)
     error(message('Wavelet:modwt:ZeroNonBoundaryCFS'));
 end
Jmax = min(Jmax,level);
w1 = w1(1:Jmax,1:N);
w2 = w2(1:Jmax,1:N);
V1 = V1(1:N);
V2 = V2(1:N);
 
 if (Jmax-level==0)
     scalingcorr = true;
 end
 
 % Remove boundary coefficients
cfs1 = removemodwtboundarycoeffs(w1,V1,N,Jmax,filtlen,scalingcorr);
[cfs2,MJ] = removemodwtboundarycoeffs(w2,V2,N,Jmax,filtlen,scalingcorr);

J = 1:Jmax;

if scalingcorr
    % If we are returning a correlation estimate of the scaling
    % coefficients
   J = [J Jmax];
end

% Adjust confidence level for symmetric distribution
ConfLevelComplement = 1-ConfLevel;
ConfLevelComplement = ConfLevelComplement/2;
ConfLevel = 1-ConfLevelComplement;

% Get critical value
qnorm = -sqrt(2)*erfcinv(2*ConfLevel);

xcseq = cell(numel(J),1);
xcseqci = cell(numel(J),1);
lags = cell(numel(J),1);


% The N in the Fisher Z-transformation for the cross-correlation
% coefficients comes from the critically-sampled DWT
NDWT = floor(size(w1,2) ./ 2.^J);

% calculate the estimator of the wavelet autocorrelation or
% cross-correlation sequence
for jj = 1:numel(J)
    cfsNoNaN1 = cfs1(jj,~isnan(cfs1(jj,:)));
    cfsNoNaN2 = cfs2(jj,~isnan(cfs2(jj,:)));
    [ccs,tau] = modwtCCS(cfsNoNaN1,cfsNoNaN2,MJ(jj));
    lowerci = real(tanh(atanh(ccs) - qnorm ./ sqrt(NDWT(jj)-3)));
    upperci = real(tanh(atanh(ccs)+qnorm./sqrt(NDWT(jj)-3)));
    xcseq{jj} = ccs';
    xcseqci{jj} = [lowerci' upperci'];
    lags{jj} = tau;
end

% For any NDWT<=3 replace confidence intervals with NaNs
if any(NDWT<=3)
 [xcseqci{NDWT<=3}] = deal(NaN(1,2));
end



%-----------------------------------------------------------------------
function [wccs,tau] = modwtCCS(cfs1,cfs2,MJ)
% Cross-correlation sequence estimation
N = size(cfs1,2);
fftpad = 2^nextpow2(2*N-1);
% floor the fftpad for the edge case that MJ = 1
zerolag = floor(fftpad/2)+1;
idxbegin = zerolag-(MJ-1);
idxend = zerolag+(MJ-1);
tau = (-(MJ-1):MJ-1)';
SSX = sum(abs(cfs1).^2);
SSY = sum(abs(cfs2).^2);
scalefactor = sqrt(SSX*SSY);
wccsDFT = fft(cfs1,fftpad).*conj(fft(cfs2,fftpad));
wccs = ifftshift(ifft(wccsDFT));
wccs = wccs(idxbegin:idxend);
wccs = wccs./scalefactor;

%-----------------------------------------------------------------------
function params = parseinputs(varargin)
% Parse inputs
         
params.boundary = 'periodic';
params.L = 8;
params.ConfLevel = 0.95;
    
tfboundary = strcmpi(varargin,'reflection');
    if any(tfboundary)
        params.boundary = 'reflection';
        varargin(tfboundary>0) = [];
    end
    
    if isempty(varargin)
        return;
    end
    
Len = length(varargin);
% The wavelet must be the first input argument in varargin
wavlen = varargin{1};
% Handle cases where the wavelet is a string, or a scalar, or a vector, or
% empty

    if ischar(wavlen)
        [~,~,Lo,~] = wfilters(wavlen);
        params.L = length(Lo);
    elseif isscalar(wavlen)
        params.L = wavlen;
    elseif (isrow(wavlen) || iscolumn(wavlen))
        params.L = length(wavlen);
    elseif isempty(wavlen)
        params.L = 8;
    else
        error(message('Wavelet:modwt:InvalidWavelet'));
    end
validateattributes(params.L,{'numeric'},{'real','positive','even'});

    if (Len>1)
        params.ConfLevel = varargin{2};
            if isempty(params.ConfLevel)
                params.ConfLevel = 0.975;
            end
    end

validateattributes(params.ConfLevel,{'double'},{'>',0,'<',1});