fft.cpp

最高1024点的FFT程序

//输入内容可在put_in内设置，最高1024点的fft
#include <conio.h>
#include <math.h>
#include <process.h>
#include <conio.h>
#include <stdio.h>
//#include <complex.h>

#define pi 3.1416

double x[1024],y[1024],am[1024];

void main()
{
	int l,N,i;
	void put_in(int N);
	void re_bit(int N);
	void butterfly(int l);

	printf("Please input the order of fft l(l<=10)=");
	scanf("%d",&l);
	N=int(pow(2,l));

	put_in(N);
	re_bit(N);
	//for(i=0;i<N;i++)
	//{printf("%f_",x[i]);//位反序后的结果
	//}
	butterfly(l);
	

	for(i=0;i<N;i++)
			printf("%f+j%f	",x[i],y[i]);
	printf("\n\n");

		for(i=0;i<N;i++)
		{
			am[i]=sqrt(x[i]*x[i]+y[i]*y[i]);
			printf("%f	",am[i]);
		}
	printf("\n\n	");
	//printf("\n%f",floor(0.1));//测试取整函数
}

void put_in(int N)			//输入信号
{
	int i;
	for(i=0;i<N;i++)
	{
		//x[i]=cos(pi*i*i/N);
		//y[i]=sin(pi*i*i/N);
		x[i]=i;
		y[i]=0;

	}	

}

void re_bit(int N)			//位反序
{
	int i,j,M,s;
	double a;
	M=N/2;
	i=0;

	for(j=1;j<N-1;j++)
	{
		s=M;
		while(i>=s)
		{i=i-s;s=s/2;}
		i=i+s;
		if(j<i)
		{
			a=x[i];x[i]=x[j];x[j]=a;
			a=y[i];y[i]=y[j];y[j]=a;
			
		}
		
	}

		for(i=0;i<N;i++)
		{
			printf("%f\n",x[i]);//输入内容
		}
}

void butterfly(int l)						//蝶形运算
{
	int i,j,r,m1,m2,m3,m4,N,k1,k2;
	double u,v;
	N=int(pow(2,l));
	for(i=1;i<l+1;i++)
	{
		m1=int(pow(2,i-1));
		m2=2*m1;
		m3=N/m2;
		for(j=0;j<m3;j++)
		{
			m4=j*m2;
			for(r=0;r<m1;r++)
			{
				k1=m4+r;
				k2=k1+m1;
				u=x[k2]*cos(2*pi*r/m2)+y[k2]*sin(2*pi*r/m2);
				v=y[k2]*cos(2*pi*r/m2)-x[k2]*sin(2*pi*r/m2);
				x[k2]=x[k1]-u;
				y[k2]=y[k1]-v;
				x[k1]=x[k1]+u;
				y[k1]=y[k1]+v;
			}
		}
	}
}