FFT.rar_SRFFT.c_site:www.pudn.com_分裂傅里叶_分裂基_分裂基FFT

//--------------------------------------------------------------------------- #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;