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
44 namespace blink {
45
46 #if OS(MACOSX)
47 const int kBufferSize = 1024;
48 #endif
49
Biquad()50 Biquad::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
~Biquad()70 Biquad::~Biquad()
71 {
72 #if USE(WEBAUDIO_IPP)
73 ippsFree(m_ippInternalBuffer);
74 #endif // USE(WEBAUDIO_IPP)
75 }
76
process(const float * sourceP,float * destP,size_t framesToProcess)77 void 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
processFast(const float * sourceP,float * destP,size_t framesToProcess)134 void 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
processSliceFast(double * sourceP,double * destP,double * coefficientsP,size_t framesToProcess)170 void 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
reset()187 void 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
setLowpassParams(double cutoff,double resonance)209 void 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
setHighpassParams(double cutoff,double resonance)245 void 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
setNormalizedCoefficients(double b0,double b1,double b2,double a0,double a1,double a2)283 void 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
setLowShelfParams(double frequency,double dbGain)307 void 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
setHighShelfParams(double frequency,double dbGain)342 void 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
setPeakingParams(double frequency,double Q,double dbGain)377 void 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
setAllpassParams(double frequency,double Q)415 void 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
setNotchParams(double frequency,double Q)451 void 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
setBandpassParams(double frequency,double Q)487 void 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
setZeroPolePairs(const Complex & zero,const Complex & pole)527 void 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
setAllpassPole(const Complex & pole)542 void Biquad::setAllpassPole(const Complex &pole)
543 {
544 Complex zero = Complex(1, 0) / pole;
545 setZeroPolePairs(zero, pole);
546 }
547
getFrequencyResponse(int nFrequencies,const float * frequency,float * magResponse,float * phaseResponse)548 void 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