content/media/webaudio/blink/HRTFElevation.cpp

Thu, 22 Jan 2015 13:21:57 +0100

author
Michael Schloh von Bennewitz <michael@schloh.com>
date
Thu, 22 Jan 2015 13:21:57 +0100
branch
TOR_BUG_9701
changeset 15
b8a032363ba2
permissions
-rw-r--r--

Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6

     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  */
    29 #include "HRTFElevation.h"
    31 #include "speex/speex_resampler.h"
    32 #include "mozilla/PodOperations.h"
    33 #include "AudioSampleFormat.h"
    35 #include "IRC_Composite_C_R0195-incl.cpp"
    37 using namespace std;
    38 using namespace mozilla;
    40 namespace WebCore {
    42 const int elevationSpacing = irc_composite_c_r0195_elevation_interval;
    43 const int firstElevation = irc_composite_c_r0195_first_elevation;
    44 const int numberOfElevations = MOZ_ARRAY_LENGTH(irc_composite_c_r0195);
    46 const unsigned HRTFElevation::NumberOfTotalAzimuths = 360 / 15 * 8;
    48 const int rawSampleRate = irc_composite_c_r0195_sample_rate;
    50 // Number of frames in an individual impulse response.
    51 const size_t ResponseFrameSize = 256;
    53 size_t HRTFElevation::sizeOfIncludingThis(mozilla::MallocSizeOf aMallocSizeOf) const
    54 {
    55     size_t amount = aMallocSizeOf(this);
    57     amount += m_kernelListL.SizeOfExcludingThis(aMallocSizeOf);
    58     for (size_t i = 0; i < m_kernelListL.Length(); i++) {
    59         amount += m_kernelListL[i]->sizeOfIncludingThis(aMallocSizeOf);
    60     }
    62     return amount;
    63 }
    65 size_t HRTFElevation::fftSizeForSampleRate(float sampleRate)
    66 {
    67     // The IRCAM HRTF impulse responses were 512 sample-frames @44.1KHz,
    68     // but these have been truncated to 256 samples.
    69     // An FFT-size of twice impulse response size is used (for convolution).
    70     // So for sample rates of 44.1KHz an FFT size of 512 is good.
    71     // We double the FFT-size only for sample rates at least double this.
    72     // If the FFT size is too large then the impulse response will be padded
    73     // with zeros without the fade-out provided by HRTFKernel.
    74     MOZ_ASSERT(sampleRate > 1.0 && sampleRate < 1048576.0);
    76     // This is the size if we were to use all raw response samples.
    77     unsigned resampledLength =
    78         floorf(ResponseFrameSize * sampleRate / rawSampleRate);
    79     // Keep things semi-sane, with max FFT size of 1024 and minimum of 4.
    80     // "size |= 3" ensures a minimum of 4 (with the size++ below) and sets the
    81     // 2 least significant bits for rounding up to the next power of 2 below.
    82     unsigned size = min(resampledLength, 1023U);
    83     size |= 3;
    84     // Round up to the next power of 2, making the FFT size no more than twice
    85     // the impulse response length.  This doubles size for values that are
    86     // already powers of 2.  This works by filling in 7 bits to right of the
    87     // most significant bit.  The most significant bit is no greater than
    88     // 1 << 9, and the least significant 2 bits were already set above.
    89     size |= (size >> 1);
    90     size |= (size >> 2);
    91     size |= (size >> 4);
    92     size++;
    93     MOZ_ASSERT((size & (size - 1)) == 0);
    95     return size;
    96 }
    98 nsReturnRef<HRTFKernel> HRTFElevation::calculateKernelForAzimuthElevation(int azimuth, int elevation, SpeexResamplerState* resampler, float sampleRate)
    99 {
   100     int elevationIndex = (elevation - firstElevation) / elevationSpacing;
   101     MOZ_ASSERT(elevationIndex >= 0 && elevationIndex <= numberOfElevations);
   103     int numberOfAzimuths = irc_composite_c_r0195[elevationIndex].count;
   104     int azimuthSpacing = 360 / numberOfAzimuths;
   105     MOZ_ASSERT(numberOfAzimuths * azimuthSpacing == 360);
   107     int azimuthIndex = azimuth / azimuthSpacing;
   108     MOZ_ASSERT(azimuthIndex * azimuthSpacing == azimuth);
   110     const int16_t (&impulse_response_data)[ResponseFrameSize] =
   111         irc_composite_c_r0195[elevationIndex].azimuths[azimuthIndex];
   113     // When libspeex_resampler is compiled with FIXED_POINT, samples in
   114     // speex_resampler_process_float are rounded directly to int16_t, which
   115     // only works well if the floats are in the range +/-32767.  On such
   116     // platforms it's better to resample before converting to float anyway.
   117 #ifdef MOZ_SAMPLE_TYPE_S16
   118 #  define RESAMPLER_PROCESS speex_resampler_process_int
   119     const int16_t* response = impulse_response_data;
   120     const int16_t* resampledResponse;
   121 #else
   122 #  define RESAMPLER_PROCESS speex_resampler_process_float
   123     float response[ResponseFrameSize];
   124     ConvertAudioSamples(impulse_response_data, response, ResponseFrameSize);
   125     float* resampledResponse;
   126 #endif
   128     // Note that depending on the fftSize returned by the panner, we may be truncating the impulse response.
   129     const size_t resampledResponseLength = fftSizeForSampleRate(sampleRate) / 2;
   131     nsAutoTArray<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);
   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,
   144                           resampled.Elements(), &out_len);
   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,
   159                     resampled.Length() - out_index);
   160         }
   162         speex_resampler_reset_mem(resampler);
   163     }
   165 #ifdef MOZ_SAMPLE_TYPE_S16
   166     nsAutoTArray<float, 2 * ResponseFrameSize> floatArray;
   167     floatArray.SetLength(resampledResponseLength);
   168     float *floatResponse = floatArray.Elements();
   169     ConvertAudioSamples(resampledResponse,
   170                         floatResponse, resampledResponseLength);
   171 #else
   172     float *floatResponse = resampledResponse;
   173 #endif
   174 #undef RESAMPLER_PROCESS
   176     return HRTFKernel::create(floatResponse, resampledResponseLength, sampleRate);
   177 }
   179 // The range of elevations for the IRCAM impulse responses varies depending on 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 };
   211 nsReturnRef<HRTFElevation> HRTFElevation::createBuiltin(int elevation, float sampleRate)
   212 {
   213     if (elevation < firstElevation ||
   214         elevation > firstElevation + numberOfElevations * elevationSpacing ||
   215         (elevation / elevationSpacing) * elevationSpacing != elevation)
   216         return nsReturnRef<HRTFElevation>();
   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,
   224                   "Not a multiple");
   225     static const unsigned InterpolationFactor =
   226         NumberOfTotalAzimuths / NumberOfRawAzimuths;
   227     static_assert(NumberOfTotalAzimuths ==
   228                   NumberOfRawAzimuths * InterpolationFactor, "Not a multiple");
   230     HRTFKernelList kernelListL;
   231     kernelListL.SetLength(NumberOfTotalAzimuths);
   233     SpeexResamplerState* resampler = sampleRate == rawSampleRate ? nullptr :
   234         speex_resampler_init(1, rawSampleRate, sampleRate,
   235                              SPEEX_RESAMPLER_QUALITY_DEFAULT, nullptr);
   237     // Load convolution kernels from HRTF files.
   238     int interpolatedIndex = 0;
   239     for (unsigned rawIndex = 0; rawIndex < NumberOfRawAzimuths; ++rawIndex) {
   240         // Don't let elevation exceed maximum for this azimuth.
   241         int maxElevation = maxElevations[rawIndex];
   242         int actualElevation = min(elevation, maxElevation);
   244         kernelListL[interpolatedIndex] = calculateKernelForAzimuthElevation(rawIndex * AzimuthSpacing, actualElevation, resampler, sampleRate);
   246         interpolatedIndex += InterpolationFactor;
   247     }
   249     if (resampler)
   250         speex_resampler_destroy(resampler);
   252     // Now go back and interpolate intermediate azimuth values.
   253     for (unsigned i = 0; i < NumberOfTotalAzimuths; i += InterpolationFactor) {
   254         int j = (i + InterpolationFactor) % NumberOfTotalAzimuths;
   256         // Create the interpolated convolution kernels and delays.
   257         for (unsigned jj = 1; jj < InterpolationFactor; ++jj) {
   258             float x = float(jj) / float(InterpolationFactor); // interpolate from 0 -> 1
   260             kernelListL[i + jj] = HRTFKernel::createInterpolatedKernel(kernelListL[i], kernelListL[j], x);
   261         }
   262     }
   264     return nsReturnRef<HRTFElevation>(new HRTFElevation(&kernelListL, elevation, sampleRate));
   265 }
   267 nsReturnRef<HRTFElevation> HRTFElevation::createByInterpolatingSlices(HRTFElevation* hrtfElevation1, HRTFElevation* hrtfElevation2, float x, float sampleRate)
   268 {
   269     MOZ_ASSERT(hrtfElevation1 && hrtfElevation2);
   270     if (!hrtfElevation1 || !hrtfElevation2)
   271         return nsReturnRef<HRTFElevation>();
   273     MOZ_ASSERT(x >= 0.0 && x < 1.0);
   275     HRTFKernelList kernelListL;
   276     kernelListL.SetLength(NumberOfTotalAzimuths);
   278     const HRTFKernelList& kernelListL1 = hrtfElevation1->kernelListL();
   279     const HRTFKernelList& kernelListL2 = hrtfElevation2->kernelListL();
   281     // Interpolate kernels of corresponding azimuths of the two elevations.
   282     for (unsigned i = 0; i < NumberOfTotalAzimuths; ++i) {
   283         kernelListL[i] = HRTFKernel::createInterpolatedKernel(kernelListL1[i], kernelListL2[i], x);
   284     }
   286     // Interpolate elevation angle.
   287     double angle = (1.0 - x) * hrtfElevation1->elevationAngle() + x * hrtfElevation2->elevationAngle();
   289     return nsReturnRef<HRTFElevation>(new HRTFElevation(&kernelListL, static_cast<int>(angle), sampleRate));
   290 }
   292 void HRTFElevation::getKernelsFromAzimuth(double azimuthBlend, unsigned azimuthIndex, HRTFKernel* &kernelL, HRTFKernel* &kernelR, double& frameDelayL, double& frameDelayR)
   293 {
   294     bool checkAzimuthBlend = azimuthBlend >= 0.0 && azimuthBlend < 1.0;
   295     MOZ_ASSERT(checkAzimuthBlend);
   296     if (!checkAzimuthBlend)
   297         azimuthBlend = 0.0;
   299     unsigned numKernels = m_kernelListL.Length();
   301     bool isIndexGood = azimuthIndex < numKernels;
   302     MOZ_ASSERT(isIndexGood);
   303     if (!isIndexGood) {
   304         kernelL = 0;
   305         kernelR = 0;
   306         return;
   307     }
   309     // Return the left and right kernels,
   310     // using symmetry to produce the right kernel.
   311     kernelL = m_kernelListL[azimuthIndex];
   312     int azimuthIndexR = (numKernels - azimuthIndex) % numKernels;
   313     kernelR = m_kernelListL[azimuthIndexR];
   315     frameDelayL = kernelL->frameDelay();
   316     frameDelayR = kernelR->frameDelay();
   318     int azimuthIndex2L = (azimuthIndex + 1) % numKernels;
   319     double frameDelay2L = m_kernelListL[azimuthIndex2L]->frameDelay();
   320     int azimuthIndex2R = (numKernels - azimuthIndex2L) % numKernels;
   321     double frameDelay2R = m_kernelListL[azimuthIndex2R]->frameDelay();
   323     // Linearly interpolate delays.
   324     frameDelayL = (1.0 - azimuthBlend) * frameDelayL + azimuthBlend * frameDelay2L;
   325     frameDelayR = (1.0 - azimuthBlend) * frameDelayR + azimuthBlend * frameDelay2R;
   326 }
   328 } // namespace WebCore

mercurial