/* * HRTF utility for producing and demonstrating the process of creating an * OpenAL Soft compatible HRIR data set. * * Copyright (C) 2011-2014 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 */ #include "config.h" #include #include #include #include #include #include #ifdef HAVE_STRINGS_H #include #endif // Rely (if naively) on OpenAL's header for the types used for serialization. #include "AL/al.h" #include "AL/alext.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 limits to the custom head radius on the command line. #define MIN_CUSTOM_RADIUS (0.05) #define MAX_CUSTOM_RADIUS (0.15) // 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) #define DEFAULT_HEAD_MODEL (HM_DATASET) #define DEFAULT_CUSTOM_RADIUS (0.0) // 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. }; // Head model used for calculating the impulse delays. enum HeadModelT { HM_NONE = 0, HM_DATASET , // Measure the onset from the dataset. HM_SPHERE // Calculate the onset using a spherical head model. }; // Desired output format from the command line. enum OutputFormatT { OF_NONE = 0, OF_MHR // OpenAL Soft MHR data set file. }; // Unsigned integer type. typedef unsigned int uint; // Serialization types. The trailing digit indicates the number of bits. typedef ALubyte uint8; typedef ALint int32; typedef ALuint uint32; typedef ALuint64SOFT uint64; typedef enum ByteOrderT ByteOrderT; typedef enum SourceFormatT SourceFormatT; typedef enum ElementTypeT ElementTypeT; typedef enum HeadModelT HeadModelT; 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 versions < 1800 (Visual Studio 2013). #if defined(_MSC_VER) && _MSC_VER < 1800 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(size_t n) { double *a; if(n == 0) n = 1; a = 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(double *a) { free(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; if (outN == 0) return; // 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, uint32 * out) { uint8 in [4]; uint32 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, uint64 * out) { uint8 in [8]; uint64 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 uint32 in, FILE * fp, const char * filename) { uint8 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 { uint32 ui; int32 i; float f; } v4; union { uint64 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) { uint32 fourCC, chunkSize; uint32 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) { uint32 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) { uint32 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 onset time of an HRIR and average it with any existing // timing for its elevation and azimuth. static void AverageHrirOnset (const double * hrir, const double f, const uint ei, const uint ai, const HrirDataT * hData) { double mag; uint n, i, j; mag = 0.0; n = hData -> mIrPoints; for (i = 0; i < n; i ++) mag = fmax (fabs (hrir [i]), mag); mag *= 0.15; for (i = 0; i < n; i ++) { if (fabs (hrir [i]) >= mag) break; } j = hData -> mEvOffset [ei] + ai; hData -> mHrtds [j] = Lerp (hData -> mHrtds [j], ((double) i) / hData -> mIrRate, f); } // 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; } // Synthesize any missing onset timings at the bottom elevations. This just // blends between slightly exaggerated known onsets. Not an accurate model. static void SynthesizeOnsets (HrirDataT * hData) { uint oi, e, a, j0, j1; double t, of, jf; oi = hData -> mEvStart; t = 0.0; for (a = 0; a < hData -> mAzCount [oi]; a ++) t += hData -> mHrtds [hData -> mEvOffset [oi] + a]; hData -> mHrtds [0] = 1.32e-4 + (t / hData -> mAzCount [oi]); for (e = 1; e < hData -> mEvStart; e ++) { of = ((double) e) / hData -> mEvStart; for (a = 0; a < hData -> mAzCount [e]; a ++) { CalcAzIndices (hData, oi, a * 2.0 * M_PI / hData -> mAzCount [e], & j0, & j1, & jf); hData -> mHrtds [hData -> mEvOffset [e] + a] = Lerp (hData -> mHrtds [0], Lerp (hData -> mHrtds [j0], hData -> mHrtds [j1], jf), of); } } } /* 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 each minimum-phase // HRIR. static void CalculateHrtds (const HeadModelT model, const double radius, 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; if (model == HM_DATASET) { t = hData -> mHrtds [j] * radius / hData -> mRadius; } else { t = CalcLTD ((-90.0 + (e * 180.0 / (hData -> mEvCount - 1))) * M_PI / 180.0, (a * 360.0 / hData -> mAzCount [e]) * M_PI / 180.0, radius, 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, (uint32) hData -> mIrRate, fp, filename)) return (0); if (! WriteBin4 (BO_LITTLE, 1, (uint32) hData -> mIrPoints, fp, filename)) return (0); if (! WriteBin4 (BO_LITTLE, 1, (uint32) hData -> mEvCount, fp, filename)) return (0); for (e = 0; e < hData -> mEvCount; e ++) { if (! WriteBin4 (BO_LITTLE, 1, (uint32) 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, (uint32) 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, (uint32) v, 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 overridden 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 (const HeadModelT model, 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)) { if (model == HM_DATASET) AverageHrirOnset (hrir, 1.0 / factor, ei, ai, hData); 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 HeadModelT model, const double radius, 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, "", & 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 (model, & 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"); if (model == HM_DATASET) SynthesizeOnsets (& hData); SynthesizeHrirs (& hData); fprintf (stdout, "Normalizing final HRIRs...\n"); NormalizeHrirs (& hData); fprintf (stdout, "Calculating impulse delays...\n"); CalculateHrtds (model, (radius > DEFAULT_CUSTOM_RADIUS) ? radius : hData . mRadius, & 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; 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; HeadModelT model; double radius; 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 [