From 929e9bfac1c30be4fb62eaa4a714c415b152f734 Mon Sep 17 00:00:00 2001 From: Chris Robinson Date: Sun, 25 Mar 2012 19:57:40 -0700 Subject: Avoid DOS-style line-endings and stdint types for the makehrtf utility --- utils/makehrtf.c | 1479 +++++++++++++++++++++++++++--------------------------- 1 file changed, 740 insertions(+), 739 deletions(-) diff --git a/utils/makehrtf.c b/utils/makehrtf.c index a2638121..53ee46dd 100644 --- a/utils/makehrtf.c +++ b/utils/makehrtf.c @@ -1,739 +1,740 @@ -/** - * HRTF utility for producing and demonstrating the process of creating an - * OpenAL Soft compatible HRIR data set. - * - * It can currently make use of the 44.1 KHz diffuse and compact KEMAR HRIRs - * provided by MIT at: - * - * http://sound.media.mit.edu/resources/KEMAR.html - */ - -#include -#include -#include -#include -#include - -// The sample rate of the MIT HRIR data sets. -#define MIT_IR_RATE (44100) - -// The total number of used impulse responses from the MIT HRIR data sets. -#define MIT_IR_COUNT (828) - -// The size (in samples) of each HRIR in the MIT data sets. -#define MIT_IR_SIZE (128) - -// The total number of elevations given a step of 10 degrees. -#define MIT_EV_COUNT (19) - -// The first elevation that the MIT data sets have HRIRs for. -#define MIT_EV_START (5) - -// The head radius (in meters) used by the MIT data sets. -#define MIT_RADIUS (0.09f) - -// The source to listener distance (in meters) used by the MIT data sets. -#define MIT_DISTANCE (1.4f) - -// The resulting size (in samples) of a mininum-phase reconstructed HRIR. -#define MIN_IR_SIZE (32) - -// The size (in samples) of the real cepstrum used in reconstruction. This -// needs to be large enough to reduce inaccuracy. -#define CEP_SIZE (8192) - -// The OpenAL Soft HRTF format marker. It stands for minimum-phase head -// response protocol 00. -#define MHR_FORMAT ("MinPHR00") - -typedef struct ComplexT ComplexT; -typedef struct HrirDataT HrirDataT; - -// A complex number type. -struct ComplexT { - float mVec [2]; -}; - -// The HRIR data definition. This can be used to add support for new HRIR -// sources in the future. -struct HrirDataT { - int mIrRate, - mIrCount, - mIrSize, - mEvCount, - mEvStart; - const int * mEvOffset, - * mAzCount; - float mRadius, - mDistance, - * mHrirs, - * mHrtds, - mMaxHrtd; -}; - -// The linear index of the first HRIR for each elevation of the MIT data set. -static const int MIT_EV_OFFSET [MIT_EV_COUNT] = { - 0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827 -}, - -// The count of distinct azimuth steps for each elevation in the MIT data -// set. - MIT_AZ_COUNT [MIT_EV_COUNT] = { - 1, 12, 24, 36, 45, 56, 60, 72, 72, 72, 72, 72, 60, 56, 45, 36, 24, 12, 1 -}; - -// Performs a forward Fast Fourier Transform. -static void FftProc (int n, const ComplexT * fftIn, ComplexT * fftOut) { - int m2, rk, k, m; - float a, b; - int i; - float wx, wy; - int j, km2; - float tx, ty, wyd; - - // Data copy and bit-reversal ordering. - m2 = (n >> 1); - rk = 0; - for (k = 0; k < n; k ++) { - fftOut [rk] . mVec [0] = fftIn [k] . mVec [0]; - fftOut [rk] . mVec [1] = fftIn [k] . mVec [1]; - if (k < (n - 1)) { - m = m2; - while (rk >= m) { - rk -= m; - m >>= 1; - } - rk += m; - } - } - // Perform the FFT. - m2 = 1; - for (m = 2; m <= n; m <<= 1) { - a = sin (M_PI / m); - a = 2.0f * a * a; - b = sin (2.0f * M_PI / m); - for (i = 0; i < n; i += m) { - wx = 1.0f; - wy = 0.0f; - for (k = i, j = 0; j < m2; k ++, j ++) { - km2 = k + m2; - tx = (wx * fftOut [km2] . mVec [0]) - (wy * fftOut [km2] . mVec [1]); - ty = (wx * fftOut [km2] . mVec [1]) + (wy * fftOut [km2] . mVec [0]); - fftOut [km2] . mVec [0] = fftOut [k] . mVec [0] - tx; - fftOut [km2] . mVec [1] = fftOut [k] . mVec [1] - ty; - fftOut [k] . mVec [0] += tx; - fftOut [k] . mVec [1] += ty; - wyd = (a * wy) - (b * wx); - wx -= (a * wx) + (b * wy); - wy -= wyd; - } - } - m2 = m; - } -} - -// Performs an inverse Fast Fourier Transform. -static void FftInvProc (int n, const ComplexT * fftIn, ComplexT * fftOut) { - int m2, rk, k, m; - float a, b; - int i; - float wx, wy; - int j, km2; - float tx, ty, wyd, invn; - - // Data copy and bit-reversal ordering. - m2 = (n >> 1); - rk = 0; - for (k = 0; k < n; k ++) { - fftOut [rk] . mVec [0] = fftIn [k] . mVec [0]; - fftOut [rk] . mVec [1] = fftIn [k] . mVec [1]; - if (k < (n - 1)) { - m = m2; - while (rk >= m) { - rk -= m; - m >>= 1; - } - rk += m; - } - } - // Perform the IFFT. - m2 = 1; - for (m = 2; m <= n; m <<= 1) { - a = sin (M_PI / m); - a = 2.0f * a * a; - b = -sin (2.0f * M_PI / m); - for (i = 0; i < n; i += m) { - wx = 1.0f; - wy = 0.0f; - for (k = i, j = 0; j < m2; k ++, j ++) { - km2 = k + m2; - tx = (wx * fftOut [km2] . mVec [0]) - (wy * fftOut [km2] . mVec [1]); - ty = (wx * fftOut [km2] . mVec [1]) + (wy * fftOut [km2] . mVec [0]); - fftOut [km2] . mVec [0] = fftOut [k] . mVec [0] - tx; - fftOut [km2] . mVec [1] = fftOut [k] . mVec [1] - ty; - fftOut [k] . mVec [0] += tx; - fftOut [k] . mVec [1] += ty; - wyd = (a * wy) - (b * wx); - wx -= (a * wx) + (b * wy); - wy -= wyd; - } - } - m2 = m; - } - // Normalize the samples. - invn = 1.0f / n; - for (i = 0; i < n; i ++) { - fftOut [i] . mVec [0] *= invn; - fftOut [i] . mVec [1] *= invn; - } -} - -// Complex absolute value. -static void ComplexAbs (const ComplexT * in, ComplexT * out) { - out -> mVec [0] = sqrt ((in -> mVec [0] * in -> mVec [0]) + (in -> mVec [1] * in -> mVec [1])); - out -> mVec [1] = 0.0f; -} - -// Complex logarithm. -static void ComplexLog (const ComplexT * in, ComplexT * out) { - float r, t; - - r = sqrt ((in -> mVec [0] * in -> mVec [0]) + (in -> mVec [1] * in -> mVec [1])); - t = atan2 (in -> mVec [1], in -> mVec [0]); - if (t < 0.0f) - t += 2.0f * M_PI; - out -> mVec [0] = log (r); - out -> mVec [1] = t; -} - -// Complex exponent. -static void ComplexExp (const ComplexT * in, ComplexT * out) { - float e; - - e = exp (in -> mVec [0]); - out -> mVec [0] = e * cos (in -> mVec [1]); - out -> mVec [1] = e * sin (in -> mVec [1]); -} - -// Calculates the real cepstrum of a given impulse response. It currently -// uses a fixed cepstrum size. To make this more robust, it should be -// rewritten to handle a variable size cepstrum. -static void RealCepstrum (int irSize, const float * ir, float cep [CEP_SIZE]) { - ComplexT in [CEP_SIZE], out [CEP_SIZE]; - int index; - - for (index = 0; index < irSize; index ++) { - in [index] . mVec [0] = ir [index]; - in [index] . mVec [1] = 0.0f; - } - for (; index < CEP_SIZE; index ++) { - in [index] . mVec [0] = 0.0f; - in [index] . mVec [1] = 0.0f; - } - FftProc (CEP_SIZE, in, out); - for (index = 0; index < CEP_SIZE; index ++) { - ComplexAbs (& out [index], & out [index]); - if (out [index] . mVec [0] < 0.000001f) - out [index] . mVec [0] = 0.000001f; - ComplexLog (& out [index], & in [index]); - } - FftInvProc (CEP_SIZE, in, out); - for (index = 0; index < CEP_SIZE; index ++) - cep [index] = out [index] . mVec [0]; -} - -// Reconstructs the minimum-phase impulse response for a given real cepstrum. -// Like the above function, this should eventually be modified to handle a -// variable size cepstrum. -static void MinimumPhase (const float cep [CEP_SIZE], int irSize, float * mpIr) { - ComplexT in [CEP_SIZE], out [CEP_SIZE]; - int index; - - in [0] . mVec [0] = cep [0]; - for (index = 1; index < (CEP_SIZE / 2); index ++) - in [index] . mVec [0] = 2.0f * cep [index]; - if ((CEP_SIZE % 2) != 1) { - in [index] . mVec [0] = cep [index]; - index ++; - } - for (; index < CEP_SIZE; index ++) - in [index] . mVec [0] = 0.0f; - for (index = 0; index < CEP_SIZE; index ++) - in [index] . mVec [1] = 0.0f; - FftProc (CEP_SIZE, in, out); - for (index = 0; index < CEP_SIZE; index ++) - ComplexExp (& out [index], & in [index]); - FftInvProc (CEP_SIZE, in, out); - for (index = 0; index < irSize; index ++) - mpIr [index] = out [index] . mVec [0]; -} - -// Calculate the left-ear time delay using a spherical head model. -static float CalcLTD (float ev, float az, float rad, float dist) { - float azp, dlp, l, al; - - azp = asin (cos (ev) * sin (az)); - dlp = sqrt ((dist * dist) + (rad * rad) + (2.0f * dist * rad * sin (azp))); - l = sqrt ((dist * dist) - (rad * rad)); - al = (0.5f * M_PI) + azp; - if (dlp > l) - dlp = l + (rad * (al - acos (rad / dist))); - return (dlp / 343.3f); -} - -// Read a 16-bit little-endian integer from a file and convert it to a 32-bit -// floating-point value in the range of -1.0 to 1.0. -static int ReadInt16LeAsFloat32 (const char * fileName, FILE * fp, float * val) { - uint8_t vb [2]; - uint16_t vw; - - if (fread (vb, 1, sizeof (vb), fp) != sizeof (vb)) { - fclose (fp); - fprintf (stderr, "Error reading from file, '%s'.\n", fileName); - return (0); - } - vw = (((uint16_t) vb [1]) << 8) | vb [0]; - (* val) = ((int16_t) vw) / 32768.0f; - return (1); -} - -// Write a string to a file. -static int WriteString (const char * val, const char * fileName, FILE * fp) { - size_t len; - - len = strlen (val); - if (fwrite (val, 1, len, fp) != len) { - fclose (fp); - fprintf (stderr, "Error writing to file, '%s'.\n", fileName); - return (0); - } - return (1); -} - -// Write a 32-bit floating-point value in the range of -1.0 to 1.0 to a file -// as a 16-bit little-endian integer. -static int WriteFloat32AsInt16Le (float val, const char * fileName, FILE * fp) { - int16_t vw; - uint8_t vb [2]; - - vw = (short) round (32767.0f * val); - vb [0] = vw & 0x00FF; - vb [1] = (vw >> 8) & 0x00FF; - if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) { - fclose (fp); - fprintf (stderr, "Error writing to file, '%s'.\n", fileName); - return (0); - } - return (1); -} - -// Write a 32-bit little-endian unsigned integer to a file. -static int WriteUInt32Le (uint32_t val, const char * fileName, FILE * fp) { - uint8_t vb [4]; - - vb [0] = val & 0x000000FF; - vb [1] = (val >> 8) & 0x000000FF; - vb [2] = (val >> 16) & 0x000000FF; - vb [3] = (val >> 24) & 0x000000FF; - if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) { - fclose (fp); - fprintf (stderr, "Error writing to file, '%s'.\n", fileName); - return (0); - } - return (1); -} - -// Write a 16-bit little-endian unsigned integer to a file. -static int WriteUInt16Le (uint16_t val, const char * fileName, FILE * fp) { - uint8_t vb [2]; - - vb [0] = val & 0x00FF; - vb [1] = (val >> 8) & 0x00FF; - if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) { - fclose (fp); - fprintf (stderr, "Error writing to file, '%s'.\n", fileName); - return (0); - } - return (1); -} - -// Write an 8-bit unsigned integer to a file. -static int WriteUInt8 (uint8_t val, const char * fileName, FILE * fp) { - if (fwrite (& val, 1, sizeof (val), fp) != sizeof (val)) { - fclose (fp); - fprintf (stderr, "Error writing to file, '%s'.\n", fileName); - return (0); - } - return (1); -} - -// Load the MIT HRIRs. This loads the entire diffuse or compact set starting -// counter-clockwise up at the bottom elevation and clockwise at the forward -// azimuth. -static int LoadMitHrirs (const char * baseName, HrirDataT * hData) { - const int EV_ANGLE [MIT_EV_COUNT] = { - -90, -80, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90 - }; - int e, a; - char fileName [1024]; - FILE * fp = NULL; - int j0, j1, i; - float s; - - for (e = MIT_EV_START; e < MIT_EV_COUNT; e ++) { - for (a = 0; a < MIT_AZ_COUNT [e]; a ++) { - // The data packs the first 180 degrees in the left channel, and - // the last 180 degrees in the right channel. - if (round ((360.0f / MIT_AZ_COUNT [e]) * a) > 180.0f) - break; - // Determine which file to open. - snprintf (fileName, 1023, "%s%d/H%de%03da.wav", baseName, EV_ANGLE [e], EV_ANGLE [e], (int) round ((360.0f / MIT_AZ_COUNT [e]) * a)); - if ((fp = fopen (fileName, "rb")) == NULL) { - fprintf (stderr, "Could not open file, '%s'.\n", fileName); - return (0); - } - // Assuming they have not changed format, skip the .WAV header. - fseek (fp, 44, SEEK_SET); - // Map the left and right channels to their appropriate azimuth - // offsets. - j0 = (MIT_EV_OFFSET [e] + a) * MIT_IR_SIZE; - j1 = (MIT_EV_OFFSET [e] + ((MIT_AZ_COUNT [e] - a) % MIT_AZ_COUNT [e])) * MIT_IR_SIZE; - // Read in the data, converting it to floating-point. - for (i = 0; i < MIT_IR_SIZE; i ++) { - if (! ReadInt16LeAsFloat32 (fileName, fp, & s)) - return (0); - hData -> mHrirs [j0 + i] = s; - if (! ReadInt16LeAsFloat32 (fileName, fp, & s)) - return (0); - hData -> mHrirs [j1 + i] = s; - } - fclose (fp); - } - } - return (1); -} - -// Performs the minimum phase reconstruction for a given HRIR data set. The -// cepstrum size should be made configureable at some point in the future. -static void ReconstructHrirs (int minIrSize, HrirDataT * hData) { - int start, end, step, j; - float cep [CEP_SIZE]; - - start = hData -> mEvOffset [hData -> mEvStart]; - end = hData -> mIrCount; - step = hData -> mIrSize; - for (j = start; j < end; j ++) { - RealCepstrum (step, & hData -> mHrirs [j * step], cep); - MinimumPhase (cep, minIrSize, & hData -> mHrirs [j * minIrSize]); - } - hData -> mIrSize = minIrSize; -} - -// Renormalize the entire HRIR data set, and attenutate it slightly. -static void RenormalizeHrirs (const HrirDataT * hData) { - int step, start, end; - float norm; - int j, i; - - step = hData -> mIrSize; - start = hData -> mEvOffset [hData -> mEvStart] * step; - end = hData -> mIrCount * step; - norm = 0.0f; - for (j = start; j < end; j += step) { - for (i = 0; i < step; i ++) { - if (fabs (hData -> mHrirs [j + i]) > norm) - norm = fabs (hData -> mHrirs [j + i]); - } - } - if (norm > 0.000001f) - norm = 1.0f / norm; - norm *= 0.95f; - for (j = start; j < end; j += step) { - for (i = 0; i < step; i ++) - hData -> mHrirs [j + i] *= norm; - } -} - -// Given an elevation offset and azimuth, calculates two offsets for -// addressing the HRIRs buffer and their interpolation factor. -static void CalcAzIndices (const HrirDataT * hData, int oi, float az, int * j0, int * j1, float * jf) { - int ai; - - az = fmod ((2.0f * M_PI) + az, 2.0f * M_PI) * hData -> mAzCount [oi] / (2.0f * M_PI); - ai = (int) az; - az -= ai; - (* j0) = hData -> mEvOffset [oi] + ai; - (* j1) = hData -> mEvOffset [oi] + ((ai + 1) % hData -> mAzCount [oi]); - (* jf) = az; -} - -// Perform a linear interpolation. -static float Lerp (float a, float b, float f) { - return (a + (f * (b - a))); -} - -// Attempt to synthesize any missing HRIRs at the bottom elevations. Right -// now this just blends the lowest elevation HRIRs together and applies some -// attenuates and high frequency damping. It's not a realistic model to use, -// but it is simple. -static void SynthesizeHrirs (HrirDataT * hData) { - int step, oi, i, a, j, e; - float of; - int j0, j1; - float jf; - float lp [4], s0, s1; - - if (hData -> mEvStart <= 0) - return; - step = hData -> mIrSize; - oi = hData -> mEvStart; - for (i = 0; i < step; i ++) - hData -> mHrirs [i] = 0.0f; - for (a = 0; a < hData -> mAzCount [oi]; a ++) { - j = (hData -> mEvOffset [oi] + a) * step; - for (i = 0; i < step; i ++) - hData -> mHrirs [i] += hData -> mHrirs [j + i] / hData -> mAzCount [oi]; - } - for (e = 1; e < hData -> mEvStart; e ++) { - of = ((float) e) / hData -> mEvStart; - for (a = 0; a < hData -> mAzCount [e]; a ++) { - j = (hData -> mEvOffset [e] + a) * step; - CalcAzIndices (hData, oi, a * 2.0f * M_PI / hData -> mAzCount [e], & j0, & j1, & jf); - j0 *= step; - j1 *= step; - lp [0] = 0.0f; - lp [1] = 0.0f; - lp [2] = 0.0f; - lp [3] = 0.0f; - for (i = 0; i < step; 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], 0.15f - (0.15f * of)); - lp [1] = Lerp (lp [0], lp [1], 0.15f - (0.15f * of)); - lp [2] = Lerp (lp [1], lp [2], 0.15f - (0.15f * of)); - lp [3] = Lerp (lp [2], lp [3], 0.15f - (0.15f * of)); - hData -> mHrirs [j + i] = lp [3]; - } - } - } - lp [0] = 0.0f; - lp [1] = 0.0f; - lp [2] = 0.0f; - lp [3] = 0.0f; - for (i = 0; i < step; i ++) { - s0 = hData -> mHrirs [i]; - lp [0] = Lerp (s0, lp [0], 0.15f); - lp [1] = Lerp (lp [0], lp [1], 0.15f); - lp [2] = Lerp (lp [1], lp [2], 0.15f); - lp [3] = Lerp (lp [2], lp [3], 0.15f); - hData -> mHrirs [i] = lp [3]; - } - hData -> mEvStart = 0; -} - -// Calculate the effective head-related time delays for the each HRIR, now -// that they are minimum-phase. -static void CalculateHrtds (HrirDataT * hData) { - float minHrtd, maxHrtd; - int e, a, j; - float t; - - minHrtd = 1000.0f; - maxHrtd = -1000.0f; - for (e = 0; e < hData -> mEvCount; e ++) { - for (a = 0; a < hData -> mAzCount [e]; a ++) { - j = hData -> mEvOffset [e] + a; - t = CalcLTD ((-90.0f + (e * 180.0f / (hData -> mEvCount - 1))) * M_PI / 180.0f, - (a * 360.0f / hData -> mAzCount [e]) * M_PI / 180.0f, - hData -> mRadius, hData -> mDistance); - hData -> mHrtds [j] = t; - if (t > maxHrtd) - maxHrtd = t; - if (t < minHrtd) - minHrtd = t; - } - } - maxHrtd -= minHrtd; - for (j = 0; j < hData -> mIrCount; j ++) - hData -> mHrtds [j] -= minHrtd; - hData -> mMaxHrtd = maxHrtd; -} - -// Save the OpenAL Soft HRTF data set. -static int SaveMhr (const HrirDataT * hData, const char * fileName) { - FILE * fp = NULL; - int e, step, end, j, i; - - if ((fp = fopen (fileName, "wb")) == NULL) { - fprintf (stderr, "Could not create file, '%s'.\n", fileName); - return (0); - } - if (! WriteString (MHR_FORMAT, fileName, fp)) - return (0); - if (! WriteUInt32Le ((uint32_t) hData -> mIrRate, fileName, fp)) - return (0); - if (! WriteUInt16Le ((uint16_t) hData -> mIrCount, fileName, fp)) - return (0); - if (! WriteUInt16Le ((uint16_t) hData -> mIrSize, fileName, fp)) - return (0); - if (! WriteUInt8 ((uint8_t) hData -> mEvCount, fileName, fp)) - return (0); - for (e = 0; e < hData -> mEvCount; e ++) { - if (! WriteUInt16Le ((uint16_t) hData -> mEvOffset [e], fileName, fp)) - return (0); - } - step = hData -> mIrSize; - end = hData -> mIrCount * step; - for (j = 0; j < end; j += step) { - for (i = 0; i < step; i ++) { - if (! WriteFloat32AsInt16Le (hData -> mHrirs [j + i], fileName, fp)) - return (0); - } - } - for (j = 0; j < hData -> mIrCount; j ++) { - i = (int) round (44100.0f * hData -> mHrtds [j]); - if (i > 127) - i = 127; - if (! WriteUInt8 ((uint8_t) i, fileName, fp)) - return (0); - } - fclose (fp); - return (1); -} - -// Save the OpenAL Soft built-in table. -static int SaveTab (const HrirDataT * hData, const char * fileName) { - FILE * fp = NULL; - int step, end, j, i; - char text [16]; - - if ((fp = fopen (fileName, "wb")) == NULL) { - fprintf (stderr, "Could not create file, '%s'.\n", fileName); - return (0); - } - if (! WriteString ("/* This data is Copyright 1994 by the MIT Media Laboratory. It is provided free\n" - " * with no restrictions on use, provided the authors are cited when the data is\n" - " * used in any research or commercial application. */\n" - "/* Bill Gardner and Keith Martin */\n" - "\n" - " /* HRIR Coefficients */\n" - " {\n", fileName, fp)) - return (0); - step = hData -> mIrSize; - end = hData -> mIrCount * step; - for (j = 0; j < end; j += step) { - if (! WriteString (" { ", fileName, fp)) - return (0); - for (i = 0; i < step; i ++) { - snprintf (text, 15, "%+d, ", (int) round (32767.0f * hData -> mHrirs [j + i])); - if (! WriteString (text, fileName, fp)) - return (0); - } - if (! WriteString ("},\n", fileName, fp)) - return (0); - } - if (! WriteString (" },\n" - "\n" - " /* HRIR Delays */\n" - " { ", fileName, fp)) - return (0); - for (j = 0; j < hData -> mIrCount; j ++) { - snprintf (text, 15, "%d, ", (int) round (44100.0f * hData -> mHrtds [j])); - if (! WriteString (text, fileName, fp)) - return (0); - } - if (! WriteString ("}\n", fileName, fp)) - return (0); - fclose (fp); - return (1); -} - -// Loads and processes an MIT data set. At present, the HRIR and HRTD data -// is loaded and processed in a static buffer. That should change to using -// heap allocated memory in the future. A cleanup function will then be -// required. -static int MakeMit(const char *baseInName, HrirDataT *hData) -{ - static float hrirs[MIT_IR_COUNT * MIT_IR_SIZE]; - static float hrtds[MIT_IR_COUNT]; - - hData->mIrRate = MIT_IR_RATE; - hData->mIrCount = MIT_IR_COUNT; - hData->mIrSize = MIT_IR_SIZE; - hData->mEvCount = MIT_EV_COUNT; - hData->mEvStart = MIT_EV_START; - hData->mEvOffset = MIT_EV_OFFSET; - hData->mAzCount = MIT_AZ_COUNT; - hData->mRadius = MIT_RADIUS; - hData->mDistance = MIT_DISTANCE; - hData->mHrirs = hrirs; - hData->mHrtds = hrtds; - fprintf(stderr, "Loading base HRIR data...\n"); - if(!LoadMitHrirs(baseInName, hData)) - return 0; - fprintf(stderr, "Performing minimum phase reconstruction and truncation...\n"); - ReconstructHrirs(MIN_IR_SIZE, hData); - fprintf(stderr, "Renormalizing minimum phase HRIR data...\n"); - RenormalizeHrirs(hData); - fprintf(stderr, "Synthesizing missing elevations...\n"); - SynthesizeHrirs(hData); - fprintf(stderr, "Calculating impulse delays...\n"); - CalculateHrtds(hData); - return 1; -} - -// Simple dispatch. Provided a command, the path to the MIT set of choice, -// and an optional output filename, this will produce an OpenAL Soft -// compatible HRTF set in the chosen format. -int main(int argc, char *argv[]) -{ - char baseName[1024]; - const char *outName = NULL; - HrirDataT hData; - - if(argc < 3 || strcmp(argv [1], "-h") == 0 || strcmp (argv [1], "--help") == 0) - { - fprintf(stderr, "Usage: %s [ ]\n\n", argv[0]); - fprintf(stderr, "Commands:\n"); - fprintf(stderr, " -m, --make-mhr Makes an OpenAL Soft compatible HRTF data set.\n"); - fprintf(stderr, " Defaults output to: ./oal_soft_hrtf_44100.mhr\n"); - fprintf(stderr, " -t, --make-tab Makes the built-in table used when compiling OpenAL Soft.\n"); - fprintf(stderr, " Defaults output to: ./hrtf_tables.inc\n"); - fprintf(stderr, " -h, --help Displays this help information.\n"); - return 0; - } - - snprintf(baseName, sizeof(baseName), "%s/elev", argv[2]); - if(strcmp(argv[1], "-m") == 0 || strcmp(argv[1], "--make-mhr") == 0) - { - if(argc > 3) - outName = argv[3]; - else - outName = "./oal_soft_hrtf_44100.mhr"; - if(!MakeMit(baseName, &hData)) - return -1; - fprintf(stderr, "Creating data set file...\n"); - if(!SaveMhr(&hData, outName)) - return -1; - } - else if(strcmp(argv[1], "-t") == 0 || strcmp(argv[1], "--make-tab") == 0) - { - if(argc > 3) - outName = argv[3]; - else - outName = "./hrtf_tables.inc"; - if(!MakeMit(baseName, &hData)) - return -1; - fprintf(stderr, "Creating table file...\n"); - if(!SaveTab(&hData, outName)) - return -1; - } - else - { - fprintf(stderr, "Invalid command '%s'\n", argv[1]); - return -1; - } - fprintf(stderr, "Done.\n"); - return 0; -} +/** + * HRTF utility for producing and demonstrating the process of creating an + * OpenAL Soft compatible HRIR data set. + * + * It can currently make use of the 44.1 KHz diffuse and compact KEMAR HRIRs + * provided by MIT at: + * + * http://sound.media.mit.edu/resources/KEMAR.html + */ + +#include +#include +#include +#include + +#include "AL/al.h" + +// The sample rate of the MIT HRIR data sets. +#define MIT_IR_RATE (44100) + +// The total number of used impulse responses from the MIT HRIR data sets. +#define MIT_IR_COUNT (828) + +// The size (in samples) of each HRIR in the MIT data sets. +#define MIT_IR_SIZE (128) + +// The total number of elevations given a step of 10 degrees. +#define MIT_EV_COUNT (19) + +// The first elevation that the MIT data sets have HRIRs for. +#define MIT_EV_START (5) + +// The head radius (in meters) used by the MIT data sets. +#define MIT_RADIUS (0.09f) + +// The source to listener distance (in meters) used by the MIT data sets. +#define MIT_DISTANCE (1.4f) + +// The resulting size (in samples) of a mininum-phase reconstructed HRIR. +#define MIN_IR_SIZE (32) + +// The size (in samples) of the real cepstrum used in reconstruction. This +// needs to be large enough to reduce inaccuracy. +#define CEP_SIZE (8192) + +// The OpenAL Soft HRTF format marker. It stands for minimum-phase head +// response protocol 00. +#define MHR_FORMAT ("MinPHR00") + +typedef struct ComplexT ComplexT; +typedef struct HrirDataT HrirDataT; + +// A complex number type. +struct ComplexT { + float mVec [2]; +}; + +// The HRIR data definition. This can be used to add support for new HRIR +// sources in the future. +struct HrirDataT { + int mIrRate, + mIrCount, + mIrSize, + mEvCount, + mEvStart; + const int * mEvOffset, + * mAzCount; + float mRadius, + mDistance, + * mHrirs, + * mHrtds, + mMaxHrtd; +}; + +// The linear index of the first HRIR for each elevation of the MIT data set. +static const int MIT_EV_OFFSET [MIT_EV_COUNT] = { + 0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827 +}, + +// The count of distinct azimuth steps for each elevation in the MIT data +// set. + MIT_AZ_COUNT [MIT_EV_COUNT] = { + 1, 12, 24, 36, 45, 56, 60, 72, 72, 72, 72, 72, 60, 56, 45, 36, 24, 12, 1 +}; + +// Performs a forward Fast Fourier Transform. +static void FftProc (int n, const ComplexT * fftIn, ComplexT * fftOut) { + int m2, rk, k, m; + float a, b; + int i; + float wx, wy; + int j, km2; + float tx, ty, wyd; + + // Data copy and bit-reversal ordering. + m2 = (n >> 1); + rk = 0; + for (k = 0; k < n; k ++) { + fftOut [rk] . mVec [0] = fftIn [k] . mVec [0]; + fftOut [rk] . mVec [1] = fftIn [k] . mVec [1]; + if (k < (n - 1)) { + m = m2; + while (rk >= m) { + rk -= m; + m >>= 1; + } + rk += m; + } + } + // Perform the FFT. + m2 = 1; + for (m = 2; m <= n; m <<= 1) { + a = sin (M_PI / m); + a = 2.0f * a * a; + b = sin (2.0f * M_PI / m); + for (i = 0; i < n; i += m) { + wx = 1.0f; + wy = 0.0f; + for (k = i, j = 0; j < m2; k ++, j ++) { + km2 = k + m2; + tx = (wx * fftOut [km2] . mVec [0]) - (wy * fftOut [km2] . mVec [1]); + ty = (wx * fftOut [km2] . mVec [1]) + (wy * fftOut [km2] . mVec [0]); + fftOut [km2] . mVec [0] = fftOut [k] . mVec [0] - tx; + fftOut [km2] . mVec [1] = fftOut [k] . mVec [1] - ty; + fftOut [k] . mVec [0] += tx; + fftOut [k] . mVec [1] += ty; + wyd = (a * wy) - (b * wx); + wx -= (a * wx) + (b * wy); + wy -= wyd; + } + } + m2 = m; + } +} + +// Performs an inverse Fast Fourier Transform. +static void FftInvProc (int n, const ComplexT * fftIn, ComplexT * fftOut) { + int m2, rk, k, m; + float a, b; + int i; + float wx, wy; + int j, km2; + float tx, ty, wyd, invn; + + // Data copy and bit-reversal ordering. + m2 = (n >> 1); + rk = 0; + for (k = 0; k < n; k ++) { + fftOut [rk] . mVec [0] = fftIn [k] . mVec [0]; + fftOut [rk] . mVec [1] = fftIn [k] . mVec [1]; + if (k < (n - 1)) { + m = m2; + while (rk >= m) { + rk -= m; + m >>= 1; + } + rk += m; + } + } + // Perform the IFFT. + m2 = 1; + for (m = 2; m <= n; m <<= 1) { + a = sin (M_PI / m); + a = 2.0f * a * a; + b = -sin (2.0f * M_PI / m); + for (i = 0; i < n; i += m) { + wx = 1.0f; + wy = 0.0f; + for (k = i, j = 0; j < m2; k ++, j ++) { + km2 = k + m2; + tx = (wx * fftOut [km2] . mVec [0]) - (wy * fftOut [km2] . mVec [1]); + ty = (wx * fftOut [km2] . mVec [1]) + (wy * fftOut [km2] . mVec [0]); + fftOut [km2] . mVec [0] = fftOut [k] . mVec [0] - tx; + fftOut [km2] . mVec [1] = fftOut [k] . mVec [1] - ty; + fftOut [k] . mVec [0] += tx; + fftOut [k] . mVec [1] += ty; + wyd = (a * wy) - (b * wx); + wx -= (a * wx) + (b * wy); + wy -= wyd; + } + } + m2 = m; + } + // Normalize the samples. + invn = 1.0f / n; + for (i = 0; i < n; i ++) { + fftOut [i] . mVec [0] *= invn; + fftOut [i] . mVec [1] *= invn; + } +} + +// Complex absolute value. +static void ComplexAbs (const ComplexT * in, ComplexT * out) { + out -> mVec [0] = sqrt ((in -> mVec [0] * in -> mVec [0]) + (in -> mVec [1] * in -> mVec [1])); + out -> mVec [1] = 0.0f; +} + +// Complex logarithm. +static void ComplexLog (const ComplexT * in, ComplexT * out) { + float r, t; + + r = sqrt ((in -> mVec [0] * in -> mVec [0]) + (in -> mVec [1] * in -> mVec [1])); + t = atan2 (in -> mVec [1], in -> mVec [0]); + if (t < 0.0f) + t += 2.0f * M_PI; + out -> mVec [0] = log (r); + out -> mVec [1] = t; +} + +// Complex exponent. +static void ComplexExp (const ComplexT * in, ComplexT * out) { + float e; + + e = exp (in -> mVec [0]); + out -> mVec [0] = e * cos (in -> mVec [1]); + out -> mVec [1] = e * sin (in -> mVec [1]); +} + +// Calculates the real cepstrum of a given impulse response. It currently +// uses a fixed cepstrum size. To make this more robust, it should be +// rewritten to handle a variable size cepstrum. +static void RealCepstrum (int irSize, const float * ir, float cep [CEP_SIZE]) { + ComplexT in [CEP_SIZE], out [CEP_SIZE]; + int index; + + for (index = 0; index < irSize; index ++) { + in [index] . mVec [0] = ir [index]; + in [index] . mVec [1] = 0.0f; + } + for (; index < CEP_SIZE; index ++) { + in [index] . mVec [0] = 0.0f; + in [index] . mVec [1] = 0.0f; + } + FftProc (CEP_SIZE, in, out); + for (index = 0; index < CEP_SIZE; index ++) { + ComplexAbs (& out [index], & out [index]); + if (out [index] . mVec [0] < 0.000001f) + out [index] . mVec [0] = 0.000001f; + ComplexLog (& out [index], & in [index]); + } + FftInvProc (CEP_SIZE, in, out); + for (index = 0; index < CEP_SIZE; index ++) + cep [index] = out [index] . mVec [0]; +} + +// Reconstructs the minimum-phase impulse response for a given real cepstrum. +// Like the above function, this should eventually be modified to handle a +// variable size cepstrum. +static void MinimumPhase (const float cep [CEP_SIZE], int irSize, float * mpIr) { + ComplexT in [CEP_SIZE], out [CEP_SIZE]; + int index; + + in [0] . mVec [0] = cep [0]; + for (index = 1; index < (CEP_SIZE / 2); index ++) + in [index] . mVec [0] = 2.0f * cep [index]; + if ((CEP_SIZE % 2) != 1) { + in [index] . mVec [0] = cep [index]; + index ++; + } + for (; index < CEP_SIZE; index ++) + in [index] . mVec [0] = 0.0f; + for (index = 0; index < CEP_SIZE; index ++) + in [index] . mVec [1] = 0.0f; + FftProc (CEP_SIZE, in, out); + for (index = 0; index < CEP_SIZE; index ++) + ComplexExp (& out [index], & in [index]); + FftInvProc (CEP_SIZE, in, out); + for (index = 0; index < irSize; index ++) + mpIr [index] = out [index] . mVec [0]; +} + +// Calculate the left-ear time delay using a spherical head model. +static float CalcLTD (float ev, float az, float rad, float dist) { + float azp, dlp, l, al; + + azp = asin (cos (ev) * sin (az)); + dlp = sqrt ((dist * dist) + (rad * rad) + (2.0f * dist * rad * sin (azp))); + l = sqrt ((dist * dist) - (rad * rad)); + al = (0.5f * M_PI) + azp; + if (dlp > l) + dlp = l + (rad * (al - acos (rad / dist))); + return (dlp / 343.3f); +} + +// Read a 16-bit little-endian integer from a file and convert it to a 32-bit +// floating-point value in the range of -1.0 to 1.0. +static int ReadInt16LeAsFloat32 (const char * fileName, FILE * fp, float * val) { + ALubyte vb [2]; + ALushort vw; + + if (fread (vb, 1, sizeof (vb), fp) != sizeof (vb)) { + fclose (fp); + fprintf (stderr, "Error reading from file, '%s'.\n", fileName); + return (0); + } + vw = (((unsigned short) vb [1]) << 8) | vb [0]; + (* val) = ((short) vw) / 32768.0f; + return (1); +} + +// Write a string to a file. +static int WriteString (const char * val, const char * fileName, FILE * fp) { + size_t len; + + len = strlen (val); + if (fwrite (val, 1, len, fp) != len) { + fclose (fp); + fprintf (stderr, "Error writing to file, '%s'.\n", fileName); + return (0); + } + return (1); +} + +// Write a 32-bit floating-point value in the range of -1.0 to 1.0 to a file +// as a 16-bit little-endian integer. +static int WriteFloat32AsInt16Le (float val, const char * fileName, FILE * fp) { + ALshort vw; + ALubyte vb [2]; + + vw = (short) round (32767.0f * val); + vb [0] = vw & 0x00FF; + vb [1] = (vw >> 8) & 0x00FF; + if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) { + fclose (fp); + fprintf (stderr, "Error writing to file, '%s'.\n", fileName); + return (0); + } + return (1); +} + +// Write a 32-bit little-endian unsigned integer to a file. +static int WriteUInt32Le (ALuint val, const char * fileName, FILE * fp) { + ALubyte vb [4]; + + vb [0] = val & 0x000000FF; + vb [1] = (val >> 8) & 0x000000FF; + vb [2] = (val >> 16) & 0x000000FF; + vb [3] = (val >> 24) & 0x000000FF; + if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) { + fclose (fp); + fprintf (stderr, "Error writing to file, '%s'.\n", fileName); + return (0); + } + return (1); +} + +// Write a 16-bit little-endian unsigned integer to a file. +static int WriteUInt16Le (ALushort val, const char * fileName, FILE * fp) { + ALubyte vb [2]; + + vb [0] = val & 0x00FF; + vb [1] = (val >> 8) & 0x00FF; + if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) { + fclose (fp); + fprintf (stderr, "Error writing to file, '%s'.\n", fileName); + return (0); + } + return (1); +} + +// Write an 8-bit unsigned integer to a file. +static int WriteUInt8 (ALubyte val, const char * fileName, FILE * fp) { + if (fwrite (& val, 1, sizeof (val), fp) != sizeof (val)) { + fclose (fp); + fprintf (stderr, "Error writing to file, '%s'.\n", fileName); + return (0); + } + return (1); +} + +// Load the MIT HRIRs. This loads the entire diffuse or compact set starting +// counter-clockwise up at the bottom elevation and clockwise at the forward +// azimuth. +static int LoadMitHrirs (const char * baseName, HrirDataT * hData) { + const int EV_ANGLE [MIT_EV_COUNT] = { + -90, -80, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90 + }; + int e, a; + char fileName [1024]; + FILE * fp = NULL; + int j0, j1, i; + float s; + + for (e = MIT_EV_START; e < MIT_EV_COUNT; e ++) { + for (a = 0; a < MIT_AZ_COUNT [e]; a ++) { + // The data packs the first 180 degrees in the left channel, and + // the last 180 degrees in the right channel. + if (round ((360.0f / MIT_AZ_COUNT [e]) * a) > 180.0f) + break; + // Determine which file to open. + snprintf (fileName, 1023, "%s%d/H%de%03da.wav", baseName, EV_ANGLE [e], EV_ANGLE [e], (int) round ((360.0f / MIT_AZ_COUNT [e]) * a)); + if ((fp = fopen (fileName, "rb")) == NULL) { + fprintf (stderr, "Could not open file, '%s'.\n", fileName); + return (0); + } + // Assuming they have not changed format, skip the .WAV header. + fseek (fp, 44, SEEK_SET); + // Map the left and right channels to their appropriate azimuth + // offsets. + j0 = (MIT_EV_OFFSET [e] + a) * MIT_IR_SIZE; + j1 = (MIT_EV_OFFSET [e] + ((MIT_AZ_COUNT [e] - a) % MIT_AZ_COUNT [e])) * MIT_IR_SIZE; + // Read in the data, converting it to floating-point. + for (i = 0; i < MIT_IR_SIZE; i ++) { + if (! ReadInt16LeAsFloat32 (fileName, fp, & s)) + return (0); + hData -> mHrirs [j0 + i] = s; + if (! ReadInt16LeAsFloat32 (fileName, fp, & s)) + return (0); + hData -> mHrirs [j1 + i] = s; + } + fclose (fp); + } + } + return (1); +} + +// Performs the minimum phase reconstruction for a given HRIR data set. The +// cepstrum size should be made configureable at some point in the future. +static void ReconstructHrirs (int minIrSize, HrirDataT * hData) { + int start, end, step, j; + float cep [CEP_SIZE]; + + start = hData -> mEvOffset [hData -> mEvStart]; + end = hData -> mIrCount; + step = hData -> mIrSize; + for (j = start; j < end; j ++) { + RealCepstrum (step, & hData -> mHrirs [j * step], cep); + MinimumPhase (cep, minIrSize, & hData -> mHrirs [j * minIrSize]); + } + hData -> mIrSize = minIrSize; +} + +// Renormalize the entire HRIR data set, and attenutate it slightly. +static void RenormalizeHrirs (const HrirDataT * hData) { + int step, start, end; + float norm; + int j, i; + + step = hData -> mIrSize; + start = hData -> mEvOffset [hData -> mEvStart] * step; + end = hData -> mIrCount * step; + norm = 0.0f; + for (j = start; j < end; j += step) { + for (i = 0; i < step; i ++) { + if (fabs (hData -> mHrirs [j + i]) > norm) + norm = fabs (hData -> mHrirs [j + i]); + } + } + if (norm > 0.000001f) + norm = 1.0f / norm; + norm *= 0.95f; + for (j = start; j < end; j += step) { + for (i = 0; i < step; i ++) + hData -> mHrirs [j + i] *= norm; + } +} + +// Given an elevation offset and azimuth, calculates two offsets for +// addressing the HRIRs buffer and their interpolation factor. +static void CalcAzIndices (const HrirDataT * hData, int oi, float az, int * j0, int * j1, float * jf) { + int ai; + + az = fmod ((2.0f * M_PI) + az, 2.0f * M_PI) * hData -> mAzCount [oi] / (2.0f * M_PI); + ai = (int) az; + az -= ai; + (* j0) = hData -> mEvOffset [oi] + ai; + (* j1) = hData -> mEvOffset [oi] + ((ai + 1) % hData -> mAzCount [oi]); + (* jf) = az; +} + +// Perform a linear interpolation. +static float Lerp (float a, float b, float f) { + return (a + (f * (b - a))); +} + +// Attempt to synthesize any missing HRIRs at the bottom elevations. Right +// now this just blends the lowest elevation HRIRs together and applies some +// attenuates and high frequency damping. It's not a realistic model to use, +// but it is simple. +static void SynthesizeHrirs (HrirDataT * hData) { + int step, oi, i, a, j, e; + float of; + int j0, j1; + float jf; + float lp [4], s0, s1; + + if (hData -> mEvStart <= 0) + return; + step = hData -> mIrSize; + oi = hData -> mEvStart; + for (i = 0; i < step; i ++) + hData -> mHrirs [i] = 0.0f; + for (a = 0; a < hData -> mAzCount [oi]; a ++) { + j = (hData -> mEvOffset [oi] + a) * step; + for (i = 0; i < step; i ++) + hData -> mHrirs [i] += hData -> mHrirs [j + i] / hData -> mAzCount [oi]; + } + for (e = 1; e < hData -> mEvStart; e ++) { + of = ((float) e) / hData -> mEvStart; + for (a = 0; a < hData -> mAzCount [e]; a ++) { + j = (hData -> mEvOffset [e] + a) * step; + CalcAzIndices (hData, oi, a * 2.0f * M_PI / hData -> mAzCount [e], & j0, & j1, & jf); + j0 *= step; + j1 *= step; + lp [0] = 0.0f; + lp [1] = 0.0f; + lp [2] = 0.0f; + lp [3] = 0.0f; + for (i = 0; i < step; 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], 0.15f - (0.15f * of)); + lp [1] = Lerp (lp [0], lp [1], 0.15f - (0.15f * of)); + lp [2] = Lerp (lp [1], lp [2], 0.15f - (0.15f * of)); + lp [3] = Lerp (lp [2], lp [3], 0.15f - (0.15f * of)); + hData -> mHrirs [j + i] = lp [3]; + } + } + } + lp [0] = 0.0f; + lp [1] = 0.0f; + lp [2] = 0.0f; + lp [3] = 0.0f; + for (i = 0; i < step; i ++) { + s0 = hData -> mHrirs [i]; + lp [0] = Lerp (s0, lp [0], 0.15f); + lp [1] = Lerp (lp [0], lp [1], 0.15f); + lp [2] = Lerp (lp [1], lp [2], 0.15f); + lp [3] = Lerp (lp [2], lp [3], 0.15f); + hData -> mHrirs [i] = lp [3]; + } + hData -> mEvStart = 0; +} + +// Calculate the effective head-related time delays for the each HRIR, now +// that they are minimum-phase. +static void CalculateHrtds (HrirDataT * hData) { + float minHrtd, maxHrtd; + int e, a, j; + float t; + + minHrtd = 1000.0f; + maxHrtd = -1000.0f; + for (e = 0; e < hData -> mEvCount; e ++) { + for (a = 0; a < hData -> mAzCount [e]; a ++) { + j = hData -> mEvOffset [e] + a; + t = CalcLTD ((-90.0f + (e * 180.0f / (hData -> mEvCount - 1))) * M_PI / 180.0f, + (a * 360.0f / hData -> mAzCount [e]) * M_PI / 180.0f, + hData -> mRadius, hData -> mDistance); + hData -> mHrtds [j] = t; + if (t > maxHrtd) + maxHrtd = t; + if (t < minHrtd) + minHrtd = t; + } + } + maxHrtd -= minHrtd; + for (j = 0; j < hData -> mIrCount; j ++) + hData -> mHrtds [j] -= minHrtd; + hData -> mMaxHrtd = maxHrtd; +} + +// Save the OpenAL Soft HRTF data set. +static int SaveMhr (const HrirDataT * hData, const char * fileName) { + FILE * fp = NULL; + int e, step, end, j, i; + + if ((fp = fopen (fileName, "wb")) == NULL) { + fprintf (stderr, "Could not create file, '%s'.\n", fileName); + return (0); + } + if (! WriteString (MHR_FORMAT, fileName, fp)) + return (0); + if (! WriteUInt32Le ((ALuint) hData -> mIrRate, fileName, fp)) + return (0); + if (! WriteUInt16Le ((ALushort) hData -> mIrCount, fileName, fp)) + return (0); + if (! WriteUInt16Le ((ALushort) hData -> mIrSize, fileName, fp)) + return (0); + if (! WriteUInt8 ((ALubyte) hData -> mEvCount, fileName, fp)) + return (0); + for (e = 0; e < hData -> mEvCount; e ++) { + if (! WriteUInt16Le ((ALushort) hData -> mEvOffset [e], fileName, fp)) + return (0); + } + step = hData -> mIrSize; + end = hData -> mIrCount * step; + for (j = 0; j < end; j += step) { + for (i = 0; i < step; i ++) { + if (! WriteFloat32AsInt16Le (hData -> mHrirs [j + i], fileName, fp)) + return (0); + } + } + for (j = 0; j < hData -> mIrCount; j ++) { + i = (int) round (44100.0f * hData -> mHrtds [j]); + if (i > 127) + i = 127; + if (! WriteUInt8 ((ALubyte) i, fileName, fp)) + return (0); + } + fclose (fp); + return (1); +} + +// Save the OpenAL Soft built-in table. +static int SaveTab (const HrirDataT * hData, const char * fileName) { + FILE * fp = NULL; + int step, end, j, i; + char text [16]; + + if ((fp = fopen (fileName, "wb")) == NULL) { + fprintf (stderr, "Could not create file, '%s'.\n", fileName); + return (0); + } + if (! WriteString ("/* This data is Copyright 1994 by the MIT Media Laboratory. It is provided free\n" + " * with no restrictions on use, provided the authors are cited when the data is\n" + " * used in any research or commercial application. */\n" + "/* Bill Gardner and Keith Martin */\n" + "\n" + " /* HRIR Coefficients */\n" + " {\n", fileName, fp)) + return (0); + step = hData -> mIrSize; + end = hData -> mIrCount * step; + for (j = 0; j < end; j += step) { + if (! WriteString (" { ", fileName, fp)) + return (0); + for (i = 0; i < step; i ++) { + snprintf (text, 15, "%+d, ", (int) round (32767.0f * hData -> mHrirs [j + i])); + if (! WriteString (text, fileName, fp)) + return (0); + } + if (! WriteString ("},\n", fileName, fp)) + return (0); + } + if (! WriteString (" },\n" + "\n" + " /* HRIR Delays */\n" + " { ", fileName, fp)) + return (0); + for (j = 0; j < hData -> mIrCount; j ++) { + snprintf (text, 15, "%d, ", (int) round (44100.0f * hData -> mHrtds [j])); + if (! WriteString (text, fileName, fp)) + return (0); + } + if (! WriteString ("}\n", fileName, fp)) + return (0); + fclose (fp); + return (1); +} + +// Loads and processes an MIT data set. At present, the HRIR and HRTD data +// is loaded and processed in a static buffer. That should change to using +// heap allocated memory in the future. A cleanup function will then be +// required. +static int MakeMit(const char *baseInName, HrirDataT *hData) +{ + static float hrirs[MIT_IR_COUNT * MIT_IR_SIZE]; + static float hrtds[MIT_IR_COUNT]; + + hData->mIrRate = MIT_IR_RATE; + hData->mIrCount = MIT_IR_COUNT; + hData->mIrSize = MIT_IR_SIZE; + hData->mEvCount = MIT_EV_COUNT; + hData->mEvStart = MIT_EV_START; + hData->mEvOffset = MIT_EV_OFFSET; + hData->mAzCount = MIT_AZ_COUNT; + hData->mRadius = MIT_RADIUS; + hData->mDistance = MIT_DISTANCE; + hData->mHrirs = hrirs; + hData->mHrtds = hrtds; + fprintf(stderr, "Loading base HRIR data...\n"); + if(!LoadMitHrirs(baseInName, hData)) + return 0; + fprintf(stderr, "Performing minimum phase reconstruction and truncation...\n"); + ReconstructHrirs(MIN_IR_SIZE, hData); + fprintf(stderr, "Renormalizing minimum phase HRIR data...\n"); + RenormalizeHrirs(hData); + fprintf(stderr, "Synthesizing missing elevations...\n"); + SynthesizeHrirs(hData); + fprintf(stderr, "Calculating impulse delays...\n"); + CalculateHrtds(hData); + return 1; +} + +// Simple dispatch. Provided a command, the path to the MIT set of choice, +// and an optional output filename, this will produce an OpenAL Soft +// compatible HRTF set in the chosen format. +int main(int argc, char *argv[]) +{ + char baseName[1024]; + const char *outName = NULL; + HrirDataT hData; + + if(argc < 3 || strcmp(argv [1], "-h") == 0 || strcmp (argv [1], "--help") == 0) + { + fprintf(stderr, "Usage: %s [ ]\n\n", argv[0]); + fprintf(stderr, "Commands:\n"); + fprintf(stderr, " -m, --make-mhr Makes an OpenAL Soft compatible HRTF data set.\n"); + fprintf(stderr, " Defaults output to: ./oal_soft_hrtf_44100.mhr\n"); + fprintf(stderr, " -t, --make-tab Makes the built-in table used when compiling OpenAL Soft.\n"); + fprintf(stderr, " Defaults output to: ./hrtf_tables.inc\n"); + fprintf(stderr, " -h, --help Displays this help information.\n"); + return 0; + } + + snprintf(baseName, sizeof(baseName), "%s/elev", argv[2]); + if(strcmp(argv[1], "-m") == 0 || strcmp(argv[1], "--make-mhr") == 0) + { + if(argc > 3) + outName = argv[3]; + else + outName = "./oal_soft_hrtf_44100.mhr"; + if(!MakeMit(baseName, &hData)) + return -1; + fprintf(stderr, "Creating data set file...\n"); + if(!SaveMhr(&hData, outName)) + return -1; + } + else if(strcmp(argv[1], "-t") == 0 || strcmp(argv[1], "--make-tab") == 0) + { + if(argc > 3) + outName = argv[3]; + else + outName = "./hrtf_tables.inc"; + if(!MakeMit(baseName, &hData)) + return -1; + fprintf(stderr, "Creating table file...\n"); + if(!SaveTab(&hData, outName)) + return -1; + } + else + { + fprintf(stderr, "Invalid command '%s'\n", argv[1]); + return -1; + } + fprintf(stderr, "Done.\n"); + return 0; +} -- cgit v1.2.3