main.c

采用128点快速傅立叶变换计算256点实信号傅立叶变换

/*------------------------FFT运算C程序--------------------
           时间抽取方式，采用N点复数FFT计算2N点实序列FFT
		                 8月24日调试通过
		               Design By SHUIJIAN
---------------------------------------------------------------*/

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

#define NFFT 32
#define NPOINT 2*NFFT

/*-----------------------------------------------------------------------------------------
函数名称:serial_init()
函数功能:串口初始化函数
入口参数:无
出口参数:无
------------------------------------------------------------------------------------------*/

void serial_init(void)                        //串口初始化函数
{
   TMOD = 0x21;                              //T1工作方式1，T0 16位定时器方式
   TH1  = 0xF3;                              //波特率4800Bps,12MHz晶振
   TL1  = 0xF3;
   PCON |= 0x80;                             //SMOD = 1,波特率加倍
   SCON = 0x50;                              //允许发送接收
   ES   = 0;                                 //禁止串口中断
   TI   = 1;
   TR1  = 1;
}

typedef unsigned char uchar;
typedef unsigned int uint;

#define PI 3.14159
//256点逆序表
/*
int Nixu_Table[]=
{
	0,128,64,192,32,160,96,224,16,
	144,80,208,48,176,112,240,8,
	136,72,200,40,168,104,232,24,
	152,88,216,56,184,120,248,4,
	132,68,196,36,164,100,228,20,
	148,84,212,52,180,116,244,12,
	140,76,204,44,172,108,236,28,
	156,92,220,60,188,124,252,2,
	130,66,194,34,162,98,226,18,
	146,82,210,50,178,114,242,10,
	138,74,202,42,170,106,234,26,
	154,90,218,58,186,122,250,6,
	134,70,198,38,166,102,230,22,
	150,86,214,54,182,118,246,14,
	142,78,206,46,174,110,238,30,
	158,94,222,62,190,126,254,1,
	129,65,193,33,161,97,225,17,
	145,81,209,49,177,113,241,9,
	137,73,201,41,169,105,233,25,
	153,89,217,57,185,121,249,5,
	133,69,197,37,165,101,229,21,
	149,85,213,53,181,117,245,13,
	141,77,205,45,173,109,237,29,
	157,93,221,61,189,125,253,3,
	131,67,195,35,163,99,227,19,
	147,83,211,51,179,115,243,11,
	139,75,203,43,171,107,235,27,
	155,91,219,59,187,123,251,7,
	135,71,199,39,167,103,231,23,
	151,87,215,55,183,119,247,15,
	143,79,207,47,175,111,239,31,
	159,95,223,63,191,127,255
};*/

//128点逆序表
int Nixu_Table[]=
{
	0,64,32,96,16,80,48,112,
	8,72,40,104,24,88,56,120,
	4,68,36,100,20,84,52,116,
	12,76,44,108,28,92,60,124,
	2,66,34,98,18,82,50,114,
	10,74,42,106,26,90,58,122,
	6,70,38,102,22,86,54,118,
	14,78,46,110,30,94,62,126,
	1,65,33,97,17,81,49,113,
	9,73,41,105,25,89,57,121,
	5,69,37,101,21,85,53,117,
	13,77,45,109,29,93,61,125,
	3,67,35,99,19,83,51,115,
	11,75,43,107,27,91,59,123,
	7,71,39,103,23,87,55,119,
	15,79,47,111,31,95,63,127
};

typedef struct compx
{
	 float real;
	 float imag;
}compx;


compx Time_Sig[NPOINT];/*,Freq_Sig[256]*/
//时域信号Time_Sig，频域信号X1,X2,Freq_Sig

float AD_DATA[NPOINT]=
{
  0.0980,    0.1951,    0.2903,    0.3827,    0.4714,   0.5556,    0.6344,    0.7071,
  0.7730,    0.8315,    0.8819,    0.9239,    0.9569,   0.9808,    0.9952,    1.0000,
  0.9952,    0.9808,    0.9569,    0.9239,    0.8819,   0.8315,    0.7730,    0.7071,
  0.6344,    0.5556,    0.4714,    0.3827,    0.2903,   0.1951,    0.0980,    0.0000,
 -0.0980,   -0.1951,   -0.2903,   -0.3827,   -0.4714,  -0.5556,   -0.6344,   -0.7071,
 -0.7730,   -0.8315,   -0.8819,   -0.9239,   -0.9569,  -0.9808,   -0.9952,   -1.0000,
 -0.9952,   -0.9808,   -0.9569,   -0.9239,   -0.8819,  -0.8315,   -0.7730,   -0.7071,
 -0.6344,   -0.5556,   -0.4714,   -0.3827,   -0.2903,   -0.1951,  -0.0980,   -0.0000,

};  //ADC采集到的数据

/*------------------------------------------------------
奇偶分组函数void Sort(void)
将AD采集到的信号进行奇偶分组，并进行倒序运算，以进行2N点到N点的FFT计算
调试通过
--------------------------------------------------------*/
void Sort(void)   //对输入信号进行奇偶分组
{
	 int i;
	 compx Sig[NFFT];
	 float even[NFFT];     //采集数据的偶部
     float odd[NFFT];      //采集数据的奇部
	 for(i=0;i<NPOINT;i++)
	 {
		  if(i%2) odd[i/2] = AD_DATA[i];
		  else    even[i/2]  = AD_DATA[i];
          if(i%2) 
		  {	  
			  Sig[i/2].real = even[i/2];  //实部等于采集信号的偶部
			  Sig[i/2].imag = odd[i/2];   //虚部等于采集信号的奇部
		  }
	 }
	 for(i=0;i<NFFT;i++)
	 {
		 Time_Sig[i] = Sig[Nixu_Table[i]];  //倒序存储
     }
}

