FFT for Real Data in Two and Three Dimensions.

1. FFT for Real Data in Two and Three Dimensions. Sung-Ju Kang Department of Physics Kangwon National University Two-dimensional FFT are particularly important in the field of image processing. An image is usually represented as a two-dimensional array of pixel intensities, real numbers. One commonly desires to filter high, or low, frequency spatial components from an image; or to convolve or deconvolve the image with some instrumental point spread function. Use of the FFT is by far the most efficient technique. In three dimensions, a common use of the FFT is to solve Poisson’s equation for a potential on a three-dimensional lattice that represents the discretization of three-dimensional space. Here the source terms and the desired potentials are also real.

2. Programming with C language for (i2=1;i2<=nn2;i2++) { if (i3 == 1) { j2=(i2 != 1 ((nn2-i2)<<1)+3 : 1); h1r=c1*(data[i1][i2][1]+speq[j1][j2]); h1i=c1*(data[i1][i2][2]-speq[j1][j2+1]); h2i=c2*(data[i1][i2][1]-speq[j1][j2]); h2r= -c2*(data[i1][i2][2]+speq[j1][j2+1]); data[i1][i2][1]=h1r+h2r; data[i1][i2][2]=h1i+h2i; speq[j1][j2]=h1r-h2r; speq[j1][j2+1]=h2i-h1i; } else { j2=(i2 != 1 nn2-i2+2 : 1); j3=nn3+3-(i3<<1); h1r=c1*(data[i1][i2][ii3]+data[j1][j2][j3]); h1i=c1*(data[i1][i2][ii3+1]-data[j1][j2][j3+1]); h2i=c2*(data[i1][i2][ii3]-data[j1][j2][j3]); h2r= -c2*(data[i1][i2][ii3+1]+data[j1][j2][j3+1]); data[i1][i2][ii3]=h1r+wr*h2r-wi*h2i; data[i1][i2][ii3+1]=h1i+wr*h2i+wi*h2r; data[j1][j2][j3]=h1r-wr*h2r+wi*h2i; data[j1][j2][j3+1]= -h1i+wr*h2i+wi*h2r;}} wr=(wtemp=wr)*wpr-wi*wpi+wr; wi=wi*wpr+wtemp*wpi+wi;}} if (isign == -1) fourn(&data[1][1][1]-1,nn,3,isign); } #include <meth.h> void rlft3(float ***data, float **speq, unsigned long nn1, unsigned long nn2,unsigned long nn3, int isign) { void fourn(float data[], unsigned long nn[], int ndim, int isign); void nrerror(char error_text[]); unsigned long i1,i2,i3,j1,j2,j3,nn[4],ii3; double theta,wi,wpi,wpr,wr,wtemp; float c1,c2,h1r,h1i,h2r,h2i; if (1+&data[nn1][nn2][nn3]-&data[1][1][1] != nn1*nn2*nn3) nrerror("rlft3: problem with dimensions or contiguity \n"); c1 = 0.5; c2 = -0.5*isign; theta=isign*(6.28318530717959/nn3); wtemp=sin(0.5*theta); wpr = -2.0*wtemp*wtemp; wpi=sin(theta); nn[1]=nn1; nn[2]=nn2; nn[3]=nn3 >> 1; if (isign == 1) { fourn(&data[1][1][1]-1,nn,3,isign); for (i1=1;i1<=nn1;i1++) for (i2=1,j2=0;i2<=nn2;i2++) { speq[i1][++j2]=data[i1][i2][1]; speq[i1][++j2]=data[i1][i2][2]; } } for (i1=1;i1<=nn1;i1++) { j1=(i1 != 1 nn1-i1+2 : 1); wr=1.0; wi=0.0; for (ii3=1,i3=1;i3<=(nn3>>2)+1;i3++,ii3+=2) { Now give some fragments from notional calling programs, to clarify the use of rlft3 for two- and three-dimensional data.

3. Note again that the routine does not actually distinguish between two and three dimensions; Two is treated like three, but with the first dimension having length 1. The first program fragment FFTs a two-dimensional data array, allows for some processing on it, e.g., filtering, and then takes the inverse transform. The second program example illustrates a three-dimensional transform, where the three dimensions have different lengths. #include <stdlib.h> #include "nrutil.h" #define N2 256 #define N3 256 int main(void) /* First Program */ { void rlft3(float ***data, float **speq, unsigned long nn1, unsigned long nn2, unsigned long nn3, int isign); float ***data, **speq; data=f3tensor(1,1,1,N2,1,N3); speq=matrix(1,1,1,2*N2); rlft3(data,speq,1,N2,N3,1); rlft3(data,speq,1,N2,N3,-1); free_matrix(speq,1,1,1,2*N2); free_f3tensor(data,1,1,1,N2,1,N3); return 0; } #define N1 32 #define N2 64 #define N3 16 int main(void) /* Second Program */ { void rlft3(float ***data, float **speq, unsigned long nn1, unsigned long nn2, unsigned long nn3, int isign); int j; float ***data,**speq; data=f3tensor(1,N1,1,N2,1,N3); speq=matrix(1,N1,1,2*N2); rlft3(data,speq,N1,N2,N3,1); free_matrix(speq,1,N1,1,2*N2); free_f3tensor(data,1,N1,1,N2,1,N3); return 0; }