fuliye.cpp

能够数据的正反傅立叶变换及选择频率进行滤波。

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