The fast Fourier transform

The method outlined in Sect. 5.7 for solving Poisson's equation in 2-d with simple Dirichlet boundary conditions in the $y$-direction requires us to perform very many Fourier-sine transforms:

$$F_j^S = \frac{2}{J}\sum_{k=1}^{J-1} f_k\,\sin(j\,k\,\pi/J)$$ (169)

for $j=0,J$, and inverse Fourier-sine transforms:

$$f_j= \sum_{k=1}^{J-1} F_k^S\,\sin(j\,k\,\pi/J).$$ (170)

Here, $f_j$ is the value of $f(y)$ at $y_j=j\,L/J$. Thus, Eq. (169) is analogous to Eqs. (153) and (156), whereas Eq. (170) can be used to reconstruct the $u(x_i, y_j)$ from the $U_{i,j}$. Likewise, the method outlined in Sect. 5.8 for solving Poisson's equation in 2-d with simple Neumann boundary conditions in the $y$-direction requires us to perform very many Fourier-cosine transforms:

$$F_j^C = \frac{f_0}{J} + \frac{2}{J}\sum_{k=1}^{J-1} f_k\,\cos(j\,k\,\pi/J) + \frac{(-1)^j\,f_J}{J}$$ (171)

for $j=0,J$, and inverse Fourier-cosine transforms:

$$f_j= \sum_{k=0}^{J} F_k^C\,\cos(j\,k\,\pi/J).$$ (172)

Unfortunately, performing such transforms directly requires $O(J^2)$ arithmetic operations, which means that they are extremely expensive in terms of cpu resources. There is, however, an ingenious algorithm for performing Fourier transforms which only takes $O(J\,\ln J)$ arithmetic operations [which is much less than $O(J^2)$ operations when $J$ is large]. This algorithm is known as the fast Fourier transform or FFT.³⁵

The details of the FFT algorithm lie beyond the scope of this course. Roughly speaking, the algorithm works by building up the transform in stages using 2, 4, 8, 16, etc. grid-points. In this course, we shall employ the best-known publicly available FFT library, called the fftw library,³⁶ to perform all of our Fourier-sine and -cosine transforms. Unfortunately, the fftw library does not directly calculate Fourier-sine and -cosine transforms.³⁷ Instead, it calculates complex Fourier transforms:

$$F_j = \frac{1}{2J}\sum_{k=0}^{2J-1} f_k\,\exp(-{\rm i}\,j\,k\,\pi/J)$$ (173)

for $j=0,2J-1$, and complex inverse Fourier transforms:

$$f_j = \sum_{k=0}^{2J-1} F_k\,\exp(\,{\rm i}\,j\,k\,\pi/J).$$ (174)

Note that $f_j$ and $F_j$ are periodic in $j$ with period $2\,J$. Note, further, that the data-sets associated with complex Fourier transforms contain twice as many elements as the data-sets associated with sine and cosine transforms. However, we can easily convert a sine or cosine transform into a complex transform by extending its data-set. Thus, for a sine transform we write:

$$f_{2J-j} = - f_j,$$ (175)

$$F_{2J-j} = -F_j,$$ (176)

for $j=1,J-1$, in which case

$$F_j^S = 2\,{\rm i}\,F_j.$$ (177)

Likewise, for a cosine transform we write:

$$f_{2J-j} = f_j,$$ (178)

$$F_{2J-j} = F_j,$$ (179)

for $j=1,J-1$, in which case

$$F_j^C = 2\,F_j.$$ (180)

Listed below are a set of wrapper routines which employ the fftw library to perform Fourier-sine and -cosine transforms.

// FFT.cpp

// Set of functions to calculate Fourier-cosine and -sine transforms
// of real data using fftw Fast-Fourier-Transform library.
// Input/ouput arrays are assumed to be of extent J+1.
// Uses version 2 of fftw library (incompatible with vs 3).

#include <fftw.h>
#include <blitz/array.h>

using namespace blitz;

