c语言搜索算法伪代码,FFT算法伪代码

FFT简介(来源：百度百科)：

快速傅立叶变换(FFT)是离散傅立叶(DFT)的快速算法，它是根据离散傅立叶变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。它对傅立叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。

设x(n)为N项的复数序列，由DFT变换，任一X(m)的计算都需要N次复数乘法和N-1次复数加法，而一次复数乘法等于四次实数乘法和两次实数加法，一次复数加法等于两次实数加法，即使把一次复数乘法和一次复数加法定义成一次“运算”(四次实数乘法和四次实数加法)，那么求出N项复数序列的X(m),即N点DFT变换大约就需要N2次运算。当N=1024点甚至更多的时候，需要N2=1048576次运算，在FFT中，利用WN的周期性和对称性，把一个N项序列(设N=2k,k为正整数)，分为两个N/2项的子序列，每个N/2点DFT变换需要(N/2)2次运算，再用N次运算把两个N/2点的DFT变换组合成一个N点的DFT变换。这样变换以后，总的运算次数就变成N+2(N/2)2=N+N2/2。继续上面的例子，N=1024时，总的运算次数就变成了525312次，节省了大约50%的运算量。而如果我们将这种“一分为二”的思想不断进行下去，直到分成两两一组的DFT运算单元，那么N点的DFT变换就只需要Nlog2N次的运算，N在1024点时，运算量仅有10240次，是先前的直接算法的1%，点数越多，运算量的节约就越大，这就是FFT的优越性。

Iterative_FFT

参数 n ——a的元素个数

a[0,1,…,(n-1)]——系数

局部变量

wd——单位元的n次原根

w——当前点的估值多项式

y——变换结果

y = BitReverse(a)

for j = 1 to logn

d =

wd =

w = 1

for k = 0 to d/2-1

for m = k to n-1 step d

t = w×y[m + d/2]

x = y[m]

y[m] = x + t

y[m + d/2] = x – t

endfor

w = w×wd

endfor

endfor

return y

C语言实现：

void Iterative_FFT(double *ir,double *ii,double *or,double

*oi,int len)

{

int i,j,n,d,k,m;

double

t;

double wdr,wdi,wr,wi;

double wtr,wti;

double tr,ti;

double xr,xi;

n=pow2(len);

for(i=0;i

{

j = BitReverse(i,len);

or[j] = ir[i];

oi[j] = ii[i];

}

for(j=1;j<=len;j++)

{

d = pow2(j);

t = 2*pi/d;

wdr = cos(t); wdi = sin(t);

wr = 1; wi = 0;

for (k=0;k<=d/2-1;k++)

{

for (m=k;m<=n-1;m+=d)

{

tr = wr*or[m + d/2] - wi*oi[m + d/2];

ti = wr*oi[m + d/2] + wi*or[m + d/2];

xr = or[m]; xi = oi[m];

or[m] = xr + tr; oi[m] = xi + ti;

or[m+d/2] = xr - tr; oi[m+d/2] = xi -

ti;

}

wtr = wr*wdr - wi*wdi;

wti = wr*wdi + wi*wdr;

wr = wtr;

wi = wti;

}

}

}

FFT的串行递归实现伪代码：

Recursive_FFT

参数 n ——a的元素个数

a[0,1,…,(n-1)]——系数

局部变量

wn——单位元的n次原根

w——当前点的估值多项式

a0——偶数项系数

a1——奇数项系数

y——变换结果

y0——a0的FFT结果

y1——a1的FFT结果

if n = =1 then return a

else

wn =

w = 1

a0 = (a[0],a[2],…a[n-2])

a1 = (a[2],a[3],…a[n-1])

y0 = Recursive_FFT(a0,n/2)

y1 = Recursive_FFT(a1,n/2)

for k = 0 to n/2-1

y[k] = y0[k] + w×y1[k]

y[k +

n/2] = y0[k]－w×y1[k]

w =

w×wn

endfor

return y

endif

C语言实现：

void Recursive_FFT(double *ir,double *ii,double *or,double

*oi,int len)

{

int i,k;

int n;

double wnr,wni;

double wr,wi;

double wtr,wti;

double temp0_ir[MAX_N], temp0_ii[MAX_N], temp1_ir[MAX_N],

temp1_ii[MAX_N];

double temp0_or[MAX_N], temp0_oi[MAX_N], temp1_or[MAX_N],

temp1_oi[MAX_N];

memset(temp0_ir, 0 ,sizeof(temp0_ir));

memset(temp0_ii, 0 ,sizeof(temp0_ii));

memset(temp1_ir, 0 ,sizeof(temp1_ir));

memset(temp1_ii, 0 ,sizeof(temp1_ii));

memset(temp0_or, 0 ,sizeof(temp0_oi));

memset(temp0_or, 0 ,sizeof(temp0_oi));

n=pow2(len);

if(len == 0)

{

or = ir;

oi = ii;

}

else

{

wnr =

cos(2*pi/n);

wni = sin(2*pi/n);

//printf("wnr:%.2f wni:%.2f\n",wnr,wni);

wr = 1;

wi = 0;

for(i=0; i

{

temp0_ir[i/2] = ir[i];

temp0_ii[i/2] = ii[i];

temp1_ir[i/2] = ir[i+1];

temp1_ii[i/2] = ii[i+1];

}

//for(i=0; i<2; i++)

//{

// printf("%.2f %.2f %.2f

%.2f\n",temp0_ir[i],temp0_ii[i],temp1_ir[i],temp1_ii[i]);

//}

Recursive_FFT(temp0_ir, temp0_ii, temp0_or, temp0_oi, len-1);

Recursive_FFT(temp1_ir, temp1_ii, temp1_or, temp1_oi, len-1);

for(k=0; k<=n/2-1; k++)

{

or[k] = temp0_or[k] + (wr*temp1_or[k] - wi*temp1_oi[k]);

oi[k] = temp0_oi[k] + (wr*temp1_oi[k] + wi*temp1_or[k]);

or[k+n/2] = temp0_or[k] - (wr*temp1_or[k] - wi*temp1_oi[k]);

oi[k+n/2] = temp0_oi[k] - (wr*temp1_oi[k] + wi*temp1_or[k]);

wtr = wr*wnr - wi*wni;

wti = wr*wni + wi*wnr;

wr = wtr;

wi = wti;

}

}

}