dllfft.cpp

实数序列的快速傅立叶变换

//---------------------------------------------------------------------------

#include <vcl.h>
#include <windows.h>
#include "math.h"
#include "stdio.h"
#include "stdlib.h"

const float pi=3.14156;

#pragma hdrstop
//---------------------------------------------------------------------------
//   Important note about DLL memory management when your DLL uses the
//   static version of the RunTime Library:
//
//   If your DLL exports any functions that pass String objects (or structs/
//   classes containing nested Strings) as parameter or function results,
//   you will need to add the library MEMMGR.LIB to both the DLL project and
//   any other projects that use the DLL.  You will also need to use MEMMGR.LIB
//   if any other projects which use the DLL will be performing new or delete
//   operations on any non-TObject-derived classes which are exported from the
//   DLL. Adding MEMMGR.LIB to your project will change the DLL and its calling
//   EXE's to use the BORLNDMM.DLL as their memory manager.  In these cases,
//   the file BORLNDMM.DLL should be deployed along with your DLL.
//
//   To avoid using BORLNDMM.DLL, pass string information using "char *" or
//   ShortString parameters.
//
//   If your DLL uses the dynamic version of the RTL, you do not need to
//   explicitly add MEMMGR.LIB as this will be done implicitly for you
//---------------------------------------------------------------------------

#pragma argsused
int WINAPI DllEntryPoint(HINSTANCE hinst, unsigned long reason, void* lpReserved)
{
  return 1;
}
//---------------------------------------------------------------------------

extern "C" __declspec(dllexport) _stdcall void rlfft(float *x,float *yyr,float *yyi,int n);


//---------------------------------------------------------------------------


//---------------------------------------------------------------------------

#pragma package(smart_init)

typedef struct {float real,imag;} complex;
/*-------------------------------------------------------------------*/
float mabs(complex a)
{
 float m;
 m=a.real*a.real+a.imag*a.imag;
 m=sqrt(m);
 return(m);
 }
/*-------------------------------------------------------------------*/
float msign(float a,float b)
{
 float z;
 if(b>=0) z=sqrt(pow(a,2));
 else z=-sqrt(pow(a,2));
 return(z);
 }
/*-------------------------------------------------------------------*/
complex cexp(complex a)
{
 complex z;
 z.real=exp(a.real)*cos(a.imag);
 z.imag=exp(a.real)*sin(a.imag);
 return(z);
 }

////////////////////////////////////////////////////////////////////////

//定义函数的外部接口

extern "C" __declspec(dllexport) _stdcall void realfft(float *x,float *yyr,float *yyi,int n);

void msplfft(complex x[],int n,int isign)
{

//分裂基算法
/*----------------------------------------------------------------------
 Routine msplfft:to perform the split-radix DIF fft algorithm.
 input parameters:
 x : complex array.input signal is stored in x(0) to x(n-1).
 n : the dimension of x.
 isign:if isign=-1 For Forward Transform
       if isign=+1 For Inverse Transform.
 output parameters:
 x : complex array. DFT result is stored in x(0) to x(n-1).
 Notes:
     n must be power of 2.
                                        in chapter 5
----------------------------------------------------------------------*/
        complex xt;
        float es,e,a,a3,cc1,ss1,cc3,ss3,r1,r2,s1,s2,s3;
        int m,n2,k,n4,j,is,id,i1,i2,i3,i0,n1,i,nn;
        for(m=1;m<=16;m++)
    	  {
    	    nn=pow(2,m);
    	    if(n==nn)break;
    	  }
	
	    n2=n*2;
        es=-isign*atan(1.0)*8.0;
   for(k=1;k<m;k++)
     {
            n2=n2/2;
            n4=n2/4;
            e=es/n2;
            a=0.0;
            for(j=0;j<n4;j++)
               {
                a3=3*a;
                cc1=cos(a);
                ss1=sin(a);
                cc3=cos(a3);
                ss3=sin(a3);
                a=(j+1)*e;
                is=j;
                id=2*n2;
                 do
                  {
                    for(i0=is;i0<n;i0+=id)
                       {
                           i1=i0+n4;
                            i2=i1+n4;
                            i3=i2+n4;
                            r1=x[i0].real-x[i2].real;
                            s1=x[i0].imag-x[i2].imag;
                            r2=x[i1].real-x[i3].real;
                            s2=x[i1].imag-x[i3].imag;
                            x[i0].real+=x[i2].real;x[i0].imag+=x[i2].imag;
                            x[i1].real+=x[i3].real;x[i1].imag+=x[i3].imag;
                		    if(isign!=1)
                    		    {
                    		      s3=r1-s2;
                    		      r1+=s2;
                    		      s2=r2-s1;
                    		      r2+=s1;
                    		    }
                    		 else
                    		    {
                    		      s3=r1+s2;
                    		      r1=r1-s2;
                    		      s2=-r2-s1;
                    		      r2=-r2+s1;
                                 }
                    		x[i2].real=r1*cc1-s2*ss1;
                    		x[i2].imag=-s2*cc1-r1*ss1;
                    		x[i3].real=s3*cc3+r2*ss3;
                    		x[i3].imag=r2*cc3-s3*ss3;
                        }
                   is=2*id-n2+j;
                  id=4*id;
		  }while(is<n-1);
      }
  }
/*   ------------ special last stage -------------------------*/
        is=0;
        id=4;
      do
      {
        for(i0=is;i0<n;i0+=id)
           {i1=i0+1;
            xt.real=x[i0].real;
            xt.imag=x[i0].imag;
            x[i0].real=xt.real+x[i1].real;
            x[i0].imag=xt.imag+x[i1].imag;
            x[i1].real=xt.real-x[i1].real;
            x[i1].imag=xt.imag-x[i1].imag;
            }
        is=2*id-2;
        id=4*id;
       }while(is<n-1);
        j=1;
        n1=n-1;
        for(i=1;i<=n1;i++)
	   {
    	   if(i<j)
    	       {
        		 xt.real=x[j-1].real;
        		 xt.imag=x[j-1].imag;
        		 x[j-1].real=x[i-1].real;
        		 x[j-1].imag=x[i-1].imag;
        		 x[i-1].real=xt.real;
        		 x[i-1].imag=xt.imag;
        		 }
    	      k=n/2;
            do
        	  {
            		if(k>=j)
                		 break;
            	    j-=k;
            	    k/=2;
        	    }while(1);
            j+=k;
        }
        if(isign==-1)
           return;
        for(i=0;i<n;i++)
    	   {
    	    x[i].real/=(float)n;
            x[i].imag/=(float)(n);
          }
 }


