fftforc51.c

51上写的浮点FFT程序

#include <REG51F.H>
#include <math.h>
#define number 16
#define N 4
unsigned int powern[N];
unsigned int bitreverse[number];
unsigned char xdata i,j,k,r,ll,M=0;
unsigned int xdata twidnumber;
unsigned int xdata butterinputstep;
unsigned int xdata butterstep;
unsigned int xdata indextemp;
unsigned int xdata index;
unsigned int xdata twidbutternumber;
double tempv;
double PI;
typedef struct
{
	double rr;
	double ii;
}Complex;

Complex  xdata tempk1;
Complex  xdata tempk2;
double xdata wcos[number];
double xdata wsin[number];
//****************************************
Complex xdata x[number];

init();
//*******************************************

//*********************************************
main()
{
	init();
//**************************************************
for(i=0;i<number;i++)
{
	tempv=(double)(1+sin(PI*i/number)+cos(PI*3*i/number));
	x[i].rr=tempv;//偶序列
	x[i].ii=0;//奇序列
}
//*******************************************
// generate Twiddle Factor for N point FFT
for(i=0;i<number;i++)
{
	tempv=(2*PI*i);
  wcos[i]=cos(tempv/number);
  wsin[i]=-sin(tempv/number) ;
}

//****************************************
	powern[0]=1;
	for(i=1;i<N;i++)
	{
		powern[i]=2*powern[i-1];
	}
////**************************************
	for(i=0;i<number;i++)
	{
		bitreverse[i]=0;
		for(j=0;j<N;j++)
		{
			if(i&powern[j])
		  bitreverse[i]+=powern[N-1-j];
		}
	}
////**********************************************
// generate bit reverse input
   for(i=0;i<number;i++)
   {
   	if(i<bitreverse[i])
   	{
   	tempk1=x[i];
    x[i]=x[bitreverse[i]];
    x[bitreverse[i]]=tempk1;
    }
  }

//*************************************************
	while(1)
	{
   for(M=1;M<(N+1);M++)//M为运算级数
   {
   	i=0;
//**********************************************
//        	twidnumber=2^(M-1);//当前运算级旋转因子个数
//        	butterinputstep=(2^(M-1));//当前运算级蝶形结输入间距
//        	butterstep=(2^M);//当前运算级蝶形结间距
//        	indextemp=r*(2^(N-M));//当前运算级旋转因子索引变化量
//   	      twidbutternumber=(2^(N-M));//当前运算级同一旋转因子对应蝶形结个数
//*************************************************
   	twidnumber=1;//当前运算级旋转因子个数
    twidbutternumber=1;//当前运算级同一旋转因子对应蝶形结个数
  	butterinputstep=1;//当前运算级蝶形结输入间距
   	butterstep=1;//当前运算级蝶形结间距
   	indextemp=1;//当前运算级旋转因子索引变化量
   	for(j=0;j<(M-1);j++)//M=1;twidnumber=1;M=2;twidnumber=2;M=3;twidnumber=4;
   	{
   	twidnumber=2*twidnumber;
  	butterinputstep=(2*butterinputstep);//M=1;butterinputstep=1;M=2;butterinputstep=2;M=3;butterinputstep=4;
    }
   	for(j=0;j<(M);j++)////M=1;butterstep=2;M=2;butterstep=4;M=3;butterstep=8;
   	{
   	butterstep=(2*butterstep);
   	indextemp=indextemp*2;
    twidbutternumber=(2*twidbutternumber);
    }
    indextemp=number/indextemp;
    twidbutternumber=(number/twidbutternumber);//M=1;twidbutternumber=4;M=2;twidbutternumber=2;M=3;twidbutternumber=1;
   	for(r=0;r<twidnumber; )//r为第M级旋转因子索引
        {
        	index=r*(indextemp);
        	for(k=0;k<twidbutternumber;k++)
        	   {

tempk1.rr=x[i+butterinputstep].rr*wcos[index]-x[i+butterinputstep].ii*wsin[index];
tempk1.ii=x[i+butterinputstep].rr*wsin[index]+x[i+butterinputstep].ii*wcos[index];
x[i+butterinputstep].rr=tempk1.rr;
x[i+butterinputstep].ii=tempk1.ii;
//****************************************************************
tempk1.rr=x[i].rr+x[i+butterinputstep].rr;
tempk1.ii=x[i].ii+x[i+butterinputstep].ii;
//************************************************************
tempk2.rr=x[i].rr-x[i+butterinputstep].rr;
tempk2.ii=x[i].ii-x[i+butterinputstep].ii;
//*************************************************************
x[i].rr=tempk1.rr;
x[i].ii=tempk1.ii;
x[i+butterinputstep].rr=tempk2.rr;
x[i+butterinputstep].ii=tempk2.ii;
//***************************************************************
        	    i=i+butterstep;
        	   }
        	i=++r;
        }
   }
   for(ll=0;ll<number;ll++)
   {
   	 tempk1.rr=x[ll].rr*x[ll].rr;
     tempk1.ii=x[ll].ii*x[ll].ii;
     x[ll].rr=(sqrt(tempk1.rr+tempk1.ii));
   }
   ll=0;
 }
}

init()
{

    TMOD=0x11;//T0,T1为16位定时器
    EA=0;
    ET0=0;//定时器0中断允许
    ET1=0;//定时器1中断允许
    EX0=0;
    PS=1;//串口中断为高级中断
    ES=0; //控制cebotelv(void)和jieshou(void) interrupt 4//串口中断，只当在查询到正确的波特率后，才在cebotelv(void)中使ES＝1，使其进入jieshou(void) interrupt 4中断服务子程序
    TH0=0xfc;//1ms中断1次**以1ms为最小单位，送20次，每次送6个数
    TL0=0x2b;
    TH1=0x63;//中断时间为40ms
    TL1=0xc0;
    TR1=0;
    IT0=0;//外部中断为边沿触发,下降沿有效
    RCLK=0;//T2用作波特率发生器，接收与发射的波特率都由t2产生
    TCLK=0;
    SCON=0x50;
    TR2=0;
    PI=atan(1)*4;
    i,j,k,r,M=0;
    twidnumber=0;
    butterinputstep=0;
    butterstep=0;
    index=0;
    twidbutternumber=0;
}
//*********************