00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include <stdlib.h>
00008 #include <stdio.h>
00009 #include <string.h>
00010 #include <math.h>
00011 
00012 #undef PI
00013 #define PI (3.141592653589793238462643)
00014 
00015 typedef struct complex { float r , i ; } complex ;
00016 
00017 static void csfft_print( int mode , int idim , complex * xc ) ;
00018 static void csfft_trigconsts( int idim ) ;
00019 
00020 int main( int argc , char * argv[] )
00021 {
00022    int len , ii ;
00023    complex * cx ;
00024 
00025    if( argc < 2 ){printf("Usage: fftprint len\n");exit(0);}
00026 
00027    len = strtol( argv[1] , NULL , 10 ) ;
00028    if( len < 4 ){ fprintf(stderr,"Illegal length\n"); exit(1); }
00029 
00030    ii = 2 ; do{ ii *= 2 ; } while( ii < len ) ;
00031    if( ii != len ){ fprintf(stderr,"Illegal length\n"); exit(1); }
00032 
00033    cx = (complex *) malloc( sizeof(complex) * len ) ;
00034    csfft_print(-1,len,cx) ;
00035    exit(0) ;
00036 }
00037 
00038 static complex * csplus = NULL , * csminus = NULL ;  /* trig consts */
00039 static int nold = -666 ;
00040 
00041 static void csfft_trigconsts( int idim )  /* internal function */
00042 {
00043    register unsigned int  m, n, i0, i1, i2, i3, k;
00044    register complex       *r0, *csp;
00045    register float         co, si, f0, f1, f2, f3, f4;
00046    double                 al;
00047 
00048    if( idim == nold ) return ;
00049 
00050    if( idim > nold ){
00051       if( csplus != NULL ){ free(csplus) ; free(csminus) ; }
00052       csplus  = (complex *) malloc( sizeof(complex) * idim ) ;
00053       csminus = (complex *) malloc( sizeof(complex) * idim ) ;
00054       if( csplus == NULL || csminus == NULL ){
00055          fprintf(stderr,"\n*** csfft cannot malloc space! ***\n"); exit(1) ;
00056       }
00057    }
00058    nold = n = idim ;
00059 
00060    f1 = 1.0 ;  /* csplus init */
00061    m  = 1; k  = 0;
00062    while (n > m) {
00063       i3 = m << 1;
00064       f2 = m; al = f1*PI/f2;
00065       co = cos(al); si = sin(al);
00066       (csplus + k)->r = 1.;
00067       (csplus + k)->i = 0.;
00068       for (i0=0; i0 < m; i0++) {
00069          k++;
00070          csp = csplus + k; r0 = csp - 1;
00071          csp->r = r0->r * co - r0->i * si;
00072          csp->i = r0->i * co + r0->r * si;
00073       }
00074       m = i3;
00075    }
00076 
00077    f1 = -1.0 ;  /* csminus init */
00078    m  = 1; k  = 0;
00079    while (n > m) {
00080       i3 = m << 1;
00081       f2 = m; al = f1*PI/f2;
00082       co = cos(al); si = sin(al);
00083       (csminus + k)->r = 1.;
00084       (csminus + k)->i = 0.;
00085       for (i0=0; i0 < m; i0++) {
00086          k++;
00087          csp = csminus + k; r0  = csp - 1;
00088          csp->r = r0->r * co - r0->i * si;
00089          csp->i = r0->i * co + r0->r * si;
00090       }
00091       m = i3;
00092    }
00093    return ;
00094 }
00095 
00096 /*****************************************************************************
00097   The routine that generates a C function for a completely unrolled FFT.
00098   Basically the same as csfft_cox() with print statements replacing
00099   the actual computations.
00100   The resulting code also needs the csfft_trigconsts() function above.
00101 ******************************************************************************/
00102 
00103 #define EPS 1.e-5  /* set EPS to 0 to disable this option */
00104 
00105 void csfft_print( int mode , int idim , complex * xc )
00106 {
00107    register unsigned int  m, n, i0, i1, i2, i3, k;
00108    register complex       *r0, *r1, *csp;
00109    register float         co, si, f0, f1, f2, f3, f4;
00110    double                 al;
00111 
00112    csfft_trigconsts( idim ) ;
00113    csp = csplus ;
00114 
00115    n   = idim;
00116    i2  = idim >> 1;
00117    i1  = 0;
00118 
00119    /*-- function header --*/
00120 
00121    printf( "/**************************************/\n"
00122            "/* FFT routine unrolled of length %3d */\n"
00123            "\n"
00124            "void fft%d( int mode , complex * xc )\n"
00125            "{\n"
00126            "   register complex * csp , * xcx=xc;\n"
00127            "   register float f1,f2,f3,f4 ;\n"
00128            "\n"
00129            "   /** perhaps initialize **/\n"
00130            "\n"
00131            "   if( nold != %d ) csfft_trigconsts( %d ) ;\n"
00132            "\n"
00133            "   csp = (mode > 0) ? csplus : csminus ;  /* choose const array */\n"
00134            "\n"
00135            "   /** data swapping part **/\n" ,
00136       idim,idim,idim,idim ) ;
00137 
00138    /*-- swap some data elements --*/
00139 
00140    for (i0=0; i0 < n; i0 ++) {
00141       if ( i1 > i0 ) {
00142 
00143 printf("\n") ;
00144 printf("   f1 = xcx[%d].r ; f2 = xcx[%d].i ;\n",i0,i0) ;
00145 printf("   xcx[%d].r = xcx[%d].r ; xcx[%d].i = xcx[%d].i ;\n" , i0,i1,i0,i1) ;
00146 printf("   xcx[%d].r = f1 ; xcx[%d].i = f2 ;\n" , i1,i1) ;
00147 
00148       }
00149       m = i2;
00150       while ( m && !(i1 < m) ) {
00151          i1 -= m;
00152          m >>= 1;
00153      }
00154      i1 += m;
00155    }
00156 
00157    /*-- the actual computations --*/
00158 
00159    printf("\n   /** butterflying part **/\n") ;
00160 
00161    m = 1;
00162    k = 0;
00163    while (n > m) {
00164       i3 = m << 1;
00165       for (i0=0; i0 < m; i0 ++) {
00166          for (i1=i0; i1 < n; i1 += i3) {
00167 
00168 printf("\n") ;
00169 if( csp[k].r == 1.0 ){
00170   printf("   f1 = xcx[%d].r ; f3 = xcx[%d].i ;  /* cos=1 sin=0 */\n", i1+m,i1+m ) ;
00171 
00172 } else if( fabs(csp[k].r) < EPS ){
00173   printf("   f1 = - xcx[%d].i * csp[%d].i ; /* cos=0 twiddles */\n", i1+m,k ) ;
00174   printf("   f3 = xcx[%d].r * csp[%d].i ;\n", i1+m,k ) ;
00175 
00176 } else {
00177   printf("   f1 = xcx[%d].r * csp[%d].r - xcx[%d].i * csp[%d].i ; /* twiddles */\n",
00178          i1+m,k , i1+m,k ) ;
00179   printf("   f3 = xcx[%d].r * csp[%d].i + xcx[%d].i * csp[%d].r ;\n",
00180          i1+m,k , i1+m,k ) ;
00181 }
00182 
00183 printf("   f2 = xcx[%d].r ; f4 = xcx[%d].i ;\n" , i1,i1 ) ;
00184 printf("   xcx[%d].r = f2-f1 ; xcx[%d].i = f4-f3 ;\n",i1+m,i1+m) ;
00185 printf("   xcx[%d].r = f2+f1 ; xcx[%d].i = f4+f3 ;\n",i1,i1) ;
00186 
00187          }
00188          k++;
00189       }
00190       m = i3;
00191    }
00192 
00193 printf("\n   return ;\n}\n") ;
00194 
00195    return ;
00196 }