Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

Emulation of Diode ladder lowpass filter

By: a guest on Apr 5th, 2012  |  syntax: C++  |  size: 3.63 KB  |  views: 620  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
This paste has a previous version, view the difference. Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1.  
  2. // This code is released under the MIT license (see below).
  3. //
  4. // The MIT License
  5. //
  6. // Copyright (c) 2012 Dominique Wurtz (www.blaukraut.info)
  7. //
  8. // Permission is hereby granted, free of charge, to any person obtaining a copy
  9. // of this software and associated documentation files (the "Software"), to deal
  10. // in the Software without restriction, including without limitation the rights
  11. // to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  12. // copies of the Software, and to permit persons to whom the Software is
  13. // furnished to do so, subject to the following conditions:
  14. //
  15. // The above copyright notice and this permission notice shall be included in
  16. // all copies or substantial portions of the Software.
  17. //
  18. // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  19. // IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  20. // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  21. // AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  22. // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  23. // OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
  24. // THE SOFTWARE.
  25.  
  26. #ifndef __DIODE_LADDER_FILTER_HPP__
  27. #define __DIODE_LADDER_FILTER_HPP__
  28.  
  29. #include <cmath>
  30. #include <algorithm>
  31.  
  32. // Emulation of Diode ladder lowpass filter as found in Roland TB303 or EMS VCS3
  33. //
  34. // Version 0.2 (04/05/2012) greatly simplified equations; add highpass filter in feedback path
  35. // Version 0.1 (04/03/2012) initial version
  36.  
  37. class DiodeLadderFilter
  38. {
  39. public:
  40.  
  41.         DiodeLadderFilter()
  42.         {
  43.                 std::fill(z, z + 5, 0);
  44.                 set_q(0);
  45.         }
  46.  
  47.         // fc: normalized cutoff frequency in the range [0..1] => 0 HZ .. Nyquist
  48.         void set_feedback_hpf_cutoff(const double fc)
  49.         {
  50.                 const double K = fc * M_PI;
  51.                 ah = (K - 2) / (K + 2);
  52.                 bh = 2 / (K + 2);
  53.         }
  54.  
  55.         void reset()
  56.         {
  57.                 if (k < 17) std::fill(z, z + 5, 0);
  58.         }
  59.  
  60.         // q: resonance in the range [0..1]
  61.         void set_q(const double q)
  62.         {
  63.                 assert(q >= 0 && q <= 1.);
  64.                 k = 20 * q;
  65.                 A = 1 + 0.5*k; // resonance gain compensation
  66.         }
  67.  
  68.         // Process one sample.
  69.         //
  70.         // x: input signal
  71.         // fc: normalized cutoff frequency in the range [0..1] => 0 HZ .. Nyquist
  72.         __forceinline double tick(const double x, const double fc)
  73.         {
  74.                 assert(fc > 0 && fc < 1);
  75.                 const double a = M_PI * fc; // PI is Nyquist frequency
  76.                 // a = 2 * tan(0.5*a); // dewarping, not required with 2x oversampling
  77.                 const double ainv = 1/a;
  78.                 const double a2 = a*a;
  79.                 const double b = 2*a + 1;
  80.                 const double b2 = b*b;
  81.                 const double c = 1 / (2*a2*a2 - 4*a2*b2 + b2*b2);
  82.                 const double g0 = 2*a2*a2*c;
  83.                 const double g = g0 * bh;
  84.  
  85.                 // current state
  86.                 const double s0 = (a2*a*z[0] + a2*b*z[1] + z[2]*(b2 - 2*a2)*a + z[3]*(b2 - 3*a2)*b) * c;
  87.                 const double s = bh*s0 - z[4];
  88.  
  89.                 // solve feedback loop (linear)
  90.                 double y5 = (g*x + s) / (1 + g*k);
  91.  
  92.                 // input clipping
  93.                 const double y0 = clip(x - k*y5);
  94.                 y5 = g*y0 + s;
  95.  
  96.                 // compute integrator outputs
  97.                 const double y4 = g0*y0 + s0;
  98.                 const double y3 = (b*y4 - z[3]) * ainv;
  99.                 const double y2 = (b*y3 - a*y4 - z[2]) * ainv;
  100.                 const double y1 = (b*y2 - a*y3 - z[1]) * ainv;
  101.  
  102.                 // update filter state
  103.                 z[0] += 4*a*(y0 - y1 + y2);
  104.                 z[1] += 2*a*(y1 - 2*y2 + y3);
  105.                 z[2] += 2*a*(y2 - 2*y3 + y4);
  106.                 z[3] += 2*a*(y3 - 2*y4);
  107.                 z[4] = bh*y4 + ah*y5;
  108.  
  109.                 return A*y4;
  110.         }
  111.        
  112. private:
  113.         double k, A;
  114.         double z[5]; // filter memory (4 integrators plus 1st order HPF)
  115.         double ah, bh; // feedback HPF coeffs
  116.  
  117.         static __forceinline double clip(const double x)
  118.         {
  119.                 return x / (1 + abs(x));
  120.         }
  121. };
  122.  
  123. #endif // __DIODE_LADDER_FILTER_HPP__