FFT.cpp

Go to the documentation of this file.

/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */

#include "FFT.h"

#include <fftw.h>

FFT::FFT(unsigned int size) :
    m_size(size),
    m_table(0),
    m_spare(0)
{
    if (m_size < 2) throw InvalidSize;
    if (m_size & (m_size-1)) throw InvalidSize;

    m_table = new int[m_size];
    
    unsigned int bits;
    unsigned int i, j, k, m;

    for (i = 0; ; ++i) {
        if (m_size & (1 << i)) {
            bits = i;
            break;
        }
    }

    for (i = 0; i < m_size; ++i) {
        
        m = i;
        
        for (j = k = 0; j < bits; ++j) {
            k = (k << 1) | (m & 1);
            m >>= 1;
        }
        
        m_table[i] = k;
    }

    m_spare = new double[m_size];
}

FFT::~FFT()
{
    delete[] m_table;
    delete[] m_spare;
}

void
FFT::tune()
{
    double *a = new double[m_size];
    double *b = new double[m_size];
    double *c = new double[m_size];
    double *d = new double[m_size];
    fftw_complex *ac = new fftw_complex[m_size];
    
    for (unsigned int i = 0; i < m_size; ++i) {
        a[i] = drand48();
        b[i] = 0;
    }

    fftw_plan fplan = fftw_plan_dft_r2c_1d(m_size, a, ac, FFTW_MEASURE);
    fftw_plan iplan = fftw_plan_dft_c2r_1d(m_size, ac, a, FFTW_MEASURE);

    Profiler *profiler = new Profiler("FFT::tune: FFTW", true);

    for (unsigned int count = 0; count < 1000; ++count) {
        fftw_execute(fplan);
        fftw_execute(iplan);
    }

    delete profiler;
    
    for (unsigned int i = 0; i < m_size; ++i) {
        a[i] = drand48();
        b[i] = 0;
    }

    profiler = new Profiler("FFT::tune: basefft", true);

    for (unsigned int count = 0; count < 1000; ++count) {
        basefft(false, a, b, c, d);
        basefft(true, a, b, c, d);
    }

    delete profiler;
}    

void
FFT::forward(double *realIn, double *realOut, double *imagOut)
{
    basefft(false, realIn, 0, realOut, imagOut);
}

void
FFT::inverse(double *realIn, double *imagIn, double *realOut)
{
    basefft(true, realIn, imagIn, realOut, m_spare);
}

void
FFT::basefft(bool inverse, double *ri, double *ii, double *ro, double *io)
{
    if (!ri || !ro || !io) return;

    unsigned int i, j, k, m;
    unsigned int blockSize, blockEnd;

    double tr, ti;

    double angle = 2.0 * M_PI;
    if (inverse) angle = -angle;

    const unsigned int n = m_size;

    if (ii) {
        ro[table[0]] = ri[0];
        io[table[0]] = ii[0];
        for (i = 1; i <= n/2; ++i) {
            ro[table[i]] = ri[i];
            io[table[i]] = ii[i];
            ro[table[n-i]] = ri[i];
            io[table[n-i]] = -ii[i];
        }
    } else {
        for (i = 0; i < n; ++i) {
            ro[table[i]] = ri[i];
            io[table[i]] = 0.0;
        }
    }

    blockEnd = 1;

    for (blockSize = 2; blockSize <= n; blockSize <<= 1) {

        double delta = angle / (double)blockSize;
        double sm2 = -sin(-2 * delta);
        double sm1 = -sin(-delta);
        double cm2 = cos(-2 * delta);
        double cm1 = cos(-delta);
        double w = 2 * cm1;
        double ar[3], ai[3];

        for (i = 0; i < n; i += blockSize) {

            ar[2] = cm2;
            ar[1] = cm1;

            ai[2] = sm2;
            ai[1] = sm1;

            for (j = i, m = 0; m < blockEnd; j++, m++) {

                ar[0] = w * ar[1] - ar[2];
                ar[2] = ar[1];
                ar[1] = ar[0];

                ai[0] = w * ai[1] - ai[2];
                ai[2] = ai[1];
                ai[1] = ai[0];

                k = j + blockEnd;
                tr = ar[0] * ro[k] - ai[0] * io[k];
                ti = ar[0] * io[k] + ai[0] * ro[k];

                ro[k] = ro[j] - tr;
                io[k] = io[j] - ti;

                ro[j] += tr;
                io[j] += ti;
            }
        }

        blockEnd = blockSize;
    }

/* fftw doesn't normalise, so nor will we

    if (inverse) {

        double denom = (double)n;

        for (i = 0; i < n; i++) {
            ro[i] /= denom;
            io[i] /= denom;
        }
    }
*/
}

---