IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 16107


Ignore:
Timestamp:
Jan 17, 2008, 8:41:08 AM (18 years ago)
Author:
eugene
Message:

converting fft and fftn to code from Jones, with Magnier version of multi-dim

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

Legend:

Unmodified
Added
Removed

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

    r16099 r16107  
    4444$(SRC)/fft1d.old.$(ARCH).o              \
    4545$(SRC)/fft1d.$(ARCH).o          \
     46$(SRC)/fft2d.old.$(ARCH).o              \
    4647$(SRC)/fft2d.$(ARCH).o          \
    4748$(SRC)/fit2d.$(ARCH).o          \

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

    r16099 r16107  
    33int fft1d (int argc, char **argv) {
    44 
    5   int Npix, Nbit, ZeroImaginary;
     5  int N, Npix, Nbit, ZeroImaginary, forward;
    66  Vector *Ire, *Iim, *Ore, *Oim;
     7
     8  forward = TRUE;
     9  if ((N = get_argument (argc, argv, "-inverse"))) {
     10    remove_argument (N, &argc, argv);
     11    forward = FALSE;
     12  }
    713
    814  if ((argc != 6) || (strcmp (argv[3], "to"))) {
    … …  
    2228  Npix = Ire[0].Nelements;
    2329  if (!ZeroImaginary && (Npix != Iim[0].Nelements)) {
    24     gprint (GP_ERR, "vector mismatch in size\n");
     30    gprint (GP_ERR, "vector size mismatch\n");
    2531    return (FALSE);
    2632  }
    … …  
    5056  }   
    5157
    52   fft (Ore[0].elements, Oim[0].elements, Npix, Nbit, TRUE);
     58  fft1D (Ore[0].elements, Oim[0].elements, Npix, Nbit, forward);
    5359 
    5460  return (TRUE);

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

    r16099 r16107  
    33int fft2d (int argc, char **argv) {
    44 
    5   int i, N, Nx, Ny, Naxis[2];
    6   int Npix, ZeroImaginary, isign;
    7   float *t1, *t2, *out, *temp;
     5  int ZeroImaginary, forward;
     6  int N, Nx, Ny;
    87  Buffer *Ire, *Iim, *Ore, *Oim;;
    98
    10   isign = 1;
     9  forward = TRUE;
    1110  if ((N = get_argument (argc, argv, "-inverse"))) {
    1211    remove_argument (N, &argc, argv);
    13     isign = -1;
     12    forward = FALSE;
    1413  }
    1514
    … …  
    2827  }   
    2928  if ((Ire = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     29
     30  Nx = Ire[0].header.Naxis[0];
     31  Ny = Ire[0].header.Naxis[1];
     32
     33  // check input image dimensions (match lengths? binary lengths?)
     34  if (!ZeroImaginary) {
     35    if ((Nx != Iim[0].header.Naxis[0]) ||
     36        (Ny != Iim[0].header.Naxis[1])) {
     37      gprint (GP_ERR, "image size mismatch\n");
     38      return (FALSE);
     39    }
     40  }
     41  if (!IsBinary (Nx, NULL)) {
     42    gprint (GP_ERR, "Nx is not a binary number!\n");
     43    return (FALSE);
     44  }
     45  if (!IsBinary (Ny, NULL)) {
     46    gprint (GP_ERR, "Ny is not a binary number!\n");
     47    return (FALSE);
     48  }
     49
    3050  if ((Ore = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    3151  if ((Oim = SelectBuffer (argv[5], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    … …  
    3656  gfits_free_matrix (&Oim[0].matrix);
    3757  gfits_free_header (&Oim[0].header);
     58 
     59  /* fix up output headers (real) & allocate data buffer */
     60  CreateBuffer (Ore, Nx, Ny, -32, 0.0, 1.0);
     61  CreateBuffer (Oim, Nx, Ny, -32, 0.0, 1.0);
    3862
    39   /* get image dimensions, check value */
    40   Npix = Ire[0].header.Naxis[0]*Ire[0].header.Naxis[1];
    41   Nx = Ire[0].header.Naxis[0];
    42   Ny = Ire[0].header.Naxis[1];
    43   Naxis[0] = Ny; Naxis[1] = Nx;
    44   if (!IsBinaryOld (Npix)) {
    45     gprint (GP_ERR, "dimensions are not binary!\n");
    46     return (FALSE);
     63  gfits_copy_header (&Ire[0].header, &Ore[0].header);
     64  gfits_copy_header (&Ire[0].header, &Oim[0].header);
     65
     66  // copy data to output buffers (fft is done in place)
     67  memcpy (Ore[0].matrix.buffer, Ire[0].matrix.buffer, Nx*Ny*sizeof(float));
     68
     69  if (ZeroImaginary) {
     70    memset (Oim[0].matrix.buffer, 0, Nx*Ny*sizeof(float));
     71  } else {
     72    memcpy (Oim[0].matrix.buffer, Iim[0].matrix.buffer, Nx*Ny*sizeof(float));
    4773  }
    48  
    49   /* create working space */
    50   t1 = (float *) Ire[0].matrix.buffer;
    51   ALLOCATE (temp, float, 2*Npix);
    52   out = temp;
    5374
    54   /* copy input to working space */
    55   if (ZeroImaginary) {
    56     for (i = 0; i < Npix; i++, t1++) {
    57       *out = *t1;
    58       out++;
    59       *out = 0;
    60       out++;
    61     }
    62   } else {
    63     t2 = (float *) Iim[0].matrix.buffer;
    64     for (i = 0; i < Npix; i++, t1++, t2++) {
    65       *out = *t1;
    66       out++;
    67       *out = *t2;
    68       out++;
    69     }
    70   }
    71    
    7275  /* run the fft */
    73   fftN (temp, Naxis, 2, isign);
    74 
    75   /* fix up output headers (real) */
    76   gfits_copy_header (&Ire[0].header, &Ore[0].header);
    77   gfits_modify (&Ore[0].header, "NAXIS1", "%d", 1, Nx);
    78   gfits_modify (&Ore[0].header, "NAXIS2", "%d", 1, Ny);
    79   Ore[0].header.Naxis[0] = Nx;
    80   Ore[0].header.Naxis[1] = Ny;
    81   Ore[0].bitpix = Ire[0].bitpix;
    82   Ore[0].unsign = Ire[0].unsign;
    83   Ore[0].bscale = Ire[0].bscale;
    84   Ore[0].bzero  = Ire[0].bzero;
    85   gfits_create_matrix (&Ore[0].header, &Ore[0].matrix);
    86 
    87   /* fix up output headers (imaginary) */
    88   gfits_copy_header (&Ire[0].header, &Oim[0].header);
    89   gfits_modify (&Oim[0].header, "NAXIS1", "%d", 1, Nx);
    90   gfits_modify (&Oim[0].header, "NAXIS2", "%d", 1, Ny);
    91   Oim[0].header.Naxis[0] = Nx;
    92   Oim[0].header.Naxis[1] = Ny;
    93   Oim[0].bitpix = Ire[0].bitpix;
    94   Oim[0].unsign = Ire[0].unsign;
    95   Oim[0].bscale = Ire[0].bscale;
    96   Oim[0].bzero  = Ire[0].bzero;
    97   gfits_create_matrix (&Oim[0].header, &Oim[0].matrix);
    98 
    99   /* move data from working space to output buffers */
    100   out = temp;
    101   t1 = (float *) Ore[0].matrix.buffer;
    102   t2 = (float *) Oim[0].matrix.buffer;
    103   for (i = 0; i < Npix; i++, t1++, t2++) {
    104     *t1 = *out / Npix;
    105     out ++;
    106     *t2 = *out / Npix;
    107     out ++;
    108   }   
    109 
    110   free (temp);
     76  fftND ((float *)Ore[0].matrix.buffer, (float *)Oim[0].matrix.buffer, 2, Ire[0].header.Naxis, forward);
    11177
    11278  return (TRUE);

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

    r16099 r16107  
    2929int fft1d            PROTO((int, char **));
    3030int fft2d            PROTO((int, char **));
     31int fft2dold         PROTO((int, char **));
    3132int fit2d            PROTO((int, char **));
    3233int fit              PROTO((int, char **));
    … …  
    149150  {"fft1d",        fft1d,            "fft on the pixel-stream in an image"},
    150151  {"fft2d",        fft2d,            "fft on an image"},
     152  {"fft2dold",     fft2dold,         "fft on an image"},
    151153  {"fit",          fit,              "fit polynomial to vector pair"},
    152154  {"fit2d",        fit2d,            "fit 2-d polynomial to vector triplet"},

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

    r16099 r16107  
    8484
    8585/* in fft.c */
    86 void fft (float *dataRe, float *dataIm, int N, int Nbit, int forward);
    87 void fftN (float *data, int *nn, int ndim, int isign);
     86void fft1D (float *dataRe, float *dataIm, int N, int Nbit, int forward);
     87int fftND (float *dataRe, float *dataIm, int Ndim, int *Nsize, int forward);
    8888int IsBinary (int N, int *Nbit);
    8989

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

    r16099 r16107  
    22
    33// fft based on code by Douglas L. Jones (see note at EOF). modified for Ohana C style
    4 void fftN (float *data, int *nn, int ndim, int isign) {
    5  
    6   fprintf (stderr, "broken\n");
    7   return;
    8 }
    9 
    10 void fft (float *x, float *y, int n, int Nbit, int forward) {
     4void fft1D (float *x, float *y, int n, int Nbit, int forward) {
    115
    126  int i,j,k,n1,n2;
    … …  
    7569}                         
    7670
     71// This implementation uses the 1-D fft above for each of the vectors in each dimension.
     72// This requires 2(Nx*Ny*...) mem copies, but the fft operations are likely to happen in
     73// cache.
     74int fftND (float *x, float *y, int Ndim, int *Nsize, int forward) {
     75
     76  int i, nIndex, minor, major, iDim;
     77  int step, Nmajor, Nminor, Nmax, Ntotal;
     78  int *Nbit;
     79  float *tmpX, *tmpY;
     80
     81  ALLOCATE (Nbit, int, Ndim);
     82
     83  // find the longest axis and allocate storage for that length
     84  Nmax = 0;
     85  Ntotal = 1;
     86  for (i = 0; i < Ndim; i++) {
     87    Nmax = MAX(Nmax, Nsize[i]);
     88    Ntotal *= Nsize[i];
     89    if (!IsBinary (Nsize[i], &Nbit[i])) {
     90      free (Nbit);
     91      return (FALSE);
     92    }
     93  }
     94  ALLOCATE (tmpX, float, Nmax);
     95  ALLOCATE (tmpY, float, Nmax);
     96 
     97  step = 1;
     98  Nminor = 1;
     99  Nmajor = Ntotal;
     100  for (iDim = 0; iDim < Ndim; iDim++) {
     101    step *= Nsize[iDim];
     102    Nmajor /= Nsize[iDim];
     103
     104    // we perform the FFT along all other dimensions
     105    for (major = 0; major < Nmajor; major++) {
     106      for (minor = 0; minor < Nminor; minor++) {
     107        // nIndex = minor + i*Nminor + major*step;
     108        // extract the data values to the temp vector
     109        nIndex = minor + major*step;
     110        for (i = 0; i < Nsize[iDim]; i++) {
     111          tmpX[i] = x[nIndex];
     112          tmpY[i] = y[nIndex];
     113          nIndex += Nminor;
     114        }
     115
     116        fft1D (tmpX, tmpY, Nsize[iDim], Nbit[iDim], forward);
     117
     118        // replace the result vectors
     119        nIndex = minor + major*step;
     120        for (i = 0; i < Nsize[iDim]; i++) {
     121          x[nIndex] = tmpX[i];
     122          y[nIndex] = tmpY[i];
     123          nIndex += Nminor;
     124        }
     125      }
     126    }
     127    Nminor *= Nsize[iDim];
     128  }
     129  free (tmpX);
     130  free (tmpY);
     131  return (TRUE);
     132}
     133
    77134// check that a number is binary (2^Nbit).  returns int(log_2(N)) in Nbit
    78135int IsBinary (int N, int *Nbit) {
    … …  
    80137  int i, Nset;
    81138
    82   *Nbit = 0;
     139  if (Nbit != NULL) *Nbit = 0;
    83140  Nset = 0;
    84141  for (i = 0; i < 8*sizeof(int); i++) {
    85142    if (N & 0x01) {
    86143      Nset ++;
    87       *Nbit = i;
     144      if (Nbit != NULL) *Nbit = i;
    88145    }
    89146    N >>= 1;
    90147  }
    91148  if (Nset > 1) return (FALSE);
     149  return (TRUE); 
    92150}
    93151

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

    r16099 r16107  
    4949    mmax = istep;
    5050  }
    51 }
    52 
    53 void fftold_Nindices (float *Data, int N, int isign) {
    54 
    55   int n,ifp1,m,j,ipf2,i;
    56   double wtemp,wr,wpr,wpi,wi,theta;
    57   float tempr,tempi, *data;
    58 
    59   data = Data;
    60   n = N << 1;
    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);
    82     wtemp = sin(0.5*theta);
    83     wpr = -2.0*wtemp*wtemp;
    84     wpi = sin(theta);
    85     wr = 1.0;
    86     wi = 0.0;
    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;
    98       }
    99       wr = (wtemp = wr)*wpr - wi*wpi+wr;
    100       wi = wi*wpr + wtemp*wpi + wi;
    101     }
    102     ifp1 = ifp2;
    103   }
    104 
    10551}
    10652
Note: See TracChangeset for help on using the changeset viewer.