src/third_party/WebKit/Source/platform/audio/FFTFrame.cpp

   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 "platform/audio/FFTFrame.h"
  34
  35 #ifndef NDEBUG
  36 #include <stdio.h>
  37 #endif
  38
  39 #include "platform/Logging.h"
  40 #include "wtf/Complex.h"
  41 #include "wtf/MathExtras.h"
  42 #include "wtf/OwnPtr.h"
  43
  44 namespace WebCore {
  45
  46 void FFTFrame::doPaddedFFT(const float* data, size_t dataSize)
  47 {
  48     // Zero-pad the impulse response
  49     AudioFloatArray paddedResponse(fftSize()); // zero-initialized
  50     paddedResponse.copyToRange(data, 0, dataSize);
  51
  52     // Get the frequency-domain version of padded response
  53     doFFT(paddedResponse.data());
  54 }
  55
  56 PassOwnPtr<FFTFrame> FFTFrame::createInterpolatedFrame(const FFTFrame& frame1, const FFTFrame& frame2, double x)
  57 {
  58     OwnPtr<FFTFrame> newFrame = adoptPtr(new FFTFrame(frame1.fftSize()));
  59
  60     newFrame->interpolateFrequencyComponents(frame1, frame2, x);
  61
  62     // In the time-domain, the 2nd half of the response must be zero, to avoid circular convolution aliasing...
  63     int fftSize = newFrame->fftSize();
  64     AudioFloatArray buffer(fftSize);
  65     newFrame->doInverseFFT(buffer.data());
  66     buffer.zeroRange(fftSize / 2, fftSize);
  67
  68     // Put back into frequency domain.
  69     newFrame->doFFT(buffer.data());
  70
  71     return newFrame.release();
  72 }
  73
  74 #if OS(WIN)
  75 // On Windows, the following pragmas are equivalent to compiling the code with /fp:fast. The
  76 // following code does not need precise FP semantics, and speed is critical here. See
  77 // crbug.com/316740 and crrev.com/116823002.
  78 #pragma float_control(push)
  79 #pragma float_control(except, off)
  80 #pragma float_control(precise, off)
  81 #pragma fp_contract(on)
  82 #pragma fenv_access(off)
  83 #endif
  84
  85 void FFTFrame::interpolateFrequencyComponents(const FFTFrame& frame1, const FFTFrame& frame2, double interp)
  86 {
  87     // FIXME : with some work, this method could be optimized
  88
  89     float* realP = realData();
  90     float* imagP = imagData();
  91
  92     const float* realP1 = frame1.realData();
  93     const float* imagP1 = frame1.imagData();
  94     const float* realP2 = frame2.realData();
  95     const float* imagP2 = frame2.imagData();
  96
  97     m_FFTSize = frame1.fftSize();
  98     m_log2FFTSize = frame1.log2FFTSize();
  99
 100     double s1base = (1.0 - interp);
 101     double s2base = interp;
 102
 103     double phaseAccum = 0.0;
 104     double lastPhase1 = 0.0;
 105     double lastPhase2 = 0.0;
 106
 107     realP[0] = static_cast<float>(s1base * realP1[0] + s2base * realP2[0]);
 108     imagP[0] = static_cast<float>(s1base * imagP1[0] + s2base * imagP2[0]);
 109
 110     int n = m_FFTSize / 2;
 111
 112     for (int i = 1; i < n; ++i) {
 113         Complex c1(realP1[i], imagP1[i]);
 114         Complex c2(realP2[i], imagP2[i]);
 115
 116         double mag1 = abs(c1);
 117         double mag2 = abs(c2);
 118
 119         // Interpolate magnitudes in decibels
 120         double mag1db = 20.0 * log10(mag1);
 121         double mag2db = 20.0 * log10(mag2);
 122
 123         double s1 = s1base;
 124         double s2 = s2base;
 125
 126         double magdbdiff = mag1db - mag2db;
 127
 128         // Empirical tweak to retain higher-frequency zeroes
 129         double threshold =  (i > 16) ? 5.0 : 2.0;
 130
 131         if (magdbdiff < -threshold && mag1db < 0.0) {
 132             s1 = pow(s1, 0.75);
 133             s2 = 1.0 - s1;
 134         } else if (magdbdiff > threshold && mag2db < 0.0) {
 135             s2 = pow(s2, 0.75);
 136             s1 = 1.0 - s2;
 137         }
 138
 139         // Average magnitude by decibels instead of linearly
 140         double magdb = s1 * mag1db + s2 * mag2db;
 141         double mag = pow(10.0, 0.05 * magdb);
 142
 143         // Now, deal with phase
 144         double phase1 = arg(c1);
 145         double phase2 = arg(c2);
 146
 147         double deltaPhase1 = phase1 - lastPhase1;
 148         double deltaPhase2 = phase2 - lastPhase2;
 149         lastPhase1 = phase1;
 150         lastPhase2 = phase2;
 151
 152         // Unwrap phase deltas
 153         if (deltaPhase1 > piDouble)
 154             deltaPhase1 -= 2.0 * piDouble;
 155         if (deltaPhase1 < -piDouble)
 156             deltaPhase1 += 2.0 * piDouble;
 157         if (deltaPhase2 > piDouble)
 158             deltaPhase2 -= 2.0 * piDouble;
 159         if (deltaPhase2 < -piDouble)
 160             deltaPhase2 += 2.0 * piDouble;
 161
 162         // Blend group-delays
 163         double deltaPhaseBlend;
 164
 165         if (deltaPhase1 - deltaPhase2 > piDouble)
 166             deltaPhaseBlend = s1 * deltaPhase1 + s2 * (2.0 * piDouble + deltaPhase2);
 167         else if (deltaPhase2 - deltaPhase1 > piDouble)
 168             deltaPhaseBlend = s1 * (2.0 * piDouble + deltaPhase1) + s2 * deltaPhase2;
 169         else
 170             deltaPhaseBlend = s1 * deltaPhase1 + s2 * deltaPhase2;
 171
 172         phaseAccum += deltaPhaseBlend;
 173
 174         // Unwrap
 175         if (phaseAccum > piDouble)
 176             phaseAccum -= 2.0 * piDouble;
 177         if (phaseAccum < -piDouble)
 178             phaseAccum += 2.0 * piDouble;
 179
 180         Complex c = complexFromMagnitudePhase(mag, phaseAccum);
 181
 182         realP[i] = static_cast<float>(c.real());
 183         imagP[i] = static_cast<float>(c.imag());
 184     }
 185 }
 186
 187 double FFTFrame::extractAverageGroupDelay()
 188 {
 189     float* realP = realData();
 190     float* imagP = imagData();
 191
 192     double aveSum = 0.0;
 193     double weightSum = 0.0;
 194     double lastPhase = 0.0;
 195
 196     int halfSize = fftSize() / 2;
 197
 198     const double kSamplePhaseDelay = (2.0 * piDouble) / double(fftSize());
 199
 200     // Calculate weighted average group delay
 201     for (int i = 0; i < halfSize; i++) {
 202         Complex c(realP[i], imagP[i]);
 203         double mag = abs(c);
 204         double phase = arg(c);
 205
 206         double deltaPhase = phase - lastPhase;
 207         lastPhase = phase;
 208
 209         // Unwrap
 210         if (deltaPhase < -piDouble)
 211             deltaPhase += 2.0 * piDouble;
 212         if (deltaPhase > piDouble)
 213             deltaPhase -= 2.0 * piDouble;
 214
 215         aveSum += mag * deltaPhase;
 216         weightSum += mag;
 217     }
 218
 219     // Note how we invert the phase delta wrt frequency since this is how group delay is defined
 220     double ave = aveSum / weightSum;
 221     double aveSampleDelay = -ave / kSamplePhaseDelay;
 222
 223     // Leave 20 sample headroom (for leading edge of impulse)
 224     if (aveSampleDelay > 20.0)
 225         aveSampleDelay -= 20.0;
 226
 227     // Remove average group delay (minus 20 samples for headroom)
 228     addConstantGroupDelay(-aveSampleDelay);
 229
 230     // Remove DC offset
 231     realP[0] = 0.0f;
 232
 233     return aveSampleDelay;
 234 }
 235
 236 void FFTFrame::addConstantGroupDelay(double sampleFrameDelay)
 237 {
 238     int halfSize = fftSize() / 2;
 239
 240     float* realP = realData();
 241     float* imagP = imagData();
 242
 243     const double kSamplePhaseDelay = (2.0 * piDouble) / double(fftSize());
 244
 245     double phaseAdj = -sampleFrameDelay * kSamplePhaseDelay;
 246
 247     // Add constant group delay
 248     for (int i = 1; i < halfSize; i++) {
 249         Complex c(realP[i], imagP[i]);
 250         double mag = abs(c);
 251         double phase = arg(c);
 252
 253         phase += i * phaseAdj;
 254
 255         Complex c2 = complexFromMagnitudePhase(mag, phase);
 256
 257         realP[i] = static_cast<float>(c2.real());
 258         imagP[i] = static_cast<float>(c2.imag());
 259     }
 260 }
 261
 262 #ifndef NDEBUG
 263 void FFTFrame::print()
 264 {
 265     FFTFrame& frame = *this;
 266     float* realP = frame.realData();
 267     float* imagP = frame.imagData();
 268     WTF_LOG(WebAudio, "**** \n");
 269     WTF_LOG(WebAudio, "DC = %f : nyquist = %f\n", realP[0], imagP[0]);
 270
 271     int n = m_FFTSize / 2;
 272
 273     for (int i = 1; i < n; i++) {
 274         double mag = sqrt(realP[i] * realP[i] + imagP[i] * imagP[i]);
 275         double phase = atan2(realP[i], imagP[i]);
 276
 277         WTF_LOG(WebAudio, "[%d] (%f %f)\n", i, mag, phase);
 278     }
 279     WTF_LOG(WebAudio, "****\n");
 280 }
 281 #endif // NDEBUG
 282
 283 } // namespace WebCore
 284
 285 #endif // ENABLE(WEB_AUDIO)