fft_c.cpp

基于C实现的fft源代码.编译通过,对于学习学习FFT算法的朋友很有帮助

/*****************fft programe*********************/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include <math.h>

#define NUM 512 

struct compx
{
	double real;
	double imag;
};


struct compx EE(struct compx b1,struct compx b2)//两个复数相乘
{
	struct compx b3 ;
	b3.real=b1.real*b2.real-b1.imag*b2.imag ;
	b3.imag=b1.real*b2.imag+b1.imag*b2.real ;
	return(b3);
}

void FFT(struct compx*xin,int N)
{
	int cal_num = 0;
    int f,m,nv2,nm1,i,k,j=1,l ;
    /*int f,m,nv2,nm1,i,k,j=N/2,l;*/
    struct compx v,w,t ;
    nv2=N/2 ;
    f=N ;
    for(m=1;(f=f/2)!=1;m++)
    {
        ;
    }
    nm1=N-1 ;
    
    /*变址运算*/
    for(i=1;i<=nm1;i++)
    {
        if(i<j)
        {
            t=xin[j];
            xin[j]=xin[i];
            xin[i]=t ;
        }
        k=nv2 ;
        while(k<j)
        {
            j=j-k ;
            k=k/2 ;
        }
        j=j+k ;
    }
    
	double le,pi;
	int lei,ip ;
    for(l=1;l<=m;l++)
    {
        le=pow(2,l);
        // 这里用的是L而不是1   !!!!
        lei=le/2 ;
        pi=3.1415926 ;
        v.real=1.0 ;
        v.imag=0.0 ;
        w.real=cos(pi/lei);
        w.imag=-sin(pi/lei);
        for(j=1;j<=lei;j++)
        {
            /*double p=pow(2,m-l)*j;
             double ps=2*pi/N*p;
             w.real=cos(ps);
             w.imag=-sin(ps);*/
			for(i=j;i<=N;i=i+le)
            {
                /*   w.real=cos(ps);
                       w.imag=-sin(ps);*/
				ip=i+lei ;
                t=EE(xin[ip],v);
                xin[ip].real=xin[i].real-t.real ;
                xin[ip].imag=xin[i].imag-t.imag ;
                xin[i].real=xin[i].real+t.real ;
                xin[i].imag=xin[i].imag+t.imag ;
				cal_num++;
            }
			v=EE(v,w);
			printf("(%f,%f)\n",v.real,v.imag);
        }
	}
	printf("cal number = %d\n",cal_num);
    return ;
}

/*****************main programe********************/

double result[NUM];
struct compx s[NUM];

const double pp=3.1415926;

int main(void)
{
     int i=1 ;
     for(;i<0x101;i++)
     {
         s[i].real=sin(pp*i/32);
         s[i].imag=0 ;
     }
    
     FFT(s,NUM);
    
     for(i=1;i<NUM;i++)
     {
         result[i]=sqrt(pow(s[i].real,2)+pow(s[i].imag,2));
     }

	 return 0;
}