Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <fftw3.h>
- #include <math.h>
- #define M_PI 3.14159265358979323846264338327
- int main(void)
- {
- int N = 16; int M = 3;
- // input data (pre transform)
- double input[N][M];
- double *in = *input;
- // result of transform (in Fourier space)
- fftw_complex output[N][M];
- fftw_complex *out = (fftw_complex *) ((void *) (*output));
- // input data (post transform) so we can compute accuracy
- double input2[N][M];
- double *in2 = *input2;
- // setup plans
- const int rank = 1; const int howmany = M;
- const int istride = M; const int ostride = M;
- const int idist = 1; const int odist = 1;
- fftw_plan fp = fftw_plan_many_dft_r2c(rank, &N, howmany,
- in, NULL, istride, idist,
- out, NULL, ostride, odist,
- FFTW_MEASURE);
- fftw_plan bp = fftw_plan_many_dft_c2r(rank, &N, howmany,
- out, NULL, istride, idist,
- in2, NULL, ostride, odist,
- FFTW_MEASURE);
- for (int j = 0; j < M; ++j)
- fprintf (stderr, "cos(%d PI x) \t\t\t", 2*(j+1));
- fprintf (stderr, "\n");
- // setup input data
- for (int i = 0; i < N; ++i)
- {
- for (int j = 0; j < M; ++j)
- {
- const double x = (double) i / (double) N;
- input[i][j] = cos(2*(j+1) * M_PI * x);
- fprintf (stderr, "input[%3d][%d] = %.9f\t", i, j, input[i][j]);
- }
- fprintf (stderr, "\n");
- }
- fprintf (stderr, "\n");
- // take the forward transform
- fftw_execute(fp); fprintf (stderr, "FFT complete!\n");
- // print output (in Fourier space)
- for (int i = 0; i <= N/2; ++i)
- {
- for (int j = 0; j < M; ++j)
- {
- fprintf (stderr, "OUT[%3d] = (%.3f + %.3fi)\t",
- i, output[i][j][0], output[i][j][1]);
- }
- fprintf (stderr, "\n");
- }
- fprintf (stderr, "\n");
- fftw_execute(bp); fprintf (stderr, "IFFT complete!\n");
- double maxError = 0.0;
- for (int i = 0; i < N; ++i)
- for (int j = 0; j < M; ++j)
- {
- // remember that FFTW doesn't normalise
- const double error = fabs(input[i][j] - input2[i][j] / (double) N);
- if (error > maxError) maxError = error;
- }
- fprintf (stderr, "maximum error after 1 cycle: %.4g\n", maxError);
- // clean up
- fftw_destroy_plan(fp);
- fftw_destroy_plan(bp);
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement