Fast Fourier transform — FFT — source code

// fft.cpp - implementation of class
// of fast Fourier transform - FFT
//
// The code is property of LIBROW
// You can use it on your own
// When utilizing credit LIBROW site
// Include declaration file
#include "fft.h"
// Include math library
#include <math.h>
// FORWARD FOURIER TRANSFORM
// Input - input data
// Output - transform result
// N - length of both input data and result
bool CFFT::Forward(const complex *const Input, complex *const Output,
 const unsigned int N)
{
 // Check input parameters
 if (!Input || !Output || N < 1 || N & (N - 1)) return false;
 // Initialize data
 Rearrange(Input, Output, N);
 // Call FFT implementation
 Perform(Output, N);
 // Succeeded
 return true;
}
// FORWARD FOURIER TRANSFORM, INPLACE VERSION
// Data - both input data and output
// N - length of both input data and result
bool CFFT::Forward(complex *const Data, const unsigned int N)
{
 // Check input parameters
 if (!Data || N < 1 || N & (N - 1)) return false;
 // Rearrange
 Rearrange(Data, N);
 // Call FFT implementation
 Perform(Data, N);
 // Succeeded
 return true;
}
// INVERSE FOURIER TRANSFORM
// Input - input data
// Output - transform result
// N - length of both input data and result
// Scale - if to scale result
bool CFFT::Inverse(const complex *const Input, complex *const Output,
 const unsigned int N, const bool Scale /* = true */)
{
 // Check input parameters
 if (!Input || !Output || N < 1 || N & (N - 1)) return false;
 // Initialize data
 Rearrange(Input, Output, N);
 // Call FFT implementation
 Perform(Output, N, true);
 // Scale if necessary
 if (Scale)
 CFFT::Scale(Output, N);
 // 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 unsigned int N,
 const bool Scale /* = true */)
{
 // Check input parameters
 if (!Data || N < 1 || N & (N - 1)) return false;
 // Rearrange
 Rearrange(Data, N);
 // Call FFT implementation
 Perform(Data, N, true);
 // Scale if necessary
 if (Scale)
 CFFT::Scale(Data, N);
 // Succeeded
 return true;
}
// Rearrange function
void CFFT::Rearrange(const complex *const Input, complex *const Output,
 const unsigned int N)
{
 // Data entry position
 unsigned int Target = 0;
 // Process all positions of input signal
 for (unsigned int Position = 0; Position < N; ++Position) { // Set data entry
 Output[Target] = Input[Position];
 // Bit mask
 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;
 }
}
// Inplace version of rearrange function
void CFFT::Rearrange(complex *const Data, const unsigned int N)
{
 // 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;
 }
 // Bit mask
 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 unsigned int N,
 const bool Inverse /* = false */)
{
 const double pi = Inverse ? 3.14159265358979323846 : -3.14159265358979323846;
 // Iteration through dyads, quadruples, octads and so on...
 for (unsigned int Step = 1; Step < N; Step <<= 1) { // Jump to the next entry of the same transform factor
 const unsigned int Jump = Step << 1; // Angle increment
 const double delta = pi / double(Step);
 // Auxiliary sin(delta / 2)
 const double Sine = sin(delta * .5);
 // Multiplier for trigonometric recurrence
 const complex Multiplier(-2. * Sine * Sine, sin(delta));
 // Start value for transform factor, fi = 0
 complex Factor(1.);
 // 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
 const complex Product(Factor * Data[Match]);
 // Transform for fi + pi
 Data[Match] = Data[Pair] - Product;
 // Transform for fi
 Data[Pair] += Product;
 }
 // Successive transform factor via trigonometric recurrence
 Factor = Multiplier * Factor + Factor;
 }
 }
}
// Scaling of inverse FFT result
void CFFT::Scale(complex *const Data, const unsigned int N)
{
 const double Factor = 1. / double(N);
 // Scale all data entries
 for (unsigned int Position = 0; Position < N; ++Position) Data[Position] *= Factor; }