/*******************************************************************************************************************************************
FFT蝶形算法,程序中包含了原创的补零函数,简单的倒位序函数等
依照按时间抽取对序列进行FFT计算.
********************************************************************************************************************************************/
#include<stdio.h>
#include<math.h>
#define swap(a,b) {a=a+b;b=a-b;a=a-b;}//数值交换
#define MAXlen 256
typedef struct complx
{
double real,imag;
}cmp;
/**************复数乘法***************/
cmp plu_mul(cmp x,cmp y)
{
cmp z;
z.real=x.real*y.real-x.imag*y.imag;
z.imag=x.real*y.imag+y.real*x.imag;
return(z);
}
/*********************************************
****************Z.Y.H原创补零程序****************/
int jud_zero(cmp *a,int len)
{
int i=1,j,len_tmp=len;
while(i<len_tmp)
i*=2;
len_tmp=i;
i=0;
while(a[i].real>0)//可知未赋值的字符数组初始值为负数,由此判断数据个数
i++;
if(len_tmp>i)
for(j=i;j<len_tmp;j++)
{
a[j].real=0.0;
a[j].imag=0.0;
}
return(len_tmp);
}
/***********Z.Y.H原创倒位序算法******************/
void reverse(cmp *a,int len)
{
int tmp,i,i_tmp,j;
printf("Reverse: ");
for(i=0;i<len;i++)
{
tmp=0;
i_tmp=i;
for(j=1;j<len;j*=2)
{
tmp=tmp*2+i_tmp%2;//类似十进制倒序
i_tmp/=2;
}
printf("%d ",tmp);//测试倒序数
swap(a[i_tmp].real,a[i].real);
swap(a[i_tmp].imag,a[i].imag);
}
printf("\n");
}
/************************************************
************倒位序 Rader算法******************/
/*void reverse(cmp *a,int len)
{
int n = len,n1,i,j,k;
cmp tmp;
j=0;
n1 = n-1;
for(i = 0;i < n1;i++)
{
if(i < j)
{
tmp=a[i];
a[i]=a[j];
a[j]=tmp;
}
k=len/2;
while(k<=j)
{
j-=k;
k/=2;
}
j+=k;
}
}*/
void FFT(cmp *a,int len)
{
int k1,N,m,j,k,M,N2,l,p;
cmp tmp,w,tmp2;
double pi=3.141592653589793;
len=jud_zero(a,len);//自动补零
reverse(a,len);//////数据倒序
l=len;
M=0;//计算蝶形总级数
while(l!=1)
{
l/=2;
M++;
}
printf("Level=%d\n",M);//测试级数
for(m=1;m<=M;m++)
{
N=pow(2,m);
N2=pow(2,m-1);//N2即N/2,第m级蝶形距离
//tmp.real=1.0;
//tmp.imag=0.0;
//w.real=cos(2*pi/len);//w为系数商,根据系数的递推公式:W(kl,N)=W((k-1)l,N)*W(l,N);
//w.imag=-sin(2*pi/len);
for(j=0;j<=N2-1;j++)
{
p=pow(2,M-m)*j;
w.real=cos(2*pi*p/len);
w.imag=-sin(2*pi*p/len);
for(k=j;k<len;k+=N)
{
k1=k+N2;
tmp2=plu_mul(a[k1],w);//蝶形运算
a[k1].real=a[k].real-tmp2.real;
a[k1].imag=a[k].imag-tmp2.imag;
a[k].real=a[k].real+tmp2.real;
a[k].imag=a[k].imag+tmp2.imag;
}
//tmp=plu_mul(tmp,w);
}
}
}
void main()
{
cmp x[MAXlen];
int i;
for(i=0;i<16;i++)
{
x[i].real=1.0;
x[i].imag=0.0;
}
FFT(x,8);//FFT(x,21);//FFT(x,7);//可测试补零程序
for(i=0;i<8;i++)
printf("%f %fi\n",x[i].real,x[i].imag);
/*printf("Please enter x(n)'s lenth");
scanf("%d",&len);
for(i=0;i<len;i++)
scanf("%f %f\n",&x[i].real,&x[i].imag);*/
}