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