1 /*
2 * Copyright (C) 2010 Google Inc. All rights reserved.
3 *
4 * Redistribution and use in source and binary forms, with or without
5 * modification, are permitted provided that the following conditions
6 * are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
13 * 3. Neither the name of Apple Computer, Inc. ("Apple") nor the names of
14 * its contributors may be used to endorse or promote products derived
15 * from this software without specific prior written permission.
16 *
17 * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
18 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
19 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
20 * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
21 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
22 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
23 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
24 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
25 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
26 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29 #include "HRTFElevation.h"
30
31 #include <speex/speex_resampler.h>
32 #include "mozilla/PodOperations.h"
33 #include "AudioSampleFormat.h"
34
35 #include "IRC_Composite_C_R0195-incl.cpp"
36
37 using namespace mozilla;
38
39 namespace WebCore {
40
41 const int elevationSpacing = irc_composite_c_r0195_elevation_interval;
42 const int firstElevation = irc_composite_c_r0195_first_elevation;
43 const int numberOfElevations = MOZ_ARRAY_LENGTH(irc_composite_c_r0195);
44
45 const unsigned HRTFElevation::NumberOfTotalAzimuths = 360 / 15 * 8;
46
47 const int rawSampleRate = irc_composite_c_r0195_sample_rate;
48
49 // Number of frames in an individual impulse response.
50 const size_t ResponseFrameSize = 256;
51
sizeOfIncludingThis(mozilla::MallocSizeOf aMallocSizeOf) const52 size_t HRTFElevation::sizeOfIncludingThis(
53 mozilla::MallocSizeOf aMallocSizeOf) const {
54 size_t amount = aMallocSizeOf(this);
55
56 amount += m_kernelListL.ShallowSizeOfExcludingThis(aMallocSizeOf);
57 for (size_t i = 0; i < m_kernelListL.Length(); i++) {
58 amount += m_kernelListL[i]->sizeOfIncludingThis(aMallocSizeOf);
59 }
60
61 return amount;
62 }
63
fftSizeForSampleRate(float sampleRate)64 size_t HRTFElevation::fftSizeForSampleRate(float sampleRate) {
65 // The IRCAM HRTF impulse responses were 512 sample-frames @44.1KHz,
66 // but these have been truncated to 256 samples.
67 // An FFT-size of twice impulse response size is used (for convolution).
68 // So for sample rates of 44.1KHz an FFT size of 512 is good.
69 // We double the FFT-size only for sample rates at least double this.
70 // If the FFT size is too large then the impulse response will be padded
71 // with zeros without the fade-out provided by HRTFKernel.
72 MOZ_ASSERT(sampleRate > 1.0 && sampleRate < 1048576.0);
73
74 // This is the size if we were to use all raw response samples.
75 unsigned resampledLength =
76 floorf(ResponseFrameSize * sampleRate / rawSampleRate);
77 // Keep things semi-sane, with max FFT size of 1024.
78 unsigned size = std::min(resampledLength, 1023U);
79 // Ensure a minimum of 2 * WEBAUDIO_BLOCK_SIZE (with the size++ below) for
80 // FFTConvolver and set the 8 least significant bits for rounding up to
81 // the next power of 2 below.
82 size |= 2 * WEBAUDIO_BLOCK_SIZE - 1;
83 // Round up to the next power of 2, making the FFT size no more than twice
84 // the impulse response length. This doubles size for values that are
85 // already powers of 2. This works by filling in alls bit to right of the
86 // most significant bit. The most significant bit is no greater than
87 // 1 << 9, and the least significant 8 bits were already set above, so
88 // there is at most one bit to add.
89 size |= (size >> 1);
90 size++;
91 MOZ_ASSERT((size & (size - 1)) == 0);
92
93 return size;
94 }
95
calculateKernelForAzimuthElevation(int azimuth,int elevation,SpeexResamplerState * resampler,float sampleRate)96 nsReturnRef<HRTFKernel> HRTFElevation::calculateKernelForAzimuthElevation(
97 int azimuth, int elevation, SpeexResamplerState* resampler,
98 float sampleRate) {
99 int elevationIndex = (elevation - firstElevation) / elevationSpacing;
100 MOZ_ASSERT(elevationIndex >= 0 && elevationIndex <= numberOfElevations);
101
102 int numberOfAzimuths = irc_composite_c_r0195[elevationIndex].count;
103 int azimuthSpacing = 360 / numberOfAzimuths;
104 MOZ_ASSERT(numberOfAzimuths * azimuthSpacing == 360);
105
106 int azimuthIndex = azimuth / azimuthSpacing;
107 MOZ_ASSERT(azimuthIndex * azimuthSpacing == azimuth);
108
109 const int16_t(&impulse_response_data)[ResponseFrameSize] =
110 irc_composite_c_r0195[elevationIndex].azimuths[azimuthIndex];
111
112 // When libspeex_resampler is compiled with FIXED_POINT, samples in
113 // speex_resampler_process_float are rounded directly to int16_t, which
114 // only works well if the floats are in the range +/-32767. On such
115 // platforms it's better to resample before converting to float anyway.
116 #ifdef MOZ_SAMPLE_TYPE_S16
117 # define RESAMPLER_PROCESS speex_resampler_process_int
118 const int16_t* response = impulse_response_data;
119 const int16_t* resampledResponse;
120 #else
121 # define RESAMPLER_PROCESS speex_resampler_process_float
122 float response[ResponseFrameSize];
123 ConvertAudioSamples(impulse_response_data, response, ResponseFrameSize);
124 float* resampledResponse;
125 #endif
126
127 // Note that depending on the fftSize returned by the panner, we may be
128 // truncating the impulse response.
129 const size_t resampledResponseLength = fftSizeForSampleRate(sampleRate) / 2;
130
131 AutoTArray<AudioDataValue, 2 * ResponseFrameSize> resampled;
132 if (sampleRate == rawSampleRate) {
133 resampledResponse = response;
134 MOZ_ASSERT(resampledResponseLength == ResponseFrameSize);
135 } else {
136 resampled.SetLength(resampledResponseLength);
137 resampledResponse = resampled.Elements();
138 speex_resampler_skip_zeros(resampler);
139
140 // Feed the input buffer into the resampler.
141 spx_uint32_t in_len = ResponseFrameSize;
142 spx_uint32_t out_len = resampled.Length();
143 RESAMPLER_PROCESS(resampler, 0, response, &in_len, resampled.Elements(),
144 &out_len);
145
146 if (out_len < resampled.Length()) {
147 // The input should have all been processed.
148 MOZ_ASSERT(in_len == ResponseFrameSize);
149 // Feed in zeros get the data remaining in the resampler.
150 spx_uint32_t out_index = out_len;
151 in_len = speex_resampler_get_input_latency(resampler);
152 out_len = resampled.Length() - out_index;
153 RESAMPLER_PROCESS(resampler, 0, nullptr, &in_len,
154 resampled.Elements() + out_index, &out_len);
155 out_index += out_len;
156 // There may be some uninitialized samples remaining for very low
157 // sample rates.
158 PodZero(resampled.Elements() + out_index, resampled.Length() - out_index);
159 }
160
161 speex_resampler_reset_mem(resampler);
162 }
163
164 #ifdef MOZ_SAMPLE_TYPE_S16
165 AutoTArray<float, 2 * ResponseFrameSize> floatArray;
166 floatArray.SetLength(resampledResponseLength);
167 float* floatResponse = floatArray.Elements();
168 ConvertAudioSamples(resampledResponse, floatResponse,
169 resampledResponseLength);
170 #else
171 float* floatResponse = resampledResponse;
172 #endif
173 #undef RESAMPLER_PROCESS
174
175 return HRTFKernel::create(floatResponse, resampledResponseLength, sampleRate);
176 }
177
178 // The range of elevations for the IRCAM impulse responses varies depending on
179 // azimuth, but the minimum elevation appears to always be -45.
180 //
181 // Here's how it goes:
182 static int maxElevations[] = {
183 // Azimuth
184 //
185 90, // 0
186 45, // 15
187 60, // 30
188 45, // 45
189 75, // 60
190 45, // 75
191 60, // 90
192 45, // 105
193 75, // 120
194 45, // 135
195 60, // 150
196 45, // 165
197 75, // 180
198 45, // 195
199 60, // 210
200 45, // 225
201 75, // 240
202 45, // 255
203 60, // 270
204 45, // 285
205 75, // 300
206 45, // 315
207 60, // 330
208 45 // 345
209 };
210
createBuiltin(int elevation,float sampleRate)211 nsReturnRef<HRTFElevation> HRTFElevation::createBuiltin(int elevation,
212 float sampleRate) {
213 if (elevation < firstElevation ||
214 elevation > firstElevation + numberOfElevations * elevationSpacing ||
215 (elevation / elevationSpacing) * elevationSpacing != elevation)
216 return nsReturnRef<HRTFElevation>();
217
218 // Spacing, in degrees, between every azimuth loaded from resource.
219 // Some elevations do not have data for all these intervals.
220 // See maxElevations.
221 static const unsigned AzimuthSpacing = 15;
222 static const unsigned NumberOfRawAzimuths = 360 / AzimuthSpacing;
223 static_assert(AzimuthSpacing * NumberOfRawAzimuths == 360, "Not a multiple");
224 static const unsigned InterpolationFactor =
225 NumberOfTotalAzimuths / NumberOfRawAzimuths;
226 static_assert(
227 NumberOfTotalAzimuths == NumberOfRawAzimuths * InterpolationFactor,
228 "Not a multiple");
229
230 HRTFKernelList kernelListL;
231 kernelListL.SetLength(NumberOfTotalAzimuths);
232
233 SpeexResamplerState* resampler =
234 sampleRate == rawSampleRate
235 ? nullptr
236 : speex_resampler_init(1, rawSampleRate, sampleRate,
237 SPEEX_RESAMPLER_QUALITY_MIN, nullptr);
238
239 // Load convolution kernels from HRTF files.
240 int interpolatedIndex = 0;
241 for (unsigned rawIndex = 0; rawIndex < NumberOfRawAzimuths; ++rawIndex) {
242 // Don't let elevation exceed maximum for this azimuth.
243 int maxElevation = maxElevations[rawIndex];
244 int actualElevation = std::min(elevation, maxElevation);
245
246 kernelListL[interpolatedIndex] = calculateKernelForAzimuthElevation(
247 rawIndex * AzimuthSpacing, actualElevation, resampler, sampleRate);
248
249 interpolatedIndex += InterpolationFactor;
250 }
251
252 if (resampler) speex_resampler_destroy(resampler);
253
254 // Now go back and interpolate intermediate azimuth values.
255 for (unsigned i = 0; i < NumberOfTotalAzimuths; i += InterpolationFactor) {
256 int j = (i + InterpolationFactor) % NumberOfTotalAzimuths;
257
258 // Create the interpolated convolution kernels and delays.
259 for (unsigned jj = 1; jj < InterpolationFactor; ++jj) {
260 float x =
261 float(jj) / float(InterpolationFactor); // interpolate from 0 -> 1
262
263 kernelListL[i + jj] = HRTFKernel::createInterpolatedKernel(
264 kernelListL[i], kernelListL[j], x);
265 }
266 }
267
268 return nsReturnRef<HRTFElevation>(
269 new HRTFElevation(std::move(kernelListL), elevation, sampleRate));
270 }
271
createByInterpolatingSlices(HRTFElevation * hrtfElevation1,HRTFElevation * hrtfElevation2,float x,float sampleRate)272 nsReturnRef<HRTFElevation> HRTFElevation::createByInterpolatingSlices(
273 HRTFElevation* hrtfElevation1, HRTFElevation* hrtfElevation2, float x,
274 float sampleRate) {
275 MOZ_ASSERT(hrtfElevation1 && hrtfElevation2);
276 if (!hrtfElevation1 || !hrtfElevation2) return nsReturnRef<HRTFElevation>();
277
278 MOZ_ASSERT(x >= 0.0 && x < 1.0);
279
280 HRTFKernelList kernelListL;
281 kernelListL.SetLength(NumberOfTotalAzimuths);
282
283 const HRTFKernelList& kernelListL1 = hrtfElevation1->kernelListL();
284 const HRTFKernelList& kernelListL2 = hrtfElevation2->kernelListL();
285
286 // Interpolate kernels of corresponding azimuths of the two elevations.
287 for (unsigned i = 0; i < NumberOfTotalAzimuths; ++i) {
288 kernelListL[i] = HRTFKernel::createInterpolatedKernel(kernelListL1[i],
289 kernelListL2[i], x);
290 }
291
292 // Interpolate elevation angle.
293 double angle = (1.0 - x) * hrtfElevation1->elevationAngle() +
294 x * hrtfElevation2->elevationAngle();
295
296 return nsReturnRef<HRTFElevation>(new HRTFElevation(
297 std::move(kernelListL), static_cast<int>(angle), sampleRate));
298 }
299
getKernelsFromAzimuth(double azimuthBlend,unsigned azimuthIndex,HRTFKernel * & kernelL,HRTFKernel * & kernelR,double & frameDelayL,double & frameDelayR)300 void HRTFElevation::getKernelsFromAzimuth(
301 double azimuthBlend, unsigned azimuthIndex, HRTFKernel*& kernelL,
302 HRTFKernel*& kernelR, double& frameDelayL, double& frameDelayR) {
303 bool checkAzimuthBlend = azimuthBlend >= 0.0 && azimuthBlend < 1.0;
304 MOZ_ASSERT(checkAzimuthBlend);
305 if (!checkAzimuthBlend) azimuthBlend = 0.0;
306
307 unsigned numKernels = m_kernelListL.Length();
308
309 bool isIndexGood = azimuthIndex < numKernels;
310 MOZ_ASSERT(isIndexGood);
311 if (!isIndexGood) {
312 kernelL = 0;
313 kernelR = 0;
314 return;
315 }
316
317 // Return the left and right kernels,
318 // using symmetry to produce the right kernel.
319 kernelL = m_kernelListL[azimuthIndex];
320 int azimuthIndexR = (numKernels - azimuthIndex) % numKernels;
321 kernelR = m_kernelListL[azimuthIndexR];
322
323 frameDelayL = kernelL->frameDelay();
324 frameDelayR = kernelR->frameDelay();
325
326 int azimuthIndex2L = (azimuthIndex + 1) % numKernels;
327 double frameDelay2L = m_kernelListL[azimuthIndex2L]->frameDelay();
328 int azimuthIndex2R = (numKernels - azimuthIndex2L) % numKernels;
329 double frameDelay2R = m_kernelListL[azimuthIndex2R]->frameDelay();
330
331 // Linearly interpolate delays.
332 frameDelayL =
333 (1.0 - azimuthBlend) * frameDelayL + azimuthBlend * frameDelay2L;
334 frameDelayR =
335 (1.0 - azimuthBlend) * frameDelayR + azimuthBlend * frameDelay2R;
336 }
337
338 } // namespace WebCore
339