annotate Alc/hrtf.c @ 29:cdf320cb5949

updated test suite
author Robert McIntyre <rlm@mit.edu>
date Sat, 10 Dec 2011 21:42:50 -0600
parents f9476ff7637e
children
rev   line source
rlm@0 1 /**
rlm@0 2 * OpenAL cross platform audio library
rlm@0 3 * Copyright (C) 2011 by Chris Robinson
rlm@0 4 * This library is free software; you can redistribute it and/or
rlm@0 5 * modify it under the terms of the GNU Library General Public
rlm@0 6 * License as published by the Free Software Foundation; either
rlm@0 7 * version 2 of the License, or (at your option) any later version.
rlm@0 8 *
rlm@0 9 * This library is distributed in the hope that it will be useful,
rlm@0 10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
rlm@0 11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
rlm@0 12 * Library General Public License for more details.
rlm@0 13 *
rlm@0 14 * You should have received a copy of the GNU Library General Public
rlm@0 15 * License along with this library; if not, write to the
rlm@0 16 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
rlm@0 17 * Boston, MA 02111-1307, USA.
rlm@0 18 * Or go to http://www.gnu.org/copyleft/lgpl.html
rlm@0 19 */
rlm@0 20
rlm@0 21 #include "config.h"
rlm@0 22
rlm@0 23 #include "AL/al.h"
rlm@0 24 #include "AL/alc.h"
rlm@0 25 #include "alMain.h"
rlm@0 26 #include "alSource.h"
rlm@0 27
rlm@0 28 /* External HRTF file format (LE byte order):
rlm@0 29 *
rlm@0 30 * ALchar magic[8] = "MinPHR00";
rlm@0 31 * ALuint sampleRate;
rlm@0 32 *
rlm@0 33 * ALushort hrirCount; // Required value: 828
rlm@0 34 * ALushort hrirSize; // Required value: 32
rlm@0 35 * ALubyte evCount; // Required value: 19
rlm@0 36 *
rlm@0 37 * ALushort evOffset[evCount]; // Required values:
rlm@0 38 * { 0, 1, 13, 37, 73, 118, 174, 234, 306, 378, 450, 522, 594, 654, 710, 755, 791, 815, 827 }
rlm@0 39 *
rlm@0 40 * ALushort coefficients[hrirCount][hrirSize];
rlm@0 41 * ALubyte delays[hrirCount]; // Element values must not exceed 127
rlm@0 42 */
rlm@0 43
rlm@0 44 static const ALchar magicMarker[8] = "MinPHR00";
rlm@0 45
rlm@0 46 #define HRIR_COUNT 828
rlm@0 47 #define ELEV_COUNT 19
rlm@0 48
rlm@0 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 };
rlm@0 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 };
rlm@0 51
rlm@0 52 static struct Hrtf {
rlm@0 53 ALuint sampleRate;
rlm@0 54 ALshort coeffs[HRIR_COUNT][HRIR_LENGTH];
rlm@0 55 ALubyte delays[HRIR_COUNT];
rlm@0 56 } Hrtf = {
rlm@0 57 44100,
rlm@0 58 #include "hrtf_tables.inc"
rlm@0 59 };
rlm@0 60
rlm@0 61 // Calculate the elevation indices given the polar elevation in radians.
rlm@0 62 // This will return two indices between 0 and (ELEV_COUNT-1) and an
rlm@0 63 // interpolation factor between 0.0 and 1.0.
rlm@0 64 static void CalcEvIndices(ALfloat ev, ALuint *evidx, ALfloat *evmu)
rlm@0 65 {
rlm@0 66 ev = (M_PI/2.0f + ev) * (ELEV_COUNT-1) / M_PI;
rlm@0 67 evidx[0] = (ALuint)ev;
rlm@0 68 evidx[1] = minu(evidx[0] + 1, ELEV_COUNT-1);
rlm@0 69 *evmu = ev - evidx[0];
rlm@0 70 }
rlm@0 71
rlm@0 72 // Calculate the azimuth indices given the polar azimuth in radians. This
rlm@0 73 // will return two indices between 0 and (azCount [ei] - 1) and an
rlm@0 74 // interpolation factor between 0.0 and 1.0.
rlm@0 75 static void CalcAzIndices(ALuint evidx, ALfloat az, ALuint *azidx, ALfloat *azmu)
rlm@0 76 {
rlm@0 77 az = (M_PI*2.0f + az) * azCount[evidx] / (M_PI*2.0f);
rlm@0 78 azidx[0] = (ALuint)az % azCount[evidx];
rlm@0 79 azidx[1] = (azidx[0] + 1) % azCount[evidx];
rlm@0 80 *azmu = az - floor(az);
rlm@0 81 }
rlm@0 82
rlm@0 83 // Calculates the normalized HRTF transition factor (delta) from the changes
rlm@0 84 // in gain and listener to source angle between updates. The result is a
rlm@0 85 // normalized delta factor than can be used to calculate moving HRIR stepping
rlm@0 86 // values.
rlm@0 87 ALfloat CalcHrtfDelta(ALfloat oldGain, ALfloat newGain, const ALfloat olddir[3], const ALfloat newdir[3])
rlm@0 88 {
rlm@0 89 ALfloat gainChange, angleChange;
rlm@0 90
rlm@0 91 // Calculate the normalized dB gain change.
rlm@0 92 newGain = maxf(newGain, 0.0001f);
rlm@0 93 oldGain = maxf(oldGain, 0.0001f);
rlm@0 94 gainChange = aluFabs(log10(newGain / oldGain) / log10(0.0001f));
rlm@0 95
rlm@0 96 // Calculate the normalized listener to source angle change when there is
rlm@0 97 // enough gain to notice it.
rlm@0 98 angleChange = 0.0f;
rlm@0 99 if(gainChange > 0.0001f || newGain > 0.0001f)
rlm@0 100 {
rlm@0 101 // No angle change when the directions are equal or degenerate (when
rlm@0 102 // both have zero length).
rlm@0 103 if(newdir[0]-olddir[0] || newdir[1]-olddir[1] || newdir[2]-olddir[2])
rlm@0 104 angleChange = aluAcos(olddir[0]*newdir[0] +
rlm@0 105 olddir[1]*newdir[1] +
rlm@0 106 olddir[2]*newdir[2]) / M_PI;
rlm@0 107
rlm@0 108 }
rlm@0 109
rlm@0 110 // Use the largest of the two changes for the delta factor, and apply a
rlm@0 111 // significance shaping function to it.
rlm@0 112 return clampf(angleChange*2.0f, gainChange*2.0f, 1.0f);
rlm@0 113 }
rlm@0 114
rlm@0 115 // Calculates static HRIR coefficients and delays for the given polar
rlm@0 116 // elevation and azimuth in radians. Linear interpolation is used to
rlm@0 117 // increase the apparent resolution of the HRIR dataset. The coefficients
rlm@0 118 // are also normalized and attenuated by the specified gain.
rlm@0 119 void GetLerpedHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat (*coeffs)[2], ALuint *delays)
rlm@0 120 {
rlm@0 121 ALuint evidx[2], azidx[2];
rlm@0 122 ALfloat mu[3];
rlm@0 123 ALuint lidx[4], ridx[4];
rlm@0 124 ALuint i;
rlm@0 125
rlm@0 126 // Claculate elevation indices and interpolation factor.
rlm@0 127 CalcEvIndices(elevation, evidx, &mu[2]);
rlm@0 128
rlm@0 129 // Calculate azimuth indices and interpolation factor for the first
rlm@0 130 // elevation.
rlm@0 131 CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
rlm@0 132
rlm@0 133 // Calculate the first set of linear HRIR indices for left and right
rlm@0 134 // channels.
rlm@0 135 lidx[0] = evOffset[evidx[0]] + azidx[0];
rlm@0 136 lidx[1] = evOffset[evidx[0]] + azidx[1];
rlm@0 137 ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
rlm@0 138 ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
rlm@0 139
rlm@0 140 // Calculate azimuth indices and interpolation factor for the second
rlm@0 141 // elevation.
rlm@0 142 CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
rlm@0 143
rlm@0 144 // Calculate the second set of linear HRIR indices for left and right
rlm@0 145 // channels.
rlm@0 146 lidx[2] = evOffset[evidx[1]] + azidx[0];
rlm@0 147 lidx[3] = evOffset[evidx[1]] + azidx[1];
rlm@0 148 ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
rlm@0 149 ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
rlm@0 150
rlm@0 151 // Calculate the normalized and attenuated HRIR coefficients using linear
rlm@0 152 // interpolation when there is enough gain to warrant it. Zero the
rlm@0 153 // coefficients if gain is too low.
rlm@0 154 if(gain > 0.0001f)
rlm@0 155 {
rlm@0 156 ALdouble scale = gain * (1.0/32767.0);
rlm@0 157 for(i = 0;i < HRIR_LENGTH;i++)
rlm@0 158 {
rlm@0 159 coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),
rlm@0 160 lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),
rlm@0 161 mu[2]) * scale;
rlm@0 162 coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),
rlm@0 163 lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),
rlm@0 164 mu[2]) * scale;
rlm@0 165 }
rlm@0 166 }
rlm@0 167 else
rlm@0 168 {
rlm@0 169 for(i = 0;i < HRIR_LENGTH;i++)
rlm@0 170 {
rlm@0 171 coeffs[i][0] = 0.0f;
rlm@0 172 coeffs[i][1] = 0.0f;
rlm@0 173 }
rlm@0 174 }
rlm@0 175
rlm@0 176 // Calculate the HRIR delays using linear interpolation.
rlm@0 177 delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),
rlm@0 178 lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),
rlm@0 179 mu[2]) * 65536.0f);
rlm@0 180 delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),
rlm@0 181 lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),
rlm@0 182 mu[2]) * 65536.0f);
rlm@0 183 }
rlm@0 184
rlm@0 185 // Calculates the moving HRIR target coefficients, target delays, and
rlm@0 186 // stepping values for the given polar elevation and azimuth in radians.
rlm@0 187 // Linear interpolation is used to increase the apparent resolution of the
rlm@0 188 // HRIR dataset. The coefficients are also normalized and attenuated by the
rlm@0 189 // specified gain. Stepping resolution and count is determined using the
rlm@0 190 // given delta factor between 0.0 and 1.0.
rlm@0 191 ALuint GetMovingHrtfCoeffs(ALfloat elevation, ALfloat azimuth, ALfloat gain, ALfloat delta, ALint counter, ALfloat (*coeffs)[2], ALuint *delays, ALfloat (*coeffStep)[2], ALint *delayStep)
rlm@0 192 {
rlm@0 193 ALuint evidx[2], azidx[2];
rlm@0 194 ALuint lidx[4], ridx[4];
rlm@0 195 ALfloat left, right;
rlm@0 196 ALfloat mu[3];
rlm@0 197 ALfloat step;
rlm@0 198 ALuint i;
rlm@0 199
rlm@0 200 // Claculate elevation indices and interpolation factor.
rlm@0 201 CalcEvIndices(elevation, evidx, &mu[2]);
rlm@0 202
rlm@0 203 // Calculate azimuth indices and interpolation factor for the first
rlm@0 204 // elevation.
rlm@0 205 CalcAzIndices(evidx[0], azimuth, azidx, &mu[0]);
rlm@0 206
rlm@0 207 // Calculate the first set of linear HRIR indices for left and right
rlm@0 208 // channels.
rlm@0 209 lidx[0] = evOffset[evidx[0]] + azidx[0];
rlm@0 210 lidx[1] = evOffset[evidx[0]] + azidx[1];
rlm@0 211 ridx[0] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[0]) % azCount[evidx[0]]);
rlm@0 212 ridx[1] = evOffset[evidx[0]] + ((azCount[evidx[0]]-azidx[1]) % azCount[evidx[0]]);
rlm@0 213
rlm@0 214 // Calculate azimuth indices and interpolation factor for the second
rlm@0 215 // elevation.
rlm@0 216 CalcAzIndices(evidx[1], azimuth, azidx, &mu[1]);
rlm@0 217
rlm@0 218 // Calculate the second set of linear HRIR indices for left and right
rlm@0 219 // channels.
rlm@0 220 lidx[2] = evOffset[evidx[1]] + azidx[0];
rlm@0 221 lidx[3] = evOffset[evidx[1]] + azidx[1];
rlm@0 222 ridx[2] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[0]) % azCount[evidx[1]]);
rlm@0 223 ridx[3] = evOffset[evidx[1]] + ((azCount[evidx[1]]-azidx[1]) % azCount[evidx[1]]);
rlm@0 224
rlm@0 225 // Calculate the stepping parameters.
rlm@0 226 delta = maxf(floor(delta*(Hrtf.sampleRate*0.015f) + 0.5), 1.0f);
rlm@0 227 step = 1.0f / delta;
rlm@0 228
rlm@0 229 // Calculate the normalized and attenuated target HRIR coefficients using
rlm@0 230 // linear interpolation when there is enough gain to warrant it. Zero
rlm@0 231 // the target coefficients if gain is too low. Then calculate the
rlm@0 232 // coefficient stepping values using the target and previous running
rlm@0 233 // coefficients.
rlm@0 234 if(gain > 0.0001f)
rlm@0 235 {
rlm@0 236 ALdouble scale = gain * (1.0/32767.0);
rlm@0 237 for(i = 0;i < HRIR_LENGTH;i++)
rlm@0 238 {
rlm@0 239 left = coeffs[i][0] - (coeffStep[i][0] * counter);
rlm@0 240 right = coeffs[i][1] - (coeffStep[i][1] * counter);
rlm@0 241
rlm@0 242 coeffs[i][0] = lerp(lerp(Hrtf.coeffs[lidx[0]][i], Hrtf.coeffs[lidx[1]][i], mu[0]),
rlm@0 243 lerp(Hrtf.coeffs[lidx[2]][i], Hrtf.coeffs[lidx[3]][i], mu[1]),
rlm@0 244 mu[2]) * scale;
rlm@0 245 coeffs[i][1] = lerp(lerp(Hrtf.coeffs[ridx[0]][i], Hrtf.coeffs[ridx[1]][i], mu[0]),
rlm@0 246 lerp(Hrtf.coeffs[ridx[2]][i], Hrtf.coeffs[ridx[3]][i], mu[1]),
rlm@0 247 mu[2]) * scale;
rlm@0 248
rlm@0 249 coeffStep[i][0] = step * (coeffs[i][0] - left);
rlm@0 250 coeffStep[i][1] = step * (coeffs[i][1] - right);
rlm@0 251 }
rlm@0 252 }
rlm@0 253 else
rlm@0 254 {
rlm@0 255 for(i = 0;i < HRIR_LENGTH;i++)
rlm@0 256 {
rlm@0 257 left = coeffs[i][0] - (coeffStep[i][0] * counter);
rlm@0 258 right = coeffs[i][1] - (coeffStep[i][1] * counter);
rlm@0 259
rlm@0 260 coeffs[i][0] = 0.0f;
rlm@0 261 coeffs[i][1] = 0.0f;
rlm@0 262
rlm@0 263 coeffStep[i][0] = step * -left;
rlm@0 264 coeffStep[i][1] = step * -right;
rlm@0 265 }
rlm@0 266 }
rlm@0 267
rlm@0 268 // Calculate the HRIR delays using linear interpolation. Then calculate
rlm@0 269 // the delay stepping values using the target and previous running
rlm@0 270 // delays.
rlm@0 271 left = delays[0] - (delayStep[0] * counter);
rlm@0 272 right = delays[1] - (delayStep[1] * counter);
rlm@0 273
rlm@0 274 delays[0] = (ALuint)(lerp(lerp(Hrtf.delays[lidx[0]], Hrtf.delays[lidx[1]], mu[0]),
rlm@0 275 lerp(Hrtf.delays[lidx[2]], Hrtf.delays[lidx[3]], mu[1]),
rlm@0 276 mu[2]) * 65536.0f);
rlm@0 277 delays[1] = (ALuint)(lerp(lerp(Hrtf.delays[ridx[0]], Hrtf.delays[ridx[1]], mu[0]),
rlm@0 278 lerp(Hrtf.delays[ridx[2]], Hrtf.delays[ridx[3]], mu[1]),
rlm@0 279 mu[2]) * 65536.0f);
rlm@0 280
rlm@0 281 delayStep[0] = (ALint)(step * (delays[0] - left));
rlm@0 282 delayStep[1] = (ALint)(step * (delays[1] - right));
rlm@0 283
rlm@0 284 // The stepping count is the number of samples necessary for the HRIR to
rlm@0 285 // complete its transition. The mixer will only apply stepping for this
rlm@0 286 // many samples.
rlm@0 287 return (ALuint)delta;
rlm@0 288 }
rlm@0 289
rlm@0 290 ALCboolean IsHrtfCompatible(ALCdevice *device)
rlm@0 291 {
rlm@0 292 if(device->FmtChans == DevFmtStereo && device->Frequency == Hrtf.sampleRate)
rlm@0 293 return ALC_TRUE;
rlm@0 294 ERR("Incompatible HRTF format: %s %uhz (%s %uhz needed)\n",
rlm@0 295 DevFmtChannelsString(device->FmtChans), device->Frequency,
rlm@0 296 DevFmtChannelsString(DevFmtStereo), Hrtf.sampleRate);
rlm@0 297 return ALC_FALSE;
rlm@0 298 }
rlm@0 299
rlm@0 300 void InitHrtf(void)
rlm@0 301 {
rlm@0 302 const char *fname;
rlm@0 303 FILE *f = NULL;
rlm@0 304
rlm@0 305 fname = GetConfigValue(NULL, "hrtf_tables", "");
rlm@0 306 if(fname[0] != '\0')
rlm@0 307 {
rlm@0 308 f = fopen(fname, "rb");
rlm@0 309 if(f == NULL)
rlm@0 310 ERR("Could not open %s\n", fname);
rlm@0 311 }
rlm@0 312 if(f != NULL)
rlm@0 313 {
rlm@0 314 const ALubyte maxDelay = SRC_HISTORY_LENGTH-1;
rlm@0 315 ALboolean failed = AL_FALSE;
rlm@0 316 struct Hrtf newdata;
rlm@0 317 ALchar magic[9];
rlm@0 318 ALsizei i, j;
rlm@0 319
rlm@0 320 if(fread(magic, 1, sizeof(magicMarker), f) != sizeof(magicMarker))
rlm@0 321 {
rlm@0 322 ERR("Failed to read magic marker\n");
rlm@0 323 failed = AL_TRUE;
rlm@0 324 }
rlm@0 325 else if(memcmp(magic, magicMarker, sizeof(magicMarker)) != 0)
rlm@0 326 {
rlm@0 327 magic[8] = 0;
rlm@0 328 ERR("Invalid magic marker: \"%s\"\n", magic);
rlm@0 329 failed = AL_TRUE;
rlm@0 330 }
rlm@0 331
rlm@0 332 if(!failed)
rlm@0 333 {
rlm@0 334 ALushort hrirCount, hrirSize;
rlm@0 335 ALubyte evCount;
rlm@0 336
rlm@0 337 newdata.sampleRate = fgetc(f);
rlm@0 338 newdata.sampleRate |= fgetc(f)<<8;
rlm@0 339 newdata.sampleRate |= fgetc(f)<<16;
rlm@0 340 newdata.sampleRate |= fgetc(f)<<24;
rlm@0 341
rlm@0 342 hrirCount = fgetc(f);
rlm@0 343 hrirCount |= fgetc(f)<<8;
rlm@0 344
rlm@0 345 hrirSize = fgetc(f);
rlm@0 346 hrirSize |= fgetc(f)<<8;
rlm@0 347
rlm@0 348 evCount = fgetc(f);
rlm@0 349
rlm@0 350 if(hrirCount != HRIR_COUNT || hrirSize != HRIR_LENGTH || evCount != ELEV_COUNT)
rlm@0 351 {
rlm@0 352 ERR("Unsupported value: hrirCount=%d (%d), hrirSize=%d (%d), evCount=%d (%d)\n",
rlm@0 353 hrirCount, HRIR_COUNT, hrirSize, HRIR_LENGTH, evCount, ELEV_COUNT);
rlm@0 354 failed = AL_TRUE;
rlm@0 355 }
rlm@0 356 }
rlm@0 357
rlm@0 358 if(!failed)
rlm@0 359 {
rlm@0 360 for(i = 0;i < HRIR_COUNT;i++)
rlm@0 361 {
rlm@0 362 ALushort offset;
rlm@0 363 offset = fgetc(f);
rlm@0 364 offset |= fgetc(f)<<8;
rlm@0 365 if(offset != evOffset[i])
rlm@0 366 {
rlm@0 367 ERR("Unsupported evOffset[%d] value: %d (%d)\n", i, offset, evOffset[i]);
rlm@0 368 failed = AL_TRUE;
rlm@0 369 }
rlm@0 370 }
rlm@0 371 }
rlm@0 372
rlm@0 373 if(!failed)
rlm@0 374 {
rlm@0 375 for(i = 0;i < HRIR_COUNT;i++)
rlm@0 376 {
rlm@0 377 for(j = 0;j < HRIR_LENGTH;j++)
rlm@0 378 {
rlm@0 379 ALshort coeff;
rlm@0 380 coeff = fgetc(f);
rlm@0 381 coeff |= fgetc(f)<<8;
rlm@0 382 newdata.coeffs[i][j] = coeff;
rlm@0 383 }
rlm@0 384 }
rlm@0 385 for(i = 0;i < HRIR_COUNT;i++)
rlm@0 386 {
rlm@0 387 ALubyte delay;
rlm@0 388 delay = fgetc(f);
rlm@0 389 newdata.delays[i] = delay;
rlm@0 390 if(delay > maxDelay)
rlm@0 391 {
rlm@0 392 ERR("Invalid delay[%d]: %d (%d)\n", i, delay, maxDelay);
rlm@0 393 failed = AL_TRUE;
rlm@0 394 }
rlm@0 395 }
rlm@0 396
rlm@0 397 if(feof(f))
rlm@0 398 {
rlm@0 399 ERR("Premature end of data\n");
rlm@0 400 failed = AL_TRUE;
rlm@0 401 }
rlm@0 402 }
rlm@0 403
rlm@0 404 fclose(f);
rlm@0 405 f = NULL;
rlm@0 406
rlm@0 407 if(!failed)
rlm@0 408 Hrtf = newdata;
rlm@0 409 else
rlm@0 410 ERR("Failed to load %s\n", fname);
rlm@0 411 }
rlm@0 412 }