fft-by skyilne.cpp

这是一个我自己写的FFT的c程序

#include <math.h>
#include <stdio.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define NULL 0
#define TYPE struct compx
#define LEN sizeof (struct compx)

struct compx 
{ double real;
  double imag;
  struct compx *next;

} 
s[16]=
{
	{0,0},
	{1,1},
	{2,2},
	{3,3},
	{4,4},
	{5,5},
	{6,6},
	{7,7},
	{8,8},
	{9,9},
	{10,10},
	{11,11},
	{12,12},
	{13,13},
	{14,14},
	{15,15},
};




float  result[257];
int   N=0,L=0,M=0;         //N采样点数   L第几级    M共多少级 N=M^2
const float PI=3.141592;

struct compx mult(struct compx b1,struct compx b2)
{
	struct compx b3;
	b3.real=b1.real*b2.real-b1.imag*b2.imag;
	b3.imag=b1.real*b2.imag+b1.imag*b2.real;
	return(b3);
}

struct compx *creat(int n)
{
struct compx *head,*pf,*pb;
int i;
for(i=0;i<n;i++)
{ 
pb=(TYPE*) malloc(LEN);
printf("input real and imag\n");
scanf("%d%d",&pb->real,&pb->imag);
if(i==0)
pf=head=pb;
else pf->next=pb;
pb->next=NULL;
pf=pb;
}
return(head);
}


int reverse(int x,int M); 
void FFT(int N);
void requeue(int M,int N);
int judge_M(int N);//求级数M  2的多少次方 
void cul(int M,int N);

void main()
{

int v=6,i=0,num=0;
float a,b;









printf("input num (  0<num<=4096  )\n");
scanf("%d",&num);
printf("num= %d\n",num);



if((num<=4096)&&(num>2048))
  N=4096;
else if((num<=2048)&&(num>1024))
  N=2048;
else if((num<=1024)&&(num>512))
  N=1024;
else if((num<=512)&&(num>256))
  N=512;
else if((num<=256)&&(num>128))
  N=256;
else if((num<=128)&&(num>64))
  N=128;
else if((num<=64)&&(num>32))
  N=64;
else if((num<=32)&&(num>16))
  N=32;
else if((num<=16)&&(num>8))
  N=16;
else if((num<=8)&&(num>4))
  N=8;
else if((num<=4)&&(num>2))
  N=4;
else if((num<=2)&&(num>0))
  N=2;




for(i=0;i<N;i++)
scanf("%f%f",&a,&b);
printf("%f\n%f\n",a,b);


FFT(16);

}

void FFT(int N)
{
   int i=0;

   M=judge_M(N);//求级数M  最大2的多少次方   最多4096个点 M=12
   requeue(M,N);//重新排队
   cul(M,N);//开始计算


   //输出结果
  printf("N=%d\n",N); 


}


void cul(int M,int N)  //开始计算
{
    struct compx temp,temp1,w;//compx.w是旋转因子
	int i=0,j=0,k=0,p,q,a,b;  //i  目前在第几级i=L    j 这一级中的第几个旋转因子  k 这个旋转因子是第几次使用  p是每级中旋转因子的个数  q是共使用几次
	for(i=1;i<=M;i++)//按L分  第几级运算  i=L   i=1~M
		{ 
		   
		   p=pow(2,(i-1));//p是每级旋转因子的个数
		   for(j=0;j<p;j++)   //旋转因子的个数分  第j个旋转因子
				{
			       w.real=cos(-(2*PI*k)/pow(2,(i+1)));//求旋转因子
		           w.imag=sin(-(2*PI*k)/pow(2,(i+1)));	   
			       q=pow(2,(M-i));//这个旋转因子的使用次数	
			       for(k=0;k<q;k++)//按每个旋转因子使用的次数分   一次蝶形运算并保存
						{
							a=(N-(pow(2,i-1)-j))-k*pow(2,i);//下面的
							b=(N-(pow(2,i-1)-j))-k*pow(2,i)-pow(2,i-1);//上面的
					        temp=mult(s[a],w);//s[(N-(pow(2,i-1)-j))-k*pow(2,i)]是 应该乘以旋转因子的那值
                            temp1=s[b];
							s[b].real=temp1.real+temp.real;
                            s[b].imag=temp1.imag+temp.imag;
                            s[a].real=temp.real-temp1.real;
                            s[a].imag=temp.imag-temp.imag;
							
						}
				}

		}

}

void requeue(int M, int N)  //重新排队程序
{
	int i=0,rev;
	struct compx temp;
	for(i=0;i<N;i++)
		{
		 rev=reverse(i,M);  //	求倒序的号
         if (rev<i)//防止 一对数据被交换两次
			 goto loop1;
		 temp=s[i];
		 s[i]=s[rev];
         s[rev]=temp;
loop1:    rev++;
		}

}



reverse(int x,int L)     //求倒序后的序列号  返回值是倒序的值   L是二进制的位数
{
	int i=0,temp,rev=0;
	temp=x;

	for(i=0;i<L;i++)
		{
		  if((temp&1)==1)
			  rev=rev|1;
		   if(i!=(L-1))
		   { rev<<=1;
			  temp>>=1;
		   }
		  else
		  {  break;}
		  
		}
    return(rev);
 


}

int judge_M(int N)//求级数M  2的多少次方 
{
if((N<=4096)&&(N>2048))
	   M=12;
   else if((N<=2048)&&(N>1024))
	   M=11;
   else if((N<=1024)&&(N>512))
	   M=10;
   else if((N<=512)&&(N>256))
       M=9;
   else if((N<=256)&&(N>128))
       M=8;
   else if((N<=128)&&(N>64))
       M=7;
   else if((N<=64)&&(N>32))
       M=6;
   else if((N<=32)&&(N>16))
       M=5;
   else if((N<=16)&&(N>8))
       M=4;
   else if((N<=8)&&(N>4))
       M=3;
   else if((N<=4)&&(N>2))
       M=2;
   else 
	   M=1;
   return (M);

}