binaries/src/mafft/core/fft.c

   1 #include "mltaln.h"
   2 #include "mtxutl.h"
   3
   4 /*
   5   from "C gengo niyoru saishin algorithm jiten" ISBN4-87408-414-1 Haruhiko Okumura
   6 */
   7 static void make_sintbl(int n, float sintbl[])
   8 {
   9         int i, n2, n4, n8;
  10         double c, s, dc, ds, t;
  11
  12         n2 = n / 2;  n4 = n / 4;  n8 = n / 8;
  13         t = sin(PI / n);
  14         dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
  15         t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[0] = 0;
  16         for (i = 1; i < n8; i++) {
  17                 c -= dc;  dc += t * c;
  18                 s += ds;  ds -= t * s;
  19                 sintbl[i] = s;  sintbl[n4 - i] = c;
  20         }
  21         if (n8 != 0) sintbl[n8] = sqrt(0.5);
  22         for (i = 0; i < n4; i++)
  23                 sintbl[n2 - i] = sintbl[i];
  24         for (i = 0; i < n2 + n4; i++)
  25                 sintbl[i + n2] = - sintbl[i];
  26 }
  27 /*
  28   {\tt fft()}.
  29 */
  30 static void make_bitrev(int n, int bitrev[])
  31 {
  32         int i, j, k, n2;
  33
  34         n2 = n / 2;  i = j = 0;
  35         for ( ; ; ) {
  36                 bitrev[i] = j;
  37                 if (++i >= n) break;
  38                 k = n2;
  39                 while (k <= j) {  j -= k;  k /= 2;  }
  40                 j += k;
  41         }
  42 }
  43 /*
  44 */
  45 int fft(int n, Fukusosuu *x, int freeflag)
  46 {
  47         static TLS int    last_n = 0;    /*  {\tt n} */
  48         static TLS int   *bitrev = NULL; /*  */
  49         static TLS float *sintbl = NULL; /*  */
  50         int i, j, k, ik, h, d, k2, n4, inverse;
  51         float t, s, c, dR, dI;
  52
  53                 if (freeflag)
  54                 {
  55                         if (bitrev) free(bitrev);
  56                         if (sintbl) free(sintbl);
  57                         return( 0 );
  58                 }
  59
  60         /*  */
  61         if (n < 0) {
  62                 n = -n;  inverse = 1;  /*  */
  63         } else inverse = 0;
  64         n4 = n / 4;
  65         if (n != last_n || n == 0) {
  66                 last_n = n;
  67 #if 0
  68                 if (sintbl != NULL) {
  69                                         free(sintbl);
  70                                         sintbl = NULL;
  71                                 }
  72                 if (bitrev != NULL) {
  73                                         free(bitrev);
  74                                         bitrev = NULL;
  75                                 }
  76                 if (n == 0) return 0;  /*  */
  77                 sintbl = (float *)malloc((n + n4) * sizeof(float));
  78                 bitrev = (int *)malloc(n * sizeof(int));
  79 #else /* by T. Nishiyama */
  80                                 sintbl = realloc(sintbl, (n + n4) * sizeof(float));
  81                                 bitrev = realloc(bitrev, n * sizeof(int));
  82 #endif
  83                 if (sintbl == NULL || bitrev == NULL) {
  84                         fprintf(stderr, "\n");  return 1;
  85                 }
  86                 make_sintbl(n, sintbl);
  87                 make_bitrev(n, bitrev);
  88         }
  89         for (i = 0; i < n; i++) {    /*  */
  90                 j = bitrev[i];
  91                 if (i < j) {
  92                         t = x[i].R;  x[i].R = x[j].R;  x[j].R = t;
  93                         t = x[i].I;  x[i].I = x[j].I;  x[j].I = t;
  94                 }
  95         }
  96         for (k = 1; k < n; k = k2) {    /*  */
  97 #if 0
  98                                 fprintf( stderr, "%d / %d\n", k, n );
  99 #endif
 100                 h = 0;  k2 = k + k;  d = n / k2;
 101                 for (j = 0; j < k; j++) {
 102 #if 0
 103                                         if( j % 1 == 0 )
 104                                                 fprintf( stderr, "%d / %d\r", j, k );
 105 #endif
 106                         c = sintbl[h + n4];
 107                         if (inverse) s = - sintbl[h];
 108                         else         s =   sintbl[h];
 109                         for (i = j; i < n; i += k2) {
 110 #if 0
 111                                                                 if( k>=4194000 ) fprintf( stderr, "in loop %d -  %d < %d, k2=%d\r", j, i, n, k2 );
 112 #endif
 113                                 ik = i + k;
 114                                 dR = s * x[ik].I + c * x[ik].R;
 115                                 dI = c * x[ik].I - s * x[ik].R;
 116                                 x[ik].R = x[i].R - dR;  x[i].R += dR;
 117                                 x[ik].I = x[i].I - dI;  x[i].I += dI;
 118                         }
 119                         h += d;
 120                 }
 121         }
 122         if (! inverse)    /* n */
 123                 for (i = 0; i < n; i++) {  x[i].R /= n;  x[i].I /= n;  }
 124         return 0;  /*  */
 125 }