home *** CD-ROM | disk | FTP | other *** search
/ Amiga Developer CD v1.2 / amidev_cd_12.iso / reference / amiga_mail_vol1 / math / fftnew next >
Text File  |  1990-01-26  |  11KB  |  288 lines
  1. (c)  Copyright 1989 Commodore-Amiga, Inc.   All rights reserved.
  2. The information contained herein is subject to change without notice, and 
  3. is provided "as is" without warranty of any kind, either expressed or implied.  
  4. The entire risk as to the use of this information is assumed by the user.
  5.  
  6.  
  7.  
  8.  
  9.                     Amiga Math: Fast Fourier Transform
  10.  
  11.                              by Dan Baker
  12.  
  13.  
  14.  
  15.      The Fourier transform is a numerical process which is used widely in
  16. digital signal analysis.  The fast Fourier transform or FFT is a special
  17. algorithm which can be used to speed up the transform.
  18.  
  19.      The example code listed below shows one method for doing an FFT on the
  20. Amiga.  The example takes a 256-point sample and computes the FFT.  The
  21. magnitude of the first 20 harmonics computed with the FFT are displayed.
  22.  
  23.      The example program can be compiled in several ways.  If you want to 
  24. use the standard, single-precision math package that comes with Lattice or 
  25. Manx, set the pre-processor flag NORMAL to a one.  To compile the program 
  26. under Lattice C (V4.01), use lc with no flags and link with c.o, lcm.lib, 
  27. lc.lib and amiga.lib.  To compile the program under Manx C (V3.4a), use 
  28. cc +L and link with ln -lm32 -lc32.  
  29.  
  30.      If you want to use the IEEE double precision math libraries with
  31. Lattice C (V4.01) then set the pre-processor flag IEEE_LATTICE to a one
  32. and set NORMAL to zero.  Then compile the program with lc (no flags) and 
  33. link with c.o and these libraries:
  34.  
  35.      mathieeedoubbas_lib.lib
  36.      mathieeedoubtrans_lib.lib
  37.      lcm.lib
  38.      lc.lib
  39.      amiga.lib
  40.  
  41.      If you compile two versions of the program, one with single-precision
  42. and the other with double precision, you will not see any difference in the
  43. results - you have to go out to the fourteenth decimal place before the extra
  44. precision of the IEEE libraries has an affect on the calculations!
  45.  
  46.      There is one other way to compile the FFT program and that is under DEC 
  47. Ultrix (V2.2).  Set the pre-processor flag NORMAL to a one and IEEE_LATTICE 
  48. to zero.  Replace the first line of the program with #include "math.h".  
  49. Use cc fft.c -lm to compile the program.
  50.  
  51.      For more information on FFTs and their application, see the article
  52. "Fast Fourier Transforms on Your Home Computer" by Stanley and Peterson from
  53. the Byte Book of Computer Music (C. P. Morgan, ed.)
  54.  
  55.  
  56.  
  57. ----------------------- Code Starts Here ---------------------
  58.  
  59. #include "exec/types.h"
  60. /*==================================*/
  61. /* 256 Point Fast Fourier Transform */
  62. /*==================================*/
  63.  
  64.  
  65. #define N     256                /* Number of points in the signal sample */
  66. #define NPOW  8                  /* Number of points as a power of 2      */
  67. #define PI    3.14159            /* A definition of PI                    */
  68.  
  69.  
  70. /*==============================*/
  71. /* Compile Flags - set only one */
  72. /*==============================*/                           
  73. #define NORMAL       1           /* Set this to a one to use the standard */
  74.                                  /* math package under Manx or Lattice    */
  75.  
  76. #define IEEE_LATTICE 0           /* Set this to a one to use the double   */
  77.                                  /* precision IEEE library with Lattice   */
  78.                                  /* (Set only one of these flags)         */
  79.  
  80. long scram();
  81.  
  82. #if NORMAL
  83. /*==================================*/
  84. /* Transcendental Math Declarations */
  85. /*==================================*/
  86.  
  87. #define SQRT      sqrt        /* We use the standard Unix(tm) function  */
  88. #define POW       pow         /* names for the transcendental math      */
  89. #define SIN       sin         /* routines used in the FFT.              */
  90. #define COS       cos
  91.  
  92. double sqrt();                /* Under both Manx and Lattice, these      */
  93. double pow();                 /* functions take and return doubles.      */
  94. double sin();
  95. double cos();
  96. /*==========================*/
  97. #endif NORMAL
  98.  
  99.  
  100.  
  101.  
  102. #if IEEE_LATTICE
  103. /*==========================*/
  104. /* IEEE Transcendental Math */
  105. /*==========================*/
  106. #define SQRT     IEEEDPSqrt       /* Here are the IEEE double precison     */
  107. #define POW(x,y) IEEEDPPow(y,x)   /* math fuction names. The parameters    */
  108. #define SIN      IEEEDPSin        /* for the power routine are backwards.  */
  109. #define COS      IEEEDPCos
  110.  
  111. double IEEEDPSqrt();            /* The IEEE transcendental functions    */
  112. double IEEEDPPow();             /* take doubles as arguments and return */
  113. double IEEEDPSin();             /* doubles.                             */
  114. double IEEEDPCos();
  115.  
  116. extern ULONG MathIeeeDoubBasBase;      /* Library bases */
  117.        ULONG MathIeeeDoubTransBase;
  118. /*==========================*/
  119. #endif IEEE_LATTICE
  120.  
  121.  
  122. void
  123. main()
  124. {                                /* Signal arrays.  X1 is the original   */
  125. double  x1[256];                 /* signal.  When done, X1 has the real  */
  126. double  x2[256];                 /* part of the transform and X2 has the */
  127.                                  /* imaginary part.                      */
  128.  
  129. double  icrou[128];        /* Imaginary part of the complex roots of unity */
  130. double  rcrou[128];        /* Real part of the complex roots of unity      */
  131. double  mag  [128];        /* Magnitude of the freq. spectrum components   */
  132.  
  133. double a1,a2,b1,b2,i1,i2,i3,i4,  /* These floats hold intermediate values */
  134.        z1,z2,max,n,v,z;
  135.  
  136. long   i,j,k,m,y;                /* Loop counters, array subscripts */
  137.  
  138. #if IEEE_LATTICE
  139. /*=====================*/
  140. /* Open Math Libraries */
  141. /*=====================*/
  142.  
  143. MathIeeeDoubBasBase=OpenLibrary("mathieeedoubbas.library",0);
  144. if(MathIeeeDoubBasBase==0)exit(0);
  145. MathIeeeDoubTransBase=OpenLibrary("mathieeedoubtrans.library",0);
  146. if(MathIeeeDoubTransBase==0)
  147.   {
  148.   printf("Can't open transcendental library\n");
  149.   CloseLibrary(MathIeeeDoubBasBase);
  150.   exit(0);
  151.   }
  152. #endif IEEE_LATTICE
  153.  
  154. printf("Calculating complex roots of unity...\n");
  155. n=N;
  156. v=2*PI/n;                          /* 256 points sampled over one cycle */
  157. for(i=0;i<n/2;i++)
  158.   {
  159.     /*====================*/
  160.     /* Calculate complex  */
  161.     /* roots of unity     */
  162.     /*====================*/
  163.     z=v*i;                          /* Now fill our arrays for the real  */
  164.     rcrou[i]=  COS(z);              /* and imaginary parts of the 128th  */
  165.     icrou[i]= -SIN(z);              /* roots of unity.                   */
  166.   }
  167.  
  168. /*=============================*/
  169. /*   Do the Set Up, 25% Duty   */
  170. /*   Cycle Square Wave Pulse   */
  171. /*=============================*/
  172.  
  173. printf("Initializing wave table...\n");
  174.  
  175. for(i=0;i<(n/4);i++)x1[i]=1.0;     /* X1 contains the original sampled   */
  176. for(i=n/4;i<n;i++)  x1[i]=0.0;     /* signal to be transformed. A square */
  177.                                    /* wave is used here for an example.  */
  178. for(i=0;i<n;i++)
  179.    {
  180.    x1[i]=x1[i]/n;                 /* Now scale the input signal and set */
  181.    x2[i]=0.0;                     /* X2[] to zeroes                     */
  182.    }
  183.  
  184. /*===================*/
  185. /*   Calculate FFT   */
  186. /*===================*/
  187. printf("FFT in progress...\n");
  188. i1=n/2.0;
  189. i2=1.0;
  190. for(i=1;i<NPOW+1;i++)
  191.   {                                      /*  Here we use the efficient FFT*/
  192.   i3=0.0;                                /* algorithm which takes on the  */
  193.   i4=i1;                                 /* order of NlogN computations   */
  194.   for( k=1 ; (k < (i2+1)) ; k++ )        /* instead of N*N computations.  */
  195.     {                                    /* 
  196.     j=scram((double)(long)((i3/i1)+0.5));/* (Truncation by casting to long)*/
  197.                                          /*  A point in the transform is   */
  198.     z1=rcrou[j];                         /* computed from the even and odd */
  199.     z2=icrou[j];                         /* points of the signal times an  */     
  200.                                          /* appopriate root of unity.      */
  201.     for(m=i3;m<i4;m++)                   /*  The FFT processing results in */
  202.       {                                  /* a scrambled order for spectrum */
  203.       a1=x1[m];                          /* componenmts.  Scram() is used  */
  204.       a2=x2[m];
  205.       b1=z1 * x1[m+(long)i1]-z2 * x2[m+(long)i1];
  206.       b2=z2 * x1[m+(long)i1]+z1 * x2[m+(long)i1];
  207.       x1[m]=a1+b1;
  208.       x2[m]=a2+b2;                       /* to get the appropriate root of*/
  209.       x1[m+(long)i1]=a1-b1;              /* unity and to unscramble the   */
  210.       x2[m+(long)i1]=a2-b2;              /* components when the FFT is    */
  211.       }                                  /* finished.                     */
  212.     i3=i3+2.0*i1;
  213.     i4=i4+2.0*i1;
  214.     }
  215.   i1=i1/2.0;
  216.   i2=i2*2.0;
  217.   }
  218. /*======================*/
  219. /* Calculate magnitudes */
  220. /*======================*/
  221. max=0.0;
  222. printf("Caluculating magnitudes...\n");  /* This calculates magnitude and */
  223. for(z=0;z<n/2;z++)                       /* finds the largest component   */
  224.   {                                      /* -max- for scaling purposes.   */
  225.   y=scram(z);
  226.   mag[y]=SQRT( POW( x1[y],2.0 ) + POW( x2[y],2.0  ) );
  227.   if (mag[y]>max) max=mag[y];
  228.   }
  229.                                       /* This shows the first 19 values */
  230. /*===========================*/       /* in the frequency spectrum.     */
  231. /*   Show tabular results    */       /* There are 128 all together.    */
  232. /*===========================*/
  233. printf("Harmonic      Real        Imaginary      Magnitude\n");
  234. for(z=0;z<20;z++)
  235.   {
  236.   y=scram(z);
  237.   printf("%3d %13.6e %13.6e %13.6e\n",(long)z,x1[y],x2[y],mag[y]);
  238.   }
  239.  
  240. /*============================*/
  241. /*  Show the graphic results  */
  242. /*============================*/
  243. printf("\nHarmonic      Magnitude\n");
  244. for(z=0;z<20;z++)
  245.   {
  246.   printf("%2ld ",(long)z);             /* This gives a poor man's graphic */
  247.   y=scram(z);                          /* display of the first 19 values  */
  248.   i=(long)((56.0*mag[y]/max)+0.5);     /* in the frequency spectrum.      */
  249.   for(j=1;j<i;j++)printf("=");         /* Plenty of room for improvement  */
  250.   printf("\n");                        /* here!                           */
  251.   }
  252.  
  253. #if IEEE_LATTICE
  254. /*============================*/
  255. /*  Close the Math Libraries  */
  256. /*============================*/
  257. CloseLibrary(MathIeeeDoubTransBase);
  258. CloseLibrary(MathIeeeDoubBasBase);
  259. #endif IEEE_LATTICE
  260.  
  261. }
  262.  
  263. /*=======================*/
  264. /*  Scramble/Unscramble  */
  265. /*=======================*/
  266. long scram(xx)
  267. double xx;                      /* Since the FFT algorithm results in */
  268. {                               /* a scrambled order for the freq.    */
  269. long yy=0;                      /* components, scram is used to fetch */
  270. long n1=N;                      /* the appropriate root of unity      */
  271. long w;                         /* at calculation time and at the end */
  272. for(w=1;w<(NPOW+1);w++)         /* to unscramble the results.         */
  273.  
  274.     {
  275.     n1=n1/2;
  276.     if(!( (long)xx < n1 ))
  277.       {
  278.       yy = yy + (long)POW( 2.0,(double)w-1 );
  279.       xx=xx-n1;
  280.       }
  281.     }
  282.  
  283. return(yy);
  284. }
  285.  
  286.  
  287.  
  288.