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}{\underset{}{\overset{\u203e}{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}{\underset{}{\overset{\u203e}{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 maximum 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™, 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.

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)); figure 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 financial 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.