Newsgroups: sci.math.num-analysis,sci.math,comp.lang.c,comp.lang.c++,comp.dsp
Path: usc!howland.reston.ans.net!zaphod.mps.ohio-state.edu!sample.eng.ohio-state.edu!purdue!yuma!hercules.lance.colostate.edu!pm860605
From: pm860605@longs.LANCE.ColoState.Edu (Peter J. McKinney)
Subject: SUMMARY Repost: FFT in C (Best?,Fastest?) Part 2/3: Codes part A
Sender: news@yuma.ACNS.ColoState.EDU (News Account)
Message-ID: <pm860605.130.732479052@longs.LANCE.ColoState.Edu>
Date: Thu, 18 Mar 1993 18:24:12 GMT
Nntp-Posting-Host: hercules.lance.colostate.edu
Organization: Colorado State University
Lines: 1273
Xref: usc sci.math.num-analysis:7465 sci.math:41863 comp.lang.c:67170 comp.lang.c++:39582 comp.dsp:5659


   This is part two of a three part summary to the net concerning FFT code 
in C.  It is a repost of the orininal two part posting (part two was to 
large and wasnot completely posted).   My Apologies.
   I had asked for the best, fastest implementation of the FFT in C or for
a "standard" piece of code if such a thing existed.  Part one of the summary 
is a compilation of the replies I recieved that point out vaious places to 
get FFT code on the internet.  Parts two and three are a compilation the 
actual codes that were sent to me via email.  Some of part two and all of 
part three is in uuencoded, compressed, tar form because that is the way I 
recieved it.  Enough of the header is left from each message so that replies 
may be made to the original poster.  Thanks to all who replied.

- Pete
____________________________________________________________________________
| Peter J. McKinney                     pm860605@longs.LANCE.ColoState.Edu |
| Electrohydrodynamic Laboratory                                           |
| Fluid Mechanics and Wind Engineering Program                             |
| Civil Engineering Department                                             |
| Colorado State University                                                |
----------------------------------------------------------------------------

*******************************************************************************
Return-Path: <herman%phage.nosc.mil@nosc.mil>
~From: herman%phage.nosc.mil@nosc.mil (John W. Herman)
~Subject: Re: Recommendations for good (Fast) FFT code in C ?


/***************************************************************
 *
 * a duhamel-holman split radix dif fft
 * ref:  electronics letters, jan. 5, 1984
 *
 * complex input and output data in double
 * precision complex array c.  length is n.
 *
 ***************************************************************/

#include <math.h>
#include <stddef.h>

/*********************************************************************/

typedef struct {double r, i;} complex;

#define   cmov(a,b)     {(a).r = (b).r; (a).i = (b).i;}

#define   cadd(a,b,c)   {(a).r = (b).r + (c).r; (a).i = (b).i + (c).i;}

#define   csub(a,b,c)   {(a).r = (b).r - (c).r; (a).i = (b).i - (c).i;}

#define   cmul(a,b,c)   {complex tt; \
  tt.r = ((b).r * (c).r) - ((b).i * (c).i); \
  tt.i = ((b).r * (c).i) + ((b).i * (c).r); \
   (a).r = tt.r; (a).i = tt.i;}

#define   conjg(a,b)    {(a).r = (b).r; (a).i = -(b).i;}

#define   rcmul(a,b,c)  {(a).r = b * (c).r; (a).i = b * (c).i;}

/*********************************************************************/

#define TWOPI (M_PI + M_PI)

int srfft(c, n)

complex c[ ];
int n;
{

 register int i0, i1, i2, i3;
 register int n1, n2, n4;
 register int i, j, k, m, id, isxx;

 static double a, e, a3, r1, r2, s1, s2, s3, cc1, ss1, cc3, ss3;

/* compute power of 2 corresponding to n */


 i = 2;
 m = 1;
 while (i < n)
  {
  i = 2 * i;
  m = m + 1;
  };

 n = i;
 n1 = n - 1;
 n2 = 2 * n;
 for(k = 1; k < m; k++)
  {
  n2 = n2 / 2;
  n4 = n2 / 4;
  e = TWOPI / (double)n2;
  a = 0.;
  for(j = 1; j <= n4; j++)
   {
   a3 = a + a + a;
   cc1 = cos(a);
   ss1 = sin(a);
   cc3 = cos(a3);
   ss3 = sin(a3);
   a = (double)j * e;
   isxx = j;
   id = 2 * n2;
   while (isxx < n)
    {
    for(i0 = isxx - 1; i0 < n1 ; i0 = i0 + id)
     {
     i1 = i0 + n4;
     i2 = i1 + n4;
     i3 = i2 + n4;
     r1 = c[i0].r - c[i2].r;
     c[i0].r = c[i0].r + c[i2].r;
     r2 = c[i1].r - c[i3].r;
     c[i1].r = c[i1].r + c[i3].r;
     s1 = c[i0].i - c[i2].i;
     c[i0].i = c[i0].i + c[i2].i;
     s2 = c[i1].i - c[i3].i;
     c[i1].i = c[i1].i + c[i3].i;

     s3 = r1 - s2;
     r1 = r1 + s2;
     s2 = r2 - s1;
     r2 = r2 + s1;

     c[i2].r = r1 * cc1 - s2 * ss1;
     c[i2].i = -s2 * cc1 - r1 * ss1;
     c[i3].r = s3 * cc3 + r2 * ss3;
     c[i3].i = r2 * cc3 - s3 * ss3;

     };

    isxx = 2 * id - n2 + j;
    id = 4 * id;
    };
   };
  };

/* last stage, length = 2 butterfly */


 isxx = 1;
 id = 4;

 while (isxx < n)
  {

  for(i1 = isxx; i1 <= n; i1 = i1 + id)
   {
   i0 = i1 - 1;
   r1 = c[i0].r;
   c[i0].r = r1 + c[i1].r;
   c[i1].r = r1 - c[i1].r;
   r1 = c[i0].i;
   c[i0].i = r1 + c[i1].i;
   c[i1].i = r1 - c[i1].i;
   };

 isxx = 2 * id - 1;
 id = 4 * id;
 };


/* bit reverse counter */


 j = 1;
 for(i = 1; i <= n1; i++)
  {
  if( i < j)
   {
   i0 = i - 1;
   i1 = j - 1;
   r1 = c[i0].r;
   c[i0].r = c[i1].r;
   c[i1].r = r1;

   r1 = c[i0].i;
   c[i0].i = c[i1].i;
   c[i1].i = r1;
   };

  k = n / 2;

  while (k < j)
   {
   j = j - k;
   k = k / 2;
   };

  j = j + k;
  };

 return(n);
}
*******************************************************************************
Return-Path: <gaw-a@minster.york.ac.uk>
~From: gaw-a@minster.york.ac.uk
~Subject: Re: Recommendations for good (Fast) FFT code in C ?

