binaries/src/mafft/core/fft2.c

   1 #include "fft.h"
   2 #include "mtxutl.h"
   3
   4 FILE *gfp;
   5 /*
   6   {\tt fft()}.
   7 */
   8 static void make_sintbl(int n, double *sintbl)
   9 {
  10         int i, n2, n4, n8;
  11         double c, s, dc, ds, t;
  12
  13         n2 = n / 2;  n4 = n / 4;  n8 = n / 8;
  14 #if 0
  15         t = sin(PI / n);
  16         dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
  17         t = 2 * dc; c = 1.0; s = 0.0;
  18         sintbl[n4] = 1.0;
  19         sintbl[0] = 0.0;
  20         for (i = 1; i < n8; i++) {
  21                 c -= dc;  dc += t * c;
  22                 s += ds;  ds -= t * s;
  23                 sintbl[i] = s;  sintbl[n4 - i] = c;
  24         }
  25         if (n8 != 0) sintbl[n8] = sqrt(0.5);
  26 #else
  27         sintbl[n4] = 1.0;
  28         sintbl[0] = 0.0;
  29         for (i = 1; i < n4; i++)
  30                         sintbl[i] = sin( (double)i / n * ( 2 * PI ) );
  31 #endif
  32         for (i = 0; i < n4; i++)
  33                 sintbl[n2 - i] = sintbl[i];
  34         for (i = 0; i < n2 + n4; i++)
  35                 sintbl[i + n2] = - sintbl[i];
  36 #if 0
  37 {
  38     FILE *fp;
  39     fp = fopen( "sintbl", "w" );
  40     for( i=0; i<n; i++ )
  41         fprintf( fp, "%d %f\n", i, sintbl[i] );
  42     fclose( fp );
  43     fp = fopen( "plot", "w" );
  44     fprintf( fp, "plot 'sintbl'\npause 5" );
  45     fclose( fp );
  46     system( "gnuplot plot" );
  47 }
  48 #endif
  49 }
  50 /*
  51   {\tt fft()}.
  52 */
  53 static void make_bitrev(int n, int *bitrev)
  54 {
  55         int i, j, k, n2;
  56
  57         n2 = n / 2;  i = j = 0;
  58         for ( ; ; ) {
  59                 bitrev[i] = j;
  60                 if (++i >= n) break;
  61                 k = n2;
  62                 while (k <= j) {  j -= k;  k /= 2;  }
  63                 j += k;
  64         }
  65 }
  66 /*
  67 */
  68 #if 0
  69 int fft(int n, Fukusosuu *x)
  70 {
  71         static int    last_n = 0;    /*  {\tt n} */
  72         static int   *bitrev = NULL; /*  */
  73         static float *sintbl = NULL; /*  */
  74         int i, j, k, ik, h, d, k2, n4, inverse;
  75         float t, s, c, dR, dI;
  76
  77         /*  */
  78         if (n < 0) {
  79                 n = -n;  inverse = 1;  /*  */
  80         } else inverse = 0;
  81         n4 = n / 4;
  82         if (n != last_n || n == 0) {
  83                 last_n = n;
  84                 if (sintbl != NULL) free((void *)sintbl);
  85                 if (bitrev != NULL) free((void *)bitrev);
  86                 if (n == 0) return 0;  /*  */
  87                 sintbl = malloc((n + n4) * sizeof(float));
  88                 bitrev = malloc(n * sizeof(int));
  89                 if (sintbl == NULL || bitrev == NULL) {
  90                         fprintf(stderr, "\n");  return 1;
  91                 }
  92                 make_sintbl(n, sintbl);
  93                 make_bitrev(n, bitrev);
  94         }
  95         for (i = 0; i < n; i++) {    /*  */
  96                 j = bitrev[i];
  97                 if (i < j) {
  98                         t = x[i].R;  x[i].R = x[j].R;  x[j].R = t;
  99                         t = x[i].I;  x[i].I = x[j].I;  x[j].I = t;
 100                 }
 101         }
 102         for (k = 1; k < n; k = k2) {    /*  */
 103                 h = 0;  k2 = k + k;  d = n / k2;
 104                 for (j = 0; j < k; j++) {
 105                         c = sintbl[h + n4];
 106                         if (inverse) s = - sintbl[h];
 107                         else         s =   sintbl[h];
 108                         for (i = j; i < n; i += k2) {
 109                                 ik = i + k;
 110                                 dR = s * x[ik].I + c * x[ik].R;
 111                                 dI = c * x[ik].I - s * x[ik].R;
 112                                 x[ik].R = x[i].R - dR;  x[i].R += dR;
 113                                 x[ik].I = x[i].I - dI;  x[i].I += dI;
 114                         }
 115                         h += d;
 116                 }
 117         }
 118         if (! inverse)    /* n */
 119                 for (i = 0; i < n; i++) {  x[i].R /= (double)n;  x[i].I /= (double)n;  }
 120         return 0;  /*  */
 121 }
 122 #else
 123 int fft(int n, Fukusosuu *x, int disp)
 124 {
 125         static int    last_n = 0;    /*  {\tt n} */
 126         static int    *bitrev = NULL; /*  */
 127         static double *sintbl = NULL; /*  */
 128         int i, j, k, ik, h, d, k2, n4, inverse;
 129         double t, s, c, dR, dI;
 130                 double tmp1, tmp2;
 131
 132         /*  */
 133         if (n < 0) {
 134                 n = -n;  inverse = 1;  /*  */
 135         } else inverse = 0;
 136 #if 0
 137 if( disp )
 138 {
 139     FILE *fp;
 140     fp = fopen( "inputOfFft", "w" );
 141     for( i=0; i<n; i++ )
 142         fprintf( fp, "%f %f\n", x[i].R, x[i].I );
 143     fclose( fp );
 144     system( "vi inputOfFft < /dev/tty > /dev/tty " );
 145 }
 146 #endif
 147         n4 = n / 4;
 148         if (n != last_n || n == 0) {
 149                 last_n = n;
 150                 if (sintbl != NULL) free((void *)sintbl);
 151                 if (bitrev != NULL) free((void *)bitrev);
 152                 if (n == 0) return 0;  /*  */
 153                 sintbl = (double *)malloc((n + n4) * sizeof(double));
 154                 bitrev = (int    *)malloc(n * sizeof(int));
 155                 if (sintbl == NULL || bitrev == NULL) {
 156                         fprintf(stderr, "\n");  return 1;
 157                 }
 158                 make_sintbl(n, sintbl);
 159                 make_bitrev(n, bitrev);
 160         }
 161         for (i = 0; i < n; i++) {    /*  */
 162                 j = bitrev[i];
 163                 if (i < j) {
 164                         t = x[i].R;  x[i].R = x[j].R;  x[j].R = t;
 165                         t = x[i].I;  x[i].I = x[j].I;  x[j].I = t;
 166                 }
 167         }
 168 #if 0
 169 if( disp ) gfp = fopen( "x", "w" );
 170 if( disp ) for( i=0; i<n; i++ )
 171         fprintf( gfp, "i=%d x=%f+i%f\n", i, x[i].R, x[i].I );
 172 #endif
 173         for (k = 1; k < n; k = k2) {    /*  */
 174                 h = 0;  k2 = k + k;  d = n / k2;
 175                 for (j = 0; j < k; j++) {
 176                         c = sintbl[h + n4];
 177                         if (inverse) s = - sintbl[h];
 178                         else         s =   sintbl[h];
 179                         for (i = j; i < n; i += k2) {
 180                                 ik = i + k;
 181                                                                 tmp1 = c * x[ik].I; tmp2 = s * x[ik].R;
 182                                                                 dI = tmp1 - tmp2;
 183                                                                 tmp1 = s * x[ik].I; tmp2 = c * x[ik].R;
 184                                                                 dR = tmp1 + tmp2;
 185 #if 0
 186 if( disp ) fprintf( gfp, "B k,j,i,ik=%d,%d,%d,s=%f,c=%f,x[%d]=%f+-i%f, dR=%f ? %f+%f \n", k, j, i, s, c, i, x[i].R, x[i].I, dR, s * x[ik].I, c * x[ik].R );
 187 if( disp ) fprintf( gfp, "B k,j,i,ik=%d,%d,%d,s=%f,c=%f,x[%d]=%f+-i%f, dR=%f\n", k, j, i, s, c, ik, x[ik].R, x[ik].I, dR );
 188 #endif
 189                                 x[ik].R = x[i].R - dR;  x[i].R += dR;
 190                                 x[ik].I = x[i].I - dI;  x[i].I += dI;
 191 #if 0
 192 if( disp ) fprintf( gfp, "A k,j,i=%d,%d,%d,s=%f,c=%f,x[%d]=%f+-i%f, dR=%f\n", k, j, i, s, c, ik, x[ik].R, x[ik].I, dR );
 193 #endif
 194                         }
 195                         h += d;
 196                 }
 197 #if 0
 198 if( disp ) fprintf( gfp, "ik=%d x=%f+i%f\n", ik, x[ik].R, x[ik].I );
 199 #endif
 200         }
 201         if (! inverse)    /* n */
 202                 for (i = 0; i < n; i++) {  x[i].R /= (double)n;  x[i].I /= (double)n;  }
 203 #if 0
 204 if( disp ) fclose( gfp );
 205 if( disp ) system( "vi x < /dev/tty > /dev/tty " );
 206 #endif
 207 #if 0
 208 if( disp )
 209 {
 210         FILE *fp;
 211     fp = fopen( "outputOfFft", "w" );
 212     for( i=0; i<n; i++ )
 213         fprintf( fp, "%f %f\n", x[i].R, x[i].I );
 214     fclose( fp );
 215     system( "vi outputOfFft < /dev/tty > /dev/tty " );
 216 }
 217 #endif
 218         return 0;  /*  */
 219 }
 220 #endif