//////////////////////////////////////////////////////////////////////

void _stdcall rlfft(float *x,float *yyr,float *yyi,int n)
{
    ////////////////////////////////////////////////////////////////
    //                                                            //
    //                 此函数用来作实数的FFT变换                 //
    //            其核心是将实数分为两半，先作n/2的fft           //
	//            再组合得到N点的FFT.采用的是分裂基算法          //
	//            x输入的数组　yyr,yyi输出的实部和虚部           //
	//            n数组长度　　　　　　　　　　　　　　　　　　　//
	//                                                           //
    ///////////////////////////////////////////////////////////////

    int n2=n/2;
    int n4=n/4;

     //存放x的偶数序数值
    complex *y=new complex[n2];

   
    for(int i=0;i<n2;i++)
    {
        y[i].real =x[2*i];
        y[i].imag=x[2*i+1];
    }
    //利用分裂基fft方法对n2点数据进行FFT变换

     msplfft(y,n2,-1);        //y保存n/2点的fft变换结果

    //对N/2点的fft数据进行组合，以产生N点的FFT变换
    float *x0=new float[n2];
    float *x1=new float[n2];
    float *yr= new float[n];
    float *yi= new float[n];
    for(int i=0;i<n2;i++)
    {
      x0[i]=y[i].real;
      x1[i]=y[i].imag;
    }

    float wr=1.0;
    float wi=0;

    float theta=-8*atan(1.0)/n;
    float ct=cos(theta);
    float st=sin(theta);
    

    yr[n2]=x0[0]-x1[0];
    yi[n2]=0.0;

    yr[0]=x0[0]+x1[0];
    yi[0]=0.0;

    for(int i=1;i<n4;i++)
    {
        float x0_1=x0[i];
        float x1_1=x1[i];
        float x0_2=x0[n2-i];
        float x1_2=x1[n2-i];
        float x0_1_=0.5*(x0_1+x0_2);
        float x0_2_=0.5*(x1_1-x1_2);
        float x1_1_=0.5*(x1_1+x1_2);
        float x1_2_=-0.5*(x0_1-x0_2);
        float  wr_t=wr;

        wr=wr_t*ct-wi*st;
        wi=wr_t*st+wi*ct;
        yr[i]=x0_1_+wr*x1_1_-wi*x1_2_;
        yi[i]=x0_2_+wr*x1_2_+wi*x1_1_;
        yr[n2-i]=x0_1_-wr*x1_1_ +wi*x1_2_;
        yi[n2-i]=-x0_2_+wr*x1_2_+wi*x1_1_;
    }

     yr[n4] =  x0[n4];
     yi[n4] = -x1[n4];

     for (int i = 1; i < n2; i++)
     {
          yr[i + n2] =  yr[n2 - i];
          yi[i + n2] = -yi[n2 - i];
     }

     //将结果保存到复数yy中
     for(int i=0;i<n;i++)
     {
       yyr[i]=yr[i];
       yyi[i]=yi[i];
     }

    delete[] x0;
    delete[] x1;
    delete[] yr;
    delete[] yi;
    delete[] y;

}