/*---------------------------------------------------
复数乘法函数Compx_Mul()
入口参数:复数b1,b2
出口参数:乘法结果b3
---------------------------------------------------*/

compx Compx_Mul(compx b1,compx b2)
{
compx b3;
b3.real=b1.real*b2.real-b1.imag*b2.imag;
b3.imag=b1.real*b2.imag+b1.imag*b2.real;
return(b3);
}


void FFT(compx *xin,int N)
{
    int i,j;     //循环控制变量
	int order,vorder,le,interval,ip;  //蝶形运算阶数,阶数控制变量,蝶形运算间隔控制
	int p;
	float rol;   //旋转因子计算变量
    compx w;           //旋转因子exp(-i*2*PI/N)
	compx temp;        //中间运算变量
	int f = N;
    for(order=1;(f=f/2)!=1;order++);      //计算蝶形运算阶数
	for(vorder=1;vorder<=order;vorder++)  //order级蝶形运算,采用原位运算
	{
         le = pow(2,vorder);  
         interval = le/2;           //各级蝶形运算的间隔点数
		 for(i=0;i<interval;i++)    //控制旋转因子相位循环
		{
		 p = pow(2,order-vorder)*i; //旋转因子w=exp(-j*2*pi/N*p)
		 rol = 2*PI/N*p;            //旋转因子相位
		 w.real = cos(rol);
		 w.imag = -sin(rol); 	  
		 for(j=i;j<N;j+=le)
			  {
				      ip = j + interval;               //当前点与相差间隔点进行蝶形运算
				      temp =  xin[j];
				      xin[ip] = Compx_Mul(xin[ip],w);  //间隔点乘上一个旋转因子
				      xin[j].real =(xin[j].real + xin[ip].real)/8;
				      xin[j].imag =(xin[j].imag + xin[ip].imag)/8;
				      xin[ip].real = (temp.real - xin[ip].real)/8;
                      xin[ip].imag = (temp.imag - xin[ip].imag)/8;
				  }
		 }
	}
}	



void Final_Cal(void)  //计算两个N点FFT,以计算整个序列的最终FFT
{
	 int i;
	 compx X1[NFFT],X2[NFFT];
	 compx w;            //旋转因子等于exp(-i*PI/N)
	 for(i=0;i<NFFT;i++)  //X1[k] = 0.5*(X[k]+conj(X[N-k])) 
	 {                   //X2[k] = -0.5j*(X[k]-conj(X[N-k]))
		 X1[i].real = 0.5*(Time_Sig[i].real+Time_Sig[128-i].real);
		 X1[i].imag = 0.5*(Time_Sig[i].imag-Time_Sig[128-i].imag);
		 X2[i].real = 0.5*(Time_Sig[i].imag+Time_Sig[128-i].imag);
		 X2[i].imag = -0.5*(Time_Sig[i].real-Time_Sig[128-i].real);
	 }
	 for(i=0;i<NFFT;i++)
	 {
		 w.real = cos(i*PI/128);
		 w.imag = -sin(i*PI/128);
		 Time_Sig[i].real = X1[i].real + Compx_Mul(X2[i],w).real;
         Time_Sig[i].imag = X1[i].imag + Compx_Mul(X2[i],w).imag;
		 Time_Sig[i+128].real = X1[i].real - Compx_Mul(X2[i],w).real;
         Time_Sig[i+128].imag = X1[i].imag - Compx_Mul(X2[i],w).imag;
     }
}
 
void main(void)
{
	 int i;
	 serial_init();
	  //float result[256];
	 // FILE *p,*q,*m;
	 // p = fopen("TimeFFT.txt","w");
	 // q = fopen("TimeSig.txt","w");
	 // m = fopen("Nixu.txt","w");
	  /*for(i=0;i<256;i++)
	  {
		    AD_DATA[i]  = sin(i*PI/128);
			//将时域信号输出到文本文档
			fprintf(q,"%.4f,",AD_DATA[i]);
      }*/
	  Sort();
	  //for(i=0;i<128;i++)
		//  fprintf(m,"[%d]=%.4f+%.4fi,",i,Time_Sig[i].real,Time_Sig[i].imag);
	  FFT(Time_Sig,NFFT);
	  Final_Cal();
	  for(i=0;i<NPOINT;i++)
	  printf("real[%d]=%d   imag[%d]=%d\n",i,(int)Time_Sig[i].real,i,(int)Time_Sig[i].imag);
      while(1);
	 /* for(i=0;i<256;i++)
	  {
		  
		  //result[i]=sqrt(pow(Freq_Sig[i].real,2)+pow(Freq_Sig[i].imag,2));
		  //将频域信号输出到TXT文档
         if(Time_Sig[i].imag>0)
		  {
			  fprintf(p,"%.6f+%.6fi,",Time_Sig[i].real,Time_Sig[i].imag);
		  }
		  else if(Time_Sig[i].imag<0)
		  {
			  fprintf(p,"%.6f%.6fi,",Time_Sig[i].real,Time_Sig[i].imag);
		  }
		  else
			  fprintf(p,"%.6f,",Time_Sig[i].real);
      }*/

}