Numerical Fourier Transforms with Mathematica

Stephen Schoun

A tutorial is presented which explains how to use Mathematica to do Discrete Fourier Transforms.  Mathematica is a powerful tool for data analysis and can be ideal for manipulating a frequency spectrum, but using it to compute Fourier Transforms can be difficult.  All of the programming code necessary already exists within Mathematica, though not in one simple package.  Here we attempt to put everything together and explain the steps involved.

Continuous Fourier Transforms

For this tutorial, we use t as the time domain variable and ω as the frequency domain variable.  The Fourier Transform ℱ[g(t)] of a function g(t) will be defined with the convention:

ℱ[g(t)] ≡ (∫_(-∞)^∞)g(t) (e^(-i ω t))dt=g(ω)

The Inverse Fourier Transform (ℱ^(-1))[g(ω)] of a function g(ω) will then be defined with the convention:

(ℱ^(-1))[g(ω)]≡ (1/(2 π))(∫_(-∞)^∞)g(ω) (e^(i ω t))dω=g(t)

These functions exist in the Mathematica language as FourierTransform[] and InverseFourierTransform[]:

?FourierTransform

?InverseFourierTransform

However, Mathematica uses a different default convention {a, b} = {0, 1}:

ℱ[g(t)] ≡(Sqrt[((|b|)/(((2π))^(1-a)))])(∫_(-∞)^∞)g(t) (e^(ibωt))dt

(ℱ^(-1))[g(ω)]≡(Sqrt[((|b|)/(((2π))^(1+a)))])(∫_(-∞)^∞)g(ω) (e^(-ibωt))dω

We can implement our convention with the option:  FourierParameters -> {1, -1}

Define a gaussian envelope gt and calculate its Fourier Transform gω.

gt=(E^(-(t^2)));

gω=FourierTransform[gt,t,ω,FourierParameters->{1,-1}]

gtPlot=Plot[gt,{t,-3,3},PlotRange->Full,AxesLabel->{"t","g(t)"}];
gωPlot=Plot[gω,{ω,-5,5},AxesLabel->{"ω","g(ω)"}];
{gtPlot,gωPlot}

Discrete Fourier Transforms

Mathematica's Discrete Fourier Transform functions are called Fourier[] and InverseFourier[]:

?Fourier

?InverseFourier

Let's sample our gaussian envelope gt.  Choose n time steps of step size dn.  This gives a frequency step size of dω = ((2 π)/(n dn)).  Fourier[] works faster if n is a power of 2.

n=(2^9);dn=0.1;dω=((2 π )/(n dn));

This allows for a total sample range of (n dn) in time.

n dn//N

Create a list of data points in time.  Remember that we want to only have n number of sample points.  If we include t=0, then our list will need to have one fewer data points on one side of zero compared to the other side. As a short example, create a time list with n=8=(2^3) and step size dn=1:

TimeExample=Table[t,{t,((-(8 1))/2)+1,(((8 1))/2),1}]

We see that this has the correct number of points (8), and it is asymmetric with one less data point on the negative time side.  This result can be generalized for any n and dn; call this TimeList.  Make a similar list which samples the gaussian envelope; call this gtList.

TimeList=Table[N[t],{t,((-(n dn))/2)+dn,(((n dn))/2),dn}];

gtList=Table[gt,{t,((-(n dn))/2)+dn,(((n dn))/2),dn}];

In the frequency domain, we need to create a similar list of frequencies, now using the frequency spacing dω; call this omega.

omega=Table[N[ω],{ω,((-(n dω))/2)+dω,(((n dω))/2),dω}];

Fourier[] requires that the first element of the list be the zero time.  Cycle the elements of the list by using the RotateRight[] and RotateLeft[] functions.  Remember that our lists are asymmetric.

?RotateRight

?RotateLeft

Let's practice on TimeExample to make sure we cycle the correct number of positions.

RotateRight[TimeExample,8/2+1]

RotateLeft[RotateRight[TimeExample,8/2+1],8/2+1]

This result can be generalized for any n, to:  RotateRight[List, n / 2 + 1].

Take the Fourier Transform of gtList; call this gωList.  Use the same FourierParameters -> {1, -1}.  We need to normalize by the time spacing dn.

gωList=dn Fourier[RotateRight[gtList,n/2+1],FourierParameters->{1,-1}];

Fourier[] outputs a list which has the zero frequency as the first element.  To cycle this list back and match it with omega, we need to RotateLeft[].  Create a nested list of gωList to incorporate the frequency spacing; call it EωList2.  A 2 dimensional nested list can be created with the Transpose[] function:

?Transpose

gωList2=Transpose[{omega,RotateLeft[gωList,n/2+1]}];

Plot this result:

ListPlot[gωList2,AxesLabel->{"ω","g(ω)"}]

