FFT 算法实现




实现功能:
DIF Radix2,  DIF Radix4
DIT Radix2,  DIT Radix4
FFT complex, IFFT complex
FFT real, IFFT real

使用方法:
lguo@lguo-enst:~$ unzip Algo.zip
lguo@lguo-enst:~$ cd Algo/DIF
lguo@lguo-enst:~/Algo/DIF$ gcc -o FFT_DIF_R2 FFT_DIF_R2.c -lm
lguo@lguo-enst:~/Algo/DIF$ ./FFT_DIF_R2

lguo@lguo-enst:~$ cd Algo/FFT_Test_Linux
lguo@lguo-enst:~/Algo/FFT_Test_Linux$ make
lguo@lguo-enst:~/Algo/FFT_Test_Linux$ ./FFT

说明:
Radix2 可以计算 4,8,16,32, 64,128, 256....点FFT
Radix4 可以计算 4,16, 64, 256, 1024...点FFT
FFT_DIT_general.c 实现了 Radix2 和Radix4 的配合使用,可以计算Radix2可以计算的所有FFT,但效率比Radix2高。

FFT complex 和 FFT real:
FFT real 和 IFFT real 是指输入信号没有虚部,只是实数序列的情况。这里就可以一个N/2点的FFT complex 来计算 N点FFT real。
特别提醒的是 IFFT real 也是指输入信号没有虚部,就是频域信号没有虚部,而输出信号,也就是时域信号有虚部的情况。这在实际应用中不多见。
之后进一步的工作就是要,实现了输入信号有虚部,输出是实序列的IFFT real。这样和FFT complex 组合起来就可以做成一个常用的,时域变到频域,频域处理后再变回时域的应用。

代码下载:
http://www.cppfrance.com/codes/ALGORITHME-FFT_44369.aspx

代码示例:
/***********************************
// DIT radix-4  FFT  complex 
//
// 1. Maximium points are 2048 because that 
//    function mylog2(int N) has the limit of maximal points
//
// 2. The last stage is 2 DFT ( for 8, 32, 128, 512, 2048...)
// or 4 DFT ( for 4, 16, 64, 256, 1024 ...).
//
//
// 24 juillet 2007
// purcharse*gmail.com
***********************************
*/


#include 
<math.h>
#include 
<stdio.h>
#include 
<stdlib.h>

typedef 
double real64;      /* floating point */
typedef 
float  real32;    /* 32 bit fixed point */

struct complex

  real32 real;
  real32 imag;
}
 complex ;

static struct complex multicomplex(struct complex,struct complex);
static int mylog2(int); 

   static void DFT_2(struct complex *,struct complex *);
static void DFT_4(struct complex *,struct complex *,struct complex *,struct complex *);

static void RunFFT(struct complex *int);
static void RunIFFT(struct complex *,int);

static void BitReverse(struct complex *int);
static void FFT_R4(struct complex *intint);
static void FFT_L2(struct complex *int);

struct  complex s[257];
int   Num = 16;
const float PI=3.1415926535898;

main()
{
    
int i;

    
/* rectangle */
      
/*
      for(i=0;i<Num+1;i++)
      {
          s[i].real=0;
          s[i].imag=0;
      }

    s[0].real = 10;
    s[0].imag = 10;
    
*/


    
/*  sinus*/
    
      
for(i=0;i<Num+1;i++)
      
{
          s[i].real
=sin(PI*i/Num);
          s[i].imag
=cos(PI*i/Num);
      }
 
    
    
/*  */
    
/*
      for(i=0;i<Num+1;i++)
      {
          s[i].real=0;
          s[i].imag=0;
      }
      for(i=Num*3/8;i<Num/2;i++)
      {
          s[i].real=(i-Num*3/8);
          s[i].imag=0;
      }
      for(i=Num/2;i<Num*5/8;i++)
      {
          s[i].real=(Num*5/8-i);
          s[i].imag=0;
      }       
*/
   

    printf(
"*********** Donnees ***************** ");
    
for(i=0;i<Num;i++)
    
{
        printf(
"%.4f ",s[i].real);
        printf(
"%.4f ",s[i].imag);
    }

    
    RunFFT(s, Num);
    printf(
"*********** FFT ***************** ");
    
for(i=0;i<Num;i++)
    
{
        printf(
"%.4f ",s[i].real);
        printf(
"%.4f ",s[i].imag);
    }

    RunIFFT(s, Num);
    printf(
"*********** IFFT ***************** ");
    
for(i=0;i<Num;i++)
    
{
        printf(
"%.4f ",s[i].real);
        printf(
"%.4f ",s[i].imag);
    }

}

/***************** functions ********************/

static struct complex multicomplex(struct complex b1,struct complex b2)  /* multiplication of complex */
{
struct complex b3;
b3.real
=b1.real*b2.real-b1.imag*b2.imag;
b3.imag
=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}


