fft.cpp

//   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 input data
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;
}