www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/WaveletCrossCorrelationExample.m

    %% Wavelet Cross-Correlation for Lead-Lag Analysis
% This example shows how to use wavelet cross-correlation to measure
% similarity between two signals at different scales.
%%
% Wavelet cross-correlation is simply a scale-localized version of the
% usual cross-correlation between two signals. In cross-correlation, you
% determine the similarity between two sequences by shifting one relative
% to the other, multiplying the shifted sequences element by element and
% summing the result. For deterministic sequences, you can write this as an
% ordinary inner product: $<x_n,y_{n-k}>_n = \sum_{n} x_n \bar{y}_{n-k}$
% where $x_n$ and $y_n$ are sequences (signals) and the bar denotes complex
% conjugation. The variable, k, is the lag variable and represents the
% shift applied to $y_n$. If both $x_n$ and $y_n$ are real, the complex
% conjugate is not necessary. Assume that $y_n$ is the same sequence as
% $x_n$ but delayed by L>0 samples, where L is an integer. For
% concreteness, assume $y_n = x_{n-10}$. If you express $y_n$ in terms of
% $x_n$ above, you obtain $<x_n,x_{n-10-k}>_n = \sum_{n} x_n
% \bar{x}_{n-10-k}$. By the Cauchy-Schwartz inequality, the above is
% maximized when $k= -10$. This means if you left shift (advance) $y_n$ by
% 10 samples, you get the maximimum cross-correlation sequence value. If
% $x_n$ is a L-delayed version of $y_n$, $x_n = y_{n-L}$, then the
% cross-correlation sequence is maximized at $k=L$.
%%
% If you have the Signal Processing Toolbox(TM), you can show this by using
% |xcorr|. Create a triangular signal consisting of 20 samples. Create a
% noisy shifted version of this signal. The shift in the peak of the
% triangle is 3 samples. Plot the |x| and |y| sequences.

% Copyright 2015 The MathWorks, Inc.

x = triang(20);
rng default;
y = [zeros(3,1);x]+0.3*randn(length(x)+3,1);

subplot(2,1,1)
stem(x,'filled')
axis([0 22 -1 2])
title('Input Sequence')

