content/media/webaudio/blink/Biquad.cpp@97036ab72558
content/media/webaudio/blink/Biquad.cpp
Tue, 06 Jan 2015 21:39:09 +0100
- author
- Michael Schloh von Bennewitz <michael@schloh.com>
- date
- Tue, 06 Jan 2015 21:39:09 +0100
- branch
- TOR_BUG_9701
- changeset 8
- 97036ab72558
- permissions
- -rw-r--r--
Conditionally force memory storage according to privacy.thirdparty.isolate;
This solves Tor bug #9701, complying with disk avoidance documented in
https://www.torproject.org/projects/torbrowser/design/#disk-avoidance.
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 */
29 // For M_PI from cmath
30 #ifdef _MSC_VER
31 # define _USE_MATH_DEFINES
32 #endif
34 #include "Biquad.h"
36 #include <cmath>
37 #include <float.h>
38 #include <algorithm>
40 namespace WebCore {
42 Biquad::Biquad()
43 {
44 // Initialize as pass-thru (straight-wire, no filter effect)
45 setNormalizedCoefficients(1, 0, 0, 1, 0, 0);
47 reset(); // clear filter memory
48 }
50 Biquad::~Biquad()
51 {
52 }
54 void Biquad::process(const float* sourceP, float* destP, size_t framesToProcess)
55 {
56 // Create local copies of member variables
57 double x1 = m_x1;
58 double x2 = m_x2;
59 double y1 = m_y1;
60 double y2 = m_y2;
62 double b0 = m_b0;
63 double b1 = m_b1;
64 double b2 = m_b2;
65 double a1 = m_a1;
66 double a2 = m_a2;
68 for (size_t i = 0; i < framesToProcess; ++i) {
69 // FIXME: this can be optimized by pipelining the multiply adds...
70 double x = sourceP[i];
71 double y = b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2;
73 destP[i] = y;
75 // Update state variables
76 x2 = x1;
77 x1 = x;
78 y2 = y1;
79 y1 = y;
80 }
82 // Avoid introducing a stream of subnormals when input is silent and the
83 // tail approaches zero.
84 if (x1 == 0.0 && x2 == 0.0 && (y1 != 0.0 || y2 != 0.0) &&
85 fabs(y1) < FLT_MIN && fabs(y2) < FLT_MIN) {
86 // Flush future values to zero (until there is new input).
87 y1 = y2 = 0.0;
88 // Flush calculated values.
89 for (int i = framesToProcess; i-- && fabsf(destP[i]) < FLT_MIN; ) {
90 destP[i] = 0.0f;
91 }
92 }
93 // Local variables back to member.
94 m_x1 = x1;
95 m_x2 = x2;
96 m_y1 = y1;
97 m_y2 = y2;
98 }
100 void Biquad::reset()
101 {
102 m_x1 = m_x2 = m_y1 = m_y2 = 0;
103 }
105 void Biquad::setLowpassParams(double cutoff, double resonance)
106 {
107 // Limit cutoff to 0 to 1.
108 cutoff = std::max(0.0, std::min(cutoff, 1.0));
110 if (cutoff == 1) {
111 // When cutoff is 1, the z-transform is 1.
112 setNormalizedCoefficients(1, 0, 0,
113 1, 0, 0);
114 } else if (cutoff > 0) {
115 // Compute biquad coefficients for lowpass filter
116 resonance = std::max(0.0, resonance); // can't go negative
117 double g = pow(10.0, 0.05 * resonance);
118 double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2);
120 double theta = M_PI * cutoff;
121 double sn = 0.5 * d * sin(theta);
122 double beta = 0.5 * (1 - sn) / (1 + sn);
123 double gamma = (0.5 + beta) * cos(theta);
124 double alpha = 0.25 * (0.5 + beta - gamma);
126 double b0 = 2 * alpha;
127 double b1 = 2 * 2 * alpha;
128 double b2 = 2 * alpha;
129 double a1 = 2 * -gamma;
130 double a2 = 2 * beta;
132 setNormalizedCoefficients(b0, b1, b2, 1, a1, a2);
133 } else {
134 // When cutoff is zero, nothing gets through the filter, so set
135 // coefficients up correctly.
136 setNormalizedCoefficients(0, 0, 0,
137 1, 0, 0);
138 }
139 }
141 void Biquad::setHighpassParams(double cutoff, double resonance)
142 {
143 // Limit cutoff to 0 to 1.
144 cutoff = std::max(0.0, std::min(cutoff, 1.0));
146 if (cutoff == 1) {
147 // The z-transform is 0.
148 setNormalizedCoefficients(0, 0, 0,
149 1, 0, 0);
150 } else if (cutoff > 0) {
151 // Compute biquad coefficients for highpass filter
152 resonance = std::max(0.0, resonance); // can't go negative
153 double g = pow(10.0, 0.05 * resonance);
154 double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2);
156 double theta = M_PI * cutoff;
157 double sn = 0.5 * d * sin(theta);
158 double beta = 0.5 * (1 - sn) / (1 + sn);
159 double gamma = (0.5 + beta) * cos(theta);
160 double alpha = 0.25 * (0.5 + beta + gamma);
162 double b0 = 2 * alpha;
163 double b1 = 2 * -2 * alpha;
164 double b2 = 2 * alpha;
165 double a1 = 2 * -gamma;
166 double a2 = 2 * beta;
168 setNormalizedCoefficients(b0, b1, b2, 1, a1, a2);
169 } else {
170 // When cutoff is zero, we need to be careful because the above
171 // gives a quadratic divided by the same quadratic, with poles
172 // and zeros on the unit circle in the same place. When cutoff
173 // is zero, the z-transform is 1.
174 setNormalizedCoefficients(1, 0, 0,
175 1, 0, 0);
176 }
177 }
179 void Biquad::setNormalizedCoefficients(double b0, double b1, double b2, double a0, double a1, double a2)
180 {
181 double a0Inverse = 1 / a0;
183 m_b0 = b0 * a0Inverse;
184 m_b1 = b1 * a0Inverse;
185 m_b2 = b2 * a0Inverse;
186 m_a1 = a1 * a0Inverse;
187 m_a2 = a2 * a0Inverse;
188 }
190 void Biquad::setLowShelfParams(double frequency, double dbGain)
191 {
192 // Clip frequencies to between 0 and 1, inclusive.
193 frequency = std::max(0.0, std::min(frequency, 1.0));
195 double A = pow(10.0, dbGain / 40);
197 if (frequency == 1) {
198 // The z-transform is a constant gain.
199 setNormalizedCoefficients(A * A, 0, 0,
200 1, 0, 0);
201 } else if (frequency > 0) {
202 double w0 = M_PI * frequency;
203 double S = 1; // filter slope (1 is max value)
204 double alpha = 0.5 * sin(w0) * sqrt((A + 1 / A) * (1 / S - 1) + 2);
205 double k = cos(w0);
206 double k2 = 2 * sqrt(A) * alpha;
207 double aPlusOne = A + 1;
208 double aMinusOne = A - 1;
210 double b0 = A * (aPlusOne - aMinusOne * k + k2);
211 double b1 = 2 * A * (aMinusOne - aPlusOne * k);
212 double b2 = A * (aPlusOne - aMinusOne * k - k2);
213 double a0 = aPlusOne + aMinusOne * k + k2;
214 double a1 = -2 * (aMinusOne + aPlusOne * k);
215 double a2 = aPlusOne + aMinusOne * k - k2;
217 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2);
218 } else {
219 // When frequency is 0, the z-transform is 1.
220 setNormalizedCoefficients(1, 0, 0,
221 1, 0, 0);
222 }
223 }
225 void Biquad::setHighShelfParams(double frequency, double dbGain)
226 {
227 // Clip frequencies to between 0 and 1, inclusive.
228 frequency = std::max(0.0, std::min(frequency, 1.0));
230 double A = pow(10.0, dbGain / 40);
232 if (frequency == 1) {
233 // The z-transform is 1.
234 setNormalizedCoefficients(1, 0, 0,
235 1, 0, 0);
236 } else if (frequency > 0) {
237 double w0 = M_PI * frequency;
238 double S = 1; // filter slope (1 is max value)
239 double alpha = 0.5 * sin(w0) * sqrt((A + 1 / A) * (1 / S - 1) + 2);
240 double k = cos(w0);
241 double k2 = 2 * sqrt(A) * alpha;
242 double aPlusOne = A + 1;
243 double aMinusOne = A - 1;
245 double b0 = A * (aPlusOne + aMinusOne * k + k2);
246 double b1 = -2 * A * (aMinusOne + aPlusOne * k);
247 double b2 = A * (aPlusOne + aMinusOne * k - k2);
248 double a0 = aPlusOne - aMinusOne * k + k2;
249 double a1 = 2 * (aMinusOne - aPlusOne * k);
250 double a2 = aPlusOne - aMinusOne * k - k2;
252 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2);
253 } else {
254 // When frequency = 0, the filter is just a gain, A^2.
255 setNormalizedCoefficients(A * A, 0, 0,
256 1, 0, 0);
257 }
258 }
260 void Biquad::setPeakingParams(double frequency, double Q, double dbGain)
261 {
262 // Clip frequencies to between 0 and 1, inclusive.
263 frequency = std::max(0.0, std::min(frequency, 1.0));
265 // Don't let Q go negative, which causes an unstable filter.
266 Q = std::max(0.0, Q);
268 double A = pow(10.0, dbGain / 40);
270 if (frequency > 0 && frequency < 1) {
271 if (Q > 0) {
272 double w0 = M_PI * frequency;
273 double alpha = sin(w0) / (2 * Q);
274 double k = cos(w0);
276 double b0 = 1 + alpha * A;
277 double b1 = -2 * k;
278 double b2 = 1 - alpha * A;
279 double a0 = 1 + alpha / A;
280 double a1 = -2 * k;
281 double a2 = 1 - alpha / A;
283 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2);
284 } else {
285 // When Q = 0, the above formulas have problems. If we look at
286 // the z-transform, we can see that the limit as Q->0 is A^2, so
287 // set the filter that way.
288 setNormalizedCoefficients(A * A, 0, 0,
289 1, 0, 0);
290 }
291 } else {
292 // When frequency is 0 or 1, the z-transform is 1.
293 setNormalizedCoefficients(1, 0, 0,
294 1, 0, 0);
295 }
296 }
298 void Biquad::setAllpassParams(double frequency, double Q)
299 {
300 // Clip frequencies to between 0 and 1, inclusive.
301 frequency = std::max(0.0, std::min(frequency, 1.0));
303 // Don't let Q go negative, which causes an unstable filter.
304 Q = std::max(0.0, Q);
306 if (frequency > 0 && frequency < 1) {
307 if (Q > 0) {
308 double w0 = M_PI * frequency;
309 double alpha = sin(w0) / (2 * Q);
310 double k = cos(w0);
312 double b0 = 1 - alpha;
313 double b1 = -2 * k;
314 double b2 = 1 + alpha;
315 double a0 = 1 + alpha;
316 double a1 = -2 * k;
317 double a2 = 1 - alpha;
319 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2);
320 } else {
321 // When Q = 0, the above formulas have problems. If we look at
322 // the z-transform, we can see that the limit as Q->0 is -1, so
323 // set the filter that way.
324 setNormalizedCoefficients(-1, 0, 0,
325 1, 0, 0);
326 }
327 } else {
328 // When frequency is 0 or 1, the z-transform is 1.
329 setNormalizedCoefficients(1, 0, 0,
330 1, 0, 0);
331 }
332 }
334 void Biquad::setNotchParams(double frequency, double Q)
335 {
336 // Clip frequencies to between 0 and 1, inclusive.
337 frequency = std::max(0.0, std::min(frequency, 1.0));
339 // Don't let Q go negative, which causes an unstable filter.
340 Q = std::max(0.0, Q);
342 if (frequency > 0 && frequency < 1) {
343 if (Q > 0) {
344 double w0 = M_PI * frequency;
345 double alpha = sin(w0) / (2 * Q);
346 double k = cos(w0);
348 double b0 = 1;
349 double b1 = -2 * k;
350 double b2 = 1;
351 double a0 = 1 + alpha;
352 double a1 = -2 * k;
353 double a2 = 1 - alpha;
355 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2);
356 } else {
357 // When Q = 0, the above formulas have problems. If we look at
358 // the z-transform, we can see that the limit as Q->0 is 0, so
359 // set the filter that way.
360 setNormalizedCoefficients(0, 0, 0,
361 1, 0, 0);
362 }
363 } else {
364 // When frequency is 0 or 1, the z-transform is 1.
365 setNormalizedCoefficients(1, 0, 0,
366 1, 0, 0);
367 }
368 }
370 void Biquad::setBandpassParams(double frequency, double Q)
371 {
372 // No negative frequencies allowed.
373 frequency = std::max(0.0, frequency);
375 // Don't let Q go negative, which causes an unstable filter.
376 Q = std::max(0.0, Q);
378 if (frequency > 0 && frequency < 1) {
379 double w0 = M_PI * frequency;
380 if (Q > 0) {
381 double alpha = sin(w0) / (2 * Q);
382 double k = cos(w0);
384 double b0 = alpha;
385 double b1 = 0;
386 double b2 = -alpha;
387 double a0 = 1 + alpha;
388 double a1 = -2 * k;
389 double a2 = 1 - alpha;
391 setNormalizedCoefficients(b0, b1, b2, a0, a1, a2);
392 } else {
393 // When Q = 0, the above formulas have problems. If we look at
394 // the z-transform, we can see that the limit as Q->0 is 1, so
395 // set the filter that way.
396 setNormalizedCoefficients(1, 0, 0,
397 1, 0, 0);
398 }
399 } else {
400 // When the cutoff is zero, the z-transform approaches 0, if Q
401 // > 0. When both Q and cutoff are zero, the z-transform is
402 // pretty much undefined. What should we do in this case?
403 // For now, just make the filter 0. When the cutoff is 1, the
404 // z-transform also approaches 0.
405 setNormalizedCoefficients(0, 0, 0,
406 1, 0, 0);
407 }
408 }
410 void Biquad::setZeroPolePairs(const Complex &zero, const Complex &pole)
411 {
412 double b0 = 1;
413 double b1 = -2 * zero.real();
415 double zeroMag = abs(zero);
416 double b2 = zeroMag * zeroMag;
418 double a1 = -2 * pole.real();
420 double poleMag = abs(pole);
421 double a2 = poleMag * poleMag;
422 setNormalizedCoefficients(b0, b1, b2, 1, a1, a2);
423 }
425 void Biquad::setAllpassPole(const Complex &pole)
426 {
427 Complex zero = Complex(1, 0) / pole;
428 setZeroPolePairs(zero, pole);
429 }
431 void Biquad::getFrequencyResponse(int nFrequencies,
432 const float* frequency,
433 float* magResponse,
434 float* phaseResponse)
435 {
436 // Evaluate the Z-transform of the filter at given normalized
437 // frequency from 0 to 1. (1 corresponds to the Nyquist
438 // frequency.)
439 //
440 // The z-transform of the filter is
441 //
442 // H(z) = (b0 + b1*z^(-1) + b2*z^(-2))/(1 + a1*z^(-1) + a2*z^(-2))
443 //
444 // Evaluate as
445 //
446 // b0 + (b1 + b2*z1)*z1
447 // --------------------
448 // 1 + (a1 + a2*z1)*z1
449 //
450 // with z1 = 1/z and z = exp(j*pi*frequency). Hence z1 = exp(-j*pi*frequency)
452 // Make local copies of the coefficients as a micro-optimization.
453 double b0 = m_b0;
454 double b1 = m_b1;
455 double b2 = m_b2;
456 double a1 = m_a1;
457 double a2 = m_a2;
459 for (int k = 0; k < nFrequencies; ++k) {
460 double omega = -M_PI * frequency[k];
461 Complex z = Complex(cos(omega), sin(omega));
462 Complex numerator = b0 + (b1 + b2 * z) * z;
463 Complex denominator = Complex(1, 0) + (a1 + a2 * z) * z;
464 Complex response = numerator / denominator;
465 magResponse[k] = static_cast<float>(abs(response));
466 phaseResponse[k] = static_cast<float>(atan2(imag(response), real(response)));
467 }
468 }
470 } // namespace WebCore
