binaries/src/mafft/core/fft3.c

   1 #include "fft.h"
   2 #include "mtxutl.h"
   3
   4 /***********************************************************
   5         fft.c -- FFT (¹âÂ®FourierÊÑ´¹)
   6 ***********************************************************/
   7 /*
   8   ´Ø¿ô{\tt fft()}¤Î²¼ÀÁ¤±¤È¤·¤Æ»°³Ñ´Ø¿ôÉ½¤òºî¤ë.
   9 */
  10 static void make_sintbl(int n, float sintbl[])
  11 {
  12         int i, n2, n4, n8;
  13         double c, s, dc, ds, t;
  14
  15         n2 = n / 2;  n4 = n / 4;  n8 = n / 8;
  16 #if 0
  17         t = sin(PI / n);
  18         dc = 2 * t * t;  ds = sqrt(dc * (2 - dc));
  19         t = 2 * dc;  c = sintbl[n4] = 1;  s = sintbl[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 }
  37 /*
  38   ´Ø¿ô{\tt fft()}¤Î²¼ÀÁ¤±¤È¤·¤Æ¥Ó¥Ã¥ÈÈ¿Å¾É½¤òºî¤ë.
  39 */
  40 static void make_bitrev(int n, int bitrev[])
  41 {
  42         int i, j, k, n2;
  43
  44         n2 = n / 2;  i = j = 0;
  45         for ( ; ; ) {
  46                 bitrev[i] = j;
  47                 if (++i >= n) break;
  48                 k = n2;
  49                 while (k <= j) {  j -= k;  k /= 2;  }
  50                 j += k;
  51         }
  52 }
  53 /*
  54   ¹âÂ®FourierÊÑ´¹ (Cooley--Tukey¤Î¥¢¥ë¥´¥ê¥º¥à).
  55   É¸ËÜÅÀ¤Î¿ô {\tt n} ¤Ï2¤ÎÀ°¿ô¾è¤Ë¸Â¤ë.
  56   {\tt x[$k$]} ¤¬¼ÂÉô, {\tt y[$k$]} ¤¬µõÉô ($k = 0$, $1$, $2$,
  57   \ldots, $|{\tt n}| - 1$).
  58   ·ë²Ì¤Ï {\tt x[]}, {\tt y[]} ¤Ë¾å½ñ¤\ad¤µ¤ì¤ë.
  59   ${\tt n} = 0$ ¤Ê¤éÉ½¤Î¥á¥â¥ê¤ò²òÊü¤¹¤ë.
  60   ${\tt n} < 0$ ¤Ê¤éµÕÊÑ´¹¤ò¹Ô¤¦.
  61   Á°²ó¤È°Û¤Ê¤ë $|{\tt n}|$ ¤ÎÃÍ¤Ç¸Æ¤Ó½Ð¤¹¤È,
  62   »°³Ñ´Ø¿ô¤È¥Ó¥Ã¥ÈÈ¿Å¾¤ÎÉ½¤òºî¤ë¤¿¤á¤ËÂ¿¾¯Í¾Ê¬¤Ë»þ´Ö¤¬¤«¤«¤ë.
  63   ¤³¤ÎÉ½¤Î¤¿¤á¤Îµ\ad²±ÎÎ°è³ÍÆÀ¤Ë¼ºÇÔ¤¹¤ë¤È1¤òÊÖ¤¹ (Àµ¾ï½ªÎ»»þ
  64   ¤ÎÌá¤êÃÍ¤Ï0).
  65   ¤³¤ì¤é¤ÎÉ½¤Îµ\ad²±ÎÎ°è¤ò²òÊü¤¹¤ë¤Ë¤Ï ${\tt n} = 0$ ¤È¤·¤Æ
  66   ¸Æ¤Ó½Ð¤¹ (¤³¤Î¤È¤\ad¤Ï {\tt x[]}, {\tt y[]} ¤ÎÃÍ¤ÏÊÑ¤ï¤é¤Ê¤¤).
  67 */
  68 int fft_hontai(int n, float x[], float y[])
  69 {
  70         static int    last_n = 0;    /* Á°²ó¸Æ½Ð¤·»þ¤Î {\tt n} */
  71         static int   *bitrev = NULL; /* ¥Ó¥Ã¥ÈÈ¿Å¾É½ */
  72         static float *sintbl = NULL; /* »°³Ñ´Ø¿ôÉ½ */
  73         int i, j, k, ik, h, d, k2, n4, inverse;
  74         float t, s, c, dx, dy;
  75
  76         /* ½àÈ÷ */
  77         if (n < 0) {
  78                 n = -n;  inverse = 1;  /* µÕÊÑ´¹ */
  79         } else inverse = 0;
  80         n4 = n / 4;
  81         if (n != last_n || n == 0) {
  82                 last_n = n;
  83                 if (sintbl != NULL) free(sintbl);
  84                 if (bitrev != NULL) free(bitrev);
  85                 if (n == 0) return 0;  /* µ\ad²±ÎÎ°è¤ò²òÊü¤·¤¿ */
  86                 sintbl = malloc((n + n4) * sizeof(float));
  87                 bitrev = malloc(n * sizeof(int));
  88                 if (sintbl == NULL || bitrev == NULL) {
  89                         fprintf(stderr, "µ\ad²±ÎÎ°èÉÔÂ\ad\n");  return 1;
  90                 }
  91                 make_sintbl(n, sintbl);
  92                 make_bitrev(n, bitrev);
  93         }
  94         for (i = 0; i < n; i++) {    /* ¥Ó¥Ã¥ÈÈ¿Å¾ */
  95                 j = bitrev[i];
  96                 if (i < j) {
  97                         t = x[i];  x[i] = x[j];  x[j] = t;
  98                         t = y[i];  y[i] = y[j];  y[j] = t;
  99                 }
 100         }
 101         for (k = 1; k < n; k = k2) {    /* ÊÑ´¹ */
 102                 h = 0;  k2 = k + k;  d = n / k2;
 103                 for (j = 0; j < k; j++) {
 104                         c = sintbl[h + n4];
 105                         if (inverse) s = - sintbl[h];
 106                         else         s =   sintbl[h];
 107                         for (i = j; i < n; i += k2) {
 108                                 ik = i + k;
 109                                 dx = s * y[ik] + c * x[ik];
 110                                 dy = c * y[ik] - s * x[ik];
 111                                 x[ik] = x[i] - dx;  x[i] += dx;
 112                                 y[ik] = y[i] - dy;  y[i] += dy;
 113                         }
 114                         h += d;
 115                 }
 116         }
 117         if (! inverse)    /* µÕÊÑ´¹¤Ç¤Ê¤¤¤Ê¤én¤Ç³ä¤ë */
 118                 for (i = 0; i < n; i++) {  x[i] /= n;  y[i] /= n;  }
 119         return 0;  /* Àµ¾ï½ªÎ» */
 120 }
 121
 122 #if 0
 123 int fft( int n, Fukusosuu *in, int disp )
 124 {
 125         int i;
 126         static int last_n = 0;
 127         static float *x, *y;
 128         if( last_n != abs( n ) )
 129         {
 130
 131                 if( last_n )
 132                 {
 133                         free( x );
 134                         free( y );
 135                 }
 136                 last_n = abs( n );
 137                 x = calloc( n, sizeof( float ) );
 138                 y = calloc( n, sizeof( float ) );
 139         }
 140         for( i=0; i<n; i++ )
 141         {
 142                 x[i] = (float)in[i].R;
 143                 y[i] = (float)in[i].I;
 144         }
 145         fft_hontai( n, x, y );
 146         for( i=0; i<n; i++ )
 147         {
 148                 (in+i)->R = (double)x[i];
 149                 (in+i)->I = (double)y[i];
 150         }
 151 #if 1
 152 if( disp )
 153 {
 154     FILE *fp;
 155     fp = fopen( "outputOfFft", "w" );
 156     for( i=0; i<n; i++ )
 157         fprintf( fp, "%f %f\n", in[i].R, in[i].I );
 158     fclose( fp );
 159     system( "vi outputOfFft < /dev/tty > /dev/tty " );
 160 }
 161 #endif
 162 }
 163 #else
 164 int fft( int n, Fukusosuu *in, int disp )
 165 {
 166         int i, m;
 167         float *x, *y;
 168
 169         m = abs( n );
 170
 171         x = calloc( m, sizeof( float ) );
 172         y = calloc( m, sizeof( float ) );
 173
 174         for( i=0; i<m; i++ )
 175         {
 176                 x[i] = (float)in[i].R;
 177                 y[i] = (float)in[i].I;
 178         }
 179
 180 #if 0
 181 if( disp )
 182 {
 183     FILE *fp;
 184     fp = fopen( "inputOfFft", "w" );
 185         fprintf( fp, "m=%d\n", m );
 186     for( i=0; i<m; i++ )
 187         fprintf( fp, "%f %f\n", x[i], y[i] );
 188     fclose( fp );
 189     system( "vi inputOfFft < /dev/tty > /dev/tty " );
 190 }
 191 #endif
 192         fft_hontai( n, x, y );
 193 #if 0
 194 if( disp )
 195 {
 196     FILE *fp;
 197     fp = fopen( "outputOfFft", "w" );
 198         fprintf( fp, "m=%d\n", m );
 199     for( i=0; i<m; i++ )
 200         fprintf( fp, "%f %f\n", x[i], y[i] );
 201     fclose( fp );
 202     system( "vi outputOfFft < /dev/tty > /dev/tty " );
 203 }
 204 #endif
 205         for( i=0; i<m; i++ )
 206         {
 207                 (in+i)->R = (double)x[i];
 208                 (in+i)->I = (double)y[i];
 209         }
 210         free( x );
 211         free( y );
 212 }
 213 #endif