Implement Biquad lowpass and highpass filters according to spec

third_party/WebKit/Source/platform/audio/Biquad.cpp

 /*
  * Copyright (C) 2010 Google Inc. All rights reserved.
  *
  * Redistribution and use in source and binary forms, with or without
  * modification, are permitted provided that the following conditions
  * are met:
  *
  * 1.  Redistributions of source code must retain the above copyright
  *     notice, this list of conditions and the following disclaimer.
  * 2.  Redistributions in binary form must reproduce the above copyright
@@ skipping 237 matching lines @@
     // Limit cutoff to 0 to 1.
     cutoff = clampTo(cutoff, 0.0, 1.0);
 
     if (cutoff == 1) {
         // When cutoff is 1, the z-transform is 1.
         setNormalizedCoefficients(index,
             1, 0, 0,
             1, 0, 0);
     } else if (cutoff > 0) {
         // Compute biquad coefficients for lowpass filter
-        resonance = std::max(0.0, resonance); // can't go negative
-        double g = pow(10.0, 0.05 * resonance);
-        double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2);
 
+        resonance = pow(10, resonance / 20);
         double theta = piDouble * cutoff;
-        double sn = 0.5 * d * sin(theta);
-        double beta = 0.5 * (1 - sn) / (1 + sn);
-        double gamma = (0.5 + beta) * cos(theta);
-        double alpha = 0.25 * (0.5 + beta - gamma);
+        double alpha = sin(theta) / (2 * resonance);
+        double cosw = cos(theta);
+        double beta = (1 - cosw) / 2;
 
-        double b0 = 2 * alpha;
-        double b1 = 2 * 2 * alpha;
-        double b2 = 2 * alpha;
-        double a1 = 2 * -gamma;
-        double a2 = 2 * beta;
+        double b0 = beta;
+        double b1 = 2 * beta;
+        double b2 = beta;
 
-        setNormalizedCoefficients(index, b0, b1, b2, 1, a1, a2);
+        double a0 = 1 + alpha;
+        double a1 = -2 * cosw;
+        double a2 = 1 - alpha;
+
+        setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
     } else {
         // When cutoff is zero, nothing gets through the filter, so set
         // coefficients up correctly.
         setNormalizedCoefficients(index,
             0, 0, 0,
             1, 0, 0);
     }
 }
 
 void Biquad::setHighpassParams(int index, double cutoff, double resonance)
 {
     // Limit cutoff to 0 to 1.
     cutoff = clampTo(cutoff, 0.0, 1.0);
 
     if (cutoff == 1) {
         // The z-transform is 0.
         setNormalizedCoefficients(index,
             0, 0, 0,
             1, 0, 0);
     } else if (cutoff > 0) {
         // Compute biquad coefficients for highpass filter
-        resonance = std::max(0.0, resonance); // can't go negative
-        double g = pow(10.0, 0.05 * resonance);
-        double d = sqrt((4 - sqrt(16 - 16 / (g * g))) / 2);
 
+        resonance = pow(10, resonance / 20);
         double theta = piDouble * cutoff;
-        double sn = 0.5 * d * sin(theta);
-        double beta = 0.5 * (1 - sn) / (1 + sn);
-        double gamma = (0.5 + beta) * cos(theta);
-        double alpha = 0.25 * (0.5 + beta + gamma);
+        double alpha = sin(theta) / (2 * resonance);
+        double cosw = cos(theta);
+        double beta = (1 + cosw) / 2;
 
-        double b0 = 2 * alpha;
-        double b1 = 2 * -2 * alpha;
-        double b2 = 2 * alpha;
-        double a1 = 2 * -gamma;
-        double a2 = 2 * beta;
+        double b0 = beta;
+        double b1 = -2 * beta;
+        double b2 = beta;
 
-        setNormalizedCoefficients(index, b0, b1, b2, 1, a1, a2);
+        double a0 = 1 + alpha;
+        double a1 = -2 * cosw;
+        double a2 = 1 - alpha;
+
+        setNormalizedCoefficients(index, b0, b1, b2, a0, a1, a2);
     } else {
         // When cutoff is zero, we need to be careful because the above
         // gives a quadratic divided by the same quadratic, with poles
         // and zeros on the unit circle in the same place. When cutoff
         // is zero, the z-transform is 1.
         setNormalizedCoefficients(index,
             1, 0, 0,
             1, 0, 0);
     }
 }
@@ skipping 286 matching lines @@
         std::complex<double> denominator =
             std::complex<double>(1, 0) + (a1 + a2 * z) * z;
         std::complex<double> response = numerator / denominator;
         magResponse[k] = static_cast<float>(abs(response));
         phaseResponse[k] = static_cast<float>(atan2(imag(response), real(response)));
     }
 }
 
 } // namespace blink