home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS 1992 December / simtel1292_SIMTEL_1292_Walnut_Creek.iso / msdos / pctech / hlsrc.arc / HLFLOAT.C < prev    next >
C/C++ Source or Header  |  1988-09-09  |  9KB  |  403 lines
  1. /*+
  2.     Name:       HLFLOAT.C
  3.     Author:     Kent J. Quirk
  4.             (c) Copyright 1988 Ziff Communications Co.
  5.     Date:       April 1988
  6.     Abstract:   Performs a floating point benchmark by generating a
  7.             waveform using Fourier synthesis, then running an FFT
  8.             over the waveform to generate Fourier Components, and
  9.             finally regenerating the waveform using the derived
  10.             parameters.
  11.     History:    09-Sep-88    kjq    Version 1.00
  12. -*/
  13.  
  14. #include <stdio.h>
  15. #include <stdlib.h>
  16. #include <conio.h>
  17. #include <float.h>
  18. #include <math.h>
  19. #include <graph.h>
  20.  
  21. #include "hl.h"
  22. #include "hltimes.h"
  23.  
  24. extern int *colorlist;
  25.  
  26. /* Macros for complex number manipulation */
  27. #define C_MULT(a,b) (a.x * (b.x + b.y) + a.y * (b.x + b.y))
  28. #define C_ADD(r,a,b) r.x = a.x + b.x; r.y = a.y + b.y
  29. #define C_SUB(r,a,b) r.x = a.x - b.x; r.y = a.y - b.y
  30. typedef struct complex COMPLEX;
  31.  
  32. COMPLEX *data;
  33. #define PI2 6.283185
  34.  
  35. /***** b i t r e v ****
  36.     Given an integer and a number of bits, this reverses the bits
  37.     in the integer.  For example, the integer 5 (0101) reverses to
  38.     the integer 10 (1010) with 4-bit numbers, but reverses to 160
  39.     (10100000) with 8-bit numbers.
  40.     This is useful in the FFT algorithm, since the values are processed
  41.     by the "butterfly" algorithm.
  42. ***********************/
  43. unsigned int bitrev(in, nbits)
  44. unsigned int in;
  45. int nbits;
  46. {
  47.     int out = 0;
  48.     while (nbits--)         /* for each bit on the input */
  49.     {
  50.         out <<= 1;          /* make a place for an output bit */
  51.         if (in & 1)         /* if the low bit in input is set */
  52.             out |= 1;       /* set the output bit */
  53.         in >>= 1;           /* and rotate the input */
  54.     }
  55.     return(out);
  56. }
  57.  
  58. /**** f f t ****
  59.     Abstract:    Performs an n-point complex FFT
  60.     Parameters: COMPLEX array[] -- a pointer to an array of COMPLEX 
  61.         numbers containing the data to be transformed.
  62.         int N -- the number of points to be processed
  63.         int NU -- the number of bits in N (ex: N=1024 -> NU=10)
  64.     Returns:    Nothing, but array[] is changed to contain the FFT data.
  65.     Comments:    NU could easily be calculated from N.
  66. ****************************/
  67. void fft(array, N, NU)
  68. COMPLEX array[];
  69. int N;
  70. int NU;
  71. {
  72.     int n2, i, l, k, NU1;
  73.     double p, c, s;
  74.     COMPLEX t;
  75.  
  76.     n2 = N/2;
  77.     NU1 = NU - 1;
  78.     k = 0;
  79.     for (l=0; l<NU; l++)
  80.     {
  81.         while (k<N)
  82.         {
  83.             for (i=0; i<n2; i++)
  84.             {
  85.                 p = (double)bitrev(k/(1 << NU1), NU);
  86.                 p = PI2 * p/(double)N;
  87.                 c = cos(p);
  88.                 s = sin(p);
  89.                 t.x = array[k+n2].x * c + array[k+n2].y * s;
  90.                 t.y = array[k+n2].y * c - array[k+n2].x * s;
  91.                 C_SUB(array[k+n2], array[k], t);
  92.                 C_ADD(array[k], array[k], t);
  93.                 k++;
  94.             }
  95.             k += n2;
  96.         }
  97.         k = 0;
  98.         NU1--;
  99.         n2 /= 2;
  100.     }
  101.  
  102.     for (k=0; k<N; k++)
  103.     {
  104.         i = bitrev(k, NU);
  105.         if (i > k)
  106.         {
  107.             t = array[k];
  108.             array[k] = array[i];
  109.             array[i] = t;
  110.         }
  111.     }
  112. }
  113.  
  114. /**** l o g 2 ****
  115.     Abstract:    Given a number n, this calculates the integer log base 2 of n
  116.     Parameters: An integer n.
  117.     Returns:    The integer log base 2 of n -- log2(1024) == 10
  118. ****************************/
  119. int log2(n)
  120. unsigned int n;
  121. {
  122.     int i = 0;
  123.     for (i = 0; n > 1; i++)
  124.         n >>= 1;
  125.     return(i);
  126. }
  127.  
  128.  
  129. /**** u s a g e ****
  130.     Abstract:    Prints usage info and exits
  131.     Parameters: None
  132.     Returns:    Never 
  133. ****************************/
  134. void usage()
  135. {
  136.     printf("Usage:  HLFLOAT [-nNBITS] [-q|-t|-m|-s] [-r] [-f] [-g] [-?]\n");
  137.     printf("   where NBITS is the number of bits in the FFT\n");
  138.     printf("   (example: 7 is 128 points, 9 is 512 points)\n");
  139.     printf("   Other options choose data set: -q is squarewave,\n");
  140.     printf("   -t is triangle, -s sinewave, -m mixture of two sines\n");
  141.     printf("   -r toggles addition of random factor to the data\n");
  142.     printf("   -h toggles hi-freq filter; filters everything above the\n");
  143.     printf("      tenth harmonic before reconstructing the inverse FFT.\n");
  144.     printf("   -b sets batch mode; no wait for user keypress.\n");
  145.     printf("   Defaults: -n11, Squarewave waveform, -f and -r ON, -b OFF");
  146.     printf("   -g prevents use of graphics mode.\n");
  147.     printf("   Note:  This could take a long time on a machine without");
  148.     printf("     an 80x87 math coprocessor.");
  149.     exit(1);
  150. }
  151.  
  152. #define SQUAREWAVE 0
  153. #define TRIANGLEWAVE 1
  154. #define SINEWAVE 2
  155. #define LFSQUAREWAVE 3
  156.  
  157.  
  158. TIME_REC timerec[] = {
  159.     { 0L, "Total: Floating Point" },
  160.     { 0L, "Forward FFT, 2048 points" },
  161.     { 0L, "Reverse FFT, 2048 points" },
  162.     { 0L, "Total Time including graphics" },
  163.     { -1L, "HLFLOAT" }
  164. };
  165.  
  166. int main(argc, argv)
  167. int argc;
  168. char *argv[];
  169. {
  170.     int nbits = 11;          /* default to 2048-point run */
  171.     int npts = 1 << nbits;
  172.     int i;
  173.     double d;
  174.     double xsc, ysc;
  175.     long ffttime;
  176.     struct videoconfig config;
  177.     int dataset = SQUAREWAVE;
  178.     int random = 1, hf_filter = 1, batch = 0, bench = 0;
  179.     int program = -1;
  180.     int graphics = 1;
  181.     char *filename = NULL;
  182.  
  183.     for (i=1; i<argc; i++)
  184.     {
  185.         if (argv[i][0] != '-')
  186.             usage();
  187.  
  188.         switch(tolower(argv[i][1])) {
  189.         case 'n':
  190.             nbits = atoi(argv[i]+2);
  191.             npts = 1 << nbits;
  192.             break;
  193.         case 'q':
  194.             dataset = SQUAREWAVE;
  195.             break;
  196.         case 's':
  197.             dataset = SINEWAVE;
  198.             break;
  199.         case 't':
  200.             dataset = TRIANGLEWAVE;
  201.             break;
  202.         case 'm':
  203.             dataset = LFSQUAREWAVE;
  204.             break;
  205.         case 'r':
  206.             random = !random;
  207.             break;
  208.         case 'h':
  209.             hf_filter = !hf_filter;
  210.             break;
  211.         case 'a':
  212.             batch = 1;
  213.             break;
  214.         case 'b':
  215.             bench = 1;
  216.             break;
  217.         case 'g':
  218.             graphics = 0;
  219.             break;
  220.         case 'f':
  221.             filename = argv[i]+2;
  222.             break;
  223.         case 'p':
  224.             program = atoi(argv[i]+2);
  225.             break;
  226.         default:
  227.         case '?':
  228.             usage();
  229.             break;
  230.         }
  231.     }
  232.  
  233.     if (argc == 1)
  234.         random = hf_filter = 1;         /* if no args, use specials */
  235.  
  236.     if ((data = malloc(sizeof(COMPLEX) * (npts+1))) == NULL)
  237.     {
  238.         printf("Unable to allocate space for %d points!\n", npts);
  239.         exit(1);
  240.     }
  241.  
  242.     d = PI2/npts;
  243.     if (graphics)
  244.     {
  245.         if (!init_video((struct videoconfig far *)&config))
  246.         {
  247.             printf("Unable to use graphics.\n");
  248.             graphics = 0;
  249.         }
  250.     }
  251.     else
  252.         printf("Not using graphics.\n");
  253.  
  254.     if (graphics)
  255.     {
  256.         ysc = (double)config.numypixels / -3.0;         /* invert the axes */
  257.         xsc = (double)config.numxpixels / (double)npts;
  258.         if (xsc < 2.0)
  259.         {
  260.             i = 1 << log2(config.numxpixels);
  261.             xsc = (double)i / (double)npts;
  262.         }
  263.  
  264.         _setcolor(colorlist[15]);          /* force white */
  265.         _moveto(0,0);
  266.         _lineto(0,config.numypixels);
  267.         _setlogorg(0, config.numypixels/2);
  268.         _moveto(0,0);
  269.         _lineto(config.numxpixels, 0);
  270.  
  271.         _settextcolor(colorlist[15]);
  272.         _outtext("Data...");
  273.     }
  274.     else
  275.         printf("Data...");
  276.  
  277.     start_timer(0);
  278.     for (i=0; i<npts; i++)
  279.     {
  280.         switch (dataset) {
  281.         case SQUAREWAVE:
  282.             data[i].x = (i > npts/2) ? -1.0 : 1.0;
  283.             break;
  284.         case SINEWAVE:
  285.             data[i].x = sin((double)i * d);
  286.             break;
  287.         case TRIANGLEWAVE:
  288.             data[i].x = 1.0 - 2.0 * modf(2*(double)i/(double)npts, &d);
  289.             break;
  290.         case LFSQUAREWAVE:
  291.             data[i].x = sin((double)i * d) + sin(11.0*i * d)/3.0;
  292.             break;
  293.         }
  294.         data[i].y = 0.0;
  295.         if (random)
  296.             data[i].x += (double)rand()/200000.0;
  297.  
  298.     }
  299.  
  300.     if (graphics)
  301.     {
  302.         _setcolor(colorlist[2]);
  303.         _moveto(0, (int)(ysc*data[0].x));
  304.         for (d=1.0; d<(float)npts; d += 1.0)
  305.             _lineto((int)(xsc*d), (int)(ysc*data[(int)d].x));
  306.  
  307.         _settextcolor(colorlist[2]);
  308.         _outtext("FFT...");
  309.     }
  310.     else
  311.         printf("FFT...");
  312.  
  313.     start_timer(1);
  314.     fft(data, npts, nbits);      /* convert it */
  315.     stop_timer();
  316.     ffttime = timerec[1].ticks = get_timer(1);
  317.  
  318.     if (graphics)
  319.     {
  320.         if (hf_filter)
  321.         {
  322.             _settextcolor(colorlist[3]);
  323.             _outtext("filtered...");
  324.         }
  325.  
  326.         for (i=0; i<=npts/2; i++)
  327.         {
  328.             d = cabs(data[i]);
  329.             if (hf_filter && ((i > 10) || (i == 0)))
  330.             {
  331.                 _setcolor(colorlist[2]);
  332.                 data[i].x = data[i].y = 0.0;      /* HF, DC Filter */
  333.                 data[npts-i].x = data[npts-i].y = 0.0;    /* on both ends */
  334.             }
  335.             else
  336.                 _setcolor(colorlist[3]);
  337.  
  338.             _moveto((int)(xsc*i), (int)(2.0*ysc*d/npts));  /* make bigger */
  339.             _lineto((int)(xsc*i), 0);
  340.             _moveto((int)(xsc*(npts-i)), (int)(2.0*ysc*d/npts));
  341.             _lineto((int)(xsc*(npts-i)), 0);
  342.         }
  343.  
  344.         _settextcolor(colorlist[4]);
  345.         _outtext("Reverse...");
  346.     }
  347.     else
  348.         printf("Reverse...");
  349.  
  350.     start_timer(1);
  351.     fft(data, npts, nbits);      /* and convert it back */
  352.     stop_timer();
  353.     ffttime += timerec[2].ticks = get_timer(1);
  354.  
  355.     if (graphics)
  356.     {
  357.         _setcolor(colorlist[4]);
  358.         _moveto(0, (int)(ysc*data[0].x/npts));
  359.         for (i=1; i<npts; i++)
  360.             _lineto((int)(xsc*i), (int)(ysc*data[i].x/-npts));
  361.  
  362.         stop_timer();
  363.         _settextcolor(colorlist[15]);
  364.         _outtext("Done.\nTest complete -- press any key...");
  365.         if (!batch)
  366.             getch();
  367.         _setvideomode(_DEFAULTMODE);
  368.     }
  369.     else
  370.     {
  371.         stop_timer();
  372.         printf("Done.\n");
  373.     }
  374.  
  375.     timerec[3].ticks = get_timer(0);
  376.     timerec[0].ticks = timerec[1].ticks + timerec[2].ticks;
  377.  
  378.     sprintf(timerec[1].desc, "Forward FFT, %d Points", npts);
  379.     sprintf(timerec[2].desc, "Reverse FFT, %d Points", npts);
  380.     
  381.     if ((!bench) && (!batch))
  382.     {
  383.         printf("FFT: %d points\n", npts);
  384.         for (i=0; timerec[i].ticks != -1; i++)
  385.             printf("%s secs:  %s\n", time_secs(timerec[i].ticks),
  386.                 timerec[i].desc);
  387.     }
  388.  
  389.     if ((filename != NULL) && (program != -1))
  390.     {
  391.         opentime(filename);
  392.         for (i=0; ; i++)
  393.         {
  394.             savetime(program, i, &timerec[i]);
  395.             if (timerec[i].ticks == -1)
  396.                 break;
  397.         }
  398.         closetime();
  399.     }
  400.  
  401.     return(0);
  402. }
  403.