FFT 算法实现

fft[/b] 算法实现 （转）
2007年10月17日 星期三 16:04
/**//***********************************
// dit radix-4   fft[/b]   complex
//
// 1. maximium points are 2048 because that
//     function mylog2(int n) has the limit of maximal points
//
// 2. the last stage is 2 dft ( for 8, 32, 128, 512, 2048)
// or 4 dft ( for 4, 16, 64, 256, 1024 ).
//
//
// 24 juillet 2007
// purcharse*gmail.com
************************************/
#include
#include
#include
typedef double real64;       /**//* floating point */
typedef float   real32;     /**//* 32 bit fixed point */
struct complex
{
        real32 real;
        real32 imag;
} complex ;
static struct complex multicomplex(struct complex,struct complex);
static int mylog2(int);
static void dft_2(struct complex *,struct complex *);
static void dft_4(struct complex *,struct complex *,struct complex *,struct complex *);
static void runfft(struct complex *, int);
static void runifft(struct complex *,int);
static void bitreverse(struct complex *, int);
static void fft[/b]_r4(struct complex *, int, int);
static void fft[/b]_l2(struct complex *, int);
struct   complex s[257];
int    num = 16;
const float pi=3.1415926535898;
main()
{
        int i;
        /**//* rectangle */
        /**//*
        for(i=0;i
        s[0].real = 10;
        s[0].imag = 10;
        */
                /**//*   sinus*/
                for(i=0;i
        /**//*   */
        /**//*
        for(i=0;i
        printf(*********** donnees ***************** \n);
        for(i=0;i
        runfft(s, num);
        printf(*********** fft[/b] ***************** \n);
        for(i=0;i
static struct complex multicomplex(struct complex b1,struct complex b2)   /**//* multiplication of complex */
{
        struct complex b3;
        b3.real=b1.real*b2.real-b1.imag*b2.imag;
        b3.imag=b1.real*b2.imag+b1.imag*b2.real;
        return(b3);
}
static int mylog2(int n)      /**//* max(n) = 4098 */
{
        int k=0;
        if (n>>12) { k+=12; n>>=12; }
        if (n>>8) { k+=8; n>>=8; }
        if (n>>4) { k+=4; n>>=4; }
        if (n>>2) { k+=2; n>>=2; }
        if (n>>1) { k+=1; n>>=1; }
        return k ;
}
static void bitreverse(struct complex *xin, int n)
{
        int lh, i, j, k;
        struct complex tmp;
        lh=n/2;    
        j = n/2;
        for( i = 1; i = k)
                {
                        j = j-k;
                        k = k/2;
                }
                j = j + k;
        }
}
static void dft_2(struct complex *b1,struct complex *b2)
{
        struct complex tmp;
        tmp = *b1;
        (*b1).real = (*b1).real + (*b2).real;
        (*b1).imag = (*b1).imag + (*b2).imag;
        (*b2).real = tmp.real - (*b2).real;
        (*b2).imag = tmp.imag - (*b2).imag;    
}
static void dft_4(struct complex* b0, struct complex* b1, struct complex* b2,struct complex* b3)
{
        /**//*variables locales*/
        struct complex temp[4];
        /**//*calcul x1*/
        temp[0].real=(*b0).real+(*b1).real;        
        temp[0].imag=(*b0).imag+(*b1).imag;
        /**//*calcul x2*/
        temp[1].real=(*b0).real-(*b1).real;        
        temp[1].imag=(*b0).imag-(*b1).imag;
        /**//*calcul x3*/
        temp[2].real=(*b2).real+(*b3).real;    
        temp[2].imag=(*b2).imag+(*b3).imag;
        /**//*calcul x4 + multiplication with -j*/
        temp[3].imag=(*b3).real-(*b2).real;        
        temp[3].real=(*b2).imag-(*b3).imag;
        /**//*the last stage*/
        (*b0).real=temp[0].real+temp[2].real;
        (*b0).imag=temp[0].imag+temp[2].imag;
        (*b1).real=temp[1].real+temp[3].real;    
        (*b1).imag=temp[1].imag+temp[3].imag;
        (*b2).real=temp[0].real-temp[2].real;    
        (*b2).imag=temp[0].imag-temp[2].imag;
        (*b3).real=temp[1].real-temp[3].real;
        (*b3).imag=temp[1].imag-temp[3].imag;
}
static void fft[/b]_r4(struct complex *xin, int n, int m)
{
        int   i, l, j;
        double   ps1, ps2, ps3, p ;
        int le,b;
        struct complex w[4];
        for( l = 1; l
                p = n/le;
                for(j = 0; j
                        ps1 = (2*pi/n)*p*2*j;
                        w[1].real = cos(ps1);
                        w[1].imag = -sin(ps1);
                        ps2 = (2*pi/n)*p*1*j;
                        w[2].real = cos(ps2);
                        w[2].imag = -sin(ps2);
                        ps3 = (2*pi/n)*p*3*j;
                        w[3].real = cos(ps3);
                        w[3].imag = -sin(ps3);
                        for(i = j; i
}//fin du fft[/b]_r4
static void fft[/b]_l2(struct complex *xin, int n)
{                 /**//* for the last stage 2 dft*/
        int j, b;
        double p, ps ;
        struct complex w;
        b = n/2;              
        for(j = 0; j
                /**//* multiple avec w */
                xin[j+ b] = multicomplex(xin[j + b], w);
                dft_2(xin + j ,xin + j + b);
        }
}//fin du fft[/b]_l2
static void runfft(struct complex *xin, int n)
{
        int m, i;
        bitreverse(xin, n);
        m = mylog2(n);
        if( (m%2) == 0 )
        {
                /**//*all the stages are radix 4*/
                fft[/b]_r4(xin, n, m/2);
        }
        else
        {
                /**//*the last stage is radix 2*/
                fft[/b]_r4(xin, n, m/2);
                fft[/b]_l2(xin, n);
        }
}
static void runifft(struct complex *xin,int n)
{             /**//* inverse fft[/b] */
        int i;
        for(i=0; i
----------------------------------------------
此篇文章从博客转发
原文地址： Http://blog.gkong.com/more.asp?id=89182&Name=mahuaqiang