Relay-Version: version B 2.10 5/3/83; site utzoo.UUCP
Path: utzoo!mnetor!uunet!seismo!lll-lcc!lll-tis!ames!sdcsvax!ucsdhub!hp-sdd!hplabs!hpcea!wunder
From: wunder@hpcea.HP.COM (Walter Underwood)
Newsgroups: comp.sources.d
Subject: Re: Re: fast fourier transform
Message-ID: <630001@hpcea.HP.COM>
Date: Wed, 15-Jul-87 21:23:27 EDT
Article-I.D.: hpcea.630001
Posted: Wed Jul 15 21:23:27 1987
Date-Received: Sat, 18-Jul-87 05:51:01 EDT
References: <407@uhnix2.UUCP>
Organization: HP Corporate Engineering - Palo Alto, CA
Lines: 155


   / 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);
}