In article <pm860605.108.731192932@longs.LANCE.ColoState.Edu> you write:
> 
>    I'm looking for c code to do an FFT on a unix box.  Is there an
> established "standard" fast code avaible on the net?  If this is an FAQ,
> please point me to the FAQ.  Otherwise I'd like suggestions as to the best
> code available - there seems to be quite a bit to choose from.  I'll check
> the news or my mail for replies and post a summary to the net if there is
> interest.  Thanks,
Hmm... I asked much the same question on our local net, and I got this reply

I haven't managed to find the time to look at it yet, but here's the thing
in its entirety.

If you get any decent replies to this, can you mail them to me, too? Thanx.

---------------------------------------------------------------------------
>From mjl-b Tue Feb 16 09:42:35 1993
~Subject: FFT code
To: gaw-a
X-Mailer: Ream 4.12b York Edition
Status: OR

Here's some code I pulled off news:


/*
 C Program: Test_Fast_Fourier_Transform_in_Double_Precision (TFFTDP.c)
 by: Dave Edelblute, edelblut@cod.nosc.mil, 05 Jan 1993
*/

#include <stdio.h>
#include <math.h>

/***************************************************************/
/* Timer options. You MUST uncomment one of the options below  */
/* or compile, for example, with the '-DUNIX' option.          */
/*                                                             */
/* Example compilation:                                        */
/* UNIX systems:                                               */
/* cc -DUNIX -O tfftdp.c -o fft -lm                            */
/***************************************************************/
/* #define Amiga       */
/* #define UNIX        */
/* #define UNIX_Old    */
/* #define VMS         */
/* #define BORLAND_C   */
/* #define MSC         */
/* #define MAC         */
/* #define IPSC        */
/* #define FORTRAN_SEC */
/* #define GTODay      */
/* #define CTimer      */
/* #define UXPM        */

#ifndef PI
#define PI  3.14159265358979323846
#endif

double xr[262144],xi[262144];

void main()
{
int i,j,np,npm,lp,n2,kr,ki;
double a,enp,t,x,y,zr,zi,zm;
double t1,t2,t3,benchtime,VAX_REF,dtime();
int dfft();

np = 8;
      printf("\nFFT benchmark - Double Precision - V1.0 -  05 Jan 1993\n\n");
      printf("FFT size   Time(sec)   max error\n");

      t3 = dtime();
      for (lp = 1; lp <= 15; lp++ ){ 
   np = np * 2; 
   printf(" %7d ",np);
   enp = np; 
   npm = np / 2  - 1;  
   t = PI / enp;
   xr[0] = (enp - 1.0) * 0.5;  
   xi[0] = 0;
   n2 = np / 2;  
   xr[n2] = -0.5; 
   xi[n2] = 0.0;
   
   for (i = 1; i <= npm; i++) {
       j = np - i;
       xr[i] = -0.5; 
       xr[j] = -0.5;
       y = t * i;  
       y = -0.5*(cos(y)/sin(y));
       xi[i] = y; xi[j] = -y;  }
   
   t1 = dtime();
   dfft(xr,xi,np);
   t2 = dtime() - t1;
   if (t2 < 0.0) t2 = 0.0;
   printf(" %10.4lf ",t2);
/*
   for (i=0; i<=15; i++) printf("%d  %g  %g\n",i,xr[i],xi[i]);
*/
   zr = 0.0; zi = 0.0; kr = 0; ki = 0; npm = np-1;
   for (i = 0; i <= npm; i++ ) { 
       a = fabs( xr[i] - i );
       if (zr < a) { zr = a; kr = i; }
       a = fabs(xi[i]);
       if (zi < a) { zi = a; ki = i; }
       }
   zm = zr;
   if ( fabs(zr) < fabs(zi) ) zm = zi;
   printf(" %10.2g %10.4\n",zm,t2);

/*
   printf("re %7d  %10.2g  im %7d  %10.2g\n",kr,zr,ki,zi);
*/
   }
   benchtime = dtime() - t3;
   VAX_REF = 140.658;
   printf("\nBenchTime (sec) = %10.4lf\n",benchtime);
   printf("VAX_FFTs = %10.3lf\n\n",VAX_REF/benchtime);

}



/*
       A Duhamel-Hollman split-radix dif fft
       Ref: Electronics Letters, Jan. 5, 1984
       Complex input and output data in arrays x and y
       Length is n.
*/

