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

第三版 295页伪码，跟一般的计算方法的顺序有些不同，这个是先计算后改变位置。简单易懂，明了。陈国良果然大家。

#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);
	}
}