FFT快速傅里叶代码

#include <math.h>　　　
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
#include <iostream>
using namespace std;
using std::ofstream;
using std::endl;
#ifndef FFT_Def
#define FFT_Def


#define PI 3.1415926535
#define uchar unsigned char


#define N 64						//FFT点数（长度）
#define BIT 16						//FFT的位数（每个数据点的类型），只能是8、16、32

#if (BIT==8)						//8-bit 有关常量
	#define BITTYPE char			//char代表16位整数
	#define OVS 52					//最大不溢出幅度，如果计算过程中超过此幅度程序将所有点除2
	#define DOVS 104				//最大不溢出幅度*2，如果计算过程中超过此幅度程序将所有点除4
	#define WNS 0x7F				//旋转因子，也就是此位数下的最大整数
	#define SHN 7					//移位因子，当乘以旋转因子后，数值被扩大了2^SHN，因此需要右移

#elif (BIT==16)
	#define BITTYPE short
	#define OVS 13572
	#define DOVS 27145
	#define WNS 0x7FFF
	#define SHN 15 

#elif (BIT==32)
	#define BITTYPE int
	#define OVS 889516850
	#define DOVS 1779033700
	#define WNS 0x7FFFFFFF
	#define SHN  31

#endif
//提示：
//OVS的值是通过设想最坏情况下，在输出不超过此类型允许数值时，加乘运算输入端的最大幅值
//此值的大小是：此类型的最大正数/(1+sqrt(2))


struct Complex{							//构造复数结构
BITTYPE real;
BITTYPE imag;
};
 Complex data[N];
 Complex wn[N/2];
 
 float sin_tab[N/4+1];
 void tab_sin()
{
	int i;	
    sin_tab[0]=0;

	for(i=1;i<=(N/4);i++)//第一象限
	{
		sin_tab[i]=sin(2*PI/N*i);
	 
	}
}

 void tab_wn()
{   
	tab_sin();
	for(int k=0;k<N/4;k++)
	{
		wn[k].real=(BITTYPE)(sin_tab[N/4-k]*WNS);
		wn[k].imag=(BITTYPE)(-sin_tab[k]*WNS);
	}
		for(int k=N/4;k<N/2;k++)
	{
		wn[k].real=(BITTYPE)(-sin_tab[k-N/4]*WNS);
		wn[k].imag=(BITTYPE)(-sin_tab[N/2-k]*WNS);
	}
		
}

 struct Complex_float{							//构造复数结构
float real;
float imag;
};
Complex_float data1[N];

typedef struct Complex COMPLEX;

/*序列的倒序*/
void ReverseIndex(COMPLEX *x)
{
int LH,f,m,nm,j,i,k;
struct Complex t;
LH=N/2;
f=N;
for(m=1;(f=f/2)!=1;m++){;}/*确定FFT运算级数M*/
nm=N-2; /*顺序数i的终止值*/  
j=N/2;/*M位二进制最高位的权值*/

for(i=1;i<=nm;i++)
{
if(i<j){t=x[j];x[j]=x[i];x[i]=t;}
k=LH;
while(j>=k){j=j-k;k=k/2;}
j=j+k;
}
}

