#include"conio.h"
#include"iostream.h"
#include"math.h"
#include"stdio.h"
#define N 64
#define PI 3.1415926
#define w0 (0.125 *PI)
#define Cmul(a,b,c) a.x=b.x*c.x-b.y*c.y;a.y=b.x*c.y+b.y*c.x;
#define Cequal(a,b) a.x=b.x;a.y=b.y;
#define Cadd(a,b,c) a.x=b.x+c.x;a.y=b.y+c.y;
#define Csub(a,b,c) a.x=b.x-c.x;a.y=b.y-c.y;
#define Wn(w,r) w.x=cos(2.0*PI*r/n);w.y=-sin(2.0*PI*r/n);
struct comp
{
 float x;
 float y;
};
void main()
{
	int i,j,nu2,nm1,n,m,le,le1,k,ip,z;
	int flag,f,n1;
	struct comp a[N],t,t1,w,d;
	float a_ipx,m1;
	cout<<"這個是按頻率抽選的FFT算法:";
	cout<<'\n';
	cout<<"請輸入N:";
	cin>>n1;
	n=n1;
	m1=log(n1)/log(2);
	m=log(n1)/log(2);
	if(m!=m1) n=pow(2,m+1);
	for(i=0;i<n;i++) 
	{
		a[i].x=a[i].y=0.0;
	}
	cout<<"\n";
	for (i=0;i<n1;i++)
	{
	cout<<"請輸入數據實部:";
    cin>>a[i].x;
	cout<<"請輸入數據虛部:";
	cin>>a[i].y;
	}
	for(z=0;z<=1;z++)
	{
		flag=-1;
		for(m=(log(n)/log(2));m>=1;m--)
		{
			flag=-1;
			for(m=(log(n)/log(2));m>=1;m--)
			{
				le=pow(2,m);
				flag++;
				le1=le/2;
				for(j=0;j<le1;j++)
					ip=i+le1;
				Cequal(t,a[ip]);
				Cequal(t1,a[ip]);
				f=(int)(i*pow(2,flag))%n;
				Wn(w,f);
				Cadd(a[i],t,t1);
				a_ipx=a[ip].x;
				if(z==1)
				{
					w.y*=-1;
				}
				a[ip].x=a[ip].x*w.x-a[ip].y*w.y;
				a[ip].y=a_ipx*w.y+a[ip].y*w.x;
			}
		}
	}
	nu2=n/2;
	nm1=n-2;
	j=0;i=0;
	while(i<=nm1)
	{
		if(i<j)
		{
			Cequal(d,a[j]);
		    Cequal(a[j],a[i]);
			Cequal(a[i],d);
		}
		k=nu2;
		while(k<=j)
		{
			j=j-k;k=k/2;
		}
		j=j+k;
		i=i+1;
		if(z==0)
		{
			cout<<"序列的FFT變換為：";
		}
		else
			cout<<"用FFT計算出的原序列為:"<<'\n';
			for(i=0;i<n;i++)
				if (z==0)
				{
					cout<<a[i].x;
					if(a[i].y>=0)
						cout<<"+j"<<a[i].y<<'\n';
					else
						cout<<"-j"<<a[i].y<<'\n';
					a[i].y=-a[i].y;
				}
				else
				{
					cout<<a[i].x/n;
					a[i].y=-a[i].y/n;
					if(a[i].y>=0)
						cout<<"+j"<<a[i].y<<'\n';
					else
						cout<<"-j"<<fabs(a[i].y)<<'\n';
				}
		}
						cout<<"\n";
}