/**** fft programe ***** Dev-C++ 4.9.9.2 ************/// 02-11-07 17:51/* 修改 radix 使循环控制完全取决于级数，为流水工作做准备 */ // 用不同的 RAM 存储中间结果 ,单此时仍旧是"原位",需要修改保证每次从不同的RAM读取一个数据 #include<math.h>#include<stdio.h>#include<stdlib.h> #define STAGE 2 #define DISPLAY 1#define Num (1<<STAGE*2)#if DISPLAY    int count = 0;#endif#define START 0    #define END 8        #define WINTH   ((1<<15)-1)#define F(t)  ((int)(WINTH*t)&0xffff)double result[Num];double x[STAGE+1][Num],y[STAGE+1][Num];double TAB[Num*2];const double PI = 3.14159265358979323846; const double EI = PI*2/Num;void rad4(double x0[],double y0[],double x1[],double y1[],double TAB[],                 int i,int n2,int adr1,int adr2,int adr3){    int i1,i2,i3;    double co1,co2,co3,si1,si2,si3;    double r1,r2,r3,r4,s1,s2,s3,s4;        i1=i+n2;i2=i1+n2;i3=i2+n2;                       /* 数据位置 */    #if DISPLAY         if(count>=START && count<END){            printf("data: (%d)%.1f+%.1f, (%d)%.1f+%.1f, (%d)%.1f+%.1f, (%d)%.1f+%.1f\n",i,x0[i],y0[i],i1,x0[i1],y0[i1],i2,x0[i2],y0[i2],i3,x0[i3],y0[i3]);        }    #endif    co1=TAB[2*adr1];co2=TAB[2*adr2];co3=TAB[2*adr3]; /* 旋转因子 */    si1=TAB[2*adr1+1];si2=TAB[2*adr2+1];si3=TAB[2*adr3+1];    #if DISPLAY         if(count>=START && count<END){            printf("factor: (%d)%4x+%xj, (%d)%4x+%xj, (%d)%4x+%xj\n",adr1,F(co1),F(si1),adr2,F(co2),F(si2),adr3,F(co3),F(si3));        }    #endif    r1=x0[i]+x0[i2];r3=x0[i]-x0[i2];    s1=y0[i]+y0[i2];s3=y0[i]-y0[i2];    r2=x0[i1]+x0[i3];r4=x0[i1]-x0[i3];    s2=y0[i1]+y0[i3];s4=y0[i1]-y0[i3];    x1[i]=r1+r2; y1[i]=s1+s2;      r2=r1-r2; r1=r3-s4;r3=r3+s4;    s2=s1-s2; s1=s3+r4;s3=s3-r4;    x1[i1]=co1*r3+si1*s3;y1[i1]=co1*s3-si1*r3;    x1[i2]=co2*r2+si2*s2;y1[i2]=co2*s2-si2*r2;    x1[i3]=co3*r1+si3*s1;y1[i3]=co3*s1-si3*r1;    #if DISPLAY         if(count>=START & count<END){            printf("result: (%d)%.1f+%.1f, (%d)%.1f+%.1f, (%d)%.1f+%.1f, (%d)%.1f+%.1f\n",i,x1[i],y1[i],i1,x1[i1],y1[i1],i2,x1[i2],y1[i2],i3,x1[i3],y1[i3]);        }    #endif}void switchadr(double x[],double y[],int n){    double tempx,tempy;    int i,j,k;    for(i=0,j=0;i<=n-1;i++)    {        if(i<j)        {            tempx=x[j]; x[j]=x[i]; x[i]=tempx;              tempy=y[j]; y[j]=y[i]; y[i]=tempy;        }        k=n/4;         while(3*k<(j+1)&&(j>0) )        {            j=j-3*k; k=k/4;        }        j=j+k;    }}void radix4(double x0[],double y0[],double x1[],double y1[],double TAB[],int n,int L) // 只处理第L级 {    int i,j,e,n1,n2;    int a,b,c;        n1 = n >> (2*L);            n2 = n1 >> 2;        e = 1 <<(2*L);    for(j=0;j<n2;j++)    {        a=j*e;  b=a+a; c=a+b;        a=a%n;  b=b%n; c=c%n;        for(i=j;i<n;i+=n1)        {            #if DISPLAY                 if(count>=START & count<END)                printf("Group : %d\n",count) ;           #endif           rad4(x0,y0,x1,y1,TAB,i,n2,a,b,c);           #if DISPLAY                count ++ ;           #endif        }    }}/***************** ShowResult programe *****************//* 显示 0 -- 4^m-1 数据 */void ShowResult(double x[],double y[],const int n){    int i;    double result[n];    printf("n     real(n)      imag(n)      fabs(n)\n") ;    for(i=0;i<n;i++)         {        result[i]=sqrt( x[i]*x[i]+y[i]*y[i] ); /* 取模 */        printf("%d     %f     %f     %f\n",i,x[i],y[i],result[i]) ;    }}/*****************main programe********************/int main(void){    int i,k;    /*******赋初值**********/    for(i=0;i<Num;i++)    {        x[0][i]=i;   y[0][i]= 0 ;    }    /*******给表格赋初值**********/    for(i=0;i<Num;i++)          /*计算NUM 点的FFT的旋转因子*/     {        TAB[i*2] = cos(EI*i);   /* 存余弦值 */        TAB[i*2+1] = sin(EI*i); /* 存正弦值 */    }    for(k=0;k<STAGE;k++)    /*内部循环没有使用K值，只是用来控制循环次数而已*/     {        radix4(x[k],y[k],x[k+1],y[k+1],TAB,Num,k);        #if !DISPLAY             ShowResult(x[k+1],y[k+1],Num);        #endif    }    switchadr(x[STAGE],y[STAGE],Num);        getchar();    return 0;}

评论