fft.cpp

频域抽取16点基4 FFT c语言代码

#include <stdio.h>
#include <math.h>
#define L 2
//定义结构体由来表示复数
struct realnumber
{
	double re;
	double im;
};
//函数功能:计算4点DFT
void dft(struct realnumber *pd,struct realnumber *px,int h)
{
	struct realnumber dtemp[4];//定义临时数组
	int i;
	for(i=0;i<4;i++)//将要进行DFT的4个点存入临时数组
	{
		dtemp[i].re=pd->re;
		dtemp[i].im=pd->im;
	    pd=pd+h;
	}
	printf("Begin DFT...\n");
	//以下是计算4点DFT的过程
	px->re=dtemp[0].re+dtemp[1].re+dtemp[2].re+dtemp[3].re;
	px->im=dtemp[0].im+dtemp[1].im+dtemp[2].im+dtemp[3].im;
	printf("%8f+%8fj\n",px->re,px->im);
	px=px+h;
	px->re=dtemp[0].re+dtemp[1].im-dtemp[2].re-dtemp[3].im;
	px->im=dtemp[0].im-dtemp[1].re-dtemp[2].im+dtemp[3].re;
	printf("%8f+%8fj\n",px->re,px->im);
	px=px+h;
	px->re=dtemp[0].re-dtemp[1].re+dtemp[2].re-dtemp[3].re;
	px->im=dtemp[0].im-dtemp[1].im+dtemp[2].im-dtemp[3].im;
	printf("%8f+%8fj\n",px->re,px->im);
	px=px+h;
	px->re=dtemp[0].re-dtemp[1].im-dtemp[2].re+dtemp[3].im;
	px->im=dtemp[0].im+dtemp[1].re-dtemp[2].im-dtemp[3].re;
	printf("%8f+%8fj\n",px->re,px->im);
}

//函数功能：计算旋转因子
void weight(struct realnumber *px,int rn)
{
	int m,n,t;
	double re1,im1,retemp;//定义临时变量
	double o;//定义旋转因子的复角
	t=rn/4;
	for(m=0;m<4;m++)
	{
		printf("\n");
		for(n=0;n<t;n++)
	     {
			 o=6.283052-6.283052*m*n/rn;//计算旋转因子的复角
			 re1=cos(o);//计算旋转因子的实部
			 im1=sin(o);//计算旋转因子的虚部
              retemp=(px->re*re1)-(px->im*im1);//计算加权后的实部，原变量保持实部不变
			 px->im=(px->re*im1)+(px->im*re1);//计算加权后的虚部
			 px->re=retemp;//将临时变量赋给原变量的实部
			 px++;//指向下一数据单元			 
		 }
	}
}

//函数功能：输出
void output(struct realnumber *ndata,int n)
{
	int i,j;
	j=0;
printf("\nCongratulations! Now,We get the final answers as follow:\n\n");
	for(i=0;i<n;i++)
	{
		if(j==4)
		{
			printf("\n");
			j=0;
		}
        j++;
		printf("%d|%8f+%8fj|;   ",i,ndata->re,ndata->im);
		ndata++;
	}
}

//函数功能：倒序
void reverse(struct realnumber *pd,struct realnumber *px,int n)
{

    int p=0;
    int i,r,j,k,t,a;
    for(i=L-1;i>0;i--)//确定层数
    {
		p=0;
	    a=(int)pow(4,i);//计算子序列个数
		t=n/a;
		for(r=0;r<a/4;r++)//确定子序列
		{
			for(j=0;j<t;j++)
		       {
			       for(k=0;k<4;k++)
				  {
				     (pd+p)->re=(px+4*t*r+t*k+j)->re;//抽取t*k+j(0<k<4;0<j<t)
				     (pd+p)->im=(px+4*t*r+t*k+j)->im;//抽取t*k+j(0<k<4;0<j<t)
					 p++;
				  }
			  }
		 }
		 for(j=0;j<n;j++)//将倒序后的值存入数据单元
		 {
			 (px+j)->re=(pd+j)->re;
			 (px+j)->im=(pd+j)->im;
		 }
	}

}

//主程序
void main()
{
	int i,j,k,n,r,rn,h;
	struct realnumber nd[16]={{1,0},{1,0},{1,0},{1,0}, 
								{0,0},{0,0},{0,0},{0,0},
								{1,0},{1,0},{1,0},{1,0}, 
								{0,0},{0,0},{0,0},{0,0}};//离散序列x[n]
  struct realnumber xd[16];//定义结果的存储单元
	struct realnumber *ndata,*xdata,*pd,*px,*t,*w;//定义指针
	ndata=nd;
	xdata=xd;
   n=(int)pow(4,L);//计算总点数N
   //input();
	for(i=1;i<=L;i++)//i 控制层数
   {
       pd=ndata;
       px=xdata;
       r=(int)pow(4,i);//计算该层子序个列数
	   r=r/4;
	   w=px; 
	   rn=n/r;//计算子序列长度
       h=rn/4;//计算抽取间隔
       j=0;
go:	   if(j==r)
		   goto next;//r 每层组数
       { 
		   for(k=0;k<h;k++)//每组做4点DFT的次数
           {
     	      // printf("j=%d\n",j);
			   dft(pd,px,h);//调用一次4点DFT函数
               pd++;
               px++;
            }
           for(k=0;k<0.75*rn;k++)//指针指向下一组子序列
		   {
			   pd++;
			   px++;
		   }
           if(i<L)
		   {
			   weight(w,rn);//每组子序列加权，w为起始地址
               w=w+rn;
		   }
	    j++;
	    goto go;
	   }

next:  t=ndata;//将计算结果作为下一层运算的输入
       ndata=xdata;
       xdata=t;
     }
	reverse(xdata,ndata,n);//最终计算结果的倒序
   output(ndata,n);//输出
}