fft.cpp

fft算法C程序实现,可以用来快速傅里叶变换

#include "math.h"
#include "stdio.h"
#include "malloc.h"
#include "time.h"
#define PI 3.1415926
int N,jie=0;						//全局变量N及N的阶数

struct complex                      //复数结构体
{
	double real;
	double imreal;
};

complex unit(int k)					//计算旋转因子函数
{
	complex a;
	a.real=cos(2*PI*k/N);
	a.imreal=-sin(2*PI*k/N);
	return a;
}

complex add(complex a,complex b)	//复数加法
{
	complex c;
	c.real=a.real+b.real;
	c.imreal=a.imreal+b.imreal;
	return c;
}

complex sub(complex a,complex b)	//复数减法
{
	complex c;
	c.real=a.real-b.real;
	c.imreal=a.imreal-b.imreal;
	return c;
}

complex mul(complex a,complex b)	//复数乘法
{
	complex c;
	c.real=a.real*b.real-a.imreal*b.imreal;
	c.imreal=a.real*b.imreal+a.imreal*b.real;
	return c;
}

int bit_reverse(int k,int n)		//比特倒序
{
	int a,b=0;
	for(int i=0;i<n;i++)
	{
		a=k%2;
		b=b+a*int(pow(2,n-1-i));
		k=k/2;
	}
	return b;
}

void dit(float *y,complex *X,complex *u)   //时间抽取fft算法
{
	complex *g=(complex *)malloc(sizeof(complex)*N);
	int temp,i,j,k;
	complex a;
	for(i=0;i<N;i++)
	{
		g[i].real=y[i];
		g[i].imreal=0;
	}
	for(i=0;i<jie;i++)
	{
		temp=int(pow(2,i));
		for(j=0;j<N/2/temp;j++)
			for(k=0;k<temp;k++)
			{
				a=mul(u[N/temp/2*k],g[k+2*temp*j+temp]);
				X[k+2*temp*j]=add(g[k+2*temp*j],a);
				X[k+2*temp*j+temp]=sub(g[k+2*temp*j],a);
			}
		for(k=0;k<N;k++)
			g[k]=X[k];
	}
}

void dif(float *y,complex *X,complex *u)   //频率抽取fft算法
{
	complex *g=(complex *)malloc(sizeof(complex)*N);
	complex a;
	int temp,i,j,k,n;
	for(i=0;i<N;i++)
	{
		g[i].real=y[i];
		g[i].imreal=0;
	}
	for(i=0;i<jie;i++)
	{
		temp=int(pow(2,i));
		for(j=0;j<N/2/temp;j++)
			for(k=0;k<temp;k++)
			{
				a=sub(g[k+2*temp*j],g[k+2*temp*j+temp]);
				n=bit_reverse(j,jie-1-i);
				X[k+2*temp*j]=add(g[k+2*temp*j],g[k+2*temp*j+temp]);
				X[k+2*temp*j+temp]=mul(a,u[temp*n]);
			}
		for(k=0;k<N;k++)
			g[k]=X[k];
	}
}

void dft(float *y,complex *X,complex *u)   //dft算法
{
	int i,j;
	complex *g=(complex *)malloc(sizeof(complex)*N);
	for(i=0;i<N;i++)
	{
		g[i].real=y[i];
		g[i].imreal=0;
	}
	for(i=0;i<N;i++)
	{
		X[i].real=0;
		X[i].imreal=0;
		for(j=0;j<N;j++)
			X[i]=add(X[i],mul(g[j],unit(j*i)));
	}
}

void main()
{
	float *x,*y;
	int i=1,j,choice;
	complex *X_dit,*X_dif,*X_dft,*u,*u_dft;
	clock_t start_dit,finish_dit,start_dif,finish_dif,start_dft,finish_dft;  //计时变量
	printf("input the number of DFT:");  //输入点数（2的某次方）
	scanf("%d",&N);
	while(i<N)						//算得N是2的多少阶
	{
		i=i*2;
		jie++;
	}
	x=(float *)malloc(sizeof(float)*N);
	for(i=0;i<N;i++)				//赋初值，三角波，周期为5（可改为别的输入）
		x[i]=i%5;
	u=(complex *)malloc(sizeof(complex)*N/2);						
	for(i=0;i<N/2;i++)				//fft旋转因子
		u[i]=unit(i);
	u_dft=(complex *)malloc(sizeof(complex)*N*N);
	for(i=0;i<N;i++)				//dft旋转因子
		for(j=0;j<N;j++)
			u_dft[i*N+j]=unit(j*i);
	y=(float *)malloc(sizeof(float)*N); //将时域抽样点比特转置，使得正序输出
	for(i=0;i<N;i++)
		y[i]=x[bit_reverse(i,jie)];
	X_dit=(complex *)malloc(sizeof(complex)*N);		//DIT基2
	start_dit=clock();
	dit(y,X_dit,u);
	finish_dit=clock();
	X_dif=(complex *)malloc(sizeof(complex)*N);		//DIF基2
	start_dif=clock();
    dif(y,X_dif,u);
	finish_dif=clock();
	X_dft=(complex *)malloc(sizeof(complex)*N);     //DFT
	start_dft=clock();
    dft(x,X_dft,u_dft);
	finish_dft=clock();
	do												//操作选择
	{
		printf("请选择:\n1、输出结果\n2、输出运行时间\n3、退出\n");
		scanf("%d",&choice);
		switch(choice)
		{
		case 1: printf("时间抽取fft算法(基2)\n");
				for(i=0;i<N;i++)
				printf("%8.6f  %8.6fi\n",X_dit[i].real,X_dit[i].imreal);
				printf("频率抽取fft算法(基2)\n");
				for(i=0;i<N;i++)
				printf("%8.6f  %8.6fi\n",X_dif[i].real,X_dif[i].imreal);
				printf("DFT\n");
				for(i=0;i<N;i++)
				printf("%8.6f  %8.6fi\n",X_dft[i].real,X_dft[i].imreal);
				break;
		case 2: printf("时间抽取fft算法(基2):%d\n",finish_dit-start_dit);
				printf("频率抽取fft算法(基2):%d\n",finish_dif-start_dif);
				printf("DFT:%d\n",finish_dft-start_dft);
		default: ;
		}
	}while(choice!=3);
}