main.c

FFT算法

#include <stdio.h>
#include <math.h>


#define pi 3.1415926
#define SIG_LEN    256

typedef struct{
    float real;
    float image;
}Complex;


Complex ComMul(Complex p1,Complex p2);
Complex ComAdd(Complex p1,Complex p2);
Complex ComSub(Complex p1,Complex p2);

int Reform(float **in,int len);
int BitRevers(int src,int size);
void Displace(float *in,int size,int M);
void CalcW(Complex *w,int len);
Complex *fft_1d(float *in,int len);

main()
{
    int slen,i;
  
    FILE *fp,*fout;
    float *data;
    Complex *out;

    
    data = malloc(SIG_LEN*sizeof(float));

    if((fp=fopen("e:\input","r"))==NULL)
    {
        printf("cannot open file");
        exit(0);
    }
    
    if((fout=fopen("e:\output","w+"))==NULL)
    {
        printf("cannot open file");
        exit(0);
    }

    for(i=0;i<SIG_LEN;i++)
        fscanf(fp,"%f\n",&data[i]);
    
    out = fft_1d(data,SIG_LEN);
    for(i=0;i<SIG_LEN;i++)
        fprintf(fout,"%f\n",sqrt(pow(out[i].real,2)+pow(out[i].image,2)));
}


Complex ComMul(Complex p1,Complex p2)
{
    Complex res;
    res.real = p1.real*p2.real - p1.image*p2.image;
    res.image = p1.real*p2.image + p1.image*p2.real;
    
    return res;
}

Complex ComAdd(Complex p1,Complex p2)
{
    Complex res;
    res.real = p1.real + p2.real;
    res.image = p1.image + p2.image;
    
    return res;
}

Complex ComSub(Complex p1,Complex p2)
{
    Complex res;
    res.real = p1.real - p2.real;
    res.image = p1.image - p2.image;
    
    return res;
}
        

int BitRevers(int src,int size)
{
    int temp=src;
    int dst=0;
    int i=0;
    for(i=size-1;i>=0;i--)
    {
        dst=((temp&0x1)<<i)|dst;
        temp=temp>>1;
    }
    return dst;
}

void Displace(float *in,int size,int M)
{
    int i;
    int new_i;
    float t;
    FILE *fp;
    for(i=1;i<size;i++)
    {
        new_i=BitRevers(i,M);
        if(new_i>i)
        {
            t=in[i];
            in[i]=in[new_i];
            in[new_i]=t;
        }
     }
    if((fp=fopen("e:\displace","w+"))==NULL)
    {
        printf("cannot open file");
        exit(0);
    }
    for(i=0;i<size;i++)
        fprintf(fp,"%f\n",in[i]);
}
void CalcW(Complex *w,int N)
{
    int i;
    FILE    *fp;
    for(i=0;i<N;i++)
    {
        w[i].real = cos(pi*i/N);
        w[i].image = -1*sin(pi*i/N);
    }
    if((fp=fopen("e:\w","w+"))==NULL)
    {
        printf("cannot open file");
        exit(0);
    }
    for(i=0;i<N;i++)
        {fprintf(fp,"        .word    %ld\n",(long)(w[i].real*16384));
         fprintf(fp,"        .word    %ld\n",(long)(w[i].image*16384));
       }

}

int Reform(float **in,int len)
{
    int i=0;
    int w=1;
    
    while(w*2<=len) {w=w*2;i++;}

    if(w<len)
    {
        
        *in = realloc(*in,w*2*sizeof(float));
        for(i=len;i<w*2;i++)
            in[i]=0;
        return i+1;
    }
    return i;
}
Complex *fft_1d(float *data,int len)
{
    int M,i=0,j,k;
    int GroupNum;
    int CellNum;
    int reallen;
    int pos1,pos2;
    float *in;
    Complex *res,*t;
   
    Complex *w,mul;

    in = data;
    M = Reform(&in,len);                    /*let the data length is pow of 2*/
    reallen = pow(2,M);
    res = malloc(reallen*sizeof(Complex));    
    
    w = malloc(reallen*sizeof(Complex)/2);  
    CalcW(w,reallen/2);

    Displace(in,reallen,M);
    while(i<reallen) {res[i].real = in[i];res[i].image=0;i++;}
    
    

    GroupNum = reallen/2;
    CellNum = 1;
    
    for(i=0;i<M;i++)
    {       
        for(j=0;j<GroupNum;j++)
        {
            for(k=0;k<CellNum;k++)
            {
                pos1 = j*CellNum*2 + k;
                pos2 = pos1 + CellNum;

                mul = ComMul(res[pos2],w[k*GroupNum]);
                
                res[pos2] = ComSub(res[pos1] ,mul);
                res[pos1] = ComAdd(res[pos1] ,mul);
            }
        }
        GroupNum = GroupNum/2;
        CellNum = CellNum*2;    
  
    }
    return res;
}