zynaddsubfx

ZynAddSubFX open source synthesizer
Log | Files | Refs | Submodules | LICENSE

MoogFilter.cpp (4502B)

      1 #include <cassert>
      2 #include <cstdio>
      3 #include <cmath>
      4 #include <stdio.h>
      5 
      6 #include "../Misc/Util.h"
      7 #include "MoogFilter.h"
      8 
      9 // theory from "THE ART OF VA FILTER DESIGN"
     10 // by Vadim Zavalishin
     11 
     12 namespace zyn{
     13 
     14 MoogFilter::MoogFilter(unsigned char Ftype, float Ffreq, float Fq,
     15     unsigned int srate, int bufsize)
     16     :Filter(srate, bufsize), sr(srate), gain(1.0f)
     17 {
     18     setfreq_and_q(Ffreq/srate, Fq);
     19     settype(Ftype); // q must be set before
     20 
     21     // initialize state but not to 0, in case of being used as denominator
     22     for (unsigned int i = 0; i<(sizeof(state)/sizeof(*state)); i++)
     23     {
     24         state[i] = 0.0001f;
     25     }
     26 }
     27 
     28 MoogFilter::~MoogFilter(void)
     29 {
     30 
     31 }
     32 
     33 inline float MoogFilter::tanhX(const float x) const
     34 {
     35     // Pade approximation of tanh(x) bound to [-1 .. +1]
     36     // https://mathr.co.uk/blog/2017-09-06_approximating_hyperbolic_tangent.html
     37     float x2 = x*x;
     38     return (x*(105.0f+10.0f*x2)/(105.0f+(45.0f+x2)*x2));
     39 }
     40 
     41 
     42 inline float MoogFilter::tanhXdivX(float x) const
     43 {
     44 
     45     //add DC offset to raise even harmonics (like transistor bias current)
     46     x+= 0.1f;
     47     // pre calc often used term
     48     float x2 = x*x;
     49     // Pade approximation for tanh(x)/x used in filter stages (5x per sample)
     50     //~ return ((15.0+x2)/(15.0+6.0*x2));
     51     // faster approximation without division
     52     // tuned for more distortion in self oscillation
     53     return (1.0f-(0.35f*x2)+(0.06f*x2*x2));
     54 }
     55 
     56 inline float MoogFilter::step(float input)
     57 {
     58     // transconductance
     59     // gM(vIn) = tanh( vIn ) / vIn
     60     float gm0 = tanhXdivX(state[0]);
     61     // to ease calculations the others are 1.0 (so daturation)
     62     // denominators
     63     float d0 = 1.0f / (1.0f + c*gm0);
     64     float dp1 = 1.0f / (1.0f + c);
     65     float dp2 = dp1*dp1;
     66     float dp3 = dp2*dp1;
     67 
     68     // instantaneous response estimate
     69     float y3Estimate =
     70         cp4 * d0 * gm0 * dp3 * input +
     71         cp3 * gm0 * dp3     * d0 * state[0] +
     72         cp2 * dp3 * state[1] +
     73         c   * dp2 * state[2] +
     74               dp1 * state[3];
     75 
     76     // mix input and gained feedback estimate for
     77     // cheaper feedback gain compensation. Idea from
     78     // Antti Huovilainen and Vesa Välimäki, "NEW APPROACHES TO DIGITAL SUBTRACTIVE SYNTHESIS"
     79     float u = input - tanhX(feedbackGain * (y3Estimate - 0.5f*input));
     80     // output of all stages
     81     float y0 = gm0 * d0 * (state[0] + c * u);
     82     float y1 = dp1 * (state[1] + c * y0);
     83     float y2 = dp1 * (state[2] + c * y1);
     84     float y3 = dp1 * (state[3] + c * y2);
     85 
     86     // update state
     87     state[0] += ct2 * (u - y0);
     88     state[1] += ct2 * (y0 - y1);
     89     state[2] += ct2 * (y1 - y2);
     90     state[3] += ct2 * (y2 - y3);
     91 
     92     // calculate multimode filter output
     93     return (a0 * u
     94           + a1 * y0
     95           + a2 * y1
     96           + a3 * y2
     97           + a4 * y3);
     98 }
     99 
    100 void MoogFilter::filterout(float *smp)
    101 {
    102     for (int i = 0; i < buffersize; i ++)
    103     {
    104         smp[i] = step(tanhX(smp[i]*gain));
    105         smp[i] *= outgain;
    106     }
    107 }
    108 
    109 void MoogFilter::setfreq_and_q(float frequency, float q_)
    110 {
    111     setfreq(frequency/sr);
    112     setq(q_);
    113 }
    114 
    115 inline float MoogFilter::tan_2(const float x) const
    116 {
    117     //Pade approximation tan(x) hand tuned to map fCutoff
    118     float x2 = x*x;
    119     //~ return ((9.54f*x*((11.08f - x2)))/(105.0f - x2*(45.0f + x2))); // more accurate but instable at high frequencies
    120     return (x+0.15f*x2+0.3f*x2*x2); // faster, no division (by zero)
    121 }
    122 
    123 void MoogFilter::setfreq(float ff)
    124 {
    125     // pre warp cutoff to map to reality
    126     c = tan_2(PI * ff);
    127     // limit cutoff to prevent overflow
    128     c = limit(c,0.0006f,1.5f);
    129     // pre calculate some stuff outside the hot zone
    130     ct2 = c * 2.0f;
    131     cp2 = c * c;
    132     cp3 = cp2 * c;
    133     cp4 = cp2 * cp2;
    134 }
    135 
    136 void MoogFilter::setq(float q)
    137 {
    138     // flattening the Q input
    139     // self oscillation begins around 4.0
    140     // mapped to match the ANALOG filter class
    141     feedbackGain = cbrtf(q/1000.0f)*4.0f + 0.3f;
    142     // compensation factor for passband reduction by the negative feedback
    143     passbandCompensation = 1.0f + limit(feedbackGain, 0.0f, 1.0f);
    144 }
    145 
    146 void MoogFilter::setgain(float dBgain)
    147 {
    148     gain = dB2rap(dBgain);
    149 }
    150 
    151 void MoogFilter::settype(unsigned char ftype)
    152 {
    153     switch (ftype)
    154     {
    155         case 0:
    156             a0 = 1.0f; a1 =-4.0f; a2 = 6.0f; a3 =-4.0f; a4 = 1.0f;
    157             break;
    158         case 1:
    159             a0 = 0.0f; a1 = 0.0f; a2 = 4.0f; a3 =-8.0f; a4 = 4.0f;
    160             break;
    161         case 2:
    162         default:
    163             a0 = 0.0f; a1 = 0.0f; a2 = 0.0f; a3 = 0.0f; a4 = passbandCompensation;
    164             break;
    165     }
    166 }
    167 
    168 };