home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 December / simtel1292_SIMTEL_1292_Walnut_Creek.iso / msdos / turbo_c / fft.arc / FFT.C next >
C/C++ Source or Header  |  1988-08-07  |  4KB  |  229 lines
  1. /*
  2.  *    fft.c
  3.  *
  4.  *    C Version 1.0 by Steve Sampson, Public Domain
  5.  *
  6.  *    This program is based on the work by W. D. Stanley
  7.  *    and S. J. Peterson, Old Dominion University.
  8.  *
  9.  *    This program produces a Frequency Domain display
  10.  *    from the Time Domain data input using the Fast Fourier Transform.
  11.  *
  12.  *    The REAL data is generated by the in-phase (I) channel and the
  13.  *    IMAGINARY data is produced by the quadrature-phase (Q) channel of
  14.  *    a Doppler Radar receiver.  The middle filter is zero Hz.  Closing
  15.  *    targets are displayed to the right, and Opening targets to the left.
  16.  *
  17.  *    Note: With IMAGINARY data set to zero the output is a mirror image.
  18.  *
  19.  *    Usage:    fft  samples  input_data  output_data
  20.  *    Where 'samples' is a power of two
  21.  *
  22.  *    Array Version for Turbo C 1.5
  23.  */
  24.  
  25. /* Includes */
  26.  
  27. #include <stdlib.h>
  28. #include <stdio.h>
  29. #include <math.h>
  30. #include <alloc.h>
  31.  
  32. /* Defines */
  33.  
  34. #define    TWO_PI    ((double)2.0 * M_PI)
  35.  
  36. /* Globals */
  37.  
  38. int    samples, power;
  39. double    *real, *imag, max;
  40. FILE    *fpi, *fpo;
  41.  
  42. /* Prototypes and forward declarations */
  43.  
  44. void    fft(void), max_amp(void), display(void);
  45. int    permute(int);
  46. double    magnitude(int);
  47.  
  48. /* The program */
  49.  
  50. main(argc, argv)
  51. int    argc;
  52. char    *argv[];
  53. {
  54.     int    n;
  55.  
  56.     if (argc != 4)  {
  57. err1:        fprintf(stderr, "Usage: fft samples input output\n");
  58.         fprintf(stderr, "Where samples is a power of 2\n");
  59.         exit(1);
  60.     }
  61.  
  62.     samples = abs(atoi(argv[1]));
  63.     power =  log10((double)samples) / log10((double)2.0);
  64.  
  65.     if ((real = (double *)malloc(samples * sizeof(double))) == NULL)
  66. err2:        fprintf(stderr, "Out of memory\n");
  67.  
  68.     if ((imag = (double *)malloc(samples * sizeof(double))) == NULL)
  69.         goto err2;
  70.  
  71.     if ((fpo = fopen(argv[3], "wb")) == (FILE *)NULL)
  72.         goto err1;
  73.  
  74.     if ((fpi = fopen(argv[2], "rb")) != (FILE *)NULL)  {
  75.         fread(real, sizeof(double), samples, fpi);
  76.         fread(imag, sizeof(double), samples, fpi);
  77.         fclose(fpi);
  78.     }
  79.     else
  80.         goto err1;
  81.  
  82.     fft();
  83.     max_amp();
  84.     display();
  85.  
  86.     fwrite(real, sizeof(double), samples, fpo);
  87.     fwrite(imag, sizeof(double), samples, fpo);
  88.  
  89.     fclose(fpo);
  90. }
  91.  
  92.  
  93. void fft()
  94. {
  95.     unsigned i1, i2, i3, i4, y;
  96.     int     loop, loop1, loop2;
  97.     double     a1, a2, b1, b2, z1, z2, v;
  98.  
  99.     /* Scale the data */
  100.  
  101.     for (loop = 0; loop < samples; loop++)  {
  102.         real[loop] /= (double)samples;
  103.         imag[loop] /= (double)samples;
  104.     }
  105.  
  106.     i1 = samples >> 1;
  107.     i2 = 1;
  108.     v = TWO_PI * ((double)1.0 / (double)samples);
  109.  
  110.     for (loop = 0; loop < power; loop++)  {
  111.         i3 = 0;
  112.         i4 = i1;
  113.  
  114.         for (loop1 = 0; loop1 < i2; loop1++)  {
  115.             y = permute(i3 / i1);
  116.             z1 =  cos(v * y);
  117.             z2 = -sin(v * y);
  118.  
  119.             for (loop2 = i3; loop2 < i4; loop2++)  {
  120.                 a1 = real[loop2];
  121.                 a2 = imag[loop2];
  122.  
  123.                 b1 = z1*real[loop2+i1] - z2*imag[loop2+i1];
  124.                 b2 = z2*real[loop2+i1] + z1*imag[loop2+i1];
  125.  
  126.                 real[loop2]      = a1 + b1;
  127.                 imag[loop2]      = a2 + b2;
  128.  
  129.                 real[loop2 + i1] = a1 - b1;
  130.                 imag[loop2 + i1] = a2 - b2;
  131.             }
  132.  
  133.             i3 += (i1 << 1);
  134.             i4 += (i1 << 1);
  135.         }
  136.  
  137.         i1 >>= 1;
  138.         i2 <<= 1;
  139.     }
  140. }
  141.  
  142. /*
  143.  *    Find maximum amplitude
  144.  */
  145.  
  146. void max_amp()
  147. {
  148.     double    mag;
  149.     int    loop;
  150.  
  151.     max = (double)0.0;
  152.     for (loop = 0; loop < samples; loop++)  {
  153.         if ((mag = magnitude(loop)) > max)
  154.             max = mag;
  155.     }
  156. }
  157.  
  158. /*
  159.  *    Display the frequency domain.
  160.  *    The filters are aranged so that DC is in the middle filter.
  161.  *    Thus -Doppler is on the left, +Doppler on the right.
  162.  */
  163.  
  164. void display()
  165. {
  166.     int    c, n, x, loop;
  167.  
  168.     n = samples / 2;
  169.     
  170.     for (loop = n; loop < samples; loop++)  {
  171.         x = (int)(magnitude(loop) * (double)56.0 / max);
  172.         printf("%d\t|", loop - n);
  173.         c = 0;
  174.         while (++c <= x)
  175.             putchar('=');
  176.  
  177.         putchar('\n');
  178.     }
  179.  
  180.     for (loop = 0; loop < n; loop++)  {
  181.         x = (int)(magnitude(loop) * (double)56.0 / max);
  182.         printf("%d\t|", loop + n);
  183.         c = 0;
  184.         while (++c <= x)
  185.             putchar('=');
  186.  
  187.         putchar('\n');
  188.     }
  189. }
  190.  
  191. /*
  192.  *    Calculate Power Magnitude
  193.  */
  194.  
  195. double magnitude(n)
  196. int    n;
  197. {
  198.     n = permute(n);
  199.     return (sqrt(real[n] * real[n] + imag[n] * imag[n]));
  200. }
  201.  
  202. /*
  203.  *    Bit reverse the number
  204.  *
  205.  *    Change 11100000b to 00000111b or vice-versa
  206.  */
  207.  
  208. int permute(index)
  209. int    index;
  210. {
  211.     int    n1, result, loop;
  212.  
  213.     n1 = samples;
  214.     result = 0;
  215.  
  216.     for (loop = 0; loop < power; loop++)  {
  217.         n1 >>= 1;            /* n1 / 2.0 */
  218.         if (index < n1)
  219.             continue;
  220.  
  221.         result += (int) pow((double)2.0, (double)loop);
  222.         index -= n1;
  223.     }
  224.  
  225.     return result;
  226. }
  227.  
  228. /* EOF */
  229.