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 "config.h"
30
31#if ENABLE(WEB_AUDIO)
32
33#include "HRTFElevation.h"
34
35#include "AudioBus.h"
36#include "AudioFileReader.h"
37#include "Biquad.h"
38#include "FFTFrame.h"
39#include "HRTFDatabaseLoader.h"
40#include "HRTFPanner.h"
41#include <algorithm>
42#include <math.h>
43#include <wtf/OwnPtr.h>
44
45using namespace std;
46
47namespace WebCore {
48
49const unsigned HRTFElevation::AzimuthSpacing = 15;
50const unsigned HRTFElevation::NumberOfRawAzimuths = 360 / AzimuthSpacing;
51const unsigned HRTFElevation::InterpolationFactor = 8;
52const unsigned HRTFElevation::NumberOfTotalAzimuths = NumberOfRawAzimuths * InterpolationFactor;
53
54// Total number of components of an HRTF database.
55const size_t TotalNumberOfResponses = 240;
56
57// Number of frames in an individual impulse response.
58const size_t ResponseFrameSize = 256;
59
60// Sample-rate of the spatialization impulse responses as stored in the resource file.
61// The impulse responses may be resampled to a different sample-rate (depending on the audio hardware) when they are loaded.
62const float ResponseSampleRate = 44100;
63
64#if PLATFORM(MAC) || USE(WEBAUDIO_GSTREAMER)
65#define USE_CONCATENATED_IMPULSE_RESPONSES
66#endif
67
68#ifdef USE_CONCATENATED_IMPULSE_RESPONSES
69// Lazily load a concatenated HRTF database for given subject and store it in a
70// local hash table to ensure quick efficient future retrievals.
71static AudioBus* getConcatenatedImpulseResponsesForSubject(const String& subjectName)
72{
73    typedef HashMap<String, AudioBus*> AudioBusMap;
74    DEFINE_STATIC_LOCAL(AudioBusMap, audioBusMap, ());
75
76    AudioBus* bus;
77    AudioBusMap::iterator iterator = audioBusMap.find(subjectName);
78    if (iterator == audioBusMap.end()) {
79        RefPtr<AudioBus> concatenatedImpulseResponses = AudioBus::loadPlatformResource(subjectName.utf8().data(), ResponseSampleRate);
80        ASSERT(concatenatedImpulseResponses);
81        if (!concatenatedImpulseResponses)
82            return 0;
83
84        bus = concatenatedImpulseResponses.release().leakRef();
85        audioBusMap.set(subjectName, bus);
86    } else
87        bus = iterator->value;
88
89    size_t responseLength = bus->length();
90    size_t expectedLength = static_cast<size_t>(TotalNumberOfResponses * ResponseFrameSize);
91
92    // Check number of channels and length. For now these are fixed and known.
93    bool isBusGood = responseLength == expectedLength && bus->numberOfChannels() == 2;
94    ASSERT(isBusGood);
95    if (!isBusGood)
96        return 0;
97
98    return bus;
99}
100#endif
101
102// Takes advantage of the symmetry and creates a composite version of the two measured versions.  For example, we have both azimuth 30 and -30 degrees
103// where the roles of left and right ears are reversed with respect to each other.
104bool HRTFElevation::calculateSymmetricKernelsForAzimuthElevation(int azimuth, int elevation, float sampleRate, const String& subjectName,
105                                                                 RefPtr<HRTFKernel>& kernelL, RefPtr<HRTFKernel>& kernelR)
106{
107    RefPtr<HRTFKernel> kernelL1;
108    RefPtr<HRTFKernel> kernelR1;
109    bool success = calculateKernelsForAzimuthElevation(azimuth, elevation, sampleRate, subjectName, kernelL1, kernelR1);
110    if (!success)
111        return false;
112
113    // And symmetric version
114    int symmetricAzimuth = !azimuth ? 0 : 360 - azimuth;
115
116    RefPtr<HRTFKernel> kernelL2;
117    RefPtr<HRTFKernel> kernelR2;
118    success = calculateKernelsForAzimuthElevation(symmetricAzimuth, elevation, sampleRate, subjectName, kernelL2, kernelR2);
119    if (!success)
120        return false;
121
122    // Notice L/R reversal in symmetric version.
123    kernelL = HRTFKernel::createInterpolatedKernel(kernelL1.get(), kernelR2.get(), 0.5f);
124    kernelR = HRTFKernel::createInterpolatedKernel(kernelR1.get(), kernelL2.get(), 0.5f);
125
126    return true;
127}
128
129bool HRTFElevation::calculateKernelsForAzimuthElevation(int azimuth, int elevation, float sampleRate, const String& subjectName,
130                                                        RefPtr<HRTFKernel>& kernelL, RefPtr<HRTFKernel>& kernelR)
131{
132    // Valid values for azimuth are 0 -> 345 in 15 degree increments.
133    // Valid values for elevation are -45 -> +90 in 15 degree increments.
134
135    bool isAzimuthGood = azimuth >= 0 && azimuth <= 345 && (azimuth / 15) * 15 == azimuth;
136    ASSERT(isAzimuthGood);
137    if (!isAzimuthGood)
138        return false;
139
140    bool isElevationGood = elevation >= -45 && elevation <= 90 && (elevation / 15) * 15 == elevation;
141    ASSERT(isElevationGood);
142    if (!isElevationGood)
143        return false;
144
145    // Construct the resource name from the subject name, azimuth, and elevation, for example:
146    // "IRC_Composite_C_R0195_T015_P000"
147    // Note: the passed in subjectName is not a string passed in via JavaScript or the web.
148    // It's passed in as an internal ASCII identifier and is an implementation detail.
149    int positiveElevation = elevation < 0 ? elevation + 360 : elevation;
150
151#ifdef USE_CONCATENATED_IMPULSE_RESPONSES
152    AudioBus* bus(getConcatenatedImpulseResponsesForSubject(subjectName));
153
154    if (!bus)
155        return false;
156
157    int elevationIndex = positiveElevation / AzimuthSpacing;
158    if (positiveElevation > 90)
159        elevationIndex -= AzimuthSpacing;
160
161    // The concatenated impulse response is a bus containing all
162    // the elevations per azimuth, for all azimuths by increasing
163    // order. So for a given azimuth and elevation we need to compute
164    // the index of the wanted audio frames in the concatenated table.
165    unsigned index = ((azimuth / AzimuthSpacing) * HRTFDatabase::NumberOfRawElevations) + elevationIndex;
166    bool isIndexGood = index < TotalNumberOfResponses;
167    ASSERT(isIndexGood);
168    if (!isIndexGood)
169        return false;
170
171    // Extract the individual impulse response from the concatenated
172    // responses and potentially sample-rate convert it to the desired
173    // (hardware) sample-rate.
174    unsigned startFrame = index * ResponseFrameSize;
175    unsigned stopFrame = startFrame + ResponseFrameSize;
176    RefPtr<AudioBus> preSampleRateConvertedResponse = AudioBus::createBufferFromRange(bus, startFrame, stopFrame);
177    RefPtr<AudioBus> response = AudioBus::createBySampleRateConverting(preSampleRateConvertedResponse.get(), false, sampleRate);
178    AudioChannel* leftEarImpulseResponse = response->channel(AudioBus::ChannelLeft);
179    AudioChannel* rightEarImpulseResponse = response->channel(AudioBus::ChannelRight);
180#else
181    String resourceName = String::format("IRC_%s_C_R0195_T%03d_P%03d", subjectName.utf8().data(), azimuth, positiveElevation);
182
183    RefPtr<AudioBus> impulseResponse(AudioBus::loadPlatformResource(resourceName.utf8().data(), sampleRate));
184
185    ASSERT(impulseResponse.get());
186    if (!impulseResponse.get())
187        return false;
188
189    size_t responseLength = impulseResponse->length();
190    size_t expectedLength = static_cast<size_t>(256 * (sampleRate / 44100.0));
191
192    // Check number of channels and length.  For now these are fixed and known.
193    bool isBusGood = responseLength == expectedLength && impulseResponse->numberOfChannels() == 2;
194    ASSERT(isBusGood);
195    if (!isBusGood)
196        return false;
197
198    AudioChannel* leftEarImpulseResponse = impulseResponse->channelByType(AudioBus::ChannelLeft);
199    AudioChannel* rightEarImpulseResponse = impulseResponse->channelByType(AudioBus::ChannelRight);
200#endif
201
202    // Note that depending on the fftSize returned by the panner, we may be truncating the impulse response we just loaded in.
203    const size_t fftSize = HRTFPanner::fftSizeForSampleRate(sampleRate);
204    kernelL = HRTFKernel::create(leftEarImpulseResponse, fftSize, sampleRate);
205    kernelR = HRTFKernel::create(rightEarImpulseResponse, fftSize, sampleRate);
206
207    return true;
208}
209
210// The range of elevations for the IRCAM impulse responses varies depending on azimuth, but the minimum elevation appears to always be -45.
211//
212// Here's how it goes:
213static int maxElevations[] = {
214        //  Azimuth
215        //
216    90, // 0
217    45, // 15
218    60, // 30
219    45, // 45
220    75, // 60
221    45, // 75
222    60, // 90
223    45, // 105
224    75, // 120
225    45, // 135
226    60, // 150
227    45, // 165
228    75, // 180
229    45, // 195
230    60, // 210
231    45, // 225
232    75, // 240
233    45, // 255
234    60, // 270
235    45, // 285
236    75, // 300
237    45, // 315
238    60, // 330
239    45 //  345
240};
241
242PassOwnPtr<HRTFElevation> HRTFElevation::createForSubject(const String& subjectName, int elevation, float sampleRate)
243{
244    bool isElevationGood = elevation >= -45 && elevation <= 90 && (elevation / 15) * 15 == elevation;
245    ASSERT(isElevationGood);
246    if (!isElevationGood)
247        return nullptr;
248
249    OwnPtr<HRTFKernelList> kernelListL = adoptPtr(new HRTFKernelList(NumberOfTotalAzimuths));
250    OwnPtr<HRTFKernelList> kernelListR = adoptPtr(new HRTFKernelList(NumberOfTotalAzimuths));
251
252    // Load convolution kernels from HRTF files.
253    int interpolatedIndex = 0;
254    for (unsigned rawIndex = 0; rawIndex < NumberOfRawAzimuths; ++rawIndex) {
255        // Don't let elevation exceed maximum for this azimuth.
256        int maxElevation = maxElevations[rawIndex];
257        int actualElevation = min(elevation, maxElevation);
258
259        bool success = calculateKernelsForAzimuthElevation(rawIndex * AzimuthSpacing, actualElevation, sampleRate, subjectName, kernelListL->at(interpolatedIndex), kernelListR->at(interpolatedIndex));
260        if (!success)
261            return nullptr;
262
263        interpolatedIndex += InterpolationFactor;
264    }
265
266    // Now go back and interpolate intermediate azimuth values.
267    for (unsigned i = 0; i < NumberOfTotalAzimuths; i += InterpolationFactor) {
268        int j = (i + InterpolationFactor) % NumberOfTotalAzimuths;
269
270        // Create the interpolated convolution kernels and delays.
271        for (unsigned jj = 1; jj < InterpolationFactor; ++jj) {
272            float x = float(jj) / float(InterpolationFactor); // interpolate from 0 -> 1
273
274            (*kernelListL)[i + jj] = HRTFKernel::createInterpolatedKernel(kernelListL->at(i).get(), kernelListL->at(j).get(), x);
275            (*kernelListR)[i + jj] = HRTFKernel::createInterpolatedKernel(kernelListR->at(i).get(), kernelListR->at(j).get(), x);
276        }
277    }
278
279    OwnPtr<HRTFElevation> hrtfElevation = adoptPtr(new HRTFElevation(kernelListL.release(), kernelListR.release(), elevation, sampleRate));
280    return hrtfElevation.release();
281}
282
283PassOwnPtr<HRTFElevation> HRTFElevation::createByInterpolatingSlices(HRTFElevation* hrtfElevation1, HRTFElevation* hrtfElevation2, float x, float sampleRate)
284{
285    ASSERT(hrtfElevation1 && hrtfElevation2);
286    if (!hrtfElevation1 || !hrtfElevation2)
287        return nullptr;
288
289    ASSERT(x >= 0.0 && x < 1.0);
290
291    OwnPtr<HRTFKernelList> kernelListL = adoptPtr(new HRTFKernelList(NumberOfTotalAzimuths));
292    OwnPtr<HRTFKernelList> kernelListR = adoptPtr(new HRTFKernelList(NumberOfTotalAzimuths));
293
294    HRTFKernelList* kernelListL1 = hrtfElevation1->kernelListL();
295    HRTFKernelList* kernelListR1 = hrtfElevation1->kernelListR();
296    HRTFKernelList* kernelListL2 = hrtfElevation2->kernelListL();
297    HRTFKernelList* kernelListR2 = hrtfElevation2->kernelListR();
298
299    // Interpolate kernels of corresponding azimuths of the two elevations.
300    for (unsigned i = 0; i < NumberOfTotalAzimuths; ++i) {
301        (*kernelListL)[i] = HRTFKernel::createInterpolatedKernel(kernelListL1->at(i).get(), kernelListL2->at(i).get(), x);
302        (*kernelListR)[i] = HRTFKernel::createInterpolatedKernel(kernelListR1->at(i).get(), kernelListR2->at(i).get(), x);
303    }
304
305    // Interpolate elevation angle.
306    double angle = (1.0 - x) * hrtfElevation1->elevationAngle() + x * hrtfElevation2->elevationAngle();
307
308    OwnPtr<HRTFElevation> hrtfElevation = adoptPtr(new HRTFElevation(kernelListL.release(), kernelListR.release(), static_cast<int>(angle), sampleRate));
309    return hrtfElevation.release();
310}
311
312void HRTFElevation::getKernelsFromAzimuth(double azimuthBlend, unsigned azimuthIndex, HRTFKernel* &kernelL, HRTFKernel* &kernelR, double& frameDelayL, double& frameDelayR)
313{
314    bool checkAzimuthBlend = azimuthBlend >= 0.0 && azimuthBlend < 1.0;
315    ASSERT(checkAzimuthBlend);
316    if (!checkAzimuthBlend)
317        azimuthBlend = 0.0;
318
319    unsigned numKernels = m_kernelListL->size();
320
321    bool isIndexGood = azimuthIndex < numKernels;
322    ASSERT(isIndexGood);
323    if (!isIndexGood) {
324        kernelL = 0;
325        kernelR = 0;
326        return;
327    }
328
329    // Return the left and right kernels.
330    kernelL = m_kernelListL->at(azimuthIndex).get();
331    kernelR = m_kernelListR->at(azimuthIndex).get();
332
333    frameDelayL = m_kernelListL->at(azimuthIndex)->frameDelay();
334    frameDelayR = m_kernelListR->at(azimuthIndex)->frameDelay();
335
336    int azimuthIndex2 = (azimuthIndex + 1) % numKernels;
337    double frameDelay2L = m_kernelListL->at(azimuthIndex2)->frameDelay();
338    double frameDelay2R = m_kernelListR->at(azimuthIndex2)->frameDelay();
339
340    // Linearly interpolate delays.
341    frameDelayL = (1.0 - azimuthBlend) * frameDelayL + azimuthBlend * frameDelay2L;
342    frameDelayR = (1.0 - azimuthBlend) * frameDelayR + azimuthBlend * frameDelay2R;
343}
344
345} // namespace WebCore
346
347#endif // ENABLE(WEB_AUDIO)
348