This example shows how to use the cross spectrum to obtain the phase lag between sinusoidal components in a bivariate time series. The example also uses the magnitude-squared coherence (MSC) to identify significant frequency-domain correlation at the sine wave frequencies.
Create the bivariate time series. The individual series consist of two sine waves with frequencies of 100 and 200 Hz in additive white Gaussian noise. The sine waves in the x-series both have amplitudes equal to 1. The 100-Hz sine wave in the y-series has amplitude 0.5 and the 200-Hz sine wave in the y-series has amplitude 0.35. The sine waves in the y-series are phase-lagged by π/4 radians (100 Hz) and π/2 radians (200 Hz). You can think of y-series as the noise-corrupted output of a linear system with input x. In the following code, set the random number generator to the default settings for reproducible results.
Fs = 1000; t = 0:1/Fs:1-1/Fs; rng default; x = cos(2*pi*100*t)+sin(2*pi*200*t)+0.5*randn(size(t)); y = 0.5*cos(2*pi*100*t-pi/4)+0.35*sin(2*pi*200*t-pi/2)+0.5*randn(size(t));
Obtain the magnitude-squared coherence (MSC) for the bivariate time series. The magnitude-squared coherence enables you to identify significant frequency-domain correlation between the two time series. Phase estimates in the cross spectrum are only useful where significant frequency-domain correlation exists.
To prevent obtaining a magnitude-squared coherence estimate, which is identically 1 for all frequencies, you must use an averaged MSC estimator. Both Welch's overlapped segment averaging (WOSA) and mulitaper techniques are appropriate. mscohere implements a WOSA estimator.
Set the window length to 100 samples. This window length contains 10 periods of the 100-Hz sine wave and 20 periods of the 200-Hz sine wave. Use an overlap of 80 samples with the default Hamming window. Plot the magnitude-squared coherence.
[Pxy,F] = mscohere(x,y,hamming(100),80,100,Fs); plot(F,Pxy,'linewidth',2); title('Magnitude-squared Coherence'); xlabel('Hz'); grid on;
You see that the magnitude-squared coherence is greater than 0.8 at 100 and 200 Hz.
Obtain the cross spectrum of x and y using cpsd. Use the same parameters to obtain the cross spectrum that you used in the MSC estimate. Plot the phase of the cross spectrum and indicate the frequencies with significant coherence between the two times. Mark the known phase lags between the sinusoidal components.
[Cxy,F] = cpsd(x,y,hamming(100),80,100,Fs); figure; plot(F,-angle(Cxy),'linewidth',2); title('Cross Spectrum Phase'); xlabel('Radians'); set(gca,'xtick',[100 200]); set(gca,'ytick',[-pi -pi/2 -pi/4 0 pi/4 pi/2 pi]); grid on;
You see that at 100 Hz and 200 Hz, the phase lags estimated from the cross spectrum are close to the true values.
In this example, the cross spectrum estimates are spaced at 1000/100= 10 Hz. You can return the phase estimates at those frequency bins. Keep in mind that the first frequency bin corresponds to 0 Hz, or DC.
phi100 = - angle(Cxy(11)); phi200 = - angle(Cxy(21));
You see that phi100 and phi200 are close to –π/4 and –π/2.