1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/content/media/webaudio/FFTBlock.cpp Wed Dec 31 06:09:35 2014 +0100
1.3 @@ -0,0 +1,227 @@
1.4 +/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
1.5 +/* vim:set ts=4 sw=4 sts=4 et cindent: */
1.6 +/*
1.7 + * Copyright (C) 2010 Google Inc. All rights reserved.
1.8 + *
1.9 + * Redistribution and use in source and binary forms, with or without
1.10 + * modification, are permitted provided that the following conditions
1.11 + * are met:
1.12 + *
1.13 + * 1. Redistributions of source code must retain the above copyright
1.14 + * notice, this list of conditions and the following disclaimer.
1.15 + * 2. Redistributions in binary form must reproduce the above copyright
1.16 + * notice, this list of conditions and the following disclaimer in the
1.17 + * documentation and/or other materials provided with the distribution.
1.18 + * 3. Neither the name of Apple Computer, Inc. ("Apple") nor the names of
1.19 + * its contributors may be used to endorse or promote products derived
1.20 + * from this software without specific prior written permission.
1.21 + *
1.22 + * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY
1.23 + * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
1.24 + * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
1.25 + * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY
1.26 + * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
1.27 + * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
1.28 + * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
1.29 + * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
1.30 + * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
1.31 + * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1.32 + */
1.33 +
1.34 +#include "FFTBlock.h"
1.35 +
1.36 +#include <complex>
1.37 +
1.38 +namespace mozilla {
1.39 +
1.40 +typedef std::complex<double> Complex;
1.41 +
1.42 +FFTBlock* FFTBlock::CreateInterpolatedBlock(const FFTBlock& block0, const FFTBlock& block1, double interp)
1.43 +{
1.44 + FFTBlock* newBlock = new FFTBlock(block0.FFTSize());
1.45 +
1.46 + newBlock->InterpolateFrequencyComponents(block0, block1, interp);
1.47 +
1.48 + // In the time-domain, the 2nd half of the response must be zero, to avoid circular convolution aliasing...
1.49 + int fftSize = newBlock->FFTSize();
1.50 + nsTArray<float> buffer;
1.51 + buffer.SetLength(fftSize);
1.52 + newBlock->GetInverseWithoutScaling(buffer.Elements());
1.53 + AudioBufferInPlaceScale(buffer.Elements(), 1.0f / fftSize, fftSize / 2);
1.54 + PodZero(buffer.Elements() + fftSize / 2, fftSize / 2);
1.55 +
1.56 + // Put back into frequency domain.
1.57 + newBlock->PerformFFT(buffer.Elements());
1.58 +
1.59 + return newBlock;
1.60 +}
1.61 +
1.62 +void FFTBlock::InterpolateFrequencyComponents(const FFTBlock& block0, const FFTBlock& block1, double interp)
1.63 +{
1.64 + // FIXME : with some work, this method could be optimized
1.65 +
1.66 + kiss_fft_cpx* dft = mOutputBuffer.Elements();
1.67 +
1.68 + const kiss_fft_cpx* dft1 = block0.mOutputBuffer.Elements();
1.69 + const kiss_fft_cpx* dft2 = block1.mOutputBuffer.Elements();
1.70 +
1.71 + MOZ_ASSERT(mFFTSize == block0.FFTSize());
1.72 + MOZ_ASSERT(mFFTSize == block1.FFTSize());
1.73 + double s1base = (1.0 - interp);
1.74 + double s2base = interp;
1.75 +
1.76 + double phaseAccum = 0.0;
1.77 + double lastPhase1 = 0.0;
1.78 + double lastPhase2 = 0.0;
1.79 +
1.80 + int n = mFFTSize / 2;
1.81 +
1.82 + dft[0].r = static_cast<float>(s1base * dft1[0].r + s2base * dft2[0].r);
1.83 + dft[n].r = static_cast<float>(s1base * dft1[n].r + s2base * dft2[n].r);
1.84 +
1.85 + for (int i = 1; i < n; ++i) {
1.86 + Complex c1(dft1[i].r, dft1[i].i);
1.87 + Complex c2(dft2[i].r, dft2[i].i);
1.88 +
1.89 + double mag1 = abs(c1);
1.90 + double mag2 = abs(c2);
1.91 +
1.92 + // Interpolate magnitudes in decibels
1.93 + double mag1db = 20.0 * log10(mag1);
1.94 + double mag2db = 20.0 * log10(mag2);
1.95 +
1.96 + double s1 = s1base;
1.97 + double s2 = s2base;
1.98 +
1.99 + double magdbdiff = mag1db - mag2db;
1.100 +
1.101 + // Empirical tweak to retain higher-frequency zeroes
1.102 + double threshold = (i > 16) ? 5.0 : 2.0;
1.103 +
1.104 + if (magdbdiff < -threshold && mag1db < 0.0) {
1.105 + s1 = pow(s1, 0.75);
1.106 + s2 = 1.0 - s1;
1.107 + } else if (magdbdiff > threshold && mag2db < 0.0) {
1.108 + s2 = pow(s2, 0.75);
1.109 + s1 = 1.0 - s2;
1.110 + }
1.111 +
1.112 + // Average magnitude by decibels instead of linearly
1.113 + double magdb = s1 * mag1db + s2 * mag2db;
1.114 + double mag = pow(10.0, 0.05 * magdb);
1.115 +
1.116 + // Now, deal with phase
1.117 + double phase1 = arg(c1);
1.118 + double phase2 = arg(c2);
1.119 +
1.120 + double deltaPhase1 = phase1 - lastPhase1;
1.121 + double deltaPhase2 = phase2 - lastPhase2;
1.122 + lastPhase1 = phase1;
1.123 + lastPhase2 = phase2;
1.124 +
1.125 + // Unwrap phase deltas
1.126 + if (deltaPhase1 > M_PI)
1.127 + deltaPhase1 -= 2.0 * M_PI;
1.128 + if (deltaPhase1 < -M_PI)
1.129 + deltaPhase1 += 2.0 * M_PI;
1.130 + if (deltaPhase2 > M_PI)
1.131 + deltaPhase2 -= 2.0 * M_PI;
1.132 + if (deltaPhase2 < -M_PI)
1.133 + deltaPhase2 += 2.0 * M_PI;
1.134 +
1.135 + // Blend group-delays
1.136 + double deltaPhaseBlend;
1.137 +
1.138 + if (deltaPhase1 - deltaPhase2 > M_PI)
1.139 + deltaPhaseBlend = s1 * deltaPhase1 + s2 * (2.0 * M_PI + deltaPhase2);
1.140 + else if (deltaPhase2 - deltaPhase1 > M_PI)
1.141 + deltaPhaseBlend = s1 * (2.0 * M_PI + deltaPhase1) + s2 * deltaPhase2;
1.142 + else
1.143 + deltaPhaseBlend = s1 * deltaPhase1 + s2 * deltaPhase2;
1.144 +
1.145 + phaseAccum += deltaPhaseBlend;
1.146 +
1.147 + // Unwrap
1.148 + if (phaseAccum > M_PI)
1.149 + phaseAccum -= 2.0 * M_PI;
1.150 + if (phaseAccum < -M_PI)
1.151 + phaseAccum += 2.0 * M_PI;
1.152 +
1.153 + dft[i].r = static_cast<float>(mag * cos(phaseAccum));
1.154 + dft[i].i = static_cast<float>(mag * sin(phaseAccum));
1.155 + }
1.156 +}
1.157 +
1.158 +double FFTBlock::ExtractAverageGroupDelay()
1.159 +{
1.160 + kiss_fft_cpx* dft = mOutputBuffer.Elements();
1.161 +
1.162 + double aveSum = 0.0;
1.163 + double weightSum = 0.0;
1.164 + double lastPhase = 0.0;
1.165 +
1.166 + int halfSize = FFTSize() / 2;
1.167 +
1.168 + const double kSamplePhaseDelay = (2.0 * M_PI) / double(FFTSize());
1.169 +
1.170 + // Remove DC offset
1.171 + dft[0].r = 0.0f;
1.172 +
1.173 + // Calculate weighted average group delay
1.174 + for (int i = 1; i < halfSize; i++) {
1.175 + Complex c(dft[i].r, dft[i].i);
1.176 + double mag = abs(c);
1.177 + double phase = arg(c);
1.178 +
1.179 + double deltaPhase = phase - lastPhase;
1.180 + lastPhase = phase;
1.181 +
1.182 + // Unwrap
1.183 + if (deltaPhase < -M_PI)
1.184 + deltaPhase += 2.0 * M_PI;
1.185 + if (deltaPhase > M_PI)
1.186 + deltaPhase -= 2.0 * M_PI;
1.187 +
1.188 + aveSum += mag * deltaPhase;
1.189 + weightSum += mag;
1.190 + }
1.191 +
1.192 + // Note how we invert the phase delta wrt frequency since this is how group delay is defined
1.193 + double ave = aveSum / weightSum;
1.194 + double aveSampleDelay = -ave / kSamplePhaseDelay;
1.195 +
1.196 + // Leave 20 sample headroom (for leading edge of impulse)
1.197 + aveSampleDelay -= 20.0;
1.198 + if (aveSampleDelay <= 0.0)
1.199 + return 0.0;
1.200 +
1.201 + // Remove average group delay (minus 20 samples for headroom)
1.202 + AddConstantGroupDelay(-aveSampleDelay);
1.203 +
1.204 + return aveSampleDelay;
1.205 +}
1.206 +
1.207 +void FFTBlock::AddConstantGroupDelay(double sampleFrameDelay)
1.208 +{
1.209 + int halfSize = FFTSize() / 2;
1.210 +
1.211 + kiss_fft_cpx* dft = mOutputBuffer.Elements();
1.212 +
1.213 + const double kSamplePhaseDelay = (2.0 * M_PI) / double(FFTSize());
1.214 +
1.215 + double phaseAdj = -sampleFrameDelay * kSamplePhaseDelay;
1.216 +
1.217 + // Add constant group delay
1.218 + for (int i = 1; i < halfSize; i++) {
1.219 + Complex c(dft[i].r, dft[i].i);
1.220 + double mag = abs(c);
1.221 + double phase = arg(c);
1.222 +
1.223 + phase += i * phaseAdj;
1.224 +
1.225 + dft[i].r = static_cast<float>(mag * cos(phase));
1.226 + dft[i].i = static_cast<float>(mag * sin(phase));
1.227 + }
1.228 +}
1.229 +
1.230 +} // namespace mozilla
