realfft.cpp

快速傅立叶变换在C＋＋上的实现

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

#include "stdafx.h"
#include "math.h"
#define wei 128
#define pi 3.14159265358979323

void rfft(double x[],int n);

int main(int argc, char* argv[])
{
	double data[wei*2+1];
	int i;

	//2K方波，采样频率为20K
	for(i=1;i<=wei;i++)
	{
		if(i%10<5)
			data[i]=5;
		else
			data[i]=-5;
	}

		printf("\n原始数据:\n");
	for(i=1;i<=wei;i++)
	{
		if(i%8==1)
			printf("\n%6d",i);
		printf("%12.6f",data[i]);	
	}

	printf("\n快速傅立叶变换:\n");

	rfft(data,wei);

	double f;
	for(i=1;i<=wei;i+=2)
	{
		f=(i/2)*20000.0/wei;
		printf("\n%8.1f: %12.6f",f,data[i]);
		//printf("\n%8.1f: %12.6f  ||%12.6f",f,data[i],data[wei*2-i+2]);	
	}

	return 0;
}

///////////////////////////
void rfft(double x[],int n)
{
	int i,j,k,m,i1,i2,i3,i4,n1,n2,n4;
	double a,e,cc,ss,xt,t1,t2;
	for(j=1,i=1;i<16;i++)
	{
		m=i;
		j=2*j;
		if(j==n)break;
	}

	n1=n-1;

	for(j=0,i=0;i<n1;i++)//倒序
	{
		if(i<j)
		{
			xt=x[j];
			x[j]=x[i];
			x[i]=xt;
		}

		k=n/2;

		while(k<(j+1))
		{
			j-=k;
			k>>=1;
		}
		j+=k;
	}

	for(i=0;i<n;i+=2)
	{
		xt=x[i];
		x[i]=xt+x[i+1];
		x[i+1]=xt-x[i+1];
	}

	n2=1;
	for(k=2;k<=m;k++)
	{
		n4=n2;
		n2=2*n4;
		n1=2*n2;
		e=6.28318530718/n1;

		for(i=0;i<n;i+=n1)
		{
			xt=x[i];
			x[i]=xt+x[i+n2];
			x[i+n2]=xt-x[i+n2];
			x[i+n2+n4]=-x[i+n2+n4];
			a=e;
			for(j=1;j<=(n4-1);j++)
			{
				i1=i+j;
				i2=i-j+n2;
				i3=i+j+n2;
				i4=i-j+n1;
				cc=cos(a);
				ss=sin(a);
				a=a+e;
				t1=cc*x[i3]+ss*x[i4];
				t2=ss*x[i3]-cc*x[i4];
				x[i4]=x[i2]-t2;
				x[i3]=-x[i2]-t2;
				x[i2]=x[i1]-t1;
				x[i1]=x[i]+t1;
			}
		}
	}

}