/* This code provided as is */

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <limits.h>

/* Add #define RGBA if your are on a big-endian machine (ie non-Intel) */

#define RGBA
#define FILENAME_SUFFIX 8
#define PI              3.14159265359
#define MFALSE          0
#define MTRUE           !MFALSE

unsigned char iiiInterpolateGrayFromXY(float xoriginal, float yoriginal,
                                       unsigned char *grayimage, int w, int h);
unsigned char *iiRotateGrayImage(int w, int h, unsigned char *grayimage,
                                 int *newwidth, int *newheight,
                                 float degreestorotate, float scale);
unsigned char *iiAverageGrayImages (int w, int h, unsigned char **grayimage,
                                    int numberofimages);
unsigned char *iiConvertToGrayscale(unsigned int *rgbimage, int w, int h);
unsigned char *iiReadPGMFile (FILE *fp, int *w, int *h);
unsigned int *iiReadPPMFile (FILE *fp, int *w, int *h);
int iiWritePPMFile (char *filename, unsigned int *pixel, int width, int height);
int iiWritePGMFileFloat (char *filename, float *pixel, int width, int height);
int iiWritePGMFileInt (char *filename, unsigned char *pixel, int width,
                       int height);
int iiWriteHistogramFile (char *filename, unsigned char *pixel, int width,
                          int height);
void iiComputeHistogram (unsigned char *pixel, int width, int height,
                         unsigned int *histogram, unsigned int histsize);
void iiComputeLookup(unsigned int *histogram, float *lookup, int width,
                     int height, unsigned int histsize);
void iiEqualizeImage(char *filename, unsigned char *grayimage, int width,
                     int height);
int iiGetRed(unsigned int pixVal);
int iiGetGreen(unsigned int pixVal);
int iiGetBlue(unsigned int pixVal);
int iiWriteFourierResult(char *outrealimagfilename, char *outmagphasefilename,
                         float *realarray, float *complexarray,
                         int streamlength);
int iiFastFourier1D(float *realfunc, float *complexfunc, int numfunc,
                    float **realarray, float **complexarray, int inverse);
int iiRead1DFourierFile(char *filename, int *numfunc, float **func);
int iiFastFourier2D(float *inputrealimage, float *inputcompleximage, int width,
                    int height, float **realimage, float **compleximage,
                    int inverse);
int iiRemoveHighFrequencies(float *realimage, float *compleximage,
                            int width, int height);
unsigned char *iiReadOffsetLinear8BitAudioFile (FILE *fp, int *streamlength);

#if 0
unsigned int *floatToColor(float *pixels, int w, int h)
{
  int r,g,b;
  int x,y;
  unsigned int *newPix;
  
  newPix = new unsigned int[w * h];

  for (y = 0; y < h; y++)
  {
    for (x = 0; x < w; x++)
    {
      r  = g = b = (int)(pixels[y * w + x]);
      newPix[y * w + x] = 0;
      newPix[y * w + x] |= (r & 0xFF);
      newPix[y * w + x] |= (g & 0xFF) << 8;
      newPix[y * w + x] |= (b & 0xFF) << 16;
      /*      printf("%f %x\n", pixels[y * w + x], newPix[y * w + x]); */

    }
  }

  return newPix;
}
#endif

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiGetRed(unsigned int pixVal)
{
  return 0xFF & pixVal;
}

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiGetGreen(unsigned int pixVal)
{
  return 0xFF & (pixVal >> 8);
}

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiGetBlue(unsigned int pixVal)
{
  return 0xFF & (pixVal >> 16);
}

/* -----------------------------------------------------------------------------
NAME:     iiReadPGMFile
PURPOSE:  Reads a PGM File into an unsigned char array
PRE:      The fp should be opened
POST:     The fp should be closed
RETURN:   An array of unsigned chars
AUTHOR:   CS Dept BYU
DATE:
MODIFIED: CMS-September 20, 1999
----------------------------------------------------------------------------- */
unsigned char *iiReadPGMFile (FILE *fp, int *w, int *h)
{
  int            max, type, x, y;
  char           buffer[250];
  unsigned char *pixel=NULL;

  if (!fp)
    return(pixel);
  fscanf(fp,"P%d\n",&type);
  if (type != 5)
    return NULL;
  buffer[0] = '#';
  while(buffer[0] == '#')
    fgets(buffer,250,fp);
  sscanf(buffer,"%d%d", w, h);
  pixel = (unsigned char *)malloc((*w) * (*h) * sizeof(unsigned char));
  if (!pixel)
    return(pixel);
  fgets(buffer,250,fp);
  sscanf(buffer,"%d",&max);
  for (y=0; y<(*h); y++) {
    for (x=0; x<(*w); x++) {
      fscanf(fp,"%c",&(pixel[y*(*w)+x]));
    }
  }
  return pixel;
} /* iiReadPGMFile */

