home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 545b.lha / SRFFT / fftg.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-09-07  |  8.4 KB  |  302 lines
  1. /* ------------------------------ fftg.c ------------------------------------ */
  2. /*                                                                            */
  3. /* Author:      Eyal Lebedinsky                                               */
  4. /* Date:        May 1990                                                      */
  5. /* Version:     9 June 1991                                                   */
  6. /*                                                                            */
  7. /* Program generator for integer DFFT. This program was developed for a very  */
  8. /* specific need and no attempt was made to make it general. Still, the basic */
  9. /* algorithm is a standard one from the literature. For maximum performance   */
  10. /* of the integer calculation a special mechanism was employed to track shift */
  11. /* operations that could be delayed. Coefficients are kept as 16 bits with 1  */
  12. /* sign bit, 1 integer bit and 14 fractional bits. So a typicall multiply     */
  13. /* will do 16bit*16bits=32bits, then shift down one bit and use the HIGH 16   */
  14. /* bits as the result. Before adding two numbers they must be scaled down one */
  15. /* bit to avoid overflow. Adding more numbers calls for more scaling. The     */
  16. /* delayed scaling allows the intermediate array to hold numbers that are     */
  17. /* scaled to a varying degree until a clash occurs. Well, this space is too   */
  18. /* short to elaborate any further and if you did similar kind of work then    */
  19. /* you should know what I am talking about.                                   */
  20. /*                                                                            */
  21. /* usage: fftgc File Size EpName                                              */
  22. /*      File - output file name.                                              */
  23. /*      Size - sample order: 8 for 256 points, 10 for 1024...                 */
  24. /*      EpName - name of the generated function ('_' added)                   */
  25. /* Now assemble the output file and then use as:                              */
  26. /*                                                                            */
  27. /* #define M 8                                                                */
  28. /* #define N (1 << M)                                                         */
  29. /* short x[N+1];                [input/output array]                          */
  30. /* short qf[(N/2)+1];           [power spectrum]                              */
  31. /* extern void far fft ()                                                     */
  32. /*                                                                            */
  33. /* [then call it as: fft();]                                                  */
  34. /*                                                                            */
  35. /* This is a modified version of Sorensen et al, IEEE ASSP-35 no 6.           */
  36. /*                                                                            */
  37. /* - uses 3+3 complex mult.                                                   */
  38. /* - no pre scrambling. descrambling done when power spectrum calculated.     */
  39. /* - moved to 0-origin as proper for c.                                       */
  40. /* - trig values in table.                                                    */
  41. /* - the case for sin==cos treated separately.                                */
  42. /* - uses delayed scaling down.                                               */
  43. /* - does length 256 fft with 642 mults.                                      */
  44. /*                                                                            */
  45. /*                                                                            */
  46. /* This program is released into the public domain.                           */
  47. /*                                                                            */
  48. /*----------------------------------------------------------------------------*/
  49.  
  50. #include <stdio.h>
  51. #include <math.h>
  52.  
  53. #define SC15 32768                /* scaling factor */
  54.  
  55. static short mm[1025];                /* fp(16, 0) */
  56. static short ad[1025];                /* fp(16, 0) */
  57. static int cntx[9];
  58.  
  59. static short x[1025];                /* fp(16,0)  */
  60. static short qf[513];                /* fp(16,0)  */
  61.  
  62. static char *file_name, *ep_name;
  63.  
  64. extern void start_fft (), end_fft ();
  65. extern void fft1 (), fft2 (), fft3 (), fft4 (), fft5 (), fft7 (), fft8 ();
  66.  
  67. static int mad;                /* maximum adjust at this level */
  68.  
  69. static void
  70. smad (i, j)                    /* set proper adjust           */
  71. int    i, j;
  72. {
  73.     while (ad[i] < (mad-j)) {
  74.         fprintf (stderr, "x[%u] >>= %d (%d-%d)\n", i, 1, mad, j);
  75.         fft1 (mm[i]);
  76.         ++ad[i];
  77.         cntx[0] += 1;
  78.         cntx[1] += 1;
  79.         cntx[2] += 1;
  80.         cntx[5] += 1;
  81.     }
  82. }
  83.  
  84. main (argc, argv)
  85. int    argc;
  86. char    *argv[];
  87. {
  88.     int n, m, j, i, k, is, id, n1, n2, n4, n8, ind,
  89.         i0, i1, i2, i3, i4, i5, i6, i7, i8,
  90.         cc1, ss1, cc3, ss3;
  91.     float e, a, a3;
  92.  
  93.     if (argc < 4) {
  94.         printf ("usage: fftg File PowerOf2 EpName [stderrFile]\n");
  95.         exit (0);
  96.     }
  97.  
  98.     file_name = argv[1];
  99.     sscanf (argv[2], "%u", &m);
  100.     ep_name = argv[3];
  101.     if (argc > 4)
  102.         freopen (argv[4], "wt", stderr);
  103.     start_fft (file_name, ep_name);
  104.  
  105. /* ---------- Initialise  ------------------------------------------------ */
  106.  
  107.     n = 1 << m;
  108.  
  109.     for (i = 0; i < 9; ++i)
  110.         cntx[i] = 0;
  111.  
  112. /* ---------- Digit reverse counter -------------------------------------- */
  113.  
  114.     for (i = 1; i <= n; ++i) {           /* '-1' for shifting to 0-origin    */
  115.         mm[i] = i-1;
  116.         ad[i] = 0;
  117.     }
  118.  
  119.     j = 1;
  120.     n1 = n - 1;
  121.     for (i = 1; i <= n1; ++i) {
  122.         if (i < j) {
  123.             mm[i] = j-1;               /* '-1' for shifting to 0-origin    */
  124.             mm[j] = i-1;               /* '-1' for shifting to 0-origin    */
  125.         }
  126.  
  127.         k = n / 2;
  128.  
  129.         while (k < j) {
  130.             j = j - k;
  131.             k = k / 2;
  132.         }
  133.  
  134.         j = j + k;
  135.     }
  136.  
  137. /* ---------- Length two butterflies ------------------------------------- */
  138.  
  139.     mad = 0;
  140.  
  141.     is = 1;
  142.     id = 4;
  143.  
  144.     do {
  145.         for (i0 = is; i0 <= n; i0 += id) {
  146.             i1 = i0 + 1;
  147.  
  148.             smad (i0, 0);
  149.             smad (i1, 0);
  150.             fft2 (mm[i0], mm[i1]);
  151.             ++ad[i0];
  152.             ++ad[i1];
  153.             cntx[0] += 2;
  154.             cntx[1] += 2;
  155.             cntx[2] += 2;
  156.             cntx[3] += 2;
  157.         }
  158.         is = 2 * id - 1;
  159.         id = 4 * id;
  160.     } while (is < n);
  161.  
  162. /* ---------- L shaped butterflies --------------------------------------- */
  163.  
  164.     n2 = 2;
  165.     for (k = 2; k <= m; ++k) {
  166.         ++mad;
  167.         n2 = n2 * 2;
  168.         n4 = n2 / 4;
  169.         n8 = n2 / 8;
  170.         e = 6.2831853071719586 / n2;
  171.         is = 0;
  172.         id = n2 * 2;
  173.  
  174.         do {
  175.             cc1 = (SC15/2) * sqrt (0.5);
  176.             for (i = is; i <= n - 1; i += id) {
  177.                 i1 = i + 1;
  178.                 i2 = i1 + n4;
  179.                 i3 = i2 + n4;
  180.                 i4 = i3 + n4;
  181.  
  182.                 smad (i1, 0);
  183.                 smad (i3, 1);
  184.                 smad (i4, 1);
  185.                 fft3 (mm[i1], mm[i3], mm[i4]);
  186.                 ad[i1] += 1;
  187.                 ad[i3] += 2;
  188.                 ad[i4] += 2;
  189.                 cntx[0] += 3;
  190.                 cntx[1] += 3;
  191.                 cntx[2] += 5;
  192.                 cntx[3] += 4;
  193.  
  194.                 if (n4 != 1) {
  195.                     i1 = i1 + n8;
  196.                     i2 = i2 + n8;
  197.                     i3 = i3 + n8;
  198.                     i4 = i4 + n8;
  199.  
  200.                     smad (i1, 1);
  201.                     smad (i2, 0);
  202.                     smad (i3, 1);
  203.                     smad (i4, 1);
  204.                     fft4 (mm[i1], mm[i2], mm[i3], mm[i4], cc1);
  205.                     ad[i1] += 2;
  206.                     ad[i2] += 1;
  207.                     ad[i3] += 2;
  208.                     ad[i4] += 2;
  209.                     cntx[0] += 4;
  210.                     cntx[1] += 4;
  211.                     cntx[2] += 2;
  212.                     cntx[3] += 6;
  213.                     cntx[4] += 2;
  214.                 }
  215.             }
  216.             is = 2 * id - n2;
  217.             id = 4 * id;
  218.         } while (is < n);
  219.  
  220.         a = e;
  221.         for (j = 2; j <= n8; ++j) {
  222.             a3 = 3 * a;
  223.             cc1 = (SC15/2) * cos (a);
  224.             ss1 = (SC15/2) * sin (a);
  225.             cc3 = (SC15/2) * cos (a3);
  226.             ss3 = (SC15/2) * sin (a3);
  227. /*            a = j * e;*/
  228.             is = 0;
  229.             id = 2 * n2;
  230.  
  231.             do {
  232.                 for (i = is; i <= n - 1; i += id) {
  233.                     i1 = i + j;
  234.                     i2 = i1 + n4;
  235.                     i3 = i2 + n4;
  236.                     i4 = i3 + n4;
  237.                     i5 = i  + n4 - j + 2;
  238.                     i6 = i5 + n4;
  239.                     i7 = i6 + n4;
  240.                     i8 = i7 + n4;
  241.  
  242.                     smad (i1, 0);
  243.                     smad (i2, 0);
  244.                     smad (i3, 2);
  245.                     smad (i4, 2);
  246.                     smad (i5, 0);
  247.                     smad (i6, 0);
  248.                     smad (i7, 1);
  249.                     smad (i8, 1);
  250.                     if (ind = (ad[i3] < (mad-1))) {
  251.                         ++ad[i3];
  252.                         ++ad[i4];
  253.                         cntx[2] += 1;
  254.                     }
  255.                     fft5 (mm[i1], mm[i2], mm[i3], mm[i4],
  256.                         mm[i5], mm[i6], mm[i7], mm[i8],
  257.                         ss1-cc1, -ss1-cc1, cc1,
  258.                         ss3-cc3, -ss3-cc3, cc3, ind);
  259.                     ad[i1] += 1;
  260.                     ad[i2] += 1;
  261.                     ad[i3] += 2;
  262.                     ad[i4] += 2;
  263.                     ad[i5] += 1;
  264.                     ad[i6] += 1;
  265.                     ad[i7] += 2;
  266.                     ad[i8] += 2;
  267.                     cntx[0] += 8;
  268.                     cntx[1] += 8;
  269.                     cntx[2] += 4;
  270.                     cntx[3] += 18;
  271.                     cntx[4] += 6;
  272.                 }
  273.                 is = 2 * id - n2;
  274.                 id = 4 * id;
  275.             } while (is < n);
  276.             a += e;
  277.         }
  278.     }
  279.  
  280.     ++mad;
  281.     for (i = 1; i <= n; ++i)
  282.         smad (i, 0);
  283.  
  284.     fft7 (0, mm[1]);
  285.  
  286.     for (i = 1; i < n/2; ++i)
  287.         fft8 (i, mm[i+1], mm[n-i+1]);
  288.  
  289.     fft7 ((n/2), mm[n/2+1]);
  290.  
  291.     end_fft (file_name, ep_name);
  292.  
  293.     for (i = 0; i < 9; ++i)
  294.         fprintf (stderr, " %5u", i);
  295.     fprintf (stderr, "\n  Load Store Shift   Add  Mult   Adj\n");
  296.     for (i = 0; i < 9; ++i)
  297.         fprintf (stderr, " %5u", cntx[i]);
  298.     fprintf (stderr, "\n");
  299.  
  300.     exit (0);
  301. }
  302.