The package VQM`FastFourier`

Bernd Thaller
Institute of Mathematics
University of Graz
Austria
04-10-2004

Abstract

Mathematica has a very efficient algorithm for calculating the Fourier transformation of a list. We use it for an approximate calculation of the Fourier transformation of a function.

This package VQM`FastFourier` provides commands for the numerical Fourier transform of complex-valued functions that are given as a list of complex values on a suitable grid of space points. This package is part of the VQM packages which can be obtained here:  http://www.uni-graz.at/imawww/vqm/software.html (free download). The VQM packages are part of the Visual Quantum Mechanics project, see http://www.uni-graz.at/imawww/vqm/. In particular, these packages are distributed with the book Advanced Visual Quantum Mechanics, Springer-Verlag New York, 2004.

This notebook describes the application of Mathematica's fast Fourier algorithm to the Fourier transform of a complex-valued function. The package VQM`FastFourier` implements this method. The

Fourier transform of a list

The discrete Fourier transform of a list of n complex numbers {(Subscript[a, 1]), (Subscript[a, 2]), ... , (Subscript[a, n])} is again a list of n complex numbers which is defined by

(1)        (Subscript[(Overscript[a, ⋀]), s]) = (1/(Sqrt[n]))(∑_(r = 1)^n)(Subscript[a, r])Exp[-i((2π)/(n ))(s-1)(r-1)],    s = 1, ... , n .

In Mathematica, this operation is implemented by the command InverseFourier:

{(Subscript[(Overscript[a, ⋀]), 1]),... , (Subscript[(Overscript[a, ⋀]), n])} = InverseFourier[{(Subscript[a, 1]),... , (Subscript[a, n])}]

As the name suggests this function is inverse to Fourier[], which is defined by a reversed sign in the exponent. (The name is only a matter of convention. For the purpose of quantum mechanics, the definition with a minus in the exponent is more common).

The formula (1) has a certain symmetry. If we multiply the whole equation (1) by the s-dependent phase factor Exp[2π i((m(s-1))/n)], this is the same as a cyclic permutation of the coefficients (Subscript[a, r]):

(Subscript[(Overscript[a, ⋀]), s])Exp[2π i((m(s-1))/n)]= (1/(Sqrt[n]))(∑_(r = 1)^n)(Subscript[a, r])Exp[-2π i(((r - 1 - m) (s - 1))/n)]= (1/(Sqrt[n]))(∑_(r = 1)^n)(Subscript[b, r])Exp[-i((2π)/(n ))(s-1)(r-1)]

where (Subscript[b, r ])= (Subscript[a, (r + m) (mod n)]), or

{(Subscript[b, 1]),  ...  , (Subscript[b, n])} = RotateLeft[{(Subscript[a, 1]),  ...  , (Subscript[a, n])}, m].

Fourier transform of a function

For an integrable function ψ(x), the Fourier transform is given by

(2)        (Overscript[ψ, ⋀])(k) = (1/(Sqrt[2  π]))(∫_(-∞)^∞)ψ(x)  (e^(-i k x))dx.

We may approximate this integral by a finite Riemann sum. To this purpose we assume that the function ψ is reasonably well localized within a finite interval [(Subscript[x, 0]), (Subscript[x, n])]. We divide this interval into n subintervals by introducing equally spaced points (Subscript[x, r]),

(Subscript[x, 0]) < (Subscript[x, 1 ])< (Subscript[x, 2 ])< ... < (Subscript[x, n])

and define

dx = (Subscript[x, r]) - (Subscript[x, r - 1]) .

The Riemann sum approximating the Fourier integral is then given by

(3)        (Overscript[ψ, ⋀])(k) ≈ (1/(Sqrt[2 π])) (∑_(r = 1)^n)ψ((Subscript[x, r]))(e^(-i k (Subscript[x, r])))dx

Applying the command InverseFourier

This sum (3) looks already similar to (1). We set

(Subscript[k, s]) = ((2  π)/(n dx))(s - 1)  and  (Subscript[x, r]) = (Subscript[x, 1]) + (r - 1)dx

and obtain

(4)        (Overscript[ψ, ⋀])((Subscript[k, s])) ≈ (1/(Sqrt[2  π])) (∑_(r = 1)^n)ψ((Subscript[x, r]))(e^(-i (Subscript[k, s])(Subscript[x, 1])))Exp[-i((2π)/(n ))(s-1)(r-1)] dx

