Page 107 - MATLAB Recipes for Earth Sciences

100 5 Time-Series Analysis

The phase shift at a frequency of f=0.2 (period τ=5) can be interpolated from
the phase spectrum

interp1(f,phase,0.2)

which produces the output
ans =
1.2568
The phase spectrum is normalized to one full period τ=2›, therefore a phase
shift of 1.2568 equals (1.2568*5)/(2*›) = 1.0001, which is the phase shift of
one that we introduced at the beginning.
We now use two sine waves with different periodicities to illustrate
crossspectral analysis. The both have a periodicity of 5, but with a phase
shift of 1, then they have both one other period, which are different, how-
ever.

clear
t = 0.01 : 0.1 : 1000;
y1 = sin(2*pi*t/15) + 0.5*sin(2*pi*t/5);
y2 = 2*sin(2*pi*t/50) + 0.5*sin(2*pi*t/5+2*pi/5);
plot(t,y1,'b-',t,y2,'r-')

Now we compute the crossspectrum, which clearly shows the common pe-
riod of τ=5 or frequency of f=0.2.

[Pxy,f] = cpsd(y1,y2,[],0,512,10);
magnitude = abs(Pxy);

plot(f,magnitude);
xlabel('Frequency')
ylabel('Power')
title('Cross PSD Estimate via Welch')

The coherence shows a large value of approximately one at f=0.2.

[Cxy,f] = mscohere(y1,y2,[],0,512,10);
plot(f,Cxy)
xlabel('Frequency')
ylabel('Magnitude Squared Coherence')
title('Coherence Estimate via Welch')
The complex part is required for calculating the phase shift between the two
sine waves.