dif_fft.cpp

按频率抽取计算FFT的C++源代码

//dif_fft.cpp
//本程序为按频率抽取法计算序列FFT（16点）

#include "iostream.h"
#include "math.h"

//定义复数类
class complex
{
public:
	complex()
	{real=imag=0.0;}
	complex(double x)
	{real=x;imag=0.0;}
	complex(double x, double y)
	{
		real=x;
		imag=y;
	}
//运算符重载
	complex operator +(const complex &c)
	{return complex(real+c.real,imag+c.imag);}
	complex operator -(const complex &c)
	{return complex(real-c.real,imag-c.imag);}
	complex operator *(const complex &c)
	{return complex(real*c.real-imag*c.imag,real*c.imag+imag*c.real);}
	complex operator /(const complex &c)
	{return complex((real*c.real+imag*c.imag)/(c.real*c.real+c.imag*c.imag),
	                (imag*c.real-real*c.imag)/(c.real*c.real+c.imag*c.imag));}
	
	void print()   
	{
		if(imag<0)
			cout<<real<<imag<<"i";
		else
			cout<<real<<"+"<<imag<<"i";
	}

private:
	double real,imag;

};


	

//DIF_FFT算法子函数
void dif_fft(complex *X,unsigned N,unsigned L)
{      //X为序列数组，n为序列长度，L=log2 N
	const double pai=3.141592653589793;    //定义常数π
	complex W,t;                         //
	unsigned Nv2, k, interval, group, nm1;
	unsigned i,j,m;
	double r;
	const double angle=2*pai/N;  
	//进行L级递推运算求X[k]
	for(m=1;m<=L;m++)           //第m级
	{
		group=pow(2,m-1);       //group为在第m级的交叉碟形单元组组数
		interval=N/pow(2,m);    //interval为蝶形单元两节点间的节距
		for(i=0;i<group;i++)         //计算第i组交叉蝶形单元组的蝶形运算
		{
			for(j=0;j<interval;j++)       //计算第i组交叉蝶形单元组的蝶形运算
			{
				k=j+i*interval*2;
				r=k<<(m-1);
				W=complex(cos(angle*r),-sin(angle*r));
				t=X[k];
				X[k]=X[k]+X[k+interval];
				X[k+interval]=(t-X[k+interval])*W;
			}
		}
	}			


	//对输出序列X[k]进行倒位序处理____雷得算法
	Nv2=N>>1;     //Nv2=N/2
	nm1=N-1;
	j=0;
	for(i=0;i<nm1;i++)
	{
		if(i<j)
		{
			t=X[j];
			X[j]=X[i];
			X[i]=t;
		}
		k=Nv2;
		while(k<=j)
		{
			j-=k;
			k>>=1;
		}
		j+=k;
	}
}


void main()
{
	const unsigned n=16, L=4;    //定义为16点
	complex X[n];
	double T=0.03125;
	for(int i=0;i<n;i++)
	{
		X[i]=complex(exp(-i*T));
	}

	cout<<"输入序列x[n]为：";
	for(i=0;i<n;i++)
	{
		cout<<"x["<<i<<"]=";
		X[i].print();
		cout<<endl;
	}


	dif_fft(X, n, L);     //进行傅立叶变换FFT

	cout<<"输出序列X[K]为：";
	for( i=0;i<n;i++)
	{
		cout<<"X["<<i<<"]=";
		X[i].print();
		cout<<endl;
	}

}