int dfft( x, y, np )
double x[2]; double y[2]; int np ;
{
double *px,*py;
int i,j,k,m,n,i0,i1,i2,i3,is,id,n1,n2,n4 ;
double  a,e,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3,xt ;
  px = x - 1; py = y - 1;
  i = 2; m = 1; while (i < np) { i = i+i; m = m+1; };
  n = i; if (n != np) {
  for (i = np+1; i <= n; i++)  { *(px + i) = 0.0; *(py + i) = 0.0; };
  printf("\nuse %d point fft",n); }

  n2 = n+n;
  for (k = 1;  k <= m-1; k++ ) {
    n2 = n2 / 2; n4 = n2 / 4; e = 2.0 * PI / n2; a = 0.0;
    for (j = 1; j<= n4 ; j++) {
      a3 = 3.0*a; cc1 = cos(a); ss1 = sin(a);
      cc3 = cos(a3); ss3 = sin(a3); a = j*e; is = j; id = 2*n2;
   while ( is < n ) {
   for (i0 = is; i0 <= n-1; i0 = i0 + id) {
      i1 = i0 + n4; i2 = i1 + n4; i3 = i2 + n4;
      r1 = *(px+i0) - *(px+i2);
      *(px+i0) = *(px+i0) + *(px+i2);
      r2 = *(px+i1) - *(px+i3);
      *(px+i1) = *(px+i1) + *(px+i3);
      s1 = *(py+i0) - *(py+i2);
      *(py+i0) = *(py+i0) + *(py+i2);
      s2 = *(py+i1) - *(py+i3);
      *(py+i1) = *(py+i1) + *(py+i3);
      s3 = r1 - s2; r1 = r1 + s2; s2 = r2 - s1; r2 = r2 + s1;
      *(px+i2) = r1*cc1 - s2*ss1; 
   *(py+i2) = -s2*cc1 - r1*ss1;
      *(px+i3) = s3*cc3 + r2*ss3;
       *(py+i3) = r2*cc3 - s3*ss3;
  }
       is = 2*id - n2 + j; id = 4*id;
 }
    }
  }

/*
---------------------Last stage, length=2 butterfly---------------------
*/
  is = 1; id = 4;
  while ( is < n) {
  for (i0 = is; i0 <= n; i0 = i0 + id) {
      i1 = i0 + 1; r1 = *(px+i0);
      *(px+i0) = r1 + *(px+i1);
      *(px+i1) = r1 - *(px+i1);
      r1 = *(py+i0);
      *(py+i0) = r1 + *(py+i1);
      *(py+i1) = r1 - *(py+i1);
    }
  is = 2*id - 1; id = 4 * id; }

/*
c--------------------------Bit reverse counter
*/
  j = 1; n1 = n - 1;
  for (i = 1; i <= n1; i++) {
    if (i < j) {
   xt = *(px+j);
      *(px+j) = *(px+i); *(px+i) = xt;
      xt = *(py+j); *(py+j) = *(py+i);
      *(py+i) = xt;
    }
    k = n / 2; while (k < j) { j = j - k; k = k / 2; }
    j = j + k;
  }

/*
  for (i = 1; i<=16; i++) printf("%d  %g   %gn",i,*(px+i),(py+i));
*/
  
  return(n);
  }


/*****************************************************/
/* Various timer routines.                           */
/* Al Aburto, aburto@marlin.nosc.mil, 20 Dec 1992    */
/*                                                   */
/* t = dtime() outputs the current time in seconds.  */
/* Use CAUTION as some of these routines will mess   */
/* up when timing across the hour mark!!!            */
/*                                                   */
/* For timing I use the 'user' time whenever         */
/* possible. Using 'user+sys' time is a separate     */
/* issue.                                            */
/*                                                   */
/* Example Usage:                                    */
/* [timer options added here]                        */
/* main()                                            */
/* {                                                 */
/* double starttime,benchtime,dtime();               */
/*                                                   */
/* starttime = dtime();                              */ 
/* [routine to time]                                 */
/* benchtime = dtime() - starttime;                  */
/* }                                                 */
/*                                                   */
/* [timer code below added here]                     */
/*****************************************************/

/*********************************/
/* Timer code.                   */
/*********************************/
/*******************/
/*  Amiga dtime()  */
/*******************/
#ifdef Amiga
#include <ctype.h>
#define HZ 50

double dtime()
{
   double q;

   struct   tt {
      long  days;
      long  minutes;
      long  ticks;
   } tt;

   DateStamp(&tt);

   q = ((double)(tt.ticks + (tt.minutes * 60L * 50L))) / (double)HZ;

   return q;
}
#endif

/*****************************************************/
/*  UNIX dtime(). This is the preferred UNIX timer.  */
/*  Provided by: Markku Kolkka, mk59200@cc.tut.fi    */
/*  HP-UX Addition by: Bo Thide', bt@irfu.se         */
/*****************************************************/
#ifdef UNIX
#include <sys/time.h>
#include <sys/resource.h>

#ifdef hpux
#include <sys/syscall.h>
#define getrusage(a,b) syscall(SYS_getrusage,a,b)
#endif

struct rusage rusage;

double dtime()
{
   double q;

   getrusage(RUSAGE_SELF,&rusage);

   q = (double)(rusage.ru_utime.tv_sec);
   q = q + (double)(rusage.ru_utime.tv_usec) * 1.0e-06;
   
   return q;
}
#endif

/***************************************************/
/*  UNIX_Old dtime(). This is the old UNIX timer.  */
/*  Use only if absolutely necessary as HZ may be  */
/*  ill defined on your system.                    */
/***************************************************/
#ifdef UNIX_Old
#include <sys/types.h>
#include <sys/times.h>
#include <sys/param.h>

#ifndef HZ
#define HZ 60
#endif

struct tms tms;

double dtime()
{
   double q;

   times(&tms);

   q = (double)(tms.tms_utime) / (double)HZ;
   
   return q;
}
#endif

/*********************************************************/
/*  VMS dtime() for VMS systems.                         */
/*  Provided by: RAMO@uvphys.phys.UVic.CA                */
/*  Some people have run into problems with this timer.  */
/*********************************************************/
#ifdef VMS
#include time

#ifndef HZ
#define HZ 100
#endif

struct tbuffer_t
       {
 int proc_user_time;
 int proc_system_time;
 int child_user_time;
 int child_system_time;
       };
struct tbuffer_t tms;

double dtime()
{
   double q;

   times(&tms);

   q = (double)(tms.proc_user_time) / (double)HZ;
   
   return q;
}
#endif

/******************************/
/*  BORLAND C dtime() for DOS */
/******************************/
#ifdef BORLAND_C
#include <ctype.h>
#include <dos.h>
#include <time.h>

#define HZ 100
struct time tnow;

double dtime()
{
   double q;

   gettime(&tnow);

   q = 60.0 * (double)(tnow.ti_min);
   q = q + (double)(tnow.ti_sec);
   q = q + (double)(tnow.ti_hund)/(double)HZ;
   
   return q;
}
#endif

/**************************************/
/*  Microsoft C (MSC) dtime() for DOS */
/**************************************/
#ifdef MSC
#include <time.h>
#include <ctype.h>