We see that the built-in discrete Fourier transform is obviously suitable to calculate (a numerical approximation of) the Fourier integral (Overscript[ψ, ⋀])(k)at the points (Subscript[k, s]) which are given by

{(Subscript[k, 1]), (Subscript[k, 2]), (Subscript[k, 3]), ... , (Subscript[k, n])} = {0, ((2  π)/(n dx)), ((4  π)/(n dx)), ... , ((2 (n - 1) π)/(n dx))} ,

These points form a grid in k-space where the distance between neighbouring points is

dk = ((2  π)/(n dx))= ((2  π)/((Subscript[x, n])-(Subscript[x, 0])))

and the length of the interval in Fourier space is  n dk = 2π/dx. It depends on the distance dx between the grid points in x-space.

Eq. (4) can be described as the (inverse) Fourier transform of a list, except that the result is multiplied by an s-dependent phase factor (e^(i (Subscript[k, s])(Subscript[x, 1]))).

(1/(Sqrt[2 π])) (∑_(r = 1)^n)ψ((Subscript[x, r]))(e^(-i (Subscript[k, s])(Subscript[x, 1])))Exp[-i((2π)/(n ))(s-1)(r-1)]dx
            = (Sqrt[(n/(2 π))]) dx (1/(Sqrt[n])) (∑_(r = 1)^n)ψ((Subscript[x, r]))(e^(-i (Subscript[k, s])(Subscript[x, 1])))Exp[-i((2π)/(n ))(s-1)(r-1)]dx

The element with index s in the list InverseFourier[{ψ((Subscript[x, 1])), ... ,ψ((Subscript[x, n]))}] has to be multiplied by the phase factor  (e^(i (Subscript[k, s])(Subscript[x, 1]))). This multiplication by an s-dependent phase factor can be achieved more effectively by a cyclic permutation of the elements {ψ((Subscript[x, 1])), ... ,ψ((Subscript[x, n]))}. We define

ind[x] := ((x-(Subscript[x, 0]))/(dx))

This is the index number of the point x in the sequence {(Subscript[x, 1]),(Subscript[x, 2]), ... ,(Subscript[x, n])}. Hence ind[(Subscript[x, 1])]=1, and ind[0] = -(Subscript[x, 0])/dx. In Mathematica, the cyclic permutation is performed by the command RotateLeft:

{(Subscript[a, 1]),   ...  , (Subscript[a, n])} = RotateLeft[{ψ((Subscript[x, 1])), ... ,ψ((Subscript[x, n]))},ind[0]]

This has precisely the same effect as the multiplication with the phase factor in (4). Hence the Fourier transform of the function is given by

