1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/content/media/webaudio/blink/HRTFElevation.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,328 @@
1.4 +/*
1.5 + * Copyright (C) 2010 Google Inc. All rights reserved.
1.6 + *
1.7 + * Redistribution and use in source and binary forms, with or without
1.8 + * modification, are permitted provided that the following conditions
1.9 + * are met:
1.10 + *
1.11 + * 1. Redistributions of source code must retain the above copyright
1.12 + * notice, this list of conditions and the following disclaimer.
1.13 + * 2. Redistributions in binary form must reproduce the above copyright
1.14 + * notice, this list of conditions and the following disclaimer in the
1.15 + * documentation and/or other materials provided with the distribution.
1.16 + * 3. Neither the name of Apple Computer, Inc. ("Apple") nor the names of
1.17 + * its contributors may be used to endorse or promote products derived
1.18 + * from this software without specific prior written permission.
1.19 + *
1.20 + * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
1.21 + * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
1.22 + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
1.23 + * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
1.24 + * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
1.25 + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
1.26 + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
1.27 + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
1.28 + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
1.29 + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1.30 + */
1.31 +
1.32 +#include "HRTFElevation.h"
1.33 +
1.34 +#include "speex/speex_resampler.h"
1.35 +#include "mozilla/PodOperations.h"
1.36 +#include "AudioSampleFormat.h"
1.37 +
1.38 +#include "IRC_Composite_C_R0195-incl.cpp"
1.39 +
1.40 +using namespace std;
1.41 +using namespace mozilla;
1.42 +
1.43 +namespace WebCore {
1.44 +
1.45 +const int elevationSpacing = irc_composite_c_r0195_elevation_interval;
1.46 +const int firstElevation = irc_composite_c_r0195_first_elevation;
1.47 +const int numberOfElevations = MOZ_ARRAY_LENGTH(irc_composite_c_r0195);
1.48 +
1.49 +const unsigned HRTFElevation::NumberOfTotalAzimuths = 360 / 15 * 8;
1.50 +
1.51 +const int rawSampleRate = irc_composite_c_r0195_sample_rate;
1.52 +
1.53 +// Number of frames in an individual impulse response.
1.54 +const size_t ResponseFrameSize = 256;
1.55 +
1.56 +size_t HRTFElevation::sizeOfIncludingThis(mozilla::MallocSizeOf aMallocSizeOf) const
1.57 +{
1.58 + size_t amount = aMallocSizeOf(this);
1.59 +
1.60 + amount += m_kernelListL.SizeOfExcludingThis(aMallocSizeOf);
1.61 + for (size_t i = 0; i < m_kernelListL.Length(); i++) {
1.62 + amount += m_kernelListL[i]->sizeOfIncludingThis(aMallocSizeOf);
1.63 + }
1.64 +
1.65 + return amount;
1.66 +}
1.67 +
1.68 +size_t HRTFElevation::fftSizeForSampleRate(float sampleRate)
1.69 +{
1.70 + // The IRCAM HRTF impulse responses were 512 sample-frames @44.1KHz,
1.71 + // but these have been truncated to 256 samples.
1.72 + // An FFT-size of twice impulse response size is used (for convolution).
1.73 + // So for sample rates of 44.1KHz an FFT size of 512 is good.
1.74 + // We double the FFT-size only for sample rates at least double this.
1.75 + // If the FFT size is too large then the impulse response will be padded
1.76 + // with zeros without the fade-out provided by HRTFKernel.
1.77 + MOZ_ASSERT(sampleRate > 1.0 && sampleRate < 1048576.0);
1.78 +
1.79 + // This is the size if we were to use all raw response samples.
1.80 + unsigned resampledLength =
1.81 + floorf(ResponseFrameSize * sampleRate / rawSampleRate);
1.82 + // Keep things semi-sane, with max FFT size of 1024 and minimum of 4.
1.83 + // "size |= 3" ensures a minimum of 4 (with the size++ below) and sets the
1.84 + // 2 least significant bits for rounding up to the next power of 2 below.
1.85 + unsigned size = min(resampledLength, 1023U);
1.86 + size |= 3;
1.87 + // Round up to the next power of 2, making the FFT size no more than twice
1.88 + // the impulse response length. This doubles size for values that are
1.89 + // already powers of 2. This works by filling in 7 bits to right of the
1.90 + // most significant bit. The most significant bit is no greater than
1.91 + // 1 << 9, and the least significant 2 bits were already set above.
1.92 + size |= (size >> 1);
1.93 + size |= (size >> 2);
1.94 + size |= (size >> 4);
1.95 + size++;
1.96 + MOZ_ASSERT((size & (size - 1)) == 0);
1.97 +
1.98 + return size;
1.99 +}
1.100 +
1.101 +nsReturnRef<HRTFKernel> HRTFElevation::calculateKernelForAzimuthElevation(int azimuth, int elevation, SpeexResamplerState* resampler, float sampleRate)
1.102 +{
1.103 + int elevationIndex = (elevation - firstElevation) / elevationSpacing;
1.104 + MOZ_ASSERT(elevationIndex >= 0 && elevationIndex <= numberOfElevations);
1.105 +
1.106 + int numberOfAzimuths = irc_composite_c_r0195[elevationIndex].count;
1.107 + int azimuthSpacing = 360 / numberOfAzimuths;
1.108 + MOZ_ASSERT(numberOfAzimuths * azimuthSpacing == 360);
1.109 +
1.110 + int azimuthIndex = azimuth / azimuthSpacing;
1.111 + MOZ_ASSERT(azimuthIndex * azimuthSpacing == azimuth);
1.112 +
1.113 + const int16_t (&impulse_response_data)[ResponseFrameSize] =
1.114 + irc_composite_c_r0195[elevationIndex].azimuths[azimuthIndex];
1.115 +
1.116 + // When libspeex_resampler is compiled with FIXED_POINT, samples in
1.117 + // speex_resampler_process_float are rounded directly to int16_t, which
1.118 + // only works well if the floats are in the range +/-32767. On such
1.119 + // platforms it's better to resample before converting to float anyway.
1.120 +#ifdef MOZ_SAMPLE_TYPE_S16
1.121 +# define RESAMPLER_PROCESS speex_resampler_process_int
1.122 + const int16_t* response = impulse_response_data;
1.123 + const int16_t* resampledResponse;
1.124 +#else
1.125 +# define RESAMPLER_PROCESS speex_resampler_process_float
1.126 + float response[ResponseFrameSize];
1.127 + ConvertAudioSamples(impulse_response_data, response, ResponseFrameSize);
1.128 + float* resampledResponse;
1.129 +#endif
1.130 +
1.131 + // Note that depending on the fftSize returned by the panner, we may be truncating the impulse response.
1.132 + const size_t resampledResponseLength = fftSizeForSampleRate(sampleRate) / 2;
1.133 +
1.134 + nsAutoTArray<AudioDataValue, 2 * ResponseFrameSize> resampled;
1.135 + if (sampleRate == rawSampleRate) {
1.136 + resampledResponse = response;
1.137 + MOZ_ASSERT(resampledResponseLength == ResponseFrameSize);
1.138 + } else {
1.139 + resampled.SetLength(resampledResponseLength);
1.140 + resampledResponse = resampled.Elements();
1.141 + speex_resampler_skip_zeros(resampler);
1.142 +
1.143 + // Feed the input buffer into the resampler.
1.144 + spx_uint32_t in_len = ResponseFrameSize;
1.145 + spx_uint32_t out_len = resampled.Length();
1.146 + RESAMPLER_PROCESS(resampler, 0, response, &in_len,
1.147 + resampled.Elements(), &out_len);
1.148 +
1.149 + if (out_len < resampled.Length()) {
1.150 + // The input should have all been processed.
1.151 + MOZ_ASSERT(in_len == ResponseFrameSize);
1.152 + // Feed in zeros get the data remaining in the resampler.
1.153 + spx_uint32_t out_index = out_len;
1.154 + in_len = speex_resampler_get_input_latency(resampler);
1.155 + out_len = resampled.Length() - out_index;
1.156 + RESAMPLER_PROCESS(resampler, 0, nullptr, &in_len,
1.157 + resampled.Elements() + out_index, &out_len);
1.158 + out_index += out_len;
1.159 + // There may be some uninitialized samples remaining for very low
1.160 + // sample rates.
1.161 + PodZero(resampled.Elements() + out_index,
1.162 + resampled.Length() - out_index);
1.163 + }
1.164 +
1.165 + speex_resampler_reset_mem(resampler);
1.166 + }
1.167 +
1.168 +#ifdef MOZ_SAMPLE_TYPE_S16
1.169 + nsAutoTArray<float, 2 * ResponseFrameSize> floatArray;
1.170 + floatArray.SetLength(resampledResponseLength);
1.171 + float *floatResponse = floatArray.Elements();
1.172 + ConvertAudioSamples(resampledResponse,
1.173 + floatResponse, resampledResponseLength);
1.174 +#else
1.175 + float *floatResponse = resampledResponse;
1.176 +#endif
1.177 +#undef RESAMPLER_PROCESS
1.178 +
1.179 + return HRTFKernel::create(floatResponse, resampledResponseLength, sampleRate);
1.180 +}
1.181 +
1.182 +// The range of elevations for the IRCAM impulse responses varies depending on azimuth, but the minimum elevation appears to always be -45.
1.183 +//
1.184 +// Here's how it goes:
1.185 +static int maxElevations[] = {
1.186 + // Azimuth
1.187 + //
1.188 + 90, // 0
1.189 + 45, // 15
1.190 + 60, // 30
1.191 + 45, // 45
1.192 + 75, // 60
1.193 + 45, // 75
1.194 + 60, // 90
1.195 + 45, // 105
1.196 + 75, // 120
1.197 + 45, // 135
1.198 + 60, // 150
1.199 + 45, // 165
1.200 + 75, // 180
1.201 + 45, // 195
1.202 + 60, // 210
1.203 + 45, // 225
1.204 + 75, // 240
1.205 + 45, // 255
1.206 + 60, // 270
1.207 + 45, // 285
1.208 + 75, // 300
1.209 + 45, // 315
1.210 + 60, // 330
1.211 + 45 // 345
1.212 +};
1.213 +
1.214 +nsReturnRef<HRTFElevation> HRTFElevation::createBuiltin(int elevation, float sampleRate)
1.215 +{
1.216 + if (elevation < firstElevation ||
1.217 + elevation > firstElevation + numberOfElevations * elevationSpacing ||
1.218 + (elevation / elevationSpacing) * elevationSpacing != elevation)
1.219 + return nsReturnRef<HRTFElevation>();
1.220 +
1.221 + // Spacing, in degrees, between every azimuth loaded from resource.
1.222 + // Some elevations do not have data for all these intervals.
1.223 + // See maxElevations.
1.224 + static const unsigned AzimuthSpacing = 15;
1.225 + static const unsigned NumberOfRawAzimuths = 360 / AzimuthSpacing;
1.226 + static_assert(AzimuthSpacing * NumberOfRawAzimuths == 360,
1.227 + "Not a multiple");
1.228 + static const unsigned InterpolationFactor =
1.229 + NumberOfTotalAzimuths / NumberOfRawAzimuths;
1.230 + static_assert(NumberOfTotalAzimuths ==
1.231 + NumberOfRawAzimuths * InterpolationFactor, "Not a multiple");
1.232 +
1.233 + HRTFKernelList kernelListL;
1.234 + kernelListL.SetLength(NumberOfTotalAzimuths);
1.235 +
1.236 + SpeexResamplerState* resampler = sampleRate == rawSampleRate ? nullptr :
1.237 + speex_resampler_init(1, rawSampleRate, sampleRate,
1.238 + SPEEX_RESAMPLER_QUALITY_DEFAULT, nullptr);
1.239 +
1.240 + // Load convolution kernels from HRTF files.
1.241 + int interpolatedIndex = 0;
1.242 + for (unsigned rawIndex = 0; rawIndex < NumberOfRawAzimuths; ++rawIndex) {
1.243 + // Don't let elevation exceed maximum for this azimuth.
1.244 + int maxElevation = maxElevations[rawIndex];
1.245 + int actualElevation = min(elevation, maxElevation);
1.246 +
1.247 + kernelListL[interpolatedIndex] = calculateKernelForAzimuthElevation(rawIndex * AzimuthSpacing, actualElevation, resampler, sampleRate);
1.248 +
1.249 + interpolatedIndex += InterpolationFactor;
1.250 + }
1.251 +
1.252 + if (resampler)
1.253 + speex_resampler_destroy(resampler);
1.254 +
1.255 + // Now go back and interpolate intermediate azimuth values.
1.256 + for (unsigned i = 0; i < NumberOfTotalAzimuths; i += InterpolationFactor) {
1.257 + int j = (i + InterpolationFactor) % NumberOfTotalAzimuths;
1.258 +
1.259 + // Create the interpolated convolution kernels and delays.
1.260 + for (unsigned jj = 1; jj < InterpolationFactor; ++jj) {
1.261 + float x = float(jj) / float(InterpolationFactor); // interpolate from 0 -> 1
1.262 +
1.263 + kernelListL[i + jj] = HRTFKernel::createInterpolatedKernel(kernelListL[i], kernelListL[j], x);
1.264 + }
1.265 + }
1.266 +
1.267 + return nsReturnRef<HRTFElevation>(new HRTFElevation(&kernelListL, elevation, sampleRate));
1.268 +}
1.269 +
1.270 +nsReturnRef<HRTFElevation> HRTFElevation::createByInterpolatingSlices(HRTFElevation* hrtfElevation1, HRTFElevation* hrtfElevation2, float x, float sampleRate)
1.271 +{
1.272 + MOZ_ASSERT(hrtfElevation1 && hrtfElevation2);
1.273 + if (!hrtfElevation1 || !hrtfElevation2)
1.274 + return nsReturnRef<HRTFElevation>();
1.275 +
1.276 + MOZ_ASSERT(x >= 0.0 && x < 1.0);
1.277 +
1.278 + HRTFKernelList kernelListL;
1.279 + kernelListL.SetLength(NumberOfTotalAzimuths);
1.280 +
1.281 + const HRTFKernelList& kernelListL1 = hrtfElevation1->kernelListL();
1.282 + const HRTFKernelList& kernelListL2 = hrtfElevation2->kernelListL();
1.283 +
1.284 + // Interpolate kernels of corresponding azimuths of the two elevations.
1.285 + for (unsigned i = 0; i < NumberOfTotalAzimuths; ++i) {
1.286 + kernelListL[i] = HRTFKernel::createInterpolatedKernel(kernelListL1[i], kernelListL2[i], x);
1.287 + }
1.288 +
1.289 + // Interpolate elevation angle.
1.290 + double angle = (1.0 - x) * hrtfElevation1->elevationAngle() + x * hrtfElevation2->elevationAngle();
1.291 +
1.292 + return nsReturnRef<HRTFElevation>(new HRTFElevation(&kernelListL, static_cast<int>(angle), sampleRate));
1.293 +}
1.294 +
1.295 +void HRTFElevation::getKernelsFromAzimuth(double azimuthBlend, unsigned azimuthIndex, HRTFKernel* &kernelL, HRTFKernel* &kernelR, double& frameDelayL, double& frameDelayR)
1.296 +{
1.297 + bool checkAzimuthBlend = azimuthBlend >= 0.0 && azimuthBlend < 1.0;
1.298 + MOZ_ASSERT(checkAzimuthBlend);
1.299 + if (!checkAzimuthBlend)
1.300 + azimuthBlend = 0.0;
1.301 +
1.302 + unsigned numKernels = m_kernelListL.Length();
1.303 +
1.304 + bool isIndexGood = azimuthIndex < numKernels;
1.305 + MOZ_ASSERT(isIndexGood);
1.306 + if (!isIndexGood) {
1.307 + kernelL = 0;
1.308 + kernelR = 0;
1.309 + return;
1.310 + }
1.311 +
1.312 + // Return the left and right kernels,
1.313 + // using symmetry to produce the right kernel.
1.314 + kernelL = m_kernelListL[azimuthIndex];
1.315 + int azimuthIndexR = (numKernels - azimuthIndex) % numKernels;
1.316 + kernelR = m_kernelListL[azimuthIndexR];
1.317 +
1.318 + frameDelayL = kernelL->frameDelay();
1.319 + frameDelayR = kernelR->frameDelay();
1.320 +
1.321 + int azimuthIndex2L = (azimuthIndex + 1) % numKernels;
1.322 + double frameDelay2L = m_kernelListL[azimuthIndex2L]->frameDelay();
1.323 + int azimuthIndex2R = (numKernels - azimuthIndex2L) % numKernels;
1.324 + double frameDelay2R = m_kernelListL[azimuthIndex2R]->frameDelay();
1.325 +
1.326 + // Linearly interpolate delays.
1.327 + frameDelayL = (1.0 - azimuthBlend) * frameDelayL + azimuthBlend * frameDelay2L;
1.328 + frameDelayR = (1.0 - azimuthBlend) * frameDelayR + azimuthBlend * frameDelay2R;
1.329 +}
1.330 +
1.331 +} // namespace WebCore
