Time-Series Analysis - MATLAB Recipes for Earth Sciences - page 159

method. h is method was invented by Arthur Schuster in 1898 for studying
the climate and calculates the power spectrum by performing a Fourier
transform directly on a sequence without requiring prior calculation of
the autocorrelation sequence. h e periodogram method can therefore be
considered a special case of the Blackman and Tukey (1958) method, applied
with the time lag k set to unity (Muller and Macdonald 2000). At the time
of its introduction in 1958, the indirect computation of the power spectrum
via an autocorrelation sequence was faster than calculating the Fourier
transformation for the full data series x(t) directly. At er the introduction
of the Fast Fourier Transform (FFT) by Cooley and Tukey (1965), and
subsequent faster computer hardware, the higher computing speed of the
Blackman-Tukey approach compared to the periodogram method became
relatively unimportant.
For this next example we again use the synthetic time series x , xn and xt
generated in Section 5.2 as the input:
clear
t = 1 : 1000; t = t';
x = 2*sin(2*pi*t/50) + sin(2*pi*t/15) + 0.5*sin(2*pi*t/5);
randn('seed',0)
n = randn(1000,1);
xn = x + n;
xt = x + 0.005*t;
We then compute the periodogram by calculating the Fourier transform of
the sequence x . h e fastest possible Fourier transform using fft computes
the Fourier transform for nfft frequencies, where nfft is the next power of
two closest to the number of data points n in the original signal x . Since the
length of the data series is n=1000 , the Fourier transform is computed for
nfft=1024 frequencies, while the signal is padded with nfft-n=24 zeros.
Xxx = fft(x,1024);
If nfft is even, as in our example, then Xxx is symmetric. For example, as the i rst
(1+nfft/2) points in Xxx are unique, the remaining points are symmetrically
redundant. h e power spectral density is dei ned as Pxx2=(abs(Xxx).^2)/Fs ,
where Fs is the sampling frequency. h e function periodogram also scales the
power spectral density by the length of the data series, i.e., it divides by Fs=1
and length(x)=1000 .
Pxx2 = abs(Xxx).^2/1000;