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