c语言递归和迭代作业答案,快速傅立叶变换（FFT）串行迭代和递归实现（C语言）...-CSDN博客

#include

#define pi 3.14159265358979

const int MAX_N = 256; // 最大DFT点数

//======================二进制位倒序排列函数==============================

//将整数src的二进制位按倒序排列后返回，其中size为整数src的二进制的长度

int BitReverse(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;

}

//=====================内部的一个求2的次幂的函数==============================

inline int pow2(int n)

{

return 1<

}

//========================FFT:Fast Fourier Transform==========================

//========================Iterative_FFT=======================================

//快速傅立叶变换算法的串行迭代实现

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;

}

//========================Recursive_FFT====================================

//快速傅立叶变换算法的串行递归实现

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

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;

}

//===================主函数========================

int main()

{

double fr[4],fi[4],gr[4],gi[4];

int i;

memset(gr,0,sizeof(gr));

memset(gi,0,sizeof(gi));

//==========一组测试用的数据==================

//==========(1,2,4,3)的傅立叶变换为(10,-3-i,0,-3+i)=========

fr[0]=1;fi[0]=0;

fr[1]=2;fi[1]=0;

fr[2]=4;fi[2]=0;

fr[3]=3;fi[3]=0;

Iterative_FFT(fr,fi,gr,gi,2);

//Recursive_FFT(fr,fi,gr,gi,2);

printf("原始数据:\n");

for(i=0;i<4;i++)

{

printf("%.2f+(%.2f)i\n",fr[i],fi[i]);

}

printf("FFT:\n");

for(i=0;i<4;i++)

{

printf("%.2f+(%.2f)i\n",gr[i],gi[i]);

}

return 0;

}