rlm@0: /**
rlm@0:  * OpenAL cross platform audio library
rlm@0:  * Copyright (C) 2011 by Chris Robinson
rlm@0:  * This library is free software; you can redistribute it and/or
rlm@0:  *  modify it under the terms of the GNU Library General Public
rlm@0:  *  License as published by the Free Software Foundation; either
rlm@0:  *  version 2 of the License, or (at your option) any later version.
rlm@0:  *
rlm@0:  * This library is distributed in the hope that it will be useful,
rlm@0:  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
rlm@0:  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
rlm@0:  *  Library General Public License for more details.
rlm@0:  *
rlm@0:  * You should have received a copy of the GNU Library General Public
rlm@0:  *  License along with this library; if not, write to the
rlm@0:  *  Free Software Foundation, Inc., 59 Temple Place - Suite 330,
rlm@0:  *  Boston, MA  02111-1307, USA.
rlm@0:  * Or go to http://www.gnu.org/copyleft/lgpl.html
rlm@0:  */
rlm@0: 
rlm@0: #include "config.h"
rlm@0: 
rlm@0: #include "AL/al.h"
rlm@0: #include "AL/alc.h"
rlm@0: #include "alMain.h"
rlm@0: #include "alSource.h"
rlm@0: 
rlm@0: /* External HRTF file format (LE byte order):
rlm@0:  *
rlm@0:  * ALchar   magic[8] = "MinPHR00";
rlm@0:  * ALuint   sampleRate;
rlm@0:  *
rlm@0:  * ALushort hrirCount; // Required value: 828
rlm@0:  * ALushort hrirSize;  // Required value: 32
rlm@0:  * ALubyte  evCount;   // Required value: 19
rlm@0:  *
rlm@0:  * ALushort evOffset[evCount]; // Required values:
rlm@0:  *   { 0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827 }
rlm@0:  *
rlm@0:  * ALushort coefficients[hrirCount][hrirSize];
rlm@0:  * ALubyte  delays[hrirCount]; // Element values must not exceed 127
rlm@0:  */
rlm@0: 
rlm@0: static const ALchar magicMarker[8] = "MinPHR00";
rlm@0: 
rlm@0: #define HRIR_COUNT 828
rlm@0: #define ELEV_COUNT 19
rlm@0: 
rlm@0: 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 };
rlm@0: 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 };
rlm@0: 
rlm@0: static struct Hrtf {
rlm@0:     ALuint sampleRate;
rlm@0:     ALshort coeffs[HRIR_COUNT][HRIR_LENGTH];
rlm@0:     ALubyte delays[HRIR_COUNT];
rlm@0: } Hrtf = {
rlm@0:     44100,
rlm@0: #include "hrtf_tables.inc"
rlm@0: };
rlm@0: 
rlm@0: // Calculate the elevation indices given the polar elevation in radians.
rlm@0: // This will return two indices between 0 and (ELEV_COUNT-1) and an
rlm@0: // interpolation factor between 0.0 and 1.0.
rlm@0: static void CalcEvIndices(ALfloat ev, ALuint *evidx, ALfloat *evmu)
rlm@0: {
rlm@0:     ev = (M_PI/2.0f + ev) * (ELEV_COUNT-1) / M_PI;
rlm@0:     evidx[0] = (ALuint)ev;
rlm@0:     evidx[1] = minu(evidx[0] + 1, ELEV_COUNT-1);
rlm@0:     *evmu = ev - evidx[0];
rlm@0: }
rlm@0: 
rlm@0: // Calculate the azimuth indices given the polar azimuth in radians.  This
rlm@0: // will return two indices between 0 and (azCount [ei] - 1) and an
rlm@0: // interpolation factor between 0.0 and 1.0.
rlm@0: static void CalcAzIndices(ALuint evidx, ALfloat az, ALuint *azidx, ALfloat *azmu)
rlm@0: {
rlm@0:     az = (M_PI*2.0f + az) * azCount[evidx] / (M_PI*2.0f);
rlm@0:     azidx[0] = (ALuint)az % azCount[evidx];
rlm@0:     azidx[1] = (azidx[0] + 1) % azCount[evidx];
rlm@0:     *azmu = az - floor(az);
rlm@0: }
rlm@0: 
rlm@0: // Calculates the normalized HRTF transition factor (delta) from the changes
rlm@0: // in gain and listener to source angle between updates.  The result is a
rlm@0: // normalized delta factor than can be used to calculate moving HRIR stepping
rlm@0: // values.
rlm@0: ALfloat CalcHrtfDelta(ALfloat oldGain, ALfloat newGain, const ALfloat olddir[3], const ALfloat newdir[3])
rlm@0: {
rlm@0:     ALfloat gainChange, angleChange;
rlm@0: 
rlm@0:     // Calculate the normalized dB gain change.
rlm@0:     newGain = maxf(newGain, 0.0001f);
rlm@0:     oldGain = maxf(oldGain, 0.0001f);
rlm@0:     gainChange = aluFabs(log10(newGain / oldGain) / log10(0.0001f));
rlm@0: 
rlm@0:     // Calculate the normalized listener to source angle change when there is
rlm@0:     // enough gain to notice it.
rlm@0:     angleChange = 0.0f;
rlm@0:     if(gainChange > 0.0001f || newGain > 0.0001f)
rlm@0:     {
rlm@0:         // No angle change when the directions are equal or degenerate (when
rlm@0:         // both have zero length).
rlm@0:         if(newdir[0]-olddir[0] || newdir[1]-olddir[1] || newdir[2]-olddir[2])
rlm@0:             angleChange = aluAcos(olddir[0]*newdir[0] +
rlm@0:                                   olddir[1]*newdir[1] +
rlm@0:                                   olddir[2]*newdir[2]) / M_PI;
rlm@0: 
rlm@0:     }
rlm@0: 
rlm@0:     // Use the largest of the two changes for the delta factor, and apply a
rlm@0:     // significance shaping function to it.
rlm@0:     return clampf(angleChange*2.0f, gainChange*2.0f, 1.0f);
rlm@0: }
rlm@0: 
rlm@0: // Calculates static HRIR coefficients and delays for the given polar
rlm@0: // elevation and azimuth in radians.  Linear interpolation is used to
rlm@0: // increase the apparent resolution of the HRIR dataset.  The coefficients
rlm@0: // are also normalized and attenuated by the specified gain.
rlm@0: void GetLerpedHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat (*coeffs)[2], ALuint *delays)
rlm@0: {
rlm@0:     ALuint evidx[2], azidx[2];
rlm@0:     ALfloat mu[3];
rlm@0:     ALuint lidx[4], ridx[4];
rlm@0:     ALuint i;
rlm@0: 
rlm@0:     // Claculate elevation indices and interpolation factor.
rlm@0:     CalcEvIndices(elevation, evidx, &mu[2]);
rlm@0: 
rlm@0:     // Calculate azimuth indices and interpolation factor for the first
rlm@0:     // elevation.
rlm@0:     CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
rlm@0: 
rlm@0:     // Calculate the first set of linear HRIR indices for left and right
rlm@0:     // channels.
rlm@0:     lidx[0] = evOffset[evidx[0]] + azidx[0];
rlm@0:     lidx[1] = evOffset[evidx[0]] + azidx[1];
rlm@0:     ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
rlm@0:     ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
rlm@0: 
rlm@0:     // Calculate azimuth indices and interpolation factor for the second
rlm@0:     // elevation.
rlm@0:     CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
rlm@0: 
rlm@0:     // Calculate the second set of linear HRIR indices for left and right
rlm@0:     // channels.
rlm@0:     lidx[2] = evOffset[evidx[1]] + azidx[0];
rlm@0:     lidx[3] = evOffset[evidx[1]] + azidx[1];
rlm@0:     ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
rlm@0:     ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
rlm@0: 
rlm@0:     // Calculate the normalized and attenuated HRIR coefficients using linear
rlm@0:     // interpolation when there is enough gain to warrant it.  Zero the
rlm@0:     // coefficients if gain is too low.
rlm@0:     if(gain > 0.0001f)
rlm@0:     {
rlm@0:         ALdouble scale = gain * (1.0/32767.0);
rlm@0:         for(i = 0;i < HRIR_LENGTH;i++)
rlm@0:         {
rlm@0:             coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),
rlm@0:                                 lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),
rlm@0:                                 mu[2]) * scale;
rlm@0:             coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),
rlm@0:                                 lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),
rlm@0:                                 mu[2]) * scale;
rlm@0:         }
rlm@0:     }
rlm@0:     else
rlm@0:     {
rlm@0:         for(i = 0;i < HRIR_LENGTH;i++)
rlm@0:         {
rlm@0:             coeffs[i][0] = 0.0f;
rlm@0:             coeffs[i][1] = 0.0f;
rlm@0:         }
rlm@0:     }
rlm@0: 
rlm@0:     // Calculate the HRIR delays using linear interpolation.
rlm@0:     delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),
rlm@0:                               lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),
rlm@0:                               mu[2]) * 65536.0f);
rlm@0:     delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),
rlm@0:                               lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),
rlm@0:                               mu[2]) * 65536.0f);
rlm@0: }
rlm@0: 
rlm@0: // Calculates the moving HRIR target coefficients, target delays, and
rlm@0: // stepping values for the given polar elevation and azimuth in radians.
rlm@0: // Linear interpolation is used to increase the apparent resolution of the
rlm@0: // HRIR dataset.  The coefficients are also normalized and attenuated by the
rlm@0: // specified gain.  Stepping resolution and count is determined using the
rlm@0: // given delta factor between 0.0 and 1.0.
rlm@0: ALuint GetMovingHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat delta, ALint counter, ALfloat (*coeffs)[2], ALuint *delays, ALfloat (*coeffStep)[2], ALint *delayStep)
rlm@0: {
rlm@0:     ALuint evidx[2], azidx[2];
rlm@0:     ALuint lidx[4], ridx[4];
rlm@0:     ALfloat left, right;
rlm@0:     ALfloat mu[3];
rlm@0:     ALfloat step;
rlm@0:     ALuint i;
rlm@0: 
rlm@0:     // Claculate elevation indices and interpolation factor.
rlm@0:     CalcEvIndices(elevation, evidx, &mu[2]);
rlm@0: 
rlm@0:     // Calculate azimuth indices and interpolation factor for the first
rlm@0:     // elevation.
rlm@0:     CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
rlm@0: 
rlm@0:     // Calculate the first set of linear HRIR indices for left and right
rlm@0:     // channels.
rlm@0:     lidx[0] = evOffset[evidx[0]] + azidx[0];
rlm@0:     lidx[1] = evOffset[evidx[0]] + azidx[1];
rlm@0:     ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
rlm@0:     ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
rlm@0: 
rlm@0:     // Calculate azimuth indices and interpolation factor for the second
rlm@0:     // elevation.
rlm@0:     CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
rlm@0: 
rlm@0:     // Calculate the second set of linear HRIR indices for left and right
rlm@0:     // channels.
rlm@0:     lidx[2] = evOffset[evidx[1]] + azidx[0];
rlm@0:     lidx[3] = evOffset[evidx[1]] + azidx[1];
rlm@0:     ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
rlm@0:     ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
rlm@0: 
rlm@0:     // Calculate the stepping parameters.
rlm@0:     delta = maxf(floor(delta*(Hrtf.sampleRate*0.015f) + 0.5), 1.0f);
rlm@0:     step = 1.0f / delta;
rlm@0: 
rlm@0:     // Calculate the normalized and attenuated target HRIR coefficients using
rlm@0:     // linear interpolation when there is enough gain to warrant it.  Zero
rlm@0:     // the target coefficients if gain is too low.  Then calculate the
rlm@0:     // coefficient stepping values using the target and previous running
rlm@0:     // coefficients.
rlm@0:     if(gain > 0.0001f)
rlm@0:     {
rlm@0:         ALdouble scale = gain * (1.0/32767.0);
rlm@0:         for(i = 0;i < HRIR_LENGTH;i++)
rlm@0:         {
rlm@0:             left = coeffs[i][0] - (coeffStep[i][0] * counter);
rlm@0:             right = coeffs[i][1] - (coeffStep[i][1] * counter);
rlm@0: 
rlm@0:             coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),
rlm@0:                                 lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),
rlm@0:                                 mu[2]) * scale;
rlm@0:             coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),
rlm@0:                                 lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),
rlm@0:                                 mu[2]) * scale;
rlm@0: 
rlm@0:             coeffStep[i][0] = step * (coeffs[i][0] - left);
rlm@0:             coeffStep[i][1] = step * (coeffs[i][1] - right);
rlm@0:         }
rlm@0:     }
rlm@0:     else
rlm@0:     {
rlm@0:         for(i = 0;i < HRIR_LENGTH;i++)
rlm@0:         {
rlm@0:             left = coeffs[i][0] - (coeffStep[i][0] * counter);
rlm@0:             right = coeffs[i][1] - (coeffStep[i][1] * counter);
rlm@0: 
rlm@0:             coeffs[i][0] = 0.0f;
rlm@0:             coeffs[i][1] = 0.0f;
rlm@0: 
rlm@0:             coeffStep[i][0] = step * -left;
rlm@0:             coeffStep[i][1] = step * -right;
rlm@0:         }
rlm@0:     }
rlm@0: 
rlm@0:     // Calculate the HRIR delays using linear interpolation.  Then calculate
rlm@0:     // the delay stepping values using the target and previous running
rlm@0:     // delays.
rlm@0:     left = delays[0] - (delayStep[0] * counter);
rlm@0:     right = delays[1] - (delayStep[1] * counter);
rlm@0: 
rlm@0:     delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),
rlm@0:                               lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),
rlm@0:                               mu[2]) * 65536.0f);
rlm@0:     delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),
rlm@0:                               lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),
rlm@0:                               mu[2]) * 65536.0f);
rlm@0: 
rlm@0:     delayStep[0] = (ALint)(step * (delays[0] - left));
rlm@0:     delayStep[1] = (ALint)(step * (delays[1] - right));
rlm@0: 
rlm@0:     // The stepping count is the number of samples necessary for the HRIR to
rlm@0:     // complete its transition.  The mixer will only apply stepping for this
rlm@0:     // many samples.
rlm@0:     return (ALuint)delta;
rlm@0: }
rlm@0: 
rlm@0: ALCboolean IsHrtfCompatible(ALCdevice *device)
rlm@0: {
rlm@0:     if(device->FmtChans == DevFmtStereo && device->Frequency == Hrtf.sampleRate)
rlm@0:         return ALC_TRUE;
rlm@0:     ERR("Incompatible HRTF format: %s %uhz (%s %uhz needed)\n",
rlm@0:         DevFmtChannelsString(device->FmtChans), device->Frequency,
rlm@0:         DevFmtChannelsString(DevFmtStereo), Hrtf.sampleRate);
rlm@0:     return ALC_FALSE;
rlm@0: }
rlm@0: 
rlm@0: void InitHrtf(void)
rlm@0: {
rlm@0:     const char *fname;
rlm@0:     FILE *f = NULL;
rlm@0: 
rlm@0:     fname = GetConfigValue(NULL, "hrtf_tables", "");
rlm@0:     if(fname[0] != '\0')
rlm@0:     {
rlm@0:         f = fopen(fname, "rb");
rlm@0:         if(f == NULL)
rlm@0:             ERR("Could not open %s\n", fname);
rlm@0:     }
rlm@0:     if(f != NULL)
rlm@0:     {
rlm@0:         const ALubyte maxDelay = SRC_HISTORY_LENGTH-1;
rlm@0:         ALboolean failed = AL_FALSE;
rlm@0:         struct Hrtf newdata;
rlm@0:         ALchar magic[9];
rlm@0:         ALsizei i, j;
rlm@0: 
rlm@0:         if(fread(magic, 1, sizeof(magicMarker), f) != sizeof(magicMarker))
rlm@0:         {
rlm@0:             ERR("Failed to read magic marker\n");
rlm@0:             failed = AL_TRUE;
rlm@0:         }
rlm@0:         else if(memcmp(magic, magicMarker, sizeof(magicMarker)) != 0)
rlm@0:         {
rlm@0:             magic[8] = 0;
rlm@0:             ERR("Invalid magic marker: \"%s\"\n", magic);
rlm@0:             failed = AL_TRUE;
rlm@0:         }
rlm@0: 
rlm@0:         if(!failed)
rlm@0:         {
rlm@0:             ALushort hrirCount, hrirSize;
rlm@0:             ALubyte  evCount;
rlm@0: 
rlm@0:             newdata.sampleRate  = fgetc(f);
rlm@0:             newdata.sampleRate |= fgetc(f)<<8;
rlm@0:             newdata.sampleRate |= fgetc(f)<<16;
rlm@0:             newdata.sampleRate |= fgetc(f)<<24;
rlm@0: 
rlm@0:             hrirCount  = fgetc(f);
rlm@0:             hrirCount |= fgetc(f)<<8;
rlm@0: 
rlm@0:             hrirSize  = fgetc(f);
rlm@0:             hrirSize |= fgetc(f)<<8;
rlm@0: 
rlm@0:             evCount = fgetc(f);
rlm@0: 
rlm@0:             if(hrirCount != HRIR_COUNT || hrirSize != HRIR_LENGTH || evCount != ELEV_COUNT)
rlm@0:             {
rlm@0:                 ERR("Unsupported value: hrirCount=%d (%d), hrirSize=%d (%d), evCount=%d (%d)\n",
rlm@0:                     hrirCount, HRIR_COUNT, hrirSize, HRIR_LENGTH, evCount, ELEV_COUNT);
rlm@0:                 failed = AL_TRUE;
rlm@0:             }
rlm@0:         }
rlm@0: 
rlm@0:         if(!failed)
rlm@0:         {
rlm@0:             for(i = 0;i < HRIR_COUNT;i++)
rlm@0:             {
rlm@0:                 ALushort offset;
rlm@0:                 offset  = fgetc(f);
rlm@0:                 offset |= fgetc(f)<<8;
rlm@0:                 if(offset != evOffset[i])
rlm@0:                 {
rlm@0:                     ERR("Unsupported evOffset[%d] value: %d (%d)\n", i, offset, evOffset[i]);
rlm@0:                     failed = AL_TRUE;
rlm@0:                 }
rlm@0:             }
rlm@0:         }
rlm@0: 
rlm@0:         if(!failed)
rlm@0:         {
rlm@0:             for(i = 0;i < HRIR_COUNT;i++)
rlm@0:             {
rlm@0:                 for(j = 0;j < HRIR_LENGTH;j++)
rlm@0:                 {
rlm@0:                     ALshort coeff;
rlm@0:                     coeff  = fgetc(f);
rlm@0:                     coeff |= fgetc(f)<<8;
rlm@0:                     newdata.coeffs[i][j] = coeff;
rlm@0:                 }
rlm@0:             }
rlm@0:             for(i = 0;i < HRIR_COUNT;i++)
rlm@0:             {
rlm@0:                 ALubyte delay;
rlm@0:                 delay = fgetc(f);
rlm@0:                 newdata.delays[i] = delay;
rlm@0:                 if(delay > maxDelay)
rlm@0:                 {
rlm@0:                     ERR("Invalid delay[%d]: %d (%d)\n", i, delay, maxDelay);
rlm@0:                     failed = AL_TRUE;
rlm@0:                 }
rlm@0:             }
rlm@0: 
rlm@0:             if(feof(f))
rlm@0:             {
rlm@0:                 ERR("Premature end of data\n");
rlm@0:                 failed = AL_TRUE;
rlm@0:             }
rlm@0:         }
rlm@0: 
rlm@0:         fclose(f);
rlm@0:         f = NULL;
rlm@0: 
rlm@0:         if(!failed)
rlm@0:             Hrtf = newdata;
rlm@0:         else
rlm@0:             ERR("Failed to load %s\n", fname);
rlm@0:     }
rlm@0: }