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, 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     ALfloat mu[3];
 117     ALuint lidx[4], ridx[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 the normalized and attenuated HRIR coefficients using linear
 146     // interpolation when there is enough gain to warrant it.  Zero the
 147     // coefficients if gain is too low.
 148     if(gain > 0.0001f)
 149     {
 150         gain *= 1.0f/32767.0f;
 151         for(i = 0;i < HRIR_LENGTH;i++)
 152         {
 153             coeffs[i][0] = lerp(lerp(Hrtf->coeffs[lidx[0]][i], Hrtf->coeffs[lidx[1]][i], mu[0]),
 154                                 lerp(Hrtf->coeffs[lidx[2]][i], Hrtf->coeffs[lidx[3]][i], mu[1]),
 155                                 mu[2]) * gain;
 156             coeffs[i][1] = lerp(lerp(Hrtf->coeffs[ridx[0]][i], Hrtf->coeffs[ridx[1]][i], mu[0]),
 157                                 lerp(Hrtf->coeffs[ridx[2]][i], Hrtf->coeffs[ridx[3]][i], mu[1]),
 158                                 mu[2]) * gain;
 159         }
 160     }
 161     else
 162     {
 163         for(i = 0;i < HRIR_LENGTH;i++)
 164         {
 165             coeffs[i][0] = 0.0f;
 166             coeffs[i][1] = 0.0f;
 167         }
 168     }
 169
 170     // Calculate the HRIR delays using linear interpolation.
 171     delays[0] = fastf2u(lerp(lerp(Hrtf->delays[lidx[0]], Hrtf->delays[lidx[1]], mu[0]),
 172                              lerp(Hrtf->delays[lidx[2]], Hrtf->delays[lidx[3]], mu[1]),
 173                              mu[2]) * 65536.0f);
 174     delays[1] = fastf2u(lerp(lerp(Hrtf->delays[ridx[0]], Hrtf->delays[ridx[1]], mu[0]),
 175                              lerp(Hrtf->delays[ridx[2]], Hrtf->delays[ridx[3]], mu[1]),
 176                              mu[2]) * 65536.0f);
 177 }
 178
 179 // Calculates the moving HRIR target coefficients, target delays, and
 180 // stepping values for the given polar elevation and azimuth in radians.
 181 // Linear interpolation is used to increase the apparent resolution of the
 182 // HRIR dataset.  The coefficients are also normalized and attenuated by the
 183 // specified gain.  Stepping resolution and count is determined using the
 184 // given delta factor between 0.0 and 1.0.
 185 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)
 186 {
 187     ALuint evidx[2], azidx[2];
 188     ALuint lidx[4], ridx[4];
 189     ALfloat left, right;
 190     ALfloat mu[3];
 191     ALfloat step;
 192     ALuint i;
 193
 194     // Claculate elevation indices and interpolation factor.
 195     CalcEvIndices(elevation, evidx, &mu[2]);
 196
 197     // Calculate azimuth indices and interpolation factor for the first
 198     // elevation.
 199     CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
 200
 201     // Calculate the first set of linear HRIR indices for left and right
 202     // channels.
 203     lidx[0] = evOffset[evidx[0]] + azidx[0];
 204     lidx[1] = evOffset[evidx[0]] + azidx[1];
 205     ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
 206     ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
 207
 208     // Calculate azimuth indices and interpolation factor for the second
 209     // elevation.
 210     CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
 211
 212     // Calculate the second set of linear HRIR indices for left and right
 213     // channels.
 214     lidx[2] = evOffset[evidx[1]] + azidx[0];
 215     lidx[3] = evOffset[evidx[1]] + azidx[1];
 216     ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
 217     ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
 218
 219     // Calculate the stepping parameters.
 220     delta = maxf(floorf(delta*(Hrtf->sampleRate*0.015f) + 0.5f), 1.0f);
 221     step = 1.0f / delta;
 222
 223     // Calculate the normalized and attenuated target HRIR coefficients using
 224     // linear interpolation when there is enough gain to warrant it.  Zero
 225     // the target coefficients if gain is too low.  Then calculate the
 226     // coefficient stepping values using the target and previous running
 227     // coefficients.
 228     if(gain > 0.0001f)
 229     {
 230         gain *= 1.0f/32767.0f;
 231         for(i = 0;i < HRIR_LENGTH;i++)
 232         {
 233             left = coeffs[i][0] - (coeffStep[i][0] * counter);
 234             right = coeffs[i][1] - (coeffStep[i][1] * counter);
 235
 236             coeffs[i][0] = lerp(lerp(Hrtf->coeffs[lidx[0]][i], Hrtf->coeffs[lidx[1]][i], mu[0]),
 237                                 lerp(Hrtf->coeffs[lidx[2]][i], Hrtf->coeffs[lidx[3]][i], mu[1]),
 238                                 mu[2]) * gain;
 239             coeffs[i][1] = lerp(lerp(Hrtf->coeffs[ridx[0]][i], Hrtf->coeffs[ridx[1]][i], mu[0]),
 240                                 lerp(Hrtf->coeffs[ridx[2]][i], Hrtf->coeffs[ridx[3]][i], mu[1]),
 241                                 mu[2]) * gain;
 242
 243             coeffStep[i][0] = step * (coeffs[i][0] - left);
 244             coeffStep[i][1] = step * (coeffs[i][1] - right);
 245         }
 246     }
 247     else
 248     {
 249         for(i = 0;i < HRIR_LENGTH;i++)
 250         {
 251             left = coeffs[i][0] - (coeffStep[i][0] * counter);
 252             right = coeffs[i][1] - (coeffStep[i][1] * counter);
 253
 254             coeffs[i][0] = 0.0f;
 255             coeffs[i][1] = 0.0f;
 256
 257             coeffStep[i][0] = step * -left;
 258             coeffStep[i][1] = step * -right;
 259         }
 260     }
 261
 262     // Calculate the HRIR delays using linear interpolation.  Then calculate
 263     // the delay stepping values using the target and previous running
 264     // delays.
 265     left = (ALfloat)(delays[0] - (delayStep[0] * counter));
 266     right = (ALfloat)(delays[1] - (delayStep[1] * counter));
 267
 268     delays[0] = fastf2u(lerp(lerp(Hrtf->delays[lidx[0]], Hrtf->delays[lidx[1]], mu[0]),
 269                              lerp(Hrtf->delays[lidx[2]], Hrtf->delays[lidx[3]], mu[1]),
 270                              mu[2]) * 65536.0f);
 271     delays[1] = fastf2u(lerp(lerp(Hrtf->delays[ridx[0]], Hrtf->delays[ridx[1]], mu[0]),
 272                              lerp(Hrtf->delays[ridx[2]], Hrtf->delays[ridx[3]], mu[1]),
 273                              mu[2]) * 65536.0f);
 274
 275     delayStep[0] = fastf2i(step * (delays[0] - left));
 276     delayStep[1] = fastf2i(step * (delays[1] - right));
 277
 278     // The stepping count is the number of samples necessary for the HRIR to
 279     // complete its transition.  The mixer will only apply stepping for this
 280     // many samples.
 281     return fastf2u(delta);
 282 }
 283
 284 const struct Hrtf *GetHrtf(ALCdevice *device)
 285 {
 286     if(device->FmtChans == DevFmtStereo)
 287     {
 288         ALuint i;
 289         for(i = 0;i < NumLoadedHrtfs;i++)
 290         {
 291             if(device->Frequency == LoadedHrtfs[i].sampleRate)
 292                 return &LoadedHrtfs[i];
 293         }
 294         if(device->Frequency == DefaultHrtf.sampleRate)
 295             return &DefaultHrtf;
 296     }
 297     ERR("Incompatible format: %s %uhz\n",
 298         DevFmtChannelsString(device->FmtChans), device->Frequency);
 299     return NULL;
 300 }
 301
 302 void InitHrtf(void)
 303 {
 304     char *fnamelist=NULL, *next=NULL;
 305     const char *val;
 306
 307     if(ConfigValueStr(NULL, "hrtf_tables", &val))
 308         next = fnamelist = strdup(val);
 309     while(next && *next)
 310     {
 311         const ALubyte maxDelay = SRC_HISTORY_LENGTH-1;
 312         struct Hrtf newdata;
 313         ALboolean failed;
 314         ALchar magic[9];
 315         ALsizei i, j;
 316         char *fname;
 317         FILE *f;
 318
 319         fname = next;
 320         next = strchr(fname, ',');
 321         if(next)
 322         {
 323             while(next != fname)
 324             {
 325                 next--;
 326                 if(!isspace(*next))
 327                 {
 328                     *(next++) = '\0';
 329                     break;
 330                 }
 331             }
 332             while(isspace(*next) || *next == ',')
 333                 next++;
 334         }
 335
 336         if(!fname[0])
 337             continue;
 338         TRACE("Loading %s\n", fname);
 339         f = fopen(fname, "rb");
 340         if(f == NULL)
 341         {
 342             ERR("Could not open %s\n", fname);
 343             continue;
 344         }
 345
 346         failed = AL_FALSE;
 347         if(fread(magic, 1, sizeof(magicMarker), f) != sizeof(magicMarker))
 348         {
 349             ERR("Failed to read magic marker\n");
 350             failed = AL_TRUE;
 351         }
 352         else if(memcmp(magic, magicMarker, sizeof(magicMarker)) != 0)
 353         {
 354             magic[8] = 0;
 355             ERR("Invalid magic marker: \"%s\"\n", magic);
 356             failed = AL_TRUE;
 357         }
 358
 359         if(!failed)
 360         {
 361             ALushort hrirCount, hrirSize;
 362             ALubyte  evCount;
 363
 364             newdata.sampleRate  = fgetc(f);
 365             newdata.sampleRate |= fgetc(f)<<8;
 366             newdata.sampleRate |= fgetc(f)<<16;
 367             newdata.sampleRate |= fgetc(f)<<24;
 368
 369             hrirCount  = fgetc(f);
 370             hrirCount |= fgetc(f)<<8;
 371
 372             hrirSize  = fgetc(f);
 373             hrirSize |= fgetc(f)<<8;
 374
 375             evCount = fgetc(f);
 376
 377             if(hrirCount != HRIR_COUNT || hrirSize != HRIR_LENGTH || evCount != ELEV_COUNT)
 378             {
 379                 ERR("Unsupported value: hrirCount=%d (%d), hrirSize=%d (%d), evCount=%d (%d)\n",
 380                     hrirCount, HRIR_COUNT, hrirSize, HRIR_LENGTH, evCount, ELEV_COUNT);
 381                 failed = AL_TRUE;
 382             }
 383         }
 384
 385         if(!failed)
 386         {
 387             for(i = 0;i < ELEV_COUNT;i++)
 388             {
 389                 ALushort offset;
 390                 offset  = fgetc(f);
 391                 offset |= fgetc(f)<<8;
 392                 if(offset != evOffset[i])
 393                 {
 394                     ERR("Unsupported evOffset[%d] value: %d (%d)\n", i, offset, evOffset[i]);
 395                     failed = AL_TRUE;
 396                 }
 397             }
 398         }
 399
 400         if(!failed)
 401         {
 402             for(i = 0;i < HRIR_COUNT;i++)
 403             {
 404                 for(j = 0;j < HRIR_LENGTH;j++)
 405                 {
 406                     ALshort coeff;
 407                     coeff  = fgetc(f);
 408                     coeff |= fgetc(f)<<8;
 409                     newdata.coeffs[i][j] = coeff;
 410                 }
 411             }
 412             for(i = 0;i < HRIR_COUNT;i++)
 413             {
 414                 ALubyte delay;
 415                 delay = fgetc(f);
 416                 newdata.delays[i] = delay;
 417                 if(delay > maxDelay)
 418                 {
 419                     ERR("Invalid delay[%d]: %d (%d)\n", i, delay, maxDelay);
 420                     failed = AL_TRUE;
 421                 }
 422             }
 423
 424             if(feof(f))
 425             {
 426                 ERR("Premature end of data\n");
 427                 failed = AL_TRUE;
 428             }
 429         }
 430
 431         fclose(f);
 432         f = NULL;
 433
 434         if(!failed)
 435         {
 436             void *temp = realloc(LoadedHrtfs, (NumLoadedHrtfs+1)*sizeof(LoadedHrtfs[0]));
 437             if(temp != NULL)
 438             {
 439                 LoadedHrtfs = temp;
 440                 TRACE("Loaded HRTF support for format: %s %uhz\n",
 441                       DevFmtChannelsString(DevFmtStereo), newdata.sampleRate);
 442                 LoadedHrtfs[NumLoadedHrtfs++] = newdata;
 443             }
 444         }
 445         else
 446             ERR("Failed to load %s\n", fname);
 447     }
 448     free(fnamelist);
 449     fnamelist = NULL;
 450 }
 451
 452 void FreeHrtf(void)
 453 {
 454     NumLoadedHrtfs = 0;
 455     free(LoadedHrtfs);
 456     LoadedHrtfs = NULL;
 457 }