1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/Alc/hrtf.c	Tue Oct 25 13:02:31 2011 -0700
     1.3 @@ -0,0 +1,412 @@
     1.4 +/**
     1.5 + * OpenAL cross platform audio library
     1.6 + * Copyright (C) 2011 by Chris Robinson
     1.7 + * This library is free software; you can redistribute it and/or
     1.8 + *  modify it under the terms of the GNU Library General Public
     1.9 + *  License as published by the Free Software Foundation; either
    1.10 + *  version 2 of the License, or (at your option) any later version.
    1.11 + *
    1.12 + * This library is distributed in the hope that it will be useful,
    1.13 + *  but WITHOUT ANY WARRANTY; without even the implied warranty of
    1.14 + *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    1.15 + *  Library General Public License for more details.
    1.16 + *
    1.17 + * You should have received a copy of the GNU Library General Public
    1.18 + *  License along with this library; if not, write to the
    1.19 + *  Free Software Foundation, Inc., 59 Temple Place - Suite 330,
    1.20 + *  Boston, MA  02111-1307, USA.
    1.21 + * Or go to http://www.gnu.org/copyleft/lgpl.html
    1.22 + */
    1.23 +
    1.24 +#include "config.h"
    1.25 +
    1.26 +#include "AL/al.h"
    1.27 +#include "AL/alc.h"
    1.28 +#include "alMain.h"
    1.29 +#include "alSource.h"
    1.30 +
    1.31 +/* External HRTF file format (LE byte order):
    1.32 + *
    1.33 + * ALchar   magic[8] = "MinPHR00";
    1.34 + * ALuint   sampleRate;
    1.35 + *
    1.36 + * ALushort hrirCount; // Required value: 828
    1.37 + * ALushort hrirSize;  // Required value: 32
    1.38 + * ALubyte  evCount;   // Required value: 19
    1.39 + *
    1.40 + * ALushort evOffset[evCount]; // Required values:
    1.41 + *   { 0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827 }
    1.42 + *
    1.43 + * ALushort coefficients[hrirCount][hrirSize];
    1.44 + * ALubyte  delays[hrirCount]; // Element values must not exceed 127
    1.45 + */
    1.46 +
    1.47 +static const ALchar magicMarker[8] = "MinPHR00";
    1.48 +
    1.49 +#define HRIR_COUNT 828
    1.50 +#define ELEV_COUNT 19
    1.51 +
    1.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 };
    1.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 };
    1.54 +
    1.55 +static struct Hrtf {
    1.56 +    ALuint sampleRate;
    1.57 +    ALshort coeffs[HRIR_COUNT][HRIR_LENGTH];
    1.58 +    ALubyte delays[HRIR_COUNT];
    1.59 +} Hrtf = {
    1.60 +    44100,
    1.61 +#include "hrtf_tables.inc"
    1.62 +};
    1.63 +
    1.64 +// Calculate the elevation indices given the polar elevation in radians.
    1.65 +// This will return two indices between 0 and (ELEV_COUNT-1) and an
    1.66 +// interpolation factor between 0.0 and 1.0.
    1.67 +static void CalcEvIndices(ALfloat ev, ALuint *evidx, ALfloat *evmu)
    1.68 +{
    1.69 +    ev = (M_PI/2.0f + ev) * (ELEV_COUNT-1) / M_PI;
    1.70 +    evidx[0] = (ALuint)ev;
    1.71 +    evidx[1] = minu(evidx[0] + 1, ELEV_COUNT-1);
    1.72 +    *evmu = ev - evidx[0];
    1.73 +}
    1.74 +
    1.75 +// Calculate the azimuth indices given the polar azimuth in radians.  This
    1.76 +// will return two indices between 0 and (azCount [ei] - 1) and an
    1.77 +// interpolation factor between 0.0 and 1.0.
    1.78 +static void CalcAzIndices(ALuint evidx, ALfloat az, ALuint *azidx, ALfloat *azmu)
    1.79 +{
    1.80 +    az = (M_PI*2.0f + az) * azCount[evidx] / (M_PI*2.0f);
    1.81 +    azidx[0] = (ALuint)az % azCount[evidx];
    1.82 +    azidx[1] = (azidx[0] + 1) % azCount[evidx];
    1.83 +    *azmu = az - floor(az);
    1.84 +}
    1.85 +
    1.86 +// Calculates the normalized HRTF transition factor (delta) from the changes
    1.87 +// in gain and listener to source angle between updates.  The result is a
    1.88 +// normalized delta factor than can be used to calculate moving HRIR stepping
    1.89 +// values.
    1.90 +ALfloat CalcHrtfDelta(ALfloat oldGain, ALfloat newGain, const ALfloat olddir[3], const ALfloat newdir[3])
    1.91 +{
    1.92 +    ALfloat gainChange, angleChange;
    1.93 +
    1.94 +    // Calculate the normalized dB gain change.
    1.95 +    newGain = maxf(newGain, 0.0001f);
    1.96 +    oldGain = maxf(oldGain, 0.0001f);
    1.97 +    gainChange = aluFabs(log10(newGain / oldGain) / log10(0.0001f));
    1.98 +
    1.99 +    // Calculate the normalized listener to source angle change when there is
   1.100 +    // enough gain to notice it.
   1.101 +    angleChange = 0.0f;
   1.102 +    if(gainChange > 0.0001f || newGain > 0.0001f)
   1.103 +    {
   1.104 +        // No angle change when the directions are equal or degenerate (when
   1.105 +        // both have zero length).
   1.106 +        if(newdir[0]-olddir[0] || newdir[1]-olddir[1] || newdir[2]-olddir[2])
   1.107 +            angleChange = aluAcos(olddir[0]*newdir[0] +
   1.108 +                                  olddir[1]*newdir[1] +
   1.109 +                                  olddir[2]*newdir[2]) / M_PI;
   1.110 +
   1.111 +    }
   1.112 +
   1.113 +    // Use the largest of the two changes for the delta factor, and apply a
   1.114 +    // significance shaping function to it.
   1.115 +    return clampf(angleChange*2.0f, gainChange*2.0f, 1.0f);
   1.116 +}
   1.117 +
   1.118 +// Calculates static HRIR coefficients and delays for the given polar
   1.119 +// elevation and azimuth in radians.  Linear interpolation is used to
   1.120 +// increase the apparent resolution of the HRIR dataset.  The coefficients
   1.121 +// are also normalized and attenuated by the specified gain.
   1.122 +void GetLerpedHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat (*coeffs)[2], ALuint *delays)
   1.123 +{
   1.124 +    ALuint evidx[2], azidx[2];
   1.125 +    ALfloat mu[3];
   1.126 +    ALuint lidx[4], ridx[4];
   1.127 +    ALuint i;
   1.128 +
   1.129 +    // Claculate elevation indices and interpolation factor.
   1.130 +    CalcEvIndices(elevation, evidx, &mu[2]);
   1.131 +
   1.132 +    // Calculate azimuth indices and interpolation factor for the first
   1.133 +    // elevation.
   1.134 +    CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
   1.135 +
   1.136 +    // Calculate the first set of linear HRIR indices for left and right
   1.137 +    // channels.
   1.138 +    lidx[0] = evOffset[evidx[0]] + azidx[0];
   1.139 +    lidx[1] = evOffset[evidx[0]] + azidx[1];
   1.140 +    ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
   1.141 +    ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
   1.142 +
   1.143 +    // Calculate azimuth indices and interpolation factor for the second
   1.144 +    // elevation.
   1.145 +    CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
   1.146 +
   1.147 +    // Calculate the second set of linear HRIR indices for left and right
   1.148 +    // channels.
   1.149 +    lidx[2] = evOffset[evidx[1]] + azidx[0];
   1.150 +    lidx[3] = evOffset[evidx[1]] + azidx[1];
   1.151 +    ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
   1.152 +    ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
   1.153 +
   1.154 +    // Calculate the normalized and attenuated HRIR coefficients using linear
   1.155 +    // interpolation when there is enough gain to warrant it.  Zero the
   1.156 +    // coefficients if gain is too low.
   1.157 +    if(gain > 0.0001f)
   1.158 +    {
   1.159 +        ALdouble scale = gain * (1.0/32767.0);
   1.160 +        for(i = 0;i < HRIR_LENGTH;i++)
   1.161 +        {
   1.162 +            coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),
   1.163 +                                lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),
   1.164 +                                mu[2]) * scale;
   1.165 +            coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),
   1.166 +                                lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),
   1.167 +                                mu[2]) * scale;
   1.168 +        }
   1.169 +    }
   1.170 +    else
   1.171 +    {
   1.172 +        for(i = 0;i < HRIR_LENGTH;i++)
   1.173 +        {
   1.174 +            coeffs[i][0] = 0.0f;
   1.175 +            coeffs[i][1] = 0.0f;
   1.176 +        }
   1.177 +    }
   1.178 +
   1.179 +    // Calculate the HRIR delays using linear interpolation.
   1.180 +    delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),
   1.181 +                              lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),
   1.182 +                              mu[2]) * 65536.0f);
   1.183 +    delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),
   1.184 +                              lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),
   1.185 +                              mu[2]) * 65536.0f);
   1.186 +}
   1.187 +
   1.188 +// Calculates the moving HRIR target coefficients, target delays, and
   1.189 +// stepping values for the given polar elevation and azimuth in radians.
   1.190 +// Linear interpolation is used to increase the apparent resolution of the
   1.191 +// HRIR dataset.  The coefficients are also normalized and attenuated by the
   1.192 +// specified gain.  Stepping resolution and count is determined using the
   1.193 +// given delta factor between 0.0 and 1.0.
   1.194 +ALuint GetMovingHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat delta, ALint counter, ALfloat (*coeffs)[2], ALuint *delays, ALfloat (*coeffStep)[2], ALint *delayStep)
   1.195 +{
   1.196 +    ALuint evidx[2], azidx[2];
   1.197 +    ALuint lidx[4], ridx[4];
   1.198 +    ALfloat left, right;
   1.199 +    ALfloat mu[3];
   1.200 +    ALfloat step;
   1.201 +    ALuint i;
   1.202 +
   1.203 +    // Claculate elevation indices and interpolation factor.
   1.204 +    CalcEvIndices(elevation, evidx, &mu[2]);
   1.205 +
   1.206 +    // Calculate azimuth indices and interpolation factor for the first
   1.207 +    // elevation.
   1.208 +    CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
   1.209 +
   1.210 +    // Calculate the first set of linear HRIR indices for left and right
   1.211 +    // channels.
   1.212 +    lidx[0] = evOffset[evidx[0]] + azidx[0];
   1.213 +    lidx[1] = evOffset[evidx[0]] + azidx[1];
   1.214 +    ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
   1.215 +    ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
   1.216 +
   1.217 +    // Calculate azimuth indices and interpolation factor for the second
   1.218 +    // elevation.
   1.219 +    CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
   1.220 +
   1.221 +    // Calculate the second set of linear HRIR indices for left and right
   1.222 +    // channels.
   1.223 +    lidx[2] = evOffset[evidx[1]] + azidx[0];
   1.224 +    lidx[3] = evOffset[evidx[1]] + azidx[1];
   1.225 +    ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
   1.226 +    ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
   1.227 +
   1.228 +    // Calculate the stepping parameters.
   1.229 +    delta = maxf(floor(delta*(Hrtf.sampleRate*0.015f) + 0.5), 1.0f);
   1.230 +    step = 1.0f / delta;
   1.231 +
   1.232 +    // Calculate the normalized and attenuated target HRIR coefficients using
   1.233 +    // linear interpolation when there is enough gain to warrant it.  Zero
   1.234 +    // the target coefficients if gain is too low.  Then calculate the
   1.235 +    // coefficient stepping values using the target and previous running
   1.236 +    // coefficients.
   1.237 +    if(gain > 0.0001f)
   1.238 +    {
   1.239 +        ALdouble scale = gain * (1.0/32767.0);
   1.240 +        for(i = 0;i < HRIR_LENGTH;i++)
   1.241 +        {
   1.242 +            left = coeffs[i][0] - (coeffStep[i][0] * counter);
   1.243 +            right = coeffs[i][1] - (coeffStep[i][1] * counter);
   1.244 +
   1.245 +            coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),
   1.246 +                                lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),
   1.247 +                                mu[2]) * scale;
   1.248 +            coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),
   1.249 +                                lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),
   1.250 +                                mu[2]) * scale;
   1.251 +
   1.252 +            coeffStep[i][0] = step * (coeffs[i][0] - left);
   1.253 +            coeffStep[i][1] = step * (coeffs[i][1] - right);
   1.254 +        }
   1.255 +    }
   1.256 +    else
   1.257 +    {
   1.258 +        for(i = 0;i < HRIR_LENGTH;i++)
   1.259 +        {
   1.260 +            left = coeffs[i][0] - (coeffStep[i][0] * counter);
   1.261 +            right = coeffs[i][1] - (coeffStep[i][1] * counter);
   1.262 +
   1.263 +            coeffs[i][0] = 0.0f;
   1.264 +            coeffs[i][1] = 0.0f;
   1.265 +
   1.266 +            coeffStep[i][0] = step * -left;
   1.267 +            coeffStep[i][1] = step * -right;
   1.268 +        }
   1.269 +    }
   1.270 +
   1.271 +    // Calculate the HRIR delays using linear interpolation.  Then calculate
   1.272 +    // the delay stepping values using the target and previous running
   1.273 +    // delays.
   1.274 +    left = delays[0] - (delayStep[0] * counter);
   1.275 +    right = delays[1] - (delayStep[1] * counter);
   1.276 +
   1.277 +    delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),
   1.278 +                              lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),
   1.279 +                              mu[2]) * 65536.0f);
   1.280 +    delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),
   1.281 +                              lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),
   1.282 +                              mu[2]) * 65536.0f);
   1.283 +
   1.284 +    delayStep[0] = (ALint)(step * (delays[0] - left));
   1.285 +    delayStep[1] = (ALint)(step * (delays[1] - right));
   1.286 +
   1.287 +    // The stepping count is the number of samples necessary for the HRIR to
   1.288 +    // complete its transition.  The mixer will only apply stepping for this
   1.289 +    // many samples.
   1.290 +    return (ALuint)delta;
   1.291 +}
   1.292 +
   1.293 +ALCboolean IsHrtfCompatible(ALCdevice *device)
   1.294 +{
   1.295 +    if(device->FmtChans == DevFmtStereo && device->Frequency == Hrtf.sampleRate)
   1.296 +        return ALC_TRUE;
   1.297 +    ERR("Incompatible HRTF format: %s %uhz (%s %uhz needed)\n",
   1.298 +        DevFmtChannelsString(device->FmtChans), device->Frequency,
   1.299 +        DevFmtChannelsString(DevFmtStereo), Hrtf.sampleRate);
   1.300 +    return ALC_FALSE;
   1.301 +}
   1.302 +
   1.303 +void InitHrtf(void)
   1.304 +{
   1.305 +    const char *fname;
   1.306 +    FILE *f = NULL;
   1.307 +
   1.308 +    fname = GetConfigValue(NULL, "hrtf_tables", "");
   1.309 +    if(fname[0] != '\0')
   1.310 +    {
   1.311 +        f = fopen(fname, "rb");
   1.312 +        if(f == NULL)
   1.313 +            ERR("Could not open %s\n", fname);
   1.314 +    }
   1.315 +    if(f != NULL)
   1.316 +    {
   1.317 +        const ALubyte maxDelay = SRC_HISTORY_LENGTH-1;
   1.318 +        ALboolean failed = AL_FALSE;
   1.319 +        struct Hrtf newdata;
   1.320 +        ALchar magic[9];
   1.321 +        ALsizei i, j;
   1.322 +
   1.323 +        if(fread(magic, 1, sizeof(magicMarker), f) != sizeof(magicMarker))
   1.324 +        {
   1.325 +            ERR("Failed to read magic marker\n");
   1.326 +            failed = AL_TRUE;
   1.327 +        }
   1.328 +        else if(memcmp(magic, magicMarker, sizeof(magicMarker)) != 0)
   1.329 +        {
   1.330 +            magic[8] = 0;
   1.331 +            ERR("Invalid magic marker: \"%s\"\n", magic);
   1.332 +            failed = AL_TRUE;
   1.333 +        }
   1.334 +
   1.335 +        if(!failed)
   1.336 +        {
   1.337 +            ALushort hrirCount, hrirSize;
   1.338 +            ALubyte  evCount;
   1.339 +
   1.340 +            newdata.sampleRate  = fgetc(f);
   1.341 +            newdata.sampleRate |= fgetc(f)<<8;
   1.342 +            newdata.sampleRate |= fgetc(f)<<16;
   1.343 +            newdata.sampleRate |= fgetc(f)<<24;
   1.344 +
   1.345 +            hrirCount  = fgetc(f);
   1.346 +            hrirCount |= fgetc(f)<<8;
   1.347 +
   1.348 +            hrirSize  = fgetc(f);
   1.349 +            hrirSize |= fgetc(f)<<8;
   1.350 +
   1.351 +            evCount = fgetc(f);
   1.352 +
   1.353 +            if(hrirCount != HRIR_COUNT || hrirSize != HRIR_LENGTH || evCount != ELEV_COUNT)
   1.354 +            {
   1.355 +                ERR("Unsupported value: hrirCount=%d (%d), hrirSize=%d (%d), evCount=%d (%d)\n",
   1.356 +                    hrirCount, HRIR_COUNT, hrirSize, HRIR_LENGTH, evCount, ELEV_COUNT);
   1.357 +                failed = AL_TRUE;
   1.358 +            }
   1.359 +        }
   1.360 +
   1.361 +        if(!failed)
   1.362 +        {
   1.363 +            for(i = 0;i < HRIR_COUNT;i++)
   1.364 +            {
   1.365 +                ALushort offset;
   1.366 +                offset  = fgetc(f);
   1.367 +                offset |= fgetc(f)<<8;
   1.368 +                if(offset != evOffset[i])
   1.369 +                {
   1.370 +                    ERR("Unsupported evOffset[%d] value: %d (%d)\n", i, offset, evOffset[i]);
   1.371 +                    failed = AL_TRUE;
   1.372 +                }
   1.373 +            }
   1.374 +        }
   1.375 +
   1.376 +        if(!failed)
   1.377 +        {
   1.378 +            for(i = 0;i < HRIR_COUNT;i++)
   1.379 +            {
   1.380 +                for(j = 0;j < HRIR_LENGTH;j++)
   1.381 +                {
   1.382 +                    ALshort coeff;
   1.383 +                    coeff  = fgetc(f);
   1.384 +                    coeff |= fgetc(f)<<8;
   1.385 +                    newdata.coeffs[i][j] = coeff;
   1.386 +                }
   1.387 +            }
   1.388 +            for(i = 0;i < HRIR_COUNT;i++)
   1.389 +            {
   1.390 +                ALubyte delay;
   1.391 +                delay = fgetc(f);
   1.392 +                newdata.delays[i] = delay;
   1.393 +                if(delay > maxDelay)
   1.394 +                {
   1.395 +                    ERR("Invalid delay[%d]: %d (%d)\n", i, delay, maxDelay);
   1.396 +                    failed = AL_TRUE;
   1.397 +                }
   1.398 +            }
   1.399 +
   1.400 +            if(feof(f))
   1.401 +            {
   1.402 +                ERR("Premature end of data\n");
   1.403 +                failed = AL_TRUE;
   1.404 +            }
   1.405 +        }
   1.406 +
   1.407 +        fclose(f);
   1.408 +        f = NULL;
   1.409 +
   1.410 +        if(!failed)
   1.411 +            Hrtf = newdata;
   1.412 +        else
   1.413 +            ERR("Failed to load %s\n", fname);
   1.414 +    }
   1.415 +}