# Numerical Fourier transform of a complicated function

Say I have a function $f(x)$ that is given explicitly in its functional form, and I want to find its Fourier transform[1]. If $f$ is too complicated to have an analytic expression for $\hat f(k)$, how do I obtain it numerically?

The naive and stupid way, which I currently use, is evaluating the integral for many $k$'s and then interpolating between them. This is outrageously inefficient, and I'm sure that FFT should be used, but I can't seem to find a good recipe for that. I keep getting the frequencies and the normalisation wrong.

This task seems to be common enough so that advanced software like Mathematica would have a built-in function for that. Any ideas?

[1] Being a physicist, I use $\hat f(k)\equiv \int e^{ikx}f(x)dx$ as the definition of FourierTransform, but you can stick in $2\pi$ wherever you feel like if you're used to a different one.

• FFT makes sense if your function has a finite support, or most of its energy is in a clearly defined interval in time. Is that your case? Also, you need to know beforehand some maximum frequency component (or again, a frequency after which your signal barely has any energy left). If you can get a list of samples of the entire support of the signal, sampled at more than twice the maximum frequency component (Nyquist theorem) then FFT is the tool. – Rojo Feb 13 '12 at 22:52
• @Rojo I think the question is: I have f(x) and can evaluate it for any x. Assume that it vanishes outside some range, say $-L<x<L$. How do I produce a list, do something to it and end up with a discretisation of its fourier transform? – acl Feb 13 '12 at 23:01
• @Rojo This is indeed the case. My function decays exponentially for large $|x|$, and is smooth (exponentially small energy at large $k$)- looks like a Gaussian, more or less. As I stated in the question - I can get a list of samples wherever I want, as I have the explicit functional form. – yohbs Feb 13 '12 at 23:01
• @yohbs, a function that does that would be useful to me too, I'll give it a shot and post it. However, I'm used to another definition of the Fourier transform. I'll try to make it work for both – Rojo Feb 14 '12 at 0:05

There is the function NFourierTransform[] (as well as NInverseFourierTransform[]) implemented in the package FourierSeries. The function, as with the related kernel functions, takes a FourierParameters option so you can adjust computations to your preferred normalization as needed. For your specific normalization, you apparently want the setting FourierParameters -> {1, 1}.

Now, since NFourierTransform[] internally uses NIntegrate[], the function also takes options that can be passed to NIntegrate[]; you can thus change Methods as seen fit. I will in particular recommend that you look into the Method choices "ExtrapolatingOscillatory" (Longman's method), "DoubleExponentialOscillatory" (Ooura-Mori double exponential quadrature) and "LevinRule" (Levin's method), as well as the general controller method "OscillatorySelection". See this scicomp.SE answer for a short description of oscillatory quadrature methods and links to references, and NIntegrate Integration Strategies in the Mathematica help file for more details on how to employ them within Mathematica.

If you take a look at the documentation, Mathematica's symbolic Fourier transform function, FourierTransform, computes

$$\hat f(k) = \frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(x)e^{ikx}\mathrm{d}x$$

You can discretize some piece of this integral by limiting $x$ and $k$ to values $x_1 + (r-1)\Delta x$ and $(s-1)\Delta k$ respectively, where $\Delta x\Delta k = 2\pi/N$, giving

\begin{align}\hat f_s &= \frac{1}{\sqrt{2\pi}}\sum_{r=1}^N f(x_1 + (r-1)\Delta x)e^{i(s-1)\Delta k[x_1 + (r-1)\Delta x]}\Delta x\\ &= \frac{1}{\sqrt{2\pi}}e^{i(s-1)\Delta k x_1/N}\sum_{r=1}^N f_r\ e^{2\pi i(r-1)(s-1)/N}\Delta x\end{align}

Now compare this with the Fourier function, which calculates

$$\frac{1}{\sqrt{N}}\sum_{r=1}^{N}u_r\ e^{2\pi i(r-1)(s-1)/N}$$

Evidently

$$u_r = \sqrt{\frac{N}{2\pi}}\ f_r\exp\biggl(\frac{2\pi i(s-1)x_1}{N\Delta x}\biggr)$$

The $2\pi/\Delta x$ in that exponent is a translation generator, and $(s-1)x_1/N$ is the corresponding parameter.

