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