Numerical Fourier transform of a complicated function

Question

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.

Answer 1

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.

Answer 2

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).

Answer 3

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.

Answer 4

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 *)  
   Internal`InheritedBlock[{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.

Answer 5

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]