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/Biquad.h" 34 35 #include <stdio.h> 36 #include <algorithm> 37 #include "core/platform/audio/DenormalDisabler.h" 38 #include "wtf/MathExtras.h" 39 40 #if OS(DARWIN) 41 #include <Accelerate/Accelerate.h> 42 #endif 43 44 namespace WebCore { 45 46 const int kBufferSize = 1024; 47 48 Biquad::Biquad() 49 { 50 #if OS(DARWIN) 51 // Allocate two samples more for filter history 52 m_inputBuffer.allocate(kBufferSize + 2); 53 m_outputBuffer.allocate(kBufferSize + 2); 54 #endif 55 56 #if USE(WEBAUDIO_IPP) 57 int bufferSize; 58 ippsIIRGetStateSize64f_BiQuad_32f(1, &bufferSize); 59 m_ippInternalBuffer = ippsMalloc_8u(bufferSize); 60 #endif // USE(WEBAUDIO_IPP) 61 62 // Initialize as pass-thru (straight-wire, no filter effect) 63 setNormalizedCoefficients(1, 0, 0, 1, 0, 0); 64 65 reset(); // clear filter memory 66 } 67 68 Biquad::~Biquad() 69 { 70 #if USE(WEBAUDIO_IPP) 71 ippsFree(m_ippInternalBuffer); 72 #endif // USE(WEBAUDIO_IPP) 73 } 74 75 void Biquad::process(const float* sourceP, float* destP, size_t framesToProcess) 76 { 77 #if OS(DARWIN) 78 // Use vecLib if available 79 processFast(sourceP, destP, framesToProcess); 80 81 #elif USE(WEBAUDIO_IPP) 82 ippsIIR64f_32f(sourceP, destP, static_cast<int>(framesToProcess), m_biquadState); 83 #else // USE(WEBAUDIO_IPP) 84 85 int n = framesToProcess; 86 87 // Create local copies of member variables 88 double x1 = m_x1; 89 double x2 = m_x2; 90 double y1 = m_y1; 91 double y2 = m_y2; 92 93 double b0 = m_b0; 94 double b1 = m_b1; 95 double b2 = m_b2; 96 double a1 = m_a1; 97 double a2 = m_a2; 98 99 while (n--) { 100 // FIXME: this can be optimized by pipelining the multiply adds... 101 float x = *sourceP++; 102 float y = b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2; 103 104 *destP++ = y; 105 106 // Update state variables 107 x2 = x1; 108 x1 = x; 109 y2 = y1; 110 y1 = y; 111 } 112 113 // Local variables back to member. Flush denormals here so we 114 // don't slow down the inner loop above. 115 m_x1 = DenormalDisabler::flushDenormalFloatToZero(x1); 116 m_x2 = DenormalDisabler::flushDenormalFloatToZero(x2); 117 m_y1 = DenormalDisabler::flushDenormalFloatToZero(y1); 118 m_y2 = DenormalDisabler::flushDenormalFloatToZero(y2); 119 120 m_b0 = b0; 121 m_b1 = b1; 122 m_b2 = b2; 123 m_a1 = a1; 124 m_a2 = a2; 125 #endif 126 } 127 128 #if OS(DARWIN) 129 130 // Here we have optimized version using Accelerate.framework 131 132 void Biquad::processFast(const float* sourceP, float* destP, size_t framesToProcess) 133 { 134 double filterCoefficients[5]; 135 filterCoefficients[0] = m_b0; 136 filterCoefficients[1] = m_b1; 137 filterCoefficients[2] = m_b2; 138 filterCoefficients[3] = m_a1; 139 filterCoefficients[4] = m_a2; 140 141 double* inputP = m_inputBuffer.data(); 142 double* outputP = m_outputBuffer.data(); 143 144 double* input2P = inputP + 2; 145 double* output2P = outputP + 2; 146 147 // Break up processing into smaller slices (kBufferSize) if necessary. 148 149 int n = framesToProcess; 150 151 while (n > 0) { 152 int framesThisTime = n < kBufferSize ? n : kBufferSize; 153 154 // Copy input to input buffer 155 for (int i = 0; i < framesThisTime; ++i) 156 input2P[i] = *sourceP++; 157 158 processSliceFast(inputP, outputP, filterCoefficients, framesThisTime); 159 160 // Copy output buffer to output (converts float -> double). 161 for (int i = 0; i < framesThisTime; ++i) 162 *destP++ = static_cast<float>(output2P[i]); 163 164 n -= framesThisTime; 165 } 166 } 167 168 void Biquad::processSliceFast(double* sourceP, double* destP, double* coefficientsP, size_t framesToProcess) 169 { 170 // Use double-precision for filter stability 171 vDSP_deq22D(sourceP, 1, coefficientsP, destP, 1, framesToProcess); 172 173 // Save history. Note that sourceP and destP reference m_inputBuffer and m_outputBuffer respectively. 174 // These buffers are allocated (in the constructor) with space for two extra samples so it's OK to access 175 // array values two beyond framesToProcess. 176 sourceP[0] = sourceP[framesToProcess - 2 + 2]; 177 sourceP[1] = sourceP[framesToProcess - 1 + 2]; 178 destP[0] = destP[framesToProcess - 2 + 2]; 179 destP[1] = destP[framesToProcess - 1 + 2]; 180 } 181 182 #endif // OS(DARWIN) 183 184 185 void Biquad::reset() 186 { 187 #if OS(DARWIN) 188 // Two extra samples for filter history 189 double* inputP = m_inputBuffer.data(); 190 inputP[0] = 0; 191 inputP[1] = 0; 192 193 double* outputP = m_outputBuffer.data(); 194 outputP[0] = 0; 195 outputP[1] = 0; 196 197 #elif USE(WEBAUDIO_IPP) 198 int bufferSize; 199 ippsIIRGetStateSize64f_BiQuad_32f(1, &bufferSize); 200 ippsZero_8u(m_ippInternalBuffer, bufferSize); 201 202 #else 203 m_x1 = m_x2 = m_y1 = m_y2 = 0; 204 #endif 205 } 206 207 void Biquad::setLowpassParams(double cutoff, double resonance) 208 { 209 // Limit cutoff to 0 to 1. 210 cutoff = std::max(0.0, std::min(cutoff, 1.0)); 211 212 if (cutoff == 1) { 213 // When cutoff is 1, the z-transform is 1. 214 setNormalizedCoefficients(1, 0, 0, 215 1, 0, 0); 216 } else if (cutoff > 0) { 217 // Compute biquad coefficients for lowpass filter 218 resonance = std::max(0.0, resonance); // can't go negative 219 double g = pow(10.0, 0.05 * resonance); 220 double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2); 221 222 double theta = piDouble * cutoff; 223 double sn = 0.5 * d * sin(theta); 224 double beta = 0.5 * (1 - sn) / (1 + sn); 225 double gamma = (0.5 + beta) * cos(theta); 226 double alpha = 0.25 * (0.5 + beta - gamma); 227 228 double b0 = 2 * alpha; 229 double b1 = 2 * 2 * alpha; 230 double b2 = 2 * alpha; 231 double a1 = 2 * -gamma; 232 double a2 = 2 * beta; 233 234 setNormalizedCoefficients(b0, b1, b2, 1, a1, a2); 235 } else { 236 // When cutoff is zero, nothing gets through the filter, so set 237 // coefficients up correctly. 238 setNormalizedCoefficients(0, 0, 0, 239 1, 0, 0); 240 } 241 } 242 243 void Biquad::setHighpassParams(double cutoff, double resonance) 244 { 245 // Limit cutoff to 0 to 1. 246 cutoff = std::max(0.0, std::min(cutoff, 1.0)); 247 248 if (cutoff == 1) { 249 // The z-transform is 0. 250 setNormalizedCoefficients(0, 0, 0, 251 1, 0, 0); 252 } else if (cutoff > 0) { 253 // Compute biquad coefficients for highpass filter 254 resonance = std::max(0.0, resonance); // can't go negative 255 double g = pow(10.0, 0.05 * resonance); 256 double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2); 257 258 double theta = piDouble * cutoff; 259 double sn = 0.5 * d * sin(theta); 260 double beta = 0.5 * (1 - sn) / (1 + sn); 261 double gamma = (0.5 + beta) * cos(theta); 262 double alpha = 0.25 * (0.5 + beta + gamma); 263 264 double b0 = 2 * alpha; 265 double b1 = 2 * -2 * alpha; 266 double b2 = 2 * alpha; 267 double a1 = 2 * -gamma; 268 double a2 = 2 * beta; 269 270 setNormalizedCoefficients(b0, b1, b2, 1, a1, a2); 271 } else { 272 // When cutoff is zero, we need to be careful because the above 273 // gives a quadratic divided by the same quadratic, with poles 274 // and zeros on the unit circle in the same place. When cutoff 275 // is zero, the z-transform is 1. 276 setNormalizedCoefficients(1, 0, 0, 277 1, 0, 0); 278 } 279 } 280 281 void Biquad::setNormalizedCoefficients(double b0, double b1, double b2, double a0, double a1, double a2) 282 { 283 double a0Inverse = 1 / a0; 284 285 m_b0 = b0 * a0Inverse; 286 m_b1 = b1 * a0Inverse; 287 m_b2 = b2 * a0Inverse; 288 m_a1 = a1 * a0Inverse; 289 m_a2 = a2 * a0Inverse; 290 291 #if USE(WEBAUDIO_IPP) 292 Ipp64f taps[6]; 293 taps[0] = m_b0; 294 taps[1] = m_b1; 295 taps[2] = m_b2; 296 taps[3] = 1; 297 taps[4] = m_a1; 298 taps[5] = m_a2; 299 m_biquadState = 0; 300 301 ippsIIRInit64f_BiQuad_32f(&m_biquadState, taps, 1, 0, m_ippInternalBuffer); 302 #endif // USE(WEBAUDIO_IPP) 303 } 304 305 void Biquad::setLowShelfParams(double frequency, double dbGain) 306 { 307 // Clip frequencies to between 0 and 1, inclusive. 308 frequency = std::max(0.0, std::min(frequency, 1.0)); 309 310 double A = pow(10.0, dbGain / 40); 311 312 if (frequency == 1) { 313 // The z-transform is a constant gain. 314 setNormalizedCoefficients(A * A, 0, 0, 315 1, 0, 0); 316 } else if (frequency > 0) { 317 double w0 = piDouble * frequency; 318 double S = 1; // filter slope (1 is max value) 319 double alpha = 0.5 * sin(w0) * sqrt((A + 1 / A) * (1 / S - 1) + 2); 320 double k = cos(w0); 321 double k2 = 2 * sqrt(A) * alpha; 322 double aPlusOne = A + 1; 323 double aMinusOne = A - 1; 324 325 double b0 = A * (aPlusOne - aMinusOne * k + k2); 326 double b1 = 2 * A * (aMinusOne - aPlusOne * k); 327 double b2 = A * (aPlusOne - aMinusOne * k - k2); 328 double a0 = aPlusOne + aMinusOne * k + k2; 329 double a1 = -2 * (aMinusOne + aPlusOne * k); 330 double a2 = aPlusOne + aMinusOne * k - k2; 331 332 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); 333 } else { 334 // When frequency is 0, the z-transform is 1. 335 setNormalizedCoefficients(1, 0, 0, 336 1, 0, 0); 337 } 338 } 339 340 void Biquad::setHighShelfParams(double frequency, double dbGain) 341 { 342 // Clip frequencies to between 0 and 1, inclusive. 343 frequency = std::max(0.0, std::min(frequency, 1.0)); 344 345 double A = pow(10.0, dbGain / 40); 346 347 if (frequency == 1) { 348 // The z-transform is 1. 349 setNormalizedCoefficients(1, 0, 0, 350 1, 0, 0); 351 } else if (frequency > 0) { 352 double w0 = piDouble * frequency; 353 double S = 1; // filter slope (1 is max value) 354 double alpha = 0.5 * sin(w0) * sqrt((A + 1 / A) * (1 / S - 1) + 2); 355 double k = cos(w0); 356 double k2 = 2 * sqrt(A) * alpha; 357 double aPlusOne = A + 1; 358 double aMinusOne = A - 1; 359 360 double b0 = A * (aPlusOne + aMinusOne * k + k2); 361 double b1 = -2 * A * (aMinusOne + aPlusOne * k); 362 double b2 = A * (aPlusOne + aMinusOne * k - k2); 363 double a0 = aPlusOne - aMinusOne * k + k2; 364 double a1 = 2 * (aMinusOne - aPlusOne * k); 365 double a2 = aPlusOne - aMinusOne * k - k2; 366 367 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); 368 } else { 369 // When frequency = 0, the filter is just a gain, A^2. 370 setNormalizedCoefficients(A * A, 0, 0, 371 1, 0, 0); 372 } 373 } 374 375 void Biquad::setPeakingParams(double frequency, double Q, double dbGain) 376 { 377 // Clip frequencies to between 0 and 1, inclusive. 378 frequency = std::max(0.0, std::min(frequency, 1.0)); 379 380 // Don't let Q go negative, which causes an unstable filter. 381 Q = std::max(0.0, Q); 382 383 double A = pow(10.0, dbGain / 40); 384 385 if (frequency > 0 && frequency < 1) { 386 if (Q > 0) { 387 double w0 = piDouble * frequency; 388 double alpha = sin(w0) / (2 * Q); 389 double k = cos(w0); 390 391 double b0 = 1 + alpha * A; 392 double b1 = -2 * k; 393 double b2 = 1 - alpha * A; 394 double a0 = 1 + alpha / A; 395 double a1 = -2 * k; 396 double a2 = 1 - alpha / A; 397 398 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); 399 } else { 400 // When Q = 0, the above formulas have problems. If we look at 401 // the z-transform, we can see that the limit as Q->0 is A^2, so 402 // set the filter that way. 403 setNormalizedCoefficients(A * A, 0, 0, 404 1, 0, 0); 405 } 406 } else { 407 // When frequency is 0 or 1, the z-transform is 1. 408 setNormalizedCoefficients(1, 0, 0, 409 1, 0, 0); 410 } 411 } 412 413 void Biquad::setAllpassParams(double frequency, double Q) 414 { 415 // Clip frequencies to between 0 and 1, inclusive. 416 frequency = std::max(0.0, std::min(frequency, 1.0)); 417 418 // Don't let Q go negative, which causes an unstable filter. 419 Q = std::max(0.0, Q); 420 421 if (frequency > 0 && frequency < 1) { 422 if (Q > 0) { 423 double w0 = piDouble * frequency; 424 double alpha = sin(w0) / (2 * Q); 425 double k = cos(w0); 426 427 double b0 = 1 - alpha; 428 double b1 = -2 * k; 429 double b2 = 1 + alpha; 430 double a0 = 1 + alpha; 431 double a1 = -2 * k; 432 double a2 = 1 - alpha; 433 434 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); 435 } else { 436 // When Q = 0, the above formulas have problems. If we look at 437 // the z-transform, we can see that the limit as Q->0 is -1, so 438 // set the filter that way. 439 setNormalizedCoefficients(-1, 0, 0, 440 1, 0, 0); 441 } 442 } else { 443 // When frequency is 0 or 1, the z-transform is 1. 444 setNormalizedCoefficients(1, 0, 0, 445 1, 0, 0); 446 } 447 } 448 449 void Biquad::setNotchParams(double frequency, double Q) 450 { 451 // Clip frequencies to between 0 and 1, inclusive. 452 frequency = std::max(0.0, std::min(frequency, 1.0)); 453 454 // Don't let Q go negative, which causes an unstable filter. 455 Q = std::max(0.0, Q); 456 457 if (frequency > 0 && frequency < 1) { 458 if (Q > 0) { 459 double w0 = piDouble * frequency; 460 double alpha = sin(w0) / (2 * Q); 461 double k = cos(w0); 462 463 double b0 = 1; 464 double b1 = -2 * k; 465 double b2 = 1; 466 double a0 = 1 + alpha; 467 double a1 = -2 * k; 468 double a2 = 1 - alpha; 469 470 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); 471 } else { 472 // When Q = 0, the above formulas have problems. If we look at 473 // the z-transform, we can see that the limit as Q->0 is 0, so 474 // set the filter that way. 475 setNormalizedCoefficients(0, 0, 0, 476 1, 0, 0); 477 } 478 } else { 479 // When frequency is 0 or 1, the z-transform is 1. 480 setNormalizedCoefficients(1, 0, 0, 481 1, 0, 0); 482 } 483 } 484 485 void Biquad::setBandpassParams(double frequency, double Q) 486 { 487 // No negative frequencies allowed. 488 frequency = std::max(0.0, frequency); 489 490 // Don't let Q go negative, which causes an unstable filter. 491 Q = std::max(0.0, Q); 492 493 if (frequency > 0 && frequency < 1) { 494 double w0 = piDouble * frequency; 495 if (Q > 0) { 496 double alpha = sin(w0) / (2 * Q); 497 double k = cos(w0); 498 499 double b0 = alpha; 500 double b1 = 0; 501 double b2 = -alpha; 502 double a0 = 1 + alpha; 503 double a1 = -2 * k; 504 double a2 = 1 - alpha; 505 506 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); 507 } else { 508 // When Q = 0, the above formulas have problems. If we look at 509 // the z-transform, we can see that the limit as Q->0 is 1, so 510 // set the filter that way. 511 setNormalizedCoefficients(1, 0, 0, 512 1, 0, 0); 513 } 514 } else { 515 // When the cutoff is zero, the z-transform approaches 0, if Q 516 // > 0. When both Q and cutoff are zero, the z-transform is 517 // pretty much undefined. What should we do in this case? 518 // For now, just make the filter 0. When the cutoff is 1, the 519 // z-transform also approaches 0. 520 setNormalizedCoefficients(0, 0, 0, 521 1, 0, 0); 522 } 523 } 524 525 void Biquad::setZeroPolePairs(const Complex &zero, const Complex &pole) 526 { 527 double b0 = 1; 528 double b1 = -2 * zero.real(); 529 530 double zeroMag = abs(zero); 531 double b2 = zeroMag * zeroMag; 532 533 double a1 = -2 * pole.real(); 534 535 double poleMag = abs(pole); 536 double a2 = poleMag * poleMag; 537 setNormalizedCoefficients(b0, b1, b2, 1, a1, a2); 538 } 539 540 void Biquad::setAllpassPole(const Complex &pole) 541 { 542 Complex zero = Complex(1, 0) / pole; 543 setZeroPolePairs(zero, pole); 544 } 545 546 void Biquad::getFrequencyResponse(int nFrequencies, 547 const float* frequency, 548 float* magResponse, 549 float* phaseResponse) 550 { 551 // Evaluate the Z-transform of the filter at given normalized 552 // frequency from 0 to 1. (1 corresponds to the Nyquist 553 // frequency.) 554 // 555 // The z-transform of the filter is 556 // 557 // H(z) = (b0 + b1*z^(-1) + b2*z^(-2))/(1 + a1*z^(-1) + a2*z^(-2)) 558 // 559 // Evaluate as 560 // 561 // b0 + (b1 + b2*z1)*z1 562 // -------------------- 563 // 1 + (a1 + a2*z1)*z1 564 // 565 // with z1 = 1/z and z = exp(j*pi*frequency). Hence z1 = exp(-j*pi*frequency) 566 567 // Make local copies of the coefficients as a micro-optimization. 568 double b0 = m_b0; 569 double b1 = m_b1; 570 double b2 = m_b2; 571 double a1 = m_a1; 572 double a2 = m_a2; 573 574 for (int k = 0; k < nFrequencies; ++k) { 575 double omega = -piDouble * frequency[k]; 576 Complex z = Complex(cos(omega), sin(omega)); 577 Complex numerator = b0 + (b1 + b2 * z) * z; 578 Complex denominator = Complex(1, 0) + (a1 + a2 * z) * z; 579 Complex response = numerator / denominator; 580 magResponse[k] = static_cast<float>(abs(response)); 581 phaseResponse[k] = static_cast<float>(atan2(imag(response), real(response))); 582 } 583 } 584 585 } // namespace WebCore 586 587 #endif // ENABLE(WEB_AUDIO) 588