fft.c

这是用C语言实现数字信号处理中的FFT变换的代码。。可以很方便的进行快速傅里叶变换

/*-----------------------------------------------------------------//
//                 Fast Fourier Transform program                  //
//-----------------------------------------------------------------//
//  Algrithem: Cooley-Tukey Algrithem ( decimation in frequency )  //
//                                                                 // 
//   xr  : Real part of original data as the input;                // 
//         Real part of the output.                                //
//   xi  : Imaginary part of original data as the input;           //
//         Imaginary part of the output.                           //
//   nr  : must be an positive integer,                            // 
//         N=pow(2,nu) is the length of input array (points).      //
//   T   : =1.0 is forward transform;  =-1.0 is reverse transform  //
//                                                                 //
//-----------------------------------------------------------------*/
#include <math.h>
#include <stdio.h>
void fft(float *xr, float *xi, int nr, float T);
void  main ()
{
 float  ra[18],ri[18];
 int i;
 for(i=0;i<8;i++)
 {
	ra[i]=0.0;
	ri[i]=0.0;
 }
   ra[4]=10.0;
   fft(ra,ri,3,1.0);
   printf("%f", ra[1]);

   fft(ra,ri,3,-1.0);
   for(i=0;i<8;i++)
   printf("\n   %f",  ra[i]);
}

//===============================================================//
void fft(float *xr, float *xi, int nr, float T)
{
	int   i, j, k, l, n, n2, nr1, i1, j1, k2, k1n2;
	float c, s, p, tr, ti, trc, tic, ars;

    n=(int)(pow(2,nr));
    n2=n/2;
    nr1=nr-1;

    k=0;
    for(l=1; l<=nr; l++)
	{
loop1:  for(j=1; j<=n2; j++)
		{
	    k2=k/(int)(pow(2,nr1));
	    p=(float)(IBIT(k2,nr));
	    ars=(float)(6.2831852*p/(float)(n));
	    c=(float)(cos(ars));
	    s=(float)(sin(ars));
            k1n2=k+n2;
            tr=xr[k1n2]*c+xi[k1n2]*s*T;
            ti=xi[k1n2]*c-xr[k1n2]*s*T;
            xr[k1n2]=xr[k]-tr;
            xi[k1n2]=xi[k]-ti;
            xr[k]=xr[k]+tr;
            xi[k]=xi[k]+ti;
            k++;		
		}
		k+=n2;
        if(k<n)
		{  goto loop1;  }
		else
		{
			k=0;
            nr1=nr1-1;
            n2 /=2;
		}
	}

    for(j=1; j<=n; j++)
	{
		i=IBIT(j-1,nr)+1;
		if(i>j)
		{
			j1=j-1;
			i1=i-1;
			trc   =xr[j1];
            tic   =xi[j1];
            xr[j1]=xr[i1];
            xi[j1]=xi[i1];
            xr[i1]=trc;
            xi[i1]=tic;
		}
	}

    if(T<0.0)
    {
		for(j=0; j<=n; j++)
		{
			xr[j]/=n;
		    xi[j]/=n;
		}
    }
}

int IBIT(int j, int m)
{
	int i, it, j1, j2;
    it=0;
    j1=j;
    for(i=1; i<=m; i++)
	{
		j2=j1/2;
		it=it*2+(j1-2*j2);
		j1=j2;
	}
    return it;
}