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.