32
$\begingroup$

In Wolfram Mathematica the function Fourier has the following declaration

Fourier[list]

And after a list is given to the function, a simple Plot gives the following result: Graphic

My question is: What do the X and Y axis represent in this result?

$\endgroup$
4

1 Answer 1

78
+50
$\begingroup$

It looks like you need some basic facts on numerical Fourier transforms. I am going to assume that you are going from the time domain to the frequency domain.

  1. The number of points in the time domain equals the number of points in the frequency domain.
  2. As the data is sampled in the time domain i.e. is a set of equally spaced points, then in the frequency domain the spectrum is periodic. Only one period is calculated.
  3. The Fourier representation of a time history gives values at positive and negative frequencies. The convention in Fourier analysis is to give the positive frequencies first and then negative. Following from 2 this is a quaint but valid representation.
  4. If in the time domain you have a cosine function that has a whole number of cycles, say n cycles, then there will be a single point at the n+1 value in the frequency domain and another value at n+1 from the end of the spectrum. All other points will be zero. These correspond to the positive and negative frequency values.
  5. For real functions in the time domain the real part of the Fourier transform is an even function and the imaginary part an odd function.
  6. The first point in the spectrum is the zero frequency value (the D.C. value).
  7. If in the time domain you have a sample rate of SR then in the frequency domain the points along the x axis go from zero to one less than the sample rate; i.e. the frequency of the nth point is f = (n-1) SR/N where N is the number of points. Here f is in Hz and the sample rate in samples per second.
  8. The ordinates of the Fourier transform are scaled in various ways but a basic theorem is that there is a scaling such that the mean square value in the time domain equals the sum of squared values in the frequency domain (Parseval's theorem). Actually as the values in the frequency domain are complex it is the square of the modulus values that must be summed.

Here are some examples

nn = 1000; (* Number of points *)
sr = 500; (* Sample rate *)
dt = 1/sr;  (* Time increment *)
df = sr/nn; (* Frequency increment *)
f1 = 30/(dt nn); (* Example frequency at 15 Hz with exactly 30 periods \
                     in time interval *)
p = 1/f1; (* Period of example frequency *)
a = 5; (* Amplitude of cosine wave *)

data = Table[N[{t, a Cos[2 \[Pi] f1 t]}], {t, 0, dt (nn - 1), dt}];
ListLinePlot[data]

Mathematica graphics

To set a desired scaling use the FourierParameters option. This is my favorite.

ft = Fourier[data[[All, 2]], FourierParameters -> {-1, -1}];
ListPlot[Abs[ft], PlotRange -> All]

Mathematica graphics

Examine the points at which we expect values at the 31st and 971th values.

ft[[{30, 31, 32}]]

(* {0., 2.5, 0.}  *)

ft[[{nn - 30, nn - 30 + 1, nn - 30 + 2}]]

(* {0., 2.5, 0.} *)

With my choice of Fourier parameters the height of each peak is half the amplitude. This comes from combining two values of Sqrt[2]. Examine the mean square value (time domain) and total square value (frequency domain).

{data[[All, 2]].data[[All, 2]]/nn, ft.Conjugate[ft]}
(* {12.5, 12.5}  *)

Start again with a frequency that does not have an exact number of periods in the time interval.

nn = 1000; (* Number of points *)
sr = 500; (* Sample rate *)
dt = 1/sr;  (* Time increment *)
df = sr/nn; (* Frequency increment *)
f1 = 14.314; (* Example frequency  *)
p = 1/f1; (* Period of example frequency *)
a = 5; (* Amplitude of cosine wave *)

data = Table[N[{t, a Cos[2 \[Pi] f1 t]}], {t, 0, dt (nn - 1), dt}];
ft = Fourier[data[[All, 2]], FourierParameters -> {-1, -1}];

If you wish to make a frequency axis then you can use the frequency increment given in point 5

freqs = Table[(n - 1) sr/nn, {n, nn}];
ListPlot[Transpose[{freqs, Abs[ft]}], PlotRange -> All]

Mathematica graphics

Note how now there are values at all frequencies.

ListPlot[Transpose[{freqs, Abs[ft]}], PlotRange -> {{0, 30}, All}]

Mathematica graphics

Again look at the scaling of the ordinates.

{data[[All, 2]].data[[All, 2]]/nn, ft.Conjugate[ft]}
(* {12.5408, 12.5408 + 0. I} *)

Keep working at simple examples like this and you will get the idea of numerical Fourier transforms.

Addendum

I have been asked "...what do each of the axes represent with standard Fourier in Mathematica". I will attempt to add some comments.

Fourier takes a list of numbers and outputs a list of numbers. In the write-up above I interpret the input as a time history and the output as a frequency spectra. Other interpretations are possible. The input could be a list of displacements and then the output would be a list of wavenumbers. You must decide what the input list means. When calculating the Fourier transform, Mathematica does not need to know the meaning of your input. The key idea is given in point 4 above; a cosine function that fits a whole number of cycles into the input list will produce two non-zero points in the output. If you use Listplot then the x-axis (abscissae) will just be the point number with the point number going from 1 to the number of points in the list. Point 7 above and the examples show how to add abscissae so that the axis is scaled. The scaling of the ordinates (y-axis) depends on what FourierParameters you have used. Point 8 above discusses this.

Putting in a frequency axis

Several options are available for constructing a frequency axis. Possibilities include a simple frequency axis that goes from 0 to the sample rate (less one increment). This is fine for the first half of the data but can be misleading for the second half (see point 3 above). Consequently it is common to remove the second half of the spectrum since typically it is uninformative being just the mirror image of the first half. Alternatively, you may wish to have negative and positive values properly represented. In this case you have to have to do slightly different things depending on whether the data is even or odd in length. Here we are going to only consider working in Hz. If you wish to work in radians per second then you will have to introduce a 2 Pi.

Having added the frequency axis there is then the issue of how to plot. As the data from Fourier is complex, it is necessary to convert it to absolute values before plotting. You may also wish to have a separate plot of phase. Here are some useful modules for performing these functions.

 ClearAll[insertFrequencies];
    insertFrequencies::usage = 
      "insertFrequencies[fd, sr] adds frequency values to a Fourier \
    spectrum. Here fd is the output from Fourier and sr is the sample \
    rate. Note that in the second half of the spectrum, containing the \
    negative frequencies, the frequency values will not be negative.   ";
    insertFrequencies[fd_, sr_] := Module[{nn},
      nn = Length[fd];
      Transpose[{Table[(n - 1) sr/nn, {n, nn}], fd}]
      ]
    
    ClearAll[toFreqMod];
    toFreqMod::usage = 
      "toFreqMod[fdata]  converts frequency data to absolute values. \
    Iinput fdata should be {{f1,y1},{f2,y2}...} where f is the frequency \
    and y are complex values. Output is {{f1,Abs[y1}},{f2,Abs[y2}},...}";
    toFreqMod[frf_] := Transpose[{frf[[All, 1]], Abs[frf[[All, 2]]]}]
    
ClearAll[insertNegativeFrequencies];
    insertNegativeFrequencies::usage = 
      "insertNegativeFrequencies[fd, sr] converts Fourier output to \
    spectra with negative and positive frequencies. fd is output from \
    Fourier and sr is the sample rate. The output is in the form {{f1, \
    y1}, {f2, y2}...}";
    insertNegativeFrequencies[fd_, sr_] := Module[{nn},
      nn = Length@fd;
      If[EvenQ[nn], 
       Transpose[{Table[sr/nn n, {n, -(nn/2) + 1, nn/2}], 
         RotateRight[fd, nn/2 - 1]}], 
       Transpose[{Table[sr/nn n, {n, -((nn - 1)/2), (nn - 1)/2}], 
         RotateRight[fd, (nn - 1)/2]}]]
      ]

Here is an example of the workflow, this time with data that has a mean value.

nn = 1000; (* Number of points *)
sr = 500; (* Sample rate *)
dt = 1/sr;  (* Time increment *)
df = sr/nn; (* Frequency increment *)
f1 = 30/(dt nn); (* Example frequency at 15 Hz with exactly 30 periods \
in time interval *)
p = 1/f1; (* Period of example frequency *)
a = 5; (* Amplitude of cosine wave *)
a0 = 1; (* Mean value *)
data = Table[
   N[{t, a0 + a Cos[2 π f1 t]}], {t, 0, dt (nn - 1), dt}];
ListLinePlot[data]

Mathematica graphics

To calculate the spectrum and add a frequency axis

b = Fourier[data[[All, 2]], FourierParameters -> {-1, -1}];
fd1 = insertFrequencies[b, sr];
fd2 = insertNegativeFrequencies[b, sr];
fd1a = toFreqMod[fd1];
fd2a = toFreqMod[fd2];
ListLinePlot[fd1a, PlotRange -> All]
ListLinePlot[fd2a, PlotRange -> All]

Mathematica graphics Mathematica graphics

Let me know if you need more.

$\endgroup$
13
  • $\begingroup$ Hugh, this is an excellent introductory writeup. Thank you! (+1) $\endgroup$
    – MarcoB
    Jun 5, 2015 at 9:22
  • $\begingroup$ @MarcoB Thank you . There are frequent requests for basic information on Fourier. I thought it would be useful to put together the above list of points. I am happy to expand further if there is a demand. $\endgroup$
    – Hugh
    Jun 5, 2015 at 10:19
  • $\begingroup$ @Hugh This is really an excellent piece of work. Very relevant to something I am working on. Thanks $\endgroup$
    – user27119
    Feb 21, 2019 at 13:39
  • $\begingroup$ @QuantumPenguin Thanks for your comments. Let me know if you need more. $\endgroup$
    – Hugh
    Feb 21, 2019 at 16:59
  • 1
    $\begingroup$ @QuantumPenguin For references some of the differences between Fourier and FourerTransform are considered here. $\endgroup$
    – Hugh
    Apr 29, 2019 at 11:56

Not the answer you're looking for? Browse other questions tagged or ask your own question.