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