Tue, 06 Jan 2015 21:39:09 +0100
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.
michael@0 | 1 | /* |
michael@0 | 2 | * Copyright (C) 2010 Google Inc. All rights reserved. |
michael@0 | 3 | * |
michael@0 | 4 | * Redistribution and use in source and binary forms, with or without |
michael@0 | 5 | * modification, are permitted provided that the following conditions |
michael@0 | 6 | * are met: |
michael@0 | 7 | * |
michael@0 | 8 | * 1. Redistributions of source code must retain the above copyright |
michael@0 | 9 | * notice, this list of conditions and the following disclaimer. |
michael@0 | 10 | * 2. Redistributions in binary form must reproduce the above copyright |
michael@0 | 11 | * notice, this list of conditions and the following disclaimer in the |
michael@0 | 12 | * documentation and/or other materials provided with the distribution. |
michael@0 | 13 | * 3. Neither the name of Apple Computer, Inc. ("Apple") nor the names of |
michael@0 | 14 | * its contributors may be used to endorse or promote products derived |
michael@0 | 15 | * from this software without specific prior written permission. |
michael@0 | 16 | * |
michael@0 | 17 | * THIS SOFTWARE IS PROVIDED BY APPLE AND ITS CONTRIBUTORS "AS IS" AND ANY |
michael@0 | 18 | * EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED |
michael@0 | 19 | * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE |
michael@0 | 20 | * DISCLAIMED. IN NO EVENT SHALL APPLE OR ITS CONTRIBUTORS BE LIABLE FOR ANY |
michael@0 | 21 | * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES |
michael@0 | 22 | * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; |
michael@0 | 23 | * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND |
michael@0 | 24 | * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
michael@0 | 25 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF |
michael@0 | 26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
michael@0 | 27 | */ |
michael@0 | 28 | |
michael@0 | 29 | // For M_PI from cmath |
michael@0 | 30 | #ifdef _MSC_VER |
michael@0 | 31 | # define _USE_MATH_DEFINES |
michael@0 | 32 | #endif |
michael@0 | 33 | |
michael@0 | 34 | #include "Biquad.h" |
michael@0 | 35 | |
michael@0 | 36 | #include <cmath> |
michael@0 | 37 | #include <float.h> |
michael@0 | 38 | #include <algorithm> |
michael@0 | 39 | |
michael@0 | 40 | namespace WebCore { |
michael@0 | 41 | |
michael@0 | 42 | Biquad::Biquad() |
michael@0 | 43 | { |
michael@0 | 44 | // Initialize as pass-thru (straight-wire, no filter effect) |
michael@0 | 45 | setNormalizedCoefficients(1, 0, 0, 1, 0, 0); |
michael@0 | 46 | |
michael@0 | 47 | reset(); // clear filter memory |
michael@0 | 48 | } |
michael@0 | 49 | |
michael@0 | 50 | Biquad::~Biquad() |
michael@0 | 51 | { |
michael@0 | 52 | } |
michael@0 | 53 | |
michael@0 | 54 | void Biquad::process(const float* sourceP, float* destP, size_t framesToProcess) |
michael@0 | 55 | { |
michael@0 | 56 | // Create local copies of member variables |
michael@0 | 57 | double x1 = m_x1; |
michael@0 | 58 | double x2 = m_x2; |
michael@0 | 59 | double y1 = m_y1; |
michael@0 | 60 | double y2 = m_y2; |
michael@0 | 61 | |
michael@0 | 62 | double b0 = m_b0; |
michael@0 | 63 | double b1 = m_b1; |
michael@0 | 64 | double b2 = m_b2; |
michael@0 | 65 | double a1 = m_a1; |
michael@0 | 66 | double a2 = m_a2; |
michael@0 | 67 | |
michael@0 | 68 | for (size_t i = 0; i < framesToProcess; ++i) { |
michael@0 | 69 | // FIXME: this can be optimized by pipelining the multiply adds... |
michael@0 | 70 | double x = sourceP[i]; |
michael@0 | 71 | double y = b0*x + b1*x1 + b2*x2 - a1*y1 - a2*y2; |
michael@0 | 72 | |
michael@0 | 73 | destP[i] = y; |
michael@0 | 74 | |
michael@0 | 75 | // Update state variables |
michael@0 | 76 | x2 = x1; |
michael@0 | 77 | x1 = x; |
michael@0 | 78 | y2 = y1; |
michael@0 | 79 | y1 = y; |
michael@0 | 80 | } |
michael@0 | 81 | |
michael@0 | 82 | // Avoid introducing a stream of subnormals when input is silent and the |
michael@0 | 83 | // tail approaches zero. |
michael@0 | 84 | if (x1 == 0.0 && x2 == 0.0 && (y1 != 0.0 || y2 != 0.0) && |
michael@0 | 85 | fabs(y1) < FLT_MIN && fabs(y2) < FLT_MIN) { |
michael@0 | 86 | // Flush future values to zero (until there is new input). |
michael@0 | 87 | y1 = y2 = 0.0; |
michael@0 | 88 | // Flush calculated values. |
michael@0 | 89 | for (int i = framesToProcess; i-- && fabsf(destP[i]) < FLT_MIN; ) { |
michael@0 | 90 | destP[i] = 0.0f; |
michael@0 | 91 | } |
michael@0 | 92 | } |
michael@0 | 93 | // Local variables back to member. |
michael@0 | 94 | m_x1 = x1; |
michael@0 | 95 | m_x2 = x2; |
michael@0 | 96 | m_y1 = y1; |
michael@0 | 97 | m_y2 = y2; |
michael@0 | 98 | } |
michael@0 | 99 | |
michael@0 | 100 | void Biquad::reset() |
michael@0 | 101 | { |
michael@0 | 102 | m_x1 = m_x2 = m_y1 = m_y2 = 0; |
michael@0 | 103 | } |
michael@0 | 104 | |
michael@0 | 105 | void Biquad::setLowpassParams(double cutoff, double resonance) |
michael@0 | 106 | { |
michael@0 | 107 | // Limit cutoff to 0 to 1. |
michael@0 | 108 | cutoff = std::max(0.0, std::min(cutoff, 1.0)); |
michael@0 | 109 | |
michael@0 | 110 | if (cutoff == 1) { |
michael@0 | 111 | // When cutoff is 1, the z-transform is 1. |
michael@0 | 112 | setNormalizedCoefficients(1, 0, 0, |
michael@0 | 113 | 1, 0, 0); |
michael@0 | 114 | } else if (cutoff > 0) { |
michael@0 | 115 | // Compute biquad coefficients for lowpass filter |
michael@0 | 116 | resonance = std::max(0.0, resonance); // can't go negative |
michael@0 | 117 | double g = pow(10.0, 0.05 * resonance); |
michael@0 | 118 | double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2); |
michael@0 | 119 | |
michael@0 | 120 | double theta = M_PI * cutoff; |
michael@0 | 121 | double sn = 0.5 * d * sin(theta); |
michael@0 | 122 | double beta = 0.5 * (1 - sn) / (1 + sn); |
michael@0 | 123 | double gamma = (0.5 + beta) * cos(theta); |
michael@0 | 124 | double alpha = 0.25 * (0.5 + beta - gamma); |
michael@0 | 125 | |
michael@0 | 126 | double b0 = 2 * alpha; |
michael@0 | 127 | double b1 = 2 * 2 * alpha; |
michael@0 | 128 | double b2 = 2 * alpha; |
michael@0 | 129 | double a1 = 2 * -gamma; |
michael@0 | 130 | double a2 = 2 * beta; |
michael@0 | 131 | |
michael@0 | 132 | setNormalizedCoefficients(b0, b1, b2, 1, a1, a2); |
michael@0 | 133 | } else { |
michael@0 | 134 | // When cutoff is zero, nothing gets through the filter, so set |
michael@0 | 135 | // coefficients up correctly. |
michael@0 | 136 | setNormalizedCoefficients(0, 0, 0, |
michael@0 | 137 | 1, 0, 0); |
michael@0 | 138 | } |
michael@0 | 139 | } |
michael@0 | 140 | |
michael@0 | 141 | void Biquad::setHighpassParams(double cutoff, double resonance) |
michael@0 | 142 | { |
michael@0 | 143 | // Limit cutoff to 0 to 1. |
michael@0 | 144 | cutoff = std::max(0.0, std::min(cutoff, 1.0)); |
michael@0 | 145 | |
michael@0 | 146 | if (cutoff == 1) { |
michael@0 | 147 | // The z-transform is 0. |
michael@0 | 148 | setNormalizedCoefficients(0, 0, 0, |
michael@0 | 149 | 1, 0, 0); |
michael@0 | 150 | } else if (cutoff > 0) { |
michael@0 | 151 | // Compute biquad coefficients for highpass filter |
michael@0 | 152 | resonance = std::max(0.0, resonance); // can't go negative |
michael@0 | 153 | double g = pow(10.0, 0.05 * resonance); |
michael@0 | 154 | double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2); |
michael@0 | 155 | |
michael@0 | 156 | double theta = M_PI * cutoff; |
michael@0 | 157 | double sn = 0.5 * d * sin(theta); |
michael@0 | 158 | double beta = 0.5 * (1 - sn) / (1 + sn); |
michael@0 | 159 | double gamma = (0.5 + beta) * cos(theta); |
michael@0 | 160 | double alpha = 0.25 * (0.5 + beta + gamma); |
michael@0 | 161 | |
michael@0 | 162 | double b0 = 2 * alpha; |
michael@0 | 163 | double b1 = 2 * -2 * alpha; |
michael@0 | 164 | double b2 = 2 * alpha; |
michael@0 | 165 | double a1 = 2 * -gamma; |
michael@0 | 166 | double a2 = 2 * beta; |
michael@0 | 167 | |
michael@0 | 168 | setNormalizedCoefficients(b0, b1, b2, 1, a1, a2); |
michael@0 | 169 | } else { |
michael@0 | 170 | // When cutoff is zero, we need to be careful because the above |
michael@0 | 171 | // gives a quadratic divided by the same quadratic, with poles |
michael@0 | 172 | // and zeros on the unit circle in the same place. When cutoff |
michael@0 | 173 | // is zero, the z-transform is 1. |
michael@0 | 174 | setNormalizedCoefficients(1, 0, 0, |
michael@0 | 175 | 1, 0, 0); |
michael@0 | 176 | } |
michael@0 | 177 | } |
michael@0 | 178 | |
michael@0 | 179 | void Biquad::setNormalizedCoefficients(double b0, double b1, double b2, double a0, double a1, double a2) |
michael@0 | 180 | { |
michael@0 | 181 | double a0Inverse = 1 / a0; |
michael@0 | 182 | |
michael@0 | 183 | m_b0 = b0 * a0Inverse; |
michael@0 | 184 | m_b1 = b1 * a0Inverse; |
michael@0 | 185 | m_b2 = b2 * a0Inverse; |
michael@0 | 186 | m_a1 = a1 * a0Inverse; |
michael@0 | 187 | m_a2 = a2 * a0Inverse; |
michael@0 | 188 | } |
michael@0 | 189 | |
michael@0 | 190 | void Biquad::setLowShelfParams(double frequency, double dbGain) |
michael@0 | 191 | { |
michael@0 | 192 | // Clip frequencies to between 0 and 1, inclusive. |
michael@0 | 193 | frequency = std::max(0.0, std::min(frequency, 1.0)); |
michael@0 | 194 | |
michael@0 | 195 | double A = pow(10.0, dbGain / 40); |
michael@0 | 196 | |
michael@0 | 197 | if (frequency == 1) { |
michael@0 | 198 | // The z-transform is a constant gain. |
michael@0 | 199 | setNormalizedCoefficients(A * A, 0, 0, |
michael@0 | 200 | 1, 0, 0); |
michael@0 | 201 | } else if (frequency > 0) { |
michael@0 | 202 | double w0 = M_PI * frequency; |
michael@0 | 203 | double S = 1; // filter slope (1 is max value) |
michael@0 | 204 | double alpha = 0.5 * sin(w0) * sqrt((A + 1 / A) * (1 / S - 1) + 2); |
michael@0 | 205 | double k = cos(w0); |
michael@0 | 206 | double k2 = 2 * sqrt(A) * alpha; |
michael@0 | 207 | double aPlusOne = A + 1; |
michael@0 | 208 | double aMinusOne = A - 1; |
michael@0 | 209 | |
michael@0 | 210 | double b0 = A * (aPlusOne - aMinusOne * k + k2); |
michael@0 | 211 | double b1 = 2 * A * (aMinusOne - aPlusOne * k); |
michael@0 | 212 | double b2 = A * (aPlusOne - aMinusOne * k - k2); |
michael@0 | 213 | double a0 = aPlusOne + aMinusOne * k + k2; |
michael@0 | 214 | double a1 = -2 * (aMinusOne + aPlusOne * k); |
michael@0 | 215 | double a2 = aPlusOne + aMinusOne * k - k2; |
michael@0 | 216 | |
michael@0 | 217 | setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); |
michael@0 | 218 | } else { |
michael@0 | 219 | // When frequency is 0, the z-transform is 1. |
michael@0 | 220 | setNormalizedCoefficients(1, 0, 0, |
michael@0 | 221 | 1, 0, 0); |
michael@0 | 222 | } |
michael@0 | 223 | } |
michael@0 | 224 | |
michael@0 | 225 | void Biquad::setHighShelfParams(double frequency, double dbGain) |
michael@0 | 226 | { |
michael@0 | 227 | // Clip frequencies to between 0 and 1, inclusive. |
michael@0 | 228 | frequency = std::max(0.0, std::min(frequency, 1.0)); |
michael@0 | 229 | |
michael@0 | 230 | double A = pow(10.0, dbGain / 40); |
michael@0 | 231 | |
michael@0 | 232 | if (frequency == 1) { |
michael@0 | 233 | // The z-transform is 1. |
michael@0 | 234 | setNormalizedCoefficients(1, 0, 0, |
michael@0 | 235 | 1, 0, 0); |
michael@0 | 236 | } else if (frequency > 0) { |
michael@0 | 237 | double w0 = M_PI * frequency; |
michael@0 | 238 | double S = 1; // filter slope (1 is max value) |
michael@0 | 239 | double alpha = 0.5 * sin(w0) * sqrt((A + 1 / A) * (1 / S - 1) + 2); |
michael@0 | 240 | double k = cos(w0); |
michael@0 | 241 | double k2 = 2 * sqrt(A) * alpha; |
michael@0 | 242 | double aPlusOne = A + 1; |
michael@0 | 243 | double aMinusOne = A - 1; |
michael@0 | 244 | |
michael@0 | 245 | double b0 = A * (aPlusOne + aMinusOne * k + k2); |
michael@0 | 246 | double b1 = -2 * A * (aMinusOne + aPlusOne * k); |
michael@0 | 247 | double b2 = A * (aPlusOne + aMinusOne * k - k2); |
michael@0 | 248 | double a0 = aPlusOne - aMinusOne * k + k2; |
michael@0 | 249 | double a1 = 2 * (aMinusOne - aPlusOne * k); |
michael@0 | 250 | double a2 = aPlusOne - aMinusOne * k - k2; |
michael@0 | 251 | |
michael@0 | 252 | setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); |
michael@0 | 253 | } else { |
michael@0 | 254 | // When frequency = 0, the filter is just a gain, A^2. |
michael@0 | 255 | setNormalizedCoefficients(A * A, 0, 0, |
michael@0 | 256 | 1, 0, 0); |
michael@0 | 257 | } |
michael@0 | 258 | } |
michael@0 | 259 | |
michael@0 | 260 | void Biquad::setPeakingParams(double frequency, double Q, double dbGain) |
michael@0 | 261 | { |
michael@0 | 262 | // Clip frequencies to between 0 and 1, inclusive. |
michael@0 | 263 | frequency = std::max(0.0, std::min(frequency, 1.0)); |
michael@0 | 264 | |
michael@0 | 265 | // Don't let Q go negative, which causes an unstable filter. |
michael@0 | 266 | Q = std::max(0.0, Q); |
michael@0 | 267 | |
michael@0 | 268 | double A = pow(10.0, dbGain / 40); |
michael@0 | 269 | |
michael@0 | 270 | if (frequency > 0 && frequency < 1) { |
michael@0 | 271 | if (Q > 0) { |
michael@0 | 272 | double w0 = M_PI * frequency; |
michael@0 | 273 | double alpha = sin(w0) / (2 * Q); |
michael@0 | 274 | double k = cos(w0); |
michael@0 | 275 | |
michael@0 | 276 | double b0 = 1 + alpha * A; |
michael@0 | 277 | double b1 = -2 * k; |
michael@0 | 278 | double b2 = 1 - alpha * A; |
michael@0 | 279 | double a0 = 1 + alpha / A; |
michael@0 | 280 | double a1 = -2 * k; |
michael@0 | 281 | double a2 = 1 - alpha / A; |
michael@0 | 282 | |
michael@0 | 283 | setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); |
michael@0 | 284 | } else { |
michael@0 | 285 | // When Q = 0, the above formulas have problems. If we look at |
michael@0 | 286 | // the z-transform, we can see that the limit as Q->0 is A^2, so |
michael@0 | 287 | // set the filter that way. |
michael@0 | 288 | setNormalizedCoefficients(A * A, 0, 0, |
michael@0 | 289 | 1, 0, 0); |
michael@0 | 290 | } |
michael@0 | 291 | } else { |
michael@0 | 292 | // When frequency is 0 or 1, the z-transform is 1. |
michael@0 | 293 | setNormalizedCoefficients(1, 0, 0, |
michael@0 | 294 | 1, 0, 0); |
michael@0 | 295 | } |
michael@0 | 296 | } |
michael@0 | 297 | |
michael@0 | 298 | void Biquad::setAllpassParams(double frequency, double Q) |
michael@0 | 299 | { |
michael@0 | 300 | // Clip frequencies to between 0 and 1, inclusive. |
michael@0 | 301 | frequency = std::max(0.0, std::min(frequency, 1.0)); |
michael@0 | 302 | |
michael@0 | 303 | // Don't let Q go negative, which causes an unstable filter. |
michael@0 | 304 | Q = std::max(0.0, Q); |
michael@0 | 305 | |
michael@0 | 306 | if (frequency > 0 && frequency < 1) { |
michael@0 | 307 | if (Q > 0) { |
michael@0 | 308 | double w0 = M_PI * frequency; |
michael@0 | 309 | double alpha = sin(w0) / (2 * Q); |
michael@0 | 310 | double k = cos(w0); |
michael@0 | 311 | |
michael@0 | 312 | double b0 = 1 - alpha; |
michael@0 | 313 | double b1 = -2 * k; |
michael@0 | 314 | double b2 = 1 + alpha; |
michael@0 | 315 | double a0 = 1 + alpha; |
michael@0 | 316 | double a1 = -2 * k; |
michael@0 | 317 | double a2 = 1 - alpha; |
michael@0 | 318 | |
michael@0 | 319 | setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); |
michael@0 | 320 | } else { |
michael@0 | 321 | // When Q = 0, the above formulas have problems. If we look at |
michael@0 | 322 | // the z-transform, we can see that the limit as Q->0 is -1, so |
michael@0 | 323 | // set the filter that way. |
michael@0 | 324 | setNormalizedCoefficients(-1, 0, 0, |
michael@0 | 325 | 1, 0, 0); |
michael@0 | 326 | } |
michael@0 | 327 | } else { |
michael@0 | 328 | // When frequency is 0 or 1, the z-transform is 1. |
michael@0 | 329 | setNormalizedCoefficients(1, 0, 0, |
michael@0 | 330 | 1, 0, 0); |
michael@0 | 331 | } |
michael@0 | 332 | } |
michael@0 | 333 | |
michael@0 | 334 | void Biquad::setNotchParams(double frequency, double Q) |
michael@0 | 335 | { |
michael@0 | 336 | // Clip frequencies to between 0 and 1, inclusive. |
michael@0 | 337 | frequency = std::max(0.0, std::min(frequency, 1.0)); |
michael@0 | 338 | |
michael@0 | 339 | // Don't let Q go negative, which causes an unstable filter. |
michael@0 | 340 | Q = std::max(0.0, Q); |
michael@0 | 341 | |
michael@0 | 342 | if (frequency > 0 && frequency < 1) { |
michael@0 | 343 | if (Q > 0) { |
michael@0 | 344 | double w0 = M_PI * frequency; |
michael@0 | 345 | double alpha = sin(w0) / (2 * Q); |
michael@0 | 346 | double k = cos(w0); |
michael@0 | 347 | |
michael@0 | 348 | double b0 = 1; |
michael@0 | 349 | double b1 = -2 * k; |
michael@0 | 350 | double b2 = 1; |
michael@0 | 351 | double a0 = 1 + alpha; |
michael@0 | 352 | double a1 = -2 * k; |
michael@0 | 353 | double a2 = 1 - alpha; |
michael@0 | 354 | |
michael@0 | 355 | setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); |
michael@0 | 356 | } else { |
michael@0 | 357 | // When Q = 0, the above formulas have problems. If we look at |
michael@0 | 358 | // the z-transform, we can see that the limit as Q->0 is 0, so |
michael@0 | 359 | // set the filter that way. |
michael@0 | 360 | setNormalizedCoefficients(0, 0, 0, |
michael@0 | 361 | 1, 0, 0); |
michael@0 | 362 | } |
michael@0 | 363 | } else { |
michael@0 | 364 | // When frequency is 0 or 1, the z-transform is 1. |
michael@0 | 365 | setNormalizedCoefficients(1, 0, 0, |
michael@0 | 366 | 1, 0, 0); |
michael@0 | 367 | } |
michael@0 | 368 | } |
michael@0 | 369 | |
michael@0 | 370 | void Biquad::setBandpassParams(double frequency, double Q) |
michael@0 | 371 | { |
michael@0 | 372 | // No negative frequencies allowed. |
michael@0 | 373 | frequency = std::max(0.0, frequency); |
michael@0 | 374 | |
michael@0 | 375 | // Don't let Q go negative, which causes an unstable filter. |
michael@0 | 376 | Q = std::max(0.0, Q); |
michael@0 | 377 | |
michael@0 | 378 | if (frequency > 0 && frequency < 1) { |
michael@0 | 379 | double w0 = M_PI * frequency; |
michael@0 | 380 | if (Q > 0) { |
michael@0 | 381 | double alpha = sin(w0) / (2 * Q); |
michael@0 | 382 | double k = cos(w0); |
michael@0 | 383 | |
michael@0 | 384 | double b0 = alpha; |
michael@0 | 385 | double b1 = 0; |
michael@0 | 386 | double b2 = -alpha; |
michael@0 | 387 | double a0 = 1 + alpha; |
michael@0 | 388 | double a1 = -2 * k; |
michael@0 | 389 | double a2 = 1 - alpha; |
michael@0 | 390 | |
michael@0 | 391 | setNormalizedCoefficients(b0, b1, b2, a0, a1, a2); |
michael@0 | 392 | } else { |
michael@0 | 393 | // When Q = 0, the above formulas have problems. If we look at |
michael@0 | 394 | // the z-transform, we can see that the limit as Q->0 is 1, so |
michael@0 | 395 | // set the filter that way. |
michael@0 | 396 | setNormalizedCoefficients(1, 0, 0, |
michael@0 | 397 | 1, 0, 0); |
michael@0 | 398 | } |
michael@0 | 399 | } else { |
michael@0 | 400 | // When the cutoff is zero, the z-transform approaches 0, if Q |
michael@0 | 401 | // > 0. When both Q and cutoff are zero, the z-transform is |
michael@0 | 402 | // pretty much undefined. What should we do in this case? |
michael@0 | 403 | // For now, just make the filter 0. When the cutoff is 1, the |
michael@0 | 404 | // z-transform also approaches 0. |
michael@0 | 405 | setNormalizedCoefficients(0, 0, 0, |
michael@0 | 406 | 1, 0, 0); |
michael@0 | 407 | } |
michael@0 | 408 | } |
michael@0 | 409 | |
michael@0 | 410 | void Biquad::setZeroPolePairs(const Complex &zero, const Complex &pole) |
michael@0 | 411 | { |
michael@0 | 412 | double b0 = 1; |
michael@0 | 413 | double b1 = -2 * zero.real(); |
michael@0 | 414 | |
michael@0 | 415 | double zeroMag = abs(zero); |
michael@0 | 416 | double b2 = zeroMag * zeroMag; |
michael@0 | 417 | |
michael@0 | 418 | double a1 = -2 * pole.real(); |
michael@0 | 419 | |
michael@0 | 420 | double poleMag = abs(pole); |
michael@0 | 421 | double a2 = poleMag * poleMag; |
michael@0 | 422 | setNormalizedCoefficients(b0, b1, b2, 1, a1, a2); |
michael@0 | 423 | } |
michael@0 | 424 | |
michael@0 | 425 | void Biquad::setAllpassPole(const Complex &pole) |
michael@0 | 426 | { |
michael@0 | 427 | Complex zero = Complex(1, 0) / pole; |
michael@0 | 428 | setZeroPolePairs(zero, pole); |
michael@0 | 429 | } |
michael@0 | 430 | |
michael@0 | 431 | void Biquad::getFrequencyResponse(int nFrequencies, |
michael@0 | 432 | const float* frequency, |
michael@0 | 433 | float* magResponse, |
michael@0 | 434 | float* phaseResponse) |
michael@0 | 435 | { |
michael@0 | 436 | // Evaluate the Z-transform of the filter at given normalized |
michael@0 | 437 | // frequency from 0 to 1. (1 corresponds to the Nyquist |
michael@0 | 438 | // frequency.) |
michael@0 | 439 | // |
michael@0 | 440 | // The z-transform of the filter is |
michael@0 | 441 | // |
michael@0 | 442 | // H(z) = (b0 + b1*z^(-1) + b2*z^(-2))/(1 + a1*z^(-1) + a2*z^(-2)) |
michael@0 | 443 | // |
michael@0 | 444 | // Evaluate as |
michael@0 | 445 | // |
michael@0 | 446 | // b0 + (b1 + b2*z1)*z1 |
michael@0 | 447 | // -------------------- |
michael@0 | 448 | // 1 + (a1 + a2*z1)*z1 |
michael@0 | 449 | // |
michael@0 | 450 | // with z1 = 1/z and z = exp(j*pi*frequency). Hence z1 = exp(-j*pi*frequency) |
michael@0 | 451 | |
michael@0 | 452 | // Make local copies of the coefficients as a micro-optimization. |
michael@0 | 453 | double b0 = m_b0; |
michael@0 | 454 | double b1 = m_b1; |
michael@0 | 455 | double b2 = m_b2; |
michael@0 | 456 | double a1 = m_a1; |
michael@0 | 457 | double a2 = m_a2; |
michael@0 | 458 | |
michael@0 | 459 | for (int k = 0; k < nFrequencies; ++k) { |
michael@0 | 460 | double omega = -M_PI * frequency[k]; |
michael@0 | 461 | Complex z = Complex(cos(omega), sin(omega)); |
michael@0 | 462 | Complex numerator = b0 + (b1 + b2 * z) * z; |
michael@0 | 463 | Complex denominator = Complex(1, 0) + (a1 + a2 * z) * z; |
michael@0 | 464 | Complex response = numerator / denominator; |
michael@0 | 465 | magResponse[k] = static_cast<float>(abs(response)); |
michael@0 | 466 | phaseResponse[k] = static_cast<float>(atan2(imag(response), real(response))); |
michael@0 | 467 | } |
michael@0 | 468 | } |
michael@0 | 469 | |
michael@0 | 470 | } // namespace WebCore |
michael@0 | 471 |