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:
My question is: What do the X and Y axis represent in this result?
Mathematica Stack Exchange is a question and answer site for users of Wolfram Mathematica. It only takes a minute to sign up.
Sign up to join this communityIn 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:
My question is: What do the X and Y axis represent in this result?
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.
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]
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]
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]
Note how now there are values at all frequencies.
ListPlot[Transpose[{freqs, Abs[ft]}], PlotRange -> {{0, 30}, All}]
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]
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]
Let me know if you need more.