//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;
	}

}