//   fft.cpp - impelementation of class
//   of fast Fourier transform - FFT
//

#include "fft.h"
#include <math.h>

//   FORWARD FOURIER TRANSFORM, INPLACE VERSION
//     Data - both input data and output
//     N    - length of input data
bool CFFT::Forward (complex * const Data) const {
  //   Check input parameters
  if (!Data || N < 1 || N & (N - 1)) return false;
  //   Rearrange
  Rearrange (Data);
  //   Call FFT implementation
  Perform (Data);
  //   Succeeded
  return true;
}
//   INVERSE FOURIER TRANSFORM, INPLACE VERSION
//     Data  - both input data and output
//     N     - length of both input data and result
//     Scale - if to scale result
bool CFFT::Inverse (complex * const Data, const bool Sc /* = true */) const {
  //   Check input parameters
  if (!Data || N < 1 || N & (N - 1)) return false;
  //   Rearrange
  Rearrange (Data);
  //   Call FFT implementation
  Perform (Data, true);
  //   Scale if necessary
  if (Sc)  Scale (Data);
  //   Succeeded
  return true;
}
//   Inplace version of rearrange function
void CFFT::Rearrange (complex * const Data) const {
  //   Swap position
  unsigned int Target = 0;
  //   Process all positions of input signal
  for (unsigned int Position = 0; Position < N; ++Position) {
    //   Only for not yet swapped entries
    if (Target > Position) {
      //   Swap entries
      const complex Temp (Data[Target]);
      Data[Target]   = Data[Position];
      Data[Position] = Temp;
    }
    //   Výpočet Target (v podstatě reverzace pořadí bitů Position)
    unsigned int Mask = N;
    //   While bit is set
    while (Target & (Mask >>= 1))
      //   Drop bit
      Target &= ~Mask;
    //   The current bit is 0 - set it
    Target |= Mask;
  }
}

//   FFT implementation
void CFFT::Perform (complex * const Data, const bool Inverse /* = false */) const {
  complex Product;
  unsigned int Step;
  const real pi = Inverse ? 3.14159265358979323846 : -3.14159265358979323846;
  //   Iteration through dyads, quadruples, octads and so on...
  for (Step = 1; Step < N; Step <<= 1) {  // 1,2,...N/2
    //   Jump to the next entry of the same transform factor Jump = Step * 2
    const unsigned int Jump = Step << 1;            // 2,4,...N
    //   Angle increment
    const real delta = pi  / real (Step);
    const real incr  = 1.0 / real (Step);
    //   Multiplier for trigonometric recurrence
    const complex Multiplier (cos (delta), sin (delta));
    //   Start value for transform factor, fi = 0
    complex Factor (1.0);
    real    rot = 0.0;
    //   Iteration through groups of different transform factor
    for (unsigned int Group = 0; Group < Step; ++Group) {
      //   Iteration within group
      for (unsigned int Pair = Group; Pair < N; Pair += Jump) {
        //   Match position
        const unsigned int Match = Pair + Step;
        //   Second term of two-point transform
        Product = Factor * Data[Match];
        //   Transform for fi + pi
        Data[Match] = Data[Pair] - Product;
        //   Transform for fi
        Data[Pair]  = Data[Pair] + Product;
      }
      //   Successive transform factor via trigonometric recurrence
      //   Factor = Multiplier * Factor + Factor;
      Factor *= Multiplier;
      rot    += incr;
    }
  }
}

//   Scaling of inverse FFT result
void CFFT::Scale (complex *const Data) const {
  const real Factor = 1.0 / real (N);
  //   Scale all data entries
  for (unsigned int Position = 0; Position < N; ++Position)
    Data[Position] *= Factor;
}