fftw3 example
/***** | |
fftw-3 example | |
Brock Palen | |
brockp@umich.edu | |
http://www.umich.edu/~brockp/ | |
http://cac.engin.umich.edu | |
take the DFT of cos(2*PI*x/10) | |
To not print results compile with -DNDEBUG | |
To change the size of the transform -DSAMPLE=2048 (1024 default) | |
module load fftw | |
icc cosfft.c -I$FFTW_INC -L$FFTW_LINK -lfftw3 -DNDEBUG | |
http://www.fftw.org/doc/One_002dDimensional-DFTs-of-Real-Data.html#One_002dDimensional-DFTs-of-Real-Data | |
*****/ | |
#include <stdio.h> | |
#include <stdlib.h> | |
#include <errno.h> | |
#include <string.h> | |
#include <math.h> | |
#include <fftw3.h> | |
#ifndef M_PI | |
#define M_PI 3.14159265 | |
#endif | |
#ifndef SAMPLE | |
#define SAMPLE 1024 | |
#endif | |
#define WHERESTR "[file %s, line %d]" | |
#define WHEREARG __FILE__,__LINE__ | |
int main(void){ | |
fftw_complex *out; | |
fftw_plan p; | |
double *in, tmp=0; | |
double add=0, mul=0, fma=0; | |
in = (double*) fftw_malloc(sizeof(double)*SAMPLE); | |
out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex)*ceil(SAMPLE/2+1)); /* see note below */ | |
/*p = fftw_plan_dft_r2c_1d(SAMPLE, in, out, FFTW_MEASURE); */ | |
p = fftw_plan_dft_r2c_1d(SAMPLE, in, out, FFTW_ESTIMATE); | |
if(out == NULL || in == NULL){ | |
fprintf(stderr, WHERESTR " ERROR: %s\n", WHEREARG, strerror(errno)); | |
} | |
int i=0; | |
#pragma omp parallel | |
#pragma omp for | |
for (i=0; i<SAMPLE; i++){ | |
in[i]=cos(2.0*M_PI*(double)i/10); | |
} /*end load */ | |
fftw_execute(p); | |
/* transforms of real data have sysmitry. Because of this fftw only computes half the transform (SAMPLE/2+1). | |
This is why 'out' is only that large as well as only printing those values. */ | |
#ifndef NDEBUG | |
for (i=0; i<ceil(SAMPLE/2+1); i++){ | |
tmp = sqrt(out[i][0]*out[i][0]+out[i][1]*out[i][1]); | |
printf("%8.8f\n",tmp); | |
} /* end print */ | |
#endif | |
fftw_flops(p, &add, &mul, &fma); | |
printf("Plan for %d samples used: %f add's, %f multiplies, %f fused Multiply-Add\n" | |
,SAMPLE, add, mul, fma); | |
fftw_destroy_plan(p); | |
fftw_free(in); | |
fftw_free(out); | |
return 0; | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment