Fourier Transformations by Steve Lamont February 20, 1993 In article <6873@pdxgate.UUCP> you write: >In article lusardi@cs.buffalo.edu (Christopher Lusardi) writes: >>Does anyone have a "simple" program in C (Modula 2, Pascal) that computes >>the 1D or 2D Fourier transform on some data and maybe the inverse of the >>transform? (How do I explictly use your implementation? What data should >>I use on it et cetera?) >> >>I just want to learn the transform; so simplicity in the program is a must >>for me! Would you also send me your mathematical formula for the transform? >>I'll accept anything that is working and is dissectable. > >I'd like to have a copy of that source too. Would it be possible to put it >up for ftp somewhere, assuming it exists (the source, that is). Thanks. I was going to post this but my nntp server seems to have indigestion right now. If you wish to repost the code, you have my permission. This is from Numerical Recipes, adapted from FORTRAN to C. Just for the heck of it, I've also attached the original FORTRAN, since I used it as a check against my C coding. For an explanation of the theory, the code, and how to use it, see Press, _et_al_, _Numerical Recipes_ (ISBN 0-521-30811-9), pp. 449-453. spl ``` #include #include /* * N dimensional DFT, filched from Numerical Recipes and translated * by hand from the original FORTRAN. */ #define PI 3.14159265358979323846 #define TWO_PI ( PI * 2.0 ) void ndfft( double *data, int *nn, int ndim, int isign ) { int ntot = 1; int nprev = 1; int idim; double i2pi = isign * TWO_PI; for ( idim = 0; idim < ndim; idim++ ) ntot *= *( nn + idim ); for ( idim = 0; idim < ndim; idim++ ) { int n = *( nn + idim ); int nrem = ntot / ( n * nprev ); int ip1 = 2 * nprev; int ip2 = ip1 * n; int ip3 = ip2 * nrem; int i2rev = 0; int i2; int ifp1; /* * Bit reversal stuff. */ for ( i2 = 0; i2 < ip2; i2 += ip1 ) { int ibit; if ( i2 < i2rev ) { int i1; for ( i1 = i2; i1 < i2 + ip1; i1 += 2 ) { register int i3; for ( i3 = i1; i3 < ip3; i3 += ip2 ) { register int i3rev = i2rev + i3 - i2; register double tempr = *( data + i3 ); register double tempi = *( data + i3 + 1 ); *( data + i3 ) = *( data + i3rev ); *( data + i3 + 1 ) = *( data + i3rev + 1 ); *( data + i3rev ) = tempr; *( data + i3rev + 1 ) = tempi; } } } ibit = ip2 / 2; while ( ( ibit > ip1 ) && ( i2rev > ibit - 1 ) ) { i2rev -= ibit; ibit /= 2; } i2rev += ibit; } /* * Danielson-Lanczos stuff. */ ifp1 = ip1; while ( ifp1 < ip2 ) { int ifp2 = 2 * ifp1; double theta = i2pi / ( ( double ) ifp2 / ip1 ); register double wpr; register double wpi; register double wr = 1.0; register double wi = 0.0; int i3; wpr = sin( 0.5 * theta ); wpr *= wpr * -2.0; wpi = sin( theta ); for ( i3 = 0; i3 < ifp1; i3 += ip1 ) { register int i1; register double wtemp; for ( i1 = i3; i1 < i3 + ip1; i1 += 2 ) { register int i2; for ( i2 = i1; i2 < ip3; i2 += ifp2 ) { register int i21 = i2 + 1; register int k2 = i2 + ifp1; register int k21 = k2 + 1; register double tempr = ( wr * *( data + k2 ) ) - ( wi * *( data + k21 ) ); register double tempi = ( wr * *( data + k21 ) ) + ( wi * *( data + k2 ) ); *( data + k2 ) = *( data + i2 ) - tempr; *( data + k21 ) = *( data + i21 ) - tempi; *( data + i2 ) += tempr; *( data + i21 ) += tempi; } } wtemp = wr; wr += ( wr * wpr ) - ( wi * wpi ); wi += ( wi * wpr ) + ( wtemp * wpi ); } ifp1 = ifp2; } nprev *= n; } } ----- cut here ------ subroutine fourn( data, nn, ndim, isign ) implicit double precision ( a-h, o-z ) double precision data(*) integer nn(ndim) integer ndim integer isign ntot = 1 do idim = 1, ndim ntot = ntot * nn(idim) enddo nprev = 1 do idim = 1, ndim n = nn(idim) nrem = ntot / ( n * nprev ) ip1 = 2 * nprev ip2 = ip1 * n ip3 = ip2 * nrem i2rev = 1 do i2 = 1, ip2, ip1 if ( i2 .lt. i2rev ) then do i1 = i2, i2 + ip1 - 2, 2 do i3 = i1, ip3, ip2 i3rev = i2rev + i3 - i2 tempr = data(i3) tempi = data(i3 + 1) data(i3) = data(i3rev) data(i3 + 1) = data(i3rev + 1) data(i3rev) = tempr data(i3rev + 1) = tempi enddo enddo endif ibit = ip2 / 2 1 continue if ( ( ibit .ge. ip1 ) .and. ( i2rev .gt. ibit ) ) then i2rev = i2rev - ibit ibit = ibit / 2 go to 1 endif i2rev = i2rev + ibit enddo ifp1 = ip1 2 continue if ( ifp1 .lt. ip2 ) then ifp2 = 2 * ifp1 theta = isign * 3.14159265358979323846d0 * 2.0 / | ( ifp2 / ip1 ) wpr = -2.0 * sin( 0.5d0 * theta )**2 wpi = sin( theta ) wr = 1.0 wi = 0.0 do i3 = 1, ifp1, ip1 do i1 = i3, i3 + ip1 - 2, 2 do i2 = i1, ip3, ifp2 k1 = i2 k2 = k1 + ifp1 tempr = ( wr * data(k2) ) - ( wi * data(k2 + 1) ) tempi = ( wr * data(k2 + 1) ) + | ( wi * data(k2) ) data(k2) = data(k1) - tempr data(k2 + 1) = data(k1 + 1) - tempi data(k1) = data(k1) + tempr data(k1 + 1) = data(k1 + 1) + tempi enddo enddo wtemp = wr wr = ( wr * wpr ) - ( wi * wpi ) + wr wi = ( wi * wpr ) + ( wtemp * wpi ) + wi enddo ifp1 = ifp2 go to 2 endif nprev = n * nprev enddo return end ``` Steve Lamont, SciViGuy -- (619) 534-7968 -- spl@szechuan.ucsd.edu UCSD Microscopy and Imaging Resource/UCSD Med School/La Jolla, CA 92093-0608 "maggotbox /mag'et-boks/ n. See Macintrash. This is even more derogatory."     - The New Hacker's Dictionary, Eric Raymond, ed. Discuss this article in the forums Date this article was posted to GameDev.net: 7/16/1999 (Note that this date does not necessarily correspond to the date the article was written) See Also: Fourier Transforms © 1999-2011 Gamedev.net. All rights reserved. Terms of Use Privacy Policy Comments? Questions? Feedback? Click here!