Question

想用C语言实现一个1024点的FFT,找到的基2-FFT的程序能实现128以内的FFT,运行结果和matlab的fft()是一样的

typedef struct{
double r;
double i;
}fcomplex;

//检查a是否为2的整数次方数
#define NOT2POW(a) (((a)-1)&(a)||(a)<=0)
//pi
#define PAI 3.14159265358979323846

int fft(fcomplex* a, unsigned int len){
unsigned int ex=0,n=len,N=len;//ex表示共有ex级蝶形运算
unsigned int nv2,nm1,k;//nv2=n/2,nm1=n-1
unsigned int le,lei,ip;//在一级蝶形运算中，le表示同种蝶形结间的距离，lei表示一个蝶形两节点的距离，
//一个蝶形的两个蝶形下标为[i]和[ip]且ip=i+lei
unsigned int i,j,m;//m:第m级蝶形运算;j:W(r,N)中的r;i:蝶形结的下标[i]
fcomplex u,w,t,u1;
double tmp;

if(NOT2POW(n)) return NULL; //如果失败，返回空指针
for(;!(N&1);N>>=1) ex++; //len应该等于2的ex次方

//用雷德算法实现倒位序
nv2=n>>1;
nm1=n-1;
j=0;
i=0;
for(i;i<nm1;i++)
{
if(i<j)
{
t.r=a[j].r;
t.i=a[j].i;

a[j].r=a[i].r;
a[j].i=a[i].i;

a[i].r=t.r;
a[i].i=t.i;
}
k=nv2;
while(k<=j)
{
j-=k;
k>>=1;
}
j+=k;
}

//ex级递推计算
le=1;
for(m=1;m<=ex;m++)
{
lei=le;
le<<=1;
u.r=1; u.i=0;
tmp=PAI/lei;
w.r=cos(tmp);w.i=-sin(tmp);
for(j=0;j<lei;j++)
{
for(i=j;i<n;i+=le)//
{
ip=i+lei;
t.r=a[ip].r*u.r-a[ip].i*u.i;//一次复乘
t.i=a[ip].r*u.i+a[ip].i*u.r;

a[ip].r=a[i].r-t.r;//先求后项的原因，a[i]的值不会被改变
a[ip].i=a[i].i-t.i;

a[i].r=a[i].r+t.r;
a[i].i=a[i].i+t.i;
}
u1=u;
u.r=u1.r*w.r-u1.i*w.i;
u.i=u1.r*w.i+u1.i*w.r;
}
}
return 1;
}
//测试主程序
int main()
{
int i, j,k,DATA_LEN,flag;
fcomplex *s；
printf("基2_FFT测试\n输入生成序列长度:");
scanf("%d",&DATA_LEN);
s=(fcomplex*)malloc(DATA_LEN*sizeof(fcomplex));
double a[8]={0.3535,0.3535,0.6464,1.0607,0.3535,-1.0607,-1.3535,-0.3535};
for(i=0;i<8;i++)
{
s[i].r=a[i];
}
for(i=8;i<DATA_LEN;i++)
{
s[i].r=0;
}

for(j=0;j<DATA_LEN;j++)
{
s[j].i=0;
}
flag=fft(s,DATA_LEN);
if(!flag){
printf("序列长度不为2的整数次方!\n");
return 0;
}
printf("\n实部\t\t虚部\n");
for(k=0;k<DATA_LEN;k++)
printf("%lf\t%lf\n",s[k].r,s[k].i);
free(s);
return 0;
}
可是512点和1024点就出了问题，哪位高手可以指点一下！！非常感谢！

Answer 1

我也写了个fft程序，256点时计算还正确，但是512和1024点结果就错了，希望楼主答案的话能告知，