Mercurial > audio-send

     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 "AL/al.h"

    24 #include "AL/alc.h"

    25 #include "alMain.h"

    26 #include "alSource.h"

    27 

    28 /* External HRTF file format (LE byte order):

    29  *

    30  * ALchar   magic[8] = "MinPHR00";

    31  * ALuint   sampleRate;

    32  *

    33  * ALushort hrirCount; // Required value: 828

    34  * ALushort hrirSize;  // Required value: 32

    35  * ALubyte  evCount;   // Required value: 19

    36  *

    37  * ALushort evOffset[evCount]; // Required values:

    38  *   { 0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827 }

    39  *

    40  * ALushort coefficients[hrirCount][hrirSize];

    41  * ALubyte  delays[hrirCount]; // Element values must not exceed 127

    42  */

    43 

    44 static const ALchar magicMarker[8] = "MinPHR00";

    45 

    46 #define HRIR_COUNT 828

    47 #define ELEV_COUNT 19

    48 

    49 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 };

    50 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 };

    51 

    52 static struct Hrtf {

    53     ALuint sampleRate;

    54     ALshort coeffs[HRIR_COUNT][HRIR_LENGTH];

    55     ALubyte delays[HRIR_COUNT];

    56 } Hrtf = {

    57     44100,

    58 #include "hrtf_tables.inc"

    59 };

    60 

    61 // Calculate the elevation indices given the polar elevation in radians.

    62 // This will return two indices between 0 and (ELEV_COUNT-1) and an

    63 // interpolation factor between 0.0 and 1.0.

    64 static void CalcEvIndices(ALfloat ev, ALuint *evidx, ALfloat *evmu)

    65 {

    66     ev = (M_PI/2.0f + ev) * (ELEV_COUNT-1) / M_PI;

    67     evidx[0] = (ALuint)ev;

    68     evidx[1] = minu(evidx[0] + 1, ELEV_COUNT-1);

    69     *evmu = ev - evidx[0];

    70 }

    71 

    72 // Calculate the azimuth indices given the polar azimuth in radians.  This

    73 // will return two indices between 0 and (azCount [ei] - 1) and an

    74 // interpolation factor between 0.0 and 1.0.

    75 static void CalcAzIndices(ALuint evidx, ALfloat az, ALuint *azidx, ALfloat *azmu)

    76 {

    77     az = (M_PI*2.0f + az) * azCount[evidx] / (M_PI*2.0f);

    78     azidx[0] = (ALuint)az % azCount[evidx];

    79     azidx[1] = (azidx[0] + 1) % azCount[evidx];

    80     *azmu = az - floor(az);

    81 }

    82 

    83 // Calculates the normalized HRTF transition factor (delta) from the changes

    84 // in gain and listener to source angle between updates.  The result is a

    85 // normalized delta factor than can be used to calculate moving HRIR stepping

    86 // values.

    87 ALfloat CalcHrtfDelta(ALfloat oldGain, ALfloat newGain, const ALfloat olddir[3], const ALfloat newdir[3])

    88 {

    89     ALfloat gainChange, angleChange;

    90 

    91     // Calculate the normalized dB gain change.

    92     newGain = maxf(newGain, 0.0001f);

    93     oldGain = maxf(oldGain, 0.0001f);

    94     gainChange = aluFabs(log10(newGain / oldGain) / log10(0.0001f));

    95 

    96     // Calculate the normalized listener to source angle change when there is

    97     // enough gain to notice it.

    98     angleChange = 0.0f;

    99     if(gainChange > 0.0001f || newGain > 0.0001f)

   100     {

   101         // No angle change when the directions are equal or degenerate (when

   102         // both have zero length).

   103         if(newdir[0]-olddir[0] || newdir[1]-olddir[1] || newdir[2]-olddir[2])

   104             angleChange = aluAcos(olddir[0]*newdir[0] +

   105                                   olddir[1]*newdir[1] +

   106                                   olddir[2]*newdir[2]) / M_PI;

   107 

   108     }

   109 

   110     // Use the largest of the two changes for the delta factor, and apply a

   111     // significance shaping function to it.

   112     return clampf(angleChange*2.0f, gainChange*2.0f, 1.0f);

   113 }

   114 

   115 // Calculates static HRIR coefficients and delays for the given polar

   116 // elevation and azimuth in radians.  Linear interpolation is used to

   117 // increase the apparent resolution of the HRIR dataset.  The coefficients

   118 // are also normalized and attenuated by the specified gain.

   119 void GetLerpedHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat (*coeffs)[2], ALuint *delays)

   120 {

   121     ALuint evidx[2], azidx[2];

   122     ALfloat mu[3];

   123     ALuint lidx[4], ridx[4];

   124     ALuint i;

   125 

   126     // Claculate elevation indices and interpolation factor.

   127     CalcEvIndices(elevation, evidx, &mu[2]);

   128 

   129     // Calculate azimuth indices and interpolation factor for the first

   130     // elevation.

   131     CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);

   132 

   133     // Calculate the first set of linear HRIR indices for left and right

   134     // channels.

   135     lidx[0] = evOffset[evidx[0]] + azidx[0];

   136     lidx[1] = evOffset[evidx[0]] + azidx[1];

   137     ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);

   138     ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);

   139 

   140     // Calculate azimuth indices and interpolation factor for the second

   141     // elevation.

   142     CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);

   143 

   144     // Calculate the second set of linear HRIR indices for left and right

   145     // channels.

   146     lidx[2] = evOffset[evidx[1]] + azidx[0];

   147     lidx[3] = evOffset[evidx[1]] + azidx[1];

   148     ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);

   149     ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);

   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         ALdouble scale = gain * (1.0/32767.0);

   157         for(i = 0;i < HRIR_LENGTH;i++)

   158         {

   159             coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),

   160                                 lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),

   161                                 mu[2]) * scale;

   162             coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),

   163                                 lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),

   164                                 mu[2]) * scale;

   165         }

   166     }

   167     else

   168     {

   169         for(i = 0;i < HRIR_LENGTH;i++)

   170         {

   171             coeffs[i][0] = 0.0f;

   172             coeffs[i][1] = 0.0f;

   173         }

   174     }

   175 

   176     // Calculate the HRIR delays using linear interpolation.

   177     delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),

   178                               lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),

   179                               mu[2]) * 65536.0f);

   180     delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),

   181                               lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),

   182                               mu[2]) * 65536.0f);

   183 }

   184 

   185 // Calculates the moving HRIR target coefficients, target delays, and

   186 // stepping values for the given polar elevation and azimuth in radians.

   187 // Linear interpolation is used to increase the apparent resolution of the

   188 // HRIR dataset.  The coefficients are also normalized and attenuated by the

   189 // specified gain.  Stepping resolution and count is determined using the

   190 // given delta factor between 0.0 and 1.0.

   191 ALuint GetMovingHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat delta, ALint counter, ALfloat (*coeffs)[2], ALuint *delays, ALfloat (*coeffStep)[2], ALint *delayStep)

   192 {

   193     ALuint evidx[2], azidx[2];

   194     ALuint lidx[4], ridx[4];

   195     ALfloat left, right;

   196     ALfloat mu[3];

   197     ALfloat step;

   198     ALuint i;

   199 

   200     // Claculate elevation indices and interpolation factor.

   201     CalcEvIndices(elevation, evidx, &mu[2]);

   202 

   203     // Calculate azimuth indices and interpolation factor for the first

   204     // elevation.

   205     CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);

   206 

   207     // Calculate the first set of linear HRIR indices for left and right

   208     // channels.

   209     lidx[0] = evOffset[evidx[0]] + azidx[0];

   210     lidx[1] = evOffset[evidx[0]] + azidx[1];

   211     ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);

   212     ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);

   213 

   214     // Calculate azimuth indices and interpolation factor for the second

   215     // elevation.

   216     CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);

   217 

   218     // Calculate the second set of linear HRIR indices for left and right

   219     // channels.

   220     lidx[2] = evOffset[evidx[1]] + azidx[0];

   221     lidx[3] = evOffset[evidx[1]] + azidx[1];

   222     ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);

   223     ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);

   224 

   225     // Calculate the stepping parameters.

   226     delta = maxf(floor(delta*(Hrtf.sampleRate*0.015f) + 0.5), 1.0f);

   227     step = 1.0f / delta;

   228 

   229     // Calculate the normalized and attenuated target HRIR coefficients using

   230     // linear interpolation when there is enough gain to warrant it.  Zero

   231     // the target coefficients if gain is too low.  Then calculate the

   232     // coefficient stepping values using the target and previous running

   233     // coefficients.

   234     if(gain > 0.0001f)

   235     {

   236         ALdouble scale = gain * (1.0/32767.0);

   237         for(i = 0;i < HRIR_LENGTH;i++)

   238         {

   239             left = coeffs[i][0] - (coeffStep[i][0] * counter);

   240             right = coeffs[i][1] - (coeffStep[i][1] * counter);

   241 

   242             coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),

   243                                 lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),

   244                                 mu[2]) * scale;

   245             coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),

   246                                 lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),

   247                                 mu[2]) * scale;

   248 

   249             coeffStep[i][0] = step * (coeffs[i][0] - left);

   250             coeffStep[i][1] = step * (coeffs[i][1] - right);

   251         }

   252     }

   253     else

   254     {

   255         for(i = 0;i < HRIR_LENGTH;i++)

   256         {

   257             left = coeffs[i][0] - (coeffStep[i][0] * counter);

   258             right = coeffs[i][1] - (coeffStep[i][1] * counter);

   259 

   260             coeffs[i][0] = 0.0f;

   261             coeffs[i][1] = 0.0f;

   262 

   263             coeffStep[i][0] = step * -left;

   264             coeffStep[i][1] = step * -right;

   265         }

   266     }

   267 

   268     // Calculate the HRIR delays using linear interpolation.  Then calculate

   269     // the delay stepping values using the target and previous running

   270     // delays.

   271     left = delays[0] - (delayStep[0] * counter);

   272     right = delays[1] - (delayStep[1] * counter);

   273 

   274     delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),

   275                               lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),

   276                               mu[2]) * 65536.0f);

   277     delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),

   278                               lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),

   279                               mu[2]) * 65536.0f);

   280 

   281     delayStep[0] = (ALint)(step * (delays[0] - left));

   282     delayStep[1] = (ALint)(step * (delays[1] - right));

   283 

   284     // The stepping count is the number of samples necessary for the HRIR to

   285     // complete its transition.  The mixer will only apply stepping for this

   286     // many samples.

   287     return (ALuint)delta;

   288 }

   289 

   290 ALCboolean IsHrtfCompatible(ALCdevice *device)

   291 {

   292     if(device->FmtChans == DevFmtStereo && device->Frequency == Hrtf.sampleRate)

   293         return ALC_TRUE;

   294     ERR("Incompatible HRTF format: %s %uhz (%s %uhz needed)\n",

   295         DevFmtChannelsString(device->FmtChans), device->Frequency,

   296         DevFmtChannelsString(DevFmtStereo), Hrtf.sampleRate);

   297     return ALC_FALSE;

   298 }

   299 

   300 void InitHrtf(void)

   301 {

   302     const char *fname;

   303     FILE *f = NULL;

   304 

   305     fname = GetConfigValue(NULL, "hrtf_tables", "");

   306     if(fname[0] != '\0')

   307     {

   308         f = fopen(fname, "rb");

   309         if(f == NULL)

   310             ERR("Could not open %s\n", fname);

   311     }

   312     if(f != NULL)

   313     {

   314         const ALubyte maxDelay = SRC_HISTORY_LENGTH-1;

   315         ALboolean failed = AL_FALSE;

   316         struct Hrtf newdata;

   317         ALchar magic[9];

   318         ALsizei i, j;

   319 

   320         if(fread(magic, 1, sizeof(magicMarker), f) != sizeof(magicMarker))

   321         {

   322             ERR("Failed to read magic marker\n");

   323             failed = AL_TRUE;

   324         }

   325         else if(memcmp(magic, magicMarker, sizeof(magicMarker)) != 0)

   326         {

   327             magic[8] = 0;

   328             ERR("Invalid magic marker: \"%s\"\n", magic);

   329             failed = AL_TRUE;

   330         }

   331 

   332         if(!failed)

   333         {

   334             ALushort hrirCount, hrirSize;

   335             ALubyte  evCount;

   336 

   337             newdata.sampleRate  = fgetc(f);

   338             newdata.sampleRate |= fgetc(f)<<8;

   339             newdata.sampleRate |= fgetc(f)<<16;

   340             newdata.sampleRate |= fgetc(f)<<24;

   341 

   342             hrirCount  = fgetc(f);

   343             hrirCount |= fgetc(f)<<8;

   344 

   345             hrirSize  = fgetc(f);

   346             hrirSize |= fgetc(f)<<8;

   347 

   348             evCount = fgetc(f);

   349 

   350             if(hrirCount != HRIR_COUNT || hrirSize != HRIR_LENGTH || evCount != ELEV_COUNT)

   351             {

   352                 ERR("Unsupported value: hrirCount=%d (%d), hrirSize=%d (%d), evCount=%d (%d)\n",

   353                     hrirCount, HRIR_COUNT, hrirSize, HRIR_LENGTH, evCount, ELEV_COUNT);

   354                 failed = AL_TRUE;

   355             }

   356         }

   357 

   358         if(!failed)

   359         {

   360             for(i = 0;i < HRIR_COUNT;i++)

   361             {

   362                 ALushort offset;

   363                 offset  = fgetc(f);

   364                 offset |= fgetc(f)<<8;

   365                 if(offset != evOffset[i])

   366                 {

   367                     ERR("Unsupported evOffset[%d] value: %d (%d)\n", i, offset, evOffset[i]);

   368                     failed = AL_TRUE;

   369                 }

   370             }

   371         }

   372 

   373         if(!failed)

   374         {

   375             for(i = 0;i < HRIR_COUNT;i++)

   376             {

   377                 for(j = 0;j < HRIR_LENGTH;j++)

   378                 {

   379                     ALshort coeff;

   380                     coeff  = fgetc(f);

   381                     coeff |= fgetc(f)<<8;

   382                     newdata.coeffs[i][j] = coeff;

   383                 }

   384             }

   385             for(i = 0;i < HRIR_COUNT;i++)

   386             {

   387                 ALubyte delay;

   388                 delay = fgetc(f);

   389                 newdata.delays[i] = delay;

   390                 if(delay > maxDelay)

   391                 {

   392                     ERR("Invalid delay[%d]: %d (%d)\n", i, delay, maxDelay);

   393                     failed = AL_TRUE;

   394                 }

   395             }

   396 

   397             if(feof(f))

   398             {

   399                 ERR("Premature end of data\n");

   400                 failed = AL_TRUE;

   401             }

   402         }

   403 

   404         fclose(f);

   405         f = NULL;

   406 

   407         if(!failed)

   408             Hrtf = newdata;

   409         else

   410             ERR("Failed to load %s\n", fname);

   411     }

   412 }
author	Robert McIntyre <rlm@mit.edu>
date	Tue, 25 Oct 2011 13:03:35 -0700
parents	f9476ff7637e
children