49f010b3cc27f7e88a0979316867e12e6a545780
[WebKit.git] / 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 <algorithm>
36 #include <float.h>
37 #include <math.h>
38 #include <stdio.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.resize(kBufferSize + 2);
53     m_outputBuffer.resize(kBufferSize + 2);
54 #endif
55
56     // Initialize as pass-thru (straight-wire, no filter effect)
57     m_a0 = 1.0;
58     m_a1 = 0.0;
59     m_a2 = 0.0;
60     m_b1 = 0.0;
61     m_b2 = 0.0;
62
63     m_g = 1.0;
64
65     reset(); // clear filter memory
66 }
67
68 void Biquad::process(const float* sourceP, float* destP, size_t framesToProcess)
69 {
70 #if OS(DARWIN)
71     // Use vecLib if available
72     processFast(sourceP, destP, framesToProcess);
73 #else
74     int n = framesToProcess;
75
76     // Create local copies of member variables
77     double x1 = m_x1;
78     double x2 = m_x2;
79     double y1 = m_y1;
80     double y2 = m_y2;
81
82     double a0 = m_a0;
83     double a1 = m_a1;
84     double a2 = m_a2;
85     double b1 = m_b1;
86     double b2 = m_b2;
87
88     while (n--) {
89         // FIXME: this can be optimized by pipelining the multiply adds...
90         float x = *sourceP++;
91         float y = a0*x + a1*x1 + a2*x2 - b1*y1 - b2*y2;
92
93         y *= m_g;
94
95         *destP++ = y;
96
97         // Update state variables
98         x2 = x1;
99         x1 = x;
100         y2 = y1;
101         y1 = y;
102     }
103
104     // Local variables back to member
105     m_x1 = x1;
106     m_x2 = x2;
107     m_y1 = y1;
108     m_y2 = y2;
109
110     m_a0 = a0;
111     m_a1 = a1;
112     m_a2 = a2;
113     m_b1 = b1;
114     m_b2 = b2;
115 #endif
116 }
117
118 #if OS(DARWIN)
119
120 // Here we have optimized version using Accelerate.framework
121
122 void Biquad::processFast(const float* sourceP, float* destP, size_t framesToProcess)
123 {
124     // Filter coefficients
125     double B[5];
126     B[0] = m_a0;
127     B[1] = m_a1;
128     B[2] = m_a2;
129     B[3] = m_b1;
130     B[4] = m_b2;
131
132     double* inputP = m_inputBuffer.data();
133     double* outputP = m_outputBuffer.data();
134
135     double* input2P = inputP + 2;
136     double* output2P = outputP + 2;
137
138     // Break up processing into smaller slices (kBufferSize) if necessary.
139
140     int n = framesToProcess;
141
142     while (n > 0) {
143         int framesThisTime = n < kBufferSize ? n : kBufferSize;
144
145         // Copy input to input buffer
146         for (int i = 0; i < framesThisTime; ++i)
147             input2P[i] = *sourceP++;
148
149         processSliceFast(inputP, outputP, B, framesThisTime);
150
151         // Copy output buffer to output (converts float -> double).
152         for (int i = 0; i < framesThisTime; ++i)
153             *destP++ = static_cast<float>(output2P[i]);
154
155         n -= framesThisTime;
156     }
157 }
158
159 void Biquad::processSliceFast(double* sourceP, double* destP, double* coefficientsP, size_t framesToProcess)
160 {
161     // Use double-precision for filter stability
162     vDSP_deq22D(sourceP, 1, coefficientsP, destP, 1, framesToProcess);
163
164     // Save history.  Note that sourceP and destP reference m_inputBuffer and m_outputBuffer respectively.
165     // These buffers are allocated (in the constructor) with space for two extra samples so it's OK to access
166     // array values two beyond framesToProcess.
167     sourceP[0] = sourceP[framesToProcess - 2 + 2];
168     sourceP[1] = sourceP[framesToProcess - 1 + 2];
169     destP[0] = destP[framesToProcess - 2 + 2];
170     destP[1] = destP[framesToProcess - 1 + 2];
171 }
172
173 #endif // OS(DARWIN)
174
175
176 void Biquad::reset()
177 {
178     m_x1 = m_x2 = m_y1 = m_y2 = 0.0;
179
180 #if OS(DARWIN)
181     // Two extra samples for filter history
182     double* inputP = m_inputBuffer.data();
183     inputP[0] = 0.0;
184     inputP[1] = 0.0;
185
186     double* outputP = m_outputBuffer.data();
187     outputP[0] = 0.0;
188     outputP[1] = 0.0;
189 #endif
190 }
191
192 void Biquad::setLowpassParams(double cutoff, double resonance)
193 {
194     resonance = std::max(0.0, resonance); // can't go negative
195
196     double g = pow(10.0, 0.05 * resonance);
197     double d = sqrt((4.0 - sqrt(16.0 - 16.0 / (g * g))) / 2.0);
198
199     // Compute biquad coefficients for lopass filter
200     double theta = M_PI * cutoff;
201     double sn = 0.5 * d * sin(theta);
202     double beta = 0.5 * (1.0 - sn) / (1.0 + sn);
203     double gamma = (0.5 + beta) * cos(theta);
204     double alpha = 0.25 * (0.5 + beta - gamma);
205
206     m_a0 = 2.0 * alpha;
207     m_a1 = 2.0 * 2.0*alpha;
208     m_a2 = 2.0 * alpha;
209     m_b1 = 2.0 * -gamma;
210     m_b2 = 2.0 * beta;
211 }
212
213 void Biquad::setHighpassParams(double cutoff, double resonance)
214 {
215     resonance = std::max(0.0, resonance); // can't go negative
216
217     double g = pow(10.0, 0.05 * resonance);
218     double d = sqrt((4.0 - sqrt(16.0 - 16.0 / (g * g))) / 2.0);
219
220     // Compute biquad coefficients for highpass filter
221     double theta = M_PI * cutoff;
222     double sn = 0.5 * d * sin(theta);
223     double beta = 0.5 * (1.0 - sn) / (1.0 + sn);
224     double gamma = (0.5 + beta) * cos(theta);
225     double alpha = 0.25 * (0.5 + beta + gamma);
226
227     m_a0 = 2.0 * alpha;
228     m_a1 = 2.0 * -2.0*alpha;
229     m_a2 = 2.0 * alpha;
230     m_b1 = 2.0 * -gamma;
231     m_b2 = 2.0 * beta;
232 }
233
234 void Biquad::setLowShelfParams(double cutoff, double dbGain)
235 {
236     double theta = M_PI * cutoff;
237
238     double A = pow(10.0, dbGain / 40.0);
239     double S = 1.0; // filter slope (1.0 is max value)
240     double alpha = 0.5 * sin(theta) * sqrt((A + 1.0 / A) * (1.0 / S - 1.0) + 2.0);
241
242     double k = cos(theta);
243     double k2 = 2.0 * sqrt(A) * alpha;
244
245     double b0 = A * ((A + 1.0) - (A - 1.0) * k + k2);
246     double b1 = 2.0 * A * ((A - 1.0) - (A + 1.0) * k);
247     double b2 = A * ((A + 1.0) - (A - 1.0) * k - k2);
248     double a0 = (A + 1.0) + (A - 1.0) * k + k2;
249     double a1 = -2.0 * ((A - 1.0) + (A + 1.0) * k);
250     double a2 = (A + 1.0) + (A - 1.0) * k - k2;
251
252     double a0Inverse = 1.0 / a0;
253     
254     m_a0 = b0 * a0Inverse;
255     m_a1 = b1 * a0Inverse;
256     m_a2 = b2 * a0Inverse;
257     m_b1 = a1 * a0Inverse;
258     m_b2 = a2 * a0Inverse;
259 }
260
261 void Biquad::setZeroPolePairs(const Complex &zero, const Complex &pole)
262 {
263     m_a0 = 1.0;
264     m_a1 = -2.0 * zero.real();
265
266     double zeroMag = abs(zero);
267     m_a2 = zeroMag * zeroMag;
268
269     m_b1 = -2.0 * pole.real();
270
271     double poleMag = abs(pole);
272     m_b2 = poleMag * poleMag;
273 }
274
275 void Biquad::setAllpassPole(const Complex &pole)
276 {
277     Complex zero = Complex(1.0, 0.0) / pole;
278     setZeroPolePairs(zero, pole);
279 }
280
281 } // namespace WebCore
282
283 #endif // ENABLE(WEB_AUDIO)