static int mylog2(int N)     /* Max(N) = 4098 */
{
  
int k=0;
  
  
if (N>>12{ k+=12; N>>=12; }
  
if (N>>8{ k+=8; N>>=8; }
  
if (N>>4{ k+=4; N>>=4; }
  
if (N>>2{ k+=2; N>>=2; }
  
if (N>>1{ k+=1; N>>=1; }
  
return k ;
}



static void BitReverse(struct complex *xin, int N)
{
    
int LH, i, j, k;
    
struct complex tmp;

            LH
=N/2;    
            j 
= N/2;

            
for( i = 1; i <= (N -2); i++)
            
{
                
if(i < j)
                
{
                    tmp 
= xin[j];
                    xin[j] 
= xin[i];
                    xin[i] 
= tmp;
                }

            k 
= LH;
                
while(j >= k)
                
{
                    j 
= j-k;
                    k 
= k/2;
                }


            j 
= j + k;

            }


}


static void DFT_2(struct complex *b1,struct complex *b2)
{
struct complex tmp;
  tmp 
= *b1;

  (
*b1).real = (*b1).real + (*b2).real;
  (
*b1).imag = (*b1).imag + (*b2).imag;

  (
*b2).real = tmp.real - (*b2).real;
  (
*b2).imag = tmp.imag - (*b2).imag;    
}


static void DFT_4(struct complex* b0, struct complex* b1, struct complex* b2,struct complex* b3)
{
  
/*variables locales*/
  
struct complex temp[4];
  
  
/*calcul x1*/
  temp[
0].real=(*b0).real+(*b1).real;        
  temp[
0].imag=(*b0).imag+(*b1).imag;

  
/*calcul x2*/
  temp[
1].real=(*b0).real-(*b1).real;        
  temp[
1].imag=(*b0).imag-(*b1).imag;

  
/*calcul x3*/
  temp[
2].real=(*b2).real+(*b3).real;    
  temp[
2].imag=(*b2).imag+(*b3).imag;

  
/*calcul x4 + multiplication with -j*/
  temp[
3].imag=(*b3).real-(*b2).real;        
  temp[
3].real=(*b2).imag-(*b3).imag; 

  
/*the last stage*/
  (
*b0).real=temp[0].real+temp[2].real;
  (
*b0).imag=temp[0].imag+temp[2].imag;
  
  (
*b1).real=temp[1].real+temp[3].real;     
  (
*b1).imag=temp[1].imag+temp[3].imag; 
  
  (
*b2).real=temp[0].real-temp[2].real;     
  (
*b2).imag=temp[0].imag-temp[2].imag;
  
  (
*b3).real=temp[1].real-temp[3].real; 
  (
*b3).imag=temp[1].imag-temp[3].imag;
}


static void FFT_R4(struct complex *xin, int N, int m)
{
        
int  i, L, j;
        
double  ps1, ps2, ps3, p ;
        
int le,B;
        
struct complex w[4];

        
for( L = 1; L <= m; L++)
        
{
            le 
= pow(4 ,L);
            B 
= le/4;       /*the distance of buttefly*/
        
            p 
= N/le;

             
for(j = 0; j <= B-1 ; j++)    
              
{
                   
//  ps0 = (2*pp/N) * 0 * j;
                   
//  w[0].real = cos(ps0);
                   
//  w[0].imag = -sin(ps0);
    
                   ps1 
= (2*PI/N)*p*2*j;
                   w[
1].real = cos(ps1);
                   w[
1].imag = -sin(ps1);

                   ps2 
= (2*PI/N)*p*1*j;
                   w[
2].real = cos(ps2);
                   w[
2].imag = -sin(ps2);

                   ps3 
= (2*PI/N)*p*3*j;
                   w[
3].real = cos(ps3);
                   w[
3].imag = -sin(ps3);

                   
for(i = j; i <= N-1; i = i + le)    /* controle those same butteflies*/
                     
{
                        
/* multiple with W */
                 
//   xin[i] = multicomplex(xin[i], w[0]);
                      xin[i + B] = multicomplex(xin[i + B], w[1]);
                      xin[i 
+ 2*B] = multicomplex(xin[i + 2*B], w[2]);
                      xin[i 
+ 3*B] = multicomplex(xin[i + 3*B], w[3]);
                        
/* DFT-4 */
                      DFT_4(xin 
+ i, xin + i + B, xin + i + 2*B, xin + i + 3*B);
                     }

                }

        
/*
        printf("*****N°%d ********** ", L);
        for(i=0;i<Num;i++)
        {
        printf("%.8f ",xin[i].real);
        printf("%.8f ",xin[i].imag);
        }
        
*/

        }

    
}
//fin du FFT_R4

static void FFT_L2(struct complex *xin, int N)
{                /* For the last stage 2 DFT*/
    
int j, B;
    
double p, ps ;
    
struct complex w;

        B 
= N/2;               
        
for(j = 0; j <= B - 1; j++)
        
{            
            ps 
= (2*PI/N)*j;
                w.real 
= cos(ps);
            w.imag 
= -sin(ps);

                 
/* multiple avec W */
                     xin[j
+ B] = multicomplex(xin[j + B], w);
             DFT_2(xin 
+ j ,xin + j + B);
        }

                
}
//fin du FFT_L2

static void RunFFT(struct complex *xin, int N)
{
    
int m, i;
    
    BitReverse(xin, N);
    m 
= mylog2(N);
    
if( (m%2== 0 )
    
{
        
/*All the stages are radix 4*/
        FFT_R4(xin, N, m
/2);
    }

    
else 
    
{
        
/*the last stage is radix 2*/
        FFT_R4(xin, N, m
/2);
        FFT_L2(xin, N);
    }

}


static void RunIFFT(struct complex *xin,int N)
{            /* inverse FFT */
      
int i;
    
      
for(i=0; i < Num + 1 ; i++)
      

    xin[i].imag 
= -xin[i].imag;
      }

      RunFFT(xin,N);
      
for(i = 0; i < Num + 1 ; i++)
      

    xin[i].real 
= xin[i].real/Num;
    xin[i].imag 
= -xin[i].imag/Num;
      }

}
阅读更多
想对作者说点什么? 我来说一句

FFT算法实现

2014年04月22日 54KB 下载

FFT算法及程序vc++程序实现

2009年04月24日 247KB 下载

没有更多推荐了,返回首页

加入CSDN,享受更精准的内容推荐,与500万程序员共同成长!
关闭
关闭