Guest User

diode_ladder_filter.hpp

a guest
Apr 3rd, 2012
1,535
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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. // Version 0.1 (04/03/2012)
  34.  
  35. class DiodeLadderFilter
  36. {
  37. public:
  38.  
  39.     DiodeLadderFilter()
  40.     {
  41.         std::fill(z, z + 4, 0);
  42.         set_q(0);
  43.     }
  44.  
  45.     void reset()
  46.     {
  47.         if (k < 17) std::fill(z, z + 4, 0);
  48.     }
  49.  
  50.     // q: resonance in the range [0..1]
  51.     void set_q(const double q)
  52.     {
  53.         assert(q >= 0 && q <= 1.);
  54.         k = 20 * q;
  55.         A = 1 + 0.5*k; // resonance gain compensation
  56.     }
  57.  
  58.     // Process one sample.
  59.     //
  60.     // x: input signal
  61.     // fc: normalized cutoff frequency in the range [0..1] => 0 HZ .. Nyquist
  62.     __forceinline double tick(const double x, const double fc)
  63.     {
  64.         assert(fc > 0 && fc < 1);
  65.         const double wc = PI_HALF * fc; // PI is Nyquist frequency
  66.         // wc = 2 * tan(0.5*wc); // dewarping, not required with 2x oversampling
  67.         const double wc2 = wc*wc;
  68.         const double wc3 = wc2*wc;
  69.         const double wc4 = wc3*wc;
  70.         const double b = 1 / (1+8*wc+20*wc2+16*wc3+2*wc4);
  71.         const double g = 2*wc4 * b;
  72.  
  73.         // current state
  74.         const double s = (z[0]*wc3 + z[1]*(wc2+2*wc3) + z[2]*(wc+4*wc2+2*wc3) + z[3]*(1+6*wc+9*wc2+2*wc3)) * b;
  75.        
  76.         // solve feedback loop (linear)
  77.         double y4 = (g*x + s) / (1 + g*k);
  78.  
  79.         // input clipping
  80.         const double y0 = fast_tanh(x - k*y4);
  81.  
  82.         // Compute all integrator outputs (y1, y2, y3, y4).
  83.         // Unlike in the well-known Moog transistor ladder, this gets quite nasty due the
  84.         // inherent coupling between filter stages.
  85.         const double y1 = (y0*(2*wc+12*wc2+20*wc3+8*wc4) + z[0]*(1+6*wc+10*wc2+4*wc3) +
  86.             z[1]*(2*wc+8*wc2+6*wc3) + z[2]*(2*wc2+4*wc3) + z[3]*2*wc3)*b;
  87.         const double y2 = (y0*(2*wc2+8*wc3+6*wc4) + z[0]*(wc+4*wc2+3*wc3) +
  88.             z[1]*(1+6*wc+11*wc2+6*wc3) + z[2]*(wc+4*wc2+4*wc3) + z[3]*(wc2+2*wc3))*b;
  89.         const double y3 = (y0*(2*wc3+4*wc4) + z[0]*(wc2+2*wc3) +
  90.             z[1]*(wc+4*wc2+4*wc3) + z[2]*(1+6*wc+10*wc2+4*wc3) + z[3]*(wc+4*wc2+2*wc3))*b;
  91.         y4 = g*y0 + s;
  92.  
  93.         // update filter state
  94.         z[0] += 4*wc*(y0 - y1 + y2);
  95.         z[1] += 2*wc*(y1 - 2*y2 + y3);
  96.         z[2] += 2*wc*(y2 - 2*y3 + y4);
  97.         z[3] += 2*wc*(y3 - 2*y4);
  98.  
  99.         return A*y4;
  100.     }
  101.    
  102. private:
  103.     double k, A;
  104.     double z[4];
  105.  
  106.     static __forceinline double fast_tanh(const double x)
  107.     {
  108.         return x / (1 + abs(x));
  109.     }
  110. };
  111.  
  112. #endif // __DIODE_LADDER_FILTER_HPP__
RAW Paste Data