The Tor Browser / file comparison
comparison: content/media/webaudio/FFTBlock.cpp
content/media/webaudio/FFTBlock.cpp
- branch
- TOR_BUG_9701
- changeset 13
- 44a2da4a2ab2
equal
deleted
inserted
replaced
| -1:000000000000 | 0:5b16a8c18c51 |
|---|---|
| 1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */ | |
| 2 /* vim:set ts=4 sw=4 sts=4 et cindent: */ | |
| 3 /* | |
| 4 * Copyright (C) 2010 Google Inc. All rights reserved. | |
| 5 * | |
| 6 * Redistribution and use in source and binary forms, with or without | |
| 7 * modification, are permitted provided that the following conditions | |
| 8 * are met: | |
| 9 * | |
| 10 * 1. Redistributions of source code must retain the above copyright | |
| 11 * notice, this list of conditions and the following disclaimer. | |
| 12 * 2. Redistributions in binary form must reproduce the above copyright | |
| 13 * notice, this list of conditions and the following disclaimer in the | |
| 14 * documentation and/or other materials provided with the distribution. | |
| 15 * 3. Neither the name of Apple Computer, Inc. ("Apple") nor the names of | |
| 16 * its contributors may be used to endorse or promote products derived | |
| 17 * from this software without specific prior written permission. | |
| 18 * | |
| 19 * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY | |
| 20 * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED | |
| 21 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE | |
| 22 * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY | |
| 23 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES | |
| 24 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; | |
| 25 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND | |
| 26 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT | |
| 27 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF | |
| 28 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | |
| 29 */ | |
| 30 | |
| 31 #include "FFTBlock.h" | |
| 32 | |
| 33 #include <complex> | |
| 34 | |
| 35 namespace mozilla { | |
| 36 | |
| 37 typedef std::complex<double> Complex; | |
| 38 | |
| 39 FFTBlock* FFTBlock::CreateInterpolatedBlock(const FFTBlock& block0, const FFTBlock& block1, double interp) | |
| 40 { | |
| 41 FFTBlock* newBlock = new FFTBlock(block0.FFTSize()); | |
| 42 | |
| 43 newBlock->InterpolateFrequencyComponents(block0, block1, interp); | |
| 44 | |
| 45 // In the time-domain, the 2nd half of the response must be zero, to avoid circular convolution aliasing... | |
| 46 int fftSize = newBlock->FFTSize(); | |
| 47 nsTArray<float> buffer; | |
| 48 buffer.SetLength(fftSize); | |
| 49 newBlock->GetInverseWithoutScaling(buffer.Elements()); | |
| 50 AudioBufferInPlaceScale(buffer.Elements(), 1.0f / fftSize, fftSize / 2); | |
| 51 PodZero(buffer.Elements() + fftSize / 2, fftSize / 2); | |
| 52 | |
| 53 // Put back into frequency domain. | |
| 54 newBlock->PerformFFT(buffer.Elements()); | |
| 55 | |
| 56 return newBlock; | |
| 57 } | |
| 58 | |
| 59 void FFTBlock::InterpolateFrequencyComponents(const FFTBlock& block0, const FFTBlock& block1, double interp) | |
| 60 { | |
| 61 // FIXME : with some work, this method could be optimized | |
| 62 | |
| 63 kiss_fft_cpx* dft = mOutputBuffer.Elements(); | |
| 64 | |
| 65 const kiss_fft_cpx* dft1 = block0.mOutputBuffer.Elements(); | |
| 66 const kiss_fft_cpx* dft2 = block1.mOutputBuffer.Elements(); | |
| 67 | |
| 68 MOZ_ASSERT(mFFTSize == block0.FFTSize()); | |
| 69 MOZ_ASSERT(mFFTSize == block1.FFTSize()); | |
| 70 double s1base = (1.0 - interp); | |
| 71 double s2base = interp; | |
| 72 | |
| 73 double phaseAccum = 0.0; | |
| 74 double lastPhase1 = 0.0; | |
| 75 double lastPhase2 = 0.0; | |
| 76 | |
| 77 int n = mFFTSize / 2; | |
| 78 | |
| 79 dft[0].r = static_cast<float>(s1base * dft1[0].r + s2base * dft2[0].r); | |
| 80 dft[n].r = static_cast<float>(s1base * dft1[n].r + s2base * dft2[n].r); | |
| 81 | |
| 82 for (int i = 1; i < n; ++i) { | |
| 83 Complex c1(dft1[i].r, dft1[i].i); | |
| 84 Complex c2(dft2[i].r, dft2[i].i); | |
| 85 | |
| 86 double mag1 = abs(c1); | |
| 87 double mag2 = abs(c2); | |
| 88 | |
| 89 // Interpolate magnitudes in decibels | |
| 90 double mag1db = 20.0 * log10(mag1); | |
| 91 double mag2db = 20.0 * log10(mag2); | |
| 92 | |
| 93 double s1 = s1base; | |
| 94 double s2 = s2base; | |
| 95 | |
| 96 double magdbdiff = mag1db - mag2db; | |
| 97 | |
| 98 // Empirical tweak to retain higher-frequency zeroes | |
| 99 double threshold = (i > 16) ? 5.0 : 2.0; | |
| 100 | |
| 101 if (magdbdiff < -threshold && mag1db < 0.0) { | |
| 102 s1 = pow(s1, 0.75); | |
| 103 s2 = 1.0 - s1; | |
| 104 } else if (magdbdiff > threshold && mag2db < 0.0) { | |
| 105 s2 = pow(s2, 0.75); | |
| 106 s1 = 1.0 - s2; | |
| 107 } | |
| 108 | |
| 109 // Average magnitude by decibels instead of linearly | |
| 110 double magdb = s1 * mag1db + s2 * mag2db; | |
| 111 double mag = pow(10.0, 0.05 * magdb); | |
| 112 | |
| 113 // Now, deal with phase | |
| 114 double phase1 = arg(c1); | |
| 115 double phase2 = arg(c2); | |
| 116 | |
| 117 double deltaPhase1 = phase1 - lastPhase1; | |
| 118 double deltaPhase2 = phase2 - lastPhase2; | |
| 119 lastPhase1 = phase1; | |
| 120 lastPhase2 = phase2; | |
| 121 | |
| 122 // Unwrap phase deltas | |
| 123 if (deltaPhase1 > M_PI) | |
| 124 deltaPhase1 -= 2.0 * M_PI; | |
| 125 if (deltaPhase1 < -M_PI) | |
| 126 deltaPhase1 += 2.0 * M_PI; | |
| 127 if (deltaPhase2 > M_PI) | |
| 128 deltaPhase2 -= 2.0 * M_PI; | |
| 129 if (deltaPhase2 < -M_PI) | |
| 130 deltaPhase2 += 2.0 * M_PI; | |
| 131 | |
| 132 // Blend group-delays | |
| 133 double deltaPhaseBlend; | |
| 134 | |
| 135 if (deltaPhase1 - deltaPhase2 > M_PI) | |
| 136 deltaPhaseBlend = s1 * deltaPhase1 + s2 * (2.0 * M_PI + deltaPhase2); | |
| 137 else if (deltaPhase2 - deltaPhase1 > M_PI) | |
| 138 deltaPhaseBlend = s1 * (2.0 * M_PI + deltaPhase1) + s2 * deltaPhase2; | |
| 139 else | |
| 140 deltaPhaseBlend = s1 * deltaPhase1 + s2 * deltaPhase2; | |
| 141 | |
| 142 phaseAccum += deltaPhaseBlend; | |
| 143 | |
| 144 // Unwrap | |
| 145 if (phaseAccum > M_PI) | |
| 146 phaseAccum -= 2.0 * M_PI; | |
| 147 if (phaseAccum < -M_PI) | |
| 148 phaseAccum += 2.0 * M_PI; | |
| 149 | |
| 150 dft[i].r = static_cast<float>(mag * cos(phaseAccum)); | |
| 151 dft[i].i = static_cast<float>(mag * sin(phaseAccum)); | |
| 152 } | |
| 153 } | |
| 154 | |
| 155 double FFTBlock::ExtractAverageGroupDelay() | |
| 156 { | |
| 157 kiss_fft_cpx* dft = mOutputBuffer.Elements(); | |
| 158 | |
| 159 double aveSum = 0.0; | |
| 160 double weightSum = 0.0; | |
| 161 double lastPhase = 0.0; | |
| 162 | |
| 163 int halfSize = FFTSize() / 2; | |
| 164 | |
| 165 const double kSamplePhaseDelay = (2.0 * M_PI) / double(FFTSize()); | |
| 166 | |
| 167 // Remove DC offset | |
| 168 dft[0].r = 0.0f; | |
| 169 | |
| 170 // Calculate weighted average group delay | |
| 171 for (int i = 1; i < halfSize; i++) { | |
| 172 Complex c(dft[i].r, dft[i].i); | |
| 173 double mag = abs(c); | |
| 174 double phase = arg(c); | |
| 175 | |
| 176 double deltaPhase = phase - lastPhase; | |
| 177 lastPhase = phase; | |
| 178 | |
| 179 // Unwrap | |
| 180 if (deltaPhase < -M_PI) | |
| 181 deltaPhase += 2.0 * M_PI; | |
| 182 if (deltaPhase > M_PI) | |
| 183 deltaPhase -= 2.0 * M_PI; | |
| 184 | |
| 185 aveSum += mag * deltaPhase; | |
| 186 weightSum += mag; | |
| 187 } | |
| 188 | |
| 189 // Note how we invert the phase delta wrt frequency since this is how group delay is defined | |
| 190 double ave = aveSum / weightSum; | |
| 191 double aveSampleDelay = -ave / kSamplePhaseDelay; | |
| 192 | |
| 193 // Leave 20 sample headroom (for leading edge of impulse) | |
| 194 aveSampleDelay -= 20.0; | |
| 195 if (aveSampleDelay <= 0.0) | |
| 196 return 0.0; | |
| 197 | |
| 198 // Remove average group delay (minus 20 samples for headroom) | |
| 199 AddConstantGroupDelay(-aveSampleDelay); | |
| 200 | |
| 201 return aveSampleDelay; | |
| 202 } | |
| 203 | |
| 204 void FFTBlock::AddConstantGroupDelay(double sampleFrameDelay) | |
| 205 { | |
| 206 int halfSize = FFTSize() / 2; | |
| 207 | |
| 208 kiss_fft_cpx* dft = mOutputBuffer.Elements(); | |
| 209 | |
| 210 const double kSamplePhaseDelay = (2.0 * M_PI) / double(FFTSize()); | |
| 211 | |
| 212 double phaseAdj = -sampleFrameDelay * kSamplePhaseDelay; | |
| 213 | |
| 214 // Add constant group delay | |
| 215 for (int i = 1; i < halfSize; i++) { | |
| 216 Complex c(dft[i].r, dft[i].i); | |
| 217 double mag = abs(c); | |
| 218 double phase = arg(c); | |
| 219 | |
| 220 phase += i * phaseAdj; | |
| 221 | |
| 222 dft[i].r = static_cast<float>(mag * cos(phase)); | |
| 223 dft[i].i = static_cast<float>(mag * sin(phase)); | |
| 224 } | |
| 225 } | |
| 226 | |
| 227 } // namespace mozilla |
