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

void   fft(double a[],double b[],int n,int flag,double c[],double d[]);

 void main()
 { 

	FILE *fp_result;
	fp_result=fopen("fuliye_result.txt","w");

	if((fp_result = fopen("fuliye_result.txt","w")) == NULL)
	{
	        printf("fp_result 文件打不開\n");
	        return;
	}

	int N,FLAG;
	double T;

	printf("Please input the 周期：");
	scanf("%lf",&T);
	printf("Please input the 采樣點數N和旗幟FLAG： ");
	printf("FLAG 表示是進行正變換還是逆變換\n.");
	scanf("%d%d",&N,&FLAG);
	double h; //h為時間域的采樣間隔。
	h=T/N;
	printf("h=%lf\n",h);
	cout<<endl;
	
	double  *pr=new double [N];
	double  *pi=new double [N];
	double  *fr=new double [N];
	double  *fi=new double [N];

    for (int i=0; i<N; i++)
    { 
		pr[i]=exp(-h*(i+0.5));  
	    pi[i]=0.0;
		fr[i]=0;
		fi[i]=0;
  	}
	// 以上循環是計算所需要變換的離散序列,可以是自輸入的離散數據。
   

	/*
	for (i=0; i<N/4; i++)
      { for (int j=0; j<=3; j++)
          printf("pr[%d]=%8.5f  ",4*i+j,pr[4*i+j]);
          printf("pi[%d]=%8.5f  ",4*i+j,pi[4*i+j]);
        printf("\n");
      }
	*/
	 //以上是打印原始離散序列。

    fft(pr,pi,N,0,fr,fi);
	for(i=0;i<N;i++)
	{
		printf("pr[%2d]=%8.5f ; pi[%2d]=%8.5f ;",i,pr[i],i,pi[i]);
		printf("fr[%2d]=%8.5f ; fi[%2d]=%8.5f \n",i,fr[i],i,fi[i]);
	}
	//調用子程序，求取離散序列的傅立葉變換

    cout<<"調用結束，已經進行傅立葉變換，開始返回"<<endl;
    cout<<"fr為返回的實部，fi為返回的虛步"<<endl;
	cout<<"pr為原始序列實部，pi為原始序列虛部"<<endl;

	double *frequency=new double [N];
    for (i=0;i<N;i++)
	{
		frequency[i]=0;
		//printf(" frequency[%d]=%8.5f   \n",i,frequency[i]);
    }
  
	for (i=0;i<N;i++)
	{
		frequency[i]=i*1/(N*h);
		//printf(" frequency[%d]=%8.5f   \n",i,frequency[i]);
    }
	//計算離散變換后,每個點所對應的頻率。采樣頻率間隔=1/（N*采樣間隔）;本例子中為：1/(N*h);
	//值得注意的是，我們顯示的頻譜的頻率是求得的最大頻率的一半。
    
	cout<<"please input the low and high frequency for filter:\n";
	cout<<"the maximum  frquency ="<<frequency[N-1]<<endl;

	double lowfrequency, highfrequency;

	cin>>lowfrequency>>highfrequency;

    double  *lhFrequency=new double [N];
    for(i=0;i<N;i++	)
	{
		lhFrequency[i]=0;
	}

   	for (i=0;i<N;i++)
	{
		if (frequency[i]>=lowfrequency && frequency[i]<=highfrequency)
		{
			lhFrequency[i]=1;
	     	lhFrequency[N-1-i]=1;
		}
	}
	for(i=0;i<N;i++)
	{
		printf("lhFrequency[%2d]=%8.5f\n",i,lhFrequency[i]);
	}
   //以上是設計帶通濾波器，注意的是：要針對傅立葉變換的對稱性，設計帶通。

    double  *Amplify=new double [N];
	double  *afterReal=new double [N];
	double  *afterImagine=new double [N];
	double  *afterAmplify=new double [N];
	double  *afterPhase=new double [N];
	for(i=0;i<N;i++)
	{
		Amplify[i]=0;
		afterReal[i]=0;
		afterImagine[i]=0;
		afterAmplify[i]=0;
		afterPhase[i]=0;
	}

    for (i=0;i<N;i++)
	{
		Amplify[i]=sqrt(fr[i]*fr[i]+fi[i]*fi[i]);
		afterReal[i]=lhFrequency[i]*fr[i];
		afterImagine[i]=lhFrequency[i]*fi[i];
		afterAmplify[i]=sqrt(afterReal[i]*afterReal[i]+afterImagine[i]*afterImagine[i] );

		if( fabs(fr[i]) <= 0.000001*fabs(fi[i]) )
		{
			if( fr[i]*fi[i]>0) afterPhase[i]=90;
			else afterPhase[i]=-90;
		}
		else 
			afterPhase[i]=atan( fi[i]/fr[i] )*360.0/6.283185306;
	}
	//以上是進行帶通濾波處理；

	double * timereal=new double [N];
	double * timeimag=new double [N];
	for(i=0;i<N;i++)
	{
		timereal[i]=0;
		timeimag[i]=0;
	}

    fft(afterReal,afterImagine,N,1,timereal,timeimag); 
	//對濾波后頻率域的信號，在反變換到時間域。故調用了傅立葉逆變換程序。
	
   	fprintf(fp_result,"%8s  %8s  %8s  %8s  %8s  %8s  %8s  %8s  %8s  %8s  %8s  \n","序列號碼",
		"實際序列","變換振幅","變換相位","變換實部","變換虛部","濾波振幅","濾波相位",
		"序列實部","序列虛部","對應頻率");

    for(i=0;i<N;i++)
	{
		fprintf(fp_result,"%8d  %8.5f  %8.5f  %8.5f  %8.5f  %8.5f  %8.5f  %8.5f  %8.5f  %8.5f  %8.5f\n",
			i,pr[i],Amplify[i],afterPhase[i],fr[i],fi[i],afterAmplify[i],afterPhase[i],timereal[i],
			timeimag[i],frequency[i]); 
	}
   
   
	
	 delete   [] pr; delete  [] pi;delete [] fr;delete [] fi; 
     delete [] frequency;  delete []  lhFrequency; delete [] Amplify;
     delete [] afterReal;
	 delete [] afterImagine;
	 delete [] afterAmplify;
	 delete [] afterPhase;
     delete [] timereal;
	 delete [] timeimag;
	 
	  fclose (fp_result);
}

//////////////////////////////////////////////////////
// fft是個傅立葉變換的子程序，當flag=0時，做正變換；//
//                             falg!=0時，做逆變換。//
//////////////////////////////////////////////////////
void   fft(double a[],double b[],int n,int flag,double c[],double d[])
{
   
    double 	p=6.283185306/(1.0*n);
	int jk=0;
	if(flag==0)
	for(int j=0;j<n;j++)
	{
		for (int k=0;k<n;k++)
		{
		    jk=-j*k;
			c[j]=c[j]+a[k]*cos(p*jk)-b[k]*sin(p*jk);
			d[j]=d[j]+a[k]*sin(p*jk)+b[k]*cos(p*jk);
		}
	}

	jk=0;
	if(flag!=0)
	for(int j=0;j<n;j++)
	{
		for (int k=0;k<n;k++)
		{
		    jk=j*k;
		    c[j]=c[j]+a[k]*cos(p*jk)-b[k]*sin(p*jk);
			d[j]=d[j]+a[k]*sin(p*jk)+b[k]*cos(p*jk);
		}
		c[j]=c[j]/n;
		d[j]=d[j]/n;
	}
}