# Using Mathematica to take Fourier transform of data

i'm trying to use Mathematica to take a discrete Fourier transform of 68 data points (0 to 67) to perform further analysis.

My data points, when plotted, look like this: Pretty straightforward.

Now, I need to take the DFT. These points are in a variable as

t={0, 0, 0.0147059, 0.0294118, 0.0441176, 0.0588235, 0.0735294,
0.0882353, 0.102941, 0.0441176, -0.0147059, -0.0735294, -0.191176,
-0.308824, -0.426471, -0.544118, -0.691176, -0.838235, -0.985294,
-1.13235, -1.27941, -1.42647, -1.57353, -1.72059, -1.86765, -1.01471,
-1.16176, -1.30882, -1.44118, -1.57353, -1.70588, -1.83824, -1.97059,
-2.10294, -2.23529, -2.36765, -2.5, -2.63235, -2.76471, -2.89706,
-3.01471, -3.13235, -3.25, -3.36765, -3.48529, -3.60294, -3.72059,
-3.83824, -3.95588, -4.07353, -4.19118, -4.30882, -4.42647, -6.54412,
-6.66176, -6.77941, -6.89706, -7.01471, -7.13235, -7.25, -7.36765,
-7.48529, -7.60294, -7.72059, -7.83824, -7.95588, -8.07353, -8.19118}


First question: should I have this in a vector of something like t={{x1,y1},{x2,y2},...}?

Because when the transform is taken, I'm not sure if mathematica will convert the time values (0-67) into the frequency domain. I do know that for this particular graph, the points should each be 2*Pi/N radians apart, where N is the fundamental period of this function (68).

Second Question: When the transform is taken, is the frequency data in radian format? So that I can plot the data from -Pi radians to Pi radians? I assume this has something to do with the FourierParameters[] function but I'm not completely sure what all the different combinations are. The documentation is not very good.

So to sum up, I guess, what is the best format for my input data to be in to take the Fourier transform, and then how is the best way to plot the output data (Abs value plot and complex phase plot)

I have put some basic information on numerical Fourier transforms here where you will find answers to all your questions. For your specific example we proceed as follows

ft = Fourier[t, FourierParameters -> {-1, -1}];
sr = 250;
ff = Table[(n - 1) sr/Length[ft], {n, Length[ft]}] // N;
ListLinePlot[Transpose[{ff, Abs[ft]}], PlotRange -> All,
Frame -> True, FrameLabel -> {"Frequency/Hz", "Absolute Value"}]
ListLinePlot[Transpose[{ff, Arg[ft]}], PlotRange -> All,
Frame -> True, FrameLabel -> {"Frequency/Hz", "Phase"}] I do not know your sample rate so I have assumed a sample rate (sr) of 250 points per second.

1. Do not include the abscissae in the input to Fourier.
2. If you want radians per second for the ordinates of the spectra then multiply my frequency values above (ff) by 2 Pi.

Hope that helps.

Edit. Frequency axis going from negative to positive values

The output of Fourier is periodic so going from negative frequencies to positive frequencies is a shift of the origin which is easily obtained using RotateRight. Care must be taken to ensure that the first point of the output, the zero frequency spectral value, ends up at the origin. The only tricky bit is that slightly different codes are needed depending on whether the number of points in the spectrum is even or odd. If the number of points is even then there is a value at SR/2. This point is missing for an odd number of points.

Below I do the general case where there can be any sample rate. In your case the sample rate is 1.

For an even number of points (which you have in your example) work as follows:

nn = Length@t;
ft = Fourier[t, FourierParameters -> {-1, -1}];
sr = 1.;
ff = Table[2 π (n - 1) sr/nn, {n, (-nn + 4)/2, (nn + 2)/2}];
ft1 = RotateRight[ft, (nn - 2)/2];
ListLinePlot[Transpose[{ff, Abs[ft1]}], PlotRange -> All,
Frame -> True, FrameLabel -> {"Frequency/Hz", "Absolute Value"}]
ListLinePlot[Transpose[{ff, Arg[ft1]}], PlotRange -> All,
Frame -> True, FrameLabel -> {"Frequency/Hz", "Phase"}] For an odd number of points the code is as follows:

Use the same data but drop the last point to give an odd number of points.

t1 = Most@t;
nn = Length@t1;
ft = Fourier[t1, FourierParameters -> {-1, -1}];
sr = 1.;
ff = Table[2 π (n - 1) sr/nn, {n, (-nn + 3)/2, (nn + 1)/2}];
ft1 = RotateRight[ft, (nn - 1)/2];
ListLinePlot[Transpose[{ff, Abs[ft1]}], PlotRange -> All,
Frame -> True, FrameLabel -> {"Frequency/Hz", "Absolute Value"}]
ListLinePlot[Transpose[{ff, Arg[ft1]}], PlotRange -> All,
Frame -> True, FrameLabel -> {"Frequency/Hz", "Phase"}] • Thank you very much for the prompt reply, very much appreciated. These look very similar to what I got on my own. I'm not sure what you mean by sampling rate -- like I said above, each point needs to occur at 2*Pi/N where N is the period but I'm not sure if that is what you are talking about. What is n in your code? In the third line? Nov 21, 2016 at 20:57
• The n in my code is just an index. When you say that the fundamental period is 68 you need to give us some units. Is that 68 points in one second or 68 seconds in one record length? In the former case it is a sample rate of 68 samples per second. In the latter it is a sample rate of 1 sample per second.
– Hugh
Nov 21, 2016 at 21:05
• Ahh I see. It is 68 samples in one period, each sample representing 1 second. Nov 21, 2016 at 21:07
• is this Fourier function anything like the FFT which shifts the first point of the transform to 0 frequency? I'm just wondering if I need to pad with 0s and all that to shift the transform. Nov 21, 2016 at 23:16
• The FFT is an algorithm for calculating the numerical Fourier transform. It requires the record length to be a power of 2 e.g. 4096. No such restrictions are required for Fourier here. Fourier will use the FFT if the record length is a power of 2. To force this to occur it one can increase the record length by padding zeros to reach a power of two. However this has unfortunate side effects and is probably not worth doing for your data.
– Hugh
Nov 22, 2016 at 10:49

