Home | Mathematics | * Maxima |     Share This Page
Fourier Analysis
P. Lutus Message Page

Copyright © 2007, P. Lutus

Overview | The Frequency -> Time Task | The Time -> Frequency Task
Conclusion

(double-click any word to see its definition)

 
Overview

Fourier analysis is a fascinating activity. It deals with the essential properties of periodic waveforms of all kinds, and it can be used to find signals lost in apparently overwhelming noise. As just one example of its usefulness, if SETI (the Search for Extraterrestrial Intelligence) should ever detect an alien signal, that discovery will be made using Fourier analysis.

It's important to say that, even though this article deals with simple waveforms, Fourier analysis is by no means limited to these classic examples — it can analyze and process images, it can efficiently compress images and video streams, and it can assist in visual pattern recognition, where a complex pattern may be efficiently and concisely described using a set of Fourier terms.

One of the principles of Fourier analysis is that any imaginable waveform can be constructed out of a carefully chosen set of sinewave components, assembled in a particular way (the frequency -> time task). And conversely, any complex periodic signal can be broken down into a series of sinewave components for analysis (the time -> frequency task). And most important, the described tasks are reciprocal operations — in the same way that integration and differentiation are reciprocal operations in Calculus, encoding and decoding signals are reciprocal operations in Fourier analysis.

Using that notion as a guide, this page has a section in which we will construct complex signals using sinewave components, and another in which we will decode complex signals into their components.

The Frequency -> Time Task
This section deals with the task of creating time-domain signals out of individual sinewave components in the frequency domain.

First, download this Fourier example file and load it into Maxima. It contains a set of example functions that show aspects of Fourier analysis.

In this section we'll be using Maxima to build complex waveforms out of simple components. To reiterate, any signal can be built out of individual sinewave components, arranged in particular ways. Remember this the next time you're listening to your favorite music — in principle, it can be created out of a mathematical description consisting only of sinewave frequencies, amplitudes and phases.

Let's begin with a table of classic wave equations and corresponding graphs:

Name Equation Graph
Sine
Square
Triangle
Sawtooth
Rectified Sine

Obviously these graphs weren't created by summing terms from 1 to ∞ as indicated by the equations in the left column. In practice one chooses a summation number that balances display quality with a reasonable running time. But by choosing different upper bounds for these summations, we can learn something important about signal generation.

Make sure you have loaded this Fourier examples file into Maxima, then enter this at the keyboard:
fpsq(1);
This will accumulate one element from the square-wave equation, which, if you think about it, must be a sine wave (because all these generators use individual sine waves as their building blocks). Now change your entry to this:
fpsq(2);

Notice the difference. Now change the number, go as high as you like, considering the issue of running time, and see how the plotted wave changes. You will begin to see a pronounced "ringing" transient at the switching points in the waveform (like the graph on this page), and the ringing will become more severe as the number of accumulations increases. There is a longstanding debate among mathematicians as to whether this ringing is an artifact of our imperfect generation methods or a real property of square waves.

In case it isn't obvious from the equation listed above, the square wave is created by accumulating the odd multiples of a specified frequency — in this way (where x = 2 π f t):

y = (4/π) (sin(x) + 1/3 sin(3x) + 1/5 sin(5x) + 1/7 sin(7x) ...)

These sine terms are combined to create the time-varying waveform in the display. And if the opposite operation were to be performed on the result (something called a Fourier transform), these individual harmonic elements would reappear. It is important to reemphasize that waveform generation, and the Fourier transform, are reciprocal operations. You can use frequency components to generate a waveform in the time domain, then transform the result back to the frequency domain and recover what you started with. This reciprocal relationship is to Fourier analysis what the Fundamental Theorem of Calculus (the idea that integration and derivation are reciprocal operations) is to Calculus.

The Fourier examples file contains a set of waveform generating functions, one for each kind of waveform in the above table. Here are their names:

  • fpsq(n); (square)
  • fptri(n); (triangle)
  • fpsaw(n); (sawtooth)
  • fprec(n); (rectified sine)
Play with these functions — each will create a graph that results from n summations of the generating function, and you choose the value for n.

The Fourier examples file also has functions to create multiple plots for each choice, with results that look like this:

Function Graph
fsqmult(8);
fsawmult(8);

Here is a list of the multiple-plot functions:

  • fsqmult(n); (square)
  • ftrimult(n); (triangle)
  • fsawmult(n); (sawtooth)
  • frecmult(n); (rectified sine)

These functions take much more time to generate their plots (because they must perform a summation for each individual trace in the plot), so don't enter large numbers unless you are a very patient person.

The Time -> Frequency Task

Converting from the time domain to the frequency domain is a bit more labor-intensive than the task of generating complex waveforms, the topic covered above. The general name for this conversion is "Fourier Transform", and because of its usefulness, much thought and ingenuity have been expended on this task.

In the 1960s, an apparently new and much more efficient method for this conversion, now called the "Fast Fourier Transform" (FFT), was devised by J. W. Cooley and John Tukey. I say "apparently" because it turns out the basic idea behind the FFT was developed by Carl Friedrich Gauss in 1805 and was rediscovered by others in various limited forms in the intervening years.

