SHOW:
|
|
- or go back to the newest paste.
| 1 | #include <fftw3.h> | |
| 2 | #include <math.h> | |
| 3 | ||
| 4 | ||
| 5 | #define M_PI 3.14159265358979323846264338327 | |
| 6 | ||
| 7 | int main(void) | |
| 8 | {
| |
| 9 | int N = 16; int M = 3; | |
| 10 | ||
| 11 | // input data (pre transform) | |
| 12 | double input[N][M]; | |
| 13 | double *in = *input; | |
| 14 | ||
| 15 | // result of transform (in Fourier space) | |
| 16 | fftw_complex output[N][M]; | |
| 17 | fftw_complex *out = (fftw_complex *) ((void *) (*output)); | |
| 18 | ||
| 19 | // input data (post transform) so we can compute accuracy | |
| 20 | double input2[N][M]; | |
| 21 | double *in2 = *input2; | |
| 22 | ||
| 23 | ||
| 24 | // setup plans | |
| 25 | const int rank = 1; const int howmany = M; | |
| 26 | const int istride = M; const int ostride = M; | |
| 27 | const int idist = 1; const int odist = 1; | |
| 28 | ||
| 29 | fftw_plan fp = fftw_plan_many_dft_r2c(rank, &N, howmany, | |
| 30 | in, NULL, istride, idist, | |
| 31 | out, NULL, ostride, odist, | |
| 32 | FFTW_MEASURE); | |
| 33 | ||
| 34 | fftw_plan bp = fftw_plan_many_dft_c2r(rank, &N, howmany, | |
| 35 | out, NULL, istride, idist, | |
| 36 | in2, NULL, ostride, odist, | |
| 37 | FFTW_MEASURE); | |
| 38 | ||
| 39 | for (int j = 0; j < M; ++j) | |
| 40 | fprintf (stderr, "cos(%d PI x) \t\t\t", 2*(j+1)); | |
| 41 | fprintf (stderr, "\n"); | |
| 42 | ||
| 43 | // setup input data | |
| 44 | for (int i = 0; i < N; ++i) | |
| 45 | {
| |
| 46 | for (int j = 0; j < M; ++j) | |
| 47 | {
| |
| 48 | const double x = (double) i / (double) N; | |
| 49 | input[i][j] = cos(2*(j+1) * M_PI * x); | |
| 50 | fprintf (stderr, "input[%3d][%d] = %.9f\t", i, j, input[i][j]); | |
| 51 | } | |
| 52 | fprintf (stderr, "\n"); | |
| 53 | } | |
| 54 | fprintf (stderr, "\n"); | |
| 55 | ||
| 56 | ||
| 57 | ||
| 58 | // take the forward transform | |
| 59 | fftw_execute(fp); fprintf (stderr, "FFT complete!\n"); | |
| 60 | ||
| 61 | // print output (in Fourier space) | |
| 62 | - | for (int i = 0; i < N/2; ++i) |
| 62 | + | for (int i = 0; i <= N/2; ++i) |
| 63 | {
| |
| 64 | for (int j = 0; j < M; ++j) | |
| 65 | {
| |
| 66 | fprintf (stderr, "OUT[%3d] = (%.3f + %.3fi)\t", | |
| 67 | i, output[i][j][0], output[i][j][1]); | |
| 68 | } | |
| 69 | fprintf (stderr, "\n"); | |
| 70 | } | |
| 71 | ||
| 72 | fprintf (stderr, "\n"); | |
| 73 | ||
| 74 | ||
| 75 | ||
| 76 | fftw_execute(bp); fprintf (stderr, "IFFT complete!\n"); | |
| 77 | ||
| 78 | double maxError = 0.0; | |
| 79 | ||
| 80 | for (int i = 0; i < N; ++i) | |
| 81 | for (int j = 0; j < M; ++j) | |
| 82 | {
| |
| 83 | // remember that FFTW doesn't normalise | |
| 84 | const double error = fabs(input[i][j] - input2[i][j] / (double) N); | |
| 85 | if (error > maxError) maxError = error; | |
| 86 | } | |
| 87 | ||
| 88 | fprintf (stderr, "maximum error after 1 cycle: %.4g\n", maxError); | |
| 89 | ||
| 90 | // clean up | |
| 91 | fftw_destroy_plan(fp); | |
| 92 | fftw_destroy_plan(bp); | |
| 93 | ||
| 94 | return 0; | |
| 95 | } |