(4')        (Overscript[ψ, ⋀])((Subscript[k, s])) ≈ (Sqrt[(n/(2 π))]) dx (1/(Sqrt[n]))  (∑_(r = 1)^n)(Subscript[a, r])Exp[-i((2π)/(n ))(s-1)(r-1)]
                             =  (Sqrt[(n/(2 π))]) dx  InverseFourier[{(Subscript[a, 1]),   ...  , (Subscript[a, n])}].

Discrete Fourier sums can be treated as periodic functions. In particular, if p = 2π/dx is the length of the interval in k-space, we can use Eq. (3) and find, using (e^(i (α ± 2 π))) = (e^(i α))

(Overscript[ψ, ⋀])((Subscript[k, s])±p)= ((dx)/((Sqrt[2 π]))) (∑_(r = 1)^n)(Subscript[a, r])(e^(-i ((Subscript[k, s]) ± p)dx (r-1))) = ((dx)/((Sqrt[2 π]))) (∑_(r = 1)^n)(Subscript[a, r])(e^(-i (Subscript[k, s]) dx (r-1) ))= (Overscript[ψ, ⋀])((Subscript[k, s])).

Usually, we are interested in an interval in momentum space, which is symmetric with respect to the origin. Using the periodicity in k-space, we can correct this again by a RotateLeft command. Since  (Overscript[ψ, ⋀])((Subscript[k, s])) is periodic, we obtain the values of (Overscript[ψ, ⋀]) at the n points:

{-(π/(dx))+ dk, ... , -dk, 0, dk, 2dk, ... , (π/( dx))}

(which are obtained from {(Subscript[k, 1]), (Subscript[k, 2]), (Subscript[k, 3]), ... , (Subscript[k, n])} by shifting these points n/2 + 1 places towards the left) by performing the cyclic permutation

RotateLeft[{ (Overscript[ψ, ⋀])((Subscript[k, 1])), (Overscript[ψ, ⋀])((Subscript[k, 2])), ... ,  (Overscript[ψ, ⋀])((Subscript[k, n]))}, n/2+1] .

Practical implementation

Assume that we are given a function f[x] which vanishes outside an interval [a, b]. We want to determine a discrete approximation for the Fourier transform, suitable for plotting it quickly.
We first choose a grid of points on the x-axis. We assume that the grid points are equally spaced by dx and given by the numbers n dx. Since it is best to choose n as a power of 2, we calculate dx:

dx = ((b - a)/n)

xgrid = Table[Floor[(a/dx)] dx + k dx, {k, 1, n}]

Of course, the size of dx depends on the variation of the given function. dx determines the size of the interval in k-space, where the values of the Fourier transform are calculated. The spacing of the points in k-space depends on the size of the interval [a,b]. If we want to evaluate the k-values at points which are closer together, we have to choose the interval in x-space larger. The points in xgrid are inside the significant interval [a,b]. We determine the function values at these points:

fvals = f /@ xgrid

Next we rotate the list of function values

ind[x_] := Floor[((x - a)/dx)] + 1

RotateLeft[fvals, ind[0] - 1]

From this, we calculate the discrete Fourier transform

(Sqrt[(n/(2 π))]) dx InverseFourier[
        RotateLeft[fvals, ind[0] - 1]]

The result has been multiplied by the appropriate factor. Finally, the result is rotated by n/2

RotateLeft[(Sqrt[(n/(2 π))]) dx InverseFourier[
          RotateLeft[fvals, ind[0] - 1]], Round[(n/2) + 1]]

Explicit example

This is the function we want to transform:

f[
          x_] = ((((1/π)))^(1/4)) Exp[-((x^2)/2) + 
              I 2 x];

We choose an interval and grid points in x-space:

n = 512;
a = -10; b = 14;

dx = ((b - a)/n);

xgrid = Table[Floor[(a/dx)] dx + k dx, {k, 1, n}];

Note that 0 is an element of xgrid. Now we evaluate the function on the grid points. This gives the discretization of the function f:

fvals = f /@ xgrid;

The following definition is useful for identifying the value of f at a particular point:

ind[x_] := Floor[((x - a)/dx)] + 1

Next, we compute the Fourier transform of this list, as described above:

kvals = 
        Chop[RotateLeft[(Sqrt[(n/(2 π))]) dx InverseFourier[
                RotateLeft[fvals, ind[0] - 1]], Round[(n/2)]]];

This list may be interpreted as a discretization of the Fourier transform of the function f. The grid in k-space is given by:

kleft = -(π/dx); kright = (π/dx); 
    dk = ((2 π)/(n dx));

In order to visualize the result, we load the package

<< "VQM`ArgColorPlot`"

This is the original function (or rather its discretization)

QListArgColorPlot[Take[fvals, {ind[-3], ind[3]}], 
        Frame -> True, 
        FrameTicks -> {QNiceTicks[-3, 3, dx], 
            Automatic}];

The following function is needed in order to extract parts of the Fourier-transformed list:

indk[k_] := Floor[((k - kleft)/dk)] + 1;

Now we can plot the interesting part of the list kvals

QListArgColorPlot[Take[kvals, {indk[-1], indk[5]}], 
        Frame -> True, 
        FrameTicks -> {QNiceTicks[-1, 5, dk], Automatic}, 
        PlotRange -> All];

We can easily combine the computation above in a single function. The following computes the Fourier transform of a numerically given function xlist and returns the list of Fourier-transformed values together with the information about the grid in k-space:

fourierList[xlist_, {xleft_, xright_}] := 
        Module[{n, dx, a, b, klist, kleft, kright, dk}, 
          n = Length[xlist] - 1; dx = ((xright - xleft)/n); 
          a = -(xleft/dx); b = (xright/dx); 
          klist = (Sqrt[(n/(2 π))]) dx RotateLeft[
                InverseFourier[RotateLeft[xlist, Round[a]]], 
                Round[(n/2) + 1]]; kleft = -(π/dx); 
          kright = (π/dx); {klist, {kleft, kright}}];

This function is implemented as QFourierList in the package VQM`FastFourier`

The package VQM`FastFourier`

The results described above form the theoretical background of the package VQM`FastFourier`. Here we show how to obtain the Fourier transform of a function f with the help of this package. The package VQM`ArgColorPlot` will be useful for the visualization of the function and its Fourier transform.

This loads all declarations of the VQM packages:

<<VQM`

The most important function defined by VQM`FastFourier` is QFourierList. (Note that all symbols defined by the VQM packages start with the letter Q.)

The first argument of QFourierList is a list flist of complex numbers. You may interpret this list as the discretization of a function defined in some interval [a,b]. The second argument of QFourierList specifies this interval. Hence the function is called like this: QFourierList[flist,{a,b}]. In order to have an example, we generate a list of complex values from a function f. The package VQM`FastFourier` provides some auxiliary functions for this.

We repeat the example above with the help of this package and consider the function

f[
          x_] = ((((1/π)))^(1/4)) Exp[-((x^2)/2) + 
              I 2 x];

in the interval [-10,14]:

QArgColorPlot[f[x],{x,-10,14},PlotRange->All]

Let us first generate a suitable grid in x-space:

a = -10.; b = 14.; n = 512;

xgrid = QGrid[a,b,n];

This command divides the interval [-10,10] by equally spaced points (Subscript[x, 0]) =a< (Subscript[x, 1 ])< (Subscript[x, 2 ])< ... < (Subscript[x, n])= b and returns the list { (Subscript[x, 1]), (Subscript[x, 2]), ... , (Subscript[x, n])}. The number n should be even. You can obtain information about the space grid via the following functions.

{QLeftBorder[a,b,n],QRightBorder[a,b,n],QSpaceStep[a,b,n]}

QIndexPosition[x,a,b,n] lets you access the element of xgrid corresponding to a special value of x

QIndexPosition[5.,a,b,n]

xgrid[[321]]

We can now define a discretization of the function f by

fvals = f /@ QGrid[a,b,n];

Usually, the function in position space will be the result of a numerical computation, and the list fvals of complex numbers is the object to start with.

QListArgColorPlot[fvals,PlotRange->All]

The option QHorizontalRange lets you specify the plot range on the horizontal axis:

QListArgColorPlot[fvals,PlotRange->All, QHorizontalRange->QSpaceInterval[a,b,n]]

QListArgColorPlot[fvals,PlotRange->All, QHorizontalRange->{QSpaceInterval[a,b,n],{-5,5}}]

The points in x-space determine the interval and the grid in k-space, as explained above. You can obtain a list of the points in k-space by QFourierGrid[-10.,10,100]. Information about the grid in Fourier space can be obtained by

{kleft,kright,dk}={QFourierLeftBorder[a,b,n],QFourierRightBorder[a,b,n],QFourierStep[a,b,n]}

Hence

QFourierGrid[a,b,n]==Table[n,{n,kleft,kright,dk}]

If the function in x-space is given on the points of QGrid, then QFourierList returns the values of the Fourier transform on the points QFourierGrid in k-space.

Ff = QFourierList[fvals,{a,b}];

QFourierList returns the list of Fourier values and the interval in k-space: {Fflist,{kleft,kright}}.

Ffvals = Ff[[1]]; kinterval = Ff[[2]];

QListArgColorPlot[Ffvals,PlotRange->All, QHorizontalRange->{kinterval,{-1,5}},Frame->True]

This is the discretization of a shifted Gaussian function.

We can (approximately) reverse the whole process by QInverseFourierList.

result = QInverseFourierList[Ff];

QListArgColorPlot[result[[1]],PlotRange->All,QHorizontalRange->result[[2]]]

Now the interval in x-space is symmetric with respect to x = 0.

Computing the Fourier transformation again gives the same result as before (up to numerical errors)

{testFfvals, testkinterval} = QFourierList[result];

QListArgColorPlot[testFfvals,PlotRange->All,QHorizontalRange->{testkinterval,{-1,5}},Frame->True]

The numerical error is:

Max[Abs[testFfvals-Ffvals]]