So here's what you need to do to compute a numerical Fourier transform: first, choose a grid on which to sample your function, consisting of equally spaced points $x_1,\ldots,x_N$.

grid = Range[x1, x1+(n-1)*deltax, deltax]


Here x1 and x1+n*deltax should be the endpoints of the region over which you are going to compute the Fourier transform, and deltax should be some interval of $x$ small enough to capture the smallest details in your function.

Then sample your function on this grid,

samples = f/@grid


Notice that the factor $\exp\bigl(\frac{2\pi i(s-1)x_1}{N\Delta x}\bigr)$ depends only on $s$, which is a frequency space index. But it's independent of $r$, which is a position space index. So you can compute that factor separately and merge it in to the result of the FFT.

factor = Sqrt[n/(2 Pi)] Exp[2 Pi I x1/(n*deltax) Range[0, n-1]]
transform = Chop[Fourier[samples] * factor]


In order to actually turn this into the Fourier transformation of your function, you need to know that the frequencies which your transformed values correspond to start at zero, increase up to some maximum value which occurs in the middle of the array, then jump down to a negative value and increase up to zero again. (Actually, the FFT assumes momentum space is periodic and calculates from $k=0$ to $k=(N-1)\Delta k$, and you then need to map the second half of this to $k=-\frac{N}{2}\Delta k$ to $k = -\Delta k$.) You can create the frequency array using

freq = #~Join~Most@Reverse[-#]&@Range[0, Pi/deltax, 2 Pi/(n*deltax)]


and then you should be able to, say, plot your Fourier transform using

ListPlot[Transpose[freq, transform]]


