uace0@uhnix2.UUCP (Univ ATARI Comp Enthusiasts) (07/11/87)
In article <1065@bloom-beacon.MIT.EDU> dquah@ares.UUCP (Danny Quah) writes: > This is my third posting on this in as many months; I ... ^^^^^ I've watched all 3 go by with no response as well. So does anyone out there have one? (A fast fourier transform algorithm written in C. -- Actually I could port it from FORTRAN, Danny.) Neal Symms, University of Houston ..!academ!uhnix1!uhnix2!uace0-- >>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<< UACE + A Smith & Wesson beats a four of a kind! + uhnix2!uace0 >>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<
wunder@hpcea.HP.COM (Walter Underwood) (07/16/87)
/ hpcea:comp.sources.d / uace0@uhnix2.UUCP (Univ ATARI Comp Enthusiasts) / 6:04 pm Jul 10, 1987 / In article <1065@bloom-beacon.MIT.EDU> dquah@ares.UUCP (Danny Quah) writes: > This is my third posting on this in as many months; I ... ^^^^^ I've watched all 3 go by with no response as well. So does anyone out there have one? Yes. --------------------------------------- /* Decimation-in-time FFT benchmark based on the Cooley-Tukey algorithm. * This C Language implementation is translated from Rabiner & Gold. * Standard C is used as defined in "The C Programming Language" by * Kernighan & Ritchie. * Paul R. Mouchon (415) 852-4770 * Ford Aerospace MS-X21 * 3939 Fabian Way * Palo Alto, CA 94303 */ #include <stdio.h> #define TWOPI 6.283185308 #define N 512 /* fft size, number of points */ #define NU 9 /* radix 2 exponent, N = 2**NU */ #define NLOOPS 100 /* number of N-point fft's to compute */ main() /* main benchmark driver */ { int i,j; /* loop counter */ double xreal[N],ximag[N]; /* data arrays */ double ccoef[N],scoef[N]; /* cosine & sine coefficient tables */ for ( j=0; j<NLOOPS; j++ ) { for ( i=0; i<N; i++ ) /* initialize data arrays to empty */ { xreal[i] = 4.01*i; ximag[i] = 5.022*i; } for ( i=0; i<N; i++ ) /* init cos & sin coefficient tables */ { ccoef[i] = 1.01*i+.005; scoef[i] = 1.02*i+.004; } fft(xreal,ximag,N,NU,ccoef,scoef); /* radix 2 fft, n=2**nu */ } } /*********************************************************************/ fft(xreal,ximag,n,nu,ccoef,scoef) /* radix 2 fft, n=2**nu */ double xreal[],ximag[]; /* real & imaginary data arrays */ double ccoef[],scoef[]; /* cosine and sine coefficient tables */ int n,nu; /* number of points & radix 2 exponent */ { double nd, /* number of data points as type double */ c,s, /* cosine & sine coefficients */ treal,timag; /* temporary storage */ int l, /* counter for nu in-place computational stages */ n2, /* number dual node pairs in computational stage l */ i, /* loop counter for n2 dual node pair computations */ k, /* node index into data arrays */ p, /* index into first node of a dual node pair */ kp, /* index into corresponding dual node of pair */ nu1; /* nu-1 */ nd = (double) n; n2 = n/2; nu1 = nu-1; for ( l=1; l<=nu; l++, nu1-=1, n2/=2 ) { for ( k=0; k<n; k+=n2 ) { for ( i=1; i<=n2; i++, k++ ) { p = bitrev((k/power(2,nu1)),nu); c = ccoef[p]; s = scoef[p]; kp = k + n2; treal = (xreal[kp] * c) + (ximag[kp] * s); timag = (ximag[kp] * c) - (xreal[kp] * s); xreal[kp] = xreal[k] - treal; ximag[kp] = ximag[k] - timag; xreal[k] = xreal[k] + treal; ximag[k] = ximag[k] + timag; } } } scramrx2(xreal,ximag,n); /* unscramble transformed data */ } /*********************************************************************/ scramrx2(x,y,n) /* radix 2 scrambler-unscrambler */ double x[],y[]; /* data arrays to be reordered */ int n; /* number of data points */ { double tx,ty; /* temporary storage */ register int i,j,k; /* loop counters and array indeces */ for ( i=0, j=0; i<(n-1); i++ ) { if ( i<j ) { tx = x[j]; ty = y[j]; x[j] = x[i]; y[j] = y[i]; x[i] = tx; y[i] = ty; } for ( k=n/2; k<=j; k=k/2 ) j = j-k; j = j+k; } } /*********************************************************************/ int bitrev(j,nbits) /* bit reverse the low nbits of j */ int j,nbits; { int i,j2,br; for ( i=1, br=0; i<=nbits; i++ ) { j2 = j/2; br = (br*2) + (j-(2*j2)); j = j2; } return(br); } /*********************************************************************/ power(x,n) /* raise x to the n-th power, n>0 */ int x,n; { int p; for ( p=1; n>0; --n ) p = p*x; return(p); }
rob@cs.vu.nl (Rob van Leeuwen) (07/16/87)
A friend of mine wrote an fft-package in C. If you're interested, send him a message: valke@cs.vu.nl. He's away for holidays now, but will be back in three weeks or so.
richard@islenet.UUCP (Richard Foulk) (07/17/87)
In article <407@uhnix2.UUCP> uace0@uhnix2.UUCP (Neal Symms) writes: > In article <1065@bloom-beacon.MIT.EDU> dquah@ares.UUCP (Danny Quah) writes: > > This is my third posting on this in as many months; I ... > ^^^^^ > I've watched all 3 go by with no response as well. > So does anyone out there have one? (A fast fourier transform algorithm > written in C. -- Actually I could port it from FORTRAN, Danny.) > Or if someone has a version in BASIC (gross, but I've seen them) I'll be glad to run it through our translator. -- Richard Foulk ...{dual,vortex,ihnp4}!islenet!richard Honolulu, Hawaii
coleman@sask.UUCP (Geoff Coleman @ College of Engineering) (07/20/87)
> Keywords: fft, general radix, C > Xref: sask comp.lang.c:2650 comp.sources.d:863 > Posted: Thu Jul 16 22:13:36 1987 >> > This is my third posting on this in as many months; I ... >> ^^^^^ Ok already I'll dig out a C version of an fft program. I have it stored somewhere on tape. It is a conversion from the fortran fft subroutine that is given in "Digital Image Processing" by Ganzalez and Wintz. Hopefully I'll get it posted to comp.sources tomorrow. -- Geoff Coleman | BITNET: Coleman@sask College of Engineering | UUCP: {utcsri,ihnp4}!sask!skul!geoff University of Saskatchewan | Compserve: 76515,1513 just a number Saskatoon, Saskatchewan | voice: (306) 966-5415
matt@jaws.UChicago.EDU (Matt Crawford) (07/21/87)
In response to those looking for a fast Fourier transformer:
Send a mail message to research!netlib, with a body consisting
of:
send index
or:
send index from fftpack
or include both lines in a single message. You'll get back a
list of available programs (mostly fortran) and a short
description of each. You can ask for any of the programs listed
by sending another message of the form:
send fooprog from barpack
Internet users can send requests to the address
netlib@anl-mcs.arpa.
________________________________________________________
Matt University matt@oddjob.uchicago.edu
Crawford of Chicago {astrovax,ihnp4}!oddjob!matt