Check that we have normalized and shifted everything correctly by plotting the DFT gωList2 with the continuous function gω.  Combine multiple plots with the function Show[].

?Show

Show[Plot[gω,{ω,-5,5},AxesLabel->{"ω","g(ω)"},PlotStyle->Red,PlotRange->{{-5,5},Full}],ListPlot[gωList2,PlotStyle->Opacity[0.7]]]

Compute the Inverse Fourier Transform to make sure we get our original gaussian envelope; call this final gaussian gtfList and its 2 dimensional list gtfList2.  We need to normalize by (1/dn).

gtfList=(1/dn)InverseFourier[gωList,FourierParameters->{1,-1}];

InverseFourier[] outputs a list that has zero time as the first element.  To cycle back, use RotateLeft[]:

gtfList2=Transpose[{TimeList,RotateLeft[gtfList,n/2 +1]}];

Check that we have normalized and shifted everything correctly by plotting the Discrete Inverse Fourier Transform gtfList2 with the continuous function gt:

Show[Plot[gt,{t,-5,5},AxesLabel->{"t","g(t)"},PlotStyle->Red],ListPlot[gtfList2,PlotStyle->Opacity[0.7]]]

Getting More Out of Your DFT

In general, functions g(t) and g(ω) will be complex.  Fourier[] and InverseFourier[] can easily handle lists of complex numbers.  To find the real part use Re[]; to find the imaginary part use Im[]; to find the magnitude use Abs[]; to find the phase use Arg[].  But be careful, do not try to apply these functions to the 2 dimensional nested lists because, for instance, we don't want to take the absolute value of omega or TimeList.

?Re

?Im

?Abs

?Arg

In order to unwrap the phase, we have to write a program (SuperDagger[like]) Unwrap[].  Be careful to only use Unwrap[] after Rotating a list back to its proper order.

Unwrap[x_List]:=Module[{diff},diff=Drop[x,1]-Drop[x,-1];x-FoldList[Plus,0,2 π Round[diff/(2 π)]]]

Test Unwrap[] with a function that has a well known phase like g(t)=(e^(i t)):

PhaseExample=Table[Arg[(E^(I t))],{t,0,20,0.1}];ListPlot[{PhaseExample,Unwrap[PhaseExample]}/π,Joined->True,AxesLabel->{"t",(("phase")/π)}, PlotLabel->"Phase Wrapped and Unwrapped"]

Example:  Linearly Chirped Gaussian

For one last example, let's consider a linearly chirped gaussian envelope gct and its Fourier Transform gcω :

gct=(E^(-(1-I) (t^2)));

gctList=Table[gct,{t,((-(n dn))/2)+dn,(((n dn))/2),dn}];

gcωList=dn Fourier[RotateRight[gctList,n/2+1],FourierParameters->{1,-1}];

Calculate the Real and Imaginary parts of gcω :

gcωList2Real=Transpose[{omega,Re[RotateLeft[gcωList,n/2+1]]}];

gcωList2Imaginary=Transpose[{omega,Im[RotateLeft[gcωList,n/2+1]]}];

ListPlot[{gcωList2Real,gcωList2Imaginary},Joined->True,PlotRange->{{-10,10},Full},AxesLabel->{"ω","Real & Imaginary gc(ω)"}]

Calculate the spectral phase of gcω :

gcωList2Phase=Transpose[{omega,Unwrap[Arg[RotateLeft[gcωList,n/2+1]]]/(2 π)}];

ListPlot[gcωList2Phase,Joined->True,PlotRange->Full,AxesLabel->{"ω","\!\(\*FractionBox[\(phase\), \(2\\ π\)]\)"}]

Other Considerations

It may be that sometime we will have to import data from an external file.  To do this, click on the Insert menu with the mouse, scroll down and select file path.  Find the file we want and hit open.  Name this string something like MyFile and use Import[] to allow Mathematica to read the file.  For example: MyData = Import[MyFile, "table"]  imports a file named MyFile and saves the data to a table named MyData.

?Import

When importing data, we are going to have to discover the number of data points n and the time spacing dn.  The data might be symmetric and have no t=0 point instead of asymmetric around t=0.  In this case we are going to have to cycle our lists differently with RotateRight[] and RotateLeft[].  We are also going to have to Table our TimeList and omega differently to account for the symmetry.

If we want to use a different convention for the Fourier Transform and Inverse Fourier Transform, we should run through this tutorial again, making sure that no spacing or normalization needs to be changed.  The safest procedure is always to check the Discrete Fourier Transform against the Continuous Fourier Transform.

†  Dick, Samuel, Alfred Riddle, and Douglas Stein.  Mathematica in the Laboratory.  Cambridge:  Cambridge University Press, 1997.  ISBN 0521499062