#define HZ CLK_TCK
clock_t tnow;

double dtime()
{
   double q;

   tnow = clock();

   q = (double)tnow / (double)HZ;
   
   return q;
}
#endif

/*************************************/
/*  Macintosh (MAC) Think C dtime()  */
/*************************************/
#ifdef MAC
#include <time.h>

#define HZ 60

double dtime()
{
   double q;

   q = (double)clock() / (double)HZ;
   
   return q;
}
#endif

/************************************************************/
/*  iPSC/860 (IPSC) dtime() for i860.                       */
/*  Provided by: Dan Yergeau, yergeau@gloworm.Stanford.EDU  */
/************************************************************/
#ifdef IPSC
extern double dclock();

double dtime()
{
   double q;

   q = dclock();
   
   return q;
}
#endif

/**************************************************/
/*  FORTRAN dtime() for Cray type systems.        */
/*  This is the preferred timer for Cray systems. */
/**************************************************/
#ifdef FORTRAN_SEC

fortran double second();

double dtime()
{
   double q;

   second(&q);
   
   return q;
}
#endif

/***********************************************************/
/*  UNICOS C dtime() for Cray UNICOS systems.  Don't use   */
/*  unless absolutely necessary as returned time includes  */
/*  'user+system' time.  Provided by: R. Mike Dority,      */
/*  dority@craysea.cray.com                                */
/***********************************************************/
#ifdef CTimer
#include <time.h>

double dtime()
{
   double    q;
   clock_t   t;

       t = clock();

       q = (double)t / (double)CLOCKS_PER_SEC;

       return q;
}
#endif

/********************************************/
/* Another UNIX timer using gettimeofday(). */
/* However, getrusage() is preferred.       */
/********************************************/
#ifdef GTODay
#include <sys/time.h>

struct timeval tnow;

double dtime()
{
   double q;

   gettimeofday(&tnow,NULL);
   q = (double)tnow.tv_sec + (double)tnow.tv_usec * 1.0e-6;

   return q;
}
#endif

/*****************************************************/
/*  Fujitsu UXP/M timer.                             */
/*  Provided by: Mathew Lim, ANUSF, M.Lim@anu.edu.au */
/*****************************************************/
#ifdef UXPM
#include <sys/types.h>
#include <sys/timesu.h>
struct tmsu rusage;

double dtime()
{
   double q;

   timesu(&rusage);

   q = (double)(rusage.tms_utime) * 1.0e-06;
   
   return q;
}
#endif


| Mathew Lodge             | "What's your name, boy?" "Kate." "Isn't that |
| mjl-b@minster.york.ac.uk |  a bit of a girl's name?" "It's short for... |
| University of York, UK   |  Bob." -- Blackadder II     |
---------------------------------------------------------------------------
Geoff

Strange how the worst loneliness is when you're surrounded by
the best of your friends...

These opinions are my own. If anyone agrees, more fool them.



*******************************************************************************
Return-Path: <p4tb18@park1.park.se>
~Subject: simplefft.doc
~From: Mikael Olsson <p4tb18@park1.park.se>


FFT.C - Spectrum analyzation of audio input.
Author : Steve Sampson, 1984.
Modifyed: Mikael Olsson, march 1993.

Usage :

FFT rate length in_file out_file

The argument "rate" is actually not necessary, it is only used to
give you a clear view in the spectrum analyzer. Set it to the same
as "length" if you're not interested.

"Length" has got to be a power of 2, i.e. 2,4,8,16,32,64,128,256,512,1024...
This is a must for the FFT to work.

"in_file" is a file of raw data in 8-byte DOUBLE form. 1024 samples would
the be a 8192 bytes long file to be read by the program. This may be
modifyed according to your demands.

"out_file" will be created by the program, containing the Real and
Imaginary components of the Components. These are the values that
the program uses to display the spectrum table. The structure is
"length" doubles of Real data, and then "length" doubles of Imaginary 
data.


I made the modifications of removing imaginary data input and enhancing
the spectrum analyzer output since I think this is what people are
usually interested in.  

   /Mikael Olsson      Email:p4tb18@park1.park.se

*******************************************************************************
Return-Path: <p4tb18@park1.park.se>
~Subject: simplefft.c
~From: Mikael Olsson <p4tb18@park1.park.se>


/*
 * fft.c
 *
 * Unix Version 2.4 by Steve Sampson, Public Domain, September 1988
 *
 * This program produces a Frequency Domain display from the Time Domain
 * data input; using the Fast Fourier Transform.
 *
 * Usage: fft rate samples input output
 *  1. sample rate
 *  2. samples is a variable power of two.
 *  3. Input is (samples * sizeof(double)) characters.
 *  4. Output is (samples * sizeof(double)) characters.
 *  5. Standard error is help or debugging messages.
 *
 *      Program modified and simplified by Mikael Olsson at 1993-03-10
 *      added: Phase offset calculation
 *             Enhanced output routine
 *             + Speeded the whole thing up app. 2 times.
 *        Email address : p4tb18@park1.park.se
 *
 */

/* Includes */



#include <stdio.h>
#ifdef AZTEC_C
#include <strings.h>
#else
#include <malloc.h>
#endif
#include <math.h>

/* Defines */

#define pi      3.14159265358979324
#define TWO_PI (2.0 * 3.14159265358979324) /* alias 360 deg */
#define Chunk (Samples * sizeof(double))

#define sine(x)  Sine[(x)]
#define cosine(x) Sine[((x) + (Samples >> 2)) % Samples]

/* Globals, Forward declarations */

static int Samples, Power;
static double *Real, *Imag, Maxn, rate, magnitude(), phase();
static void scale(), fft(), max_amp(), display(), err_report();

static int permute();
static double *Sine;
static void build_trig();

static FILE *Fpi, *Fpo;

extern double atof();

/* The program */

