home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_100 / 158_01 / fwt.c < prev    next >
Text File  |  1987-10-12  |  4KB  |  163 lines
  1. /* 
  2. HEADER:     CUG
  3. TITLE:        FWT.C - Fast Walsh Transform
  4. VERSION:    1.00
  5. DATE:        05/27/85
  6. DESCRIPTION:    "A Fast Walsh Transform implementation based on
  7.         Cooley's successive-doubling method. See the
  8.         September '77 issue of BYTE for a description of
  9.         this alternative to the Fourier transform."
  10. KEYWORDS:    fwt, fft, filter, walsh, fourier, transform
  11. SYSTEM:        Any
  12. FILENAME:    FWT.C
  13. WARNINGS:    "Data must be presented in multiples of two."
  14. CRC:        xxxx
  15. SEE-ALSO:    FFT.C
  16. AUTHORS:    Ian Ashdown - byHeart Software
  17. COMPILERS:    Any C compiler
  18. REFERENCES:    AUTHORS: R.F. Gonzalez & P. Wintz;
  19.         TITLE:     Digital Image Processing;
  20.         AUTHORS: B. Jacobi;
  21.         TITLE:     Walsh Functions: A Digital Fourier
  22.              Series
  23.              pp. 190 - 198, BYTE, September 1977;
  24.         AUTHORS: J. Cooley, P. Lewis & P. Welch;{
  25.         TITLE:     The Fast Fourier Transform and Its
  26.              Applications
  27.              IEEE Trans. Educ.
  28.              pp. 27 - 34, Vol. E-12, No. 1, 1969;
  29. ENDREF
  30. */
  31.  
  32. /*-------------------------------------------------------------*/
  33.  
  34. #include "math.h"
  35.  
  36. typedef double REAL;
  37.  
  38. /* FWT() - Fast Walsh Transform accepts a one-dimensional
  39.  *       "double" array of one row of N+1 elements, where the
  40.  *       zeroth column element is not used and N is equal to
  41.  *       an integer power of two, and a variable n, set equal
  42.  *       to N above. The discrete Walsh transform is computed
  43.  *       in place and returned in the same array.
  44.  */
  45.  
  46. void fwt(n,point)
  47. int n;
  48. REAL point[];
  49. {
  50.   register int a,
  51.            b,
  52.            e;
  53.   int c = 1,
  54.       d;
  55.   REAL temp;
  56.   void r_bitrev();
  57.  
  58.   /* Reorder real data vector by bit reversal rule */
  59.  
  60.   r_bitrev(n,point);
  61.  
  62.   /* Perform successive doubling calculations */
  63.  
  64.   while(c < n)
  65.   {
  66.     d = c;
  67.     c *= 2;
  68.     for(b = 1; b <= d; b++)
  69.       for(a = b; a <= n ; a += c)
  70.       {
  71.     e = a + d;
  72.     temp = point[e];
  73.     point[e] = point[a] - temp;
  74.     point[a] = point[a] + temp;
  75.       }
  76.   }
  77.  
  78.   /* Normalize transformed data vector */
  79.  
  80.   for(a = 1; a <= n; a++)
  81.     point[a] = point[a]/n;
  82. }
  83.  
  84. /* IFWT() - Inverse Fast Walsh Transform accepts a one-dimensional
  85.  *        "double" array of one row of N+1 elements, where the
  86.  *        zeroth column element is not used and N is equal to
  87.  *        an integer power of two, and a variable n, set equal
  88.  *        to N above. The discrete Walsh transform is computed
  89.  *        in place and returned in the same array. (The
  90.  *        algorithm is identical to the Fast Walsh Transform,
  91.  *        except that the normalization of the transformed data
  92.  *        vector is not performed.)
  93.  */
  94.  
  95. void ifwt(n,point)
  96. int n;
  97. REAL point[];
  98. {
  99.   register int a,
  100.            b,
  101.            e;
  102.   int c = 1,
  103.       d;
  104.   REAL temp;
  105.   void r_bitrev();
  106.  
  107.   /* Reorder real data vector by bit reversal rule */
  108.  
  109.   r_bitrev(n,point);
  110.  
  111.   /* Perform successive doubling calculations */
  112.  
  113.   while(c < n)
  114.   {
  115.     d = c;
  116.     c *= 2;
  117.     for(b = 1; b <= d; b++)
  118.       for(a = b; a <= n ; a += c)
  119.       {
  120.     e = a + d;
  121.     temp = point[e];
  122.     point[e] = point[a] - temp;
  123.     point[a] = point[a] + temp;
  124.       }
  125.   }
  126. }
  127.  
  128. /* R_BITREV() - Reorder real data vector by bit reversal rule */
  129.  
  130. void r_bitrev(n,point)
  131. int n;
  132. REAL point[];
  133. {
  134.   register int i,
  135.            j,
  136.            k;
  137.   REAL temp;
  138.  
  139.   j = n/2 + 1;
  140.   for(i = 2; i < n; i++)
  141.   {
  142.     if(i < j)
  143.     {
  144.       temp = point[i];
  145.       point[i] = point[j];
  146.       point[j] = temp;
  147.     }
  148.     k = n/2;
  149.     while(k < j)
  150.     {
  151.       j = j - k;
  152.       k /= 2;
  153.     }
  154.     j = j + k;
  155.   }
  156. }
  157.  
  158. /* End of FWT.C */
  159. ansformed data
  160.  *        vector is not performed.)
  161.  */
  162.  
  163. void ifwt(n,point)