fourier.C

#include "random.h"
#include <cmath>
#include <utility>

/*
    GenRan: non-Gaussian random field generator.
    Copyright (C) 2006 Peter Grassl
    Distributed under the MIT licence.

    Radix-2 Cooley-Tukey FFT, both directions, in place. The signature uses
    a 1-based real array of length 2*nn that encodes nn complex numbers as
    interleaved (real, imag) doubles: data[1..2*nn]. Callers therefore pass
    `buffer - 1` for a normal 0-based buffer of length 2*nn.

    Sign convention:
        isign = +1 :  X_k = sum_{j=0}^{nn-1} x_j * exp(+i 2*pi * j * k / nn)
        isign = -1 :  X_k = sum_{j=0}^{nn-1} x_j * exp(-i 2*pi * j * k / nn)

    Neither direction normalises by nn. nn must be a positive power of two
    (not checked here).
*/

void four1(double data[], int nn, int isign)
{
  // Switch to a 0-based view of the buffer.
  double *d = data + 1;
  const double two_pi = 6.283185307179586476925286766559;

  // Bit-reversal permutation of the nn complex elements (XOR algorithm).
  for (int i = 1, j = 0; i < nn; ++i) {
    int bit = nn >> 1;
    while (j & bit) {
      j ^= bit;
      bit >>= 1;
    }
    j ^= bit;
    if (i < j) {
      std::swap(d[2 * i],     d[2 * j]);
      std::swap(d[2 * i + 1], d[2 * j + 1]);
    }
  }

  // Iterative Cooley-Tukey butterflies, doubling stage size each pass.
  for (int len = 2; len <= nn; len <<= 1) {
    const int half = len >> 1;
    const double theta = isign * two_pi / static_cast<double>(len);
    const double w_step_re = std::cos(theta);
    const double w_step_im = std::sin(theta);

    for (int block = 0; block < nn; block += len) {
      double w_re = 1.0;
      double w_im = 0.0;
      for (int k = 0; k < half; ++k) {
        const int p = 2 * (block + k);
        const int q = 2 * (block + k + half);
        const double t_re = w_re * d[q]     - w_im * d[q + 1];
        const double t_im = w_re * d[q + 1] + w_im * d[q];
        d[q]     = d[p]     - t_re;
        d[q + 1] = d[p + 1] - t_im;
        d[p]     += t_re;
        d[p + 1] += t_im;
        const double w_re_new = w_re * w_step_re - w_im * w_step_im;
        w_im = w_re * w_step_im + w_im * w_step_re;
        w_re = w_re_new;
      }
    }
  }
}