(c)  Copyright 1989 Commodore-Amiga, Inc.   All rights reserved.
The information contained herein is subject to change without notice, and 
is provided "as is" without warranty of any kind, either expressed or implied.  
The entire risk as to the use of this information is assumed by the user.




                    Amiga Math: Fast Fourier Transform

                             by Dan Baker



     The Fourier transform is a numerical process which is used widely in
digital signal analysis.  The fast Fourier transform or FFT is a special
algorithm which can be used to speed up the transform.

     The example code listed below shows one method for doing an FFT on the
Amiga.  The example takes a 256-point sample and computes the FFT.  The
magnitude of the first 20 harmonics computed with the FFT are displayed.

     The example program can be compiled in several ways.  If you want to 
use the standard, single-precision math package that comes with Lattice or 
Manx, set the pre-processor flag NORMAL to a one.  To compile the program 
under Lattice C (V4.01), use lc with no flags and link with c.o, lcm.lib, 
lc.lib and amiga.lib.  To compile the program under Manx C (V3.4a), use 
cc +L and link with ln -lm32 -lc32.  

     If you want to use the IEEE double precision math libraries with
Lattice C (V4.01) then set the pre-processor flag IEEE_LATTICE to a one
and set NORMAL to zero.  Then compile the program with lc (no flags) and 
link with c.o and these libraries:

     mathieeedoubbas_lib.lib
     mathieeedoubtrans_lib.lib
     lcm.lib
     lc.lib
     amiga.lib

     If you compile two versions of the program, one with single-precision
and the other with double precision, you will not see any difference in the
results - you have to go out to the fourteenth decimal place before the extra
precision of the IEEE libraries has an affect on the calculations!

     There is one other way to compile the FFT program and that is under DEC 
Ultrix (V2.2).  Set the pre-processor flag NORMAL to a one and IEEE_LATTICE 
to zero.  Replace the first line of the program with #include "math.h".  
Use cc fft.c -lm to compile the program.

     For more information on FFTs and their application, see the article
"Fast Fourier Transforms on Your Home Computer" by Stanley and Peterson from
the Byte Book of Computer Music (C. P. Morgan, ed.)



----------------------- Code Starts Here ---------------------

#include "exec/types.h"
/*==================================*/
/* 256 Point Fast Fourier Transform */
/*==================================*/


#define N     256                /* Number of points in the signal sample */
#define NPOW  8                  /* Number of points as a power of 2      */
#define PI    3.14159            /* A definition of PI                    */


/*==============================*/
/* Compile Flags - set only one */
/*==============================*/                           
#define NORMAL       1           /* Set this to a one to use the standard */
                                 /* math package under Manx or Lattice    */

#define IEEE_LATTICE 0           /* Set this to a one to use the double   */
                                 /* precision IEEE library with Lattice   */
                                 /* (Set only one of these flags)         */

long scram();

#if NORMAL
/*==================================*/
/* Transcendental Math Declarations */
/*==================================*/

#define SQRT      sqrt        /* We use the standard Unix(tm) function  */
#define POW       pow         /* names for the transcendental math      */
#define SIN       sin         /* routines used in the FFT.              */
#define COS       cos

double sqrt();                /* Under both Manx and Lattice, these      */
double pow();                 /* functions take and return doubles.      */
double sin();
double cos();
/*==========================*/
#endif NORMAL




#if IEEE_LATTICE
/*==========================*/
/* IEEE Transcendental Math */
/*==========================*/
#define SQRT     IEEEDPSqrt       /* Here are the IEEE double precison     */
#define POW(x,y) IEEEDPPow(y,x)   /* math fuction names. The parameters    */
#define SIN      IEEEDPSin        /* for the power routine are backwards.  */
#define COS      IEEEDPCos

double IEEEDPSqrt();            /* The IEEE transcendental functions    */
double IEEEDPPow();             /* take doubles as arguments and return */
double IEEEDPSin();             /* doubles.                             */
double IEEEDPCos();

extern ULONG MathIeeeDoubBasBase;      /* Library bases */
       ULONG MathIeeeDoubTransBase;
/*==========================*/
#endif IEEE_LATTICE


