串行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 dataa[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) arrayPI=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);}}