fft_test.cpp

快速傅立叶变换的源代码

// FFT_Test.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#include <stdio.h>
#include <process.h>
#include <math.h>
struct compx {float real,image;};
struct xia{char name;int age;};
struct compx EE(struct compx,struct compx);
struct compx xx[257];
void FFT(struct compx *,int);
void DIT_FFT(int log2n,float *xr,float *xi,int ntype);
#define pi 3.14159265
FILE *stream;
struct compx EE(struct compx b1,struct compx b2)
{
	struct compx b3;
	b3.real=b1.real*b2.real-b1.image*b2.image;
	b3.image=b1.real*b2.image+b1.image*b2.real;
	return(b3);

}
void main(void)
{
	int i;
	char   s[] = "***************";
	char   c = '\n';

	double xlable[257];
	double avera=0;
	double cep[256];
	double t,r;
	float xxr[256],xxi[256];
	for(int j=0;j<257;j++)
	{
		xlable[j]=avera;
		avera+=pi/20;
	}
	
	stream = fopen( "fftx.out", "w" );
	for(i=0;i<257;i++)
	{
		if(i<15)
		xx[i].real=5.0;
		//xx[i].real=0;
		xx[i].image=0;
		//fprintf(stream,"%f",xx[i].real);
		//fprintf( stream, "%s%c", s, c );
		//fprintf(stream,"%f",xx[i].real);
	}
	
	//fprintf(stream,"%-10d\t",100);
	FFT(xx,256);
	for(i=0;i<256;i++)
	{
		t=pow(xx[i].real,2);
		r=pow(xx[i].image,2);
		cep[i]=pow(t+r,0.5);
	}
		
	for(i=0;i<256;i++)
	{
		//fprintf(stream,"%-10f\t",sin(xlable[i])+cos(10*xlable[i]));
		fprintf(stream,"%-10f\t",xx[i].real);
		fprintf(stream,"%-10f\t",xx[i].image);
		fprintf(stream,"%-10f\n",cep[i]);
	}
	/*for(i=0;i<256;i++)
	{
		xxr[i]=xx[i].real;
		xxi[i]=xx[i].image;
	}
		DIT_FFT(256,xxr,xxi,-1);
	for(i=0;i<256;i++)
	{
		
		fprintf(stream,"%-10f\t",xxr[i]);
		fprintf(stream,"%-10f\n",xxi[i]);
	
	}*/
	fclose( stream );
   system( "type fftx.out" );
	
}

/* 输入：xin（实部，虚部），输出：xout （实部，虚部），N：FFT点数*/
void FFT(struct compx *xin,int N)
{
	int f,m,nv2,nm1,i,j=1,k,l;
	struct compx v,w,t;
	nv2=N/2;
	f=N;
	for(m=1;(f=f/2)!=1;m++);
	nm1=N-1;
	//变址运算
	for(i=1;i<=nm1;i++)
	{
		if(i<j)
		{
			t.real=xin[j].real;
			t.image=xin[j].image;
			xin[j].real=xin[i].real;
			xin[j].image=xin[i].image;
			xin[i].real=t.real;
			xin[i].image=t.image;
		}
		k=nv2;
		while(k<j)
		{
			j=j-k;
			k=k/2;
		}
		j=j+k;
	}
/*	//*fft*/
	{
		int le,lei,ip;
	
		for(l=1;l<=m;l++)
		{
			le=pow(2,l);
			lei=le/2;
			
			v.real=1.0;
			v.image=0.0;
			w.real=cos(pi/lei);
			w.image=-sin(pi/lei);
			for(j=1;j<=lei;j++)
			{
				for(i=j;i<=N;i=i+le)
				{
					ip=i+lei;
					t=EE(xin[ip],v);
					xin[ip].real=xin[i].real-t.real;
					xin[ip].image=xin[i].image-t.image;
					xin[i].real=xin[i].real+t.real;
					xin[i].image=xin[i].image+t.image;
				}
				v=EE(v,w);
			}
		}
	}
/**/
}
/***************************************************/

/***************************************************/
void DIT_FFT(int log2n,float *xr,float *xi,int ntype)
{
	int i,j,k,n,nv2,nm1,l,le,le1,ip,sign=-1;
	float tr,ti,ur,ui,wr,wi,ain;
	if(ntype<0)sign=1;
	n=1<<log2n;
	nv2=n>>1;
	nm1=n-1;
	j=0;
	for(i=0;i<nm1;i++)
	{
		if(i<j)
		{
			tr=xr[j];
			ti=xi[j];
			xr[j]=xr[i];
			xi[j]=xi[i];
			xr[i]=tr;
			xi[i]=ti;

		}
		k=nv2;
		while(k<=j)
		{
			j-=k;
			k=k>>1;
		}
		j+=k;
	}
	//pi=4.0*atan(1.0);
	for(l=1;l<log2n;l++)
	{
		le=1<<1;
		le1=le/2;
		ur=1.;
		ui=0.;
		wr=cos(pi/le1);
		wi=sign*sin(pi/le1);
		for(j=0;j<le1;j++)
		{
			for(i=j;i<n;i+=le)
			{
				ip=i+le1;
				tr=xr[ip]*ur-xi[ip]*ui;
				ti=xr[ip]*ui+xi[ip]*ur;
				xr[ip]=xr[i]-tr;
				xi[ip]=xi[i]-ti;
				xr[i]=xr[i]+tr;
				xi[i]=xi[i]+ti;
			}
			tr=ur*wr-ui*wi;
			ti=ur*wi+ui*wr;
			ur=tr;
			ui=ti;
		}
	}
	if(ntype>0)return;
	ain=1./n;
	for(i=0;i<n;i++)
	{
		xr[i]=xr[i]*ain;
		xi[i]=xi[i]*ain;
	}
	return;
}