Logo Search packages:      
Sourcecode: mafft version File versions  Download package

fft_NC.c

#include "fft.h"
#include "mtxutl.h"

#define SWAP( a, b ) tempr=( a ); ( a ) = ( b ); ( b ) = tempr;

void four1( float *data, int nn, int isign )
{
      unsigned long n, mmax, m, j, istep, i;
      double wtemp, wr, wpr, wpi, wi, theta;
      float tempr, tempi;

      n = nn << 1;
      j = 1;
      for( i=1; i<n; i+=2 )
      {
            if( j > i )
            {
                  SWAP( data[j], data[i] );
                  SWAP( data[j+1], data[i+1] );
            }
            m = n >> 1;
            while( m >= 2 && j > m )
            {
                  j -= m;
                  m >>= 1;
            }
            j += m;
      }
      mmax = 2;
      while( n > mmax )
      {
            istep = mmax << 1;
            theta = isign * ( 2 * PI / mmax );
            wtemp = sin( 0.5 * theta );
            wpr = -2.0 * wtemp * wtemp;
            wpi = sin( theta );
            wr = 1.0;
            wi = 0.0;
            for( m=1; m<mmax; m+=2 )
            {
                  for( i=m; i<=n; i+=istep )
                  {
                        j = i + mmax;
                        tempr = wr * data[j] - wi * data[j+1];
                        tempi = wr * data[j+1] + wi * data[j];
                        data[j] = data[i] - tempr;
                        data[j+1] = data[i+1] - tempi;
                        data[i] += tempr;
                        data[i+1] += tempi;
                  }
                  wr = ( wtemp = wr ) * wpr - wi * wpi + wr;
                  wi = wi * wpr + wtemp * wpi + wi;
            }
            mmax = istep;
      }
}

int fft( int n, Fukusosuu *in, int disp )
{
      int i, m, sign;
      float *x;
      m = abs( n );
      sign = n / m;
      x = calloc( m*2+1, sizeof( float ) );

      for( i=0; i<m; i++ )
      {
            x[i*2+1] = in[i].R;
            x[i*2+2] = in[i].I;
      }
#if 0
if( disp )                                                      
{                                                                    
    FILE *fp;        
    fp = fopen( "inputOfFft", "w" );                                           
    for( i=0; i<m; i++ )                
        fprintf( fp, "%f %f\n", x[i*2+1], x[i*2+2] );        
    fclose( fp );        
    system( "vi inputOfFft < /dev/tty > /dev/tty " );        
}                
#endif        
      four1( x, m, sign );
#if 0
if( disp )                                                      
{                                                                    
    FILE *fp;        
    fp = fopen( "outputOfFft", "w" );                                           
    for( i=0; i<m; i++ )                
        fprintf( fp, "%f %f\n", x[i*2+1], x[i*2+2] );        
    fclose( fp );        
    system( "vi outputOfFft < /dev/tty > /dev/tty " );        
}                
#endif        
      for( i=0; i<m; i++ )
      {
            in[i].R = x[i*2+1];
            in[i].I = x[i*2+2];
      }
}
            

#if 0


int fft( int n, Fukusosuu *in, int disp )
{
      int i, m;
      float *x, *y;

      m = abs( n );

      x = calloc( 100000, sizeof( float ) );
      y = calloc( 100000, sizeof( float ) );

      for( i=0; i<m; i++ )
      {
            x[i] = (float)in[i].R;
            y[i] = (float)in[i].I;
      }     

      fft_hontai( n, x, y );
      for( i=0; i<m; i++ )
      {
            (in+i)->R = (double)x[i];
            (in+i)->I = (double)y[i];
      }     
      free( x );
      free( y );
}
#endif

Generated by  Doxygen 1.6.0   Back to index