IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16099


Ignore:
Timestamp:
Jan 16, 2008, 1:10:16 PM (18 years ago)
Author:
eugene
Message:

working on the fft function

Location:
trunk/Ohana/src/opihi
Files:
2 added
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.data/Makefile

    r16013 r16099  
    4242$(SRC)/erase.$(ARCH).o          \
    4343$(SRC)/extract.$(ARCH).o        \
     44$(SRC)/fft1d.old.$(ARCH).o              \
    4445$(SRC)/fft1d.$(ARCH).o          \
    4546$(SRC)/fft2d.$(ARCH).o          \

  • trunk/Ohana/src/opihi/cmd.data/fft1d.c

    r16094 r16099  
    5050  }   
    5151
    52   fft (Ore[0].elements, Oim[0].elements, Npix, Nbit);
     52  fft (Ore[0].elements, Oim[0].elements, Npix, Nbit, TRUE);
    5353 
    5454  return (TRUE);

  • trunk/Ohana/src/opihi/cmd.data/fft2d.c

    r7917 r16099  
    4242  Ny = Ire[0].header.Naxis[1];
    4343  Naxis[0] = Ny; Naxis[1] = Nx;
    44   if (!IsBinary (Npix)) {
     44  if (!IsBinaryOld (Npix)) {
    4545    gprint (GP_ERR, "dimensions are not binary!\n");
    4646    return (FALSE);

  • trunk/Ohana/src/opihi/cmd.data/init.c

    r16013 r16099  
    2626int erase            PROTO((int, char **));
    2727int extract          PROTO((int, char **));
     28int fft1dold         PROTO((int, char **));
    2829int fft1d            PROTO((int, char **));
    2930int fft2d            PROTO((int, char **));
    … …  
    145146  {"erase",        erase,            "erase objects on an overlay"},
    146147  {"extract",      extract,          "extract a portion of a buffer into another buffer"},
     148  {"fft1dold",     fft1dold,         "fft on the pixel-stream in an image"},
    147149  {"fft1d",        fft1d,            "fft on the pixel-stream in an image"},
    148150  {"fft2d",        fft2d,            "fft on an image"},

  • trunk/Ohana/src/opihi/include/data.h

    r16094 r16099  
    8484
    8585/* in fft.c */
    86 void fft (float *dataRe, float *dataIm, int N, int Nbit);
     86void fft (float *dataRe, float *dataIm, int N, int Nbit, int forward);
    8787void fftN (float *data, int *nn, int ndim, int isign);
    88 int IsBinary (int N);
    89 void fourn (float *data, int *nn, int ndim, int isign);
     88int IsBinary (int N, int *Nbit);
     89
     90// fftold.c
     91void fftold (float *Data, int N, int isign);
     92void fftNold (float *data, int *nn, int ndim, int isign);
     93int IsBinaryOld (int N);
    9094
    9195/* in gaussj.c */

  • trunk/Ohana/src/opihi/lib.data/Makefile

    r16040 r16099  
    1919$(SDIR)/page.$(ARCH).o                  \
    2020$(SDIR)/fft.$(ARCH).o                   \
     21$(SDIR)/fftold.$(ARCH).o                        \
    2122$(SDIR)/svdcmp.$(ARCH).o                \
    2223$(SDIR)/convert.$(ARCH).o               \

  • trunk/Ohana/src/opihi/lib.data/fft.c

    r16094 r16099  
    88}
    99
    10 void fft (double *x, double *y, int n, int Nbit) {
     10void fft (float *x, float *y, int n, int Nbit, int forward) {
    1111
    1212  int i,j,k,n1,n2;
    13   double c,s,e,a,t1,t2;       
     13  float c,s,e,a,t1,t2;       
     14  double factor;
    1415         
    1516  // bit-reverse
    … …  
    1920    n1 = n2;
    2021    while ( j >= n1 ) {
    21       j = j - n1;
    22       n1 = n1/2;
     22      j -= n1;
     23      n1 /= 2;
    2324    }
    24     j = j + n1;
     25    j += n1;
    2526               
    2627    if (i < j) {
    … …  
    3334    }
    3435  }
    35  
    3636                                         
    3737  n1 = 0; /* FFT */
    3838  n2 = 1;
    3939                                             
     40  if (forward) {
     41    factor = +2.0*M_PI;
     42  } else {
     43    factor = -2.0*M_PI;
     44  }
     45
    4046  for (i=0; i < Nbit; i++) {
    4147    n1 = n2;
    4248    n2 = n2 + n2;
    43     e = -6.283185307179586/n2;
     49    e = factor/n2;
    4450    a = 0.0;
    4551                                             
    … …  
    6066  }
    6167                                     
     68  // re-normalize
     69  for (i = 0; i < n; i++) {
     70    x[i] /= n;
     71    y[i] /= n;
     72  }
     73
    6274  return;
    6375}                         
    6476
    65 int IsBinary (int N, int *nbin) {
     77// check that a number is binary (2^Nbit).  returns int(log_2(N)) in Nbit
     78int IsBinary (int N, int *Nbit) {
    6679
    67   // check that a number is binary (2^nbit)
     80  int i, Nset;
    6881
    69   Nbit = 0;
     82  *Nbit = 0;
    7083  Nset = 0;
    7184  for (i = 0; i < 8*sizeof(int); i++) {
    7285    if (N & 0x01) {
    7386      Nset ++;
    74       Nbin = i;
     87      *Nbit = i;
    7588    }
    7689    N >>= 1;
    7790  }
    78   *nbin = Nbit;
    7991  if (Nset > 1) return (FALSE);
    8092}
    … …  
    92104/* m: n = 2**m                                            */
    93105/*   input/output                                         */
    94 /* x: double array of length n with real part of data     */
    95 /* y: double array of length n with imag part of data     */
     106/* x: float array of length n with real part of data     */
     107/* y: float array of length n with imag part of data     */
    96108/*                                                        */
    97109/*   Permission to copy and use this program is granted   */

  • trunk/Ohana/src/opihi/lib.data/fftold.c

    r16094 r16099  
    33#define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
    44
    5 void fft (float *Data, int N, int isign) {
     5void fftold (float *Data, int N, int isign) {
    66
    77  int n,mmax,m,j,istep,i;
    … …  
    5151}
    5252
    53 void fftold (float *Data, int N, int isign) {
    54 
    55   int n,mmax,m,j,istep,i;
     53void fftold_Nindices (float *Data, int N, int isign) {
     54
     55  int n,ifp1,m,j,ipf2,i;
    5656  double wtemp,wr,wpr,wpi,wi,theta;
    5757  float tempr,tempi, *data;
    5858
    59   data = Data - 1;
     59  data = Data;
    6060  n = N << 1;
    61   j = 1;
    62 
    63   for (i = 1; i < n; i+=2) {
    64     if (j > i) {
    65       FSWAP (data[j], data[i]);
    66       FSWAP (data[j+1], data[i+1]);
    67     }
    68     m = n >> 1;
    69     while (m >= 2 && j > m) {
    70       j -= m;
    71       m >>= 1;
    72     }
    73     j += m;
    74   }
    75   mmax = 2;
    76   while (n > mmax) {
    77     istep = 2*mmax;
    78     theta = 6.28318530717959 / (isign*mmax);
     61
     62  ip1 = 2;
     63  ip2 = n;
     64  i2rev = 0;
     65  for (i2 = 0; i2 < ip2; i2 += ip1) {
     66    if (i2 < i2ref) {
     67      FSWAP (data[i2rev], data[i2]);
     68      FSWAP (data[i2rev+1], data[i2+1]);
     69    }
     70    ibit = ip2 >> 1;
     71    while (ibit >= ip1 && i2rev >= ibit) {
     72      i2rev -= ibit;
     73      ibit >>= 1;
     74    }
     75    i2rev += ibit;
     76  }
     77
     78  ifp1 = 2;
     79  while (ifp1 < ip2) {
     80    ifp2 = ifp1 << 1;
     81    theta = isign*6.28318530717959 / (2*ifp1/ip1);
    7982    wtemp = sin(0.5*theta);
    8083    wpr = -2.0*wtemp*wtemp;
    … …  
    8285    wr = 1.0;
    8386    wi = 0.0;
    84     for (m = 1; m < mmax; m+=2) {
    85       for (i = m; i <= n; i+=istep) {
    86         j = i + mmax;
    87         tempr = wr*data[j] - wi*data[j+1];
    88         tempi = wr*data[j+1] + wi*data[j];
    89         data[j] = data[i] - tempr;
    90         data[j+1] = data[i+1] - tempi;
    91         data[i] += tempr;
    92         data[i+1] += tempi;
     87    for (i3 = 0; i3 < ifp1; i3+=ip1) {
     88      // in the N-d version, this loop is broken into two loops
     89      for (i1 = i3; i1 < n; i1+=ipf2) {
     90        k1 = i1;
     91        k2 = k1 + ifp1;
     92        tempr = wr*data[k2] - wi*data[k2+1];
     93        tempi = wr*data[k2+1] + wi*data[k2];
     94        data[k2] = data[k1] - tempr;
     95        data[k2+1] = data[k1+1] - tempi;
     96        data[k1] += tempr;
     97        data[k1+1] += tempi;
    9398      }
    9499      wr = (wtemp = wr)*wpr - wi*wpi+wr;
    95100      wi = wi*wpr + wtemp*wpi + wi;
    96101    }
    97     mmax = istep;
    98   }
     102    ifp1 = ifp2;
     103  }
     104
    99105}
    100106
    101107/* convert indices to zero reference */
    102 void fftN (float *data, int *nn, int ndim, int isign) {
     108void fftNold (float *data, int *nn, int ndim, int isign) {
    103109
    104110  int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
    … …  
    138144    while (ifp1 < ip2) {
    139145      ifp2 = ifp1 << 1;
    140       theta = isign*6.28318530717959/(ifp2/ip1);
     146      theta = isign*6.28318530717959/(2*ifp1/ip1);
    141147      wtemp = sin(0.5*theta);
    142148      wpr = -2.0*wtemp*wtemp;
    … …  
    182188*/
    183189
    184 int IsBinary (int N) {
     190int IsBinaryOld (int N) {
    185191
    186192  int i, nbit;
    … …  
    199205
    200206} 
    201 
    202 #define FSWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
    203 
    204 void fourn (float *data, int *nn, int ndim, int isign) {
    205 
    206   int i1,i2,i3,i2rev,i3rev,ip1,ip2,ip3,ifp1,ifp2;
    207   int ibit,idim,k1,k2,n,nprev,nrem,ntot;
    208   float tempi,tempr;
    209   double theta,wi,wpi,wpr,wr,wtemp;
    210 
    211   ntot = 1;
    212   for (idim = 1; idim <= ndim; idim++)
    213     ntot *= nn[idim];
    214   nprev = 1;
    215   for (idim = ndim; idim >= 1; idim--) {
    216     n = nn[idim];
    217     nrem = ntot/(n*nprev);
    218     ip1 = nprev << 1;
    219     ip2 = ip1*n;
    220     ip3 = ip2*nrem;
    221     i2rev = 1;
    222     for (i2 = 1; i2 <= ip2; i2+=ip1) {
    223       if (i2 < i2rev) {
    224         for (i1 = i2; i1 <= i2+ip1-2; i1+=2) {
    225           for (i3 = i1; i3 <= ip3;i3+=ip2) {
    226             i3rev = i2rev+i3-i2;
    227             FSWAP(data[i3],data[i3rev]);
    228             FSWAP(data[i3+1],data[i3rev+1]);
    229           }
    230         }
    231       }
    232       ibit = ip2 >> 1;
    233       while (ibit >= ip1 && i2rev > ibit) {
    234         i2rev -=  ibit;
    235         ibit >>= 1;
    236       }
    237       i2rev += ibit;
    238     }
    239     ifp1 = ip1;
    240     while (ifp1 < ip2) {
    241       ifp2 = ifp1 << 1;
    242       theta = isign*6.28318530717959/(ifp2/ip1);
    243       wtemp = sin(0.5*theta);
    244       wpr  =  -2.0*wtemp*wtemp;
    245       wpi = sin(theta);
    246       wr = 1.0;
    247       wi = 0.0;
    248       for (i3 = 1;i3<=ifp1;i3+=ip1) {
    249         for (i1 = i3;i1<=i3+ip1-2;i1+=2) {
    250           for (i2 = i1;i2<=ip3;i2+=ifp2) {
    251             k1 = i2;
    252             k2 = k1+ifp1;
    253             tempr = wr*data[k2]-wi*data[k2+1];
    254             tempi = wr*data[k2+1]+wi*data[k2];
    255             data[k2] = data[k1]-tempr;
    256             data[k2+1] = data[k1+1]-tempi;
    257             data[k1] += tempr;
    258             data[k1+1] += tempi;
    259           }
    260         }
    261         wr = (wtemp = wr)*wpr-wi*wpi+wr;
    262         wi = wi*wpr+wtemp*wpi+wi;
    263       }
    264       ifp1 = ifp2;
    265     }
    266     nprev *= n;
    267   }
    268 }
    269 
    270 #undef FSWAP
Note: See TracChangeset for help on using the changeset viewer.