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:

terms describing waveform

Figure 15-1

Using this figure we can introduce some definitions: 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
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: eq 16-1

(recall that eq 16-1 so the Fourier transform has real and imaginary components).

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

Eq 15-2: eq 16-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: eq 16-3

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

Eq 15-4: eq 16-4

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

For such discrete time series, we consider measurements of h taken at discrete intervals, Delta:

Eq 15-5: eq 16-5

The quantity 1/Delta is the sampling rate. For any particular sampling interval there is an frequency fsub c called the Nyquist critical frequency:

Eq 15-6: eq 16-6

What is the significance of the critical frequency? There are two related issues.

To summarize:

  1. Unless H(f) is zero once f exceeds 1/(2Delta), 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 ucsub b. A mass balance can be written:

Eq 15-7: eq 16-7

or rewriting in terms of a periodic variation in the delivery of material:

Eq 15-8: eq 16-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: eq 16-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: eq 16-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: eq 16-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 Delta then these occur at times tsub k=kDelta with k=0,...,N-1. We introduce the notation for the values at these times: hsub kidenth(tsub k) 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: eq 16-12

With these conventions the discrete Fourier transform is written:

Eq 15-13: eq 16-13

The power at each of these frequencies is:

Eq 15-14: eq 16-14

This information can be plotted in the form of a periodiogram:



Figure 15-5



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: eq 16-15

Then the power density is given by:

Eq 15-16: eq 16-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:

windowing function 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.


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: eq 16-17

Thus the autocorrelation of a time series onto itself is given by:

Eq 15-18: eq 16-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.
Fourier wavelet comparison

Figure 15-9

We can review these ideas with this Matlab demonstration.
Some resources for more detail on these topics:

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

Valid HTML 4.0!

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