Explanation

I've always found the Fourier function in Mathematica to be a bit lacking and non-intuitive. Like you seem to imply, a function which would take time-domain data in the form u={{t1,y1},{t2,y2},...} and return frequency domain data of the form v={{f1,Y1},{f2,Y2},...{fmax,Yn}} would be much more useful in many cases. So I've gone ahead and written my own that can simply be used as follows:

DiscreteFourier[u,fmax] (*Code provided at bottom of this answer*)


Here u is the list containing your time-domain points and fmax is the maximum frequency you want the output list to go up to.

Examples

For example, to compute the fourier transform of the function $$\sin\left( 2\pi 5 t\right)$$ from sampled data we could do something like this:

data = Table[{t, Sin[2 \[Pi] 5 t]}, {t, 0, 5, 0.001}];
ft = DiscreteFourier[data, 20];
ListPlot[ft /. c_Complex :> Abs[c], PlotRange -> All]


Which produces the following output Additionally, if you have data that isn't sampled at a constant sample rate, the function will automatically detect it and re-scale it using interpolation like this:

data = Table[{t, Sin[2 \[Pi] 5 t]}, {t, RandomReal[{0, 5}, 5/0.001]}];
ft = DiscreteFourier[data, 20];
ListPlot[ft /. c_Complex :> Abs[c], PlotRange -> All] Notice, however, that when the samples aren't evenly spaced you loose resolution in the frequency domain, despite having the same number of samples as the first example.

Code

Clear[DiscreteFourier];
DiscreteFourier[pts_, fmax_, opts:OptionsPattern[]] := Module[{\[Omega]max, pt, \[CapitalDelta]ts, \[CapitalDelta]t, u, n, b, \[Omega], ifun},
\[Omega]max = 2*Pi*fmax;
pt = SortBy[pts, First];
pt = pt /. {t_ /;  !ListQ[t], x_ /;  !ListQ[x]} :> {t - pts[][], x}; n = Length[pts];
\[CapitalDelta]ts = Transpose[pt][][[2 ;; -1]] - Transpose[pt][][[1 ;; -2]] /. t_NumberQ :> 0.0002; \[CapitalDelta]t = Mean[\[CapitalDelta]ts];
If[Max[\[CapitalDelta]ts - \[CapitalDelta]t] > OptionValue[Tolerance],
ifun = Interpolation[pt, InterpolationOrder -> OptionValue[InterpolationOrder]];
pt = Table[{((r - 1)/(n - 1))*pt[[-1]][], ifun[((r - 1)/(n - 1))*pt[[-1]][]]}, {r, 1, n}];
Message[DiscreteFourier::Rescale];
];
u = Switch[OptionValue[Method],
"LeftRiemann",
Table[pt[[r]][]*\[CapitalDelta]t, {r, 1, n - 1}],
"RightRiemann",
Table[pt[[r + 1]][]*\[CapitalDelta]t, {r, 1, n - 1}],
"Trapezoidal",
Table[((pt[[r + 1]][] + pt[[r]][])/2)*\[CapitalDelta]t, {r, 1, n - 1}],
_,
Message[DiscreteFourier::InvalidMethod, OptionValue[Method]];
Abort[];
];
b = (\[Omega]max*pt[[-1]][]*n)/(2*Pi*(n - 1)^2);
v = Fourier[u, FourierParameters -> {1, b}];
Table[{(1/(2*Pi))*((s - 1)/(n - 1))*\[Omega]max, v[[s]]}, {s, 1, n - 1}]
]
DiscreteFourier::InvalidMethod = "1 is not a valid method";
DiscreteFourier::Rescale = "The provided datapoints are not evenly spaced and so interpolation will be used to re-sample them at a constant sample rate. Note that this can introduce considerable error in some cases.";
Options[DiscreteFourier] = {Method -> "Trapezoidal", Tolerance -> 10^(-10), InterpolationOrder -> 1};


Code explained

First the function takes the data points and sorts them and shifts them so that they run from t=0 to tmax. Then it checks to see if the data points are evenly spaced in time by seeing if the difference between any 2 consecutive points is more than Tolerance away from the mean of all the other consecutive points. If the data is not evenly spaced it will use an interpolating function to re-sample them as even spacing is necessary for the FFT algorithm to work. I could probably get around this by using SparseArrays, but that's work for a future date. Anyway, after we have evenly spaced data, the function then re-arranges the points into a list u and a parameter b which are fed into Mathematica's Fourier function. Both u and b are found so that the summation preformed by Fourier[] mimics a Riemann sum. The specific type of Riemann sum being used can also be specified using Method->"Trapezoidal", Method->"LeftRiemann" or Method->"RightRiemann". Finally, we rescale the output of Fourier[] so that it returns a list of points of the form {f,Y(f)} where f is a frequency and Y(f) is the Fourier transform of the signal at that frequency.