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