Frequency Analysis of Time Series

from Essential Mathematica for Students of Science
© James J. Kelly, 1998
Last revised: January, 2006

Introduction

An important method for evaluating the dynamics of a system is Fourier analysis of time series.  Suppose that a sequence of n measurements (Subscript[x, m]) = x[(Subscript[t, m])] = x[(m - 1) τ] is made at equally spaced times covering an observation period T = (n - 1) τ.  The Fourier series is then

\(TraditionalForm\`x\&~\_k = \(1\/\@n\) \(∑\+\(m = 1\)\%n\( x\_m\) 
            Exp[\(\(2  π\ i\)\/n\) \((m - 1)\) \((k - 
                    1)\)]\)\)

such that

\(TraditionalForm\`x\_m = \(1\/\@n\) \(∑\+\(k = 1\)\%n\( x\&~\_k\) 
            Exp[\(-\(\(2  π\ i\)\/n\)\) \((k - 1)\) \((m - 
                    1)\)]\)\)

is represented as a sum over Fourier components with amplitudes \(TraditionalForm\`x\&~\_k\) and frequencies (Subscript[ν, k]) = ((k - 1)/T).  Note that the Fourier series can be defined without explicit reference to the time scale by choosing units in which τ -> 1.  The sampling rate is given by (τ^(-1)) and the maximum frequency for which meaningful information can be extracted is (((2  τ))^(-1)) because a minimum of two samples per period are needed to determine a Fourier amplitude.  This critical or maximum frequency is known in the literature as the Nyquist frequency.

If the input signal is periodic and the measurements cover an integral number of periods, the Fourier series will be strongly dominated by that frequency.  The Fourier spectrum of a quasi-periodic signal typically exhibits a strong peak and several smaller peaks at multiples of the fundamental frequency.  Random or chaotic signals, on the hand, tend to produce broad, nearly featureless, spectral distributions.  Therefore, Fourier analysis provides a valuable technique for isolating periodic components from background noise.  Signal processing can be performed automatically by specialized electronic equipment or numerically by analyzing data for sampled signals.  Electronic or numerical algorithms can be used to filter the signals and to compensate for distortion or bandwidth limitations.  In this notebook we will examine the spectral properties of several types of time signals, but will not explore the more sophisticated techniques available to modern signal processing.

Although frequency distributions are generally complex functions, the relative strength of various frequency components can be assessed using the power spectrum

\(TraditionalForm\`P\_k = \(\(2\/n\) |x\&~\_k|\^2 = \(2\/n\) 
          x\&~\_k\ \(x\&~\_k\^*\)\)\)

which has been normalized so that a sinusoidal function with frequency (Subscript[ν, k]) and unit amplitude has unit power divided equally between the peak at (Subscript[ν, k]) and its reflection at (Subscript[ν, 1 + n - k]).  [The reflection symmetry is discussed later.]  This definition is similar to that used for mechanical or electrical oscillations for which the energy is proportional to the square of the displacement, hence the description of this quantity as power, but perhaps it should be more accurately named spectral intensity because it is can also be used for quantities without an associated energy.

The function Fourier[TimeSeries] accepts a list of n equally spaced time samples (Subscript[x, m]) and returns a list of Fourier amplitudes \(x\&~\_k\) where (Subscript[t, m]) = (m - 1) τ are times and (Subscript[ν, k]) = ((k - 1)/T)are frequencies expressed in terms of the unspecified time interval T = (n - 1) τ.  The function InverseFourier[FrequencySeries] accepts a list of equally spaced Fourier amplitudes and returns a time series constructed by inverse Fourier transformation.

Obviously, the Fourier series is closely related to the Fourier transform and is often called the discrete Fourier transform.  Although the mathematical properties of Fourier integrals are simpler in many respects than those of Fourier series, the series is directly applicable to data or to approximate numerical calculations for which observations or functional values are available only for discrete times.  The Fourier transformation, on the other hand, requires continuous functions that can be evaluated for any value of their arguments.  Another important difference is that Fourier series implicitly assume that signals are periodic in time whereas Fourier transformations over an infinite domain apply equally well to nonperiodic signals.  Nevertheless, Fourier series can be used to analyze nonperiodic signals within the observation period T provided that one does not expect that representation to apply to times outside the observed range.

A valuable reference which provides greater mathematical and technical detail and extensive references to the technical literature is:
	Numerical Recipes — the Art of Scientific Computing, by W.H. Press, B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling (Cambridge University Press, Cambridge, 1989).

Initialization

ClearAll["Global`*"];
Off[General::spell,  General::spell1];
$TextStyle = {FontFamily -> "Times", 
          FontSize -> 12};

Needs["Graphics`Graphics`"];
Needs["Utilities`FilterOptions`"]

Implementation

It is useful to define functions which evaluate and plot the power spectrum for a given time series.  The function PowerSpectrum evaluates the power spectrum.  The option ScaleType->Symmetric transforms the frequency scale to the range ±(1/2), whereas ScaleType->Unity maps frequency onto the range {0,1}.  The function PlotPowerSpectrum produces the plot.  The option PlotType->ListPlot yields a linear plot while the default PlotType->LogListPlot yields a semi-log plot.

spectrum[TimeSeries_List, opts___Rule] := 
      Module[{FrequencyDistribution, PowerDistribution}, 
        FrequencyDistribution = Fourier[TimeSeries]; 
        PowerDistribution = 
          FrequencyDistribution Conjugate[FrequencyDistribution] 
            2/Length[TimeSeries]; 
        LogListPlot[ PowerDistribution,  opts ] ]

\(\(PowerSpectrum::usage = \*""PowerSpectrum[TimeSeries,opts___] evaluates the spectral intensity for a list of samples equally spaced in time.  The option ScaleType can be used to transform the frequency index. \n ScaleType->Symmetric displays frequencies on the scale {-(1/2),(1/2)}. \n ScaleType->Unity displays frequencies {0,1}."";\)
	\)
\(\(\(Options[PowerSpectrum] = {ScaleType -> Natural};\)\(
\)
    \)\)
\(PowerSpectrum[TimeSeries_List, opts___Rule] := 
      Module[{FrequencyDistribution, PowerDistribution, n, frequencies, 
          type}, n = Length[TimeSeries]; 
        FrequencyDistribution = Fourier[TimeSeries]; 
        type = \(ScaleType /. {opts}\) /. Options[PowerSpectrum]; 
        frequencies = 
          Switch[ToString[type], "Unity", 
            Table[i\/n, {i, 0, n - 1}], "Symmetric", 
            Table[If[
                i ≤ 
                  1 + n\/2, \(i - 1. \)\/\(n - 1. \), \(i - 1.  - n\)\/\(n - 1. \)], {i, 1, n}], _, Range[n]]; PowerDistribution = {}; 
        Do[AppendTo[
            PowerDistribution, {frequencies[[
                i]], \(2\/n\) 
                FrequencyDistribution[[
                  i]]\ Conjugate[
                  FrequencyDistribution[[
                    i]]]}], {i, 1, n}]; 
        Return[PowerDistribution];\ ]\)

\(\(\(PlotPowerSpectrum::usage = "PlotPowerSpectrum[PowerData,opts] plots a list of data for a power spectrum.  The option PlotType->LogListPlot or PlotType->ListPlot distinguishes between semi-log and linear plots, with LogListPlot being the default.  Additional options are passed to the plotting function.";\)\(
\)
    \)\)
\(\(\(Options[
          PlotPowerSpectrum] = {PlotType -> LogListPlot};\)\(
\)
    \)\)
PlotPowerSpectrum[data_List, opts___Rule] := 
      Module[{type, PlotFunction, optPlot}, 
        PlotFunction = PlotType /. {opts} /. Options[PlotPowerSpectrum]; 
        optPlot = FilterOptions[PlotFunction, opts]; 
        PlotFunction[data, optPlot, Frame -> True, 
          FrameLabel -> {"frequency", "power"}]]

General characteristics of Fourier series

To interpret power spectra, one must be aware of several important properties of discrete Fourier series.  Although they can be derived rigorously, we prefer to provide examples which illustrate these properties in action.  When evaluating these examples, please identify the relevant property.

Superposition: if (Subscript[f, m]) = (Subscript[x, m]) + (Subscript[y, m]) is the sum of two input signals, its Fourier transform \(TraditionalForm\`f\&~\_m = x\&~\_m + y\&~\_m\)  is the sum of the Fourier transforms of each signal independently.

If {(Subscript[x, k])} is purely real, the Fourier components satisfy the reflection symmetry \(TraditionalForm\`\(x\&~\_k\^*\) = \ x\&~\_\(n + 1 - k\)\).  Consequently, the power spectrum for a real function is symmetric about its midpoint and the maximum frequency that contains independent information is n/(2  T).

The inverse Fourier transform

\(TraditionalForm\`x\_m = \(1\/\@n\) \(∑\+\(k = 1\)\%n\( x\&~\_k\) 
            Exp[\(-\(\(2  π\ i\)\/n\)\) \((k - 1)\) \((m - 
                    1)\)]\)\)

is a periodic function, such that (Subscript[x, n + m]) = (Subscript[x, m]).  Therefore, if the input signal is composed of harmonic functions and the data cover an integral number of periods, the Fourier series consists of spikes at each input frequency and its reflection.  However, if any period is truncated a broad background is produced because suppression of the background involves delicate cancellations which require the presence of full periods sampled uniformly.

General properties

Use the simple gaussian function

gaussian[t_, τ_, σ_]  :=  
      Exp[-(((t - τ))^2)/(2  (σ^2))]

for the following exercises.  When evaluating list plots, it will often be necessary to select the appropriate ranges and sometimes it is useful to connect the dots.  Often it is necessary to use PlotRange->All in order to identify the region of interest.  Also, make sure that your data are approximate real numbers instead of exact symbolic numbers.

Construct a table of 100 input data using τ->70 and σ->20.  Use ListPlot to view the data.

Evaluate the discrete Fourier transform of these data.  Compare the first element of the transform with the sum of the input data and explain the ratio between them.  [Hint: look at the symbolic expression for \(x\&~\_1\).]

Add elements 2 through 100 to a list with elements 100 to 2 in reversed order.  What happens to the imaginary parts?  What happens to the real parts?

Construct the inverse Fourier transform explicitly without using the built-in function.  Plot the inverse Fourier series for {t, 1, 200} and compare with the original function.  [Hint: construct a table and use ListPlot.]

How does the Fourier transform depend upon the width σ?  The option ScaleType->Symmetric is useful here.

Periodic functions

Consider the simple sinusoidal function with period λ.

f[λ_, t_] := Sin[((2  π)/λ) t]

Suppose that the period is λ->20.  Compare frequency spectra for observation times of T equal to 1, 2, 5, and 10 periods with 20 points sampled per period.  Interpret the result.  [Hint: remember that if the signal is real, the power spectrum must be symmetric.]

Again suppose that the period is λ->20 but now the observation period T is  2.25, 5.25, or 10.25 periods.  How is the spectrum affected by truncation of a cycle?  Show that a single missing point matters also.  What is the effect of increasing T?

Evaluate the Fourier series for the function

\(TraditionalForm\`g[
        t]\  = \ ∑\+\(p = 1\)\%10 p\^\(-2\)\ Sin[\(\(2  π\)\/λ\) p\ t]\)

where λ->20 and n -> 50  λ.  Interpret the result.  What is the effect of truncation?  Compare the inverse Fourier series with the input function.

Wave packets

Consider a simple wave packet constructed by modulating a sinusoidal signal with a gaussian envelope.

packet[t_, λ_, τ_, σ_] := 
      f[λ, t] gaussian[t, τ, σ]

Construct and plot 1000 input data using {λ->100, τ->500, σ->100}.

Evaluate and plot the power spectrum.  Interpret the result.

Study the dependence of the power spectrum upon σ.  What is the effect of decreasing σ?

Show that changing τ does not affect the power spectrum (provided that the function is not cut off) — the position of the peak affects the phase but not the intensity of the Fourier transform.

Noise

Now suppose that our instrument receives noise in addition to the signal of interest.  Although the correlation between noise signals at neighboring time intervals depends upon details of instrumentation and mechanism far beyond the scope of this study, it is of interest to consider the effect of the simplest type of noise.  White noise is characterized by a uniform power distribution and corresponds to random, uncorrelated, signals in each time interval.  Random signals usually follow a normal or Gaussian distribution.  It is useful to load a couple of packages.

Needs["Statistics`ContinuousDistributions`"];
Needs["Statistics`DataManipulation`"]

Suppose that the noise is characterized by a width σ = 0.1 compared to a unit signal.

NoiseSignal[σ_] := 
      Random[ NormalDistribution[0. , σ] ]

It is instructive to use a histogram of the noise signal to verify that our noise function is distributed normally.

histo = 
        BinCounts[ 
          Table[NoiseSignal[0.1], {1000}],  {-0.4, 0.4, 
            0.05}];
bins = Range[-0.375, 0.375, 0.05];
histodata = 
        Transpose[{bins, histo}] /. {x_, y_} -> {x, y, 
              0.05};

GeneralizedBarChart[histodata, PlotRange -> All, 
        Frame -> True ];

Plot the power spectrum for the noise signal alone.  How does the amplitude depend upon the number of samples?

Now suppose that the noise signal is added to a sinusoidal input signal.  What is the effect upon the power spectrum?  Why do long observation periods clean up the spectrum?

Sampling errors

Even if noise is suppressed, experimental observations are still subject to random errors.  Suppose that the random error is proportional to the size of the signal itself.

Construct 1000 perfect data for a wave packet using {λ->100, τ->500, σ->100}.

Apply random ±10% errors using the normal distribution.  Compare the imprecise data with the perfect data.

Plot the interesting range of the power spectrum and interpret the result.

Repeat previous two steps several times and comment.  How would averaging several measurements help?

Spectral analysis of van der Pol oscillator

As an example of an anharmonic periodic signal, consider the van der Pol oscillator whose phase-space trajectory is given by the following function:

trajectory[τ_, ε_, x0_, y0_, tmax_]:={x[τ], y[τ]}/.NDSolve[{x'[t]==y[t],y'[t]==-x[
                          t] + ε (1 - x[t]^2) y[t],x[0] == x0,y[0] == y0},{x, y},{t, 0, tmax},MaxSteps -> (10^5)][[1]]

It is useful to perform a Fourier analysis of the motion on the limit cycle.  Because we have already determined that for a given ε the system will reach the same limit cycle for any initial conditions, we choose fixed initial conditions.  The time interval {tmin,tmax,tstep} should be chosen to begin after the limit cycle has been reached and to include a large enough number of cycles with enough points per cycle to produce a stable power spectrum.

vdPolPower[ε_, tmin_, tmax_, tstep_, opts___Rule] := 
      Module[{sample, spectrum, optPowerSpectrum}, 
        optPowerSpectrum = FilterOptions[PowerSpectrum, opts]; 
        sample = 
          Table[Evaluate[trajectory[t, ε, 1. , 0. , 
                  tmax][[1]]], {t, 
              tmin, tmax, tstep}]; 
        spectrum = PowerSpectrum[sample, optPowerSpectrum]; 
        PlotPowerSpectrum[spectrum, opts, Frame -> True, 
          FrameLabel -> {"frequency", "power"}, 
          PlotLabel -> "van der Pol oscillator, ε=" <> 
              ToString[ε]]]

The figure below shows the spectrum for a large damping parameter.  We have used a large integration range and offset to ensure that the periodicity emerges clearly.

vdPolPower[10. , 100, 500, 0.1, 
        PlotRange -> {{1, 300}, {(10^(-6)), 10}}];

Notice that peaks occur at the odd harmonics.  Why?  [Hint: consider the Fourier transform of the differential equation.]  Observations of this kind can give insight into the dynamics responsible for observed motions.

Compare spectra for ε -> 0.1 and ε -> 10.  Interpret you results.

Spectral analysis of logistic map

Recall that the logistic map, (Subscript[x, n + 1]) = μ (Subscript[x, n])(1 - (Subscript[x, n])), is a population model for species with nonoverlapping generations.  The behavior of this model is strongly affected by the growth parameter μ.  For some ranges of μ the behavior is periodic in the sense that a sequence of points is visited in each cycle, while for other ranges the system is chaotic.  Fourier analysis of sequences {(Subscript[x, n])} shows strikingly different characteristics for these types of behavior.

The following function is designed to display a power spectrum for the logistic map after the initial transients have been allowed to decay.  The zero-frequency mode is omitted also.

logisticPower[μ_, n_, ndrop_, opts___Rule] := 
      Module[{sample}, 
        sample = 
          Drop[Fourier[ 
              Drop[NestList[μ  # (1 - #) &, 0.5, n + ndrop], 
                ndrop] ], 1]; 
        LogListPlot[ sample Conjugate[sample] , opts, 
          PlotLabel -> "μ=" <> ToString[μ]]]

We studied the behavior of the logistic map in logistic map.nb and cobweb.nb.  The next few examples illustrate the spectral properties of some of the characteristic behaviors.  Be sure to use enough points to sample the frequencies well and to drop enough starting points to suppress transients — some experimentation may be needed.

We found that for 1 < μ < 3 there is a stable fixed point.  Frequency analysis of a fixed point is not very interesting!

For 3 < μ < 3.45 there is a 2-cycle — display a spectrum for μ~3.2.  Why is it peaked in the center?

For 3.45 < μ < 3.54 there is a 4-cycle — display a spectrum for μ≈3.5.

Display a spectrum for μ ≈ 3.55 where there is an 8-cycle.

As μ continues to increase, period doublings become increasingly dense until chaos reigns.  Display a spectrum for μ=3.6.

However, order miraculously reappears near μ ≈ 3.83 when a 3-cycle becomes possible, but dissolves unto chaos again quite quickly.  Display spectra for μ = 3.83, 3.845, 3.85, and 3.86.

Within chaotic regions there are also some narrow slices of order.  Inspection of the logistic map figure suggest that μ=3.74 is one such value — verify with a frequency spectrum and determine the cycle length.  How far away is chaos from there?  Find another orderly slice and display its spectrum.

Some practical considerations

Normalization and phase conventions

Many conventions for the normalization and phase of Fourier series and for the normalization of the power spectrum exist in the technical literature.  Some authors prefer to use positive imaginary arguments of the exponential in the series for (Subscript[x, m]) so that \(TraditionalForm\`\(x\&~\)\_k\) is the complex conjugate of the present definition.  Others omit the (n^(-1/2)) factor in (Subscript[x, m]) so that \(TraditionalForm\`\(x\&~\)\_k\) has a factor of (n^(-1)), rather than the more symmetric choice of (n^(-1/2)) for both used by Mathematica.  Still others put a factor of (n^(-1)) in (Subscript[x, m]) while omitting (n^(-1/2)) from \(TraditionalForm\`\(x\&~\)\_k\).  Many employ frequencies \(TraditionalForm
      \`{ν\_k, \(-n\)\/\(2\ T\), \ n\/\(2  T\), \ 1\/T}\) instead of the Mathematica choice {(Subscript[ν, k]), 0, ((n - 1)/T), (1/T)}.  All possible choices of normalization, phase, and grid conventions appear with similar frequency — you just have to cope with it!

Similarly, there are several possible normalization of the power spectrum.  The most appropriate definition for applications to physical systems depends upon the physical meaning of the displacement and time variables and their relationship to the energy and power associated with a period oscillation.  For other applications the definition is more a matter of convenience.  The fundamental relationship between power in the time and frequency domains is Parseval's theorem, which for Fourier series states

\(TraditionalForm\`∑\+\(m = 1\)\%n|x\_m|\^2 = ∑\+\(k = 1\)\%n|x\&~\_k|\^2\)

using the symmetric normalization convention.  Recognizing that for a sinusoidal oscillation with unit amplitude the average of |(x^2)| over a period is 1/2, the sum on the left-hand side reduces to n/2 for a unit sinusoid.  Thus, if we employ units of power for which a unit sinusoid has average power ρ,  we should define \(TraditionalForm\`P\_k = \(2\/n\) ρ\ |x\&~\_k|\^2\).  Of course, in the absence of a physical relationship or after nondimensionalization of expressions, we normally choose ρ->1.

Fast Fourier Transform

Direct computation of Fourier series using naive implementation of the definitions given in the introduction is quite inefficient, scaling with (n^2).  You probably noticed that the built-in Mathematica function was much faster than the inverse transformation you performed in one of the first exercises in this notebook.  All modern Fourier transformations use some variation of the fast Fourier transform (FFT) algorithm, pioneered by Cooley and Tukey in the mid-1960's, which reduces the computation time to n Log[2, n] in the same units.  For large n the reduction factor

n Log[2, n]/(n^2)  /. n -> (10^6)  // N

can be critically important.  FFT implementations are generally most efficient when n -> (2^k) is an integral power of two and will usually subdivide a data set that does not satisfy this criterion so that the transforms of subsets can be computed and combined with maximum efficiency.

Construct a data set containing 8192 points of a sinusoidal function and then evaluate Timing[Fourier[data];].  [We recommend a sinusoidal function to ensure reasonable function values and to avoid complications due to extremely small or extremely large numbers.]  Next, use Drop to make a slightly shorter data set containing 8191 points and repeat the timing evaluation.  You will probably find that Mathematica requires approximately 300 times as much time to compute the Fourier transform of 8191 points as it needs for 8192 points!  Not only is 8191 not an integral power of two, it is prime, which defeats the FFT algorithm and makes it hardly more efficient than naive summation.  Finally, you will find that transforming 8190 points is about 150 times faster than 8191 points but remains 2 times slower than 8192 points.  Therefore, although Mathematica does not restrict the number of points, its Fourier transforms are most efficient for integral powers of two and it is often advantageous to pad shorter data lists with zeros.  Above all, avoid data sets with prime length.

Discretization

The discrete Fourier transform is often, if not usually, employed as an approximation to the continuous Fourier transform of nonperiodic functions.  In addition to problems associated with the implicit periodicity of Fourier series, the discretization of a continuous function distorts the Fourier transform and its power spectrum.  More sophisticated estimates of the power spectrum attempt to correct for discretization and edge effects using windowing functions and digital filters.  We will return to these issues after studying Fourier transforms and convolution techniques.

Aliasing

Parseval's theorem tells us that the total power in the frequency domain is equal to the total power in the time domain — in other words, that the Fourier transform conserves power.  If we know that our signal is bandwidth limited, meaning that all frequency components outside the sampling range vanish, then the total Fourier in the discrete Fourier series will be the same as that of the continuous Fourier transform.  For electronic signals bandwidth limitations can be imposed using electronic filters, but it is more difficult to impose such limitations on time series for dynamical systems.  However, if the input signal is not bandwidth limited, any power at frequencies greater than the Nyquist critical frequency, (Subscript[ν, c]) = (1/(2  T)), will appear in the discrete power spectrum at lower frequencies.  This artificial difference between the continuous and discrete power spectra is an error in the discrete approximation known as the aliasing problem because high frequencies are aliased to lower frequencies.  Often there is little one can do about the aliasing problem except to employ greater sampling rates for signals whose Fourier amplitudes remain appreciable near the Nyquist frequency.