home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / procssng / ccs / ccs-11.lha / ccs-lib / lib / fourier.c < prev    next >
Encoding:
C/C++ Source or Header  |  1991-12-04  |  5.1 KB  |  216 lines
  1. /*
  2.  * "fourier.c"
  3.  */
  4.  
  5. #include    "complex.h"
  6. #include    <math.h>
  7.  
  8. COMPLEX *W_factors = 0;        /* array of W-factors */
  9. unsigned Nfactors = 0;        /* number of entries in W-factors */
  10.  
  11. /*
  12.  * load_w puts Wn ^ k (= e ^ (2pi * i * k / n)) in W_factors [k], 0 <= k < n.
  13.  * If n is equal to Nfactors then nothing is done, so the same W_factors
  14.  * array can used for several transforms of the same number of samples.
  15.  * Notice the explicit calculation of sines and cosines, an iterative approach
  16.  * introduces substantial errors.
  17.  */
  18. int load_w(n)
  19. unsigned n;
  20. {
  21. #define pi_2    3.1415926535897932384626434 * 2
  22. register unsigned    k;
  23. register FType    scale = pi_2 / n;
  24.  
  25.     if (n == Nfactors)
  26.         return 0;
  27.     if (Nfactors != 0 && W_factors != 0)
  28.         free ((char *) W_factors);
  29.     if ((Nfactors=n) == 0)
  30.         return 0;
  31.     if ((W_factors = (COMPLEX *) malloc (n * sizeof (COMPLEX))) == 0)
  32.         return -1;
  33.     {
  34.     register FType    pos=0;
  35.         for (k=0; k < n; k++) {
  36.         c_re (W_factors [k]) = cos (pos);
  37.         c_im (W_factors [k]) = sin (pos);
  38.         pos += scale;
  39.         }
  40.     }
  41. return 0;
  42. }
  43.  
  44. /*
  45.  * Recursive (reverse) complex fast Fourier transform on the n
  46.  * complex samples of array in, with the Cooley-Tukey method.
  47.  * The result is placed in out.  The number of samples, n, is arbitrary.
  48.  * The algorithm costs O (n * (r1 + .. + rk)), where k is the number
  49.  * of factors in the prime-decomposition of n (also the maximum
  50.  * depth of the recursion), and ri is the i-th primefactor.
  51.  */
  52. Fourier (in, n, out)
  53. register COMPLEX    *in;
  54. register unsigned    n;
  55. register COMPLEX    *out;
  56. {
  57. register unsigned    r;
  58. unsigned radix ();
  59.  
  60.     if ((r = radix (n)) < n)
  61.         split (in, r, n / r, out);
  62.     join (in, n / r, n, out);
  63. }
  64.  
  65. /*
  66.  * Give smallest possible radix for n samples.
  67.  * Determines (in a rude way) the smallest primefactor of n.
  68.  */
  69. static unsigned radix (n)
  70. register unsigned    n;
  71. {
  72. register unsigned    r=2, q;
  73.  
  74.     if (n < r)    return 1;
  75.  
  76.     if (n & 1){
  77.         q = n << 1 + (n>16);
  78.         while (++r < q && n % r);
  79.         if (r == q)    r = n;
  80.     }
  81.     return r;
  82. }
  83.  
  84. /*
  85.  * Split array in of r * m samples in r parts of each m samples,
  86.  * such that in [i] goes to out [(i % r) * m + (i / r)].
  87.  * Then call for each part of out Fourier, so the r recursively
  88.  * transformed parts will go back to in.
  89.  */
  90. static split (in, r, m, out)
  91. COMPLEX *in;
  92. register unsigned r, m;
  93. COMPLEX *out;
  94. {
  95.     register unsigned k, s, i, j;
  96.  
  97.     for (k = j = 0; k < r; k++)
  98.         for (s = 0, i = k; s < m; s++, i += r, j++)
  99.             out [j] = in [i];
  100.  
  101.     for (k = 0; k < r; k++, out += m, in += m)
  102.         Fourier (out, m, in);
  103. }
  104.  
  105. /*
  106.  * Sum the n / m parts of each m samples of in to n samples in out.
  107.  *            r - 1
  108.  * Out [j] becomes  sum  in [j % m] * W (j * k).  Here in is the k-th
  109.  *            k = 0   k           n         k
  110.  * part of in (indices k * m ... (k + 1) * m - 1), and r is the radix.
  111.  * For k = 0, a complex multiplication with W (0) is avoided.
  112.  */
  113. static join (in, m, n, out)
  114. register COMPLEX    *in;
  115. register unsigned    m, n;
  116. register COMPLEX    *out;
  117. {
  118. register unsigned i, j, jk, s;
  119.  
  120.     for (s = 0; s < m; s++)
  121.         for (j = s; j < n; j += m) {
  122.         out [j] = in [s];
  123.         for (i = s + m, jk = j; i < n; i += m, jk += j)
  124.         { COMPLEX C1, C2; C1 = in[i]; C2 = W (n, jk);
  125.               c_re (out[j]) += C1.re * C2.re - C1.im * C2.im;
  126.               c_im (out[j]) += C1.re * C2.im + C1.im * C2.re;
  127.         }
  128.         }
  129. }
  130.  
  131. /*
  132.  * "Double Fourier.c"
  133.  */
  134.  
  135. DBCOMPLEX *DBW_factors = 0;        /* array of DBW-factors */
  136. unsigned DBNfactors = 0;        /* number of entries in DBW-factors */
  137.  
  138. /*
  139.  * load_w puts Wn ^ k (= e ^ (2pi * i * k / n)) in DBW_factors [k], 0 <= k < n.
  140.  * If n is equal to DBNfactors then nothing is done, so the same DBW_factors
  141.  * array can used for several transforms of the same number of samples.
  142.  * Notice the explicit calculation of sines and cosines, an iterative approach
  143.  * introduces substantial errors.
  144.  */
  145. int load_DBw(n)
  146. unsigned n;
  147. {
  148. #define pi_2    3.1415926535897932384626434 * 2
  149. register unsigned    k;
  150. register double    scale = pi_2 / n;
  151.  
  152.     if (n == DBNfactors)
  153.         return 0;
  154.     if (DBNfactors != 0 && DBW_factors != 0)
  155.         free ((char *) DBW_factors);
  156.     if ((DBNfactors = n) == 0)
  157.         return 0;
  158.     if ((DBW_factors = (DBCOMPLEX *) malloc (n * sizeof (DBCOMPLEX))) == 0)
  159.         return -1;
  160.     {
  161.     register double    pos=0;
  162.         for (k=0; k < n; k++) {
  163.         c_re (DBW_factors [k]) = cos (pos);
  164.         c_im (DBW_factors [k]) = sin (pos);
  165.         pos += scale;
  166.         }
  167.     }
  168. return 0;
  169. }
  170.  
  171. DBFourier (in, n, out)
  172. register DBCOMPLEX    *in;
  173. register unsigned    n;
  174. register DBCOMPLEX    *out;
  175. {
  176. register unsigned    r;
  177. unsigned radix ();
  178.  
  179.     if ((r = radix (n)) < n)
  180.         DBsplit (in, r, n / r, out);
  181.     DBjoin (in, n / r, n, out);
  182. }
  183.  
  184. static DBsplit (in, r, m, out)
  185. DBCOMPLEX *in;
  186. register unsigned r, m;
  187. DBCOMPLEX *out;
  188. {
  189. register unsigned k, s, i, j;
  190.  
  191.     for (k = j = 0; k < r; k++)
  192.         for (s = 0, i = k; s < m; s++, i += r, j++)
  193.             out [j] = in [i];
  194.  
  195.     for (k = 0; k < r; k++, out += m, in += m)
  196.         DBFourier (out, m, in);
  197. }
  198.  
  199. static DBjoin (in, m, n, out)
  200. register DBCOMPLEX    *in;
  201. register unsigned    m, n;
  202. register DBCOMPLEX    *out;
  203. {
  204. register unsigned i, j, jk, s;
  205.  
  206.     for (s = 0; s < m; s++)
  207.         for (j = s; j < n; j += m) {
  208.         out [j] = in [s];
  209.         for (i = s + m, jk = j; i < n; i += m, jk += j)
  210.         { DBCOMPLEX C1, C2; C1 = in[i]; C2 = DBW (n, jk);
  211.               c_re (out[j]) += C1.re * C2.re - C1.im * C2.im;
  212.               c_im (out[j]) += C1.re * C2.im + C1.im * C2.re;
  213.         }
  214.         }
  215. }
  216.