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