I have a real matrix represented by the structure double input[N][M]
where I want to take a 1D FT for each column j = 0..M along the N direction. My attempt at this is as follows:
#include <fftw3.h>
#include <math.h>
#define M_PI 3.14159265358979323846264338327
int main(void)
{
int N = 16; int M = 2;
double input[N][M]; // input data
double *in = (double *) input;
fftw_complex output[N][M]; // output data
fftw_complex *out = (fftw_complex *) output;
for (int i = 0; i < N; ++i)
{
for (int j = 0; j < M; ++j)
{
double x = (double) i / (double) N;
input[i][j] = sin(2*(j+1) * M_PI * x);
fprintf (stderr, "input[%d][%d] = %.9f, ", i, j, input[i][j]);
}
fprintf (stderr, "\n");
}
fprintf (stderr, "\n");
// setup plans
int rank = 1; int howmany = M;
int istride = M; int ostride = M; int idist = 1; int odist = 1;
fftw_plan fp = fftw_plan_many_dft_r2c(rank, &N, howmany,
in, NULL, istride, idist,
out, NULL, ostride, odist,
FFTW_MEASURE);
// take the forward transform
fftw_execute(fp); fprintf (stderr, "FFT complete\n");
for (int j = 0; j < M; ++j)
for (int i = 0; i < N/2; ++i)
fprintf (stderr, "OUT[%3d] = (%.4f + %.4fi)\n",
i, output[i][j][0], output[i][j][1]);
fprintf (stderr, "\n");
fftw_destroy_plan(fp);
return 0;
}
which I'm compiling with gcc fft.c -std=c99 -g -lfftw3 -lm
. However the code doesn't appear to be working: FT output is all zeros (see here)
Documentation for the function is here.
EDIT: update, so it only seems to work for FFTW_ESTIMATE and not any of the other flags. Any idea why this might be?