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