Enable IPP for Biquad filter
[WebKit-https.git] / Source / WebCore / platform / audio / Biquad.cpp
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 "Biquad.h"
34
35 #include "DenormalDisabler.h"
36 #include <algorithm>
37 #include <stdio.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)