main(argc, argv)
int argc;
char **argv;
{
    int loop;

 if (argc != 5)
  err_report(0);

 rate = atof( *++argv );
 Samples = abs(atoi(*++argv));
 Power = log10((double)Samples) / log10(2.0);

 /* Allocate memory for the variable arrays */

 if ((Real = (double *)malloc(Chunk)) == NULL)
  err_report(1);

 if ((Imag = (double *)malloc(Chunk)) == NULL)
  err_report(2);

 if ((Sine = (double *)malloc(Chunk)) == NULL)
  err_report(3);

 /* open the data files */

 if ((Fpi = fopen(*++argv, "rb")) == NULL)
  err_report(4);

 if ((Fpo = fopen(*++argv, "wb")) == NULL)
  err_report(5);

 /* read in the data from the input file */

 fread(Real, sizeof(double), Samples, Fpi);
 fclose(Fpi);

        /* Zero out imaginary numbers for pure audio data input*/ 
        for(loop=0;loop<Samples;loop++)
         Imag[loop]=0;

 build_trig();
 scale();
 fft();

 display();

 /* write the frequency domain data to the standard output */

 fwrite(Real, sizeof(double), Samples, Fpo);
 fwrite(Imag, sizeof(double), Samples, Fpo);
 fclose(Fpo);

 /* free up memory and let's get back to our favorite shell */

 free((char *)Real);
 free((char *)Imag);
 free((char *)Sine);

 exit(0);
}


static void err_report(n)
int n;
{
 switch (n)  {
 case 0:
  fprintf(stderr, "Usage: fft rate samples in_file out_file\n");
  fprintf(stderr, "Where samples is a power of two\n");
  break;
 case 1:
  fprintf(stderr, "fft: Out of memory getting real space\n");
  break;
 case 2:
  fprintf(stderr, "fft: Out of memory getting imag space\n");
  free((char *)Real);
  break;
 case 3:
  fprintf(stderr, "fft: Out of memory getting sine space\n");
  free((char *)Real);
  free((char *)Imag);
  break;
 case 4:
  fprintf(stderr,"fft: Unable to open data input file\n");
  free((char *)Real);
  free((char *)Imag);
  free((char *)Sine);
  break;
 case 5:
  fprintf(stderr,"fft: Unable to open data output file\n");
  fclose(Fpi);
  free((char *)Real);
  free((char *)Imag);
  free((char *)Sine);
  break;
 }

 exit(1);
}


static void scale()
{
 register int loop;

 for (loop = 0; loop < Samples; loop++)  {
  Real[loop] /= Samples;
  Imag[loop] /= Samples;
 }
}


static void fft()
{
 register int loop, loop1, loop2;
 unsigned i1;   /* going to right shift this */
 int  i2, i3, i4, y, n, t;
 double  a1, a2, b1, b2, z1, z2;

 i1 = Samples >> 1;
 i2 = 1;

 /* perform the butterfly's */

 for (loop = 0; loop < Power; loop++)
 {
  i3 = 0;
  i4 = i1;

  for (loop1 = 0; loop1 < i2; loop1++)
  {
   y = permute(i3 / (int)i1);

   z1 =  cosine(y);
   z2 = -sine(y);

   for (loop2 = i3; loop2 < i4; loop2++)
   {

    a1 = Real[loop2];
    a2 = Imag[loop2];

    b1  = z1*Real[loop2+i1] - z2*Imag[loop2+i1];
    b2  = z2*Real[loop2+i1] + z1*Imag[loop2+i1];

    Real[loop2] = a1 + b1;
    Imag[loop2] = a2 + b2;

    Real[loop2+i1] = a1 - b1;
    Imag[loop2+i1] = a2 - b2;
   }

   i3 += (i1 << 1);
   i4 += (i1 << 1);
  }

  i1 >>= 1;
  i2 <<= 1;
 }
}


/*
 * Display the frequency domain.
 *
 */

static void display()
{
 register int loop, offset;
 int  n, x, i;
 double range,step;
 double p,a;

 max_amp();

 step=rate/Samples;
   range=step*(int)(rate/2/step);

 printf("\nDC Component : %10.3lf\n",magnitude(0));
   printf(  "Freq. Range 0- %10.3lf\n",range);
 printf(  "Freq. Step   : %10.3lf\n",step);
 printf("\n  Freq.  Phase   Amplitude\n");

   n=Samples>>1; i=n;

 for (loop = 0; loop < n; loop++, i++)  {
  a = magnitude(loop);
  x = (int)(a * 40.0 / Maxn);

      if(x>=1)  /* Don't care about the tiny ones... */
                /* May be removed if you want to see ALL frequencies */
   {

    if(loop==0)
     {
    printf("    DC Offset  %8.3lf |", a);
    x/=2;
         }
    else
     {
      p = phase(loop);
      printf("%7.0f %6.3lf %8.3lf |", (0.5 * (i - n) * rate / n), p, a*2);
   }
    offset = 0;
    while (++offset <= x)
     putchar('=');

    putchar('\n');
   }
 }
}

/*
 * Find maximum amplitude
 */

static void max_amp()
{
 register int loop;
 double  mag;

 Maxn = 0.0;
 for (loop = 0; loop < Samples; loop++)  {
  if ((mag = magnitude(loop)) > Maxn)
   Maxn = mag;
 }
}

/*
 * Calculate Power Magnitude
 */

static double magnitude(n)
int n;
{
 n = permute(n);
 return hypot(Real[n], Imag[n]);
}

/*
 * Bit reverse the number
 *
 * Change 11100000b to 00000111b or vice-versa
 */

 /* This was speeded up HEAVILY. Using Logarithms to do this is 
    unnecessary, or not ? =) */
static int permute(int n)
{
  int t;int r=0;
  t=Power;

  while(t--)
   {
  r<<=1;
  r|=(n&1);
  n>>=1;
   }

 return (r);
}

/*
 * Pre-compute the sine and cosine tables
 */

static void build_trig()
{
 register int loop;
 double  angle, increment;

 angle = 0.0;
 increment = TWO_PI / (double)Samples;

 for (loop = 0; loop < Samples; loop++)  {
  Sine[loop] = sin(angle);
  angle += increment;
 }
}

/* Compute the phase offset for two coordinates. */

