串行fft 源于并行计算——结构。算法。编程中伪码

#include <stdio.h>
#include <math.h>
#include <stdlib.h>

typedef struct {
double real;
double img;
} com;

double PI;

void cp(com f,com *t);
void add(com a,com b,com* r);
void mul(com a,com b,com* r);
void sub(com a,com b,com* r);
int br(int src,int size);
void init(com* w,int size);
void fft(com* a,com* w,int size);
void show(com *a,int size);
int main() {
int size=4;
com a[size];
com ws[size];
for(int i=0;i<size;++i) {//initialize test data
a[i].real=i+1;
a[i].img=0;
}
init(ws,size);
fft(a,ws,size);
show(a,size);
}
void cp(com f,com *t) {
t->real=f.real;
t->img=f.img;
}
void add(com a,com b,com *c)
{
c->real=a.real+b.real;
c->img=a.img+b.img;
}

void mul(com a,com b,com *c)
{
c->real=a.real*b.real-a.img*b.img;
c->img=a.real*b.img+a.img*b.real;
}
void sub(com a,com b,com *c)
{
c->real=a.real-b.real;
c->img=a.img-b.img;
}
int br(int src,int size) {
int tmp=src;
int des=0;
for(int i=size-1;i>=0;--i) {
des=((tmp&1)<<i)|des;
tmp=tmp>>1;
}
return des;
}
void init(com* w,int size) {//initialize the power(w,i) array
PI=atan(1)*4;
for(int i=0;i<size;++i) {
w[i].real=cos(2*PI*i/size);
w[i].img=sin(2*PI*i/size);//?
}
}
void fft(com* a,com* w,int size) {
com c[size];
com up,down;
for(int i=0;i<size;++i) {
cp(a[i],&c[i]);
}
for(int h=log(size)/log(2)-1;h>=0;--h) {
int p=pow(2,h);
int q=size/p;
int z=q/2;
for(int k=0;k<size;++k) {
if(k%p==k%(2*p)) {
add(c[k],c[k+p],&up);
sub(c[k],c[k+p],&down);
cp(up,&c[k]);
int time=z*(k%p);
mul(down,w[time],&c[k+p]);
}
}
}
int h=log(size)/log(2);
for(int i=0;i<size;++i) {
cp(c[i],&a[br(i,h)]);
}
}
void show(com *a,int size) {
for(int i=0;i<size;++i) {
printf("%.4f %.4f \n",a[i].real,a[i].img);
}
}