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 <math.h>
00008 #include <stdlib.h>
00009 #include <stdio.h>
00010 
00011 #ifndef PI
00012 #  define PI 3.1415926536
00013 #endif
00014 
00015 /*
00016    mode  : -1 forward, 1 inverse.
00017    idim  : FFT dimension, max 1024.
00018    xr    : real array.
00019    xi    : imaginary part of the array.
00020 */
00021 
00022 /* --------------------------------- */
00023   void cfft( int mode , int idim , float *xr , float *xi )
00024 /* --------------------------------- */
00025 {
00026 #define IDMAX  1024
00027 #define LIDMAX 10
00028    static int idold = -999 ;
00029    register int i0,i1,i2,i4,i5,m0, id = idim;
00030    float    f1,f3,f4,f5,al,co,si,md = mode;
00031    static   int    n,m[LIDMAX];
00032    static   float  f2,c[IDMAX/2],s[IDMAX/2];
00033 
00034    /**** preparations if id has changed since last call ****/
00035 
00036    if (idold != id) {
00037      idold = id ;
00038      /* check for power of 2 */
00039      for( i4=4 ; i4 <= IDMAX ; i4 *= 2 ){
00040         if( id == i4 )break ;
00041      }
00042      if( id != i4 ){
00043        fprintf(stderr,"\n In cfft : illegal idim=%d\n",idim);
00044        exit(1) ;
00045      }
00046      f2     = id;
00047      n      = log(f2)/log(2.) + .5;
00048      m[n-1] = 1;
00049      al     = 2.*PI/f2;
00050      co     = cos(al);
00051      si     = sin(al);
00052      c[0]   = 1.;
00053      s[0]   = 0.;
00054      for(i4 = 1; i4 < 512; i4++) {
00055        c[i4] = c[i4-1]*co - s[i4-1]*si;
00056        s[i4] = s[i4-1]*co + c[i4-1]*si;
00057      }
00058      for(i4 = n-1; i4 >= 1; i4--) m[i4-1] = m[i4]*2;
00059    }
00060 
00061    /**** Main loop starts here ****/
00062 
00063    for(i0 = 0; i0 < n; i0++) {
00064      i1 = 0;
00065      m0 = m[i0];
00066      for(i2 = 0; i2 < m[n-i0-1]; i2++) {
00067        f4 = c[i1];
00068        f5 = s[i1]*md;
00069          for(i4 = 2*m0*i2; i4 < m0*(2*i2+1); i4++) {
00070            f3 = xr[i4+m0]*f4 - xi[i4+m0]*f5;
00071            f1 = xi[i4+m0]*f4 + xr[i4+m0]*f5;
00072            xr[i4+m0] = xr[i4] - f3;
00073            xr[i4] += f3;
00074            xi[i4+m0] = xi[i4] - f1;
00075            xi[i4] += f1;
00076          }
00077         for(i4 = 1; i4 < n; i4++) {
00078            i5 = i4;
00079            if (i1 < m[i4]) goto i1_plus1;
00080            i1 -= m[i4];
00081          }
00082 i1_plus1: i1 += m[i5];
00083      }
00084    }
00085    i1 = 0;
00086    for(i4 = 0; i4 < id; i4++) {
00087      if (i1 > i4) {
00088        f3 = xr[i4];
00089        f1 = xi[i4];
00090        xr[i4] = xr[i1];
00091        xi[i4] = xi[i1];
00092        xr[i1] = f3;
00093        xi[i1] = f1;
00094      }
00095      for(i2 = 0; i2 < n; i2++) {
00096        i5 = i2;
00097        if (i1 < m[i2]) goto i1_plus2;
00098        i1 -= m[i2];
00099      }
00100 i1_plus2:   i1 += m[i5];
00101    }
00102 
00103    if (md > 0.) {
00104      f1 = 1./f2;
00105      for(i4 = 0; i4 < id; i4++) {
00106        xr[i4] *= f1;
00107        xi[i4] *= f1;
00108      }
00109    }
00110    return;
00111 }
00112 
00113 /*********************************************************************/
00114 
00115 /*----------------------------------*/
00116 void cfft2d( int mode , int nx,int ny , float *xr, float *xi )
00117 /*----------------------------------*/
00118 {
00119    float *rbuf , *ibuf ;
00120    register int ii , jj , jbase ;
00121 
00122    rbuf = (float *)malloc( ny * sizeof(float) ) ;
00123    ibuf = (float *)malloc( ny * sizeof(float) ) ;
00124    if( rbuf == NULL || ibuf == NULL ){
00125       fprintf(stderr,"malloc error in cfft2d\n") ;
00126       exit(1) ;
00127    }
00128 
00129    for( jj=0 ; jj < ny ; jj++ ){
00130       jbase = nx * jj ;
00131       cfft( mode , nx , &xr[jbase] , &xi[jbase] ) ;
00132    }
00133 
00134    for( ii=0 ; ii < nx ; ii++ ){
00135       for( jj=0 ; jj < ny ; jj++ ){
00136          rbuf[jj] = xr[ii + jj*nx] ;
00137          ibuf[jj] = xi[ii + jj*nx] ;
00138       }
00139       cfft( mode , ny , rbuf , ibuf ) ;
00140       for( jj=0 ; jj < ny ; jj++ ){
00141          xr[ii+jj*nx] = rbuf[jj] ;
00142          xi[ii+jj*nx] = ibuf[jj] ;
00143       }
00144    }
00145 
00146    free(rbuf) ; free(ibuf) ;
00147    return ;
00148 }