/* * This file contains a Fortran to C translation of the 1D transformations * based on the original FFTPACK, written by paul n swarztrauber * at the national center for atmospheric research and available * at www.netlib.org. FFTPACK is in the public domain. * * Higher-dimension transforms written by Erik Lindahl, 2008-2009. * Just as FFTPACK, this file may be redistributed freely, and can be * considered to be in the public domain. * * Any errors in this (threadsafe, but not threaded) C version * are probably due to the f2c translator, or hacks by Erik Lindahl, * rather than FFTPACK. If you find a bug, it would be great if you could * drop a line to lindahl@cbr.su.se and let me know about it! */ #include #include #include #include #include #include "fftpack.h" /** Contents of the FFTPACK fft datatype. * * FFTPACK only does 1d transforms, so we use a pointers to another fft for * the transform in the next dimension. * Thus, a 3d-structure contains a pointer to a 2d one, which in turns contains * a pointer to a 1d. The 1d structure has next==NULL. */ struct fftpack { int ndim; /**< Dimensions, including our subdimensions. */ int n; /**< Number of points in this dimension. */ int ifac[15]; /**< 15 bytes needed for cfft and rfft */ struct fftpack * next; /**< Pointer to next dimension, or NULL. */ RealOpenMM * work; /**< 1st 4n reserved for cfft, 1st 2n for rfft */ }; static void fftpack_passf2(int ido, int l1, RealOpenMM cc[], RealOpenMM ch[], RealOpenMM wa1[], int isign) { int i, k, ah, ac; RealOpenMM ti2, tr2; if (ido <= 2) { for (k=0; k= l1) { for (j=1; j idp) idlj -= idp; war = wa[idlj - 2]; wai = wa[idlj-1]; for (ik=0; ik 5) { wa[i1-1] = wa[i-1]; wa[i1] = wa[i]; } } l1 = l2; } } static int fftpack_transpose_2d(t_complex * in_data, t_complex * out_data, int nx, int ny) { t_complex * src; int i,j; if(nx<2 || ny<2) { if(in_data != out_data) { memcpy(out_data,in_data,sizeof(t_complex)*nx*ny); } return 0; } if(in_data == out_data) { src = (t_complex *)malloc(sizeof(t_complex)*nx*ny); memcpy(src,in_data,sizeof(t_complex)*nx*ny); } else { src = in_data; } for(i=0;inext = NULL; fft->n = nx; /* Need 4*n storage for 1D complex FFT */ if( (fft->work = (RealOpenMM *)malloc(sizeof(RealOpenMM)*(4*nx))) == NULL) { free(fft); return ENOMEM; } if(fft->n>1) fftpack_cffti1(nx,fft->work,fft->ifac); *pfft = fft; return 0; }; int fftpack_init_2d(fftpack_t * pfft, int nx, int ny) { fftpack_t fft; int rc; if(pfft==NULL) { fprintf(stderr,"Fatal error - Invalid FFT opaque type pointer."); return EINVAL; } *pfft = NULL; /* Create the X transform */ if( (rc = fftpack_init_1d(&fft,nx)) != 0) { return rc; } /* Create Y transform as a link from X */ if( (rc=fftpack_init_1d(&(fft->next),ny)) != 0) { free(fft); return rc; } *pfft = fft; return 0; }; int fftpack_init_3d(fftpack_t * pfft, int nx, int ny, int nz) { fftpack_t fft; int rc; if(pfft==NULL) { fprintf(stderr,"Fatal error - Invalid FFT opaque type pointer."); return EINVAL; } *pfft = NULL; /* Create the X transform */ if( (fft = (struct fftpack *)malloc(sizeof(struct fftpack))) == NULL) { return ENOMEM; } fft->n = nx; /* Need 4*nx storage for 1D complex FFT. */ if( (fft->work = (RealOpenMM *)malloc(sizeof(RealOpenMM)*(4*nx))) == NULL) { free(fft); return ENOMEM; } fftpack_cffti1(nx,fft->work,fft->ifac); /* Create 2D Y/Z transforms as a link from X */ if( (rc=fftpack_init_2d(&(fft->next),ny,nz)) != 0) { free(fft); return rc; } *pfft = fft; return 0; }; int fftpack_exec_1d (fftpack_t fft, enum fftpack_direction dir, t_complex * in_data, t_complex * out_data) { int i,n; RealOpenMM * p1; RealOpenMM * p2; n=fft->n; if(n==1) { p1 = (RealOpenMM *)in_data; p2 = (RealOpenMM *)out_data; p2[0] = p1[0]; p2[1] = p1[1]; } /* FFTPACK only does in-place transforms, so emulate out-of-place * by copying data to the output array first. */ if( in_data != out_data ) { p1 = (RealOpenMM *)in_data; p2 = (RealOpenMM *)out_data; /* n complex = 2*n RealOpenMM elements */ for(i=0;i<2*n;i++) { p2[i] = p1[i]; } } /* Elements 0 .. 2*n-1 in work are used for ffac values, * Elements 2*n .. 4*n-1 are internal FFTPACK work space. */ if(dir == FFTPACK_FORWARD) { fftpack_cfftf1(n,(RealOpenMM *)out_data,fft->work+2*n,fft->work,fft->ifac, -1); } else if(dir == FFTPACK_BACKWARD) { fftpack_cfftf1(n,(RealOpenMM *)out_data,fft->work+2*n,fft->work,fft->ifac, 1); } else { fprintf(stderr,"FFT plan mismatch - bad plan or direction."); return EINVAL; } return 0; } int fftpack_exec_2d (fftpack_t fft, enum fftpack_direction dir, t_complex * in_data, t_complex * out_data) { int i,nx,ny; t_complex * data; nx = fft->n; ny = fft->next->n; /* FFTPACK only does in-place transforms, so emulate out-of-place * by copying data to the output array first. * For 2D there is likely enough data to benefit from memcpy(). */ if( in_data != out_data ) { memcpy(out_data,in_data,sizeof(t_complex)*nx*ny); } /* Much easier to do pointer arithmetic when base has the correct type */ data = (t_complex *)out_data; /* y transforms */ for(i=0;inext,dir,data+i*ny,data+i*ny); } /* Transpose in-place to get data in place for x transform now */ fftpack_transpose_2d(data,data,nx,ny); /* x transforms */ for(i=0;in; ny=fft->next->n; nz=fft->next->next->n; /* FFTPACK only does in-place transforms, so emulate out-of-place * by copying data to the output array first. * For 3D there is likely enough data to benefit from memcpy(). */ if( in_data != out_data ) { memcpy(out_data,in_data,sizeof(t_complex)*nx*ny*nz); } /* Much easier to do pointer arithmetic when base has the correct type */ data = (t_complex *)out_data; /* Perform z transforms */ for(i=0;inext->next,dir,data+i*nz,data+i*nz); /* For each X slice, transpose the y & z dimensions inside the slice */ for(i=0;inext,dir,data+i*ny,data+i*ny); } /* Transpose back to (nx,ny,nz) */ for(i=0;iwork); if(fft->next != NULL) fftpack_destroy(fft->next); free(fft); } }