• Home
  • History
  • Annotate
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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