fft.txt

这是一个实现快速傅里叶变换的C语言的源程序。

#include<graphics.h>
#include<stdlib.h>
 
#include<stdio.h>
#include<math.h>
#define PI 3.141593
int nu,i,n=1;
float xr[1024],xi[1024],w[1024];
main()
{
	float x;
	int v;
	int graphdrv=VGA;
	int graphmode=VGAHI;
	initgraph(&graphdrv,
&graphmode,"\\tc\\bgi");
	printf("\n\n n=2**r\n\n");
	printf("r=5,6,7,8,9 or 10 \n\n");
	printf("INPUT ONE OF 
r's VALUES: \n\n r=");
	scanf("%d",&nu);
	for(i=0;i<nu;i++)
	n=2*n;
	printf("Waiting,Data are processed...\n");
	x=4*PI/n;
	for(i=0;i<n;i++)
	{
		xr[i]=cos(i*x)+0.5*cos(2*i*x)+0.8*cos(5*i*x);
		xr[i]=xr[i]/n;
		xi[i]=0;
w[i]=cos(i*x)+0.5*cos(2*i*x)+0.8*cos(5*i*x);
	}
	cleardevice();
	settextstyle(1,0,0);
	outtextxy(100,3,"Here os the curve f(t)-t:");
	showgraph(220,50,n);
	v=n/2;
	fft(v);
	settextstyle(1,0,0);
	outtextxy(100,3,"Here is the curve F(f)-f:");
	showgraph(220,300,n);
}
	/*用于画图的函数*/
	showgraph(y0,r,m)
	int y0,r,m;
	{
		float t,l;
		int y1,i,x1;
		line(0,20,0,420);
		for(i=0;i<41;i++)
		line(0,10*i+20,2,10*i+20);
		line(0,y0,639,y0);
		for(i=0;i<65;i++)
		line(10*i,y0,10*i,218);
		moveto(0,y0);
		for(i=0;i<m;i++)
		{
			l=640.0/m*i;
			x1=(int)l;
			t=y0-w[i]*r;
			y1=(int)t;
			lineto(x1,y1);	
		}
		settextstyle(3,0,0);
		outtextxy(100,400,"Press any key to continue!");
		getch();
		cleardevice();		
	}
	/*实现傅里叶变换的函数*/
	fft(n2)
	int n2;
	{
		int nu1,kn,l,i,g,k,p;
		float arg,c,s;
		int j,h=1;
		float tr,ti;
		nu1=nu-1;
		k=0;
		for(l=1;l<=nu;l++)
		{
			for(g=1;g<=nu1;g++)	
				h=2*h;
			do
			{
				for(i=1;i<=n2;i++)
				{
					j=k/h;
					p=ibitr(j,nu);
					arg=6.283185*p/n;
					c=cos(arg);
					s=sin(arg);
					kn=k+n2;
					tr=xr[kn]*c+xi[kn]*s;
					ti=xi[kn]*c-xr[kn]*s;
					xr[kn]=xr[k]-tr;
					xi[kn]=xi[k]-ti;
					xr[k]=xr[k]+tr;
					xi[k]=xi[k]+ti;
					k=k+1;	
				}
				k=k+n2;
			}
			while((k+n2)<n);
			k=0;
			nu1=nu1-1;
			h=1;
			n2=n2/2;
		}
		for(k=0;k<n;k++)
		{
			j=k;
			i=ibitr(j,nu);
			if(i<=k) continue;
			tr=xr[k];
			xr[i]=tr;
			ti=xi[k];
			xi[k]=xi[i];
			xi[i]=ti;
			
		}
		for(i=0;i<n/2;i++)
		{
			w[i]=xr[i]*xr[i]+xi[i]*xi[i];
			w[n-i-1]=0;
		}
	}
	/*实现序数恢复的函数*/
	int ibitr(j,nu)
	int j,nu;
	{
		int m,j1,j2,ibr;
		ibr=0;
		j1=j;
		for(m=1;m<=nu;m++)
		{
			j2=j1/2;
			ibr=2*ibr+(j1-j2*2);
			j1=j2;
		}
		return(ibr);
	}