105 lines
3.4 KiB
C++
105 lines
3.4 KiB
C++
// 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;
|
|
}
|