fft.cpp

用c语言写的Fft

//128 dot fft
//2004-9-12
//liaoxingen

#include "stdio.h"
#include "math.h"

#define pi  3.1415926

#define N  256
int iTemp;
int level;

void Fft128(double dataR[],double dataI[]);

void main(void)
{
	int i,j;
	double step;
	double DataR[N];
	double DataI[N];
	double Amplitude[N];
	double Phase[N];
	//输入函数10*sin(t)+5cos(3t)
	step=(2.0*pi)/N;
	for(i=0;i<N;i++)
	{
		DataR[i]=10.0*sin(step*i);
		DataI[i]=0.0;
		Amplitude[i]=0.0;
		Phase[i]=0.0;
	}
	printf("\n");
	for(i=0;i<N;i++)//输入数据
	{
		for(j=0;j<1;j++)
		{
			printf("%f     ",DataR[i*1+j]);
		}
		printf("\n");
	}
	iTemp=N;
	level=0;
	while(iTemp>1)
	{
		level++;
		iTemp=iTemp/2;
	}

	//位倒序
	int x;
	int b[8];
	for(i=0;i<N;i++)
	{
		b[7]=b[6]=b[5]=b[4]=b[3]=b[2]=b[1]=b[0]=0;
		b[level-1]=0x01&(i>>0);
		b[level-2]=0x01&(i>>1);
		b[level-3]=0x01&(i>>2);
		b[level-4]=0x01&(i>>3);
		b[level-5]=0x01&(i>>4);
		b[level-6]=0x01&(i>>5);
		b[level-7]=0x01&(i>>6);
		b[level-8]=0x01&(i>>7);
		x=b[7]*128+b[6]*64+b[5]*32+b[4]*16+b[3]*8+b[2]*4+b[1]*2+b[0]*1;
		DataI[x]=DataR[i];
	}
	for(i=0;i<N;i++)
	{
		DataR[i]=DataI[i];
		DataI[i]=0.0;
	}
	Fft128(DataR,DataI);
	for(i=0;i<N;i++)	//计算幅值和相位
	{
		Amplitude[i]=sqrt(DataR[i]*DataR[i]+DataI[i]*DataI[i]);
		Phase[i]=atan(DataR[i]/DataI[i]);
	}
	printf("\n");
	printf("////////////////////////////幅值:");
	printf("\n");
	for(i=0;i<N;i++)
	{
		for(j=0;j<1;j++)
		{
			printf("%f     ",Amplitude[i*1+j]);
		}
		printf("\n");
	}
	printf("////////////////////////////相位:");
	printf("\n");
	for(i=0;i<N;i++)
	{
		for(j=0;j<1;j++)
		{
			printf("%f     ",Phase[i*1+j]);
		}
		printf("\n");
	}
}

//fft模块
void Fft128(double dataR[],double dataI[])
{
	int i,j,L,b,s,k,nk;
	double TR,TI,temp;
	printf("\n");
	printf("calculating...\n");
	for(L=1;L<=level;L++)  //共有L级碟型运算
	{
		//每一级的群数为 b=N/(2^L)
		i=level-L;
		b=1;
		while(i>0)  
		{
			b=b*2;
			i--;
		}
		//每一级碟型因子数 s=2^(L-1)
		i=L-1;
		s=1;
		while(i>0)
		{
			s=s*2;
			i--;
		}
		for(i=1;i<=b;i++)//第L级的i个群
		{
			k=0;
			nk=1;
			for(j=1;j<=s;j++) //第L级的第i个群的碟型运算，共有s个碟型因子
			{
				nk=k*b;//求碟型因子
				TR=dataR[j+(i-1)*s*2]; TI=dataI[j+(i-1)*s*2]; temp=dataR[j+(i-1)*s*2+s];
				dataR[j+(i-1)*s*2]=dataR[j+(i-1)*s*2]+dataR[j+(i-1)*s*2+s]*cos(2.0*pi*nk/N)+dataI[j+(i-1)*s*2+s]*sin(2.0*pi*nk/N);
				dataI[j+(i-1)*s*2]=dataI[j+(i-1)*s*2]-dataR[j+(i-1)*s*2+s]*sin(2.0*pi*nk/N)+dataI[j+(i-1)*s*2+s]*cos(2.0*pi*nk/N);
				dataR[j+(i-1)*s*2+s]=TR-dataR[j+(i-1)*s*2+s]*cos(2*pi*nk/N)-dataI[j+(i-1)*s*2+s]*sin(2.0*pi*nk/N);
				dataI[j+(i-1)*s*2+s]=TI+temp*sin(2.0*pi*nk/N)-dataI[j+(i-1)*s*2+s]*cos(2.0*pi*nk/N);
				k++;     
			}
		}
		
	}
}

/*
for(i=0;i<32;i++) // 如只需要32次以下的谐波进行分析 
{ 
	w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);
	w[i]=w[i]/64；
}
	w[0]=w[0]/2;
*/