/* -----------------------------------------------------------------------------
NAME:     iiReadPPMFile
PURPOSE:  Reads a PPM file into an unsigned int array
PRE:      The fp should be opened
POST:     The fp should be closed
RETURN:   An array of 
AUTHOR:   CS Dept BYU
DATE:
MODIFIED: CMS-September 10, 1999
----------------------------------------------------------------------------- */
unsigned int *iiReadPPMFile (FILE *fp, int *w, int *h)
{
  int           max, type;
  char          buffer[250];
  unsigned int *pixel = NULL;
  int           rint, gint, bint, x, y;
  unsigned char ruchar, guchar, buchar;

  if (!fp)
    return(pixel);

  fscanf(fp,"P%d\n",&type);
  if (type != 5 && type != 6 && type != 2 && type != 3)
    return NULL;
  buffer[0] = '#';
  while(buffer[0] == '#')
    fgets(buffer,250,fp);
  sscanf(buffer, "%d%d", w, h);
  pixel = (unsigned int *)malloc((*w) * (*h) * sizeof(unsigned int));
  if (!pixel)
    return(pixel);
  fgets(buffer,250,fp);
  sscanf(buffer,"%d",&max);
  for (y=0; y<(*h); y++) {
    for (x=0; x<(*w); x++) {
      if (type == 2 || type == 3) {
        if (type == 3)
          fscanf(fp,"%d%d%d",&rint,&gint,&bint);
        else {
          fscanf(fp,"%d",&rint);
          gint=bint=rint;
        }
        pixel[y*(*w)+x] = 0;
#ifdef RGBA
        pixel[y*(*w)+x] |= (rint&0xFF) << 24;
        pixel[y*(*w)+x] |= (gint&0xFF) << 16;
        pixel[y*(*w)+x] |= (bint&0xFF) << 8;
        pixel[y*(*w)+x] |= 0xFF;
#else
        pixel[y*(*w)+x] |= (rint&0xFF);
        pixel[y*(*w)+x] |= (gint&0xFF) << 8;
        pixel[y*(*w)+x] |= (bint&0xFF) << 16;
        pixel[y*(*w)+x] |= 0xFF << 24;
#endif
      }
      else {
        if (type == 6)
          fscanf(fp,"%c%c%c",&ruchar,&guchar,&buchar);
        else {
          fscanf(fp,"%c",&ruchar);
          guchar=buchar=ruchar;
        }
        pixel[y*(*w)+x] = 0;
#ifdef RGBA
        pixel[y*(*w)+x] |= (ruchar&0xFF) << 24;
        pixel[y*(*w)+x] |= (guchar&0xFF) << 16;
        pixel[y*(*w)+x] |= (buchar&0xFF) << 8;
        pixel[y*(*w)+x] |= 0xFF;
#else
        pixel[y*(*w)+x] |= (ruchar&0xFF);
        pixel[y*(*w)+x] |= (guchar&0xFF) << 8;
        pixel[y*(*w)+x] |= (buchar&0xFF) << 16;
        pixel[y*(*w)+x] |= 0xFF << 24;
#endif
      }
    }
  }
  return pixel;
} /* iiReadPPMFile */

/* -----------------------------------------------------------------------------
NAME:     iiWritePPMFile
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiWritePPMFile (char *filename, unsigned int *pixel, int width, int height)
{
  FILE        *fp = NULL;
  int          x, y;
  unsigned int ind;

  /* NOTE All of your values better be <= 255! */

  if (!pixel)
    return 0;
  fp = fopen(filename,"wt");
  if (fp == NULL)
    return -1;
  fprintf(fp,"P6\n");
  fprintf(fp,"%d %d\n",width,height);
  fprintf(fp,"255\n");
  for (y=0; y < height; y++) {
    for (x=0; x<width; x++) {
      ind = y*width+x;
      fprintf(fp,"%c",(char)((pixel[ind])&0xFF));
      fprintf(fp,"%c",(char)((pixel[ind]>>8)&0xFF));
      fprintf(fp,"%c",(char)((pixel[ind]>>16)&0xFF));
    }
  }
  fclose(fp);
  return 0;
} /* iiWritePPMFile */

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiWritePGMFileFloat (char *filename, float *pixel, int width, int height)
{
  int          x, y;
  unsigned int ind;
  FILE        *fp = fopen(filename,"wt");

  /* NOTE All of your values better be <= 255! */

  if (!fp)
    return -1;
  fprintf(fp,"P5\n");
  fprintf(fp,"%d %d\n",width,height);
  fprintf(fp,"255\n");
  for (y=0; y < height; y++) {
    for (x=0; x< width; x++) {
      ind = y*width+x;
      fprintf(fp,"%c",(unsigned char)(pixel[ind]));
    }
  }
  fclose(fp);
  return 0;
} /* iiWritePGMFileFloat */

/* -----------------------------------------------------------------------------
NAME:
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiWritePGMFileInt (char *filename, unsigned char *pixel, int width,
                       int height)
{
  FILE        *fp = NULL;
  int          x, y;
  unsigned int ind;

  /* NOTE All of your values better be <= 255! */

  fp = fopen(filename,"wt");
  if (!fp)
    return(MFALSE);
  fprintf(fp,"P5\n");
  fprintf(fp,"%d %d\n",width,height);
  fprintf(fp,"255\n");
  for (y=0; y < height; y++) {
    for (x=0; x< width; x++) {
      ind = y*width+x;
      fprintf(fp,"%c",pixel[ind]);
    }
  }
  fclose(fp);
  return(MTRUE);
} /* iiWritePGMFileInt */