//按时间抽取法的fft变换，输入反序，输出正序
void fft(COMPLEX *x){				//x为欲变换数组，也当作输出数组，N为点数，支持1024（其它的点数请自行修改反转下标函数即可）
	int i,j,k,u=0,l=0,wi=0;				//j第二层循环（子块中的每个蝶形的循环计数）
										//k第一层循环（横向fft变换阶数，为log2（N）
										//u 蝶形上标x[upper],l 蝶形下标x[lower]，wi旋转因子下标wn[wi]
	int SubBlockNum,SubBlockStep=1;		//SubBlockNum当前k层子块数量，SubBlockStep当前k层不同子块的相同位置元素间间隔
	bool ov,ovd;							//溢出标志
	COMPLEX XkWn;				//wn为旋转因子数组，XkWn为临时存储当前蝶形的旋转因子的临时变量			
	tab_wn();//初始化wn数组，在嵌入式设备中，此数组请预先计算并写成常量数组，以便于节省时间
	
	ReverseIndex(x);					//输入反序
	for(k=N;k>1;k=(k>>1)){				//第一个循环，代表log2(k)阶的变换
		ov=0;ovd=0;						//溢出缩放检查
		for(i=0;i<N;i++)				//遍历查找是否溢出
			if (x[i].real>DOVS || x[i].imag>DOVS || x[i].real<-DOVS || x[i].imag<-DOVS)				//超出两倍允许值，所有值右移2个bit
			{
				ovd=1;
				break;
			}
			else if (x[i].real>OVS || x[i].imag>OVS || x[i].real<-OVS || x[i].imag<-OVS)			//超出允许值，所有值右移1个bit
			{
				ov=1;
			}
		if (ovd==1)
			for(i=0;i<N;i++)
			{
				x[i].real=x[i].real>>2;
				x[i].imag=x[i].imag>>2;
			}
		else if(ov==1)
			for(i=0;i<N;i++)
			{
				x[i].real=x[i].real>>1;
				x[i].imag=x[i].imag>>1;
			}
		SubBlockNum=k>>1;				//子块个数为所做点数的一半
		SubBlockStep=SubBlockStep<<1;	//子块间同等地位的元素间隔以2为倍数递增
		wi=0;							//旋转因子初始化
		for(j=0;j<SubBlockStep>>1;j++){	//第二层循环，更新j值，j表示各个子块的第j个蝶形。因为每个子块的同地位蝶形具有相同的wn，所以用第二层循环控制wn
			for(u=j;u<N;u+=SubBlockStep){		//第三层循环，循环于各个子块间的第j个蝶形，计算所有蝶形。直到下标u越界。(u>N)		
						l=u+(SubBlockStep>>1);	//下标l计算
						#if (BIT==32)			//注：__int64加入的目的是保留32位相乘结果的高32位，不加入则高32位会忽略。此句需要针对平台修改。
						XkWn.real=(((__int64)x[l].real*wn[wi].real)>>SHN)-(((__int64)x[l].imag*wn[wi].imag)>>SHN);//蝶形x[u]=x[u]+x[l]*Wn,x[l]=x[u]-x[l]*Wn的复数计算
						XkWn.imag=(((__int64)x[l].imag*wn[wi].real)>>SHN)+(((__int64)x[l].real*wn[wi].imag)>>SHN);
						#else
						XkWn.real=((x[l].real*wn[wi].real)>>SHN)-((x[l].imag*wn[wi].imag)>>SHN);//蝶形x[u]=x[u]+x[l]*Wn,x[l]=x[u]-x[l]*Wn的复数计算
						XkWn.imag=((x[l].imag*wn[wi].real)>>SHN)+((x[l].real*wn[wi].imag)>>SHN);
						#endif
						x[l].real=x[u].real-XkWn.real;
						x[l].imag=x[u].imag-XkWn.imag;
						x[u].real=x[u].real+XkWn.real;
						x[u].imag=x[u].imag+XkWn.imag;
			}
			wi+=SubBlockNum;			//第二层循环更新wi值
		}
	}

}

#endif

/*幂函数*/
int power(int a,int b) 
{ 
int i,j;  
j=1;
for(i=0;i<b;i++)
j*=a;
return j;
}

//开平方整型表示
unsigned int sqrt_16(unsigned long M) 
{ 
     unsigned int Nm, i; 
     unsigned long tmp, ttp;    // 结果、循环计数 
     if (M == 0)                // 被开方数，开方结果也为0 
         return 0; 

     Nm = 0; 

     tmp = (M >> 30);           // 获取最高位：B[m-1] 
     M <<= 2; 
     if (tmp > 1)               // 最高位为1 
     { 
         Nm ++;                  // 结果当前位为1，否则为默认的0 
         tmp -= Nm; 
     } 

     for (i=15; i>0; i--)       // 求剩余的15位 
     { 
         Nm <<= 1;               // 左移一位 

         tmp <<= 2; 
         tmp += (M >> 30);      // 假设 

         ttp = Nm; 
         ttp = (ttp<<1)+1; 

         M <<= 2; 
         if (tmp >= ttp)        // 假设成立 
         { 
             tmp -= ttp; 
             Nm ++; 
         } 

     } 
    return Nm;
} 


int ReadFile(void)//读取待处理数据
{
  ifstream DataFileIn;
  
  DataFileIn.open("E:\\FFT\\1-5.txt",ios::in);
  
  for(int i=0; i<N; i++)
  {
    DataFileIn>>data[i].real;
	data[i].imag=0;
	
  }

 DataFileIn.close();
  return 1;
}

void main()
{
ReadFile();	
	
fft(data);
for(int j=0;j<N;j++)
{
	
	data[j].real=sqrt_16(power(data[j].real,2)+power(data[j].imag,2));
printf("%d\n",data[j].real);}
ofstream ofs("E:\\FFT\\result.txt");
for( int j=0;j<N;j++)
{
	ofs<<data[j].real<<endl;
}
ofs.close();
}