static double phase (n)
int n;
{
  double x,y,a;

  n = permute(n);
  x=Imag[n];y=Real[n];

  if((x==0)&&(y>=0))
   return(0);
  if((x==0)&&(y<0))
   return(pi);
  if((y==0)&&(x>0))
   return(pi*0.5);
  if((y==0)&&(y<0))
   return(pi*1.5);

  a=atan(x/y);
  if(y<0)a+=pi;
  return(a);
}

*******************************************************************************
Return-Path: <mount@onyx.dartmouth.edu>
~From: mount@onyx.dartmouth.edu (Martin McKinney)

begin 666 FFT.uuenc
M'YV01HQ0^3*%3ATS9EX`6,BPH<.'$"-*G$BQ8D40&&_4J`$"`$:,-6!T_`A"
M!HV1)%.*!%%C!HT;-&;`T&CCX\8;-CQ:W,FSI\^?0(,*'4JTJ-&C2),F#3BP
MX,&$0]ZT@<.F#)XF8>B@<3%&J=>?'VW0..G1ILBR&$VB3(E3QLJ6+V/"L#'C
M!L88,&2TU/FUK]^_@`,+'DRX\`L5/1(K7LRXL>/'D"-+5JS`2)JJ;L*T*:,C
M0=2I5:]FW3I&`9,P;L[4"7.&<P(J:-*X60-BB`(%4][4D3.F#(@Q;\CX-O-&
MSF^I5*V":#,:A)DZ;L;02?/&S9S;38*G,9-F3%;JUG4H2)#WA8P9+W+`2'"D
MC!PY>4`(*1-F=YDS8>20\2R'/ITR9("@51IS.'=9&>.5)\,-Z<7`GGOPR4>?
M??CIE\`4=Z1!QQAH"/B&<VR\D56!9L@A%0ADZ"9&5==-YN*+,,:XF`HOW':8
MC#CF**,"243'1AW"M:CCD$1"1J,"("@P@FQC_"@<"`GPP)Q6+J#A@Y),.EE&
M`B)\EIQH5*(A@HTJM-!";4&008:9;+;IYIMPQBGGG'$JX&5H((2AIH;@*4#G
MGX`&*JB;-":0@)W(X3E$FF2@<*=R8;!06Z+*B9&"`GN,]R@>%];1Q@[CC3>%
MIRY(X5L/>9;JVPH@B*$JJ)VVX4(2;8"`:ABSULJJJ[3"VI]!<K@!PJB?*M`'
MF68.,:H8@S;K+)V(@J;<''6(08<<84C7Y[/<=CMGH8=N6MNRCE**1YZ2BFLI
MIIJ:FP`1VYGA7AG1E0'K>/`B-&^]JMJ:JJD@G.FJJ;#F*V]__-+J+ZX*"YRK
MKV4`*ZS!^_8&ZK$*')9L$W6P08>W(#<;[9?+=3P=5=U]5YV?(;?,+;@C*\JQ
MQ^5*>VZDD]K<ZJ69>N8N%":248=T]R8`='!#T]'OK:JJ,+!O9S+<AM,/CW>T
MT-+ENG#3O.J:*JU4$SS>K[L)>W72%X>J<0M#P&N'RW#+&;-R9*1A!X';QJWW
MGS"+V[;=-9.,L[H\6_T&&_D9?7A^4=3QQG3TT@'KI@DT_G@:D1<-Q>)R6`ZY
M&TIC=8:_FXJ.0A@IO%`Z:RA8"NOFB'?N^.=*0X%&&'.<FO.7MN->QNDI!+Q[
M:+WGWGH*L'J..>C^PIX?%6]LBH+SLE\>.?)C1UPV",IG;JRAR+)=G1I[E]_"
MW.<"YX8:=>#W'\OFQ\]FW^9.NG[@>*+.KL\V^[Q^^UFQ5ZC\QS[WZ2X,DQL?
M`/^C-<0$+`8N@$'1R!8L^Q4P@!<#7\;*)#Z;08]Z\G,6HMQ@!_?,`3S.,5&M
M^B,=U*@F=AX"`1PXYYSBM`%^(30?S*@W/*M\D'/XLTH"\%`XQ<4N`::BEL>*
MED23N4!T'T'5ZLZ``B+"JHD><T'Q=+>I+581>PF@H+"P*#ECA8]ZT-M4#@<U
MPA+*X8352>&)9@A#.GR(A71P8<?R4T,YW'"-\:.?SM`8/7--CX9$W!_ER,C$
M,BA1:0!#%1Z>R!H0J.`X<ZBB%F^7.S"246N2I.3H+GE"-VARBV`4(P@8:<8-
M)LM<H@-DH$96A_\LAS5NT!"0?/,&,QQ'9V[PE!C<@T-9P@UF*:I657HH&BJ*
M*Y$]2^:*MB2Z7!I$.*`:3S5U^214S2$.<J"#)DVE@DD"C%63!%LZVY!*[55P
MF]<4(,;6UD5.EL&8?Z*E+;7B&SC8$P2]S-,O21;,-@Q3#L7$9\B0J:)EUM-W
M03P7-,<CS:H8S9[9O*CO%I9'-\A`D[22E#G+T$Z)@6"+:7-E"Z2@)^@1X3X*
MG5,;W4,'$&"K;J@ID!U1=)_^.#*A,>T60Y7I&Y:2P:7W04%%BZJG-*`&E?M;
MZKMZ6H8RH/0V4SV#3[?HKQG8``81M.0JF_K4?[X`!"B005@O"84DE'1[+]5J
M5:\Z3P[&%7I&#6J=HD)"FO)4KHZ,X4V=:AV@ZM590YTF".[Z!J,JM:&^B>M6
M[5E$J1J5L%<=SV7+NE%4J54D;$V"6"4[5[."P*M@E6#V3+I9-]#U-J$JC&R_
MPA2"&`0A+V`*5V;+V[",92TL.0M)U(*6C\3`)"$!B4M@(A.QU`0CSLU)1WA+
MW>I:][K8S:Y?-EBD[GJW,@?*S&8ZH]O2G"8UJVE-9VQ3&8'DR0T!PMT<WC`&
MI_XG0,^)SG2J,P<78$<[W/'.?L-SF_*4)P<.:L][XC.?^O2G0OO!D(8X!*`8
M5J<%_&R!&=@0'^!,I98J$U9`K[`&-[S`">,QE(H-989LV?&-@@70T'RSF384
M)S[]B4,=TM"?S8"NOPI0$'K4\Z`%2\C!]\G/?H;0GP`&:$`D.I`"O$ME'!V)
MNU7.\F1XY*-="DG+8'[,E9=D!OB6P0R&(@@5B)"$)WP!"5CJLG"B-`<ZU.T-
M5;K2".A5-S/<ALQF1G,"OM"$(%`!"6^.<Y-V&:4I;<5*2N+S=OZ\G4"G>0A/
M:`(4F%`$+!3ZT(E>DIRWU"58CJ9*8]HS?">-I3(+1]!O_@)3%*TE+NE63)%>
MM9^Q'.9>.T8!1[,C'?(`AY_Z^MB.N7*=OS.&!-CA#6G8CQ`T9"HWYBX([PE#
M'B)Z22Q(*D2I`0$3WG`&)V!OV=-I-KC/D(!ITZ':)OS=ND'@!&&ZY]O5&1T4
M$->;.6"/U\A&M@*.$"(QA($-7PZXPI6=QW2C"+()N$(04E.5'HVAQYE+P&%`
M4(7<!4@V`P7Q@`'J2Q*;V`G.<7%Q"G2D)*V-*8<5&13<0QP_%B@,/*WOE%`(
M\A*50<?T&D-\F`)0XR2!Z'%$39ZR'9]>*L##!+7W&UV`$2<\SC=:R8J`T,#+
M6L*AEB`@4-B%)08-V;0,UJYP<80CAQ`H`#:^R8]J?$R'F_='/`G`PA:ZH.*)
M+QU;3?>EN`IZT/YRW.-]'/O7:XJ:`*784+JAP^*I/IYQE]M0EH^AP7-7$I+_
M7=L@.*$>RN#?=_&X#-JJSFNXCJ+3IU[$ON0G"'2;I+=SG?/J6_:/E^Z;.N3N
M.6Q(?.[H,)UPPPN/X,&[$9X@!8E+@0AHA4$*#$6$#]7\#DJ>O4!*GP0G6*$(
M4IA"$=`:@^F_ZT.R2;OVJ>`"#5K][":UP\'KX$BJ@P#N`/4ZV/,#^+`7".2R
M-P>:X1OR)P=.I5BXLW6^(1N39U@Q-RA'\FS1MG[<!@+>!@+S9GE.("FR45/'
MAWH#5D3K9B@%P1KVHF+C,8()8`1L@#ES,'-RP!0[H'&75&\&Y1Z>)P:U]!]R
ML&$N*$,X2"UBL'Z65",)H()3``>H,8.&LG%7$&U:X7FRIX/$1W,<5H0IF&^&
MPA0E&$Y,N'%=6'RC(QO"<2X!)81$>"1'J(4K*!!)F"VR<09?>$GPHGN]T2H1
M<P=5)2S+%DZ>AX9,P7)&J((LF`=1`1V2TX27=(B@@X,U!P*$AX,!184\Z(.!
MI88J"#UPP`)"\#AV]"F*>']O``?O%2!BX(DG0H;*$5`X1XE6&!^8R(9KP(2A
MV!YNX!X')T.[,4.<IXKG$HNI82A.0(LT"`+8UG^BYQO`R&X1UR-E2(L;YXS*
M\8@F=V(I)QTKAX5K&(S-:&9X4!!E``=S"`(6AW',\XB^.'8[=054$`:*I892
MI6)8H(D$4XRHHG>:V`7:&(^&,H^CV"N'$D;W02`\^'!$U8^=^&XGN''WN`4)
M^8GZJ(8/1I`X&(]8D)``F3UG0)'&$8\N5!4SN'$3=P;+-&(E9HTMAHW&(9$#
M66<5"7$15X].>)(HEY(O5H1)(I`;Z9(="9-7`)"',AXHARHHL&XI@`(S=`=I
M%4$L\%A$E0(:&#RLDEHP4'X9M1UHM5BN-W*),7O,YWS0=RF&(G$45P;E6`9T
MYR\M\%DJT%8OX)33E`+#.!YE@'`(,I8C67'188XUA2HKP)9N"9=5(9>P$I`L
MZ((PR!0)@"HHYP,^``(Q$)(U:![:^(;"8BB,"0*."9FPPH5*6%_!R)BP\I/>
M"(YPL)B<"248@54;!P5@5P>DB',M2'S+-`?>L4P[Q1QKL(`U11]S$!^Y61T:
M4ASK=P(%`@<F(E=S((CC@9P=:`8H(`),H0,@(`+($Y"/B`+#DD>MX2\2M)TF
M"`(\(&[D-HP@L`(K4(+=*98)X)R@`YTB8`?6"2ON20?PR05N,)_-:8#O&9TI
MH9\)@!`_,@=H@`)UEDQT\&\!2IS:J9ZZ\YT.*I[D66X[<)[IR9V^(98]TX94
M$(;>69C9R85YY(>H\ITBFA\U-9[F>:)^N`*HXIEPF!KLN:'=6(8?^GB:2!(O
MZH8C6D8J]I`G@BHYRBJ669@+:APH4(B,V)<@8*(<MJ02>IB.E)@"4:'HJ:2Z
M`3KLF0`TF@`?J3MD29)FN9=H&3EB19K/^'@QJ3O`D4E?BCRKN6(_62O>)!NG
M4Y9P6GLKYH]P$$D6N`7Y^"IQ*H^:J#`-&:B]HJ87^3A^JG=`V@9=\"J*BI%T
M^J>/&JF)JJB`.HJ1VJCTN"H6F)!BLZ>;"@>86JE\JC"LLJATD*A[*JJ^$446
M^*G"PZJC*H^4JJ.S^H^U<B:LZJKRZ)"H"*E+$ZJ,6@8J<`4`XZN4FJS`VH_"
M2GQ2<:K^8JO(.J?G::RM.C7**D!JBIZ::*0)@)Z/*JYHJAPN"@+G^HW_(8Z/
M=RR&`J\KV()3ZAZ*F0".B2J121+%Z#;1YAMB$!\RT"HXE@9G@`8U1:#;45-J
M>"%*Z`:&`B7YRIDI!J.@.3J&,K&1F6+K:IH)`"4\P`/Z.H[^^B0!VWDG6Q5F
MD+"QL;)8F&)@B`9O<`<RE)S]L9S:V)[\:9_1J0(`&J`;YGL%>J"1!Z?C`:_U
M>9_YJ:`;AW\'=P;%H2%H4"M5$08EI%.L=[,F,W8"99O88E!5$2!KYQY4!QMB
M!QU>JQG3]'\,:X3N!F]O5`;'J&THX&U1630;5[*^@7;ND0<#$FXG6Y/$>73N
MI898B0(?^'H@$`(O^I5!\'SLF9VT4:(52ALJ6KGH.:./IW=KT*FZP[F>^Y9+
M19B;NP6="TI_>KJT(KJ01;KQ>BB'H;//&9TN4+N6D4L$6F%,@9_SB85]@(+A
M\[9\&[=SFP</""=O%P:[:7?]5Y=ENGN-IXMUIX"U`G)E]S']D79BJQ]DZX#'
M*R<1"&T!(KQI5[P5>($96)XBJ(5+ZAZP,H+D&V_M*P?OJX5SR3]?\AIHZ:X!
M.91H991(.;-+"0--.;I1F:U4695@E)WS^Z&3@H@X>+D/W(ARD+GCL:'Q&[<-
MC"H9;#SSRP)WFV*'V\%ET,"/.;^:JV)4L+_5N@7SVP5&JG<OW,(D_,(Q[)#4
M-KRY,\-"NK^%":_'4E=F0L+?ZR8*4'VM8G;9:T*Y>(;%5E_<46&1:!P!1179
M$EBL*"SKYKU%_"9',F\D7)3Y1F]2AV_AMF]7[&_[,X)I4+_!"+>YDU%P_*"P
MDIUIX,!W/)YHW&\5F@86S*4I-L?^(L@A2[&&(LBH(LA\@%8V>%`@8`)-B@<P
M,,E6F6*-C(.,*76:^9@;FP#PJDISG%)<G$]>4DLTADO<Q$N^!'5X,L6CW,5>
M_`+@LU3,9#K/5%D0!T^[E%&ZW$VA!T[B-%+EU"^L4EM.@5NK$R;:A12^118D
MD5S%15PID18WH5QQT5QU<1=Y(19\L<S>_,W@',[BO%TWHG"]!EZ8,8#K96IA
M8AIZ9(*=8;:S41NW@03TP78&LDR/R,K*X6C.`1VO=QT*D!U]EC(#-@=X)V3I
ML1X*%B$-1B%*QA_^46%0EL]W697FP2`(5F0./2$/%M$$O1V8\V0?XGO#$2(C
M,G8NJ2=JNF()$%!+!63F[&MC5FFO-FB8IFF<YFF&AFAP-@*O)AM;\@4YO6F=
M]FD^328S#68CI'MU-V5+#69C%M2W:#1)8"@H,`,N$`,T$`,UD`,R8`,M40,X
MD`,WD`,S<!XX0`,V<"D`%]5%\G;$!B!GEG!PW5U7-FS%]FJA=RU)$U40-ZI2
ME:A]P$S9I-=T[4MU)@=_'4VYS!JP(E5T=5*<DTWE?-=#`FPF(FQS;=>8K2-7
MYC>,$E'HPDP[`RI^0R[B,CCUXSKH4QLS(TZKG2ZMC3U^XS:DS=HZX]J4PU=J
MD-O8PT/X&QH_%#ND;46O34C2PT/(+553=-S8X]SFXD6W#"J6U5)O$%>"Z1O-
M#7&,Y5BTC-RYUF>P588>"[(BRYFAJ+<$B[+QH;(LN[`O&XI3(+,TBYSD=K/,
M:2A)V[,_*Z!":Z!V5K05>K3[.;LBP+M,>TE.RP90:X!:0;7&?%L)<6OCG!3-
M#%S0/%S./,U>/1/6S%QSP=:^50,TT,T7GN(JON(L_LUO_=E:AL[THLX<6B7N
MC%[P[!G_5=`"!AX(76`R\`('EF`0PF`>G606$M)13-*_<7!CL$>V5(TUJ7(P
M%D>RMV%YT-(J)ALK+;:^!')Z(G_UPK4XET?39'_5=XLQ)(`EI.6&4F,WYE\*
M360-;>1(!F$2[61;)W;<4150#>-A5M.N=F:#5EL_3=6^(6L"D6@O#N@OTM0=
M]=2.#MHUHB2(OH6.^WR&`@.6?F9";2C=]WWA5P2&$@-*/>DXHMF/\SB=_>>H
M+B-7)H'[P13G:\:CHX$<R#R)&X+9I&KD#5LM'NS"/NS$7NS&?NS(GNS*ONS,
MWNS._NS0OA#A%(=C9Z-G6"UI.(AL:+%Q2+($B1IW.$QTH(?TTM<H^H?8'HC:
M2(A/FJ6)6(P-_(A3G(,[:(4_N(R&HHF<.*SCF*/1>XK26KW>Z'FM6.\]>(7X
MG@"SN&(;9XNX&'Q?)P>\N(`#G_##R/"75+RAEP:CM^YLN*[0>$G2>"[42)/7
M^&+[K8(=VZ[C>)9IB8X#WX$?PH[NN$SP")/]^*F2V<+YN(\XGW>%"HH:R9$&
7.4T(>:P[WY"7JHT3R9-%;U%Y1ZDSJ)'P
`
end

*******************************************************************************
____________________________________________________________________________
| Peter J. McKinney                     pm860605@longs.LANCE.ColoState.Edu |
| Electrohydrodynamic Laboratory                                           |
| Fluid Mechanics and Wind Engineering Program                             |
| Civil Engineering Department                                             |
| Colorado State University                                                |
----------------------------------------------------------------------------