void
main()
{                                /* Signal arrays.  X1 is the original   */
double  x1[256];                 /* signal.  When done, X1 has the real  */
double  x2[256];                 /* part of the transform and X2 has the */
                                 /* imaginary part.                      */

double  icrou[128];        /* Imaginary part of the complex roots of unity */
double  rcrou[128];        /* Real part of the complex roots of unity      */
double  mag  [128];        /* Magnitude of the freq. spectrum components   */

double a1,a2,b1,b2,i1,i2,i3,i4,  /* These floats hold intermediate values */
       z1,z2,max,n,v,z;

long   i,j,k,m,y;                /* Loop counters, array subscripts */

#if IEEE_LATTICE
/*=====================*/
/* Open Math Libraries */
/*=====================*/

MathIeeeDoubBasBase=OpenLibrary("mathieeedoubbas.library",0);
if(MathIeeeDoubBasBase==0)exit(0);
MathIeeeDoubTransBase=OpenLibrary("mathieeedoubtrans.library",0);
if(MathIeeeDoubTransBase==0)
  {
  printf("Can't open transcendental library\n");
  CloseLibrary(MathIeeeDoubBasBase);
  exit(0);
  }
#endif IEEE_LATTICE

printf("Calculating complex roots of unity...\n");
n=N;
v=2*PI/n;                          /* 256 points sampled over one cycle */
for(i=0;i<n/2;i++)
  {
    /*====================*/
    /* Calculate complex  */
    /* roots of unity     */
    /*====================*/
    z=v*i;                          /* Now fill our arrays for the real  */
    rcrou[i]=  COS(z);              /* and imaginary parts of the 128th  */
    icrou[i]= -SIN(z);              /* roots of unity.                   */
  }

/*=============================*/
/*   Do the Set Up, 25% Duty   */
/*   Cycle Square Wave Pulse   */
/*=============================*/

printf("Initializing wave table...\n");

for(i=0;i<(n/4);i++)x1[i]=1.0;     /* X1 contains the original sampled   */
for(i=n/4;i<n;i++)  x1[i]=0.0;     /* signal to be transformed. A square */
                                   /* wave is used here for an example.  */
for(i=0;i<n;i++)
   {
   x1[i]=x1[i]/n;                 /* Now scale the input signal and set */
   x2[i]=0.0;                     /* X2[] to zeroes                     */
   }

/*===================*/
/*   Calculate FFT   */
/*===================*/
printf("FFT in progress...\n");
i1=n/2.0;
i2=1.0;
for(i=1;i<NPOW+1;i++)
  {                                      /*  Here we use the efficient FFT*/
  i3=0.0;                                /* algorithm which takes on the  */
  i4=i1;                                 /* order of NlogN computations   */
  for( k=1 ; (k < (i2+1)) ; k++ )        /* instead of N*N computations.  */
    {                                    /* 
    j=scram((double)(long)((i3/i1)+0.5));/* (Truncation by casting to long)*/
                                         /*  A point in the transform is   */
    z1=rcrou[j];                         /* computed from the even and odd */
    z2=icrou[j];                         /* points of the signal times an  */     
                                         /* appopriate root of unity.      */
    for(m=i3;m<i4;m++)                   /*  The FFT processing results in */
      {                                  /* a scrambled order for spectrum */
      a1=x1[m];                          /* componenmts.  Scram() is used  */
      a2=x2[m];
      b1=z1 * x1[m+(long)i1]-z2 * x2[m+(long)i1];
      b2=z2 * x1[m+(long)i1]+z1 * x2[m+(long)i1];
      x1[m]=a1+b1;
      x2[m]=a2+b2;                       /* to get the appropriate root of*/
      x1[m+(long)i1]=a1-b1;              /* unity and to unscramble the   */
      x2[m+(long)i1]=a2-b2;              /* components when the FFT is    */
      }                                  /* finished.                     */
    i3=i3+2.0*i1;
    i4=i4+2.0*i1;
    }
  i1=i1/2.0;
  i2=i2*2.0;
  }
/*======================*/
/* Calculate magnitudes */
/*======================*/
max=0.0;
printf("Caluculating magnitudes...\n");  /* This calculates magnitude and */
for(z=0;z<n/2;z++)                       /* finds the largest component   */
  {                                      /* -max- for scaling purposes.   */
  y=scram(z);
  mag[y]=SQRT( POW( x1[y],2.0 ) + POW( x2[y],2.0  ) );
  if (mag[y]>max) max=mag[y];
  }
                                      /* This shows the first 19 values */
/*===========================*/       /* in the frequency spectrum.     */
/*   Show tabular results    */       /* There are 128 all together.    */
/*===========================*/
printf("Harmonic      Real        Imaginary      Magnitude\n");
for(z=0;z<20;z++)
  {
  y=scram(z);
  printf("%3d %13.6e %13.6e %13.6e\n",(long)z,x1[y],x2[y],mag[y]);
  }

/*============================*/
/*  Show the graphic results  */
/*============================*/
printf("\nHarmonic      Magnitude\n");
for(z=0;z<20;z++)
  {
  printf("%2ld ",(long)z);             /* This gives a poor man's graphic */
  y=scram(z);                          /* display of the first 19 values  */
  i=(long)((56.0*mag[y]/max)+0.5);     /* in the frequency spectrum.      */
  for(j=1;j<i;j++)printf("=");         /* Plenty of room for improvement  */
  printf("\n");                        /* here!                           */
  }

#if IEEE_LATTICE
/*============================*/
/*  Close the Math Libraries  */
/*============================*/
CloseLibrary(MathIeeeDoubTransBase);
CloseLibrary(MathIeeeDoubBasBase);
#endif IEEE_LATTICE

}

/*=======================*/
/*  Scramble/Unscramble  */
/*=======================*/
long scram(xx)
double xx;                      /* Since the FFT algorithm results in */
{                               /* a scrambled order for the freq.    */
long yy=0;                      /* components, scram is used to fetch */
long n1=N;                      /* the appropriate root of unity      */
long w;                         /* at calculation time and at the end */
for(w=1;w<(NPOW+1);w++)         /* to unscramble the results.         */

    {
    n1=n1/2;
    if(!( (long)xx < n1 ))
      {
      yy = yy + (long)POW( 2.0,(double)w-1 );
      xx=xx-n1;
      }
    }

return(yy);
}



