IMSE Chapter 8

IMSE-08(Fourier).nb
Fourier Analysis and Power Spectra
RCH 1/05/02
revised 6/27/02 - rjerwin

This notebook explores the calculation of Fourier power spectra using Mathematica.

It is often useful to characterize experimental data by calculating the Fourier power spectrum of that data.  If the data are the time series of some variable, then the Fourier power spectrum tells us how the signal can be decomposed into signals of various frequencies.  If the data represent the spatial variation of some variable, the Fourier power spectrum gives us information about the so-called spatial frequencies (wavevectors) that make up the original data.

I assume that you are familiar with the basics of Fourier analysis and Fourier power spectra.  But, for the sake of completeness, I give a brief review here.

Fourier's theorem states that any periodic function with period T can be written as an infinite series of terms, each associated with a frequency (Subscript[ω, n])which is an integer multiple of (Subscript[ω, o])= 2 π/T.  In more formal terms, we have

			f(t) = (∑_(n  = -∞)^(+∞))(Subscript[c, n])(e^(in(Subscript[ω, o]) t)),
			
where we have used a complex exponential form for the time dependence.  (It is possible to use sines and cosines as well).
The Fourier power spectrum is defined to be

			P((Subscript[ω, n])) = |(Subscript[c, n])\( | \^\(\(\ \)\(2\)\)\).

For non-periodic functions, the so-called Fourier transform is used, which is based on a continuous parameter ω

			g(t)=(1/(Sqrt[2  π]))(∫_(-∞)^(+∞))c(ω)(e^(iω t))dω.
			
The power spectrum is then a continuous function of ω:

			P(ω) = \(TraditionalForm\`\(\(|\)\(c(ω)\)\( | \^2\)\)\).

Mathematica uses the command Fourier[dataList] to compute the amplitude of the Fourier coefficients from the list of data given in dataList.  The result of Fourier[dataList] is itself a list of the Fourier coefficients (Subscript[c, n]), which, in general, are complex numbers.  To get the power spectrum, we need to use (Abs[Fourier[dataList]])^2.

Let's generate some data to try out these notions.

ndata = 
        1000;  (*The number of data points to generate.*) 
dataList = 
        Table[Sin[2*π*k/120]*Cos[2*π*k/63], {k, 0, 
            ndata}];
ListPlot[dataList, PlotJoined -> True];

fpsList = (Abs[
              Fourier[dataList]])^2;
ListPlot[fpsList, AxesLabel -> {"frequency", "power"}, 
        PlotRange -> All, PlotJoined -> True];

There are several features to be noted in the plot of the power spectrum.  First, all of the plotted points seem to be scrunched up near very low frequencies and very high frequencies.  Second, the high frequency part appears to be the mirror image of the low frequency part.  In our case, the low frequencies reflect the original data, which repeats itself only a few times during the time interval being plotted.

To understand the mirror image part, we need to know some of the conventions that Mathematica uses to computer the Fourier coefficients.  Mathematica assumes that the list of data span one unit of time - a second, for example.  If we use ListPlot to give us a graph of the power spectrum, the frequency axis (the abscissa) will run from 0 to the number of data (ndata) in the list.  You will observe, however, that the spectrum from 0 to N/2 is just the mirror image of the spectrum from ndata/2 + 1 to ndata.  This is an artifact of the numerical procedures used to calculate the spectrum.  The mirror image also reflects (no pun) the fact that the highest possible frequency that can be reproduced in a sequence of ndata data points per second is 2/ndata.  This result is known as the Nyquist criterion.  Simple picture: the highest frequency representable corresponds to an alternation between high and low values for successive data points.  For ndata points, we get ndata/2 as the largest number of alternations that can occur.  Hence, the highest frequency that can be "captured" with ndata data points per unit time is 2/ndata.

Here we restrict the range to the frequencies from 0 to ndata/100.

ListPlot[fpsList, AxesLabel -> {"frequency", "power"}, 
        PlotRange -> {{0, ndata/10}, All}, 
        PlotJoined -> True];

If you want to scale the frequency axis to match the actual time scale of the data, you will need to create another list of frequency values.  We denote the time interval between successive data points as  (Subscript[t, o]).  Then the frequency scale should run from 0 to 1/(Subscript[t, o]) in steps of 1/(N(Subscript[t, o])).

The following input cell generates the list of frequencies, combines the frequency list and the Fourier coefficients squared list and then replots the Fourier power spectrum.  We have assumed that the time interval between data points is 0.01 second.

to  =  0.01;
step  =  1/(ndata*to);
freqList = Table[k, {k, 0, 1/to, 1/(ndata*to)}];

combinedList = 
        Transpose[{freqList, fpsList}];
\(\(ListPlot[combinedList, 
        AxesLabel -> {\*""\!\(sec\^\(-1\)\)"", "power"}, 
        PlotRange -> {{0, 1/\((to*10)\)}, All}, 
        PlotJoined -> True];\)\)

There appear to be two main frequency peaks, one between 0.5 and 1.0 (sec^(-1)), and another one around 2.3 (sec^(-1)).  Now, the first graph we plotted does look like a mess, and so you're probably wondering why the spectrum only show two peaks.  If you look closely, however, you'll notice that the distances between peaks fall into two categories: data point intervals around 150, and much smaller intervals around 40.  (These are only approximate values, remember.)  Since the time interval between data points is 0.01 seconds, these come out to 0.7 (sec^(-1)) and 2.0 (sec^(-1)), respectively, and match our frequency peaks pretty well!

But don't get too smug, not just yet.  This example was an easy one.  First of all, despite the messiness of the orginal graph, we were still dealing with a nice regular signal - that is, there was no noisy component to add a lot of distracting information into the waveform and the Fourier spectrum.  Second, since the mean value of dataList is very close to zero (try it), we don't need to remove any "DC component."  By DC component, I mean a non-zero mean value, so that Fourier now includes a periodic component of infinite period.  This then shows up as a spike at zero frequency.

Standard practice is to remove the zero frequency spike by removing the average value from the waveform.  In the following example, I've simply squared the original equation.  The mean value is no longer close to zero.

dataList = 
        Table[(((Sin[2*π*k/120]*Cos[2*π*k/63]))^2), {k, 0, 
            ndata}];
ListPlot[dataList, PlotJoined -> True];

fpsList = (Abs[
              Fourier[dataList]])^2;
to  =  0.01;
step  =  1/(ndata*to);
freqList = Table[k, {k, 0, 1/to, 1/(ndata*to)}];

combinedList = 
        Transpose[{freqList, fpsList}];
\(\(ListPlot[combinedList, 
        AxesLabel -> {\*""\!\(sec\^\(-1\)\)"", "power"}, 
        PlotRange -> {{0, 1/\((to*10)\)}, All}, 
        PlotJoined -> True];\)\)

As you can see, there's a sharp jump in the power spectrum at zero frequency, which you can easily zoom in on for better view.  Let's now remove the mean value.

<< Statistics`DescriptiveStatistics`

dataminusmean = 
        dataList - N[Mean[dataList]];
fpsminusmean = (Abs[
              Fourier[dataminusmean]])^2;
combinedmm = Transpose[{freqList, fpsminusmean}];

\(\(ListPlot[combinedmm, 
        AxesLabel -> {\*""\!\(sec\^\(-1\)\)"", "power"}, 
        PlotRange -> {{0, 1/\((to*10)\)}, All}, 
        PlotJoined -> True];\)\)

Surprise!  We have the same power spectrum as before, but without the non-zero mean peak.