NCL Home>
Application examples>
Data Analysis ||
Data files for some examples
Example pages containing:
tips |
resources |
functions/procedures

# NCL: Spectral Analysis and Complex Demodulation

**Spectral Analysis**
Spectral analysis of time series is the process of partitioning the temporal variance information
into frequency variance information. The latter is called the *spectrum*.
The spectrum breaks the sample variance of time series into discret components,
each of which is associated with a particular frequency.

A nice simple example of the concept and process is provided at
**Introduction to Spectral Analysis** (D. Percival, U. Washington).

The NCL functions **specx_anal** and **specxy_anal**
perform the temporal-to-frequency transformation via the
**Fast Fourier Transform (FFT)**.
The **FFT** implicitly assume the time series are 'cyclic'. Generally, this is not true.
Hence, the source time series must be pre-treated (*aka*, "pre-whitening").
In most applications, this means the data must be tapered in the time or frequency domains.
NCL's functions perform the tapering in the time domain. Further, removing any linear
trend in the source data is recommended. The trend removal is performed first, then the
trend-removed data are tapered. This combination minimizes "ringing" and "leakage".

The FFT performed on the pre-whitened series create periodogram estimates. Nominally,
each periodogram estimate has 2 degrees of freedom for each frequency band.
("Nominally" because the source data has been pre-whitened which slightly reduces
the degrees of freedom..) To provide some level
of statistical reliability, the periodogram estimates must be smoothed.
For example, a 5-point smoother would nominally result in 10 degrees of freedom (5*2).
The side effect of the smoothing process is that the effective band-width is wider
(think, 5*periodogram-band-width). A nice description is at:
**Spectral Analysis: Smoothed Periodogram Method** (DMeko, U. Arizona).

spec_1.ncl: Reads in a time series from
a netCDF file, calculates the spectrum and plots the results.

**specx_anal**: Calculates the spectra of a time series with options on
smoothing, demeaning, detrending, and tapering.

**dimsizes**: calculates the dimension sizes of any array.

spec_2.ncl: Reads in two time series
and calculates their spectra, cospectra, coherence, quadrature spectra,
and phase and creates a panel plot of all values.

**specxy_anal**: Calculates the cospectra, coherence-squared, quadrature
spectra, and phase as well as the spectra of two time series with
options on smoothing, demeaning, detrending, and tapering.

The coherence-squared is a statistic that can be used to examine the relation between
two signals or data sets. It allows identification of significant frequency-domain
correlation between the two time series. **Phase estimates in the cross spectrum are
only useful where significant frequency-domain coherence-squared exists.** The
6.4.0 release contains two functions that can be used to test for significance:
**cohsq_c2p** and **cohsq_p2c**.

**gsn_panel** is the plot interface that panels plots together.
More examples of panel plots.

spec_3.ncl: Calculates red noise
confidence intervals and creates several plot variations.

**specx_ci**
calculates a red noise confidence interval. This function
returns four curves: the calculated spectrum, the "red noise" curve
and the curves indicating the upper and lower confidence bounds.

spec_4.ncl: Plots spectra and
"red noise" confidence intervals on a log scale and draws the bandwidth.

*trYLog* = True

*trYMinF* =
0.10

*trYMaxF* =
30.0, Changes the Y-axis to log and manually
sets the lower adn uppe limits.

**gsn_polyline** and

**gsn_text**
are used to add the bandwidth line and label to the figure.

spec_5.ncl: A variation of Example 4.
The leftmost plot uses frequency as the abscissa. The two rightmost plots
use period (1/fequency) as the abscissa. Since the period is not linear,
plotting a subset of all the periods may be more useful. The

**ind**
function is used to select the desired indices. The rightmost plot illustrates
how to reverse the abscissa order.

**Complex Demodulation:**
As noted by Bloomfield (1976, 2000) "not all 'periodic' phenomena have simple representations in
terms of cosine functions. Complex demodulation is a more flexible approach to the analysis
of such data. By trading off some *frequency resolution* for *time resolution*, complex
demodulation can describe features of data that would be missed by harmonic analysis.
The price of this flexibility is loss of precision in describing pure frequencies,
for which harmonic analysis is most exact."

Complex demodulation may be viewed as a *local* version of harmonic analysis. It enables the
amplitude and phase of a particular frequency component of a time series to be described
as functions of time. It is *local* in that these components are allowed to change slowly over time.

The complex demodulation process, implemented by **demod_cmplx**, consists of several steps:

(1) A demodulation frequency (*frqdem*) is specified.
(2) Compute the real (cosine) and imaginary (sine) components at each time step.
(3) Low pass filter each component using a specified low pass
cutoff frequency (*frqcut*).
(4) Compute the amplitudes [ **A**(t) ] and phases [ **P**(t) ]
using the low pass filtered components (3).

The phases returned in step (4) are calculated via **atan2** and range from -pi to +pi.
These are called 'wrapped phases.' Unfortunately, the resulting phase plots may appear quite noisy and
are difficult to interpret. For graphical and interpretation purposes, the wrapped phases are often
unwrapped. NCL's **unwrap_phase** can be used to perform this task.

Further, as suggested by the
**S-Plus demod** documentation:

**To better understand the results of complex demodulation
several lowpass filters should be tried:** the smaller the pass band,
the less instantaneous in time but more specific in frequency is the result.

In fact, complex demodulation (like spectral analysis), could be viewed to be as much 'art' as 'science'.
Different

*frqcut* and

*nwt* and, in some cases (Example demod_cmplx_2), different

References:

**Fourier Analysis of Time Series**
Peter Bloomfield
Wiley, 2000
WWW:
http://faculty.washington.edu/kessler/generals/complex-demodulation.pdf
http://www.pmel.noaa.gov/maillists/tmap/ferret_users/fu_2007/pdfHMANlJLqCE.pdf

demod_cmplx_1.ncl:
Demodulate Wolf's annual sun spot series using NCL's

**demod_cmplx** function.
The modulation period used is 11 years [

*frqdem*=(1.0/11)].

As noted above, several different low pass (*frqcut*, *nwt*)
filters should be tried. The left figure uses *frqcut*=0.045 (=0.5**frqdem*);
the right figure uses *frqcut*=0.025. In both figures, the *nwt*=41 is kept constant.
That was what Bloomfield (1976) used for his least squares filter.
However, that too could be changed to get different filter effects.

The unwrapped phases created by **unwrap_phase**
are at the bottom of the figure.

**Note: ***frqcut* and *frqdem* are not related. The use of *frqcut*=0.5**frqdem*
in the example (left figure) was an experiment.

demod_cmplx_2.ncl:
Read and unpack daily sea-level-pressures (slp) from grid points surrounding
Kanton Island (2.8N, 188.325E). Perform a simple areal average at each time step.
Calculate daily anomalies from a daily climatology. Perform a complex
demodulation on the anomaly time series about a period of 40 days (

*frqdem*=(1/40).
Commonly, the Madden-Julian Oscillation (MJO) is described as having a 30-60 day period.

Since the MJO spans 30-60 days several different demodulation (*frqdem*) frequencies
should be examined. Again, different *frqcut* and *nwt* should be
examined.

The unwrapped phases created by **unwrap_phase**
are at the bottom of the middle figure.