Geoscience Reference
In-Depth Information
h e resulting spectrum no longer shows the low-frequency peak.
[Pxxt,f] = periodogram(xt,[],1024,1);
[Pxxdt,f] = periodogram(xdt,[],1024,1);
subplot(1,2,1)
plot(f,Pxx), grid
xlabel('Frequency')
ylabel('Power')
subplot(1,2,2)
plot(f,Pxxdt), grid
xlabel('Frequency')
ylabel('Power')
Some data contain a high-order trend that can be removed by i tting a higher-
order polynomial to the data and subtracting the corresponding
x
(
t
) values.
We now use two sine waves with identical periodicities ˄=5 (equivalent
to
f
=0.2) and amplitudes equal to two to compute the cross-spectrum of two
time series. h e sine waves show a relative phase shit of
t
=1. In the argument
of the second sine wave this corresponds to 2ˀ/5, which is one i t h of the full
wavelength of ˄=5.
clear
t = 1 : 1000;
x = 2*sin(2*pi*t/5);
y = 2*sin(2*pi*t/5 + 2*pi/5);
plot(t,x,'b-',t,y,'r-')
axis([0 50 -2 2]), grid
h e cross-spectrum is computed by using the function
cpsd
, which uses
Welch's method for computing power spectra (Fig. 5.8).
Pxy
is complex and
contains both amplitude and phase information.
[Pxy,f] = cpsd(x,y,[],0,1024,1);
plot(f,abs(Pxy)), grid
xlabel('Frequency')
ylabel('Power')
title('Cross-Spectrum')
h e function
cpsd(x,y,window,noverlap,nfft,fs)
specii es the number of
FFT points
nfft
used to calculate the cross power spectral density, which is
1024 in our example. h e parameter
window
is empty in our example and the
default rectangular window is therefore used to obtain eight sections of
x
and
y
. h e parameter
noverlap
dei nes the number of overlapping samples, which
is zero in our example. h e sampling frequency
fs
is 1 in this example. h e
