Alc/hrtf.c

   1 /**
   2  * OpenAL cross platform audio library
   3  * Copyright (C) 2011 by Chris Robinson
   4  * This library is free software; you can redistribute it and/or
   5  *  modify it under the terms of the GNU Library General Public
   6  *  License as published by the Free Software Foundation; either
   7  *  version 2 of the License, or (at your option) any later version.
   8  *
   9  * This library is distributed in the hope that it will be useful,
  10  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12  *  Library General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU Library General Public
  15  *  License along with this library; if not, write to the
  16  *  Free Software Foundation, Inc., 59 Temple Place - Suite 330,
  17  *  Boston, MA  02111-1307, USA.
  18  * Or go to http://www.gnu.org/copyleft/lgpl.html
  19  */
  20
  21 #include "config.h"
  22
  23 #include <stdlib.h>
  24 #include <ctype.h>
  25
  26 #include "AL/al.h"
  27 #include "AL/alc.h"
  28 #include "alMain.h"
  29 #include "alSource.h"
  30
  31 /* External HRTF file format (LE byte order):
  32  *
  33  * ALchar   magic[8] = "MinPHR00";
  34  * ALuint   sampleRate;
  35  *
  36  * ALushort hrirCount; // Required value: 828
  37  * ALushort hrirSize;  // Required value: 32
  38  * ALubyte  evCount;   // Required value: 19
  39  *
  40  * ALushort evOffset[evCount]; // Required values:
  41  *   { 0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827 }
  42  *
  43  * ALshort coefficients[hrirCount][hrirSize];
  44  * ALubyte delays[hrirCount]; // Element values must not exceed 127
  45  */
  46
  47 static const ALchar magicMarker[8] = "MinPHR00";
  48
  49 #define HRIR_COUNT 828
  50 #define ELEV_COUNT 19
  51
  52 static const ALushort evOffset[ELEV_COUNT] = { 0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827 };
  53 static const ALubyte azCount[ELEV_COUNT] = { 1, 12, 24, 36, 45, 56, 60, 72, 72, 72, 72, 72, 60, 56, 45, 36, 24, 12, 1 };
  54
  55
  56 static const struct Hrtf {
  57     ALuint sampleRate;
  58     ALshort coeffs[HRIR_COUNT][HRIR_LENGTH];
  59     ALubyte delays[HRIR_COUNT];
  60 } DefaultHrtf = {
  61     44100,
  62 #include "hrtf_tables.inc"
  63 };
  64
  65 static struct Hrtf *LoadedHrtfs = NULL;
  66 static ALuint NumLoadedHrtfs = 0;
  67
  68
  69 // Calculate the elevation indices given the polar elevation in radians.
  70 // This will return two indices between 0 and (ELEV_COUNT-1) and an
  71 // interpolation factor between 0.0 and 1.0.
  72 static void CalcEvIndices(ALfloat ev, ALuint *evidx, ALfloat *evmu)
  73 {
  74     ev = (F_PI_2 + ev) * (ELEV_COUNT-1) / F_PI;
  75     evidx[0] = (ALuint)ev;
  76     evidx[1] = minu(evidx[0] + 1, ELEV_COUNT-1);
  77     *evmu = ev - evidx[0];
  78 }
  79
  80 // Calculate the azimuth indices given the polar azimuth in radians.  This
  81 // will return two indices between 0 and (azCount [ei] - 1) and an
  82 // interpolation factor between 0.0 and 1.0.
  83 static void CalcAzIndices(ALuint evidx, ALfloat az, ALuint *azidx, ALfloat *azmu)
  84 {
  85     az = (F_PI*2.0f + az) * azCount[evidx] / (F_PI*2.0f);
  86     azidx[0] = (ALuint)az % azCount[evidx];
  87     azidx[1] = (azidx[0] + 1) % azCount[evidx];
  88     *azmu = az - floor(az);
  89 }
  90
  91 // Calculates the normalized HRTF transition factor (delta) from the changes
  92 // in gain and listener to source angle between updates.  The result is a
  93 // normalized delta factor than can be used to calculate moving HRIR stepping
  94 // values.
  95 ALfloat CalcHrtfDelta(ALfloat oldGain, ALfloat newGain, const ALfloat olddir[3], const ALfloat newdir[3])
  96 {
  97     ALfloat gainChange, angleChange;
  98
  99     // Calculate the normalized dB gain change.
 100     newGain = maxf(newGain, 0.0001f);
 101     oldGain = maxf(oldGain, 0.0001f);
 102     gainChange = aluFabs(log10(newGain / oldGain) / log10(0.0001f));
 103
 104     // Calculate the normalized listener to source angle change when there is
 105     // enough gain to notice it.
 106     angleChange = 0.0f;
 107     if(gainChange > 0.0001f || newGain > 0.0001f)
 108     {
 109         // No angle change when the directions are equal or degenerate (when
 110         // both have zero length).
 111         if(newdir[0]-olddir[0] || newdir[1]-olddir[1] || newdir[2]-olddir[2])
 112             angleChange = aluAcos(olddir[0]*newdir[0] +
 113                                   olddir[1]*newdir[1] +
 114                                   olddir[2]*newdir[2]) / F_PI;
 115
 116     }
 117
 118     // Use the largest of the two changes for the delta factor, and apply a
 119     // significance shaping function to it.
 120     return clampf(angleChange*2.0f, gainChange*2.0f, 1.0f);
 121 }
 122
 123 // Calculates static HRIR coefficients and delays for the given polar
 124 // elevation and azimuth in radians.  Linear interpolation is used to
 125 // increase the apparent resolution of the HRIR dataset.  The coefficients
 126 // are also normalized and attenuated by the specified gain.
 127 void GetLerpedHrtfCoeffs(const struct Hrtf *Hrtf, ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat (*coeffs)[2], ALuint *delays)
 128 {
 129     ALuint evidx[2], azidx[2];
 130     ALfloat mu[3];
 131     ALuint lidx[4], ridx[4];
 132     ALuint i;
 133
 134     // Claculate elevation indices and interpolation factor.
 135     CalcEvIndices(elevation, evidx, &mu[2]);
 136
 137     // Calculate azimuth indices and interpolation factor for the first
 138     // elevation.
 139     CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
 140
 141     // Calculate the first set of linear HRIR indices for left and right
 142     // channels.
 143     lidx[0] = evOffset[evidx[0]] + azidx[0];
 144     lidx[1] = evOffset[evidx[0]] + azidx[1];
 145     ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
 146     ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
 147
 148     // Calculate azimuth indices and interpolation factor for the second
 149     // elevation.
 150     CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
 151
 152     // Calculate the second set of linear HRIR indices for left and right
 153     // channels.
 154     lidx[2] = evOffset[evidx[1]] + azidx[0];
 155     lidx[3] = evOffset[evidx[1]] + azidx[1];
 156     ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
 157     ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
 158
 159     // Calculate the normalized and attenuated HRIR coefficients using linear
 160     // interpolation when there is enough gain to warrant it.  Zero the
 161     // coefficients if gain is too low.
 162     if(gain > 0.0001f)
 163     {
 164         ALdouble scale = gain * (1.0/32767.0);
 165         for(i = 0;i < HRIR_LENGTH;i++)
 166         {
 167             coeffs[i][0] = lerp(lerp(Hrtf->coeffs[lidx[0]][i], Hrtf->coeffs[lidx[1]][i], mu[0]),
 168                                 lerp(Hrtf->coeffs[lidx[2]][i], Hrtf->coeffs[lidx[3]][i], mu[1]),
 169                                 mu[2]) * scale;
 170             coeffs[i][1] = lerp(lerp(Hrtf->coeffs[ridx[0]][i], Hrtf->coeffs[ridx[1]][i], mu[0]),
 171                                 lerp(Hrtf->coeffs[ridx[2]][i], Hrtf->coeffs[ridx[3]][i], mu[1]),
 172                                 mu[2]) * scale;
 173         }
 174     }
 175     else
 176     {
 177         for(i = 0;i < HRIR_LENGTH;i++)
 178         {
 179             coeffs[i][0] = 0.0f;
 180             coeffs[i][1] = 0.0f;
 181         }
 182     }
 183
 184     // Calculate the HRIR delays using linear interpolation.
 185     delays[0] = (ALuint)(lerp(lerp(Hrtf->delays[lidx[0]], Hrtf->delays[lidx[1]], mu[0]),
 186                               lerp(Hrtf->delays[lidx[2]], Hrtf->delays[lidx[3]], mu[1]),
 187                               mu[2]) * 65536.0f);
 188     delays[1] = (ALuint)(lerp(lerp(Hrtf->delays[ridx[0]], Hrtf->delays[ridx[1]], mu[0]),
 189                               lerp(Hrtf->delays[ridx[2]], Hrtf->delays[ridx[3]], mu[1]),
 190                               mu[2]) * 65536.0f);
 191 }
 192
 193 // Calculates the moving HRIR target coefficients, target delays, and
 194 // stepping values for the given polar elevation and azimuth in radians.
 195 // Linear interpolation is used to increase the apparent resolution of the
 196 // HRIR dataset.  The coefficients are also normalized and attenuated by the
 197 // specified gain.  Stepping resolution and count is determined using the
 198 // given delta factor between 0.0 and 1.0.
 199 ALuint GetMovingHrtfCoeffs(const struct Hrtf *Hrtf, ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat delta, ALint counter, ALfloat (*coeffs)[2], ALuint *delays, ALfloat (*coeffStep)[2], ALint *delayStep)
 200 {
 201     ALuint evidx[2], azidx[2];
 202     ALuint lidx[4], ridx[4];
 203     ALfloat left, right;
 204     ALfloat mu[3];
 205     ALfloat step;
 206     ALuint i;
 207
 208     // Claculate elevation indices and interpolation factor.
 209     CalcEvIndices(elevation, evidx, &mu[2]);
 210
 211     // Calculate azimuth indices and interpolation factor for the first
 212     // elevation.
 213     CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
 214
 215     // Calculate the first set of linear HRIR indices for left and right
 216     // channels.
 217     lidx[0] = evOffset[evidx[0]] + azidx[0];
 218     lidx[1] = evOffset[evidx[0]] + azidx[1];
 219     ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
 220     ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
 221
 222     // Calculate azimuth indices and interpolation factor for the second
 223     // elevation.
 224     CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
 225
 226     // Calculate the second set of linear HRIR indices for left and right
 227     // channels.
 228     lidx[2] = evOffset[evidx[1]] + azidx[0];
 229     lidx[3] = evOffset[evidx[1]] + azidx[1];
 230     ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
 231     ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
 232
 233     // Calculate the stepping parameters.
 234     delta = maxf(floor(delta*(Hrtf->sampleRate*0.015f) + 0.5), 1.0f);
 235     step = 1.0f / delta;
 236
 237     // Calculate the normalized and attenuated target HRIR coefficients using
 238     // linear interpolation when there is enough gain to warrant it.  Zero
 239     // the target coefficients if gain is too low.  Then calculate the
 240     // coefficient stepping values using the target and previous running
 241     // coefficients.
 242     if(gain > 0.0001f)
 243     {
 244         ALdouble scale = gain * (1.0/32767.0);
 245         for(i = 0;i < HRIR_LENGTH;i++)
 246         {
 247             left = coeffs[i][0] - (coeffStep[i][0] * counter);
 248             right = coeffs[i][1] - (coeffStep[i][1] * counter);
 249
 250             coeffs[i][0] = lerp(lerp(Hrtf->coeffs[lidx[0]][i], Hrtf->coeffs[lidx[1]][i], mu[0]),
 251                                 lerp(Hrtf->coeffs[lidx[2]][i], Hrtf->coeffs[lidx[3]][i], mu[1]),
 252                                 mu[2]) * scale;
 253             coeffs[i][1] = lerp(lerp(Hrtf->coeffs[ridx[0]][i], Hrtf->coeffs[ridx[1]][i], mu[0]),
 254                                 lerp(Hrtf->coeffs[ridx[2]][i], Hrtf->coeffs[ridx[3]][i], mu[1]),
 255                                 mu[2]) * scale;
 256
 257             coeffStep[i][0] = step * (coeffs[i][0] - left);
 258             coeffStep[i][1] = step * (coeffs[i][1] - right);
 259         }
 260     }
 261     else
 262     {
 263         for(i = 0;i < HRIR_LENGTH;i++)
 264         {
 265             left = coeffs[i][0] - (coeffStep[i][0] * counter);
 266             right = coeffs[i][1] - (coeffStep[i][1] * counter);
 267
 268             coeffs[i][0] = 0.0f;
 269             coeffs[i][1] = 0.0f;
 270
 271             coeffStep[i][0] = step * -left;
 272             coeffStep[i][1] = step * -right;
 273         }
 274     }
 275
 276     // Calculate the HRIR delays using linear interpolation.  Then calculate
 277     // the delay stepping values using the target and previous running
 278     // delays.
 279     left = delays[0] - (delayStep[0] * counter);
 280     right = delays[1] - (delayStep[1] * counter);
 281
 282     delays[0] = (ALuint)(lerp(lerp(Hrtf->delays[lidx[0]], Hrtf->delays[lidx[1]], mu[0]),
 283                               lerp(Hrtf->delays[lidx[2]], Hrtf->delays[lidx[3]], mu[1]),
 284                               mu[2]) * 65536.0f);
 285     delays[1] = (ALuint)(lerp(lerp(Hrtf->delays[ridx[0]], Hrtf->delays[ridx[1]], mu[0]),
 286                               lerp(Hrtf->delays[ridx[2]], Hrtf->delays[ridx[3]], mu[1]),
 287                               mu[2]) * 65536.0f);
 288
 289     delayStep[0] = (ALint)(step * (delays[0] - left));
 290     delayStep[1] = (ALint)(step * (delays[1] - right));
 291
 292     // The stepping count is the number of samples necessary for the HRIR to
 293     // complete its transition.  The mixer will only apply stepping for this
 294     // many samples.
 295     return (ALuint)delta;
 296 }
 297
 298 const struct Hrtf *GetHrtf(ALCdevice *device)
 299 {
 300     if(device->FmtChans == DevFmtStereo)
 301     {
 302         ALuint i;
 303         for(i = 0;i < NumLoadedHrtfs;i++)
 304         {
 305             if(device->Frequency == LoadedHrtfs[i].sampleRate)
 306                 return &LoadedHrtfs[i];
 307         }
 308         if(device->Frequency == DefaultHrtf.sampleRate)
 309             return &DefaultHrtf;
 310     }
 311     ERR("Incompatible format: %s %uhz\n",
 312         DevFmtChannelsString(device->FmtChans), device->Frequency);
 313     return NULL;
 314 }
 315
 316 void InitHrtf(void)
 317 {
 318     char *fnamelist=NULL, *next=NULL;
 319     const char *val;
 320
 321     if(ConfigValueStr(NULL, "hrtf_tables", &val))
 322         next = fnamelist = strdup(val);
 323     while(next && *next)
 324     {
 325         const ALubyte maxDelay = SRC_HISTORY_LENGTH-1;
 326         struct Hrtf newdata;
 327         ALboolean failed;
 328         ALchar magic[9];
 329         ALsizei i, j;
 330         char *fname;
 331         FILE *f;
 332
 333         fname = next;
 334         next = strchr(fname, ',');
 335         if(next)
 336         {
 337             while(next != fname)
 338             {
 339                 next--;
 340                 if(!isspace(*next))
 341                 {
 342                     *(next++) = '\0';
 343                     break;
 344                 }
 345             }
 346             while(isspace(*next) || *next == ',')
 347                 next++;
 348         }
 349
 350         if(!fname[0])
 351             continue;
 352         TRACE("Loading %s\n", fname);
 353         f = fopen(fname, "rb");
 354         if(f == NULL)
 355         {
 356             ERR("Could not open %s\n", fname);
 357             continue;
 358         }
 359
 360         failed = AL_FALSE;
 361         if(fread(magic, 1, sizeof(magicMarker), f) != sizeof(magicMarker))
 362         {
 363             ERR("Failed to read magic marker\n");
 364             failed = AL_TRUE;
 365         }
 366         else if(memcmp(magic, magicMarker, sizeof(magicMarker)) != 0)
 367         {
 368             magic[8] = 0;
 369             ERR("Invalid magic marker: \"%s\"\n", magic);
 370             failed = AL_TRUE;
 371         }
 372
 373         if(!failed)
 374         {
 375             ALushort hrirCount, hrirSize;
 376             ALubyte  evCount;
 377
 378             newdata.sampleRate  = fgetc(f);
 379             newdata.sampleRate |= fgetc(f)<<8;
 380             newdata.sampleRate |= fgetc(f)<<16;
 381             newdata.sampleRate |= fgetc(f)<<24;
 382
 383             hrirCount  = fgetc(f);
 384             hrirCount |= fgetc(f)<<8;
 385
 386             hrirSize  = fgetc(f);
 387             hrirSize |= fgetc(f)<<8;
 388
 389             evCount = fgetc(f);
 390
 391             if(hrirCount != HRIR_COUNT || hrirSize != HRIR_LENGTH || evCount != ELEV_COUNT)
 392             {
 393                 ERR("Unsupported value: hrirCount=%d (%d), hrirSize=%d (%d), evCount=%d (%d)\n",
 394                     hrirCount, HRIR_COUNT, hrirSize, HRIR_LENGTH, evCount, ELEV_COUNT);
 395                 failed = AL_TRUE;
 396             }
 397         }
 398
 399         if(!failed)
 400         {
 401             for(i = 0;i < HRIR_COUNT;i++)
 402             {
 403                 ALushort offset;
 404                 offset  = fgetc(f);
 405                 offset |= fgetc(f)<<8;
 406                 if(offset != evOffset[i])
 407                 {
 408                     ERR("Unsupported evOffset[%d] value: %d (%d)\n", i, offset, evOffset[i]);
 409                     failed = AL_TRUE;
 410                 }
 411             }
 412         }
 413
 414         if(!failed)
 415         {
 416             for(i = 0;i < HRIR_COUNT;i++)
 417             {
 418                 for(j = 0;j < HRIR_LENGTH;j++)
 419                 {
 420                     ALshort coeff;
 421                     coeff  = fgetc(f);
 422                     coeff |= fgetc(f)<<8;
 423                     newdata.coeffs[i][j] = coeff;
 424                 }
 425             }
 426             for(i = 0;i < HRIR_COUNT;i++)
 427             {
 428                 ALubyte delay;
 429                 delay = fgetc(f);
 430                 newdata.delays[i] = delay;
 431                 if(delay > maxDelay)
 432                 {
 433                     ERR("Invalid delay[%d]: %d (%d)\n", i, delay, maxDelay);
 434                     failed = AL_TRUE;
 435                 }
 436             }
 437
 438             if(feof(f))
 439             {
 440                 ERR("Premature end of data\n");
 441                 failed = AL_TRUE;
 442             }
 443         }
 444
 445         fclose(f);
 446         f = NULL;
 447
 448         if(!failed)
 449         {
 450             void *temp = realloc(LoadedHrtfs, (NumLoadedHrtfs+1)*sizeof(LoadedHrtfs[0]));
 451             if(temp != NULL)
 452             {
 453                 LoadedHrtfs = temp;
 454                 TRACE("Loaded HRTF support for format: %s %uhz\n",
 455                       DevFmtChannelsString(DevFmtStereo), newdata.sampleRate);
 456                 LoadedHrtfs[NumLoadedHrtfs++] = newdata;
 457             }
 458         }
 459         else
 460             ERR("Failed to load %s\n", fname);
 461     }
 462     free(fnamelist);
 463     fnamelist = NULL;
 464 }
 465
 466 void FreeHrtf(void)
 467 {
 468     NumLoadedHrtfs = 0;
 469     free(LoadedHrtfs);
 470     LoadedHrtfs = NULL;
 471 }