Main Page | Modules | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Class Members | File Members | Related Pages

FEAPIExamplePluginZeroCrossings/zplFFTGenCore.cpp

Go to the documentation of this file.
00001 //
00002 //                   FFT library
00003 //
00004 //  (one-dimensional complex and real FFTs for array 
00005 //  lengths of 2^n)
00006 //
00007 //  Author: Toth Laszlo (tothl@inf.u-szeged.hu)
00008 //  
00009 //  Research Group on Artificial Intelligence
00010 //  H-6720 Szeged, Aradi vertanuk tere 1, Hungary
00011 //  
00012 //  Last modified: 97.05.29
00014 
00015 #include "zplVecLib.h"
00016 
00017 #if !defined(_INVSQRT2)
00018 #define _INVSQRT2               (float)(0.70710678118654752440084436210485)    
00019 #endif
00020 
00022 // Sorensen in-place split-radix FFT for real values
00023 // data: array of doubles:
00024 // re(0),re(1),re(2),...,re(size-1)
00025 // 
00026 // output:
00027 // re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1)
00028 // normalized by array length
00029 //
00030 // Source: 
00031 // Sorensen et al: Real-Valued Fast Fourier Transform Algorithms,
00032 // IEEE Trans. ASSP, ASSP-35, No. 6, June 1987
00033 
00034 void  ZCDECL zRealFFTSplitGeneric(float *data,int n,float *pSinTable,float *pCosTable)
00035 {
00036     
00037     int i,j,k,i5,i6,i7,i8,i0,id,i1,i2,i3,i4,n2,n4,n8;
00038     float t1,t2,t3,t4,t5,t6,ss1,ss3,cc1,cc3;//,a,a3,e;
00039     
00040     int   iStepsize,index,Step3;
00041     
00042     Step3 = (n>>3)-2;
00043     n4=n-1;
00044     
00045     //data shuffling
00046     for (i=0,j=0,n2=n/2; i<n4 ; i++){
00047         if (i<j){
00048             t1=data[j];
00049             data[j]=data[i];
00050             data[i]=t1;
00051         }
00052         k=n2;
00053         while (k<=j){
00054             j-=k;
00055             k>>=1;  
00056         }
00057         j+=k;
00058     }
00059     
00060     /*----------------------*/
00061     
00062     //length two butterflies    
00063     i0=0;
00064     id=4;
00065     do{
00066         for (; i0<n4; i0+=id){ 
00067             i1=i0+1;
00068             t1=data[i0];
00069             data[i0]=t1+data[i1];
00070             data[i1]=t1-data[i1];
00071         }
00072         id<<=1;
00073         i0=id-2;
00074         id<<=1;
00075     } while ( i0<n4 );
00076     
00077     /*----------------------*/
00078     //L shaped butterflies
00079     n2=2;
00080     
00081     iStepsize=n/4;
00082     
00083     for(k=n;k>2;k>>=1){  
00084         n2<<=1;   
00085         n4=n2>>2; 
00086         n8=n2>>3; 
00087         
00088         //e = (float)(_2_PI/(n2));
00089         i1=0;
00090         id=n2<<1;
00091         do{ 
00092             for (; i1<n; i1+=id){
00093                 i2=i1+n4;
00094                 i3=i2+n4;
00095                 i4=i3+n4;
00096                 t1=data[i4]+data[i3];
00097                 data[i4]-=data[i3];
00098                 data[i3]=data[i1]-t1;
00099                 data[i1]+=t1;
00100                 if (n4!=1){
00101                     i0=i1+n8;
00102                     i2+=n8;
00103                     i3+=n8;
00104                     i4+=n8;
00105                     t1=(data[i3]+data[i4])*_INVSQRT2;
00106                     t2=(data[i3]-data[i4])*_INVSQRT2;
00107                     data[i4]=data[i2]-t1;
00108                     data[i3]=-data[i2]-t1;
00109                     data[i2]=data[i0]-t2;
00110                     data[i0]+=t2;
00111                 }
00112             }
00113             id<<=1;
00114             i1=id-n2;
00115             id<<=1;
00116         } while ( i1<n );
00117         //a=e;
00118         index=0;
00119         
00120         for (j=2; j<=n8; j++){  
00121             //a3=3*a;
00122             index+=iStepsize;
00123             
00124             //cc1=(float)cos(a);
00125             cc1=pCosTable[index-1];
00126             //ss1=(float)sin(a);
00127             ss1=pSinTable[index-1];
00128             //cc3=(float)cos(a3);
00129             cc3=pCosTable[index+Step3];
00130             //ss3=(float)sin(a3);
00131             ss3=pSinTable[index+Step3];
00132             //a=j*e;
00133             
00134             i=0;
00135             id=n2<<1;
00136             do{
00137                 for (; i<n; i+=id){  
00138                     i1=i+j-1;
00139                     i2=i1+n4;
00140                     i3=i2+n4;
00141                     i4=i3+n4;
00142                     i5=i+n4-j+1;
00143                     i6=i5+n4;
00144                     i7=i6+n4;
00145                     i8=i7+n4;
00146                     t1=data[i3]*cc1+data[i7]*ss1;
00147                     t2=data[i7]*cc1-data[i3]*ss1;
00148                     t3=data[i4]*cc3+data[i8]*ss3;
00149                     t4=data[i8]*cc3-data[i4]*ss3;
00150                     t5=t1+t3;
00151                     t6=t2+t4;
00152                     t3=t1-t3;
00153                     t4=t2-t4;
00154                     t2=data[i6]+t6;
00155                     data[i3]=t6-data[i6];
00156                     data[i8]=t2;
00157                     t2=data[i2]-t3;
00158                     data[i7]=-data[i2]-t3;
00159                     data[i4]=t2;
00160                     t1=data[i1]+t5;
00161                     data[i6]=data[i1]-t5;
00162                     data[i1]=t1;
00163                     t1=data[i5]+t4;
00164                     data[i5]-=t4;
00165                     data[i2]=t1;
00166                 }
00167                 id<<=1;
00168                 i=id-n2;
00169                 id<<=1;
00170             } while(i<n);
00171         }
00172         iStepsize>>=1;
00173     }
00174     
00175 }
00176 
00177 
00178 
00180 // Sorensen in-place inverse split-radix FFT for real values
00181 // data: array of doubles:
00182 // re(0),re(1),re(2),...,re(size/2),im(size/2-1),...,im(1)
00183 // 
00184 // output:
00185 // re(0),re(1),re(2),...,re(size-1)
00186 // NOT normalized by array length
00187 //
00188 // Source: 
00189 // Sorensen et al: Real-Valued Fast Fourier Transform Algorithms,
00190 // IEEE Trans. ASSP, ASSP-35, No. 6, June 1987
00191 
00192 void ZCDECL zRealIFFTSplitGeneric(float *data,int n,float *pSinTable,float* pCosTable)
00193 {
00194     
00195     int i,j,k,i5,i6,i7,i8,i0,id,i1,i2,i3,i4,n2,n4,n8,n1;
00196     float t1,t2,t3,t4,t5,ss1,ss3,cc1,cc3;//,a,a3,e;
00197     
00198     int   iStepsize,index,Step3;
00199     
00200     Step3 = (n>>3)-2;
00201     
00202     n1=n-1;
00203     n2=n<<1;
00204     
00205     iStepsize=1;
00206     
00207     for(k=n;k>2;k>>=1){  
00208         id=n2;
00209         n2>>=1;
00210         n4=n2>>2;
00211         n8=n2>>3;
00212         //e = (float)(_2_PI/(n2));
00213         i1=0;
00214         do{ 
00215             for (; i1<n; i1+=id){
00216                 i2=i1+n4;
00217                 i3=i2+n4;
00218                 i4=i3+n4;
00219                 t1=data[i1]-data[i3];
00220                 data[i1]+=data[i3];
00221                 data[i2]*=2;
00222                 data[i3]=t1-2*data[i4];
00223                 data[i4]=t1+2*data[i4];
00224                 if (n4!=1){
00225                     i0=i1+n8;
00226                     i2+=n8;
00227                     i3+=n8;
00228                     i4+=n8;
00229                     t1=(data[i2]-data[i0])*_INVSQRT2;
00230                     t2=(data[i4]+data[i3])*_INVSQRT2;
00231                     data[i0]+=data[i2];
00232                     data[i2]=data[i4]-data[i3];
00233                     data[i3]=2*(-t2-t1);
00234                     data[i4]=2*(-t2+t1);
00235                 }
00236             }
00237             id<<=1;
00238             i1=id-n2;
00239             id<<=1;
00240         } while ( i1<n1 );
00241         //a=e;
00242         index=0;
00243         
00244         for (j=2; j<=n8; j++){  
00245             //a3=3*a;
00246             
00247             index+=iStepsize;
00248             
00249             //cc1=(float)cos(a);
00250             cc1=pCosTable[index-1];
00251             //ss1=(float)sin(a);
00252             ss1=pSinTable[index-1];
00253             //cc3=(float)cos(a3);
00254             cc3=pCosTable[index+Step3];
00255             //ss3=(float)sin(a3);
00256             ss3=pSinTable[index+Step3];
00257             
00258             //a=j*e;
00259             i=0;
00260             id=n2<<1;
00261             do{
00262                 for (; i<n; i+=id){  
00263                     i1=i+j-1;
00264                     i2=i1+n4;
00265                     i3=i2+n4;
00266                     i4=i3+n4;
00267                     i5=i+n4-j+1;
00268                     i6=i5+n4;
00269                     i7=i6+n4;
00270                     i8=i7+n4;
00271                     t1=data[i1]-data[i6];
00272                     data[i1]+=data[i6];
00273                     t2=data[i5]-data[i2];
00274                     data[i5]+=data[i2];
00275                     t3=data[i8]+data[i3];
00276                     data[i6]=data[i8]-data[i3];
00277                     t4=data[i4]+data[i7];
00278                     data[i2]=data[i4]-data[i7];
00279                     t5=t1-t4;
00280                     t1+=t4;
00281                     t4=t2-t3;
00282                     t2+=t3;
00283                     data[i3]=t5*cc1+t4*ss1;
00284                     data[i7]=-t4*cc1+t5*ss1;
00285                     data[i4]=t1*cc3-t2*ss3;
00286                     data[i8]=t2*cc3+t1*ss3;
00287                 }
00288                 id<<=1;
00289                 i=id-n2;
00290                 id<<=1;
00291             } while(i<n1);
00292         }
00293         iStepsize<<=1;
00294     }   
00295     
00296     /*----------------------*/
00297     i0=0;
00298     id=4;
00299     do{
00300         for (; i0<n1; i0+=id){ 
00301             i1=i0+1;
00302             t1=data[i0];
00303             data[i0]=t1+data[i1];
00304             data[i1]=t1-data[i1];
00305         }
00306         id<<=1;
00307         i0=id-2;
00308         id<<=1;
00309     } while ( i0<n1 );
00310     
00311     /*----------------------*/
00312     
00313     //data shuffling
00314     for (i=0,j=0,n2=n/2; i<n1 ; i++){
00315         if (i<j){
00316             t1=data[j];
00317             data[j]=data[i];
00318             data[i]=t1;
00319         }
00320         k=n2;
00321         while (k<=j){
00322             j-=k;
00323             k>>=1;  
00324         }
00325         j+=k;
00326     }   
00327 }
00328 
00329 
00330 

Generated on Fri Mar 23 10:28:54 2007 for FEAPI Plugin Documentation by  doxygen 1.3.9.1