fft-dif.txt

自己编写的C语言基于频率抽取的基2FFT算法

#include"conio.h"
#include"math.h"
#include"stdio.h"
#define N 8
#define PI 3.1415926
int n;

/*复数定义*/
typedef struct
{
	float re;
	float im;
}COMPLEX;

/*复数加法运算*/
COMPLEX Add(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re+c2.re;
	c.im=c1.im+c2.im;
	return c;
}

/*复数减法运算*/
COMPLEX Sub(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re-c2.re;
	c.im=c1.im-c2.im;
	return c;
}

/*复数乘运算*/
COMPLEX Mul(COMPLEX c1, COMPLEX c2)
{
	COMPLEX c;
	c.re=c1.re*c2.re-c1.im*c2.im;
	c.im=c1.re*c2.im+c2.re*c1.im;
	return c;
}
/*复数赋值运算*/
COMPLEX Equal(COMPLEX c1)
{
        COMPLEX c;
        c.re=c1.re;
        c.im=c1.im;
        return c;
}
/* 旋转因子运算*/
COMPLEX Wn(int r)
{       COMPLEX w;
        w.re=cos(PI/r);
        w.im=-sin(PI/r);
        return w;
}
/*序列的初始化，如果序列点数不是2的整数幂，补零*/
int CSH(COMPLEX *a)
{
       int n1,i;
       int m;
       float m1;
       printf("\n This programme is about FFT by DIF way. ");
       printf("\n please enter N: ");
       scanf("%d",&n1);
       n=n1;
       m1=log(n1)/log(2);
       m=log(n1)/log(2);
       if(m!=m1) n=pow(2,m+1);
       for(i=0;i<n;i++){a[i].re=a[i].im=0.0;} /* bu ling*/
       printf("\n");
       for(i=0;i<n1;i++)
       {
        printf("\nplease enter data(%d)_[Re]:",i);
        scanf("%f",&a[i].re);
        printf("\nplease enter data(%d)_[Im]:",i);
        scanf("%f",&a[i].im); 
       } 
       printf("\n输入序列为：");
       for(i=0;i<n;i++)
       { 
       printf(" %7.3f",a[i].re);
       if(a[i].im>=0)
        printf("+%7.3fj",a[i].im); 
       else
        printf("-%7.3fj",fabs(a[i].im));      
       }     
}
/*快速傅立叶变换*/
int FFT(COMPLEX *a)
{
       int i,j,nu2,nm1,m,le,le1,k,ip;
       int n1;
       COMPLEX u,w,t;       
        /*蝶形运算*/      
        for(m=1;m<=(log(n)/log(2));m++)
        {
         le=n/(pow(2,m-1));
         le1=le/2;
         u.re=1;
         u.im=0;
         w=Wn(le1);
         for(j=0;j<le1;j++)
         {
          for(i=j;i<n;i+=le)
          {
           ip=i+le1;
           t=Sub(a[i],a[ip]);
           a[i]=Add(a[i],a[ip]);
           a[ip]=Mul(t,u);
          }
          u=Mul(u,w);
         }
        }
		/*倒位序*/
       nm1=n-2;
       nu2=n/2;
       i=j=0;
       while(i<=nm1)
        {
          if(i<j)
          {
           t=Equal(a[j]);
           a[j]=Equal(a[i]);
           a[i]=Equal(t);
          }
          k=nu2;
          while(k<=j)
          {
          j-=k;
          k=k/2;
          }
          j+=k;
          i++;       
         } 
        
}

/*快速傅立叶反变换，利用快速傅立叶变换*/
void IFFT(COMPLEX *a)
{
    int i;
	/*求频域点的共轭*/
	for(i=0;i<n;i++)
	{
		a[i].im=-a[i].im;
	}
	/*调用快速傅立叶变换*/
	FFT(a);
	/*求时域点的共轭*/
	for(i=0;i<n;i++)
	{
		a[i].re/=n;
		a[i].im=-a[i].im/n;
	}
     
}

void main()
{
	int i;
    COMPLEX a[N];
    CSH(a);
    printf("\n**********************************************\n");
    FFT(a);
    printf("\nfft变换后序列为：");
       for(i=0;i<n;i++)
       { 
       printf(" %7.3f",a[i].re);
       if(a[i].im>=0)
        printf("+%7.3fj",a[i].im); 
       else
        printf("-%7.3fj",fabs(a[i].im));      
       } 
	   printf("\n**********************************************\n");   
    IFFT(a); 
    printf("\n用ifft计算出的原序列为：");
       for(i=0;i<n;i++)
       { 
       printf(" %7.3f",a[i].re);
       if(a[i].im>=0)
        printf("+%7.3fj",a[i].im); 
       else
        printf("-%7.3fj",fabs(a[i].im));  
       }      
}