r2fft_test.cc

#include <stdio.h>
#include <stdlib.h>

#include "r2fft.h"
#include "fft.h"

#define Float double

template <class Type>
double Power(std::complex<Type> X) {
    return real(X) * real(X) + imag(X) * imag(X);
}

int main(int argc, char *argv[]) {
    r2FFT<Float> FFT;
    DFT1d<Float> RefFFT;

    int FFTsize = 1024;
    std::complex<Float> Buffer[FFTsize];

    FFT.Preset(FFTsize);
    RefFFT.PresetForward(FFTsize);

    srand(123456);
    for (int Idx = 0; Idx < FFTsize; Idx++) {
        Buffer[Idx] = std::complex<Float>((rand()&3) - 1.5, (rand()&3) - 1.5);
        RefFFT.Buffer[Idx] = Buffer[Idx];
    }

    FFT.Process(Buffer);
    RefFFT.Execute();

    double AverPower = 0;
    for (int Idx = 0; Idx < FFTsize; Idx++) {
        AverPower += Power(Buffer[Idx]);
    }
    AverPower /= FFTsize;
    printf("FFT: <Power> = %8.3f\n", AverPower);

    double RefAverPower = 0;
    for (int Idx = 0; Idx < FFTsize; Idx++) {
        RefAverPower += Power(RefFFT.Buffer[Idx]);
    }
    RefAverPower /= FFTsize;
    printf("RefFFT: <Power> = %8.3f\n", AverPower);

    double DiffAverPower = 0;
    for (int Idx = 0; Idx < FFTsize; Idx++) {
        DiffAverPower += Power(Buffer[Idx] - RefFFT.Buffer[Idx]);
    }
    DiffAverPower /= FFTsize;
    printf("FFT-RefFFT: <Power> = %8.3f\n", DiffAverPower);

    return 0;
}