FFTBlock.cpp (7498B)
1 /* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- */ 2 /* vim:set ts=4 sw=2 sts=2 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 #include "FFVPXRuntimeLinker.h" 36 37 namespace mozilla { 38 39 FFmpegFFTFuncs FFTBlock::sFFTFuncs = {}; 40 41 using Complex = std::complex<double>; 42 43 static double fdlibm_cabs(const Complex& z) { 44 return fdlibm_hypot(real(z), imag(z)); 45 } 46 47 static double fdlibm_carg(const Complex& z) { 48 return fdlibm_atan2(imag(z), real(z)); 49 } 50 51 FFTBlock* FFTBlock::CreateInterpolatedBlock(const FFTBlock& block0, 52 const FFTBlock& block1, 53 double interp) { 54 uint32_t fftSize = block0.FFTSize(); 55 FFTBlock* newBlock = 56 new FFTBlock(fftSize, 1.0f / AssertedCast<float>(fftSize)); 57 58 newBlock->InterpolateFrequencyComponents(block0, block1, interp); 59 60 // In the time-domain, the 2nd half of the response must be zero, to avoid 61 // circular convolution aliasing... 62 AlignedTArray<float> buffer(fftSize); 63 newBlock->GetInverse(buffer.Elements()); 64 PodZero(buffer.Elements() + fftSize / 2, fftSize / 2); 65 66 // Put back into frequency domain. 67 newBlock->PerformFFT(buffer.Elements()); 68 69 return newBlock; 70 } 71 72 void FFTBlock::InterpolateFrequencyComponents(const FFTBlock& block0, 73 const FFTBlock& block1, 74 double interp) { 75 // FIXME : with some work, this method could be optimized 76 77 ComplexU* dft = mOutputBuffer.Elements(); 78 79 const ComplexU* dft1 = block0.mOutputBuffer.Elements(); 80 const ComplexU* dft2 = block1.mOutputBuffer.Elements(); 81 82 MOZ_ASSERT(mFFTSize == block0.FFTSize()); 83 MOZ_ASSERT(mFFTSize == block1.FFTSize()); 84 double s1base = (1.0 - interp); 85 double s2base = interp; 86 87 double phaseAccum = 0.0; 88 double lastPhase1 = 0.0; 89 double lastPhase2 = 0.0; 90 91 int n = mFFTSize / 2; 92 93 dft[0].r = static_cast<float>(s1base * dft1[0].r + s2base * dft2[0].r); 94 dft[n].r = static_cast<float>(s1base * dft1[n].r + s2base * dft2[n].r); 95 96 for (int i = 1; i < n; ++i) { 97 Complex c1(dft1[i].r, dft1[i].i); 98 Complex c2(dft2[i].r, dft2[i].i); 99 100 double mag1 = fdlibm_cabs(c1); 101 double mag2 = fdlibm_cabs(c2); 102 103 // Interpolate magnitudes in decibels 104 double mag1db = 20.0 * fdlibm_log10(mag1); 105 double mag2db = 20.0 * fdlibm_log10(mag2); 106 107 double s1 = s1base; 108 double s2 = s2base; 109 110 double magdbdiff = mag1db - mag2db; 111 112 // Empirical tweak to retain higher-frequency zeroes 113 double threshold = (i > 16) ? 5.0 : 2.0; 114 115 if (magdbdiff < -threshold && mag1db < 0.0) { 116 s1 = fdlibm_pow(s1, 0.75); 117 s2 = 1.0 - s1; 118 } else if (magdbdiff > threshold && mag2db < 0.0) { 119 s2 = fdlibm_pow(s2, 0.75); 120 s1 = 1.0 - s2; 121 } 122 123 // Average magnitude by decibels instead of linearly 124 double magdb = s1 * mag1db + s2 * mag2db; 125 double mag = fdlibm_pow(10.0, 0.05 * magdb); 126 127 // Now, deal with phase 128 double phase1 = fdlibm_carg(c1); 129 double phase2 = fdlibm_carg(c2); 130 131 double deltaPhase1 = phase1 - lastPhase1; 132 double deltaPhase2 = phase2 - lastPhase2; 133 lastPhase1 = phase1; 134 lastPhase2 = phase2; 135 136 // Unwrap phase deltas 137 if (deltaPhase1 > M_PI) deltaPhase1 -= 2.0 * M_PI; 138 if (deltaPhase1 < -M_PI) deltaPhase1 += 2.0 * M_PI; 139 if (deltaPhase2 > M_PI) deltaPhase2 -= 2.0 * M_PI; 140 if (deltaPhase2 < -M_PI) deltaPhase2 += 2.0 * M_PI; 141 142 // Blend group-delays 143 double deltaPhaseBlend; 144 145 if (deltaPhase1 - deltaPhase2 > M_PI) 146 deltaPhaseBlend = s1 * deltaPhase1 + s2 * (2.0 * M_PI + deltaPhase2); 147 else if (deltaPhase2 - deltaPhase1 > M_PI) 148 deltaPhaseBlend = s1 * (2.0 * M_PI + deltaPhase1) + s2 * deltaPhase2; 149 else 150 deltaPhaseBlend = s1 * deltaPhase1 + s2 * deltaPhase2; 151 152 phaseAccum += deltaPhaseBlend; 153 154 // Unwrap 155 if (phaseAccum > M_PI) phaseAccum -= 2.0 * M_PI; 156 if (phaseAccum < -M_PI) phaseAccum += 2.0 * M_PI; 157 158 dft[i].r = static_cast<float>(mag * fdlibm_cos(phaseAccum)); 159 dft[i].i = static_cast<float>(mag * fdlibm_sin(phaseAccum)); 160 } 161 } 162 163 double FFTBlock::ExtractAverageGroupDelay() { 164 ComplexU* dft = mOutputBuffer.Elements(); 165 166 double aveSum = 0.0; 167 double weightSum = 0.0; 168 double lastPhase = 0.0; 169 170 int halfSize = FFTSize() / 2; 171 172 const double kSamplePhaseDelay = (2.0 * M_PI) / double(FFTSize()); 173 174 // Remove DC offset 175 dft[0].r = 0.0f; 176 177 // Calculate weighted average group delay 178 for (int i = 1; i < halfSize; i++) { 179 Complex c(dft[i].r, dft[i].i); 180 double mag = fdlibm_cabs(c); 181 double phase = fdlibm_carg(c); 182 183 double deltaPhase = phase - lastPhase; 184 lastPhase = phase; 185 186 // Unwrap 187 if (deltaPhase < -M_PI) deltaPhase += 2.0 * M_PI; 188 if (deltaPhase > M_PI) deltaPhase -= 2.0 * M_PI; 189 190 aveSum += mag * deltaPhase; 191 weightSum += mag; 192 } 193 194 // Note how we invert the phase delta wrt frequency since this is how group 195 // delay is defined 196 double ave = aveSum / weightSum; 197 double aveSampleDelay = -ave / kSamplePhaseDelay; 198 199 // Leave 20 sample headroom (for leading edge of impulse) 200 aveSampleDelay -= 20.0; 201 if (aveSampleDelay <= 0.0) return 0.0; 202 203 // Remove average group delay (minus 20 samples for headroom) 204 AddConstantGroupDelay(-aveSampleDelay); 205 206 return aveSampleDelay; 207 } 208 209 void FFTBlock::AddConstantGroupDelay(double sampleFrameDelay) { 210 int halfSize = FFTSize() / 2; 211 212 ComplexU* dft = mOutputBuffer.Elements(); 213 214 const double kSamplePhaseDelay = (2.0 * M_PI) / double(FFTSize()); 215 216 double phaseAdj = -sampleFrameDelay * kSamplePhaseDelay; 217 218 // Add constant group delay 219 for (int i = 1; i < halfSize; i++) { 220 Complex c(dft[i].r, dft[i].i); 221 double mag = fdlibm_cabs(c); 222 double phase = fdlibm_carg(c); 223 224 phase += i * phaseAdj; 225 226 dft[i].r = static_cast<float>(mag * fdlibm_cos(phase)); 227 dft[i].i = static_cast<float>(mag * fdlibm_sin(phase)); 228 } 229 } 230 231 } // namespace mozilla