00001 #include "mrilib.h"
00002 
00003 /*--------------------------------------------------------------------
00004   This is the older FFT routine by AJ and RWC,
00005   from which the routines in csfft.c were derived.
00006 ----------------------------------------------------------------------*/
00007 
00008 static complex * csplus = NULL , * csminus = NULL ;  /* trig consts */
00009 static int nold = -666 ;
00010 
00011 #undef PI
00012 #define PI (3.141592653589793238462643)  /* should be close enough */
00013 
00014 static void csfft_trigconsts( int idim )
00015 {
00016    register unsigned int  m, n, i0, i1, i2, i3, k;
00017    register complex       *r0, *csp;
00018    register float         co, si, f0, f1, f2, f3, f4;
00019    double                 al;
00020 
00021    if( idim == nold ) return ;
00022 
00023    if( idim > nold ){
00024       if( csplus != NULL ){ free(csplus) ; free(csminus) ; }
00025       csplus  = (complex *) malloc( sizeof(complex) * idim ) ;
00026       csminus = (complex *) malloc( sizeof(complex) * idim ) ;
00027       if( csplus == NULL || csminus == NULL ){
00028          fprintf(stderr,"\n*** fft cannot malloc space! ***\n"); exit(-1) ;
00029       }
00030    }
00031    nold = n = idim ;
00032 
00033    f1 = 1.0 ;  /* csplus init */
00034    m  = 1;
00035    k  = 0;
00036    while (n > m) {
00037       i3 = m << 1;
00038       f2 = m;
00039       al = f1*PI/f2;
00040       co = cos(al);
00041       si = sin(al);
00042       (csplus + k)->r = 1.;
00043       (csplus + k)->i = 0.;
00044          for (i0=0; i0 < m; i0++) {
00045          k++;
00046          csp = csplus + k;
00047          r0 = csp - 1;
00048          csp->r = r0->r * co - r0->i * si;
00049          csp->i = r0->i * co + r0->r * si;
00050       }
00051       m = i3;
00052    }
00053 
00054    f1 = -1.0 ;  /* csminus init */
00055    m  = 1;
00056    k  = 0;
00057    while (n > m) {
00058       i3 = m << 1;
00059       f2 = m;
00060       al = f1*PI/f2;
00061       co = cos(al);
00062       si = sin(al);
00063       (csminus + k)->r = 1.;
00064       (csminus + k)->i = 0.;
00065          for (i0=0; i0 < m; i0++) {
00066          k++;
00067          csp = csminus + k;
00068          r0  = csp - 1;
00069          csp->r = r0->r * co - r0->i * si;
00070          csp->i = r0->i * co + r0->r * si;
00071       }
00072       m = i3;
00073    }
00074    return ;
00075 }
00076 
00077 /*--------------------------------------------------------------------
00078    Complex-to-complex FFT in place:
00079      mode = -1 or +1 (NO SCALING ON INVERSE!)
00080      idim = dimension (power of 2)
00081      xc   = input/output array
00082    Automagically re-initializes itself when idim changes from
00083    previous call.  By AJJesmanowicz, modified by RWCox.
00084 ----------------------------------------------------------------------*/
00085 
00086 void csfft( int mode , int idim , complex * xc )
00087 {
00088    register unsigned int  m, n, i0, i1, i2, i3, k;
00089    register complex       *r0, *r1, *csp;
00090    register float         co, si, f0, f1, f2, f3, f4;
00091    double                 al;
00092 
00093    /** perhaps initialize **/
00094 
00095    if( nold != idim ) csfft_trigconsts( idim ) ;
00096 
00097    /** Main loop starts here **/
00098 
00099    n   = idim;
00100    i2  = idim >> 1;
00101    i1  = 0;
00102    csp = (mode > 0) ? csplus : csminus ;  /* choose const array */
00103 
00104    for (i0=0; i0 < n; i0 ++) {
00105       if ( i1 > i0 ) {
00106          r0    = xc + i0;
00107          r1    = xc + i1;
00108          f1    = r0->r;
00109          f2    = r0->i;
00110          r0->r = r1->r;
00111          r0->i = r1->i;
00112          r1->r = f1;
00113          r1->i = f2;
00114       }
00115       m = i2;
00116       while ( m && !(i1 < m) ) {
00117          i1 -= m;
00118          m >>= 1;
00119      }
00120      i1 += m;
00121    }
00122 
00123    m = 1;
00124    k = 0;
00125    while (n > m) {
00126       i3 = m << 1;
00127       for (i0=0; i0 < m; i0 ++) {
00128          co = (csp + k)->r;
00129          si = (csp + k)->i;
00130          for (i1=i0; i1 < n; i1 += i3) {
00131             r0    = xc + i1;
00132             r1    = r0 + m;
00133             f1    = r1->r * co;
00134             f2    = r1->i * si;
00135             f3    = r1->r * si;
00136             f4    = r1->i * co;
00137             f1   -= f2;
00138             f3   += f4;
00139             f2    = r0->r;
00140             f4    = r0->i;
00141             r1->r = f2 - f1;
00142             r1->i = f4 - f3;
00143             r0->r = f2 + f1;
00144             r0->i = f4 + f3;
00145          }
00146          k++;
00147       }
00148       m = i3;
00149    }
00150 
00151 #ifdef SCALE_INVERSE
00152    if (mode > 0) {
00153       f0 = 1.0 / idim ;
00154       i0 = 0;
00155       i1 = 1;
00156       while (i0 < n) {
00157          r0 = xc + i0;
00158          r1 = xc + i1;
00159          f1 = r0->r;
00160          f2 = r0->i;
00161          f3 = r1->r;
00162          f4 = r1->i;
00163          f1 *= f0;
00164          f2 *= f0;
00165          f3 *= f0;
00166          f4 *= f0;
00167          r0->r = f1;
00168          r0->i = f2;
00169          r1->r = f3;
00170          r1->i = f4;
00171          i0 += 2;
00172          i1 += 2;
00173       }
00174    }
00175 #endif
00176 
00177    return ;
00178 }