text: /*
   * HRTF utility for producing and demonstrating the process of creating an
   * OpenAL Soft compatible HRIR data set.
   *
   * Copyright (C) 2011-2012  Christopher Fitzgerald
   *
   * This program is free software; you can redistribute it and/or modify
   * it under the terms of the GNU General Public License as published by
   * the Free Software Foundation; either version 2 of the License, or
   * (at your option) any later version.
   *
   * This program is distributed in the hope that it will be useful,
   * but WITHOUT ANY WARRANTY; without even the implied warranty of
   * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   * GNU General Public License for more details.
   *
   * You should have received a copy of the GNU General Public License along
   * with this program; if not, write to the Free Software Foundation, Inc.,
   * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
   *
   * Or visit:  http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
   *
   * --------------------------------------------------------------------------
   *
   * A big thanks goes out to all those whose work done in the field of
   * binaural sound synthesis using measured HRTFs makes this utility and the
   * OpenAL Soft implementation possible.
   *
   * The algorithm for diffuse-field equalization was adapted from the work
   * done by Rio Emmanuel and Larcher Veronique of IRCAM and Bill Gardner of
   * MIT Media Laboratory.  It operates as follows:
   *
   *  1.  Take the FFT of each HRIR and only keep the magnitude responses.
   *  2.  Calculate the diffuse-field power-average of all HRIRs weighted by
   *      their contribution to the total surface area covered by their
   *      measurement.
   *  3.  Take the diffuse-field average and limit its magnitude range.
   *  4.  Equalize the responses by using the inverse of the diffuse-field
   *      average.
   *  5.  Reconstruct the minimum-phase responses.
   *  5.  Zero the DC component.
   *  6.  IFFT the result and truncate to the desired-length minimum-phase FIR.
   *
   * The spherical head algorithm for calculating propagation delay was adapted
   * from the paper:
   *
   *  Modeling Interaural Time Difference Assuming a Spherical Head
   *  Joel David Miller
   *  Music 150, Musical Acoustics, Stanford University
   *  December 2, 2001
   *
   * The formulae for calculating the Kaiser window metrics are from the
   * the textbook:
   *
   *  Discrete-Time Signal Processing
   *  Alan V. Oppenheim and Ronald W. Schafer
   *  Prentice-Hall Signal Processing Series
   *  1999
   */
  
  /* Needed for 64-bit unsigned integer. */
  #include "config.h"
  
  #include <stdio.h>
  #include <stdlib.h>
  #include <stdarg.h>
  #include <string.h>
  #include <ctype.h>
  #include <math.h>
  
  // Rely (if naively) on OpenAL's header for the types used for serialization.
  #include "AL/al.h"
  
  #ifndef M_PI
  #define M_PI                         (3.14159265358979323846)
  #endif
  
  #ifndef HUGE_VAL
  #define HUGE_VAL                     (1.0 / 0.0)
  #endif
  
  // The epsilon used to maintain signal stability.
  #define EPSILON                      (1e-15)
  
  // Constants for accessing the token reader's ring buffer.
  #define TR_RING_BITS                 (16)
  #define TR_RING_SIZE                 (1 << TR_RING_BITS)
  #define TR_RING_MASK                 (TR_RING_SIZE - 1)
  
  // The token reader's load interval in bytes.
  #define TR_LOAD_SIZE                 (TR_RING_SIZE >> 2)
  
  // The maximum identifier length used when processing the data set
  // definition.
  #define MAX_IDENT_LEN                (16)
  
  // The maximum path length used when processing filenames.
  #define MAX_PATH_LEN                 (256)
  
  // The limits for the sample 'rate' metric in the data set definition and for
  // resampling.
  #define MIN_RATE                     (32000)
  #define MAX_RATE                     (96000)
  
  // The limits for the HRIR 'points' metric in the data set definition.
  #define MIN_POINTS                   (16)
  #define MAX_POINTS                   (8192)
  
  // The limits to the number of 'azimuths' listed in the data set definition.
  #define MIN_EV_COUNT                 (5)
  #define MAX_EV_COUNT                 (128)
  
  // The limits for each of the 'azimuths' listed in the data set definition.
  #define MIN_AZ_COUNT                 (1)
  #define MAX_AZ_COUNT                 (128)
  
  // The limits for the listener's head 'radius' in the data set definition.
  #define MIN_RADIUS                   (0.05)
  #define MAX_RADIUS                   (0.15)
  
  // The limits for the 'distance' from source to listener in the definition
  // file.
  #define MIN_DISTANCE                 (0.5)
  #define MAX_DISTANCE                 (2.5)
  
  // The maximum number of channels that can be addressed for a WAVE file
  // source listed in the data set definition.
  #define MAX_WAVE_CHANNELS            (65535)
  
  // The limits to the byte size for a binary source listed in the definition
  // file.
  #define MIN_BIN_SIZE                 (2)
  #define MAX_BIN_SIZE                 (4)
  
  // The minimum number of significant bits for binary sources listed in the
  // data set definition.  The maximum is calculated from the byte size.
  #define MIN_BIN_BITS                 (16)
  
  // The limits to the number of significant bits for an ASCII source listed in
  // the data set definition.
  #define MIN_ASCII_BITS               (16)
  #define MAX_ASCII_BITS               (32)
  
  // The limits to the FFT window size override on the command line.
  #define MIN_FFTSIZE                  (512)
  #define MAX_FFTSIZE                  (16384)
  
  // The limits to the equalization range limit on the command line.
  #define MIN_LIMIT                    (2.0)
  #define MAX_LIMIT                    (120.0)
  
  // The limits to the truncation window size on the command line.
  #define MIN_TRUNCSIZE                (8)
  #define MAX_TRUNCSIZE                (128)
  
  // The truncation window size must be a multiple of the below value to allow
  // for vectorized convolution.
  #define MOD_TRUNCSIZE                (8)
  
  // The defaults for the command line options.
  #define DEFAULT_EQUALIZE             (1)
  #define DEFAULT_SURFACE              (1)
  #define DEFAULT_LIMIT                (24.0)
  #define DEFAULT_TRUNCSIZE            (32)
  
  // The four-character-codes for RIFF/RIFX WAVE file chunks.
  #define FOURCC_RIFF                  (0x46464952) // 'RIFF'
  #define FOURCC_RIFX                  (0x58464952) // 'RIFX'
  #define FOURCC_WAVE                  (0x45564157) // 'WAVE'
  #define FOURCC_FMT                   (0x20746D66) // 'fmt '
  #define FOURCC_DATA                  (0x61746164) // 'data'
  #define FOURCC_LIST                  (0x5453494C) // 'LIST'
  #define FOURCC_WAVL                  (0x6C766177) // 'wavl'
  #define FOURCC_SLNT                  (0x746E6C73) // 'slnt'
  
  // The supported wave formats.
  #define WAVE_FORMAT_PCM              (0x0001)
  #define WAVE_FORMAT_IEEE_FLOAT       (0x0003)
  #define WAVE_FORMAT_EXTENSIBLE       (0xFFFE)
  
  // The maximum propagation delay value supported by OpenAL Soft.
  #define MAX_HRTD                     (63.0)
  
  // The OpenAL Soft HRTF format marker.  It stands for minimum-phase head
  // response protocol 01.
  #define MHR_FORMAT                   ("MinPHR01")
  
  // Byte order for the serialization routines.
  enum ByteOrderT {
    BO_NONE   = 0,
    BO_LITTLE    ,
    BO_BIG
  };
  
  // Source format for the references listed in the data set definition.
  enum SourceFormatT {
    SF_NONE   = 0,
    SF_WAVE      ,  // RIFF/RIFX WAVE file.
    SF_BIN_LE    ,  // Little-endian binary file.
    SF_BIN_BE    ,  // Big-endian binary file.
    SF_ASCII        // ASCII text file.
  };
  
  // Element types for the references listed in the data set definition.
  enum ElementTypeT {
    ET_NONE = 0,
    ET_INT     ,    // Integer elements.
    ET_FP           // Floating-point elements.
  };
  
  // Desired output format from the command line.
  enum OutputFormatT {
    OF_NONE  = 0,
    OF_MHR      ,   // OpenAL Soft MHR data set file.
    OF_TABLE        // OpenAL Soft built-in table file (used when compiling).
  };
  
  // Unsigned integer type.
  typedef unsigned int                 uint;
  
  // Serialization types.  The trailing digit indicates the number of bytes.
  typedef ALubyte                      uint1;
  
  typedef ALint                        int4;
  typedef ALuint                       uint4;
  
  #if defined (HAVE_STDINT_H)
  #include <stdint.h>
  
  typedef uint64_t                     uint8;
  #elif defined (HAVE___INT64)
  typedef unsigned __int64             uint8;
  #elif (SIZEOF_LONG == 8)
  typedef unsigned long                uint8;
  #elif (SIZEOF_LONG_LONG == 8)
  typedef unsigned long long           uint8;
  #endif
  
  typedef enum ByteOrderT              ByteOrderT;
  typedef enum SourceFormatT           SourceFormatT;
  typedef enum ElementTypeT            ElementTypeT;
  typedef enum OutputFormatT           OutputFormatT;
  
  typedef struct TokenReaderT          TokenReaderT;
  typedef struct SourceRefT            SourceRefT;
  typedef struct HrirDataT             HrirDataT;
  typedef struct ResamplerT            ResamplerT;
  
  // Token reader state for parsing the data set definition.
  struct TokenReaderT {
    FILE                             * mFile;
    const char                       * mName;
    uint                               mLine,
                                       mColumn;
    char                               mRing [TR_RING_SIZE];
    size_t                             mIn,
                                       mOut;
  };
  
  // Source reference state used when loading sources.
  struct SourceRefT {
    SourceFormatT                      mFormat;
    ElementTypeT                       mType;
    uint                               mSize;
    int                                mBits;
    uint                               mChannel,
                                       mSkip,
                                       mOffset;
    char                               mPath [MAX_PATH_LEN + 1];
  };
  
  // The HRIR metrics and data set used when loading, processing, and storing
  // the resulting HRTF.
  struct HrirDataT {
    uint                               mIrRate,
                                       mIrCount,
                                       mIrSize,
                                       mIrPoints,
                                       mFftSize,
                                       mEvCount,
                                       mEvStart,
                                       mAzCount [MAX_EV_COUNT],
                                       mEvOffset [MAX_EV_COUNT];
    double                             mRadius,
                                       mDistance,
                                     * mHrirs,
                                     * mHrtds,
                                       mMaxHrtd;
  };
  
  // The resampler metrics and FIR filter.
  struct ResamplerT {
    uint                               mP,
                                       mQ,
                                       mM,
                                       mL;
    double                           * mF;
  };
  
  /* Token reader routines for parsing text files.  Whitespace is not
   * significant.  It can process tokens as identifiers, numbers (integer and
   * floating-point), strings, and operators.  Strings must be encapsulated by
   * double-quotes and cannot span multiple lines.
   */
  
  // Setup the reader on the given file.  The filename can be NULL if no error
  // output is desired.
  static void TrSetup (FILE * fp, const char * filename, TokenReaderT * tr) {
    const char * name = NULL;
    char ch;
  
    tr -> mFile = fp;
    name = filename;
    // If a filename was given, store a pointer to the base name.
    if (filename != NULL) {
       while ((ch = (* filename)) != '\0') {
         if ((ch == '/') || (ch == '\\'))
            name = filename + 1;
         filename ++;
       }
    }
    tr -> mName = name;
    tr -> mLine = 1;
    tr -> mColumn = 1;
    tr -> mIn = 0;
    tr -> mOut = 0;
  }
  
  // Prime the reader's ring buffer, and return a result indicating that there
  // is text to process.
  static int TrLoad (TokenReaderT * tr) {
    size_t toLoad, in, count;
  
    toLoad = TR_RING_SIZE - (tr -> mIn - tr -> mOut);
    if ((toLoad >= TR_LOAD_SIZE) && (! feof (tr -> mFile))) {
       // Load TR_LOAD_SIZE (or less if at the end of the file) per read.
       toLoad = TR_LOAD_SIZE;
       in = tr -> mIn & TR_RING_MASK;
       count = TR_RING_SIZE - in;
       if (count < toLoad) {
          tr -> mIn += fread (& tr -> mRing [in], 1, count, tr -> mFile);
          tr -> mIn += fread (& tr -> mRing [0], 1, toLoad - count, tr -> mFile);
       } else {
          tr -> mIn += fread (& tr -> mRing [in], 1, toLoad, tr -> mFile);
       }
       if (tr -> mOut >= TR_RING_SIZE) {
          tr -> mOut -= TR_RING_SIZE;
          tr -> mIn -= TR_RING_SIZE;
       }
    }
    if (tr -> mIn > tr -> mOut)
       return (1);
    return (0);
  }
  
  // Error display routine.  Only displays when the base name is not NULL.
  static void TrErrorVA (const TokenReaderT * tr, uint line, uint column, const char * format, va_list argPtr) {
    if (tr -> mName != NULL) {
       fprintf (stderr, "Error (%s:%u:%u):  ", tr -> mName, line, column);
       vfprintf (stderr, format, argPtr);
    }
  }
  
  // Used to display an error at a saved line/column.
  static void TrErrorAt (const TokenReaderT * tr, uint line, uint column, const char * format, ...) {
    va_list argPtr;
  
    va_start (argPtr, format);
    TrErrorVA (tr, line, column, format, argPtr);
    va_end (argPtr);
  }
  
  // Used to display an error at the current line/column.
  static void TrError (const TokenReaderT * tr, const char * format, ...) {
    va_list argPtr;
  
    va_start (argPtr, format);
    TrErrorVA (tr, tr -> mLine, tr -> mColumn, format, argPtr);
    va_end (argPtr);
  }
  
  // Skips to the next line.
  static void TrSkipLine (TokenReaderT * tr) {
    char ch;
  
    while (TrLoad (tr)) {
      ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
      tr -> mOut ++;
      if (ch == '\n') {
         tr -> mLine ++;
         tr -> mColumn = 1;
         break;
      }
      tr -> mColumn ++;
    }
  }
  
  // Skips to the next token.
  static int TrSkipWhitespace (TokenReaderT * tr) {
    char ch;
  
    while (TrLoad (tr)) {
      ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
      if (isspace (ch)) {
         tr -> mOut ++;
         if (ch == '\n') {
            tr -> mLine ++;
            tr -> mColumn = 1;
         } else {
            tr -> mColumn ++;
         }
      } else if (ch == '#') {
         TrSkipLine (tr);
      } else {
         return (1);
      }
    }
    return (0);
  }
  
  // Get the line and/or column of the next token (or the end of input).
  static void TrIndication (TokenReaderT * tr, uint * line, uint * column) {
    TrSkipWhitespace (tr);
    if (line != NULL)
       (* line) = tr -> mLine;
    if (column != NULL)
       (* column) = tr -> mColumn;
  }
  
  // Checks to see if a token is the given operator.  It does not display any
  // errors and will not proceed to the next token.
  static int TrIsOperator (TokenReaderT * tr, const char * op) {
    size_t out, len;
    char ch;
  
    if (! TrSkipWhitespace (tr))
       return (0);
    out = tr -> mOut;
    len = 0;
    while ((op [len] != '\0') && (out < tr -> mIn)) {
      ch = tr -> mRing [out & TR_RING_MASK];
      if (ch != op [len])
         break;
      len ++;
      out ++;
    }
    if (op [len] == '\0')
       return (1);
    return (0);
  }
  
  /* The TrRead*() routines obtain the value of a matching token type.  They
   * display type, form, and boundary errors and will proceed to the next
   * token.
   */
  
  // Reads and validates an identifier token.
  static int TrReadIdent (TokenReaderT * tr, const uint maxLen, char * ident) {
    uint col, len;
    char ch;
  
    col = tr -> mColumn;
    if (TrSkipWhitespace (tr)) {
       col = tr -> mColumn;
       ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
       if ((ch == '_') || isalpha (ch)) {
          len = 0;
          do {
            if (len < maxLen)
               ident [len] = ch;
            len ++;
            tr -> mOut ++;
            if (! TrLoad (tr))
               break;
            ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
          } while ((ch == '_') || isdigit (ch) || isalpha (ch));
          tr -> mColumn += len;
          if (len > maxLen) {
             TrErrorAt (tr, tr -> mLine, col, "Identifier is too long.\n");
             return (0);
          }
          ident [len] = '\0';
          return (1);
       }
    }
    TrErrorAt (tr, tr -> mLine, col, "Expected an identifier.\n");
    return (0);
  }
  
  // Reads and validates (including bounds) an integer token.
  static int TrReadInt (TokenReaderT * tr, const int loBound, const int hiBound, int * value) {
    uint col, digis, len;
    char ch, temp [64 + 1];
  
    col = tr -> mColumn;
    if (TrSkipWhitespace (tr)) {
       col = tr -> mColumn;
       len = 0;
       ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
       if ((ch == '+') || (ch == '-')) {
          temp [len] = ch;
          len ++;
          tr -> mOut ++;
       }
       digis = 0;
       while (TrLoad (tr)) {
         ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
         if (! isdigit (ch))
            break;
         if (len < 64)
            temp [len] = ch;
         len ++;
         digis ++;
         tr -> mOut ++;
       }
       tr -> mColumn += len;
       if ((digis > 0) && (ch != '.') && (! isalpha (ch))) {
          if (len > 64) {
             TrErrorAt (tr, tr -> mLine, col, "Integer is too long.");
             return (0);
          }
          temp [len] = '\0';
          (* value) = strtol (temp, NULL, 10);
          if (((* value) < loBound) || ((* value) > hiBound)) {
             TrErrorAt (tr, tr -> mLine, col, "Expected a value from %d to %d.\n", loBound, hiBound);
             return (0);
          }
          return (1);
       }
    }
    TrErrorAt (tr, tr -> mLine, col, "Expected an integer.\n");
    return (0);
  }
  
  // Reads and validates (including bounds) a float token.
  static int TrReadFloat (TokenReaderT * tr, const double loBound, const double hiBound, double * value) {
    uint col, digis, len;
    char ch, temp [64 + 1];
  
    col = tr -> mColumn;
    if (TrSkipWhitespace (tr)) {
       col = tr -> mColumn;
       len = 0;
       ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
       if ((ch == '+') || (ch == '-')) {
          temp [len] = ch;
          len ++;
          tr -> mOut ++;
       }
       digis = 0;
       while (TrLoad (tr)) {
         ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
         if (! isdigit (ch))
            break;
         if (len < 64)
            temp [len] = ch;
         len ++;
         digis ++;
         tr -> mOut ++;
       }
       if (ch == '.') {
          if (len < 64)
             temp [len] = ch;
          len ++;
          tr -> mOut ++;
       }
       while (TrLoad (tr)) {
         ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
         if (! isdigit (ch))
            break;
         if (len < 64)
            temp [len] = ch;
         len ++;
         digis ++;
         tr -> mOut ++;
       }
       if (digis > 0) {
          if ((ch == 'E') || (ch == 'e')) {
             if (len < 64)
                temp [len] = ch;
             len ++;
             digis = 0;
             tr -> mOut ++;
             if ((ch == '+') || (ch == '-')) {
                if (len < 64)
                   temp [len] = ch;
                len ++;
                tr -> mOut ++;
             }
             while (TrLoad (tr)) {
               ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
               if (! isdigit (ch))
                  break;
               if (len < 64)
                  temp [len] = ch;
               len ++;
               digis ++;
               tr -> mOut ++;
             }
          }
          tr -> mColumn += len;
          if ((digis > 0) && (ch != '.') && (! isalpha (ch))) {
             if (len > 64) {
                TrErrorAt (tr, tr -> mLine, col, "Float is too long.");
                return (0);
             }
             temp [len] = '\0';
             (* value) = strtod (temp, NULL);
             if (((* value) < loBound) || ((* value) > hiBound)) {
                TrErrorAt (tr, tr -> mLine, col, "Expected a value from %f to %f.\n", loBound, hiBound);
                return (0);
             }
             return (1);
          }
       } else {
          tr -> mColumn += len;
       }
    }
    TrErrorAt (tr, tr -> mLine, col, "Expected a float.\n");
    return (0);
  }
  
  // Reads and validates a string token.
  static int TrReadString (TokenReaderT * tr, const uint maxLen, char * text) {
    uint col, len;
    char ch;
  
    col = tr -> mColumn;
    if (TrSkipWhitespace (tr)) {
       col = tr -> mColumn;
       ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
       if (ch == '\"') {
          tr -> mOut ++;
          len = 0;
          while (TrLoad (tr)) {
            ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
            tr -> mOut ++;
            if (ch == '\"')
               break;
            if (ch == '\n') {
               TrErrorAt (tr, tr -> mLine, col, "Unterminated string at end of line.\n");
               return (0);
            }
            if (len < maxLen)
               text [len] = ch;
            len ++;
          }
          if (ch != '\"') {
             tr -> mColumn += 1 + len;
             TrErrorAt (tr, tr -> mLine, col, "Unterminated string at end of input.\n");
             return (0);
          }
          tr -> mColumn += 2 + len;
          if (len > maxLen) {
             TrErrorAt (tr, tr -> mLine, col, "String is too long.\n");
             return (0);
          }
          text [len] = '\0';
          return (1);
       }
    }
    TrErrorAt (tr, tr -> mLine, col, "Expected a string.\n");
    return (0);
  }
  
  // Reads and validates the given operator.
  static int TrReadOperator (TokenReaderT * tr, const char * op) {
    uint col, len;
    char ch;
  
    col = tr -> mColumn;
    if (TrSkipWhitespace (tr)) {
       col = tr -> mColumn;
       len = 0;
       while ((op [len] != '\0') && TrLoad (tr)) {
         ch = tr -> mRing [tr -> mOut & TR_RING_MASK];
         if (ch != op [len])
            break;
         len ++;
         tr -> mOut ++;
       }
       tr -> mColumn += len;
       if (op [len] == '\0')
          return (1);
    }
    TrErrorAt (tr, tr -> mLine, col, "Expected '%s' operator.\n", op);
    return (0);
  }
  
  /* Performs a string substitution.  Any case-insensitive occurrences of the
   * pattern string are replaced with the replacement string.  The result is
   * truncated if necessary.
   */
  static int StrSubst (const char * in, const char * pat, const char * rep, const size_t maxLen, char * out) {
    size_t inLen, patLen, repLen;
    size_t si, di;
    int truncated;
  
    inLen = strlen (in);
    patLen = strlen (pat);
    repLen = strlen (rep);
    si = 0;
    di = 0;
    truncated = 0;
    while ((si < inLen) && (di < maxLen)) {
      if (patLen <= (inLen - si)) {
         if (strncasecmp (& in [si], pat, patLen) == 0) {
            if (repLen > (maxLen - di)) {
               repLen = maxLen - di;
               truncated = 1;
            }
            strncpy (& out [di], rep, repLen);
            si += patLen;
            di += repLen;
         }
      }
      out [di] = in [si];
      si ++;
      di ++;
    }
    if (si < inLen)
       truncated = 1;
    out [di] = '\0';
    return (! truncated);
  }
  
  // Provide missing math routines for MSVC.
  #undef _MSC_VER
  #ifdef _MSC_VER
  static double round (double val) {
    if (val < 0.0)
       return (ceil (val - 0.5));
    return (floor (val + 0.5));
  }
  
  static double fmin (double a, double b) {
    return ((a < b) ? a : b);
  }
  
  static double fmax (double a, double b) {
    return ((a > b) ? a : b);
  }
  #endif
  
  // Simple clamp routine.
  static double Clamp (const double val, const double lower, const double upper) {
    return (fmin (fmax (val, lower), upper));
  }
  
  // Performs linear interpolation.
  static double Lerp (const double a, const double b, const double f) {
    return (a + (f * (b - a)));
  }
  
  // Performs a high-passed triangular probability density function dither from
  // a double to an integer.  It assumes the input sample is already scaled.
  static int HpTpdfDither (const double in, int * hpHist) {
    const double PRNG_SCALE = 1.0 / (RAND_MAX + 1.0);
    int prn;
    double out;
  
    prn = rand ();
    out = round (in + (PRNG_SCALE * (prn - (* hpHist))));
    (* hpHist) = prn;
    return ((int) out);
  }
  
  // Allocates an array of doubles.
  static double * CreateArray (const size_t n) {
    double * a = NULL;
  
    a = (double *) calloc (n, sizeof (double));
    if (a == NULL) {
       fprintf (stderr, "Error:  Out of memory.\n");
       exit (-1);
    }
    return (a);
  }
  
  // Frees an array of doubles.
  static void DestroyArray (const double * a) {
    free ((void *) a);
  }
  
  // Complex number routines.  All outputs must be non-NULL.
  
  // Magnitude/absolute value.
  static double ComplexAbs (const double r, const double i) {
    return (sqrt ((r * r) + (i * i)));
  }
  
  // Multiply.
  static void ComplexMul (const double aR, const double aI, const double bR, const double bI, double * outR, double * outI) {
    (* outR) = (aR * bR) - (aI * bI);
    (* outI) = (aI * bR) + (aR * bI);
  }
  
  // Base-e exponent.
  static void ComplexExp (const double inR, const double inI, double * outR, double * outI) {
    double e;
  
    e = exp (inR);
    (* outR) = e * cos (inI);
    (* outI) = e * sin (inI);
  }
  
  /* Fast Fourier transform routines.  The number of points must be a power of
   * two.  In-place operation is possible only if both the real and imaginary
   * parts are in-place together.
   */
  
  // Performs bit-reversal ordering.
  static void FftArrange (const uint n, const double * inR, const double * inI, double * outR, double * outI) {
    uint rk, k, m;
    double tempR, tempI;
  
    if ((inR == outR) && (inI == outI)) {
       // Handle in-place arrangement.
       rk = 0;
       for (k = 0; k < n; k ++) {
           if (rk > k) {
              tempR = inR [rk];
              tempI = inI [rk];
              outR [rk] = inR [k];
              outI [rk] = inI [k];
              outR [k] = tempR;
              outI [k] = tempI;
           }
           m = n;
           while (rk & (m >>= 1))
             rk &= ~m;
           rk |= m;
       }
    } else {
       // Handle copy arrangement.
       rk = 0;
       for (k = 0; k < n; k ++) {
           outR [rk] = inR [k];
           outI [rk] = inI [k];
           m = n;
           while (rk & (m >>= 1))
             rk &= ~m;
           rk |= m;
       }
    }
  }
  
  // Performs the summation.
  static void FftSummation (const uint n, const double s, double * re, double * im) {
    double pi;
    uint m, m2;
    double vR, vI, wR, wI;
    uint i, k, mk;
    double tR, tI;
  
    pi = s * M_PI;
    for (m = 1, m2 = 2; m < n; m <<= 1, m2 <<= 1) {
        // v = Complex (-2.0 * sin (0.5 * pi / m) * sin (0.5 * pi / m), -sin (pi / m))
        vR = sin (0.5 * pi / m);
        vR = -2.0 * vR * vR;
        vI = -sin (pi / m);
        // w = Complex (1.0, 0.0)
        wR = 1.0;
        wI = 0.0;
        for (i = 0; i < m; i ++) {
            for (k = i; k < n; k += m2) {
                mk = k + m;
                // t = ComplexMul (w, out [km2])
                tR = (wR * re [mk]) - (wI * im [mk]);
                tI = (wR * im [mk]) + (wI * re [mk]);
                // out [mk] = ComplexSub (out [k], t)
                re [mk] = re [k] - tR;
                im [mk] = im [k] - tI;
                // out [k] = ComplexAdd (out [k], t)
                re [k] += tR;
                im [k] += tI;
            }
            // t = ComplexMul (v, w)
            tR = (vR * wR) - (vI * wI);
            tI = (vR * wI) + (vI * wR);
            // w = ComplexAdd (w, t)
            wR += tR;
            wI += tI;
        }
    }
  }
  
  // Performs a forward FFT.
  static void FftForward (const uint n, const double * inR, const double * inI, double * outR, double * outI) {
    FftArrange (n, inR, inI, outR, outI);
    FftSummation (n, 1.0, outR, outI);
  }
  
  // Performs an inverse FFT.
  static void FftInverse (const uint n, const double * inR, const double * inI, double * outR, double * outI) {
    double f;
    uint i;
  
    FftArrange (n, inR, inI, outR, outI);
    FftSummation (n, -1.0, outR, outI);
    f = 1.0 / n;
    for (i = 0; i < n; i ++) {
        outR [i] *= f;
        outI [i] *= f;
    }
  }
  
  /* Calculate the complex helical sequence (or discrete-time analytical
   * signal) of the given input using the Hilbert transform.  Given the
   * negative natural logarithm of a signal's magnitude response, the imaginary
   * components can be used as the angles for minimum-phase reconstruction.
   */
  static void Hilbert (const uint n, const double * in, double * outR, double * outI) {
    uint i;
  
    if (in == outR) {
       // Handle in-place operation.
       for (i = 0; i < n; i ++)
           outI [i] = 0.0;
    } else {
       // Handle copy operation.
       for (i = 0; i < n; i ++) {
           outR [i] = in [i];
           outI [i] = 0.0;
       }
    }
    FftForward (n, outR, outI, outR, outI);
    /* Currently the Fourier routines operate only on point counts that are
     * powers of two.  If that changes and n is odd, the following conditional
     * should be:  i < (n + 1) / 2.
     */
    for (i = 1; i < (n / 2); i ++) {
        outR [i] *= 2.0;
        outI [i] *= 2.0;
    }
    // If n is odd, the following increment should be skipped.
    i ++;
    for (; i < n; i ++) {
        outR [i] = 0.0;
        outI [i] = 0.0;
    }
    FftInverse (n, outR, outI, outR, outI);
  }
  
  /* Calculate the magnitude response of the given input.  This is used in
   * place of phase decomposition, since the phase residuals are discarded for
   * minimum phase reconstruction.  The mirrored half of the response is also
   * discarded.
   */
  static void MagnitudeResponse (const uint n, const double * inR, const double * inI, double * out) {
    const uint m = 1 + (n / 2);
    uint i;
  
    for (i = 0; i < m; i ++)
        out [i] = fmax (ComplexAbs (inR [i], inI [i]), EPSILON);
  }
  
  /* Apply a range limit (in dB) to the given magnitude response.  This is used
   * to adjust the effects of the diffuse-field average on the equalization
   * process.
   */
  static void LimitMagnitudeResponse (const uint n, const double limit, const double * in, double * out) {
    const uint m = 1 + (n / 2);
    double halfLim;
    uint i, lower, upper;
    double ave;
  
    halfLim = limit / 2.0;
    // Convert the response to dB.
    for (i = 0; i < m; i ++)
        out [i] = 20.0 * log10 (in [i]);
    // Use six octaves to calculate the average magnitude of the signal.
    lower = ((uint) ceil (n / pow (2.0, 8.0))) - 1;
    upper = ((uint) floor (n / pow (2.0, 2.0))) - 1;
    ave = 0.0;
    for (i = lower; i <= upper; i ++)
        ave += out [i];
    ave /= upper - lower + 1;
    // Keep the response within range of the average magnitude.
    for (i = 0; i < m; i ++)
        out [i] = Clamp (out [i], ave - halfLim, ave + halfLim);
    // Convert the response back to linear magnitude.
    for (i = 0; i < m; i ++)
        out [i] = pow (10.0, out [i] / 20.0);
  }
  
  /* Reconstructs the minimum-phase component for the given magnitude response
   * of a signal.  This is equivalent to phase recomposition, sans the missing
   * residuals (which were discarded).  The mirrored half of the response is
   * reconstructed.
   */
  static void MinimumPhase (const uint n, const double * in, double * outR, double * outI) {
    const uint m = 1 + (n / 2);
    double * mags = NULL;
    uint i;
    double aR, aI;
  
    mags = CreateArray (n);
    for (i = 0; i < m; i ++) {
        mags [i] = fmax (in [i], EPSILON);
        outR [i] = -log (mags [i]);
    }
    for (; i < n; i ++) {
        mags [i] = mags [n - i];
        outR [i] = outR [n - i];
    }
    Hilbert (n, outR, outR, outI);
    // Remove any DC offset the filter has.
    outR [0] = 0.0;
    outI [0] = 0.0;
    for (i = 1; i < n; i ++) {
        ComplexExp (0.0, outI [i], & aR, & aI);
        ComplexMul (mags [i], 0.0, aR, aI, & outR [i], & outI [i]);
    }
    DestroyArray (mags);
  }
  
  /* This is the normalized cardinal sine (sinc) function.
   *
   *   sinc(x) = { 1,                   x = 0
   *             { sin(pi x) / (pi x),  otherwise.
   */
  static double Sinc (const double x) {
    if (fabs (x) < EPSILON)
       return (1.0);
    return (sin (M_PI * x) / (M_PI * x));
  }
  
  /* The zero-order modified Bessel function of the first kind, used for the
   * Kaiser window.
   *
   *   I_0(x) = sum_{k=0}^inf (1 / k!)^2 (x / 2)^(2 k)
   *          = sum_{k=0}^inf ((x / 2)^k / k!)^2
   */
  static double BesselI_0 (const double x) {
    double term, sum, x2, y, last_sum;
    int k;
  
    // Start at k=1 since k=0 is trivial.
    term = 1.0;
    sum = 1.0;
    x2 = x / 2.0;
    k = 1;
    // Let the integration converge until the term of the sum is no longer
    // significant.
    do {
      y = x2 / k;
      k ++;
      last_sum = sum;
      term *= y * y;
      sum += term;
    } while (sum != last_sum);
    return (sum);
  }
  
  /* Calculate a Kaiser window from the given beta value and a normalized k
   * [-1, 1].
   *
   *   w(k) = { I_0(B sqrt(1 - k^2)) / I_0(B),  -1 <= k <= 1
   *          { 0,                              elsewhere.
   *
   * Where k can be calculated as:
   *
   *   k = i / l,         where -l <= i <= l.
   *
   * or:
   *
   *   k = 2 i / M - 1,   where 0 <= i <= M.
   */
  static double Kaiser (const double b, const double k) {
    double k2;
  
    k2 = Clamp (k, -1.0, 1.0);
    if ((k < -1.0) || (k > 1.0))
       return (0.0);
    k2 *= k2;
    return (BesselI_0 (b * sqrt (1.0 - k2)) / BesselI_0 (b));
  }
  
  // Calculates the greatest common divisor of a and b.
  static uint Gcd (const uint a, const uint b) {
    uint x, y, z;
  
    x = a;
    y = b;
    while (y > 0) {
      z = y;
      y = x % y;
      x = z;
    }
    return (x);
  }
  
  /* Calculates the size (order) of the Kaiser window.  Rejection is in dB and
   * the transition width is normalized frequency (0.5 is nyquist).
   *
   *   M = { ceil((r - 7.95) / (2.285 2 pi f_t)),  r > 21
   *       { ceil(5.79 / 2 pi f_t),                r <= 21.
   *
   */
  static uint CalcKaiserOrder (const double rejection, const double transition) {
    double w_t;
  
    w_t = 2.0 * M_PI * transition;
    if (rejection > 21.0)
       return ((uint) ceil ((rejection - 7.95) / (2.285 * w_t)));
    return ((uint) ceil (5.79 / w_t));
  }
  
  // Calculates the beta value of the Kaiser window.  Rejection is in dB.
  static double CalcKaiserBeta (const double rejection) {
    if (rejection > 50.0)
       return (0.1102 * (rejection - 8.7));
    else if (rejection >= 21.0)
       return ((0.5842 * pow (rejection - 21.0, 0.4)) +
               (0.07886 * (rejection - 21.0)));
    else
       return (0.0);
  }
  
  /* Calculates a point on the Kaiser-windowed sinc filter for the given half-
   * width, beta, gain, and cutoff.  The point is specified in non-normalized
   * samples, from 0 to M, where M = (2 l + 1).
   *
   *   w(k) 2 p f_t sinc(2 f_t x)
   *
   *   x    -- centered sample index (i - l)
   *   k    -- normalized and centered window index (x / l)
   *   w(k) -- window function (Kaiser)
   *   p    -- gain compensation factor when sampling
   *   f_t  -- normalized center frequency (or cutoff; 0.5 is nyquist)
   */
  static double SincFilter (const int l, const double b, const double gain, const double cutoff, const int i) {
    return (Kaiser (b, ((double) (i - l)) / l) * 2.0 * gain * cutoff * Sinc (2.0 * cutoff * (i - l)));
  }
  
  /* This is a polyphase sinc-filtered resampler.
   *
   *              Upsample                      Downsample
   *
   *              p/q = 3/2                     p/q = 3/5
   *
   *          M-+-+-+->                     M-+-+-+->
   *         -------------------+          ---------------------+
   *   p  s * f f f f|f|        |    p  s * f f f f f           |
   *   |  0 *   0 0 0|0|0       |    |  0 *   0 0 0 0|0|        |
   *   v  0 *     0 0|0|0 0     |    v  0 *     0 0 0|0|0       |
   *      s *       f|f|f f f   |       s *       f f|f|f f     |
   *      0 *        |0|0 0 0 0 |       0 *         0|0|0 0 0   |
   *         --------+=+--------+       0 *          |0|0 0 0 0 |
   *          d . d .|d|. d . d            ----------+=+--------+
   *                                        d . . . .|d|. . . .
   *          q->
   *                                        q-+-+-+->
   *
   *   P_f(i,j) = q i mod p + pj
   *   P_s(i,j) = floor(q i / p) - j
   *   d[i=0..N-1] = sum_{j=0}^{floor((M - 1) / p)} {
   *                   { f[P_f(i,j)] s[P_s(i,j)],  P_f(i,j) < M
   *                   { 0,                        P_f(i,j) >= M. }
   */
  
  // Calculate the resampling metrics and build the Kaiser-windowed sinc filter
  // that's used to cut frequencies above the destination nyquist.
  static void ResamplerSetup (ResamplerT * rs, const uint srcRate, const uint dstRate) {
    uint gcd, l;
    double cutoff, width, beta;
    int i;
  
    gcd = Gcd (srcRate, dstRate);
    rs -> mP = dstRate / gcd;
    rs -> mQ = srcRate / gcd;
    /* The cutoff is adjusted by half the transition width, so the transition
     * ends before the nyquist (0.5).  Both are scaled by the downsampling
     * factor.
     */
    if (rs -> mP > rs -> mQ) {
       cutoff = 0.45 / rs -> mP;
       width = 0.1 / rs -> mP;
    } else {
       cutoff = 0.45 / rs -> mQ;
       width = 0.1 / rs -> mQ;
    }
    // A rejection of -180 dB is used for the stop band.
    l = CalcKaiserOrder (180.0, width) / 2;
    beta = CalcKaiserBeta (180.0);
    rs -> mM = (2 * l) + 1;
    rs -> mL = l;
    rs -> mF = CreateArray (rs -> mM);
    for (i = 0; i < ((int) rs -> mM); i ++)
        rs -> mF [i] = SincFilter ((int) l, beta, rs -> mP, cutoff, i);
  }
  
  // Clean up after the resampler.
  static void ResamplerClear (ResamplerT * rs) {
    DestroyArray (rs -> mF);
    rs -> mF = NULL;
  }
  
  // Perform the upsample-filter-downsample resampling operation using a
  // polyphase filter implementation.
  static void ResamplerRun (ResamplerT * rs, const uint inN, const double * in, const uint outN, double * out) {
    const uint p = rs -> mP, q = rs -> mQ, m = rs -> mM, l = rs -> mL;
    const double * f = rs -> mF;
    double * work = NULL;
    uint i;
    double r;
    uint j_f, j_s;
  
    // Handle in-place operation.
    if (in == out)
       work = CreateArray (outN);
    else
       work = out;
    // Resample the input.
    for (i = 0; i < outN; i ++) {
        r = 0.0;
        // Input starts at l to compensate for the filter delay.  This will
        // drop any build-up from the first half of the filter.
        j_f = (l + (q * i)) % p;
        j_s = (l + (q * i)) / p;
        while (j_f < m) {
          // Only take input when 0 <= j_s < inN.  This single unsigned
          // comparison catches both cases.
          if (j_s < inN)
             r += f [j_f] * in [j_s];
          j_f += p;
          j_s --;
        }
        work [i] = r;
    }
    // Clean up after in-place operation.
    if (in == out) {
       for (i = 0; i < outN; i ++)
           out [i] = work [i];
       DestroyArray (work);
    }
  }
  
  // Read a binary value of the specified byte order and byte size from a file,
  // storing it as a 32-bit unsigned integer.
  static int ReadBin4 (FILE * fp, const char * filename, const ByteOrderT order, const uint bytes, uint4 * out) {
    uint1 in [4];
    uint4 accum;
    uint i;
  
    if (fread (in, 1, bytes, fp) != bytes) {
       fprintf (stderr, "Error:  Bad read from file '%s'.\n", filename);
       return (0);
    }
    accum = 0;
    switch (order) {
      case BO_LITTLE :
        for (i = 0; i < bytes; i ++)
            accum = (accum << 8) | in [bytes - i - 1];
      break;
      case BO_BIG :
        for (i = 0; i < bytes; i ++)
            accum = (accum << 8) | in [i];
      break;
      default :
      break;
    }
    (* out) = accum;
    return (1);
  }
  
  // Read a binary value of the specified byte order from a file, storing it as
  // a 64-bit unsigned integer.
  static int ReadBin8 (FILE * fp, const char * filename, const ByteOrderT order, uint8 * out) {
    uint1 in [8];
    uint8 accum;
    uint i;
  
    if (fread (in, 1, 8, fp) != 8) {
       fprintf (stderr, "Error:  Bad read from file '%s'.\n", filename);
       return (0);
    }
    accum = 0ULL;
    switch (order) {
      case BO_LITTLE :
        for (i = 0; i < 8; i ++)
            accum = (accum << 8) | in [8 - i - 1];
      break;
      case BO_BIG :
        for (i = 0; i < 8; i ++)
            accum = (accum << 8) | in [i];
      break;
      default :
      break;
    }
    (* out) = accum;
    return (1);
  }
  
  // Write an ASCII string to a file.
  static int WriteAscii (const char * out, FILE * fp, const char * filename) {
    size_t len;
  
    len = strlen (out);
    if (fwrite (out, 1, len, fp) != len) {
       fclose (fp);
       fprintf (stderr, "Error:  Bad write to file '%s'.\n", filename);
       return (0);
    }
    return (1);
  }
  
  // Write a binary value of the given byte order and byte size to a file,
  // loading it from a 32-bit unsigned integer.
  static int WriteBin4 (const ByteOrderT order, const uint bytes, const uint4 in, FILE * fp, const char * filename) {
    uint1 out [4];
    uint i;
  
    switch (order) {
      case BO_LITTLE :
        for (i = 0; i < bytes; i ++)
            out [i] = (in >> (i * 8)) & 0x000000FF;
      break;
      case BO_BIG :
        for (i = 0; i < bytes; i ++)
            out [bytes - i - 1] = (in >> (i * 8)) & 0x000000FF;
      break;
      default :
      break;
    }
    if (fwrite (out, 1, bytes, fp) != bytes) {
       fprintf (stderr, "Error:  Bad write to file '%s'.\n", filename);
       return (0);
    }
    return (1);
  }
  
  /* Read a binary value of the specified type, byte order, and byte size from
   * a file, converting it to a double.  For integer types, the significant
   * bits are used to normalize the result.  The sign of bits determines
   * whether they are padded toward the MSB (negative) or LSB (positive).
   * Floating-point types are not normalized.
   */
  static int ReadBinAsDouble (FILE * fp, const char * filename, const ByteOrderT order, const ElementTypeT type, const uint bytes, const int bits, double * out) {
    union {
      uint4 ui;
      int4 i;
      float f;
    } v4;
    union {
      uint8 ui;
      double f;
    } v8;
  
    (* out) = 0.0;
    if (bytes > 4) {
       if (! ReadBin8 (fp, filename, order, & v8 . ui))
          return (0);
       if (type == ET_FP)
          (* out) = v8 . f;
    } else {
       if (! ReadBin4 (fp, filename, order, bytes, & v4 . ui))
          return (0);
       if (type == ET_FP) {
          (* out) = (double) v4 . f;
       } else {
          if (bits > 0)
             v4 . ui >>= (8 * bytes) - ((uint) bits);
          else
             v4 . ui &= (0xFFFFFFFF >> (32 + bits));
          if (v4 . ui & ((uint) (1 << (abs (bits) - 1))))
             v4 . ui |= (0xFFFFFFFF << abs (bits));
          (* out) = v4 . i / ((double) (1 << (abs (bits) - 1)));
       }
    }
    return (1);
  }
  
  /* Read an ascii value of the specified type from a file, converting it to a
   * double.  For integer types, the significant bits are used to normalize the
   * result.  The sign of the bits should always be positive.  This also skips
   * up to one separator character before the element itself.
   */
  static int ReadAsciiAsDouble (TokenReaderT * tr, const char * filename, const ElementTypeT type, const uint bits, double * out) {
    int v;
  
    if (TrIsOperator (tr, ","))
       TrReadOperator (tr, ",");
    else if (TrIsOperator (tr, ":"))
       TrReadOperator (tr, ":");
    else if (TrIsOperator (tr, ";"))
       TrReadOperator (tr, ";");
    else if (TrIsOperator (tr, "|"))
       TrReadOperator (tr, "|");
    if (type == ET_FP) {
       if (! TrReadFloat (tr, -HUGE_VAL, HUGE_VAL, out)) {
          fprintf (stderr, "Error:  Bad read from file '%s'.\n", filename);
          return (0);
       }
    } else {
       if (! TrReadInt (tr, -(1 << (bits - 1)), (1 << (bits - 1)) - 1, & v)) {
          fprintf (stderr, "Error:  Bad read from file '%s'.\n", filename);
          return (0);
       }
       (* out) = v / ((double) ((1 << (bits - 1)) - 1));
    }
    return (1);
  }
  
  // Read the RIFF/RIFX WAVE format chunk from a file, validating it against
  // the source parameters and data set metrics.
  static int ReadWaveFormat (FILE * fp, const ByteOrderT order, const uint hrirRate, SourceRefT * src) {
    uint4 fourCC, chunkSize;
    uint4 format, channels, rate, dummy, block, size, bits;
  
    chunkSize = 0;
    do {
      if (chunkSize > 0)
         fseek (fp, (long) chunkSize, SEEK_CUR);
      if ((! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC)) ||
          (! ReadBin4 (fp, src -> mPath, order, 4, & chunkSize)))
         return (0);
    } while (fourCC != FOURCC_FMT);
    if ((! ReadBin4 (fp, src -> mPath, order, 2, & format)) ||
        (! ReadBin4 (fp, src -> mPath, order, 2, & channels)) ||
        (! ReadBin4 (fp, src -> mPath, order, 4, & rate)) ||
        (! ReadBin4 (fp, src -> mPath, order, 4, & dummy)) ||
        (! ReadBin4 (fp, src -> mPath, order, 2, & block)))
       return (0);
    block /= channels;
    if (chunkSize > 14) {
       if (! ReadBin4 (fp, src -> mPath, order, 2, & size))
          return (0);
       size /= 8;
       if (block > size)
          size = block;
    } else {
       size = block;
    }
    if (format == WAVE_FORMAT_EXTENSIBLE) {
       fseek (fp, 2, SEEK_CUR);
       if (! ReadBin4 (fp, src -> mPath, order, 2, & bits))
          return (0);
       if (bits == 0)
          bits = 8 * size;
       fseek (fp, 4, SEEK_CUR);
       if (! ReadBin4 (fp, src -> mPath, order, 2, & format))
          return (0);
       fseek (fp, (long) (chunkSize - 26), SEEK_CUR);
    } else {
       bits = 8 * size;
       if (chunkSize > 14)
          fseek (fp, (long) (chunkSize - 16), SEEK_CUR);
       else
          fseek (fp, (long) (chunkSize - 14), SEEK_CUR);
    }
    if ((format != WAVE_FORMAT_PCM) && (format != WAVE_FORMAT_IEEE_FLOAT)) {
       fprintf (stderr, "Error:  Unsupported WAVE format in file '%s'.\n", src -> mPath);
       return (0);
    }
    if (src -> mChannel >= channels) {
       fprintf (stderr, "Error:  Missing source channel in WAVE file '%s'.\n", src -> mPath);
       return (0);
    }
    if (rate != hrirRate) {
       fprintf (stderr, "Error:  Mismatched source sample rate in WAVE file '%s'.\n", src -> mPath);
       return (0);
    }
    if (format == WAVE_FORMAT_PCM) {
       if ((size < 2) || (size > 4)) {
          fprintf (stderr, "Error:  Unsupported sample size in WAVE file '%s'.\n", src -> mPath);
          return (0);
       }
       if ((bits < 16) || (bits > (8 * size))) {
          fprintf (stderr, "Error:  Bad significant bits in WAVE file '%s'.\n", src -> mPath);
          return (0);
       }
       src -> mType = ET_INT;
    } else {
       if ((size != 4) && (size != 8)) {
          fprintf (stderr, "Error:  Unsupported sample size in WAVE file '%s'.\n", src -> mPath);
          return (0);
       }
       src -> mType = ET_FP;
    }
    src -> mSize = size;
    src -> mBits = (int) bits;
    src -> mSkip = channels;
    return (1);
  }
  
  // Read a RIFF/RIFX WAVE data chunk, converting all elements to doubles.
  static int ReadWaveData (FILE * fp, const SourceRefT * src, const ByteOrderT order, const uint n, double * hrir) {
    int pre, post, skip;
    uint i;
  
    pre = (int) (src -> mSize * src -> mChannel);
    post = (int) (src -> mSize * (src -> mSkip - src -> mChannel - 1));
    skip = 0;
    for (i = 0; i < n; i ++) {
        skip += pre;
        if (skip > 0)
           fseek (fp, skip, SEEK_CUR);
        if (! ReadBinAsDouble (fp, src -> mPath, order, src -> mType, src -> mSize, src -> mBits, & hrir [i]))
           return (0);
        skip = post;
    }
    if (skip > 0)
       fseek (fp, skip, SEEK_CUR);
    return (1);
  }
  
  // Read the RIFF/RIFX WAVE list or data chunk, converting all elements to
  // doubles.
  static int ReadWaveList (FILE * fp, const SourceRefT * src, const ByteOrderT order, const uint n, double * hrir) {
    uint4 fourCC, chunkSize, listSize, count;
    uint block, skip, offset, i;
    double lastSample;
  
    for (;;) {
        if ((! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC)) ||
            (! ReadBin4 (fp, src -> mPath, order, 4, & chunkSize)))
           return (0);
        if (fourCC == FOURCC_DATA) {
           block = src -> mSize * src -> mSkip;
           count = chunkSize / block;
           if (count < (src -> mOffset + n)) {
              fprintf (stderr, "Error:  Bad read from file '%s'.\n", src -> mPath);
              return (0);
           }
           fseek (fp, (long) (src -> mOffset * block), SEEK_CUR);
           if (! ReadWaveData (fp, src, order, n, & hrir [0]))
              return (0);
           return (1);
        } else if (fourCC == FOURCC_LIST) {
           if (! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC))
              return (0);
           chunkSize -= 4;
           if (fourCC == FOURCC_WAVL)
              break;
        }
        if (chunkSize > 0)
           fseek (fp, (long) chunkSize, SEEK_CUR);
    }
    listSize = chunkSize;
    block = src -> mSize * src -> mSkip;
    skip = src -> mOffset;
    offset = 0;
    lastSample = 0.0;
    while ((offset < n) && (listSize > 8)) {
      if ((! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC)) ||
          (! ReadBin4 (fp, src -> mPath, order, 4, & chunkSize)))
         return (0);
      listSize -= 8 + chunkSize;
      if (fourCC == FOURCC_DATA) {
         count = chunkSize / block;
         if (count > skip) {
            fseek (fp, (long) (skip * block), SEEK_CUR);
            chunkSize -= skip * block;
            count -= skip;
            skip = 0;
            if (count > (n - offset))
               count = n - offset;
            if (! ReadWaveData (fp, src, order, count, & hrir [offset]))
               return (0);
            chunkSize -= count * block;
            offset += count;
            lastSample = hrir [offset - 1];
         } else {
            skip -= count;
            count = 0;
         }
      } else if (fourCC == FOURCC_SLNT) {
         if (! ReadBin4 (fp, src -> mPath, order, 4, & count))
            return (0);
         chunkSize -= 4;
         if (count > skip) {
            count -= skip;
            skip = 0;
            if (count > (n - offset))
               count = n - offset;
            for (i = 0; i < count; i ++)
                hrir [offset + i] = lastSample;
            offset += count;
         } else {
            skip -= count;
            count = 0;
         }
      }
      if (chunkSize > 0)
         fseek (fp, (long) chunkSize, SEEK_CUR);
    }
    if (offset < n) {
       fprintf (stderr, "Error:  Bad read from file '%s'.\n", src -> mPath);
       return (0);
    }
    return (1);
  }
  
  // Load a source HRIR from a RIFF/RIFX WAVE file.
  static int LoadWaveSource (FILE * fp, SourceRefT * src, const uint hrirRate, const uint n, double * hrir) {
    uint4 fourCC, dummy;
    ByteOrderT order;
  
    if ((! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC)) ||
        (! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & dummy)))
       return (0);
    if (fourCC == FOURCC_RIFF) {
       order = BO_LITTLE;
    } else if (fourCC == FOURCC_RIFX) {
       order = BO_BIG;
    } else {
       fprintf (stderr, "Error:  No RIFF/RIFX chunk in file '%s'.\n", src -> mPath);
       return (0);
    }
    if (! ReadBin4 (fp, src -> mPath, BO_LITTLE, 4, & fourCC))
       return (0);
    if (fourCC != FOURCC_WAVE) {
       fprintf (stderr, "Error:  Not a RIFF/RIFX WAVE file '%s'.\n", src -> mPath);
       return (0);
    }
    if (! ReadWaveFormat (fp, order, hrirRate, src))
       return (0);
    if (! ReadWaveList (fp, src, order, n, hrir))
       return (0);
    return (1);
  }
  
  // Load a source HRIR from a binary file.
  static int LoadBinarySource (FILE * fp, const SourceRefT * src, const ByteOrderT order, const uint n, double * hrir) {
    uint i;
  
    fseek (fp, (long) src -> mOffset, SEEK_SET);
    for (i = 0; i < n; i ++) {
        if (! ReadBinAsDouble (fp, src -> mPath, order, src -> mType, src -> mSize, src -> mBits, & hrir [i]))
           return (0);
        if (src -> mSkip > 0)
           fseek (fp, (long) src -> mSkip, SEEK_CUR);
    }
    return (1);
  }
  
  // Load a source HRIR from an ASCII text file containing a list of elements
  // separated by whitespace or common list operators (',', ';', ':', '|').
  static int LoadAsciiSource (FILE * fp, const SourceRefT * src, const uint n, double * hrir) {
    TokenReaderT tr;
    uint i, j;
    double dummy;
  
    TrSetup (fp, NULL, & tr);
    for (i = 0; i < src -> mOffset; i ++) {
        if (! ReadAsciiAsDouble (& tr, src -> mPath, src -> mType, (uint) src -> mBits, & dummy))
           return (0);
    }
    for (i = 0; i < n; i ++) {
        if (! ReadAsciiAsDouble (& tr, src -> mPath, src -> mType, (uint) src -> mBits, & hrir [i]))
           return (0);
        for (j = 0; j < src -> mSkip; j ++) {
            if (! ReadAsciiAsDouble (& tr, src -> mPath, src -> mType, (uint) src -> mBits, & dummy))
               return (0);
        }
    }
    return (1);
  }
  
  // Load a source HRIR from a supported file type.
  static int LoadSource (SourceRefT * src, const uint hrirRate, const uint n, double * hrir) {
    FILE * fp = NULL;
    int result;
  
    if (src -> mFormat == SF_ASCII)
       fp = fopen (src -> mPath, "r");
    else
       fp = fopen (src -> mPath, "rb");
    if (fp == NULL) {
       fprintf (stderr, "Error:  Could not open source file '%s'.\n", src -> mPath);
       return (0);
    }
    if (src -> mFormat == SF_WAVE)
       result = LoadWaveSource (fp, src, hrirRate, n, hrir);
    else if (src -> mFormat == SF_BIN_LE)
       result = LoadBinarySource (fp, src, BO_LITTLE, n, hrir);
    else if (src -> mFormat == SF_BIN_BE)
       result = LoadBinarySource (fp, src, BO_BIG, n, hrir);
    else
       result = LoadAsciiSource (fp, src, n, hrir);
    fclose (fp);
    return (result);
  }
  
  // Calculate the magnitude response of an HRIR and average it with any
  // existing responses for its elevation and azimuth.
  static void AverageHrirMagnitude (const double * hrir, const double f, const uint ei, const uint ai, const HrirDataT * hData) {
    double * re = NULL, * im = NULL;
    uint n, m, i, j;
  
    n = hData -> mFftSize;
    re = CreateArray (n);
    im = CreateArray (n);
    for (i = 0; i < hData -> mIrPoints; i ++) {
        re [i] = hrir [i];
        im [i] = 0.0;
    }
    for (; i < n; i ++) {
        re [i] = 0.0;
        im [i] = 0.0;
    }
    FftForward (n, re, im, re, im);
    MagnitudeResponse (n, re, im, re);
    m = 1 + (n / 2);
    j = (hData -> mEvOffset [ei] + ai) * hData -> mIrSize;
    for (i = 0; i < m; i ++)
        hData -> mHrirs [j + i] = Lerp (hData -> mHrirs [j + i], re [i], f);
    DestroyArray (im);
    DestroyArray (re);
  }
  
  /* Calculate the contribution of each HRIR to the diffuse-field average based
   * on the area of its surface patch.  All patches are centered at the HRIR
   * coordinates on the unit sphere and are measured by solid angle.
   */
  static void CalculateDfWeights (const HrirDataT * hData, double * weights) {
    uint ei;
    double evs, sum, ev, up_ev, down_ev, solidAngle;
  
    evs = 90.0 / (hData -> mEvCount - 1);
    sum = 0.0;
    for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++) {
        // For each elevation, calculate the upper and lower limits of the
        // patch band.
        ev = -90.0 + (ei * 2.0 * evs);
        if (ei < (hData -> mEvCount - 1))
           up_ev = (ev + evs) * M_PI / 180.0;
        else
           up_ev = M_PI / 2.0;
        if (ei > 0)
           down_ev = (ev - evs) * M_PI / 180.0;
        else
           down_ev = -M_PI / 2.0;
        // Calculate the area of the patch band.
        solidAngle = 2.0 * M_PI * (sin (up_ev) - sin (down_ev));
        // Each weight is the area of one patch.
        weights [ei] = solidAngle / hData -> mAzCount [ei];
        // Sum the total surface area covered by the HRIRs.
        sum += solidAngle;
    }
    // Normalize the weights given the total surface coverage.
    for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
        weights [ei] /= sum;
  }
  
  /* Calculate the diffuse-field average from the given magnitude responses of
   * the HRIR set.  Weighting can be applied to compensate for the varying
   * surface area covered by each HRIR.  The final average can then be limited
   * by the specified magnitude range (in positive dB; 0.0 to skip).
   */
  static void CalculateDiffuseFieldAverage (const HrirDataT * hData, const int weighted, const double limit, double * dfa) {
    double * weights = NULL;
    uint ei, ai, count, step, start, end, m, j, i;
    double weight;
  
    weights = CreateArray (hData -> mEvCount);
    if (weighted) {
       // Use coverage weighting to calculate the average.
       CalculateDfWeights (hData, weights);
    } else {
       // If coverage weighting is not used, the weights still need to be
       // averaged by the number of HRIRs.
       count = 0;
       for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
           count += hData -> mAzCount [ei];
       for (ei = hData -> mEvStart; ei < hData -> mEvCount; ei ++)
           weights [ei] = 1.0 / count;
    }
    ei = hData -> mEvStart;
    ai = 0;
    step = hData -> mIrSize;
    start = hData -> mEvOffset [ei] * step;
    end = hData -> mIrCount * step;
    m = 1 + (hData -> mFftSize / 2);
    for (i = 0; i < m; i ++)
        dfa [i] = 0.0;
    for (j = start; j < end; j += step) {
        // Get the weight for this HRIR's contribution.
        weight = weights [ei];
        // Add this HRIR's weighted power average to the total.
        for (i = 0; i < m; i ++)
            dfa [i] += weight * hData -> mHrirs [j + i] * hData -> mHrirs [j + i];
        // Determine the next weight to use.
        ai ++;
        if (ai >= hData -> mAzCount [ei]) {
           ei ++;
           ai = 0;
        }
    }
    // Finish the average calculation and keep it from being too small.
    for (i = 0; i < m; i ++)
        dfa [i] = fmax (sqrt (dfa [i]), EPSILON);
    // Apply a limit to the magnitude range of the diffuse-field average if
    // desired.
    if (limit > 0.0)
       LimitMagnitudeResponse (hData -> mFftSize, limit, dfa, dfa);
    DestroyArray (weights);
  }
  
  // Perform diffuse-field equalization on the magnitude responses of the HRIR
  // set using the given average response.
  static void DiffuseFieldEqualize (const double * dfa, const HrirDataT * hData) {
    uint step, start, end, m, j, i;
  
    step = hData -> mIrSize;
    start = hData -> mEvOffset [hData -> mEvStart] * step;
    end = hData -> mIrCount * step;
    m = 1 + (hData -> mFftSize / 2);
    for (j = start; j < end; j += step) {
        for (i = 0; i < m; i ++)
            hData -> mHrirs [j + i] /= dfa [i];
    }
  }
  
  // Perform minimum-phase reconstruction using the magnitude responses of the
  // HRIR set.
  static void ReconstructHrirs (const HrirDataT * hData) {
    double * re = NULL, * im = NULL;
    uint step, start, end, n, j, i;
  
    step = hData -> mIrSize;
    start = hData -> mEvOffset [hData -> mEvStart] * step;
    end = hData -> mIrCount * step;
    n = hData -> mFftSize;
    re = CreateArray (n);
    im = CreateArray (n);
    for (j = start; j < end; j += step) {
        MinimumPhase (n, & hData -> mHrirs [j], re, im);
        FftInverse (n, re, im, re, im);
        for (i = 0; i < hData -> mIrPoints; i ++)
            hData -> mHrirs [j + i] = re [i];
    }
    DestroyArray (im);
    DestroyArray (re);
  }
  
  // Resamples the HRIRs for use at the given sampling rate.
  static void ResampleHrirs (const uint rate, HrirDataT * hData) {
    ResamplerT rs;
    uint n, step, start, end, j;
  
    ResamplerSetup (& rs, hData -> mIrRate, rate);
    n = hData -> mIrPoints;
    step = hData -> mIrSize;
    start = hData -> mEvOffset [hData -> mEvStart] * step;
    end = hData -> mIrCount * step;
    for (j = start; j < end; j += step)
        ResamplerRun (& rs, n, & hData -> mHrirs [j], n, & hData -> mHrirs [j]);
    ResamplerClear (& rs);
    hData -> mIrRate = rate;
  }
  
  /* Given an elevation index and an azimuth, calculate the indices of the two
   * HRIRs that bound the coordinate along with a factor for calculating the
   * continous HRIR using interpolation.
   */
  static void CalcAzIndices (const HrirDataT * hData, const uint ei, const double az, uint * j0, uint * j1, double * jf) {
    double af;
    uint ai;
  
    af = ((2.0 * M_PI) + az) * hData -> mAzCount [ei] / (2.0 * M_PI);
    ai = ((uint) af) % hData -> mAzCount [ei];
    af -= floor (af);
    (* j0) = hData -> mEvOffset [ei] + ai;
    (* j1) = hData -> mEvOffset [ei] + ((ai + 1) % hData -> mAzCount [ei]);
    (* jf) = af;
  }
  
  /* Attempt to synthesize any missing HRIRs at the bottom elevations.  Right
   * now this just blends the lowest elevation HRIRs together and applies some
   * attenuation and high frequency damping.  It is a simple, if inaccurate
   * model.
   */
  static void SynthesizeHrirs (HrirDataT * hData) {
    uint oi, a, e, step, n, i, j;
    double of, b;
    uint j0, j1;
    double jf;
    double lp [4], s0, s1;
  
    if (hData -> mEvStart <= 0)
       return;
    step = hData -> mIrSize;
    oi = hData -> mEvStart;
    n = hData -> mIrPoints;
    for (i = 0; i < n; i ++)
        hData -> mHrirs [i] = 0.0;
    for (a = 0; a < hData -> mAzCount [oi]; a ++) {
        j = (hData -> mEvOffset [oi] + a) * step;
        for (i = 0; i < n; i ++)
            hData -> mHrirs [i] += hData -> mHrirs [j + i] / hData -> mAzCount [oi];
    }
    for (e = 1; e < hData -> mEvStart; e ++) {
        of = ((double) e) / hData -> mEvStart;
        b = (1.0 - of) * (3.5e-6 * hData -> mIrRate);
        for (a = 0; a < hData -> mAzCount [e]; a ++) {
            j = (hData -> mEvOffset [e] + a) * step;
            CalcAzIndices (hData, oi, a * 2.0 * M_PI / hData -> mAzCount [e], & j0, & j1, & jf);
            j0 *= step;
            j1 *= step;
            lp [0] = 0.0;
            lp [1] = 0.0;
            lp [2] = 0.0;
            lp [3] = 0.0;
            for (i = 0; i < n; i ++) {
                s0 = hData -> mHrirs [i];
                s1 = Lerp (hData -> mHrirs [j0 + i], hData -> mHrirs [j1 + i], jf);
                s0 = Lerp (s0, s1, of);
                lp [0] = Lerp (s0, lp [0], b);
                lp [1] = Lerp (lp [0], lp [1], b);
                lp [2] = Lerp (lp [1], lp [2], b);
                lp [3] = Lerp (lp [2], lp [3], b);
                hData -> mHrirs [j + i] = lp [3];
            }
        }
    }
    b = 3.5e-6 * hData -> mIrRate;
    lp [0] = 0.0;
    lp [1] = 0.0;
    lp [2] = 0.0;
    lp [3] = 0.0;
    for (i = 0; i < n; i ++) {
        s0 = hData -> mHrirs [i];
        lp [0] = Lerp (s0, lp [0], b);
        lp [1] = Lerp (lp [0], lp [1], b);
        lp [2] = Lerp (lp [1], lp [2], b);
        lp [3] = Lerp (lp [2], lp [3], b);
        hData -> mHrirs [i] = lp [3];
    }
    hData -> mEvStart = 0;
  }
  
  // The following routines assume a full set of HRIRs for all elevations.
  
  // Normalize the HRIR set and slightly attenuate the result.
  static void NormalizeHrirs (const HrirDataT * hData) {
    uint step, end, n, j, i;
    double maxLevel;
  
    step = hData -> mIrSize;
    end = hData -> mIrCount * step;
    n = hData -> mIrPoints;
    maxLevel = 0.0;
    for (j = 0; j < end; j += step) {
        for (i = 0; i < n; i ++)
            maxLevel = fmax (fabs (hData -> mHrirs [j + i]), maxLevel);
    }
    maxLevel = 1.01 * maxLevel;
    for (j = 0; j < end; j += step) {
        for (i = 0; i < n; i ++)
            hData -> mHrirs [j + i] /= maxLevel;
    }
  }
  
  // Calculate the left-ear time delay using a spherical head model.
  static double CalcLTD (const double ev, const double az, const double rad, const double dist) {
    double azp, dlp, l, al;
  
    azp = asin (cos (ev) * sin (az));
    dlp = sqrt ((dist * dist) + (rad * rad) + (2.0 * dist * rad * sin (azp)));
    l = sqrt ((dist * dist) - (rad * rad));
    al = (0.5 * M_PI) + azp;
    if (dlp > l)
       dlp = l + (rad * (al - acos (rad / dist)));
    return (dlp / 343.3);
  }
  
  // Calculate the effective head-related time delays for the each HRIR, now
  // that they are minimum-phase.
  static void CalculateHrtds (HrirDataT * hData) {
    double minHrtd, maxHrtd;
    uint e, a, j;
    double t;
  
    minHrtd = 1000.0;
    maxHrtd = -1000.0;
    for (e = 0; e < hData -> mEvCount; e ++) {
        for (a = 0; a < hData -> mAzCount [e]; a ++) {
            j = hData -> mEvOffset [e] + a;
            t = CalcLTD ((-90.0 + (e * 180.0 / (hData -> mEvCount - 1))) * M_PI / 180.0,
                         (a * 360.0 / hData -> mAzCount [e]) * M_PI / 180.0,
                         hData -> mRadius, hData -> mDistance);
            hData -> mHrtds [j] = t;
            maxHrtd = fmax (t, maxHrtd);
            minHrtd = fmin (t, minHrtd);
        }
    }
    maxHrtd -= minHrtd;
    for (j = 0; j < hData -> mIrCount; j ++)
        hData -> mHrtds [j] -= minHrtd;
    hData -> mMaxHrtd = maxHrtd;
  }
  
  // Store the OpenAL Soft HRTF data set.
  static int StoreMhr (const HrirDataT * hData, const char * filename) {
    FILE * fp = NULL;
    uint e, step, end, n, j, i;
    int hpHist, v;
  
    if ((fp = fopen (filename, "wb")) == NULL) {
       fprintf (stderr, "Error:  Could not open MHR file '%s'.\n", filename);
       return (0);
    }
    if (! WriteAscii (MHR_FORMAT, fp, filename))
       return (0);
    if (! WriteBin4 (BO_LITTLE, 4, (uint4) hData -> mIrRate, fp, filename))
       return (0);
    if (! WriteBin4 (BO_LITTLE, 1, (uint4) hData -> mIrPoints, fp, filename))
       return (0);
    if (! WriteBin4 (BO_LITTLE, 1, (uint4) hData -> mEvCount, fp, filename))
       return (0);
    for (e = 0; e < hData -> mEvCount; e ++) {
        if (! WriteBin4 (BO_LITTLE, 1, (uint4) hData -> mAzCount [e], fp, filename))
           return (0);
    }
    step = hData -> mIrSize;
    end = hData -> mIrCount * step;
    n = hData -> mIrPoints;
    srand (0x31DF840C);
    for (j = 0; j < end; j += step) {
        hpHist = 0;
        for (i = 0; i < n; i ++) {
            v = HpTpdfDither (32767.0 * hData -> mHrirs [j + i], & hpHist);
            if (! WriteBin4 (BO_LITTLE, 2, (uint4) v, fp, filename))
               return (0);
        }
    }
    for (j = 0; j < hData -> mIrCount; j ++) {
        v = (int) fmin (round (hData -> mIrRate * hData -> mHrtds [j]), MAX_HRTD);
        if (! WriteBin4 (BO_LITTLE, 1, (uint4) v, fp, filename))
           return (0);
    }
    fclose (fp);
    return (1);
  }
  
  // Store the OpenAL Soft built-in table.
  static int StoreTable (const HrirDataT * hData, const char * filename) {
    FILE * fp = NULL;
    uint step, end, n, j, i;
    int hpHist, v;
    char text [128 + 1];
  
    if ((fp = fopen (filename, "wb")) == NULL) {
       fprintf (stderr, "Error:  Could not open table file '%s'.\n", filename);
       return (0);
    }
    snprintf (text, 128, "/* Elevation metrics */\n"
                         "static const ALubyte defaultAzCount[%u] = { ", hData -> mEvCount);
    if (! WriteAscii (text, fp, filename))
       return (0);
    for (i = 0; i < hData -> mEvCount; i ++) {
        snprintf (text, 128, "%u, ", hData -> mAzCount [i]);
        if (! WriteAscii (text, fp, filename))
           return (0);
    }
    snprintf (text, 128, "};\n"
                         "static const ALushort defaultEvOffset[%u] = { ", hData -> mEvCount);
    if (! WriteAscii (text, fp, filename))
       return (0);
    for (i = 0; i < hData -> mEvCount; i ++) {
        snprintf (text, 128, "%u, ", hData -> mEvOffset [i]);
        if (! WriteAscii (text, fp, filename))
           return (0);
    }
    step = hData -> mIrSize;
    end = hData -> mIrCount * step;
    n = hData -> mIrPoints;
    snprintf (text, 128, "};\n\n"
                         "/* HRIR Coefficients */\n"
                         "static const ALshort defaultCoeffs[%u] =\n{\n", hData -> mIrCount * n);
    if (! WriteAscii (text, fp, filename))
       return (0);
    srand (0x31DF840C);
    for (j = 0; j < end; j += step) {
        if (! WriteAscii ("   ", fp, filename))
           return (0);
        hpHist = 0;
        for (i = 0; i < n; i ++) {
            v = HpTpdfDither (32767.0 * hData -> mHrirs [j + i], & hpHist);
            snprintf (text, 128, " %+d,", v);
            if (! WriteAscii (text, fp, filename))
               return (0);
        }
        if (! WriteAscii ("\n", fp, filename))
           return (0);
    }
    snprintf (text, 128, "};\n\n"
                         "/* HRIR Delays */\n"
                         "static const ALubyte defaultDelays[%u] =\n{\n"
                         "   ", hData -> mIrCount);
    if (! WriteAscii (text, fp, filename))
       return (0);
    for (j = 0; j < hData -> mIrCount; j ++) {
        v = (int) fmin (round (hData -> mIrRate * hData -> mHrtds [j]), MAX_HRTD);
        snprintf (text, 128, " %d,", v);
        if (! WriteAscii (text, fp, filename))
           return (0);
    }
    if (! WriteAscii ("\n};\n\n"
                      "/* Default HRTF Definition */\n", fp, filename))
       return (0);
    snprintf (text, 128, "static const struct Hrtf DefaultHrtf = {\n"
                         "    %u, %u, %u, defaultAzCount, defaultEvOffset,\n",
                         hData -> mIrRate, hData -> mIrPoints, hData -> mEvCount);
    if (! WriteAscii (text, fp, filename))
       return (0);
    if (! WriteAscii ("    defaultCoeffs, defaultDelays, NULL\n"
                      "};\n", fp, filename))
       return (0);
    fclose (fp);
    return (1);
  }
  
  // Process the data set definition to read and validate the data set metrics.
  static int ProcessMetrics (TokenReaderT * tr, const uint fftSize, const uint truncSize, HrirDataT * hData) {
    char ident [MAX_IDENT_LEN + 1];
    uint line, col;
    int intVal;
    uint points;
    double fpVal;
    int hasRate = 0, hasPoints = 0, hasAzimuths = 0;
    int hasRadius = 0, hasDistance = 0;
  
    while (! (hasRate && hasPoints && hasAzimuths && hasRadius && hasDistance)) {
      TrIndication (tr, & line, & col);
      if (! TrReadIdent (tr, MAX_IDENT_LEN, ident))
         return (0);
      if (strcasecmp (ident, "rate") == 0) {
         if (hasRate) {
            TrErrorAt (tr, line, col, "Redefinition of 'rate'.\n");
            return (0);
         }
         if (! TrReadOperator (tr, "="))
            return (0);
         if (! TrReadInt (tr, MIN_RATE, MAX_RATE, & intVal))
            return (0);
         hData -> mIrRate = (uint) intVal;
         hasRate = 1;
      } else if (strcasecmp (ident, "points") == 0) {
         if (hasPoints) {
            TrErrorAt (tr, line, col, "Redefinition of 'points'.\n");
            return (0);
         }
         if (! TrReadOperator (tr, "="))
            return (0);
         TrIndication (tr, & line, & col);
         if (! TrReadInt (tr, MIN_POINTS, MAX_POINTS, & intVal))
            return (0);
         points = (uint) intVal;
         if ((fftSize > 0) && (points > fftSize)) {
            TrErrorAt (tr, line, col, "Value exceeds the overriden FFT size.\n");
            return (0);
         }
         if (points < truncSize) {
            TrErrorAt (tr, line, col, "Value is below the truncation size.\n");
            return (0);
         }
         hData -> mIrPoints = points;
         hData -> mFftSize = fftSize;
         if (fftSize <= 0) {
            points = 1;
            while (points < (4 * hData -> mIrPoints))
              points <<= 1;
            hData -> mFftSize = points;
            hData -> mIrSize = 1 + (points / 2);
         } else {
            hData -> mFftSize = fftSize;
            hData -> mIrSize = 1 + (fftSize / 2);
            if (points > hData -> mIrSize)
               hData -> mIrSize = points;
         }
         hasPoints = 1;
      } else if (strcasecmp (ident, "azimuths") == 0) {
         if (hasAzimuths) {
            TrErrorAt (tr, line, col, "Redefinition of 'azimuths'.\n");
            return (0);
         }
         if (! TrReadOperator (tr, "="))
            return (0);
         hData -> mIrCount = 0;
         hData -> mEvCount = 0;
         hData -> mEvOffset [0] = 0;
         for (;;) {
             if (! TrReadInt (tr, MIN_AZ_COUNT, MAX_AZ_COUNT, & intVal))
                return (0);
             hData -> mAzCount [hData -> mEvCount] = (uint) intVal;
             hData -> mIrCount += (uint) intVal;
             hData -> mEvCount ++;
             if (! TrIsOperator (tr, ","))
                break;
             if (hData -> mEvCount >= MAX_EV_COUNT) {
                TrError (tr, "Exceeded the maximum of %d elevations.\n", MAX_EV_COUNT);
                return (0);
             }
             hData -> mEvOffset [hData -> mEvCount] = hData -> mEvOffset [hData -> mEvCount - 1] + ((uint) intVal);
             TrReadOperator (tr, ",");
         }
         if (hData -> mEvCount < MIN_EV_COUNT) {
            TrErrorAt (tr, line, col, "Did not reach the minimum of %d azimuth counts.\n", MIN_EV_COUNT);
            return (0);
         }
         hasAzimuths = 1;
      } else if (strcasecmp (ident, "radius") == 0) {
         if (hasRadius) {
            TrErrorAt (tr, line, col, "Redefinition of 'radius'.\n");
            return (0);
         }
         if (! TrReadOperator (tr, "="))
            return (0);
         if (! TrReadFloat (tr, MIN_RADIUS, MAX_RADIUS, & fpVal))
            return (0);
         hData -> mRadius = fpVal;
         hasRadius = 1;
      } else if (strcasecmp (ident, "distance") == 0) {
         if (hasDistance) {
            TrErrorAt (tr, line, col, "Redefinition of 'distance'.\n");
            return (0);
         }
         if (! TrReadOperator (tr, "="))
            return (0);
         if (! TrReadFloat (tr, MIN_DISTANCE, MAX_DISTANCE, & fpVal))
            return (0);
         hData -> mDistance = fpVal;
         hasDistance = 1;
      } else {
         TrErrorAt (tr, line, col, "Expected a metric name.\n");
         return (0);
      }
      TrSkipWhitespace (tr);
    }
    return (1);
  }
  
  // Parse an index pair from the data set definition.
  static int ReadIndexPair (TokenReaderT * tr, const HrirDataT * hData, uint * ei, uint * ai) {
    int intVal;
  
    if (! TrReadInt (tr, 0, (int) hData -> mEvCount, & intVal))
       return (0);
    (* ei) = (uint) intVal;
    if (! TrReadOperator (tr, ","))
       return (0);
    if (! TrReadInt (tr, 0, (int) hData -> mAzCount [(* ei)], & intVal))
       return (0);
    (* ai) = (uint) intVal;
    return (1);
  }
  
  // Match the source format from a given identifier.
  static SourceFormatT MatchSourceFormat (const char * ident) {
    if (strcasecmp (ident, "wave") == 0)
       return (SF_WAVE);
    else if (strcasecmp (ident, "bin_le") == 0)
       return (SF_BIN_LE);
    else if (strcasecmp (ident, "bin_be") == 0)
       return (SF_BIN_BE);
    else if (strcasecmp (ident, "ascii") == 0)
       return (SF_ASCII);
    return (SF_NONE);
  }
  
  // Match the source element type from a given identifier.
  static ElementTypeT MatchElementType (const char * ident) {
    if (strcasecmp (ident, "int") == 0)
       return (ET_INT);
    else if (strcasecmp (ident, "fp") == 0)
       return (ET_FP);
    return (ET_NONE);
  }
  
  // Parse and validate a source reference from the data set definition.
  static int ReadSourceRef (TokenReaderT * tr, SourceRefT * src) {
    uint line, col;
    char ident [MAX_IDENT_LEN + 1];
    int intVal;
  
    TrIndication (tr, & line, & col);
    if (! TrReadIdent (tr, MAX_IDENT_LEN, ident))
       return (0);
    src -> mFormat = MatchSourceFormat (ident);
    if (src -> mFormat == SF_NONE) {
       TrErrorAt (tr, line, col, "Expected a source format.\n");
       return (0);
    }
    if (! TrReadOperator (tr, "("))
       return (0);
    if (src -> mFormat == SF_WAVE) {
       if (! TrReadInt (tr, 0, MAX_WAVE_CHANNELS, & intVal))
          return (0);
       src -> mType = ET_NONE;
       src -> mSize = 0;
       src -> mBits = 0;
       src -> mChannel = (uint) intVal;
       src -> mSkip = 0;
    } else {
       TrIndication (tr, & line, & col);
       if (! TrReadIdent (tr, MAX_IDENT_LEN, ident))
          return (0);
       src -> mType = MatchElementType (ident);
       if (src -> mType == ET_NONE) {
          TrErrorAt (tr, line, col, "Expected a source element type.\n");
          return (0);
       }
       if ((src -> mFormat == SF_BIN_LE) || (src -> mFormat == SF_BIN_BE)) {
          if (! TrReadOperator (tr, ","))
             return (0);
          if (src -> mType == ET_INT) {
             if (! TrReadInt (tr, MIN_BIN_SIZE, MAX_BIN_SIZE, & intVal))
                return (0);
             src -> mSize = (uint) intVal;
             if (TrIsOperator (tr, ",")) {
                TrReadOperator (tr, ",");
                TrIndication (tr, & line, & col);
                if (! TrReadInt (tr, -2147483647 - 1, 2147483647, & intVal))
                   return (0);
                if ((abs (intVal) < MIN_BIN_BITS) || (((uint) abs (intVal)) > (8 * src -> mSize))) {
                   TrErrorAt (tr, line, col, "Expected a value of (+/-) %d to %d.\n", MIN_BIN_BITS, 8 * src -> mSize);
                   return (0);
                }
                src -> mBits = intVal;
             } else {
                src -> mBits = (int) (8 * src -> mSize);
             }
          } else {
             TrIndication (tr, & line, & col);
             if (! TrReadInt (tr, -2147483647 - 1, 2147483647, & intVal))
                return (0);
             if ((intVal != 4) && (intVal != 8)) {
                TrErrorAt (tr, line, col, "Expected a value of 4 or 8.\n");
                return (0);
             }
             src -> mSize = (uint) intVal;
             src -> mBits = 0;
          }
       } else if ((src -> mFormat == SF_ASCII) && (src -> mType == ET_INT)) {
          if (! TrReadOperator (tr, ","))
             return (0);
          if (! TrReadInt (tr, MIN_ASCII_BITS, MAX_ASCII_BITS, & intVal))
             return (0);
          src -> mSize = 0;
          src -> mBits = intVal;
       } else {
          src -> mSize = 0;
          src -> mBits = 0;
       }
       if (TrIsOperator (tr, ";")) {
          TrReadOperator (tr, ";");
          if (! TrReadInt (tr, 0, 0x7FFFFFFF, & intVal))
             return (0);
          src -> mSkip = (uint) intVal;
       } else {
          src -> mSkip = 0;
       }
    }
    if (! TrReadOperator (tr, ")"))
       return (0);
    if (TrIsOperator (tr, "@")) {
       TrReadOperator (tr, "@");
       if (! TrReadInt (tr, 0, 0x7FFFFFFF, & intVal))
          return (0);
       src -> mOffset = (uint) intVal;
    } else {
       src -> mOffset = 0;
    }
    if (! TrReadOperator (tr, ":"))
       return (0);
    if (! TrReadString (tr, MAX_PATH_LEN, src -> mPath))
       return (0);
    return (1);
  }
  
  // Process the list of sources in the data set definition.
  static int ProcessSources (TokenReaderT * tr, HrirDataT * hData) {
    uint * setCount = NULL, * setFlag = NULL;
    double * hrir = NULL;
    uint line, col, ei, ai;
    SourceRefT src;
    double factor;
  
    setCount = (uint *) calloc (hData -> mEvCount, sizeof (uint));
    setFlag = (uint *) calloc (hData -> mIrCount, sizeof (uint));
    hrir = CreateArray (hData -> mIrPoints);
    while (TrIsOperator (tr, "[")) {
      TrIndication (tr, & line, & col);
      TrReadOperator (tr, "[");
      if (ReadIndexPair (tr, hData, & ei, & ai)) {
         if (TrReadOperator (tr, "]")) {
            if (! setFlag [hData -> mEvOffset [ei] + ai]) {
               if (TrReadOperator (tr, "=")) {
                  factor = 1.0;
                  for (;;) {
                      if (ReadSourceRef (tr, & src)) {
                         if (LoadSource (& src, hData -> mIrRate, hData -> mIrPoints, hrir)) {
                            AverageHrirMagnitude (hrir, 1.0 / factor, ei, ai, hData);
                            factor += 1.0;
                            if (! TrIsOperator (tr, "+"))
                               break;
                            TrReadOperator (tr, "+");
                            continue;
                         }
                      }
                      DestroyArray (hrir);
                      free (setFlag);
                      free (setCount);
                      return (0);
                  }
                  setFlag [hData -> mEvOffset [ei] + ai] = 1;
                  setCount [ei] ++;
                  continue;
               }
            } else {
               TrErrorAt (tr, line, col, "Redefinition of source.\n");
            }
         }
      }
      DestroyArray (hrir);
      free (setFlag);
      free (setCount);
      return (0);
    }
    ei = 0;
    while ((ei < hData -> mEvCount) && (setCount [ei] < 1))
      ei ++;
    if (ei < hData -> mEvCount) {
       hData -> mEvStart = ei;
       while ((ei < hData -> mEvCount) && (setCount [ei] == hData -> mAzCount [ei]))
         ei ++;
       if (ei >= hData -> mEvCount) {
          if (! TrLoad (tr)) {
             DestroyArray (hrir);
             free (setFlag);
             free (setCount);
             return (1);
          } else {
             TrError (tr, "Errant data at end of source list.\n");
          }
       } else {
          TrError (tr, "Missing sources for elevation index %d.\n", ei);
       }
    } else {
       TrError (tr, "Missing source references.\n");
    }
    DestroyArray (hrir);
    free (setFlag);
    free (setCount);
    return (0);
  }
  
  /* Parse the data set definition and process the source data, storing the
   * resulting data set as desired.  If the input name is NULL it will read
   * from standard input.
   */
  static int ProcessDefinition (const char * inName, const uint outRate, const uint fftSize, const int equalize, const int surface, const double limit, const uint truncSize, const OutputFormatT outFormat, const char * outName) {
    FILE * fp = NULL;
    TokenReaderT tr;
    HrirDataT hData;
    double * dfa = NULL;
    char rateStr [8 + 1], expName [MAX_PATH_LEN];
  
    hData . mIrRate = 0;
    hData . mIrPoints = 0;
    hData . mFftSize = 0;
    hData . mIrSize = 0;
    hData . mIrCount = 0;
    hData . mEvCount = 0;
    hData . mRadius = 0;
    hData . mDistance = 0;
    fprintf (stdout, "Reading HRIR definition...\n");
    if (inName != NULL) {
       fp = fopen (inName, "r");
       if (fp == NULL) {
          fprintf (stderr, "Error:  Could not open definition file '%s'\n", inName);
          return (0);
       }
       TrSetup (fp, inName, & tr);
    } else {
       fp = stdin;
       TrSetup (fp, "<stdin>", & tr);
    }
    if (! ProcessMetrics (& tr, fftSize, truncSize, & hData)) {
       if (inName != NULL)
          fclose (fp);
       return (0);
    }
    hData . mHrirs = CreateArray (hData . mIrCount * hData . mIrSize);
    hData . mHrtds = CreateArray (hData . mIrCount);
    if (! ProcessSources (& tr, & hData)) {
       DestroyArray (hData . mHrtds);
       DestroyArray (hData . mHrirs);
       if (inName != NULL)
          fclose (fp);
       return (0);
    }
    if (inName != NULL)
       fclose (fp);
    if (equalize) {
       dfa = CreateArray (1 + (hData . mFftSize / 2));
       fprintf (stdout, "Calculating diffuse-field average...\n");
       CalculateDiffuseFieldAverage (& hData, surface, limit, dfa);
       fprintf (stdout, "Performing diffuse-field equalization...\n");
       DiffuseFieldEqualize (dfa, & hData);
       DestroyArray (dfa);
    }
    fprintf (stdout, "Performing minimum phase reconstruction...\n");
    ReconstructHrirs (& hData);
    if ((outRate != 0) && (outRate != hData . mIrRate)) {
       fprintf (stdout, "Resampling HRIRs...\n");
       ResampleHrirs (outRate, & hData);
    }
    fprintf (stdout, "Truncating minimum-phase HRIRs...\n");
    hData . mIrPoints = truncSize;
    fprintf (stdout, "Synthesizing missing elevations...\n");
    SynthesizeHrirs (& hData);
    fprintf (stdout, "Normalizing final HRIRs...\n");
    NormalizeHrirs (& hData);
    fprintf (stdout, "Calculating impulse delays...\n");
    CalculateHrtds (& hData);
    snprintf (rateStr, 8, "%u", hData . mIrRate);
    StrSubst (outName, "%r", rateStr, MAX_PATH_LEN, expName);
    switch (outFormat) {
      case OF_MHR :
        fprintf (stdout, "Creating MHR data set file...\n");
        if (! StoreMhr (& hData, expName))
           return (0);
      break;
      case OF_TABLE :
        fprintf (stderr, "Creating OpenAL Soft table file...\n");
        if (! StoreTable (& hData, expName))
           return (0);
      break;
      default :
      break;
    }
    DestroyArray (hData . mHrtds);
    DestroyArray (hData . mHrirs);
    return (1);
  }
  
  // Standard command line dispatch.
  int main (const int argc, const char * argv []) {
    const char * inName = NULL, * outName = NULL;
    OutputFormatT outFormat;
    int argi;
    uint outRate, fftSize;
    int equalize, surface;
    double limit;
    uint truncSize;
    char * end = NULL;
  
    if (argc < 2) {
       fprintf (stderr, "Error:  No command specified.  See '%s -h' for help.\n", argv [0]);
       return (-1);
    }
    if ((strcmp (argv [1], "--help") == 0) || (strcmp (argv [1], "-h") == 0)) {
       fprintf (stdout, "HRTF Processing and Composition Utility\n\n");
       fprintf (stdout, "Usage:  %s <command> [<option>...]\n\n", argv [0]);
       fprintf (stdout, "Commands:\n");
       fprintf (stdout, " -m, --make-mhr  Makes an OpenAL Soft compatible HRTF data set.\n");
       fprintf (stdout, "                 Defaults output to: ./oalsoft_hrtf_%%r.mhr\n");
       fprintf (stdout, " -t, --make-tab  Makes the built-in table used when compiling OpenAL Soft.\n");
       fprintf (stdout, "                 Defaults output to: ./hrtf_tables.inc\n");
       fprintf (stdout, " -h, --help      Displays this help information.\n\n");
       fprintf (stdout, "Options:\n");
       fprintf (stdout, " -r=<rate>       Change the data set sample rate to the specified value and\n");
       fprintf (stdout, "                 resample the HRIRs accordingly.\n");
       fprintf (stdout, " -f=<points>     Override the FFT window size (defaults to the first power-\n");
       fprintf (stdout, "                 of-two that fits four times the number of HRIR points).\n");
       fprintf (stdout, " -e={on|off}     Toggle diffuse-field equalization (default: %s).\n", (DEFAULT_EQUALIZE ? "on" : "off"));
       fprintf (stdout, " -s={on|off}     Toggle surface-weighted diffuse-field average (default: %s).\n", (DEFAULT_SURFACE ? "on" : "off"));
       fprintf (stdout, " -l={<dB>|none}  Specify a limit to the magnitude range of the diffuse-field\n");
       fprintf (stdout, "                 average (default: %.2f).\n", DEFAULT_LIMIT);
       fprintf (stdout, " -w=<points>     Specify the size of the truncation window that's applied\n");
       fprintf (stdout, "                 after minimum-phase reconstruction (default: %u).\n", DEFAULT_TRUNCSIZE);
       fprintf (stdout, " -i=<filename>   Specify an HRIR definition file to use (defaults to stdin).\n");
       fprintf (stdout, " -o=<filename>   Specify an output file.  Overrides command-selected default.\n");
       fprintf (stdout, "                 Use of '%%r' will be substituted with the data set sample rate.\n");
       return (0);
    }
    if ((strcmp (argv [1], "--make-mhr") == 0) || (strcmp (argv [1], "-m") == 0)) {
       if (argc > 3)
          outName = argv [3];
       else
          outName = "./oalsoft_hrtf_%r.mhr";
       outFormat = OF_MHR;
    } else if ((strcmp (argv [1], "--make-tab") == 0) || (strcmp (argv [1], "-t") == 0)) {
       if (argc > 3)
          outName = argv [3];
       else
          outName = "./hrtf_tables.inc";
       outFormat = OF_TABLE;
    } else {
       fprintf (stderr, "Error:  Invalid command '%s'.\n", argv [1]);
       return (-1);
    }
    argi = 2;
    outRate = 0;
    fftSize = 0;
    equalize = DEFAULT_EQUALIZE;
    surface = DEFAULT_SURFACE;
    limit = DEFAULT_LIMIT;
    truncSize = DEFAULT_TRUNCSIZE;
    while (argi < argc) {
      if (strncmp (argv [argi], "-r=", 3) == 0) {
         outRate = strtoul (& argv [argi] [3], & end, 10);
         if ((end [0] != '\0') || (outRate < MIN_RATE) || (outRate > MAX_RATE)) {
            fprintf (stderr, "Error:  Expected a value from %u to %u for '-r'.\n", MIN_RATE, MAX_RATE);
            return (-1);
         }
      } else if (strncmp (argv [argi], "-f=", 3) == 0) {
         fftSize = strtoul (& argv [argi] [3], & end, 10);
         if ((end [0] != '\0') || (fftSize & (fftSize - 1)) || (fftSize < MIN_FFTSIZE) || (fftSize > MAX_FFTSIZE)) {
            fprintf (stderr, "Error:  Expected a power-of-two value from %u to %u for '-f'.\n", MIN_FFTSIZE, MAX_FFTSIZE);
            return (-1);
         }
      } else if (strncmp (argv [argi], "-e=", 3) == 0) {
         if (strcmp (& argv [argi] [3], "on") == 0) {
            equalize = 1;
         } else if (strcmp (& argv [argi] [3], "off") == 0) {
            equalize = 0;
         } else {
            fprintf (stderr, "Error:  Expected 'on' or 'off' for '-e'.\n");
            return (-1);
         }
      } else if (strncmp (argv [argi], "-s=", 3) == 0) {
         if (strcmp (& argv [argi] [3], "on") == 0) {
            surface = 1;
         } else if (strcmp (& argv [argi] [3], "off") == 0) {
            surface = 0;
         } else {
            fprintf (stderr, "Error:  Expected 'on' or 'off' for '-s'.\n");
            return (-1);
         }
      } else if (strncmp (argv [argi], "-l=", 3) == 0) {
         if (strcmp (& argv [argi] [3], "none") == 0) {
            limit = 0.0;
         } else {
            limit = strtod (& argv [argi] [3], & end);
            if ((end [0] != '\0') || (limit < MIN_LIMIT) || (limit > MAX_LIMIT)) {
               fprintf (stderr, "Error:  Expected 'none' or a value from %.2f to %.2f for '-l'.\n", MIN_LIMIT, MAX_LIMIT);
               return (-1);
            }
         }
      } else if (strncmp (argv [argi], "-w=", 3) == 0) {
         truncSize = strtoul (& argv [argi] [3], & end, 10);
         if ((end [0] != '\0') || (truncSize < MIN_TRUNCSIZE) || (truncSize > MAX_TRUNCSIZE) || (truncSize % MOD_TRUNCSIZE)) {
            fprintf (stderr, "Error:  Expected a value from %u to %u in multiples of %u for '-w'.\n", MIN_TRUNCSIZE, MAX_TRUNCSIZE, MOD_TRUNCSIZE);
            return (-1);
         }
      } else if (strncmp (argv [argi], "-i=", 3) == 0) {
         inName = & argv [argi] [3];
      } else if (strncmp (argv [argi], "-o=", 3) == 0) {
         outName = & argv [argi] [3];
      } else {
         fprintf (stderr, "Error:  Invalid option '%s'.\n", argv [argi]);
         return (-1);
      }
      argi ++;
    }
    if (! ProcessDefinition (inName, outRate, fftSize, equalize, surface, limit, truncSize, outFormat, outName))
       return (-1);
    fprintf (stdout, "Operation completed.\n");
    return (0);
  }
  
  
  
  
  