all_radix_fft.c

此程式具有三种FFT演算法,分别是radix-2 radix-4 radix-8,可选择其中一种演算法,接着输入FFT点数,程式就会执行FFT/IFFT,可以藉此知道是否正确

#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>
#define  size 2048
#define  pi   3.14159265
#define  FFT  1
#define  IFFT 0
int fft_pt,bit_len;
double xr[size],xi[size];
double sr[size],si[size];
double T_xr[size],T_xi[size];
void get_data();
void show_data();
void but_cell(int fir_index,int sec_index);
int bit_reversal(int tr_index,int bit_length);
void do_bit_reversal();
void do_R2_FFT(int IsFFT,int R8_stage,int R4_stage,int R2_stage);
void do_R4_FFT(int IsFFT,int R8_stage,int R4_stage,int R2_stage);
void do_R8_FFT(int IsFFT,int R8_stage,int R4_stage,int R2_stage);
int main()
{
    int i;
    int FFT_type,R8_level,R4_level,R2_level;
    float tmp;
    double error=0;
    while(1)
    {
      printf("Choose FFT algorithm-(1)Radix-2 (2)Radix-4 (3)Radix-8 :");// choose FFT algorithm
      scanf("%d",&FFT_type);
      if(FFT_type>0 && FFT_type<4)
      break;
      printf("Input number is NOT between 1 and 3 !!\nPlease input again,Press any key to continue~\n");
      getch();
    }
    while(1)
    {
      printf("Input FFT point:");                            // input count of FFT point
      scanf("%d",&fft_pt);
      tmp=log2(fft_pt);
      if(tmp-floor(tmp)==0)
      break;
      printf("FFT point is NOT the power of 2 !!\nPlease input again,Press any key to continue~\n");
      getch();
    }
    bit_len=tmp;  
    switch(FFT_type)                                       // count of do Radix-2 or Radix-4 or Radix-8 FFT
    {
      case 1 :R8_level=0;
              R4_level=0;
              R2_level=bit_len;
              break;
      case 2 :R8_level=0;
              R4_level=bit_len/2;
              R2_level=bit_len%2;
              break;
      case 3 :R8_level=bit_len/3;
              R4_level=(bit_len%3)/2;                                    
              R2_level=(bit_len%3)%2;
              break;
    }
    get_data();                                            // get random data
    switch(FFT_type)                                       // Radix-4 or Radix-4 or Radix-8 FFT
    {
      case 1 :do_R2_FFT(FFT,R8_level,R4_level,R2_level);
              break;
      case 2 :do_R4_FFT(FFT,R8_level,R4_level,R2_level);
              break;
      case 3 :do_R8_FFT(FFT,R8_level,R4_level,R2_level);
              break;
    }    
    do_bit_reversal();                                     // bit reversal
  //  printf("\n\n  after FFT  \n");                         // show FFT result 
   // show_data();                                           
    for(i=0; i<fft_pt;i++)                                 //  conjugate data in frequency
    {
      xi[i]=-xi[i];
    }
    switch(FFT_type)                                       // Radix-4 or Radix-4 or Radix-8 IFFT
    {
      case 1 :do_R2_FFT(IFFT,R8_level,R4_level,R2_level);
              break;
      case 2 :do_R4_FFT(IFFT,R8_level,R4_level,R2_level);
              break;
      case 3 :do_R8_FFT(IFFT,R8_level,R4_level,R2_level);
              break;
    }
    do_bit_reversal();                                     // bit reversal
    //printf("  after IFFT  \n");                            // show IFFT result
   // show_data();
    for(i=0; i<fft_pt;i++)                                 // show total error
    {
      error+=sqrt(pow(sr[i]-xr[i],2)+pow(si[i]-xi[i],2));
    }
    printf("total error= %f \n",error);
    getch();
}
void get_data()                                            // get randomly data function
{
     int i;
     //printf("\nx:\n");
     for(i=0; i<fft_pt;i++)
     {
		xr[i]=sr[i]=(rand()%10);
		xi[i]=si[i]=(rand()%10);
		//printf("%f   %f\n",xr[i],xi[i]);
	 }
}
void show_data()                                           // show FFT output data function
{
     int i;
     for(i=0; i<fft_pt;i++)
     {
       printf("%f   %f\n",T_xr[i],T_xi[i]);
     }
     printf("\n\n");
}
void but_cell(int fir_index,int sec_index)                 // butterfly function
{
     double tr,ti;
     tr=xr[fir_index]+xr[sec_index];
     ti=xi[fir_index]+xi[sec_index];
     xr[sec_index]=xr[fir_index]-xr[sec_index];
     xi[sec_index]=xi[fir_index]-xi[sec_index];
     xr[fir_index]=tr;
     xi[fir_index]=ti;
}
int bit_reversal(int tr_index,int bit_length)
{
  int i,inv; 
  for(i=0,inv=0;i<bit_length;i++) 
  { 
    inv=(inv<<1)|(tr_index&0x1); 
    tr_index>>=1; 
   } 
   return inv; 
}
void do_bit_reversal()                                     // bit reversal function
{
    int i,j,inv,tmp; 
    for(i=0;i<(fft_pt);i++)
    {
     /* tmp=i;
      for(j=0,inv=0;j<bit_len;j++) 
      { 
        inv=(inv<<1)|(tmp&0x1); 
        tmp>>=1; 
      } */
      T_xr[bit_reversal(i,bit_len)]=xr[i];
      T_xi[bit_reversal(i,bit_len)]=xi[i];
    }
    for(i=1;i<(fft_pt-1);i++)
    {
      xr[i]=T_xr[i];
      xi[i]=T_xi[i];
    }
}
void do_R2_FFT(int IsFFT,int R8_stage,int R4_stage,int R2_stage)// Radix-2 FFT function
{
  int i,j,k;
  int x,y,z,pt_num;
  double a,b,tmp;
  for(x=0;x<R2_stage;x++)                                
  {
    if(R8_stage>0)                                         // judge the last stage of Radix-4(8) FFT 
    {
      z=pow(8,R8_stage);
      pt_num=fft_pt/z;
    }
    else if(R4_stage>0)
    {
      z=pow(4,R4_stage);
      pt_num=fft_pt/z;
    }
    else
    {
      z=pow(2,x);
      pt_num=fft_pt/z;
    }
    for(y=0;y<z;y++)
    {
      for(i=y*pt_num;i<(pt_num/2+y*pt_num);i++)            // first stage of Radix-2 FFT 
      {
        but_cell(i,i+(pt_num/2));                          // butterfly stage
      }
      for(j=1;j<(pt_num/2);j++)                            // complex multiplication
      {
        k=j+(pt_num/2)+y*pt_num;
        if((k >= pt_num/2+y*pt_num)&&(k < pt_num+y*pt_num))
        {
          a=cos(2*pi*j/pt_num);
          b=sin(2*pi*j/pt_num);
        }
        tmp=xr[k];
        xr[k]=xr[k]*a+xi[k]*b;
        xi[k]=xi[k]*a-tmp*b;
      }
    }
  }
  if(IsFFT==0)                                             // judge FFT or IFFT 
  {
    for(i=0;i<(fft_pt/2);i++)
    {
      xr[i*2]=xr[i*2]/fft_pt;
      xi[i*2]=-xi[i*2]/fft_pt;
      xr[i*2+1]=xr[i*2+1]/fft_pt;
      xi[i*2+1]=-xi[i*2+1]/fft_pt;
    }
  }            
}
void do_R4_FFT(int IsFFT,int R8_stage,int R4_stage,int R2_stage)// Radix-4 FFT function
{
  int i,j,k,l;
  int x,y,z,pt_num;
  double a,b,tmp;
  for(x=0;x<R4_stage;x++)                                
  {
    if(R8_stage>0)                                         // judge the last stage of Radix-8 FFT 
    {
      z=pow(8,R8_stage);
      pt_num=fft_pt/z;
    }
    else
    {
      z=pow(4,x);
      pt_num=fft_pt/z;
    }
    for(y=0;y<z;y++)
    {
      for(i=y*pt_num;i<(pt_num/2+y*pt_num);i++)            // first stage of Radix-4 FFT 
      {
        but_cell(i,i+(pt_num/2));                          // butterfly stage
        if((i+(pt_num/2))>=(3*(pt_num/4))+y*pt_num)        // multiply -j
        {
          tmp= xi[i+(pt_num/2)];
          xi[i+(pt_num/2)]= -xr[i+(pt_num/2)];
          xr[i+(pt_num/2)]= tmp;
        } 
      }
      for(i=0;i<2;i++)                                     // second stage of Radix-4 FFT
      {
        for(j=y*pt_num;j<(pt_num/4+y*pt_num);j++)
        {
          but_cell(j+(i*(pt_num/2)),j+(pt_num/4)+(i*(pt_num/2)));// butterfly stage
        }
      }
      for(i=1;i<4;i++)                                     // complex multiplication
      {
        for(j=1;j<(pt_num/4);j++)
        {
          k=j+(pt_num/4)*i+y*pt_num;
          if((k >= pt_num/4+y*pt_num)&&(k < pt_num/2+y*pt_num))
          {
            a=cos(2*pi*2*j/pt_num);
            b=sin(2*pi*2*j/pt_num);
          }
          else if((k >= pt_num/2+y*pt_num)&&(k < pt_num*3/4+y*pt_num))
          {
            a=cos(2*pi*j/pt_num);
            b=sin(2*pi*j/pt_num);
          }
          else if((k >= pt_num*3/4+y*pt_num)&&(k < pt_num+y*pt_num))
          {
            a=cos(2*pi*3*j/pt_num);
            b=sin(2*pi*3*j/pt_num);
          }
          tmp=xr[k];
          xr[k]=xr[k]*a+xi[k]*b;
          xi[k]=xi[k]*a-tmp*b;
        }
      }
    }
  }
  if(R2_stage==1)                                          // judge do Radix-2 FFT algrithm or not
  {
    do_R2_FFT(IsFFT,R8_stage,R4_stage,R2_stage);
  }
  else if(IsFFT==0)                                        // judge FFT or IFFT 
  {
    for(i=0;i<(fft_pt/2);i++)
    {
      xr[i*2]=xr[i*2]/fft_pt;
      xi[i*2]=-xi[i*2]/fft_pt;
      xr[i*2+1]=xr[i*2+1]/fft_pt;
      xi[i*2+1]=-xi[i*2+1]/fft_pt;
    }
  }           
}
void do_R8_FFT(int IsFFT,int R8_stage,int R4_stage,int R2_stage)// Radix-8 FFT function
{
  
  int i,j,k;
  int x,y,z,pt_num,l;
  double a,b,tmp;
  for(x=0;x<R8_stage;x++)                                
  {
    z=pow(8,x);
    pt_num=fft_pt/z;
    for(y=0;y<z;y++)
    {
      for(i=y*pt_num;i<(pt_num/2+y*pt_num);i++)            // first stage of Radix-8 FFT 
      {
        but_cell(i,i+(pt_num/2));                          // butterfly stage
        if((i+(pt_num/2))>=(3*(pt_num/4))+y*pt_num)        // multiply -j
        {
          tmp= xi[i+(pt_num/2)];
          xi[i+(pt_num/2)]= -xr[i+(pt_num/2)];
          xr[i+(pt_num/2)]= tmp;
        } 
      }
      for(i=0;i<2;i++)                                     // second stage of Radix-8 FFT
      {
        l=i*pt_num/2;
        for(j=y*pt_num+l;j<(pt_num/4+y*pt_num+l);j++)
        {
          but_cell(j,j+(pt_num/4));                        // butterfly stage
        }
      }
      for(i=1;i<4;i++)                                     // complex multiplication
      {
        for(j=pt_num/8;j<(pt_num/4);j++)
        {
          
          k=j+(pt_num/4)*i+y*pt_num;
          if((k >= pt_num*3/8+y*pt_num)&&(k < pt_num/2+y*pt_num))
          {
            a=cos(2*pi*2/8);
            b=sin(2*pi*2/8);
          }
          else if((k >= pt_num*5/8+y*pt_num)&&(k < pt_num*3/4+y*pt_num))
          {
            a=cos(2*pi/8);
            b=sin(2*pi/8);
          }
          else if((k >= pt_num*7/8+y*pt_num)&&(k < pt_num+y*pt_num))
          {
            a=cos(2*pi*3/8);
            b=sin(2*pi*3/8);
          }
          tmp=xr[k];
          xr[k]=xr[k]*a+xi[k]*b;
          xi[k]=xi[k]*a-tmp*b;
        }
      }
      for(i=0;i<4;i++)                                     // third stage of Radix-8 FFT
      {
        l=i*pt_num/4;
        for(j=l+y*pt_num;j<(pt_num/8+l+y*pt_num);j++)
        {
          but_cell(j,j+pt_num/8);                          // butterfly stage
        }
      }
      for(i=1;i<8;i++)                                     // complex multiplication
      {
        for(j=1;j<(pt_num/8);j++)
        {
          k=j+(pt_num/8)*i+y*pt_num;
          if((k >= pt_num/8+y*pt_num)&&(k < pt_num/4+y*pt_num))//e
          {
            a=cos(2*pi*4*j/pt_num);
            b=sin(2*pi*4*j/pt_num);
          }
          else if((k >= pt_num/4+y*pt_num)&&(k < pt_num*3/8+y*pt_num))
          {
            a=cos(2*pi*2*j/pt_num);
            b=sin(2*pi*2*j/pt_num);
          }
          else if((k >= pt_num*3/8+y*pt_num)&&(k < pt_num/2+y*pt_num))
          {
            a=cos(2*pi*6*j/pt_num);
            b=sin(2*pi*6*j/pt_num);
          }
          else if((k >= pt_num/2+y*pt_num)&&(k < pt_num*5/8+y*pt_num))
          {
            a=cos(2*pi*j/pt_num);
            b=sin(2*pi*j/pt_num);
          }
          else if((k >= pt_num*5/8+y*pt_num)&&(k < pt_num*3/4+y*pt_num))
          {
            a=cos(2*pi*5*j/pt_num);
            b=sin(2*pi*5*j/pt_num);
          }
          else if((k >= pt_num*3/4+y*pt_num)&&(k < pt_num*7/8+y*pt_num))
          {
            a=cos(2*pi*3*j/pt_num);
            b=sin(2*pi*3*j/pt_num);
          }
          else if((k >= pt_num*7/8+y*pt_num)&&(k < pt_num+y*pt_num))
          {
            a=cos(2*pi*7*j/pt_num);
            b=sin(2*pi*7*j/pt_num);
          }
          tmp=xr[k];
          xr[k]=xr[k]*a+xi[k]*b;
          xi[k]=xi[k]*a-tmp*b;
        }
      }
    }
  }  
  if(R4_stage==1)                                          // judge do Radix-4(2) FFT algrithm or not
  {
    do_R4_FFT(IsFFT,R8_stage,R4_stage,R2_stage);
  }
  else if(R2_stage==1)
  {
    do_R2_FFT(IsFFT,R8_stage,R4_stage,R2_stage);
  }
  else if(IsFFT==0)                                        // judge FFT or IFFT 
  {
    for(i=0;i<(fft_pt/2);i++)
    {
      xr[i*2]=xr[i*2]/fft_pt;
      xi[i*2]=-xi[i*2]/fft_pt;
      xr[i*2+1]=xr[i*2+1]/fft_pt;
      xi[i*2+1]=-xi[i*2+1]/fft_pt;
    }
  }           
}