utils/makehrtf.c

   1 /**
   2  * HRTF utility for producing and demonstrating the process of creating an
   3  * OpenAL Soft compatible HRIR data set.
   4  *
   5  * It can currently make use of the 44.1 KHz diffuse and compact KEMAR HRIRs
   6  * provided by MIT at:
   7  *
   8  *   http://sound.media.mit.edu/resources/KEMAR.html
   9  */
  10
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <math.h>
  14 #include <string.h>
  15
  16 #include "AL/al.h"
  17
  18 #ifndef M_PI
  19 #define M_PI 3.14159265358979323846
  20 #endif
  21
  22 #ifdef _MSC_VER
  23 static double round(double val)
  24 {
  25     if(val < 0.0)
  26         return ceil(val - 0.5);
  27     return floor(val + 0.5);
  28 }
  29 #endif
  30
  31 // The sample rate of the MIT HRIR data sets.
  32 #define MIT_IR_RATE                  (44100)
  33
  34 // The total number of used impulse responses from the MIT HRIR data sets.
  35 #define MIT_IR_COUNT                 (828)
  36
  37 // The size (in samples) of each HRIR in the MIT data sets.
  38 #define MIT_IR_SIZE                  (128)
  39
  40 // The total number of elevations given a step of 10 degrees.
  41 #define MIT_EV_COUNT                 (19)
  42
  43 // The first elevation that the MIT data sets have HRIRs for.
  44 #define MIT_EV_START                 (5)
  45
  46 // The head radius (in meters) used by the MIT data sets.
  47 #define MIT_RADIUS                   (0.09f)
  48
  49 // The source to listener distance (in meters) used by the MIT data sets.
  50 #define MIT_DISTANCE                 (1.4f)
  51
  52 // The resulting size (in samples) of a mininum-phase reconstructed HRIR.
  53 #define MIN_IR_SIZE                  (32)
  54
  55 // The size (in samples) of the real cepstrum used in reconstruction.  This
  56 // needs to be large enough to reduce inaccuracy.
  57 #define CEP_SIZE                     (8192)
  58
  59 // The OpenAL Soft HRTF format marker.  It stands for minimum-phase head
  60 // response protocol 00.
  61 #define MHR_FORMAT                   ("MinPHR00")
  62
  63 typedef struct ComplexT              ComplexT;
  64 typedef struct HrirDataT             HrirDataT;
  65
  66 // A complex number type.
  67 struct ComplexT {
  68   float                              mVec [2];
  69 };
  70
  71 // The HRIR data definition.  This can be used to add support for new HRIR
  72 // sources in the future.
  73 struct HrirDataT {
  74   int                                mIrRate,
  75                                      mIrCount,
  76                                      mIrSize,
  77                                      mEvCount,
  78                                      mEvStart;
  79   const int                        * mEvOffset,
  80                                    * mAzCount;
  81   float                              mRadius,
  82                                      mDistance,
  83                                    * mHrirs,
  84                                    * mHrtds,
  85                                      mMaxHrtd;
  86 };
  87
  88 // The linear index of the first HRIR for each elevation of the MIT data set.
  89 static const int                     MIT_EV_OFFSET [MIT_EV_COUNT] = {
  90   0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827
  91 },
  92
  93 // The count of distinct azimuth steps for each elevation in the MIT data
  94 // set.
  95                                      MIT_AZ_COUNT [MIT_EV_COUNT] = {
  96   1, 12, 24, 36, 45, 56, 60, 72, 72, 72, 72, 72, 60, 56, 45, 36, 24, 12, 1
  97 };
  98
  99 // Performs a forward Fast Fourier Transform.
 100 static void FftProc (int n, const ComplexT * fftIn, ComplexT * fftOut) {
 101   int m2, rk, k, m;
 102   float a, b;
 103   int i;
 104   float wx, wy;
 105   int j, km2;
 106   float tx, ty, wyd;
 107
 108   // Data copy and bit-reversal ordering.
 109   m2 = (n >> 1);
 110   rk = 0;
 111   for (k = 0; k < n; k ++) {
 112       fftOut [rk] . mVec [0] = fftIn [k] . mVec [0];
 113       fftOut [rk] . mVec [1] = fftIn [k] . mVec [1];
 114       if (k < (n - 1)) {
 115          m = m2;
 116          while (rk >= m) {
 117            rk -= m;
 118            m >>= 1;
 119          }
 120          rk += m;
 121       }
 122   }
 123   // Perform the FFT.
 124   m2 = 1;
 125   for (m = 2; m <= n; m <<= 1) {
 126       a = sin (M_PI / m);
 127       a = 2.0f * a * a;
 128       b = sin (2.0f * M_PI / m);
 129       for (i = 0; i < n; i += m) {
 130           wx = 1.0f;
 131           wy = 0.0f;
 132           for (k = i, j = 0; j < m2; k ++, j ++) {
 133               km2 = k + m2;
 134               tx = (wx * fftOut [km2] . mVec [0]) - (wy * fftOut [km2] . mVec [1]);
 135               ty = (wx * fftOut [km2] . mVec [1]) + (wy * fftOut [km2] . mVec [0]);
 136               fftOut [km2] . mVec [0] = fftOut [k] . mVec [0] - tx;
 137               fftOut [km2] . mVec [1] = fftOut [k] . mVec [1] - ty;
 138               fftOut [k] . mVec [0] += tx;
 139               fftOut [k] . mVec [1] += ty;
 140               wyd = (a * wy) - (b * wx);
 141               wx -= (a * wx) + (b * wy);
 142               wy -= wyd;
 143           }
 144       }
 145       m2 = m;
 146   }
 147 }
 148
 149 // Performs an inverse Fast Fourier Transform.
 150 static void FftInvProc (int n, const ComplexT * fftIn, ComplexT * fftOut) {
 151   int m2, rk, k, m;
 152   float a, b;
 153   int i;
 154   float wx, wy;
 155   int j, km2;
 156   float tx, ty, wyd, invn;
 157
 158   // Data copy and bit-reversal ordering.
 159   m2 = (n >> 1);
 160   rk = 0;
 161   for (k = 0; k < n; k ++) {
 162       fftOut [rk] . mVec [0] = fftIn [k] . mVec [0];
 163       fftOut [rk] . mVec [1] = fftIn [k] . mVec [1];
 164       if (k < (n - 1)) {
 165          m = m2;
 166          while (rk >= m) {
 167            rk -= m;
 168            m >>= 1;
 169          }
 170          rk += m;
 171       }
 172   }
 173   // Perform the IFFT.
 174   m2 = 1;
 175   for (m = 2; m <= n; m <<= 1) {
 176       a = sin (M_PI / m);
 177       a = 2.0f * a * a;
 178       b = -sin (2.0f * M_PI / m);
 179       for (i = 0; i < n; i += m) {
 180           wx = 1.0f;
 181           wy = 0.0f;
 182           for (k = i, j = 0; j < m2; k ++, j ++) {
 183               km2 = k + m2;
 184               tx = (wx * fftOut [km2] . mVec [0]) - (wy * fftOut [km2] . mVec [1]);
 185               ty = (wx * fftOut [km2] . mVec [1]) + (wy * fftOut [km2] . mVec [0]);
 186               fftOut [km2] . mVec [0] = fftOut [k] . mVec [0] - tx;
 187               fftOut [km2] . mVec [1] = fftOut [k] . mVec [1] - ty;
 188               fftOut [k] . mVec [0] += tx;
 189               fftOut [k] . mVec [1] += ty;
 190               wyd = (a * wy) - (b * wx);
 191               wx -= (a * wx) + (b * wy);
 192               wy -= wyd;
 193           }
 194       }
 195       m2 = m;
 196   }
 197   // Normalize the samples.
 198   invn = 1.0f / n;
 199   for (i = 0; i < n; i ++) {
 200       fftOut [i] . mVec [0] *= invn;
 201       fftOut [i] . mVec [1] *= invn;
 202   }
 203 }
 204
 205 // Complex absolute value.
 206 static void ComplexAbs (const ComplexT * in, ComplexT * out) {
 207   out -> mVec [0] = sqrt ((in -> mVec [0] * in -> mVec [0]) + (in -> mVec [1] * in -> mVec [1]));
 208   out -> mVec [1] = 0.0f;
 209 }
 210
 211 // Complex logarithm.
 212 static void ComplexLog (const ComplexT * in, ComplexT * out) {
 213   float r, t;
 214
 215   r = sqrt ((in -> mVec [0] * in -> mVec [0]) + (in -> mVec [1] * in -> mVec [1]));
 216   t = atan2 (in -> mVec [1], in -> mVec [0]);
 217   if (t < 0.0f)
 218      t += 2.0f * M_PI;
 219   out -> mVec [0] = log (r);
 220   out -> mVec [1] = t;
 221 }
 222
 223 // Complex exponent.
 224 static void ComplexExp (const ComplexT * in, ComplexT * out) {
 225   float e;
 226
 227   e = exp (in -> mVec [0]);
 228   out -> mVec [0] = e * cos (in -> mVec [1]);
 229   out -> mVec [1] = e * sin (in -> mVec [1]);
 230 }
 231
 232 // Calculates the real cepstrum of a given impulse response.  It currently
 233 // uses a fixed cepstrum size.  To make this more robust, it should be
 234 // rewritten to handle a variable size cepstrum.
 235 static void RealCepstrum (int irSize, const float * ir, float cep [CEP_SIZE]) {
 236   ComplexT in [CEP_SIZE], out [CEP_SIZE];
 237   int index;
 238
 239   for (index = 0; index < irSize; index ++) {
 240       in [index] . mVec [0] = ir [index];
 241       in [index] . mVec [1] = 0.0f;
 242   }
 243   for (; index < CEP_SIZE; index ++) {
 244       in [index] . mVec [0] = 0.0f;
 245       in [index] . mVec [1] = 0.0f;
 246   }
 247   FftProc (CEP_SIZE, in, out);
 248   for (index = 0; index < CEP_SIZE; index ++) {
 249       ComplexAbs (& out [index], & out [index]);
 250       if (out [index] . mVec [0] < 0.000001f)
 251          out [index] . mVec [0] = 0.000001f;
 252       ComplexLog (& out [index], & in [index]);
 253   }
 254   FftInvProc (CEP_SIZE, in, out);
 255   for (index = 0; index < CEP_SIZE; index ++)
 256       cep [index] = out [index] . mVec [0];
 257 }
 258
 259 // Reconstructs the minimum-phase impulse response for a given real cepstrum.
 260 // Like the above function, this should eventually be modified to handle a
 261 // variable size cepstrum.
 262 static void MinimumPhase (const float cep [CEP_SIZE], int irSize, float * mpIr) {
 263   ComplexT in [CEP_SIZE], out [CEP_SIZE];
 264   int index;
 265
 266   in [0] . mVec [0] = cep [0];
 267   for (index = 1; index < (CEP_SIZE / 2); index ++)
 268       in [index] . mVec [0] = 2.0f * cep [index];
 269   if ((CEP_SIZE % 2) != 1) {
 270      in [index] . mVec [0] = cep [index];
 271      index ++;
 272   }
 273   for (; index < CEP_SIZE; index ++)
 274       in [index] . mVec [0] = 0.0f;
 275   for (index = 0; index < CEP_SIZE; index ++)
 276       in [index] . mVec [1] = 0.0f;
 277   FftProc (CEP_SIZE, in, out);
 278   for (index = 0; index < CEP_SIZE; index ++)
 279       ComplexExp (& out [index], & in [index]);
 280   FftInvProc (CEP_SIZE, in, out);
 281   for (index = 0; index < irSize; index ++)
 282       mpIr [index] = out [index] . mVec [0];
 283 }
 284
 285 // Calculate the left-ear time delay using a spherical head model.
 286 static float CalcLTD (float ev, float az, float rad, float dist) {
 287   float azp, dlp, l, al;
 288
 289   azp = asin (cos (ev) * sin (az));
 290   dlp = sqrt ((dist * dist) + (rad * rad) + (2.0f * dist * rad * sin (azp)));
 291   l = sqrt ((dist * dist) - (rad * rad));
 292   al = (0.5f * M_PI) + azp;
 293   if (dlp > l)
 294      dlp = l + (rad * (al - acos (rad / dist)));
 295   return (dlp / 343.3f);
 296 }
 297
 298 // Read a 16-bit little-endian integer from a file and convert it to a 32-bit
 299 // floating-point value in the range of -1.0 to 1.0.
 300 static int ReadInt16LeAsFloat32 (const char * fileName, FILE * fp, float * val) {
 301   ALubyte vb [2];
 302   ALushort vw;
 303
 304   if (fread (vb, 1, sizeof (vb), fp) != sizeof (vb)) {
 305      fclose (fp);
 306      fprintf (stderr, "Error reading from file, '%s'.\n", fileName);
 307      return (0);
 308   }
 309   vw = (((unsigned short) vb [1]) << 8) | vb [0];
 310   (* val) = ((short) vw) / 32768.0f;
 311   return (1);
 312 }
 313
 314 // Write a string to a file.
 315 static int WriteString (const char * val, const char * fileName, FILE * fp) {
 316   size_t len;
 317
 318   len = strlen (val);
 319   if (fwrite (val, 1, len, fp) != len) {
 320      fclose (fp);
 321      fprintf (stderr, "Error writing to file, '%s'.\n", fileName);
 322      return (0);
 323   }
 324   return (1);
 325 }
 326
 327 // Write a 32-bit floating-point value in the range of -1.0 to 1.0 to a file
 328 // as a 16-bit little-endian integer.
 329 static int WriteFloat32AsInt16Le (float val, const char * fileName, FILE * fp) {
 330   ALshort vw;
 331   ALubyte vb [2];
 332
 333   vw = (short) round (32767.0f * val);
 334   vb [0] =  vw       & 0x00FF;
 335   vb [1] = (vw >> 8) & 0x00FF;
 336   if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) {
 337      fclose (fp);
 338      fprintf (stderr, "Error writing to file, '%s'.\n", fileName);
 339      return (0);
 340   }
 341   return (1);
 342 }
 343
 344 // Write a 32-bit little-endian unsigned integer to a file.
 345 static int WriteUInt32Le (ALuint val, const char * fileName, FILE * fp) {
 346   ALubyte vb [4];
 347
 348   vb [0] =  val        & 0x000000FF;
 349   vb [1] = (val >>  8) & 0x000000FF;
 350   vb [2] = (val >> 16) & 0x000000FF;
 351   vb [3] = (val >> 24) & 0x000000FF;
 352   if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) {
 353      fclose (fp);
 354      fprintf (stderr, "Error writing to file, '%s'.\n", fileName);
 355      return (0);
 356   }
 357   return (1);
 358 }
 359
 360 // Write a 16-bit little-endian unsigned integer to a file.
 361 static int WriteUInt16Le (ALushort val, const char * fileName, FILE * fp) {
 362   ALubyte vb [2];
 363
 364   vb [0] =  val        & 0x00FF;
 365   vb [1] = (val >>  8) & 0x00FF;
 366   if (fwrite (vb, 1, sizeof (vb), fp) != sizeof (vb)) {
 367      fclose (fp);
 368      fprintf (stderr, "Error writing to file, '%s'.\n", fileName);
 369      return (0);
 370   }
 371   return (1);
 372 }
 373
 374 // Write an 8-bit unsigned integer to a file.
 375 static int WriteUInt8 (ALubyte val, const char * fileName, FILE * fp) {
 376   if (fwrite (& val, 1, sizeof (val), fp) != sizeof (val)) {
 377      fclose (fp);
 378      fprintf (stderr, "Error writing to file, '%s'.\n", fileName);
 379      return (0);
 380   }
 381   return (1);
 382 }
 383
 384 // Load the MIT HRIRs.  This loads the entire diffuse or compact set starting
 385 // counter-clockwise up at the bottom elevation and clockwise at the forward
 386 // azimuth.
 387 static int LoadMitHrirs (const char * baseName, HrirDataT * hData) {
 388   const int EV_ANGLE [MIT_EV_COUNT] = {
 389     -90, -80, -70, -60, -50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50, 60, 70, 80, 90
 390   };
 391   int e, a;
 392   char fileName [1024];
 393   FILE * fp = NULL;
 394   int j0, j1, i;
 395   float s;
 396
 397   for (e = MIT_EV_START; e < MIT_EV_COUNT; e ++) {
 398       for (a = 0; a < MIT_AZ_COUNT [e]; a ++) {
 399           // The data packs the first 180 degrees in the left channel, and
 400           // the last 180 degrees in the right channel.
 401           if (round ((360.0f / MIT_AZ_COUNT [e]) * a) > 180.0f)
 402              break;
 403           // Determine which file to open.
 404           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));
 405           if ((fp = fopen (fileName, "rb")) == NULL) {
 406              fprintf (stderr, "Could not open file, '%s'.\n", fileName);
 407              return (0);
 408           }
 409           // Assuming they have not changed format, skip the .WAV header.
 410           fseek (fp, 44, SEEK_SET);
 411           // Map the left and right channels to their appropriate azimuth
 412           // offsets.
 413           j0 = (MIT_EV_OFFSET [e] + a) * MIT_IR_SIZE;
 414           j1 = (MIT_EV_OFFSET [e] + ((MIT_AZ_COUNT [e] - a) % MIT_AZ_COUNT [e])) * MIT_IR_SIZE;
 415           // Read in the data, converting it to floating-point.
 416           for (i = 0; i < MIT_IR_SIZE; i ++) {
 417               if (! ReadInt16LeAsFloat32 (fileName, fp, & s))
 418                  return (0);
 419               hData -> mHrirs [j0 + i] = s;
 420               if (! ReadInt16LeAsFloat32 (fileName, fp, & s))
 421                  return (0);
 422               hData -> mHrirs [j1 + i] = s;
 423           }
 424           fclose (fp);
 425       }
 426   }
 427   return (1);
 428 }
 429
 430 // Performs the minimum phase reconstruction for a given HRIR data set.  The
 431 // cepstrum size should be made configureable at some point in the future.
 432 static void ReconstructHrirs (int minIrSize, HrirDataT * hData) {
 433   int start, end, step, j;
 434   float cep [CEP_SIZE];
 435
 436   start = hData -> mEvOffset [hData -> mEvStart];
 437   end = hData -> mIrCount;
 438   step = hData -> mIrSize;
 439   for (j = start; j < end; j ++) {
 440       RealCepstrum (step, & hData -> mHrirs [j * step], cep);
 441       MinimumPhase (cep, minIrSize, & hData -> mHrirs [j * minIrSize]);
 442   }
 443   hData -> mIrSize = minIrSize;
 444 }
 445
 446 // Renormalize the entire HRIR data set, and attenutate it slightly.
 447 static void RenormalizeHrirs (const HrirDataT * hData) {
 448   int step, start, end;
 449   float norm;
 450   int j, i;
 451
 452   step = hData -> mIrSize;
 453   start = hData -> mEvOffset [hData -> mEvStart] * step;
 454   end = hData -> mIrCount * step;
 455   norm = 0.0f;
 456   for (j = start; j < end; j += step) {
 457       for (i = 0; i < step; i ++) {
 458           if (fabs (hData -> mHrirs [j + i]) > norm)
 459              norm = fabs (hData -> mHrirs [j + i]);
 460       }
 461   }
 462   if (norm > 0.000001f)
 463      norm = 1.0f / norm;
 464   norm *= 0.95f;
 465   for (j = start; j < end; j += step) {
 466       for (i = 0; i < step; i ++)
 467           hData -> mHrirs [j + i] *= norm;
 468   }
 469 }
 470
 471 // Given an elevation offset and azimuth, calculates two offsets for
 472 // addressing the HRIRs buffer and their interpolation factor.
 473 static void CalcAzIndices (const HrirDataT * hData, int oi, float az, int * j0, int * j1, float * jf) {
 474   int ai;
 475
 476   az = fmod ((2.0f * M_PI) + az, 2.0f * M_PI) * hData -> mAzCount [oi] / (2.0f * M_PI);
 477   ai = (int) az;
 478   az -= ai;
 479   (* j0) = hData -> mEvOffset [oi] + ai;
 480   (* j1) = hData -> mEvOffset [oi] + ((ai + 1) % hData -> mAzCount [oi]);
 481   (* jf) = az;
 482 }
 483
 484 // Perform a linear interpolation.
 485 static float Lerp (float a, float b, float f) {
 486   return (a + (f * (b - a)));
 487 }
 488
 489 // Attempt to synthesize any missing HRIRs at the bottom elevations.  Right
 490 // now this just blends the lowest elevation HRIRs together and applies some
 491 // attenuates and high frequency damping.  It's not a realistic model to use,
 492 // but it is simple.
 493 static void SynthesizeHrirs (HrirDataT * hData) {
 494   int step, oi, i, a, j, e;
 495   float of;
 496   int j0, j1;
 497   float jf;
 498   float lp [4], s0, s1;
 499
 500   if (hData -> mEvStart <= 0)
 501      return;
 502   step = hData -> mIrSize;
 503   oi = hData -> mEvStart;
 504   for (i = 0; i < step; i ++)
 505       hData -> mHrirs [i] = 0.0f;
 506   for (a = 0; a < hData -> mAzCount [oi]; a ++) {
 507       j = (hData -> mEvOffset [oi] + a) * step;
 508       for (i = 0; i < step; i ++)
 509           hData -> mHrirs [i] += hData -> mHrirs [j + i] / hData -> mAzCount [oi];
 510   }
 511   for (e = 1; e < hData -> mEvStart; e ++) {
 512       of = ((float) e) / hData -> mEvStart;
 513       for (a = 0; a < hData -> mAzCount [e]; a ++) {
 514           j = (hData -> mEvOffset [e] + a) * step;
 515           CalcAzIndices (hData, oi, a * 2.0f * M_PI / hData -> mAzCount [e], & j0, & j1, & jf);
 516           j0 *= step;
 517           j1 *= step;
 518           lp [0] = 0.0f;
 519           lp [1] = 0.0f;
 520           lp [2] = 0.0f;
 521           lp [3] = 0.0f;
 522           for (i = 0; i < step; i ++) {
 523               s0 = hData -> mHrirs [i];
 524               s1 = Lerp (hData -> mHrirs [j0 + i], hData -> mHrirs [j1 + i], jf);
 525               s0 = Lerp (s0, s1, of);
 526               lp [0] = Lerp (s0, lp [0], 0.15f - (0.15f * of));
 527               lp [1] = Lerp (lp [0], lp [1], 0.15f - (0.15f * of));
 528               lp [2] = Lerp (lp [1], lp [2], 0.15f - (0.15f * of));
 529               lp [3] = Lerp (lp [2], lp [3], 0.15f - (0.15f * of));
 530               hData -> mHrirs [j + i] = lp [3];
 531           }
 532       }
 533   }
 534   lp [0] = 0.0f;
 535   lp [1] = 0.0f;
 536   lp [2] = 0.0f;
 537   lp [3] = 0.0f;
 538   for (i = 0; i < step; i ++) {
 539       s0 = hData -> mHrirs [i];
 540       lp [0] = Lerp (s0, lp [0], 0.15f);
 541       lp [1] = Lerp (lp [0], lp [1], 0.15f);
 542       lp [2] = Lerp (lp [1], lp [2], 0.15f);
 543       lp [3] = Lerp (lp [2], lp [3], 0.15f);
 544       hData -> mHrirs [i] = lp [3];
 545   }
 546   hData -> mEvStart = 0;
 547 }
 548
 549 // Calculate the effective head-related time delays for the each HRIR, now
 550 // that they are minimum-phase.
 551 static void CalculateHrtds (HrirDataT * hData) {
 552   float minHrtd, maxHrtd;
 553   int e, a, j;
 554   float t;
 555
 556   minHrtd = 1000.0f;
 557   maxHrtd = -1000.0f;
 558   for (e = 0; e < hData -> mEvCount; e ++) {
 559       for (a = 0; a < hData -> mAzCount [e]; a ++) {
 560           j = hData -> mEvOffset [e] + a;
 561           t = CalcLTD ((-90.0f + (e * 180.0f / (hData -> mEvCount - 1))) * M_PI / 180.0f,
 562                        (a * 360.0f / hData -> mAzCount [e]) * M_PI / 180.0f,
 563                        hData -> mRadius, hData -> mDistance);
 564           hData -> mHrtds [j] = t;
 565           if (t > maxHrtd)
 566              maxHrtd = t;
 567           if (t < minHrtd)
 568              minHrtd = t;
 569       }
 570   }
 571   maxHrtd -= minHrtd;
 572   for (j = 0; j < hData -> mIrCount; j ++)
 573       hData -> mHrtds [j] -= minHrtd;
 574   hData -> mMaxHrtd = maxHrtd;
 575 }
 576
 577 // Save the OpenAL Soft HRTF data set.
 578 static int SaveMhr (const HrirDataT * hData, const char * fileName) {
 579   FILE * fp = NULL;
 580   int e, step, end, j, i;
 581
 582   if ((fp = fopen (fileName, "wb")) == NULL) {
 583      fprintf (stderr, "Could not create file, '%s'.\n", fileName);
 584      return (0);
 585   }
 586   if (! WriteString (MHR_FORMAT, fileName, fp))
 587      return (0);
 588   if (! WriteUInt32Le ((ALuint) hData -> mIrRate, fileName, fp))
 589      return (0);
 590   if (! WriteUInt16Le ((ALushort) hData -> mIrCount, fileName, fp))
 591      return (0);
 592   if (! WriteUInt16Le ((ALushort) hData -> mIrSize, fileName, fp))
 593      return (0);
 594   if (! WriteUInt8 ((ALubyte) hData -> mEvCount, fileName, fp))
 595      return (0);
 596   for (e = 0; e < hData -> mEvCount; e ++) {
 597       if (! WriteUInt16Le ((ALushort) hData -> mEvOffset [e], fileName, fp))
 598          return (0);
 599   }
 600   step = hData -> mIrSize;
 601   end = hData -> mIrCount * step;
 602   for (j = 0; j < end; j += step) {
 603       for (i = 0; i < step; i ++) {
 604           if (! WriteFloat32AsInt16Le (hData -> mHrirs [j + i], fileName, fp))
 605              return (0);
 606       }
 607   }
 608   for (j = 0; j < hData -> mIrCount; j ++) {
 609       i = (int) round (44100.0f * hData -> mHrtds [j]);
 610       if (i > 127)
 611          i = 127;
 612       if (! WriteUInt8 ((ALubyte) i, fileName, fp))
 613          return (0);
 614   }
 615   fclose (fp);
 616   return (1);
 617 }
 618
 619 // Save the OpenAL Soft built-in table.
 620 static int SaveTab (const HrirDataT * hData, const char * fileName) {
 621   FILE * fp = NULL;
 622   int step, end, j, i;
 623   char text [16];
 624
 625   if ((fp = fopen (fileName, "wb")) == NULL) {
 626      fprintf (stderr, "Could not create file, '%s'.\n", fileName);
 627      return (0);
 628   }
 629   if (! WriteString ("/* This data is Copyright 1994 by the MIT Media Laboratory. It is provided free\n"
 630                      " * with no restrictions on use, provided the authors are cited when the data is\n"
 631                      " * used in any research or commercial application. */\n"
 632                      "/* Bill Gardner <billg@media.mit.edu> and Keith Martin <kdm@media.mit.edu> */\n"
 633                      "\n"
 634                      "    /* HRIR Coefficients */\n"
 635                      "    {\n", fileName, fp))
 636      return (0);
 637   step = hData -> mIrSize;
 638   end = hData -> mIrCount * step;
 639   for (j = 0; j < end; j += step) {
 640       if (! WriteString ("        { ", fileName, fp))
 641          return (0);
 642       for (i = 0; i < step; i ++) {
 643           snprintf (text, 15, "%+d, ", (int) round (32767.0f * hData -> mHrirs [j + i]));
 644           if (! WriteString (text, fileName, fp))
 645              return (0);
 646       }
 647       if (! WriteString ("},\n", fileName, fp))
 648          return (0);
 649   }
 650   if (! WriteString ("    },\n"
 651                      "\n"
 652                      "    /* HRIR Delays */\n"
 653                      "    { ", fileName, fp))
 654      return (0);
 655   for (j = 0; j < hData -> mIrCount; j ++) {
 656       snprintf (text, 15, "%d, ", (int) round (44100.0f * hData -> mHrtds [j]));
 657       if (! WriteString (text, fileName, fp))
 658          return (0);
 659   }
 660   if (! WriteString ("}\n", fileName, fp))
 661      return (0);
 662   fclose (fp);
 663   return (1);
 664 }
 665
 666 // Loads and processes an MIT data set.  At present, the HRIR and HRTD data
 667 // is loaded and processed in a static buffer.  That should change to using
 668 // heap allocated memory in the future.  A cleanup function will then be
 669 // required.
 670 static int MakeMit(const char *baseInName, HrirDataT *hData)
 671 {
 672     static float hrirs[MIT_IR_COUNT * MIT_IR_SIZE];
 673     static float hrtds[MIT_IR_COUNT];
 674
 675     hData->mIrRate = MIT_IR_RATE;
 676     hData->mIrCount = MIT_IR_COUNT;
 677     hData->mIrSize = MIT_IR_SIZE;
 678     hData->mEvCount = MIT_EV_COUNT;
 679     hData->mEvStart = MIT_EV_START;
 680     hData->mEvOffset = MIT_EV_OFFSET;
 681     hData->mAzCount = MIT_AZ_COUNT;
 682     hData->mRadius = MIT_RADIUS;
 683     hData->mDistance = MIT_DISTANCE;
 684     hData->mHrirs = hrirs;
 685     hData->mHrtds = hrtds;
 686     fprintf(stderr, "Loading base HRIR data...\n");
 687     if(!LoadMitHrirs(baseInName, hData))
 688         return 0;
 689     fprintf(stderr, "Performing minimum phase reconstruction and truncation...\n");
 690     ReconstructHrirs(MIN_IR_SIZE, hData);
 691     fprintf(stderr, "Renormalizing minimum phase HRIR data...\n");
 692     RenormalizeHrirs(hData);
 693     fprintf(stderr, "Synthesizing missing elevations...\n");
 694     SynthesizeHrirs(hData);
 695     fprintf(stderr, "Calculating impulse delays...\n");
 696     CalculateHrtds(hData);
 697     return 1;
 698 }
 699
 700 // Simple dispatch.  Provided a command, the path to the MIT set of choice,
 701 // and an optional output filename, this will produce an OpenAL Soft
 702 // compatible HRTF set in the chosen format.
 703 int main(int argc, char *argv[])
 704 {
 705     char baseName[1024];
 706     const char *outName = NULL;
 707     HrirDataT hData;
 708
 709     if(argc < 3 || strcmp(argv [1], "-h") == 0 || strcmp (argv [1], "--help") == 0)
 710     {
 711         fprintf(stderr, "Usage:  %s <command> <path of MIT set> [ <output file> ]\n\n", argv[0]);
 712         fprintf(stderr, "Commands:\n");
 713         fprintf(stderr, " -m, --make-mhr   Makes an OpenAL Soft compatible HRTF data set.\n");
 714         fprintf(stderr, "                  Defaults output to:  ./oal_soft_hrtf_44100.mhr\n");
 715         fprintf(stderr, " -t, --make-tab   Makes the built-in table used when compiling OpenAL Soft.\n");
 716         fprintf(stderr, "                  Defaults output to:  ./hrtf_tables.inc\n");
 717         fprintf(stderr, " -h, --help       Displays this help information.\n");
 718         return 0;
 719     }
 720
 721     snprintf(baseName, sizeof(baseName), "%s/elev", argv[2]);
 722     if(strcmp(argv[1], "-m") == 0 || strcmp(argv[1], "--make-mhr") == 0)
 723     {
 724         if(argc > 3)
 725             outName = argv[3];
 726         else
 727             outName = "./oal_soft_hrtf_44100.mhr";
 728         if(!MakeMit(baseName, &hData))
 729             return -1;
 730         fprintf(stderr, "Creating data set file...\n");
 731         if(!SaveMhr(&hData, outName))
 732             return -1;
 733     }
 734     else if(strcmp(argv[1], "-t") == 0 || strcmp(argv[1], "--make-tab") == 0)
 735     {
 736         if(argc > 3)
 737             outName = argv[3];
 738         else
 739             outName = "./hrtf_tables.inc";
 740         if(!MakeMit(baseName, &hData))
 741             return -1;
 742         fprintf(stderr, "Creating table file...\n");
 743         if(!SaveTab(&hData, outName))
 744             return -1;
 745     }
 746     else
 747     {
 748         fprintf(stderr, "Invalid command '%s'\n", argv[1]);
 749         return -1;
 750     }
 751     fprintf(stderr, "Done.\n");
 752     return 0;
 753 }