afftfilter.c

自己编的数值计算方法的c语言实现源程序。

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

void ChangeOrder(double xr[],double xi[],int N)
	{
		int LH,N1,I,J,K;
		double T;
		LH=N/2; J=LH; N1=N-2;
		for(I=1;I<=N1;I++)
		{
			if(I<J)
			{
				T=xr[I]; xr[I]=xr[J]; xr[J]=T;
				T=xi[I]; xi[I]=xi[J]; xi[J]=T;
			}
			K=LH;
			while(J>=K)
			{
				J=J-K;
				K=(int)(K/2+0.5);
			}
			J=J+K;
	/*	printf(" I=%d,J=%d,K=%d\n", I,J,K);  */
		}
	/*	getch(); */
	}

void Fft(double xr[], double xi[], int N, int M)
{
	int L, B, J, P, k,i;
	double rpartKB, ipartKB;
	double rcf[256],icf[256];
	double PI2=8.0*atan(1.0);
	for(i=0;i<N;i++)
	{
		rcf[i]=cos(i*PI2/N);
		icf[i]=sin(i*PI2/N);
	}
	
	ChangeOrder(xr,xi,N);
	
	for(L=1;L<=M;L++)
	{
		B=(int)(pow(2,(L-1))+0.5);
		for(J=0;J<=B-1;J++)
		{
			P=J*((int)(pow(2,(M-L))+0.5));
			for(k=J;k<=N-1;k+=(int)(pow(2,L)+0.5))
			{
				rpartKB=xr[k+B]*rcf[P]-xi[k+B]*icf[P];
				ipartKB=xr[k+B]*rcf[P]+xr[k+B]*icf[P];
				xr[k+B]=xr[k]-rpartKB;
				xi[k+B]=xi[k]-ipartKB;
				xr[k]=xr[k]+rpartKB;
				xi[k]=xi[k]+ipartKB;
			}
		}
	}
}


void Ifft(double xr[], double xi[], int N, int M)
{
	int L, B, J, P, k,i;
	double rpartKB, ipartKB;
	double rcf[256],icf[256];

	double PI2=8.0*atan(1.0);
	for(i=0;i<N;i++)
	{
		rcf[i]=cos(i*PI2/N);
		icf[i]=sin(i*PI2/N);
	}

	ChangeOrder(xr,xi,N);

	for(L=1;L<=M;L++)
	{
		B=(int)(pow(2,(L-1))+0.5);
		for(J=0;J<=B-1;J++)
		{
			P=J*((int)(pow(2,(M-L))+0.5));
			for(k=J;k<=N-1;k+=(int)(pow(2,L)+0.5))
			{
				rpartKB=xr[k+B]*rcf[P]+xi[k+B]*icf[P];
				ipartKB=xr[k+B]*rcf[P]-xr[k+B]*icf[P];
				xr[k+B]=xr[k]-rpartKB;
				xi[k+B]=xi[k]-ipartKB;
				xr[k]=xr[k]+rpartKB;
				xi[k]=xi[k]+ipartKB;
			}
		}
		
	}
	
	for(i=0;i<N;i++)
	{
		xr[i]=xr[i]/N;
		xi[i]=xi[i]/N;
	}
}	
	

	
	
	void main(void)
	{
		double y[256],yi[256];
		double t=0;
		double pace=atan(1.0)/32;
		int i;

		FILE *out;
		for(i=0;i<256;i++)
		{
			t+=pace;
			y[i]=exp(-(cos(t)*cos(t))*(sin(2*t)+2*cos(4*t)+0.4*sin(t)*sin(50*t)));
			yi[i]=0;
			printf(" %d",i);
		}
          getch();
		
		Fft(y,yi,256,8);
		for(i=6;i<250;i++)
		{	y[i]=0; }
		Ifft(y,yi,256,8);
		
		
		for(i=0;i<256;i++)
		{
			printf(" %f", y[i]);
		}
		getch();  
	}