subplot(2,1,2)
stem(y,'filled')
axis([0 22 -1 2])
title('Output Sequence')
%%
% Use |xcorr| to obtain the cross-correlation sequence and determine the
% lag where the maximum is obtained.
[xc,lags] = xcorr(x,y);
[~,I] = max(abs(xc));
figure
stem(lags,xc,'filled')
legend(sprintf('Maximum at lag %d',lags(I)))
title('Sample Cross-Correlation Sequence')
%%
% The maximum is found at a lag of -3. The signal |y| is the
% second input to |xcorr| and it is a delayed version of |x|. You have to
% shift |y| 3 samples to the left (a negative shift) to maximize the cross
% correlation. If you reverse the roles of |x| and |y| as inputs to
% |xcorr|, the maximum lag now occurs at a positive lag. 
[xc,lags] = xcorr(y,x);
[~,I] = max(abs(xc));
figure
stem(lags,xc,'filled')
legend(sprintf('Maximum at lag %d',lags(I)))
title('Sample Cross-Correlation Sequence')
%%
% |x| is an advanced version of |y| and you delay |x| by three samples
% to maximize the cross correlation. 
%%
% |modwtxcorr| is the scale-based version of |xcorr|. To use |modwtxcorr|,
% you first obtain the nondecimated wavelet transforms. 
%%
% Apply wavelet cross-correlation to two signals that are shifted versions
% of each other. Construct two exponentially-damped 200-Hz sine waves with
% additive noise. The |x| signal has its time center at $t=0.2$ seconds
% while |y| is centered at $t=0.5$ seconds.
t = 0:1/2000:1-1/2000;
x = sin(2*pi*200*t).*exp(-50*pi*(t-0.2).^2)+0.1*randn(size(t));
y = sin(2*pi*200*t).*exp(-50*pi*(t-0.5).^2)+0.1*randn(size(t));
plot(t,x)
hold on
plot(t,y)
xlabel('Seconds')
ylabel('Amplitude')
grid on
legend('x','y')
%%
% You see that |x| and |y| are very similar except that |y| is delayed by
% 0.3 seconds. Obtain the nondecimated wavelet transform of |x| and |y|
% down to level 5 using the Fejer-Korovkin (14) wavelet. The wavelet
% coefficients at level 3 with a sampling frequency of 2 kHz are an
% approximate $[2000/2^4, 2000/2^3)$ bandpass filtering of the inputs.
% The frequency localization of the Fejer-Korovkin filters ensures that
% this bandpass approximation is quite good.
wx = modwt(x,'fk14',5);
wy = modwt(y,'fk14',5);
%%
% Obtain the wavelet cross-correlation sequences for the wavelet transforms
% of |x| and |y|. Plot the level 3 wavelet cross-correlation sequence for
% 2000 lags centered at zero lag. Multiply the lags by the sampling period
% to obtain a meaningful time axis.
[xc,~,lags] = modwtxcorr(wx,wy,'fk14');
lev = 3;
zerolag = floor(numel(xc{lev})/2+1);
tlag = lags{lev}(zerolag-999:zerolag+1000).*(1/2000);
figure;
plot(tlag,xc{lev}(zerolag-999:zerolag+1000));
title('Wavelet Cross-Correlation Sequence (level 3)');
xlabel('Time')
ylabel('Cross-Correlation Coefficient');
%%
% The cross-correlation sequence peaks at a delay of -0.3 seconds. The
% wavelet transform of |y| is the second input to |modwtxcorr|. Because the
% second input of |modwtxcorr| is shifted relative to the first, the peak
% correlation occurs at a negative delay. You have to left shift (advance)
% the cross-correlation sequence to align the time series. If you reverse
% the roles of the inputs to |modwtxcorr|, you obtain the peak correlation
% at a positive lag.
[xc,~,lags] = modwtxcorr(wy,wx,'fk14');
lev = 3;
zerolag = floor(numel(xc{lev})/2+1);
tlag = lags{lev}(zerolag-999:zerolag+1000).*(1/2000);
figure;
plot(tlag,xc{lev}(zerolag-999:zerolag+1000));
title('Wavelet Cross-Correlation Sequence (level 3)');
xlabel('Time')
ylabel('Cross-Correlation Coefficient');
%%
% To show that wavelet cross-correlation enables scale(frequency)-localized
% correlation, plot the cross-correlation sequences at levels 1 and 5.
lev = 1;
zerolag = floor(numel(xc{lev})/2+1);
tlag = lags{lev}(zerolag-999:zerolag+1000).*(1/2000);
plot(tlag,xc{lev}(zerolag-999:zerolag+1000));
title('Wavelet Cross-Correlation Sequence (level 1)');
xlabel('Time')
ylabel('Cross-Correlation Coefficient');
figure;
lev = 5;
zerolag = floor(numel(xc{lev})/2+1);
tlag = lags{lev}(zerolag-999:zerolag+1000).*(1/2000);
plot(tlag,xc{lev}(zerolag-999:zerolag+1000));
title('Wavelet Cross-Correlation Sequence (level 5)');
xlabel('Time')
ylabel('Cross-Correlation Coefficient');
%%
% The wavelet cross-correlation sequences at levels 1 and 5 do not show any
% evidence of the exponentially-weighted sinusoids due to the bandpass
% nature of the wavelet transform.
%%
% With finanical data, there is often a leading or lagging relationship
% between variables. In those cases, it is useful to examine the
% cross-correlation sequence to determine if lagging one variable with
% respect to another maximizes their cross-correlation. To illustrate this,
% consider the correlation between two components of the GDP -- personal
% consumption expenditures and gross private domestic investment. The data
% is quarterly chain-weighted U.S. real GDP data for 1974Q1 to 2012Q4. The
% data were transformed by first taking the natural logarithm and then
% calculating the year-over-year difference. Look at the correlation
% between two components of the GDP -- personal consumption expenditures,
% |pc|, and gross private domestic investment, |privateinvest|.
load GDPcomponents;
piwt = modwt(privateinvest,'fk8',5);
pcwt = modwt(pc,'fk8',5);
figure;
modwtcorr(piwt,pcwt,'fk8')
%%
% Personal expenditure and personal investment are negatively correlated
% over a period of 2-4 quarters. At longer scales, there is a strong
% positive correlation between personal expenditure and personal
% investment. Examine the wavelet cross-correlation sequence at the scale
% representing 2-4 quarter cycles. Plot the cross-correlation sequence
% along with 95% confidence intervals.

[xcseq,xcseqci,lags] = modwtxcorr(piwt,pcwt,'fk8');
zerolag = floor(numel(xcseq{1})/2)+1;
figure;
plot(lags{1}(zerolag:zerolag+20),xcseq{1}(zerolag:zerolag+20));
hold on;
plot(lags{1}(zerolag:zerolag+20),xcseqci{1}(zerolag:zerolag+20,:),'r--');
xlabel('Lag (Quarters)');
ylabel('Cross-Correlation');
grid on;
title({'Wavelet Cross-Correlation Sequence -- [2Q,4Q)'; ...
    'Personal Consumption and Private Investment'});
%%
% The finest-scale wavelet cross-correlation sequence shows a peak positive
% correlation at a lag of one quarter. This indicates that personal
% investment lags personal expenditures by one quarter. If you take that
% lagging relationship into account, then there is a positive correlation
% between the GDP components at all scales.