/* -----------------------------------------------------------------------------
NAME:     iiComputeHistogram
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
void iiComputeHistogram (unsigned char *pixel, int width, int height,
                         unsigned int *histogram, unsigned int histsize)
{
  unsigned int arraysize;
  unsigned int i;

  memset(histogram, 0, histsize*sizeof(unsigned int));

  arraysize = width*height;
  for (i=0;i<arraysize;i++) {
    if (pixel[i] < histsize)
      (histogram[pixel[i]])++;
  }
} /* iiComputeHistogram */
/* -----------------------------------------------------------------------------
NAME:     iiWriteHistogramFile
PURPOSE:
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int iiWriteHistogramFile (char *filename, unsigned char *pixel, int width,
                          int height)
{
  FILE        *fp = NULL;
  unsigned int arraysize;
  unsigned int i;
  unsigned int histogram[256];

  /* NOTE All of your values better be <= 255! */

  fp = fopen(filename,"w");
  if (!fp)
    return(MFALSE);

  iiComputeHistogram(pixel, width, height, histogram, 256);

  for (i=0;i<256;i++)
    fprintf(fp, "%d\n", histogram[i]);

  fclose(fp);
  return(MTRUE);
} /* iiWriteHistogramFile */

/* -----------------------------------------------------------------------------
NAME:     iiConvertToGrayscale
PURPOSE:  Converts an RGB image to grayscale.
PRE:      The image must be read into the image array.
POST:     The grayscale image should be freed.
RETURN:   The image in grayscale
AUTHOR:   CMS
DATE:     September 11, 1999
MODIFIED:
----------------------------------------------------------------------------- */
unsigned char *iiConvertToGrayscale(unsigned int *rgbimage, int w, int h)
{
  unsigned char *grayimage = NULL;
  unsigned int   ind;
  unsigned char  r, g, b;
  int            x, y;

  if (!rgbimage)
    return(NULL);

  grayimage = (unsigned char *)malloc(w*h*sizeof(unsigned char));
  if (!grayimage)
    return(NULL);

  for (y=0; y<h; y++) {
    for (x=0; x<w; x++) {
      ind = y*w+x;
      r = (unsigned char)(rgbimage[ind])&0xFF;
      g = (unsigned char)(rgbimage[ind]>>8)&0xFF;
      b = (unsigned char)(rgbimage[ind]>>16)&0xFF;
      grayimage[ind] = (unsigned char)((float)r*0.299+(float)g*0.587+
                                       (float)b*0.114);
      if (grayimage[ind] > 255)
        grayimage[ind] = 255;
    }
  }
  return(grayimage);
} /* iiConvertToGrayscale */
/* -----------------------------------------------------------------------------
NAME:     iiComputeLookup
PURPOSE:  Computes a lookup table for equalizing the giving histogram.
PRE:      The histogram must be computed.
POST:     None.
RETURN:   None.
AUTHOR:   CMS
DATE:     September 18, 1999
MODIFIED:
----------------------------------------------------------------------------- */
void iiComputeLookup(unsigned int *histogram, float *lookup, int width,
                     int height, unsigned int histsize)
{
  unsigned int i;
  float        cumprob;
  float        imagearea;

  cumprob = 0.0;
  memset(lookup, 0, histsize*sizeof(float));
  imagearea = (float)width*(float)height;
  for (i=0; i<histsize; i++) {
    cumprob += (float)histogram[i]/imagearea;
    lookup[i] = cumprob*((float)histsize-1.001);
  }
} /* iiComputeLookup */
/* -----------------------------------------------------------------------------
NAME:     iiEqualizeImage
PURPOSE:  Equalizes the histogram of the given image and modifies the image
          to represent this equalization.
PRE:      The image must exist.
POST:     None.
RETURN:   None.
AUTHOR:   CMS
DATE:     September 18, 1999
MODIFIED:
----------------------------------------------------------------------------- */
void iiEqualizeImage(char *filename, unsigned char *grayimage, int width,
                     int height)
{
  unsigned int      histogram[256];
  float             lookup[256];
  unsigned long int imagearea;
  unsigned long int i;
  FILE             *fp;

  fp = fopen(filename,"w");

  imagearea = width*height;
  iiComputeHistogram(grayimage, width, height, histogram, 256);
  iiComputeLookup(histogram, lookup, width, height, 256);
  if (fp) {
    for (i=0;i<256;i++)
      fprintf(fp, "%1.3lf\n", lookup[i]);
  }
  for (i=0;i<imagearea;i++)
    grayimage[i] = (unsigned char)lookup[grayimage[i]];

  if (fp)
    fclose(fp);
} /* iiEqualizeImage */
/* -----------------------------------------------------------------------------
NAME:     iiAverageGrayImages
PURPOSE:  Averages an array of grayscale images
PRE:
POST:
RETURN:   An array of unsigned chars
AUTHOR:
DATE:
MODIFIED: CMS-September 26, 1999
----------------------------------------------------------------------------- */
unsigned char *iiAverageGrayImages (int w, int h, unsigned char **grayimage,
                                    int numberofimages)
{
  unsigned char     *pixel = NULL;
  float             *imagesum = NULL, myfloat;
  unsigned long int  i, j, imagearea;

  imagearea = w*h;
  pixel = (unsigned char *)malloc(imagearea*sizeof(unsigned char));
  if (!pixel)
    return(pixel);
  memset(pixel, 0, imagearea*sizeof(unsigned char));
  imagesum = (float *)malloc(imagearea*sizeof(float));
  if (!imagesum)
    return(pixel);
  memset(imagesum, 0, imagearea*sizeof(float));

  for (i=0;i<numberofimages;i++) {
    if (grayimage[i]) {
      for (j=0;j<imagearea;j++) {
        imagesum[j] += (float)grayimage[i][j];
      }
    }
  }

  for (i=0;i<imagearea;i++) {
    myfloat = imagesum[i]/(float)numberofimages;
    pixel[i] = (unsigned char)(myfloat > (255.0-1.0e-30)) ? 255.0 : myfloat;
  }

  if (imagesum)
    free(imagesum);

  return pixel;
} /* iiAverageGrayImages */
/* -----------------------------------------------------------------------------
NAME:     iiiInterpolateGrayFromXY
PURPOSE:  Given x and y locations on an image and the image, this function
          interpolates the grayscale value for the given x and y on the image.
          If the x and y are off the image, 0 is returned.
PRE:      None.
POST:     None.
RETURN:   The grayscale value.
AUTHOR:   CMS
DATE:     September 28, 1999
MODIFIED:
----------------------------------------------------------------------------- */
unsigned char iiiInterpolateGrayFromXY(float xoriginal, float yoriginal,
                                       unsigned char *grayimage, int w, int h)
{
  long int ind1 = -1, ind2 = -1, ind3 = -1, ind4 = -1;
  int      i, j;
  float    x, y;

  if ((xoriginal < 0) || (xoriginal >= w) || (yoriginal < 0) ||
      (yoriginal >= h)) {
    return(0);
  }
  i = (int)xoriginal;
  j = (int)yoriginal;
  x = xoriginal - (float)i;
  y = yoriginal - (float)j;
  ind1 = j*w+i;
  if ((j+1) < h)
    ind2 = (j+1)*w+i;
  if ((i+1) < w)
    ind3 = j*w+i+1;
  if (((j+1) < h) && ((i+1) < w))
    ind4 = (j+1)*w+i+1;
  if ((ind2 != -1) && (ind3 != -1) && (ind4 != -1))
    return(((grayimage[ind3]-grayimage[ind1])*x)+
           ((grayimage[ind2]-grayimage[ind1])*y)+
           ((grayimage[ind4]+grayimage[ind1]-
             grayimage[ind2]-grayimage[ind3])*x*y)+grayimage[ind1]);
  else if (ind2 != -1)
    return(grayimage[ind2]);
  else if (ind3 != -1)
    return(grayimage[ind3]);
  else if (ind4 != -1)
    return(grayimage[ind4]);
  else
    return(0);
} /* iiiInterpolateGrayFromXY */
/* -----------------------------------------------------------------------------
NAME:     iiRotateGrayImage
PURPOSE:  Rotates a gray image.
PRE:      The image must be set up.
POST:     Must free the finalimage.
RETURN:   the rotated image, its size, and width.
AUTHOR:   CMS
DATE:     September 28, 1999
MODIFIED:
----------------------------------------------------------------------------- */
unsigned char *iiRotateGrayImage(int w, int h, unsigned char *grayimage,
                                 int *newwidth, int *newheight,
                                 float degreestorotate, float scale)
{
  float              xcenter, ycenter, curx, cury, xoriginal, yoriginal,
                     cosdeg, sindeg;
  int                u, v;
  unsigned long int  imagearea, ind;
  unsigned char     *pixel = NULL;

  *newwidth = (int)((float)w*scale);
  *newheight = (int)((float)h*scale);
  imagearea = (*newwidth)*(*newheight);
  pixel = (unsigned char *)malloc(imagearea*sizeof(unsigned char));
  if (!pixel)
    return(pixel);
  memset(pixel, 0, imagearea*sizeof(unsigned char));

  xcenter = (float)w/2.0;
  ycenter = (float)h/2.0;
  cosdeg = -cos((PI*degreestorotate)/180.0);
  sindeg = -sin((PI*degreestorotate)/180.0);

  for (u=0;u<(*newwidth);u++) {
    curx = (float)(*newwidth-1-u);
    for (v=0;v<(*newheight);v++) {
      ind = v*(*newwidth)+u;
      cury = (float)(*newheight-1-v);
      xoriginal = ((curx-xcenter)*cosdeg)-((cury-ycenter)*sindeg)+xcenter;
      yoriginal = ((curx-xcenter)*sindeg)+((cury-ycenter)*cosdeg)+ycenter;
      pixel[ind] = iiiInterpolateGrayFromXY(xoriginal, yoriginal, grayimage,
                                            w, h);
    }
  }
  return(pixel);
} /* iiRotateGrayImage */
/* -----------------------------------------------------------------------------
NAME:     iiRead1DFourierFile
PURPOSE:  Reads a 1-D fourier file
PRE:      None.
POST:     The fourier file is read and the array is allocated.  The array
          must be free'd when it is done being used.
RETURN:   The function array and the number of elements in that array.  MTRUE if
          successful, MFALSE otherwise.
AUTHOR:   CMS
DATE:     Oct 12, 1999
MODIFIED:
----------------------------------------------------------------------------- */
int iiRead1DFourierFile(char *filename, int *numfunc, float **func)
{
  FILE *fp;
  int   retval;
  float curfloat;

  fp = fopen(filename, "r");
  if (!fp)
    return(MFALSE);

  *func = NULL;
  *numfunc = 0;

  while ((retval = fscanf(fp, "%f", &curfloat)) != EOF) {
    if (retval == 1) {
      if (*func)
        (*func) = (float *)realloc((*func), (++(*numfunc))*sizeof(float));
      else
        (*func) = (float *)malloc((++(*numfunc))*sizeof(float));
      if (!(*func)) {
        *numfunc = 0;
        fclose(fp);
        return(MFALSE);
      }
      (*func)[(*numfunc)-1] = curfloat;
    }
  }
  fclose(fp);
  return(MTRUE);
} /* iiRead1DFourierFile */
/* -----------------------------------------------------------------------------
NAME:     iiFastFourier1D
PURPOSE:  Takes the fourier transform of a function and returns the complex
          and real parts.
PRE:      The function must exist.
POST:     The complex and real parts should be free'd after they are done
          being used.
RETURN:   MTRUE or MFALSE, also the fourier transform is returned.
AUTHOR:   CMS
DATE:     October 28, 1999
MODIFIED:
----------------------------------------------------------------------------- */
int iiFastFourier1D(float *realfunc, float *complexfunc, int numfunc,
                    float **realarray, float **complexarray, int inverse)
{
  float *oddrealarray = NULL, *oddcomplexarray = NULL,
        *oddcomplex = NULL, *oddreal = NULL;
  float *evenrealarray = NULL, *evencomplexarray = NULL,
        *evencomplex = NULL, *evenreal = NULL;
  float  arg, acoeff, bcoeff;
  int    numoddeven, i, index;

  (*realarray) = (float *)malloc(numfunc*sizeof(float));
  (*complexarray) = (float *)malloc(numfunc*sizeof(float));

  if (!(*realarray) || !(*complexarray)) {
    if (*realarray)
      free(*realarray);
    (*realarray) = NULL;
    if (*complexarray)
      free(*complexarray);
    (*complexarray) = NULL;
    return(MFALSE);
  }

  memset((*realarray), 0, numfunc*sizeof(float));
  memset((*complexarray), 0, numfunc*sizeof(float));

  if (numfunc == 2) {
    (*realarray)[0] = realfunc[0] + realfunc[1];
    (*realarray)[1] = realfunc[0] - realfunc[1];
    (*complexarray)[0] = complexfunc[0] + complexfunc[1];
    (*complexarray)[1] = complexfunc[0] - complexfunc[1];
    return(MTRUE);
  }
  else {
    numoddeven = numfunc/2;
    oddrealarray = (float *)malloc(numoddeven*sizeof(int));
    oddcomplexarray = (float *)malloc(numoddeven*sizeof(int));
    evenrealarray = (float *)malloc(numoddeven*sizeof(int));
    evencomplexarray = (float *)malloc(numoddeven*sizeof(int));
    if (!oddrealarray || !oddcomplexarray || !evenrealarray ||
        !evencomplexarray) {
      if (*realarray)
        free(*realarray);
      (*realarray) = NULL;
      if (*complexarray)
        free(*complexarray);
      (*complexarray) = NULL;
      if (oddrealarray)
        free(oddrealarray);
      oddrealarray = NULL;
      if (oddcomplexarray)
        free(oddcomplexarray);
      oddcomplexarray = NULL;
      if (evenrealarray)
        free(evenrealarray);
      evenrealarray = NULL;
      if (evencomplexarray)
        free(evencomplexarray);
      evencomplexarray = NULL;
      return(MFALSE);
    }
    memset(oddrealarray, 0, numoddeven*sizeof(float));
    memset(oddcomplexarray, 0, numoddeven*sizeof(float));
    memset(evenrealarray, 0, numoddeven*sizeof(float));
    memset(evencomplexarray, 0, numoddeven*sizeof(float));

    for (i=0;i<numfunc;i++) {
      if (i%2) {
        oddrealarray[i/2] = realfunc[i];
        oddcomplexarray[i/2] = complexfunc[i];
      }
      else {
        evenrealarray[i/2] = realfunc[i];
        evencomplexarray[i/2] = complexfunc[i];
      }
    }

    if (!iiFastFourier1D(evenrealarray, evencomplexarray, (numfunc/2),
                         &evenreal, &evencomplex, inverse)) {
      if (*realarray)
        free(*realarray);
      (*realarray) = NULL;
      if (*complexarray)
        free(*complexarray);
      (*complexarray) = NULL;
      if (oddrealarray)
        free(oddrealarray);
      oddrealarray = NULL;
      if (oddcomplexarray)
        free(oddcomplexarray);
      oddcomplexarray = NULL;
      if (evenrealarray)
        free(evenrealarray);
      evenrealarray = NULL;
      if (evencomplexarray)
        free(evencomplexarray);
      evencomplexarray = NULL;
      if (oddreal)
        free(oddreal);
      oddreal = NULL;
      if (oddcomplex)
        free(oddcomplex);
      oddcomplex = NULL;
      if (evenreal)
        free(evenreal);
      evenreal = NULL;
      if (evencomplex)
        free(evencomplex);
      evencomplex = NULL;
      return(MFALSE);
    }
    if (!iiFastFourier1D(oddrealarray, oddcomplexarray, (numfunc/2), &oddreal,
                         &oddcomplex, inverse)) {
      if (*realarray)
        free(*realarray);
      (*realarray) = NULL;
      if (*complexarray)
        free(*complexarray);
      (*complexarray) = NULL;
      if (oddrealarray)
        free(oddrealarray);
      oddrealarray = NULL;
      if (oddcomplexarray)
        free(oddcomplexarray);
      oddcomplexarray = NULL;
      if (evenrealarray)
        free(evenrealarray);
      evenrealarray = NULL;
      if (evencomplexarray)
        free(evencomplexarray);
      evencomplexarray = NULL;
      if (oddreal)
        free(oddreal);
      oddreal = NULL;
      if (oddcomplex)
        free(oddcomplex);
      oddcomplex = NULL;
      if (evenreal)
        free(evenreal);
      evenreal = NULL;
      if (evencomplex)
        free(evencomplex);
      evencomplex = NULL;
      return(MFALSE);
    }
    for (i=0;i<numfunc;i++) {
      index = i%(numfunc/2);
      arg = (2.0*PI*(float)i)/(float)numfunc;
      acoeff = cos(arg);
      bcoeff = sin(arg);
      if (inverse)
        (*realarray)[i] = evenreal[index]+(acoeff*oddreal[index])-
                          (bcoeff*oddcomplex[index]);
      else
        (*realarray)[i] = evenreal[index]+(acoeff*oddreal[index])+
                          (bcoeff*oddcomplex[index]);
      if (inverse)
        (*complexarray)[i] = evencomplex[index]+(acoeff*oddcomplex[index])+
                             (bcoeff*oddreal[index]);
      else
        (*complexarray)[i] = evencomplex[index]+(acoeff*oddcomplex[index])-
                             (bcoeff*oddreal[index]);
    }
  }
  if (oddrealarray)
    free(oddrealarray);
  oddrealarray = NULL;
  if (oddcomplexarray)
    free(oddcomplexarray);
  oddcomplexarray = NULL;
  if (evenrealarray)
    free(evenrealarray);
  evenrealarray = NULL;
  if (evencomplexarray)
    free(evencomplexarray);
  evencomplexarray = NULL;
  if (oddreal)
    free(oddreal);
  oddreal = NULL;
  if (oddcomplex)
    free(oddcomplex);
  oddcomplex = NULL;
  if (evenreal)
    free(evenreal);
  evenreal = NULL;
  if (evencomplex)
    free(evencomplex);
  evencomplex = NULL;
  return(MTRUE);
} /* iiFastFourier1D */
/* -----------------------------------------------------------------------------
NAME:     iiRemoveHighFrequencies
PURPOSE:  Removes the high-magnitude frequencies from a frequency-domain signal.
PRE:
POST:
RETURN:   none.
AUTHOR:   CMS
DATE:     November 12, 1999
MODIFIED:
----------------------------------------------------------------------------- */
int iiRemoveHighFrequencies(float *realimage, float *compleximage,
                            int width, int height)
{
  unsigned long int  imagearea, i;
  float             *magnitude = NULL, minmag, maxmag, myfloat, compare;
  unsigned char     *pixel = NULL;

  imagearea = width*height;

  magnitude = (float *)malloc(imagearea*sizeof(float));
  if (!magnitude) {
    if (magnitude)
      free(magnitude);
    return(MFALSE);
  }
  memset(magnitude, 0, imagearea*sizeof(float));
  minmag = 1.0e30;
  maxmag = -1.0e30;

  for (i=0;i<imagearea;i++) {
    magnitude[i] = sqrt((realimage[i]*realimage[i])+
                        (compleximage[i]*compleximage[i]));
    if (magnitude[i] > maxmag)
      maxmag = magnitude[i];
    if (magnitude[i] < minmag)
      minmag = magnitude[i];
  }

  compare = (maxmag - minmag)*0.15;
  for (i=0;i<imagearea;i++) {
    /* If this pixel's magnitude is in the upper 90% of the magnitudes,
       set its magnitude to 0 */
    if ((magnitude[i]-minmag) > compare) {
      realimage[i] = compleximage[i] = 0.0;
    }
  }

  if (magnitude)
    free(magnitude);
  return(MTRUE);
} /* iiRemoveHighFrequencies */
/* -----------------------------------------------------------------------------
NAME:     iiWriteFourierResult
PURPOSE:  Writes out the real and complex parts, magnitude and phase of the
          fourier transform result to a PGM image file.
PRE:      The fourier transform must exist.
POST:
RETURN:   MTRUE or MFALSE.
AUTHOR:   CMS
DATE:     November 30, 1999
MODIFIED:
----------------------------------------------------------------------------- */
int iiWriteFourierResult(char *outrealimagfilename, char *outmagphasefilename,
                         float *realarray, float *complexarray,
                         int streamlength)
{
  float magnitude, phase;
  int   i;
  FILE *fp1, *fp2;

  if (!(fp1 = fopen(outrealimagfilename, "w"))) {
    printf("Cannot open file %s for writing", outrealimagfilename);
    return(MFALSE);
  }
  if (!(fp2 = fopen(outmagphasefilename, "w"))) {
    printf("Cannot open file %s for writing", outmagphasefilename);
    return(MFALSE);
  }

  for (i=0;i<streamlength;i++) {
    magnitude = (realarray[i]*realarray[i])+(complexarray[i]*complexarray[i]);
    fprintf(fp1, "%1.9lf\n", realarray[i]);
    fprintf(fp2, "%1.9lf\n", magnitude);
  }

  if (fp1)
    fclose(fp1);
  if (fp2)
    fclose(fp2);

  return(MTRUE);
} /* iiWriteFourierResult */
/* -----------------------------------------------------------------------------
NAME:     iiFastFourier2D
PURPOSE:  Computes the 2D Fourier transform of an image.
PRE:
POST:
RETURN:   Returns the complex and imaginary parts of the Fourier transform
          of the image.
AUTHOR:   CMS
DATE:     November 3, 1999
MODIFIED:
----------------------------------------------------------------------------- */
int iiFastFourier2D(float *inputrealimage, float *inputcompleximage, int width,
                    int height, float **realimage, float **compleximage,
                    int inverse)
{
  unsigned long int  i, j, imagearea, index;
  float             *realarray = NULL, *complexarray = NULL,
                    *realfloatarray = NULL, *complexfloatarray;

  imagearea = width*height;
  (*realimage) = (float *)malloc(imagearea*sizeof(float));
  (*compleximage) = (float *)malloc(imagearea*sizeof(float));
  realfloatarray = (float *)malloc(((width>height)?width:height)*sizeof(float));
  complexfloatarray = (float *)malloc(((width>height)?width:height)*
                                        sizeof(float));
  if (!(*realimage) || !(*compleximage) || !realfloatarray ||
      !complexfloatarray) {
    if (*realimage)
      free(*realimage);
    (*realimage) = NULL;
    if (*compleximage)
      free(*compleximage);
    (*compleximage) = NULL;
    if (realfloatarray)
      free(realfloatarray);
    (realfloatarray) = NULL;
    if (complexfloatarray)
      free(complexfloatarray);
    (complexfloatarray) = NULL;
    return(MFALSE);
  }
  memset(*realimage, 0, imagearea*sizeof(float));
  memset(*compleximage, 0, imagearea*sizeof(float));
  memset(realfloatarray, 0, ((width>height)?width:height)*sizeof(float));
  memset(complexfloatarray, 0, ((width>height)?width:height)*sizeof(float));
  for (i=0;i<height;i++) {
    for (j=0;j<width;j++) {
      index = i*width+j;
      realfloatarray[j] = inputrealimage[index];
      complexfloatarray[j] = inputcompleximage[index];
    }
    if (!iiFastFourier1D(realfloatarray, complexfloatarray, width, &realarray,
                         &complexarray, inverse)) {
      if (*realimage)
        free(*realimage);
      (*realimage) = NULL;
      if (*compleximage)
        free(*compleximage);
      (*compleximage) = NULL;
      if (realfloatarray)
        free(realfloatarray);
      (realfloatarray) = NULL;
      if (complexfloatarray)
        free(complexfloatarray);
      (complexfloatarray) = NULL;
      return(MFALSE);
    }
    memcpy((*realimage)+(i*width), realarray, width*sizeof(float));
    memcpy((*compleximage)+(i*width), complexarray, width*sizeof(float));
    free(realarray);
    realarray = NULL;
    free(complexarray);
    complexarray = NULL;
  }
  for (i=0;i<width;i++) {
    for (j=0;j<height;j++) {
      index = j*width+i;
      realfloatarray[j] = (*realimage)[index];
      complexfloatarray[j] = (*compleximage)[index];
    }
    if (!iiFastFourier1D(realfloatarray, complexfloatarray, height, &realarray,
                         &complexarray, inverse)) {
      if (*realimage)
        free(*realimage);
      (*realimage) = NULL;
      if (*compleximage)
        free(*compleximage);
      (*compleximage) = NULL;
      if (realfloatarray)
        free(realfloatarray);
      (realfloatarray) = NULL;
      if (complexfloatarray)
        free(complexfloatarray);
      (complexfloatarray) = NULL;
      return(MFALSE);
    }
    for (j=0;j<height;j++) {
      index = j*width+i;
      (*realimage)[index] = realarray[j];
      (*compleximage)[index] = complexarray[j];
    }
    free(realarray);
    realarray = NULL;
    free(complexarray);
    complexarray = NULL;
  }
  free(realfloatarray);
  realfloatarray = NULL;
  free(complexfloatarray);
  complexfloatarray = NULL;
  return(MTRUE);
} /* iiFastFourier2D */
/* -----------------------------------------------------------------------------
NAME:     iiReadOffsetLinear8BitAudioFile
PURPOSE:  Reads an offset linear 8-bit audio file into an unsigned char array
PRE:      The fp should be opened
POST:     The fp should be closed
RETURN:   An array of unsigned chars
AUTHOR:   CMS
DATE:
MODIFIED: CMS-November 30, 1999
----------------------------------------------------------------------------- */
unsigned char *iiReadOffsetLinear8BitAudioFile (FILE *fp, int *streamlength)
{
  int            i, done, blocksize = 10000, oldblocksize, numread;
  unsigned char *audiostream = NULL;

  if (!fp)
    return(audiostream);
  audiostream = (unsigned char *)malloc(blocksize * sizeof(unsigned char));
  if (!audiostream)
    return(audiostream);
  for (i=0; i<blocksize; i++) {
    numread = fscanf(fp, "%c", audiostream+i);
    if (numread != 1) {
      blocksize = i;
      audiostream = (unsigned char *)realloc(audiostream, blocksize *
                                             sizeof(unsigned char));
      if (!audiostream)
        return(audiostream);
    }
    else if ((numread == 1) && (i == blocksize-1)) {
      blocksize+=blocksize;
      audiostream = (unsigned char *)realloc(audiostream, blocksize *
                                             sizeof(unsigned char));
      if (!audiostream)
        return(audiostream);
    }
  }
  /* Pad with 0's if size is not a power of 2 */
  done = MFALSE;
  oldblocksize = blocksize;
  for (i=0;(i<50) && !done;i++) {
    if (blocksize == (int)(pow(2.0, (double)i)+.001))
      done = MTRUE;
    else if ((double)blocksize > pow(2.0, (double)i) &&
             (double)blocksize < pow(2.0, (double)(i+1))) {
      blocksize = (int)(pow(2.0, (double)(i+1))+.001);
      audiostream = (unsigned char *)realloc(audiostream, blocksize *
                                             sizeof(unsigned char));
      memset(audiostream+oldblocksize, 0, (blocksize-oldblocksize)*
             sizeof(unsigned char));
      done = MTRUE;
    }
  }
  *streamlength = blocksize;
  return audiostream;
} /* iiReadOffsetLinear8BitAudioFile */
/* -----------------------------------------------------------------------------
NAME:     main
PURPOSE:  This is the main function, buckeroo!
PRE:
POST:
RETURN:
AUTHOR:
DATE:
MODIFIED:
----------------------------------------------------------------------------- */
int main(int argc, char *argv[])
{
  unsigned char    *inputaudiostream = NULL;
  float            *realfloatarray = NULL, *complexfloatarray = NULL,
                   *realarray = NULL, *complexarray = NULL;
  int               audiostreamlength = 0, i;
  FILE             *fp = NULL;
  unsigned long int imagesize;

  if (argc != 4) {
    printf("Usage: fouriersound <infilename> <outrealimagfilename> "
           "<outmagphasefilename>\n");
    return(1);
  }

  if (!(fp = fopen(argv[1], "r"))) {
    printf("Cannot open file %s", argv[1]);
    return(2);
  }

  if (!(inputaudiostream = iiReadOffsetLinear8BitAudioFile(fp,
        &audiostreamlength)))
    return(3);

  if (fp)
    fclose(fp);

  realfloatarray = (float *)malloc(audiostreamlength*sizeof(float));
  complexfloatarray = (float *)malloc(audiostreamlength*sizeof(float));
  memset(realfloatarray, 0, audiostreamlength*sizeof(float));
  memset(complexfloatarray, 0, audiostreamlength*sizeof(float));

  for (i=0;i<audiostreamlength;i++) {
    realfloatarray[i] = (float)inputaudiostream[i];
  }

  if (!iiFastFourier1D(realfloatarray, complexfloatarray, audiostreamlength,
                       &realarray, &complexarray, MFALSE)) {
    if (realfloatarray)
      free(realfloatarray);
    if (complexfloatarray)
      free(complexfloatarray);
    if (realarray);
      free(realarray);
    if (complexarray)
      free(complexarray);
  }

  if (!iiWriteFourierResult(argv[2], argv[3], realarray, complexarray,
                            audiostreamlength)) {
    if (realfloatarray)
      free(realfloatarray);
    if (complexfloatarray)
      free(complexfloatarray);
    if (realarray);
      free(realarray);
    if (complexarray)
      free(complexarray);
  }

  if (realfloatarray)
    free(realfloatarray);
  if (complexfloatarray)
    free(complexfloatarray);
  if (realarray);
    free(realarray);
  if (complexarray)
    free(complexarray);

  return(0);
} /* main */