// Calculates Fourier-cosine transform of array f in array F
void fft_forward_cos (Array<double,1> f, Array<double,1>& F)
{
  // Find J. Declare local arrays.
  int J = f.extent(0) - 1;
  int N = 2 * J;
  fftw_complex ff[N], FF[N];

  // Load and extend data
  c_re (ff[0]) = f(0); c_im (ff[0]) = 0.;
  c_re (ff[J]) = f(J); c_im (ff[J]) = 0.;
  for (int j = 1; j < J; j++)
    {
      c_re (ff[j]) = f(j); c_im (ff[j]) = 0.;
      c_re (ff[2*J-j]) = f(j); c_im (ff[2*J-j]) = 0.;
    }

  // Call fftw routine
  fftw_plan p = fftw_create_plan (N, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_one (p, ff, FF);
  fftw_destroy_plan (p); 

  // Unload data
  F(0) = c_re (FF[0]); F(J) = c_re (FF[J]); 
  for (int j = 1; j < J; j++)
    {
      F(j) = 2. * c_re (FF[j]);
    }

  // Normalize data
  F /= 2. * double (J);
}

// Calculates inverse Fourier-cosine transform of array F in array f
void fft_backward_cos (Array<double,1> F, Array<double,1>& f)
{    
  // Find J. Declare local arrays.
  int J = f.extent(0) - 1;
  int N = 2 * J;
  fftw_complex ff[N], FF[N];

  // Load and extend data
  c_re (FF[0]) = F(0); c_im (FF[0]) = 0.;
  c_re (FF[J]) = F(J); c_im (FF[J]) = 0.;
  for (int j = 1; j < J; j++)
    {
      c_re (FF[j]) = F(j) / 2.; c_im (FF[j]) = 0.;
      FF[2*J-j] = FF[j];
    }

  // Call fftw routine
  fftw_plan p = fftw_create_plan (N, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_one (p, FF, ff);
  fftw_destroy_plan (p); 

  // Unload data 
  f(0) = c_re (ff[0]); f(J) = c_re (ff[J]); 
  for (int j = 1; j < J; j++)
    {
      f(j) = c_re (ff[j]);
    }
}

// Calculates Fourier-sine transform of array f in array F
void fft_forward_sin (Array<double,1> f, Array<double,1>& F)
{  
  // Find J. Declare local arrays.
  int J = f.extent(0) - 1;
  int N = 2 * J;
  fftw_complex ff[N], FF[N];
  
  // Load and extend data 
  c_re (ff[0]) = 0.; c_im (ff[0]) = 0.; 
  c_re (ff[J]) = 0.; c_im (ff[J]) = 0.;
  for (int j = 1; j < J; j++)
    {
      c_re (ff[j]) = f(j); c_im (ff[j]) = 0.;
      c_re (ff[2*J-j]) = - f(j); c_im (ff[2*J-j]) = 0.;
    }

  // Call fftw routine
  fftw_plan p = fftw_create_plan (N, FFTW_FORWARD, FFTW_ESTIMATE);
  fftw_one (p, ff, FF);
  fftw_destroy_plan (p); 

  // Unload data
  F(0) = 0.; F(J) = 0.;
  for (int j = 1; j < J; j++)
    {
      F(j) = - 2. * c_im (FF[j]);
    }  

  // Normalize data
  F /= 2. * double (J);
}

// Calculates inverse Fourier-sine transform of array F in array f
void fft_backward_sin (Array<double,1> F, Array<double,1>& f)
{ 
  // Find J. Declare local arrays.
  int J = f.extent(0) - 1;
  int N = 2 * J;
  fftw_complex ff[N], FF[N];

  // Load and extend data
  c_re (FF[0]) = 0.; c_im (FF[0]) = 0.; 
  c_re (FF[J]) = 0.; c_im (FF[J]) = 0.;
  for (int j = 1; j < J; j++)
    {
      c_re (FF[j]) = 0.; c_im (FF[j]) = - F(j) / 2.;
      c_re (FF[2*J-j]) = 0.; c_im (FF[2*J-j]) = F(j) / 2.;
    }

  // Call fftw routine
  fftw_plan p = fftw_create_plan (N, FFTW_BACKWARD, FFTW_ESTIMATE);
  fftw_one (p, FF, ff);
  fftw_destroy_plan (p); 

  // Unload data
  f(0) = 0.; f(J) = 0.;
  for (int j = 1; j < J; j++)
    {
      f(j) = c_re (ff[j]);
    }
}