mfinegan@uceng.UC.EDU (michael k finegan) (07/30/90)
3ksnn64@cadlab1.ecn.purdue.edu (Joe Cychosz) writes: >In article <1990Jul22.175529.8882@midway.uchicago.edu> lipman@fciads.bsd.uchicago.edu (Everett Lipman) writes: >>I am looking for source code for a 2 dimensional Fast Fourier >>Transform routine. If anyone has one or knows where to get >>one, please send me e-mail. C is preferred, but any language >>is fine. Thanks. >> >> Everett Lipman (lipman@fciads.bsd.uchicago.edu) lots of code deleted ... Since the original poster wanted C code, here is a 2D fft that works with 2d float arrays (where each column is 2*rows w/ real,imag,real,imag ... format). While verbose compared to the fortran version - fast, even though unoptimized. A very basic test routine is also included ... - Mike mfinegan@uceng.UC.EDU /*----------------------- cut here for file 1 of 2 --------------------------*/ /* * testfft.c * test 2D fft subroutine, using trivial example * compiled under SysV using: cc testfft.c twodfft.c -lm -o testfft */ #include <stdio.h> float **testarray; main() { int i, j; /* allocate 2D float array - use size = power of 2 */ if((testarray = (float **)malloc(sizeof(float *)*8)) == (float **)NULL) { fprintf(stderr,"Error: can't allocate float**\n"); exit(1); } for(j=0;j<8;j++) { if((testarray[j] = (float *)malloc(sizeof(float)*2*8)) == (float *)NULL) { fprintf(stderr,"Error: can't allocate float* #%d\n",j); exit(1); } } /* try a real, constant, input */ for(i=0;i<8;i++) for(j=0;j<8;j++) { testarray[i][j*2] = 1.0; testarray[i][j*2 + 1] = 0.0; } /* show contents ... */ printf("\n input data:\n"); for(i=0;i<8;i++) { for(j=0;j<8;j++) { printf("%+05.1f%+05.1f", testarray[i][j*2],testarray[i][j*2 + 1]); } printf("\n"); } twodfft(8,8,(float **)testarray,-1); /* show forward (uncentered) FFT */ printf("\n <CR> to display forward results ..."); getc(stdin); printf("\n"); for(i=0;i<8;i++) { for(j=0;j<8;j++) { printf("%+05.1f%+05.1f", testarray[i][j*2],testarray[i][j*2 + 1]); } printf("\n"); } twodfft(8,8,(float **)testarray,1); /* show twice trasformed (recentered) aproximation to original */ printf("\n <CR> to display inverse results (== original) ..."); fgetc(stdin); printf("\n"); for(i=0;i<8;i++) { for(j=0;j<8;j++) { printf("%+05.1f%+05.1f", testarray[i][j*2],testarray[i][j*2 + 1]); } printf("\n"); } } /*----------------------- cut here for file 2 of 2 --------------------------*/ /* * twodfft.c */ #include <stdio.h> #include <math.h> /* for MSC V5.1 - use float.h and huge or large model */ #define TRUE 1 #define FALSE 0 /* * Perform 2-D FFT by calling multiple vector * subroutine twice: call once for columns, * transpose results, then call again for * rows. Result will be FFT, with D.C. @ 0,0. * * NOTE: Expected: M,N power of 2, and M = N. * The caller can zeropad to a square array, * or could be added here. (I do it when allocating * the 2d array ...). Forward transform if sign = -1 * (= sign of exponent), else inverse transform is * assumed with the inverse result being divided through * by M*N. * * ALL THREE subroutines are used. (Also note inclusion of matherr()) * * 2D complex array = 2D float array with twice as many columns as rows. * It is ASSUMED that the array is dynamically allocated (some compilers * have problems accessing static 2d arrays in the subroutine as **array. * - internal structure of array must be different! */ int twodfft(nrows, ncols, array, sign) int nrows, ncols; float **array; int sign; { int i, j; int offset; float normalization; void transpose(); /* if inverse FFT, conjugate */ if(sign > 0) { for(i = 0L; i < nrows; i++) for(j = 0L; j < ncols; j++) { array[i][j*2L + 1L] *= -1.0; } } if(!(ffmvtno(nrows, ncols, array))) return(FALSE); transpose(nrows, ncols, array); if(!(ffmvtno(nrows, ncols, array))) return(FALSE); /* divide inverse transform by nrows*ncols */ if(sign > 0) { normalization = 1.0/((float)nrows*(float)ncols); for(i = 0L; i < nrows; i++) for(j = 0L; j < ncols; j++) { array[i][j*2L] *= normalization; array[i][j*2L + 1L] *= normalization; } } return(TRUE); } /* * Fast Fourier Multiple Vector Transform, with 'Natural Output' * derived from fortran subroutine 'ffatno' in appendix of * "An Efficient Two-Dimensional FFT Algorithm", L. R. Johnson and * A. K. Jain, IEEE T-PAMI Vol. PAMI-3, No. 6, pp. 698-701, 1981. * * ffatno - possibly stands for: * Fast Fourier Algorithm multiple vector Transform Natural Order * Fast Fourier Algorithm decimation in Time Natural Order * Fast Fourier Algorithm Two-dimensional Natural Order * * At any rate, converted to C from Fortran - all indices shifted down * by one (1 to N becomes 0 to N - 1), complex replaced by float, where * the even elements are real and the odd are imaginary. Complex addition * and subtraction are done component by component, while multiplication * and exponentiation are done explicitly, using components. A 2-D array, * 'a', is implemented as a column order 2-D array; elements a[0][0], a[0][1] * are the real, and imaginary, parts of the first element in the array. * * The "shuffling" was added at the beginning of the routine, as mentioned * in the paper, and its appendix. The other alternative would be to write * two versions of the subroutine, one producing shuffled output and one * expecting shuffled input, but producing "natural ordered" output. * * One final implementation detail: * The Fortran subroutine uses a 'trick', which I assume is legitimate * Fortran (it compiles fine with one Unix F77, but not with another). * The trick is the access of a 2-D array as an 1-D array, using row indices * that exceed the 2-D range, but stay within a corresponding 1-D range. * While convenient, and more efficient than copying 2-D arrays into * intermediate 1-D arrays, etc., it should have been documented. * (Authors have assumed it to be a commonly used technique ?). */ /* * Fast Fourier Multiple Vector Transform w/ Natural Output, * from subroutine 'ffatno' (c) 1978 by Lawerence R. Johnson ; * changes (c) 1989 by Michael K. Finegan, Jr. */ int bitrev[1024]; float row_copy[1024]; int ffmvtno(nrows, ncols, array) int nrows, ncols; float **array; /* origin of 2-D 'complex' array */ { int i, j, k, l, m; int npoint, nexpon, ifirst, isecnd, eq_ifirst_row, eq_ifirst_col; int eq_isecnd_row, eq_isecnd_col; float w[2L], temp[2L], first[2L], secnd[2L]; /* * create bit reverse 'look-up table'; taken from (2-D via sequential 1-D) * another 'C' FFT program (Purdue origin? - cf. Gonzalez & Wintz D.I.P.). */ m = (int)(floor(log((double)nrows)/log(2.0))); for(i=0L; i < nrows; i++) { k = 0L; l = i; for(j=0L; j < m; j++) { k = k*2L + l%2L; l = l/2L; } bitrev[i] = k; } /* * shuffle each column using 'bit reversed' row indices */ for(j=0L; j < ncols; j++) { for(i=0L; i < nrows; i++) { row_copy[i*2L] = array[bitrev[i]][j*2L]; row_copy[i*2L + 1L] = array[bitrev[i]][j*2L + 1L]; } for(i=0L; i < nrows; i++) { array[i][j*2L] = row_copy[i*2L]; array[i][j*2L + 1L] = row_copy[i*2L + 1L]; } } /* * start of 'C' coded ffatno subroutine: */ npoint = 1L; do { npoint = npoint * 2L; for(nexpon = 0L; nexpon < npoint/2L; nexpon++) { ifirst = nexpon; /* start @ array[0] */ isecnd = ifirst + npoint/2L; /* * 'C' implementation of W = cexp(Z): * 1) exp(j*theta) = cos(theta) + j*sin(theta) (Euler's formula) * 2) exp(Z1 + Z2) = exp(Z1) * exp(Z2) * let Z = Z1 + Z2, Z1 = X + j0, Z2 = 0 + jY, and we have: * W = cexp(X+jY) = cexp(X) * cexp(jY) = exp(X) * (cos(Y) + jsin(Y)) * with X = 0.0, W = 1.0 * (cos(Y) + jsin(Y)) = cos(Y) + j*sin(Y) * and so W has nonzero real, and imaginary, components. */ w[0L] = (float)cos(-2.0*3.1415926535*(double)nexpon/(double)npoint); w[1L] = (float)sin(-2.0*3.1415926535*(double)nexpon/(double)npoint); for(i=0L; i < nrows*ncols/npoint; i++) { /* NOTE: fortran apparently allows 1D access to 2D arrays!; * as a result, adjust indices to reflect `1D' access. * Since we use col order, more than addition involved * to get row, col equivalent offset. Calculate equivalent * offset once per iteration to save computation ... * eq_offset = (offset to start of row) + (offset to col) */ eq_ifirst_row = ifirst%nrows; eq_ifirst_col = (ifirst/nrows)*2L; eq_isecnd_row = isecnd%nrows; eq_isecnd_col = (isecnd/nrows)*2L; /* temp = complex multiply of w and array[isecnd][0] */ temp[0L] = w[0L]*array[eq_isecnd_row][eq_isecnd_col] - w[1L]*array[eq_isecnd_row][eq_isecnd_col + 1L]; temp[1L] = w[1L]*array[eq_isecnd_row][eq_isecnd_col] - w[0L]*array[eq_isecnd_row][eq_isecnd_col + 1L]; /* complex subtraction */ array[eq_isecnd_row][eq_isecnd_col] = array[eq_ifirst_row][eq_ifirst_col] - temp[0L]; array[eq_isecnd_row][eq_isecnd_col + 1L] = array[eq_ifirst_row][eq_ifirst_col + 1L] - temp[1L]; /* complex addition */ array[eq_ifirst_row][eq_ifirst_col] += temp[0L]; array[eq_ifirst_row][eq_ifirst_col + 1L] += temp[1L]; ifirst = ifirst + npoint; isecnd = isecnd + npoint; } } } while(npoint < nrows); /* return with success */ return(TRUE); } /* * Transpose a square vector, in place. */ void transpose(nrows, ncols, array) int nrows, ncols; float **array; { int i, j; float temp[2L]; /* 0 = real, 1 = imag */ if(nrows != ncols) { fprintf(stderr,"Error: only know how to transpose square matrices!\n"); exit(1); } for(j=0L; j < nrows - 1L; j++) { for(i=j+1L; i < nrows; i++) { /* swap row element with column element (complex) */ temp[0L] = array[i][j*2L]; temp[1L] = array[i][j*2L + 1L]; array[i][j*2L] = array[j][i*2L]; array[i][j*2L + 1L] = array[j][i*2L + 1L]; array[j][i*2L] = temp[0L]; array[j][i*2L + 1L] = temp[1L]; } } } /* * Optional, with need being compiler dependent! */ int matherr(x) struct exception *x; { int return_value = 1; switch(x->type) { case DOMAIN: fprintf(stderr,"MATH: DOMAIN error in %s: %e, %e\n",x->name,x->arg1,x->arg2); break; case SING: fprintf(stderr,"MATH: SING error in %s: %e, %e\n",x->name,x->arg1,x->arg2); break; case OVERFLOW: fprintf(stderr,"MATH: OVERFLOW error in %s: %e, %e\n",x->name,x->arg1,x->arg2); break; case PLOSS: fprintf(stderr,"MATH: PLOSS error in %s: %e, %e\n",x->name,x->arg1,x->arg2); break; case TLOSS: fprintf(stderr,"MATH: TLOSS error in %s: %e, %e\n",x->name,x->arg1,x->arg2); break; case UNDERFLOW: fprintf(stderr,"MATH: UNDERFLOW error in %s: %e, %e\n",x->name,x->arg1,x->arg2); break; default: fprintf(stderr,"MATH: unknown error in %s: %e, %e\n",x->name,x->arg1,x->arg2); return_value = 1; break; } x->retval = 1.0; return(return_value); } /*--------------------------------- EOF ---------------------------------*/