Maxima has an FFT package available, made accessible by loading a library named "fft". Here is an outline of the steps needed to perform an FFT in Maxima:

Code
Explanation
load(fft); Load the required library.
n:2048; Establish an array size, which must be:
  • A power of 2 (512,1024,2048, etc.).
  • Twice the size of the highest frequency of interest in Hertz, e.g. to resolve frequencies between 0 and N Hertz, the array size must be the next higher power of 2 above N * 2.
array(ra,float,n-1); Declare a real-value array.
array(ia,float,n-1); Declare an imaginary-value array.
for i:0 thru n-1 do ( Fill the array with the desired data.
   ra[i]:my_function(i), This is a user-defined function that is able to create suitable data.
   ia[i]:0.0 Must set the imaginary array to zero if there are no data.
);
fft(ra,ia); Perform the FFT conversion in-place (the original array data are overwritten).
recttopolar(ra,ia); Optional: convert complex a+bi form to polar magnitude and angle form.
lira:makelist(ra[i],i,0,(n/2)); For plotting, create a list out of the first half of the data array.
xIndexList:makelist(i,i,0,length(lira)); Create an index list required by the plotting routine.
wxplot2d([discrete, xIndexList,lira], [gnuplot_preamble, "unset key;set xrange [0:1024];set xlabel 'Frequency';set ylabel 'Amplitude';set title 'Amplitude Modulation — Frequency Domain';"]); Create an inline plot (one that is integrated into wxMaxima's display).

There are example functions in the Fourier examples file that show a typical conversion:

Function
Graph
Comment
fpt();

This time-domain plot shows a waveform that was once very familiar to electrical engineers: an amplitude-modulated (AM) radio signal. In the AM scheme, information is added to a radio carrier wave by changing its amplitude. Speaking mathematically, this modulation method creates sidebands whose distance from the carrier frequency is equal to the modulation frequency. Typically the modulation of a real-world AM radio signal is complex and broad. In this example, we're modulating the carrier with a simple sinewave.

This plot is the input data for the FFT. If you have loaded the Fourier examples file into Maxima, and if you enter the function name at the left, you will get this plot.

fpf();

Here is the frequency-domain equivalent of the above time-domain plot (with some random noise added for a touch of realism). The central peak is the carrier frequency, and there are two sidebands, one above and one below the carrier frequency. The sideband frequencies are located at c+m and c-m, where c = carrier frequency and m = modulation frequency.

This plot is the output data from the FFT. If you have loaded the Fourier examples file into Maxima, and if you enter the function name at the left, you will get this plot.

There is one more experiment in the Fourier examples file. It is meant to clearly show the equivalence of time-domain and frequency-domain representations of spectral data. In this experiment, we will:

  • Create a square wave signal in the time domain using the methods established above,
  • Convert the time-domain data back to the frequency domain using the FFT,
  • Compare the resulting harmonic components to the individual terms in the square wave generator equation.

Actually, this process is automated for you by two functions that are part of the the Fourier examples file. They are:

Function
Comment
fph();

Predicts (and prints a list of) frequencies and amplitudes for square-wave harmonics based on the mathematical identity:

y = (4/π) (sin(x) + 1/3 sin(3x) + 1/5 sin(5x) + 1/7 sin(7x) ...)
fch();

Generates a square-wave data table in the time domain, converts it to the frequency domain using the FFT, then prints out the frequencies and amplitudes of the resulting harmonic components.

Make sure you have loaded the Fourier examples, then enter the two function names listed at the left above. One ("fph()") will predict the expected harmonic frequencies and amplitudes based on the mathematical identity of a square wave, the other ("fch()") will use the FFT to convert a time-domain square-wave data list into its harmonic components, then search the results and print out the harmonics it finds.

If everything is in order on your system and there is nothing wrong with your Maxima installation, the predicted and result data lists should be nearly identical. Because the FFT is a numerical procedure performed on machines with finite floating-pint resolution and a finite array size, all FFT results are approximate.

I want to emphasize that the numerical, approximate nature of the FFT doesn't imply that all Fourier transforms are approximate. The FFT is a specific, very fast, embodiment of a well-understood mathematical transformation that can be carried out analytically in some circumstances. What I am saying here is ... don't confuse the FFT with the field itself.

Conclusion

Apart from being rather intriguing in an aesthetic sense, Fourier analysis is a useful thing to know in the modern world. There are applications of this field in virtually all scientific and technological disciplines.

I want to reemphasize one point about Fourier analysis, and it is the single fact essential to know about this field. If my readers have forgotten all the content of this page a week from now, I hope they remember this:

Waveform creation in the time domain, and harmonic analysis in the frequency domain, are reciprocal operations. One can take a list of harmonic components and use them to create a time-domain waveform, then one can carry out a Fourier transform on the time-domain waveform to recapture the original harmonic components. It turns out that these two representations are equivalent and interchangeable.

 

Home | Mathematics | * Maxima |     Share This Page