fft.c

Go to the documentation of this file.

#include <grass/config.h>

#if defined(HAVE_FFTW_H) || defined(HAVE_DFFTW_H) || defined(HAVE_FFTW3_H)

#ifdef HAVE_FFTW_H
#include <fftw.h>
#endif

#ifdef HAVE_DFFTW_H
#include <dfftw.h>
#endif

#ifdef HAVE_FFTW3_H
#include <fftw3.h>
#define c_re(c) ((c)[0])
#define c_im(c) ((c)[1])
#endif

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <grass/gmath.h>
#include <grass/gis.h>


int fft2(int i_sign, double (*data)[2], int NN, int dimc, int dimr)
{
#ifdef HAVE_FFTW3_H
    fftw_plan plan;
#else
    fftwnd_plan plan;
#endif
    double norm;
    int i;

    norm = 1.0 / sqrt(NN);

#ifdef HAVE_FFTW3_H
    plan = fftw_plan_dft_2d(dimr, dimc, data, data,
                            (i_sign < 0) ? FFTW_FORWARD : FFTW_BACKWARD,
                            FFTW_ESTIMATE);

    fftw_execute(plan);

    fftw_destroy_plan(plan);
#else
    plan = fftw2d_create_plan(dimc, dimr,
                              (i_sign < 0) ? FFTW_FORWARD : FFTW_BACKWARD,
                              FFTW_ESTIMATE | FFTW_IN_PLACE);

    fftwnd_one(plan, data, data);

    fftwnd_destroy_plan(plan);
#endif

    for (i = 0; i < NN; i++) {
        data[i][0] *= norm;
        data[i][1] *= norm;
    }

    return 0;
}

int fft(int i_sign, double *DATA[2], int NN, int dimc, int dimr)
{
    fftw_complex *data;
    int i;

    data = (fftw_complex *) G_malloc(NN * sizeof(fftw_complex));

    for (i = 0; i < NN; i++) {
        c_re(data[i]) = DATA[0][i];
        c_im(data[i]) = DATA[1][i];
    }

    fft2(i_sign, data, NN, dimc, dimr);

    for (i = 0; i < NN; i++) {
        DATA[0][i] = c_re(data[i]);
        DATA[1][i] = c_im(data[i]);
    }

    G_free(data);

    return 0;
}

#endif /* HAVE_FFT */