/**** 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;
}

评论