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
43namespace WebCore {
44
45const int kBufferSize = 1024;
46
47Biquad::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
67void 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
121void 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
158void 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
175void 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
191void 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
212void 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
233void 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
260void 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
274void 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