home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / zpp / fft.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1992-04-03  |  15.7 KB  |  482 lines
  1. /* -------------------------------------------------------------------- */
  2. /* fft.c        by Carl W. Moreland   Version 2.3     04/03/92          */
  3. /* -------------------------------------------------------------------- */
  4. /*                                                                      */
  5. /* Fast Fourier Transform                                               */
  6. /*                                                                      */
  7. /* The following program is intended for Fourier analysis of circuit    */
  8. /* simulations. It can be invoked from the command line with no inter-  */
  9. /* action so it is ideal for batch files. The syntax is:                */
  10. /*                                                                      */
  11. /* fft infile.ext <outfile.ext </c> <'Comments'>> </n=x> </e> </s>      */
  12. /*                                                                      */
  13. /* infile.ext is the name of the file containing data. There must be    */
  14. /*    2^n data points and must represent an integer number of cycles.   */
  15. /*                                                                      */
  16. /* outfile.ext is an optional file to store the harmonic information.   */
  17. /*                                                                      */
  18. /* /c allows the user to enter comments for each analysis. This is in   */
  19. /*    addition to the command line comment.  Each comment line is limi- */
  20. /*    ted to 80 characters; there is no line limit.                     */
  21. /*                                                                      */
  22. /* /n=x specifies the highest order harmonic to be printed. By default  */
  23. /*    the program will print up to the 10th, or numpts/2-1, whatever    */
  24. /*    is larger. If x is larger than numpts/2-1 it will be ignored.     */
  25. /*                                                                      */
  26. /* /e specifies exponential format instead of decimal. This is for      */
  27. /*    those working on 14+ bit projects.                                */
  28. /*                                                                      */
  29. /* /s suppresses screen output                                          */
  30. /*                                                                      */
  31. /* Version 2 now supports comments from the source file.  They must     */
  32. /* begin with the ""#"" character.  DC offset and fundamental amplitude */
  33. /* are now printed.                                                     */
  34. /*                                                                      */
  35. /* Version 2.1 adds dynamic memory allocation for the data array so max */
  36. /* number of points is now limited by RAM. (05/04/90)                   */
  37. /*                                                                      */
  38. /* Version 2.2 fixes a bug that gives incorrect results when the DC     */
  39. /* offset is greater than the fundamental amplitude. Harmonic bins are  */
  40. /* now sorted wrt the fundamental. (04/17/91)                           */
  41. /*                                                                      */
  42. /* Version 2.3 uses the C++ complex class. (04/03/92)                   */
  43. /*                                                                      */
  44. /* Sun compiler: cc fft.c -lm -offt (Version 2.2 or before)             */
  45. /* Note: Sun C does not recognize # conditionals or modern parameter    */
  46. /* passing.                                                             */
  47. /* -------------------------------------------------------------------- */
  48.  
  49. #include <math.h>
  50. #include <ctype.h>
  51. #include <malloc.h>
  52. #include <stdio.h>
  53. #include <stdlib.h>
  54. #include <string.h>
  55.  
  56. #include "complex.h"
  57.  
  58. void get_args(int argc, char *argv[]);
  59. int  get_data(void);
  60. void put_data(complex x[]);
  61. void fft(int N, complex *x);
  62.  
  63. struct arg_struct        /* Holds the command line arguments    */
  64. {
  65.   char input[20];
  66.   char output[20];
  67.   int  num_harm;
  68.   int  com_flag;
  69.   int  exp_flag;
  70.   int  sup_flag;
  71.   int  out_flag;
  72.   char comment[81];
  73. } args;
  74.  
  75. char *message[] =
  76. {
  77.   "+----------------------------------------------------------------------+",
  78.   "| Fast Fourier Transform   Version 2.3   04/03/92    Carl W. Moreland  |",
  79.   "|                                                                      |",
  80.   "| Syntax:                                                              |",
  81.   "|                                                                      |",
  82.   "| fft infile.ext <outfile.ext </c> <'Comments'>> </n=x> </e> </s>      |",
  83.   "|                                                                      |",
  84.   "| infile.ext is the name of the file containing data. There must be    |",
  85.   "|    2^n data points and must represent an integer number of cycles.   |",
  86.   "|                                                                      |",
  87.   "| outfile.ext is an optional file to store the harmonic information.   |",
  88.   "|                                                                      |",
  89.   "| /c allows interactive comments, 80 characters max.                   |",
  90.   "| /n=x specifies the highest order harmonic to be printed. By default  |",
  91.   "|      the program will print up to the 10th or numpts/2-1.            |",
  92.   "| /e specifies exponential format instead of decimal for the output.   |",
  93.   "| /s suppresses screen output                                          |",
  94.   "+----------------------------------------------------------------------+"
  95. };
  96.  
  97. complex *x;
  98.  
  99. /* ----- main ---------------------------------------------------------    */
  100.  
  101. main(int argc, char *argv[])
  102. {
  103.   register int i;
  104.   int numpts;
  105.  
  106.   if(argc < 2)            /* Check for proper syntax        */
  107.   {
  108.     for(i=0; i<18; i++)
  109.       printf("%s\n", message[i]);
  110.     exit(0);
  111.   }
  112.  
  113.   get_args(argc, argv);        /* Parse the command line arguments    */
  114.  
  115.   x = (complex *)malloc(8*sizeof(complex));
  116.   if(x == NULL)
  117.   {
  118.     printf("Insufficient memory for 8 data points\n");
  119.     exit(0);
  120.   }
  121.  
  122.   numpts = get_data();        /* Read the input file            */
  123.   if(numpts/2 < args.num_harm)    /* See if there's less than 10 harmonics*/
  124.     args.num_harm = numpts/2 - 1;
  125.  
  126.   fft(numpts, x);        /* Perform FFT on data, return Mag & Ph */
  127.   put_data(x);            /* Write the output to screen/file    */
  128.   free(x);
  129.  
  130.   return 0;
  131. }
  132.  
  133. /* ----- get_args -----------------------------------------------------    */
  134.  
  135. void get_args(int argc, char *argv[])
  136. {
  137.   register int i, j, k;
  138.   char temp[81];            /* Holds the comment string    */
  139.  
  140.   args.output[0]  = '\0';        /* Initialize all parameters    */
  141.   args.num_harm   = 10;
  142.   args.com_flag   = 0;
  143.   args.exp_flag   = 0;
  144.   args.sup_flag   = 0;
  145.   args.out_flag   = 0;
  146.   args.comment[0] = '\0';
  147.  
  148.   strcpy(args.input, argv[1]);        /* store the input filename    */
  149.  
  150.   for(i=2; i<argc; i++)            /* Parse thru each argument    */
  151.   {
  152.                     /* Check for "/" or "-" switch    */
  153.     if(argv[i][0] == '/' || argv[i][0] == '-')
  154.     {
  155.       if(argv[i][1] == 'c' || argv[i][1] == 'C') args.com_flag = 1;
  156.       if(argv[i][1] == 'e' || argv[i][1] == 'E') args.exp_flag = 1;
  157.       if(argv[i][1] == 's' || argv[i][1] == 'S') args.sup_flag = 1;
  158.       if(argv[i][1] == 'n' || argv[i][1] == 'N')
  159.       {
  160.     k = 0;
  161.     for(j=2;;j++)
  162.     {
  163.       if(argv[i][j] == '\0')
  164.       {
  165.         temp[k++] = argv[i][j];
  166.         break;
  167.       }
  168.       if(isdigit(argv[i][j])) temp[k++] = argv[i][j];
  169.     }
  170.     args.num_harm = atoi(temp);
  171.       }
  172.       continue;
  173.     }
  174.  
  175.     if(argv[i][0] == 0x27)        /* Check for "'" for comment    */
  176.     {
  177.       j = 1;
  178.       for(k=0; k<81; k++)        /* up to 80 characters        */
  179.       {
  180.     if(argv[i][j] == 0x27)        /* Check for "'" end of comment    */
  181.     {
  182.       temp[k] = '\0';
  183.       break;
  184.     }
  185.     if(argv[i][j] == '\0')        /* Check for new word        */
  186.     {
  187.       temp[k] = ' ';
  188.       i++;
  189.       j=0;
  190.       continue;
  191.     }
  192.     if(k == 80)            /* 80 is the max so terminate    */
  193.     {
  194.       temp[k] = '\0';
  195.       break;
  196.     }
  197.     temp[k] = argv[i][j++];        /* Copy character to temp string*/
  198.       }
  199.       strcpy(args.comment, temp);
  200.       continue;
  201.     }
  202.                     /* If it's nothing else then it */
  203.     strcpy(args.output, argv[i]);    /* must be the output file name    */
  204.     args.out_flag = 1;
  205.   }
  206. }
  207.  
  208. /* ----- get_data -----------------------------------------------------    */
  209.  
  210. int get_data()
  211. {
  212.   register int i;
  213.   int points = 8;                   /* start with 8 data points    */
  214.   char temp[81];
  215.   FILE *infile, *outfile;
  216.  
  217.   if((infile = fopen(args.input, "r")) == NULL)    /* Open input file    */
  218.   {
  219.     puts(" File not found\n");
  220.     free(x);
  221.     exit(0);
  222.   }
  223.  
  224.   if(args.out_flag) outfile = fopen(args.output, "a");
  225.  
  226.   i = 0;
  227.  
  228.   /* Read file until EOF    */
  229.  
  230.   while(fgets(temp, 81, infile) != NULL)
  231.   {
  232.     if(temp[0] == '#')
  233.     {
  234.       if(args.out_flag)
  235.       {
  236.         temp[0] = 0x20;
  237.         fprintf(outfile, "%s", temp);
  238.       }
  239.       continue;
  240.     }
  241.     if(i == points)
  242.     {
  243.       points *= 2;
  244.       x = (complex *)realloc(x, points*sizeof(complex));
  245.       if(x == NULL)
  246.       {
  247.         printf("Insufficient memory for %d data points\n", points);
  248.         exit(0);
  249.       }
  250.     }
  251.     x[i].re = atof(temp);
  252.     x[i++].im = 0;
  253.   }
  254.  
  255.   fclose(infile);
  256.   if(args.out_flag) fclose(outfile);
  257.   return(i);
  258. }
  259.  
  260. /* ----- put_data -----------------------------------------------------    */
  261.  
  262. void put_data(complex x[])
  263. {
  264.   register int i;
  265.   char comment[81];
  266.   double db, bits;
  267.   double sum = 0;
  268.   FILE *outfile;
  269.  
  270.   if(args.out_flag)            /* Write data to an output file    */
  271.   {
  272.     outfile = fopen(args.output, "a");    /* Open file for appending    */
  273.  
  274.     if(args.com_flag)            /* Allow interactive comments    */
  275.     {
  276.       printf("Enter comment(s), <RET> when through.\n");
  277.       for(;;)
  278.       {
  279.         gets(comment);
  280.         if(comment[0] == NULL) break;
  281.         fprintf(outfile, " %s\n", comment);
  282.       }
  283.     }
  284.   }
  285.  
  286.   if(!args.sup_flag)            /* Write data to the screen    */
  287.   {
  288.     printf(" DC offset is %.2f\n", x[0].re);
  289.     printf(" Fundamental amplitude is %.4f\n\n", x[1].re);
  290.     printf("  #    Magnitude    Phase        dB       Bits \n");
  291.     printf(" ---   ---------   --------   --------   ------\n");
  292.   }
  293.   if(args.out_flag)
  294.   {
  295.     fprintf(outfile, "\n");
  296.     fprintf(outfile, " DC offset is %.2f\n", x[0].re);
  297.     fprintf(outfile, " Fundamental amplitude is %.4f\n\n", x[1].re);
  298.     fprintf(outfile, "  #    Magnitude    Phase        dB       Bits \n");
  299.     fprintf(outfile, " ---   ---------   --------   --------   ------\n");
  300.   }
  301.  
  302.   for(i=2; i<=args.num_harm; i++)    /* Loop thru each harmonic    */
  303.   {
  304.     if(x[i].re == 0)
  305.     {
  306.       db  = -999;
  307.       bits = 999;
  308.     }
  309.     else
  310.     {
  311.       db   = 20*log10(x[i].re);        /* Convert to dB         */
  312.       bits = -db/6.0206;        /* Calculate number of bits    */
  313.     }
  314.  
  315.     if(args.exp_flag)            /* Check for /e switch        */
  316.     {
  317.       if(!args.sup_flag)
  318.         printf(" %2d    %.4e   %8.3f   %8.3f   %6.3f\n", i, x[i].re, x[i].im, db, bits);
  319.       if(args.out_flag)
  320.     fprintf(outfile, " %2d    %.4e   %8.3f   %8.3f   %6.3f\n", i, x[i].re, x[i].im, db, bits);
  321.     }
  322.     else
  323.     {
  324.       if(!args.sup_flag)
  325.         printf(" %2d    %9.7f   %8.3f   %8.3f   %6.3f\n", i, x[i].re, x[i].im, db, bits);
  326.       if(args.out_flag)
  327.     fprintf(outfile, " %2d    %9.7f   %8.3f   %8.3f   %6.3f\n", i, x[i].re, x[i].im, db, bits);
  328.     }
  329.  
  330.     sum += x[i].re*x[i].re;        /* Keep a running total        */
  331.   }
  332.  
  333.   sum  = sqrt(sum);            /* Calculate THD        */
  334.   if(sum == 0)
  335.   {
  336.     bits = 99;
  337.   }
  338.   else
  339.   {
  340.     db   = 20*log10(sum);        /* Convert to dB        */
  341.     bits = -db/6.0206;            /* Calculate number of bits    */
  342.   }
  343.  
  344.   if(!args.sup_flag)
  345.     printf(" ----------------------------------------------\n");
  346.   if(args.out_flag)
  347.     fprintf(outfile, " ----------------------------------------------\n");
  348.  
  349.   if(args.exp_flag)
  350.   {
  351.     if(!args.sup_flag)
  352.       printf(" THD = %.4e%% or %6.3f bits\n\n", 100*sum, bits);
  353.     if(args.out_flag)
  354.       fprintf(outfile, " THD = %.4e%% or %6.3f bits\n\n\n", 100*sum, bits);
  355.   }
  356.   else
  357.   {
  358.     if(!args.sup_flag)
  359.       printf(" THD = %7.4f%% or %6.3f bits\n\n", 100*sum, bits);
  360.     if(args.out_flag)
  361.       fprintf(outfile, " THD = %7.4f%% or %6.3f bits\n\n\n", 100*sum, bits);
  362.   }
  363.  
  364.   if(args.out_flag) fclose(outfile);
  365.  
  366.   return;
  367. }
  368.  
  369. /* --------------------------------------------------------------------    */
  370. /*                                                                      */
  371. /* >>>>>>>>>>>>>>>>> Decimation in Time FFT Algorithm <<<<<<<<<<<<<<<<< */
  372. /*                                                                      */
  373. /* The following subroutine is taken from Digital Signal Processing by    */
  374. /* A. Oppenheim and R. Schafer, p.331. The original FORTRAN code was     */
  375. /* translated into C by Carl Moreland. Since FORTRAN arrays are indexed */
  376. /* from 1 to N several modifications were made to accomodate C's 0 to   */
  377. /* N-1 index. Usage:                                                    */
  378. /*                                                                      */
  379. /* void fft(N, x)                                                       */
  380. /*                                                                      */
  381. /* integer N = # of samples                                             */
  382. /* complex x is a pointer to an array of complex of size N, where       */
  383. /*   complex is defined as a structure of double re, double im.         */
  384. /*                                                                      */
  385. /* --------------------------------------------------------------------    */
  386.  
  387. void fft(int N, complex *x)
  388. {
  389.   register int i, j, k;
  390.   int M, LE, LE1, ip, fund=1, num;
  391.   complex u, w, t, tm;
  392.   double PI = 3.14159265358;
  393.  
  394.   M = log((double)N)/log(2.0);    /* Convert numpts to log base 2        */
  395.   j = -N/2;            /* Initial value of j for bit reversal    */
  396.  
  397.   for(i=0; i<N; i++)
  398.   {
  399.  
  400.     k = N/2;            /* The next few lines calculate the     */
  401.                 /* proper j value for bit reversing    */
  402.     while(k<=j)            /* x[i] and x[j]. Some modifications    */
  403.     {                /* were made here to deal with C's       */
  404.       j = j - k;        /* array indices such as the initial      */
  405.       k = k/2;            /* j value above. A Polish way do it     */
  406.     }                /* but it works...            */
  407.     j = j + k;
  408.  
  409.  
  410.     if(i<j)            /* If the time is right swap x[i] and    */
  411.     {                /* x[j]. The complex number t is a tem- */
  412.       t = x[j];            /* porary holder.            */
  413.       x[j] = x[i];
  414.       x[i] = t;
  415.     }
  416.   }
  417.  
  418.   for(k=1; k<=M; k++)       /* Run through M stages of computations */
  419.   {
  420.     LE  = pow(2.0,(double)k);
  421.     LE1 = LE/2;
  422.     u   = Z1;
  423.     w   = complex(cos(PI/(double)LE1), -sin(PI/(double)LE1));
  424.  
  425.     for(j=0; j<LE1; j++)        /* Do N/2 butterfly computations    */
  426.     {
  427.  
  428.       for(i=j; i<N; i+=LE)      /* x[i] is one of the butterfly terms...*/
  429.       {
  430.     ip = i+LE1;        /* x[ip] is the other term...        */
  431.  
  432.     /* Butterfly algorithm calculates x[i]' and x[ip]' from x[i]     */
  433.     /* and x[ip] as follows:                    */
  434.     /*    x[i]'  = x[i] + x[ip]*w                    */
  435.     /*      x[ip]' = x[i] - x[ip]*w                    */
  436.     /* where w = exp(j2cr/LE) = W(LE)^r                 */
  437.  
  438.     t = x[ip] * u;        /* Temporary holder            */
  439.         x[ip] = x[i] - t;
  440.         x[i] += t;
  441.       }
  442.  
  443.       u *= w;
  444.     }
  445.   }
  446.  
  447.   Complex.SetArgMode(Z_DEGREES);
  448.  
  449.   x[0].topolar();
  450.  
  451.   for(i=1; i<N/2; i++)
  452.   {
  453.     x[i].topolar();             /* Convert to magnitude & phase        */
  454.  
  455.     if(x[i].re > x[fund].re)    /* Find the fundamental         */
  456.       fund = i;
  457.   }
  458.  
  459.   x[N/2] = x[0];        /* Sort the bins. Use the upper half of */
  460.                              /*   of the fft array to temporarily    */
  461.   for(i=1; i<N/2; i++)        /*   hold the harmonics.        */
  462.   {
  463.     num = fund*i;
  464.     while(num > N/2)
  465.       num = abs(N - num);
  466.     x[N/2+i] = x[num];
  467.   }
  468.  
  469.   x[0] = x[N/2];
  470.   x[1] = x[N/2+1];        /* Now copy the newly sorted harmonics  */
  471.                                /*   back down to the lower half of the */
  472.   for(i=2; i<N/2; i++)        /*   array and calculate the relative     */
  473.   {
  474.     x[i] = x[N/2+i];        /*   magnitudes.            */
  475.     x[i].re /= x[1].re;
  476.   }
  477.  
  478.   x[0].re /= N;
  479.   x[1].re /= N/2;
  480. }
  481.  
  482.