(actually my expression for freq seems a little off in tests, but I'll see if I can fix it up).

• So that one doesn't have to remember to carry around multiplicative factors in both Fourier[] and FourierTransform[], one could always set the option FourierParameters beforehand, either explicitly in the functions or through SetOptions[]. – J. M. will be back soon Feb 14 '12 at 2:25
• I don't think it makes much of a difference whether the Sqrt[n] is explicit or implicitly provided through FourierParameters, since we already have the variable n for other purposes, and there is going to be a constant factor anyway whether it includes Sqrt[n] or not. I figured since Daniel already posted the compact version, I would aim for being more educational, which means making the constant factor explicit, and plus there is something to be said for the simplicity of omitting an option. – David Z Feb 14 '12 at 3:15
• Sure; as I said, if one is always using a convention different from Mathematica's, there is always SetOptions[] to get Mathematica to always use your convention instead of having to carry around factors or explicitly specify options with each call to a Fourier function. :) – J. M. will be back soon Feb 14 '12 at 3:19
• Nice exposition. Wish I'd tried to cast it that way. – Daniel Lichtblau Feb 14 '12 at 16:43
• could it be that the expression for $u_r$ is missing a $\Delta_x$? Comparing the second line for $\hat f_s$ with the Fourier expression it seems that $u_r = \Delta_x\sqrt{\frac{N}{2\pi}}\ f_r\exp\biggl(\frac{2\pi i(s-1)x_1}{N\Delta x}\biggr)$ – ftiaronsem Jul 11 '17 at 13:45

This is borrowed from comp.soft-sys.math.mathematica posts, primarily by Szesi Mukasa. The outside factor I cribbed from a text book.

fft[ll_] :=
Exp[I*Pi*(Range[Length[ll]] + Boole[OddQ[Length[ll]]])]*
RotateRight[Fourier[ll, FourierParameters -> {0, 1}],
Quotient[Length[ll], 2]]


For an example, you could use a Gaussian.

f[x_] := Exp[-x^2]
tabl = Table[f[10*x], {x, -3., 3., .1}];


There will be smallish imaginary parts so I'll just retain the real parts.

ListPlot[Chop[fft[tabl]], PlotRange -> All]


Best to test before use; I'm a bit out of my depth here.

• Thanks, but that still leaves the problem of finding the frequencies. – yohbs Feb 14 '12 at 7:29
• Omigosh, that has to be the most complicated way to write (-1)^k that I've ever conjured. I'm fortunate I can't downvote my own post. For frequencies, David Zaslavky explained it better anyway. – Daniel Lichtblau Feb 14 '12 at 18:56
• I once spent one month doing a complicated calculation of a normal mode frequency for a BEC. I obtained the result $\omega=\sqrt{2}$ in the canonical units. turned out this was correct and there is a theorem ensuring it... Taught me to search the literature. – acl Feb 16 '12 at 1:02

Coming a little late... Anyway

I'm used to the FourierParameters that MMA describes as "signal processing" in FourierTransform's More Information {0, -2 Pi}. I'll edit if I have a more general version

Set up the option SamplingPeriod, which if it's left at it's default value it chooses the Nyquist period (1/(2* maximum frequency))

Options[FourierEstimate] ^= {SamplingPeriod -> Automatic};


The function gets the expression, the variable, the maximum time and the maximum frequency component. If you tell it the max frequency is 10Hz for instance, it assumes the negative part too (-10, 10). It also assumes the signal starts at time 0. It can be easily extended to deal with signals that start at any time, by using the shift property and calling this function. It can also be easily extended to deal with functions with it's frequencies in any range, not necessarily symmetrical (complex functions). Given that it is interpolated using the default polynomial interpolation which is quite different to the ideal sinc interpolation, it may be wise to set the maximum time some factor greater than the real maximum time, or increasing the interpolation order. Probably I'll soon edit this function to make it more versatile, since it is useful for me too

FourierEstimate[expr_, var_, maxTime_, maxFreq_,
op : OptionsPattern[]] :=
With[{sp =
OptionValue[SamplingPeriod] //
If[# === Automatic, 1./(2. maxFreq), #] &},
With[{numPoints = Quotient[maxTime, sp]},
(* I protect the options so that I don't change them to the user *)
InternalInheritedBlock[{Fourier},
(* This are the parameters I am used to using with DFT *)
SetOptions[Fourier, FourierParameters -> {1, -1}];

(* Time values where I'm sampling *)
Range[0., maxTime, sp] //
(* Getting the samples *)
Table[expr, {var, #}] & //
(* Obtaining its DFT and rescaling *)
Chop@Fourier[#] sp & //
RotateRight[#, Floor[(numPoints - 1)/2]] & //
Transpose@{Range @@ ({
(* minimum frequency *) -Floor[(numPoints - 1)/2],
(* maximum frequency *) Ceiling[(numPoints + 1)/2],
1}/(sp numPoints)),
#} & //
(* Generating the interpolated function *)
Interpolation
]
]
]


Hope it helps.

Since I recently had to deal with this issue too, I wanted to provide full working example, comparing the FourierTransform result to the Fourier result. This code is based on David Z's excellent answer.

First let's define our function (in this case the ground state of the quantum harmonic oscillator):

ClearAll["Global*"]
fun[x_] := E^(-(x^2/2))/\[Pi]^(1/4)


This function is normalized to 1

Integrate[Conjugate[fun[k]]*fun[k], {k, -Infinity, Infinity}]


And its fourier transform is also normalized to 1

ftfun[k_] = FourierTransform[fun[x], x, k]
Integrate[Conjugate[ftfun[k]]*ftfun[k], {k, -Infinity, Infinity}]


Let's see whether we can reproduce this using Fourier and the details presented in David Z's answer.

First we evaluate our function on a grid and double check the normalization

n = 50;
deltax = 0.2;
x1 = (-n/2 + 1)*deltax;
grid = Range[x1, x1 + (n - 1)*deltax, deltax];
samples = fun /@ grid;
Total[Conjugate[samples]*samples]*(grid[[2]] - grid[[1]])


Now we fourier transform

factor = Sqrt[n/(2 Pi)] Exp[2 Pi I x1/(n*deltax) Range[0, n - 1]]*deltax;
transform = Chop[Fourier[samples]*factor];


and create a frequency grid

freq=RotateRight[Range[-Ceiling[n/2]+1,Floor[n/2]]/(n*deltax),Floor[n/2]+1]*2*Pi;


And the fourier transformed function is normalized

Total[Conjugate[transform]*transform]*(freq[[2]] - freq[[1]])


Let's compare some function values directly

N[ftfun[0]] - N[transform[[1]]]
N[ftfun[freq[[2]]]] - N[transform[[2]]]
N[ftfun[freq[[10]]]] - N[transform[[10]]]


and they are identical within numerical errors. Finally let's plot the fourier transform

ListPlot[Transpose[{freq, transform}], PlotRange -> Full]
`