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