binaries/src/mafft/core/fft_NC.c

   1 #include "fft.h"
   2 #include "mtxutl.h"
   3
   4 #define SWAP( a, b ) tempr=( a ); ( a ) = ( b ); ( b ) = tempr;
   5
   6 void four1( float *data, int nn, int isign )
   7 {
   8         unsigned long n, mmax, m, j, istep, i;
   9         double wtemp, wr, wpr, wpi, wi, theta;
  10         float tempr, tempi;
  11
  12         n = nn << 1;
  13         j = 1;
  14         for( i=1; i<n; i+=2 )
  15         {
  16                 if( j > i )
  17                 {
  18                         SWAP( data[j], data[i] );
  19                         SWAP( data[j+1], data[i+1] );
  20                 }
  21                 m = n >> 1;
  22                 while( m >= 2 && j > m )
  23                 {
  24                         j -= m;
  25                         m >>= 1;
  26                 }
  27                 j += m;
  28         }
  29         mmax = 2;
  30         while( n > mmax )
  31         {
  32                 istep = mmax << 1;
  33                 theta = isign * ( 2 * PI / mmax );
  34                 wtemp = sin( 0.5 * theta );
  35                 wpr = -2.0 * wtemp * wtemp;
  36                 wpi = sin( theta );
  37                 wr = 1.0;
  38                 wi = 0.0;
  39                 for( m=1; m<mmax; m+=2 )
  40                 {
  41                         for( i=m; i<=n; i+=istep )
  42                         {
  43                                 j = i + mmax;
  44                                 tempr = wr * data[j] - wi * data[j+1];
  45                                 tempi = wr * data[j+1] + wi * data[j];
  46                                 data[j] = data[i] - tempr;
  47                                 data[j+1] = data[i+1] - tempi;
  48                                 data[i] += tempr;
  49                                 data[i+1] += tempi;
  50                         }
  51                         wr = ( wtemp = wr ) * wpr - wi * wpi + wr;
  52                         wi = wi * wpr + wtemp * wpi + wi;
  53                 }
  54                 mmax = istep;
  55         }
  56 }
  57
  58 int fft( int n, Fukusosuu *in, int disp )
  59 {
  60         int i, m, sign;
  61         float *x;
  62         m = abs( n );
  63         sign = n / m;
  64         x = calloc( m*2+1, sizeof( float ) );
  65
  66         for( i=0; i<m; i++ )
  67         {
  68                 x[i*2+1] = in[i].R;
  69                 x[i*2+2] = in[i].I;
  70         }
  71 #if 0
  72 if( disp )
  73 {
  74     FILE *fp;
  75     fp = fopen( "inputOfFft", "w" );
  76     for( i=0; i<m; i++ )
  77         fprintf( fp, "%f %f\n", x[i*2+1], x[i*2+2] );
  78     fclose( fp );
  79     system( "vi inputOfFft < /dev/tty > /dev/tty " );
  80 }
  81 #endif
  82         four1( x, m, sign );
  83 #if 0
  84 if( disp )
  85 {
  86     FILE *fp;
  87     fp = fopen( "outputOfFft", "w" );
  88     for( i=0; i<m; i++ )
  89         fprintf( fp, "%f %f\n", x[i*2+1], x[i*2+2] );
  90     fclose( fp );
  91     system( "vi outputOfFft < /dev/tty > /dev/tty " );
  92 }
  93 #endif
  94         for( i=0; i<m; i++ )
  95         {
  96                 in[i].R = x[i*2+1];
  97                 in[i].I = x[i*2+2];
  98         }
  99 }
 100
 101
 102 #if 0
 103
 104
 105 int fft( int n, Fukusosuu *in, int disp )
 106 {
 107         int i, m;
 108         float *x, *y;
 109
 110         m = abs( n );
 111
 112         x = calloc( 100000, sizeof( float ) );
 113         y = calloc( 100000, sizeof( float ) );
 114
 115         for( i=0; i<m; i++ )
 116         {
 117                 x[i] = (float)in[i].R;
 118                 y[i] = (float)in[i].I;
 119         }
 120
 121         fft_hontai( n, x, y );
 122         for( i=0; i<m; i++ )
 123         {
 124                 (in+i)->R = (double)x[i];
 125                 (in+i)->I = (double)y[i];
 126         }
 127         free( x );
 128         free( y );
 129 }
 130 #endif