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