fft 修改 1.c

同样是浮点型的fft算法一样为vc开发的

//FFT.c  C callable RADIX2 (DIF) FFT function in C

#define PTS 256			    //# of points for FFT
typedef struct {float real,imag;} COMPLEX;
extern COMPLEX w[PTS];       	    //twiddle constants stored in w

void FFT(COMPLEX *Y, int N)      //input sample array, # of points     N=256&512,(8,16,32,64,128,1024)
{
 COMPLEX temp1,temp2;             //temporary storage variables          
 int i,j,k;                       //loop counter variables               
 int upper_leg, lower_leg;   	    //index of upper/lower butterfly leg   
 int leg_diff;                    //difference between upper/lower leg   
 int num_stages = 0;              //number of FFT stages (iterations)
 int index, step;                 //index/step through twiddle constant
 i = 1;                           //log(base2) of N points= # of stages  
 
 do
  {
   num_stages +=1;
   i = i*2;
  }while (i!=N);
 leg_diff = N>>2; 		     //difference between upper&lower legs “modifie
    step=1;   		            // 512/N ?? step between values in "twiddle.h "  “modified”??  
      
 for (i = 0;i < num_stages; i++)     //carve up "num_stages" for N-point FFT   one by one              
  {
   index = 0;
   for (j = 0; j < leg_diff; j++)
    {
     for (upper_leg = j; upper_leg < N; upper_leg += (2*leg_diff))
      {   
       lower_leg = upper_leg+leg_diff;
       temp1.real = (Y[upper_leg]).real + (Y[lower_leg]).real;
       temp1.imag = (Y[upper_leg]).imag + (Y[lower_leg]).imag;
       temp2.real = (Y[upper_leg]).real - (Y[lower_leg]).real;
       temp2.imag = (Y[upper_leg]).imag - (Y[lower_leg]).imag;
       (Y[lower_leg]).real = temp2.real*(w[index]).real
				    -temp2.imag*(w[index]).imag;
       (Y[lower_leg]).imag = temp2.real*(w[index]).imag
				    +temp2.imag*(w[index]).real;
       (Y[upper_leg]).real = temp1.real;
       (Y[upper_leg]).imag = temp1.imag;
      }
   index += step;
    }
     leg_diff = leg_diff/2;
     step *= 2;
      }
       
      
   
	   
	   //**************bit reversal for resequencing data************//
	   j = 0;
	   for (i = 1; i < (N-1); i++)    
           {
	     k = N/2;
	     while (k <= j)
	      {
	       j = j - k;
	       k = k/2;
	      }
	     j = j + k;
	     if (i<j)
	      {
	       temp1.real = (Y[j]).real;
	       temp1.imag = (Y[j]).imag;
	       (Y[j]).real = (Y[i]).real;
	       (Y[j]).imag = (Y[i]).imag;
	       (Y[i]).real = temp1.real;
	       (Y[i]).imag = temp1.imag;
	      } 
           }
                      return;
}
//***************************other bit_reverse*************************

/*		 j=1;
		 for ( i=1; i<N; ++i )       //npt change to N
		 {
			 if ( i<j )
			 {
				 tr=x[j].real;
				 ti=x[j].imag;
				 x[j].real=x[i].real;
				 x[j].imag=x[i].imag;
				 x[i].real=tr;
				 x[i].imag=ti;
				 k=N/2;
				 while( k<j )
				 {
					 j=j-k;
					 k=k/2;
				 }
			 }
			 else
			 {
				  k=N/2;
				  while( k<j )
                                 {
					  j=j-k;
					  k=k/2;
				  }
			 }
			 j=j+k;
	          }


   */
    /********************************************************************
	*********************************************************************
	*********************************************************************/
    // *************************the another way of fft***********
 
 /*   void fft ()
	{
		 int	sign;
		 long	m, irem, l, le, le1, k, ip, i, j;
		 double	ur, ui, wr, wi,tr,ti,temp;
		 extern	long npt;
		 extern	int inv;
		 extern complex x[size];

	
     /****************fft or ifft**************************/
	        
     /* 				   m=0;
					irem=npt;
					while ( irem>1 )
					{
						irem=irem/2;
						   m=m+1   ;
					}
					if ( inv==1 )
						sign=1;
					 else
						sign=-1;

	/******************fft for each stage*******************/
               
	/*				 for ( l=1; l<=m; l++ )
					 {
						le =pow ( 2,l );
						 le1=le/2      ;
						  ur=1.0       ;
						  ui=0         ;
						  wr=cos( pi/le1 );
						  wi=sign*sin(pi/le1);
						  for ( j=1; j<=le1;++j )
                          {
							  i=j;
							  while ( i<=npt )
							  {
								  ip=i+le1;
								  tr=x[ip].real*ur-x[ip].imag*ui;
								  ti=x[ip].imag*ur+x[ip].real*ui;
								  x[ip].real=x[i].real-tr;
								  x[ip].imag=x[i].imag-ti;
								  x[i].real=x[i].real+tr;
								  x[i].imag=x[i].imag+ti;
								  i=i+le;
							  }
							  temp=ur*wr-ui*wi;
							  ui=ui*wr-ur*wi  ;
							  ur=temp;

						  }
					 }

		/*******************if inverse fft is desired divide each coefficient by npt**********/
                /*				 if(inv==1)
					 {
						 for ( i=1; i<=npt; ++i)
						 {
							 x[i].real=x[i].real/npt;
							 x[i].imag=x[i],iamg/npt;
						 }
					 }
		   
         }