Oceanography 540--Marine Geological Processes--Winter Quarter 2001

### Time Series Analysis

Time series of oceanic phenomena often contain periodic components related to forcing at a wide range of time scales: waves, tides and tidal currents, diurnal and annual cycles, ENSO, Pacific Decadal Oscillation, and orbital geometry, its influence on incoming solar radiation and Pleistocene climate. Time series analysis involves the extraction of information about these periodic components from time series data.

A periodic signal can be represented by its intensity or amplitude as a function of time: Figure 15-1

Using this figure we can introduce some definitions:
• The peak intensity of the signal (relative to the mean) is the amplitude, A.
• The waveform repeats with period T.
• Time can be expressed in either time units, t, or as an angle (expressed in radians).
• The angular scale is related to the time scale via the relationship: where is the angular frequency and is the phase angle.
• The period is related to the angular frequency by the relationship: • The frequency f is the inverse of the period f=1/T.
• Thus =2 f.
• With these definitions, the signal can be represented as: A physical process can generally be described by the values of some parameter h that varies with time, h=h(t). We refer to this representation as being in the time domain. An alternative and equivalent representation of the same signal would be to express it as its amplitude, H, as a function of frequency, f: H=H(f) . We refer to this representation as being in the frequency domain.

Time DomainFrequency Domain  Figure 15-2. The period of the waveform is 36 time units, so that the frequency is 1/36=.0278. (The vertical axis on the right hand side is power, not amplitude, see below).
These two representations of the same signal can be related by the Fourier transform:

Eq 15-1:  (recall that so the Fourier transform has real and imaginary components).

The original time series can be obtained from the inverse transform:

Eq 15-2:  The transforms can be thought of as weighted averages of waveforms of a given frequency; the forward transform to the frequency domain weighted by the signal amplitude h(t) and the inverse transform to the time domain weighted by the frequency amplitude H(f). In other words, when the signal h has the same frequency as the waveform, the the integral H is large.

If h(t) is real (we are generally considering real signals) then H(f) is generally a complex number. However under this condition H(-f)=H(f)*. (the notation * denotes the complex conjugate). This symmetry means that for real signals there is no unique information in the negative frequencies.

The power of the signal, P, is a measure of its squared amplitude. The power at a frequency f is defined as:

Eq 15-3:  and so the total power in a signal is (fix equation, no (f) on lhs):

Eq 15-4:  Equations 16-1 through 16-4 apply to a continuous function spanning all times. In general though real signals will:

• be represented by discrete measurements
• these measurements will span a finite interval
For such discrete time series, we consider measurements of h taken at discrete intervals, :

Eq 15-5:  The quantity 1/ is the sampling rate. For any particular sampling interval there is an frequency f called the Nyquist critical frequency:

Eq 15-6:  What is the significance of the critical frequency? There are two related issues.

• The first issue involves a concept called the sampling theorem which says:

If the spectrum of h(t) is zero outside of the frequency interval -fc<f<fc, then all spectral information is obtained if the signal is sampled at times separated by no more than 1/(2 ). This means that we gain no spectral information by sampling at a rate greater than the Nyquist frequency. We can see why by examining the periodic signal in this figure:  Figure 15-3 The period of the signal is 36 time units and so the Nyquist frequency is 1/18 time units. If we sample time unit by time unit, we get a smooth curve. Even if we sample by 18 time units, the Nyquist period, we get (provided we are not exactly in phase with the zero crossings) the highs and the lows so that we still capture the period of the signal. If we sample by 36 time units though, we only get a single value--the signal is constant through time.

Insuring that H(f)=0 beyond the Nyquist frequency is achieved in practice by applying a low pass filter (a filter that lets low frequencies pass the filter, but not high frequency components) to the signal before sampling (as we will see, nature applies a low pass filter through bioturbation in the surface-most sediments).

• The second issue is a phenomenon called aliasing. If there is power outside the frequency band -f  f f then discrete sampling folds that power into this frequency band. Schematically:  Figure 15-4. The width of the true transform is 2fc; the band from fc to 3fc maps to fc to -fc, the band from 3 fc to 5fc maps to -fc to fc, and so on.) To summarize:

1. Unless H(f) is zero once f exceeds 1/(2 ), there will be aliasing by higher frequency components.
2. To avoid aliasing, either the sampling interval must be chosen to insure that this condition is met or the original signal must be filtered to remove high frequency components before sampling.

To understand the nature of the filtering effect of bioturbation consider as an example a core taken in sediments accumulating at 3 cm/ky and sampled at 10 cm intervals. Thus, the sampling interval is 3000 y and so the Nyquist frequency is 1/6000 y.

To model bioturbation, consider a well mixed box of height h at the sediment surface with material delivered at rate F and removed as the product of sedimentation rate times concentration uc . A mass balance can be written:

Eq 15-7:  or rewriting in terms of a periodic variation in the delivery of material:

Eq 15-8:  This differential equation can be solved for some arbitrary initial condition to find that at large times, when initial conditions have decayed, there is a stationary solution:

Eq 15-9:  This equation can be simplified at the limits of high and low frequency. When f is large (the supply of material to the sediment varies at high frequency) it reduces to:

Eq 15-10:  In this case mixing averages the fluctations in the supply of rain to the sediment. When f is small (the supply to the sediment varies slowly) it reduces to:

Eq 15-11:  In this case the periodic signal is preserved. In other words, low frequency signals are preserved while high frequency ones are not; the system acts as a low pass filter. The boundary between "low" and "high" frequency is effectively f~u/h. For a depth of bioturbation of 2-10 cm in our sample core, this frequency is somewhere between 1/600 y and 1/3000 y. By carefully choosing the sampling interval, little aliasing is possible. (And of course we can also apply a low pass filter to be sure.) Conversely, a limit exists for the highest frequency signal that can be resolved in a core.

Mathematica is a useful tool for solving Eq 15-8 | Download Notebook The time series analysis of a record consists of a series of computational steps. Depending on the application these might include:

1. detrending and whitening
2. autocorrelation
3. windowing
4. spectral estimation
We will discuss these in the reverse order from which they are performed.

Spectral estimation

If we have N measurements of a parameter h sampled at an interval then these occur at times t =k with k=0,...,N-1. We introduce the notation for the values at these times: h  h(t ) From these N measurements we can derive independent information at only N discrete frequencies--these are normally spaced evenly within the Nyquist limits:

Eq 15-12:  With these conventions the discrete Fourier transform is written:

Eq 15-13:  The power at each of these frequencies is:

Eq 15-14:  This information can be plotted in the form of a periodiogram:  Figure 15-5 Windowing

A common problem is that these estimates are not within the bin implied by the geometry; they actually have a spread several bins wide. A process termed data windowing minimizes this leakage between bins. It involves applying a windowing function to the data before forming the transform:

Eq 15-15:  Then the power density is given by:

Eq 15-16:  The origin of the leakage between bins arises in the process of taking a portion out of a long time series. Consider this time series h and window function w: Figure 15-6

Extracting a discrete portion out of the time series h is equivalent to multiplying h by a window function w where w=1 in the region of interest and 0 everywhere else. The convolution theorem say that for two functions h(t) and w(t), if p=wh then P(f)=W(f)H(f). When the discrete transform is taken, the high frequency components associated with the steps in the function w are aliased into the Nyquist interval. The characteristic of a good window is to smoothly turn "on and off", thereby minimizing the introduction of high frequency components into the transform.

Autocorrelation

In performing time series analysis, one can work with autocorrelated records to produce a correlogram. If the Fourier transform of h(t) is H(f), and the autocorrelation function of h(t) is r, then the Fourier transform R(f)=|G(f)|2. For N data points taken from two times series g and h, the correlation as a function of the lag j is:

Eq 15-17:  Thus the autocorrelation of a time series onto itself is given by:

Eq 15-18:  Imagine a periodic function that is shifted progressively to the right:  Figure 15-7 As it is shifted the correlation will drop, then rise to the same value when the function has been shifted exactly one cycle:  Figure 15-8 Thus the autocorrelation produces a different time series, but with the same spectral components.

Detrending, whitening

There are generally long period and d.c. components that give rise to power in a spectra being concentrated at low frequencies. High pass filters can be applied to remove the d.c. baseline (detrending) and long period, low frequency components (whitening). The effect is to emphasize peaks at intermediate frequency, which no longer sit on the shoulder of a larger low frequency peak.

Thus far we have emphasized frequency characteristics of a time series. Time series analysis also includes consideration of the frequency components within the time domain. These are separated by:

1. transforming to the frequency domain
2. multiplying the transform by a bandpass filter. This product contains only the information contained within this frequency band
3. applying the inverse transform back to the time domain
Thus far we have been considering time series that are stationary, i.e., whose frequency content is not changing with time. We may also be interested in phenomena that do not satisfy this constraint. A traditional approach is to apply Fourier analysis to overlapping subsets of the time series in order to examine the evolution of the frequency content. A more recently developed technique is wavelet analysis (for a good review see (42)). Figure 15-9 is adapted from that reference and compares Fourier and wavelet analysis of a stationary time series with two superimposed sine waves and a time series containing signals of the same frequency content but in different segments of the observation period. Figure 15-9

We can review these ideas with this Matlab demonstration.
Some resources for more detail on these topics:
• W.H. Press, et al., Numerical Recipes in <C | Fortran | Pascal>
chapter on Fourier and Spectral Applications develops background for some basic algorithms
• J.S. Bendat and A.G. Piersol, Random Data: Analysis and Measurement Procedures
survey of data analysis, including spectral methods
• L.H. Koopmans, The Spectral Analysis of Time Series
a comprehensive treatment, used for Charlie Eriksen's Data Analysis course

Next Lecture | Lecture Index | Search the Ocean 540 Pages | Ocean 540 Home

 Oceanography 540 Pages Pages Maintained by Russ McDuff (mcduff@ocean.washington.edu) Copyright (©) 1994-2001 Russell E. McDuff and G. Ross Heath; Copyright Notice Content Last Modified 1/19/2001 | Page Last Built 1/19/2001