summaryrefslogtreecommitdiffstats
path: root/utils
diff options
context:
space:
mode:
authorChris Robinson <[email protected]>2012-03-25 19:57:40 -0700
committerChris Robinson <[email protected]>2012-03-25 19:57:40 -0700
commit929e9bfac1c30be4fb62eaa4a714c415b152f734 (patch)
treecf79a3d3e5f098f2c8510adb535f429807037512 /utils
parentb51afb74da12c790fde6804d99c982685c781d4e (diff)
Avoid DOS-style line-endings and stdint types for the makehrtf utility
Diffstat (limited to 'utils')
-rw-r--r--utils/makehrtf.c1479
1 files 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 <stdint.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <math.h>
-#include <string.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) {
- 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 <[email protected]> and Keith Martin <[email protected]> */\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 <command> <path of MIT set> [ <output file> ]\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 <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include <string.h>
+
+#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 <[email protected]> and Keith Martin <[email protected]> */\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 <command> <path of MIT set> [ <output file> ]\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;
+}