|
Fourier Transformations
In article <6873@pdxgate.UUCP> you write:
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 <stdio.h>
#include <math.h>
/*
* 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
Discuss this article in the forums
See Also: © 1999-2011 Gamedev.net. All rights reserved. Terms of Use Privacy Policy
|