Page 99 - MATLAB Recipes for Earth Sciences

92 5 Time-Series Analysis

where M is the maximum lag and f the sampling frequency. The Blackman-
s
Tukey power spectral density PSD is estimated by

The actual computation of PSD can be performed only at a fi nite number
of frequency points by employing a Fast Fourier Transformation (FFT).
The FFT is a method to compute a discrete Fourier Transform with reduced
execution time. Most FFT algorithms divide the transform into two pieces
of size N/2 at each step. It is therefore limited to blocks of power of two.
In practice, the PSD is computed by using N squared number of frequen-
cies. The actual number of frequencies used lies close to the number of data
points in the original signal x(t).
The discrete Fourier transform is an approximation of the continu-
ous Fourier transform. The Fourier transform expects an infi nite signal.
However, real data are limited at both ends, i.e., the signal amplitude is zero
beyond the limits of the time series. In the time domain, a finite signal cor-

responds to an infinite signal multiplied by a rectangular window that is one
within the limits of the signal and zero elsewhere. In the frequency domain,
the multiplication of the time series with this window equals to a convolu-
tion of the power spectrum of the signal with the spectrum of the rectangular
window. The spectrum of the window, however, equals a sin(x)/x function,
which has a main lobe and several side lobes at both sides of the main peak.
Therefore all maxima in a power spectrum leak, i.e., they lose power with
respect to the minor peaks (Fig. 5.4).
A popular way to overcome the problem of spectral leakage is windowing.
The sequence of data is simply multiplied by a window with smooth ends.
Several window shapes are available, e.g., Bartlett (triangular), Hamming
(cosinusoidal) and Hanning (slightly different cosinusoidal). The use of
these windows slightly modifies the equation of the power spectral density.

where M is the maximum lag considered and window length, and w(k) is the
windowing function. The Blackman-Tukey method therefore performs au-
tospectral analysis in three steps, calculation of the autocorrelation sequence
corr (k), windowing and fi nally computation of the discrete fourier trans-
xx
form. MATLAB allows to perform power spectral analysis with a number of
modifi cations of the above method. A useful modifi cation is the method by