src/resid-fp/filter.h

//  This file is part of reSID, a MOS6581 SID emulator engine.
//  Copyright (C) 2004  Dag Lem <resid@nimrod.no>
//
//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software
//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//  ---------------------------------------------------------------------------
//  Filter distortion code written by Antti S. Lankila 2007 - 2008.

#ifndef __FILTER_H__
#define __FILTER_H__

#include <math.h>
#include "siddefs-fp.h"

// ----------------------------------------------------------------------------
// The SID filter is modeled with a two-integrator-loop biquadratic filter,
// which has been confirmed by Bob Yannes to be the actual circuit used in
// the SID chip.
//
// Measurements show that excellent emulation of the SID filter is achieved,
// except when high resonance is combined with high sustain levels.
// In this case the SID op-amps are performing less than ideally and are
// causing some peculiar behavior of the SID filter. This however seems to
// have more effect on the overall amplitude than on the color of the sound.
//
// The theory for the filter circuit can be found in "Microelectric Circuits"
// by Adel S. Sedra and Kenneth C. Smith.
// The circuit is modeled based on the explanation found there except that
// an additional inverter is used in the feedback from the bandpass output,
// allowing the summer op-amp to operate in single-ended mode. This yields
// inverted filter outputs with levels independent of Q, which corresponds with
// the results obtained from a real SID.
//
// We have been able to model the summer and the two integrators of the circuit
// to form components of an IIR filter.
// Vhp is the output of the summer, Vbp is the output of the first integrator,
// and Vlp is the output of the second integrator in the filter circuit.
//
// According to Bob Yannes, the active stages of the SID filter are not really
// op-amps. Rather, simple NMOS inverters are used. By biasing an inverter
// into its region of quasi-linear operation using a feedback resistor from
// input to output, a MOS inverter can be made to act like an op-amp for
// small signals centered around the switching threshold.
//
// Qualified guesses at SID filter schematics are depicted below.
//
// SID filter
// ----------
//
//     -----------------------------------------------
//    |                                               |
//    |            ---Rq--                            |
//    |           |       |                           |
//    |  ------------<A]-----R1---------              |
//    | |                               |             |
//    | |                        ---C---|      ---C---|
//    | |                       |       |     |       |
//    |  --R1--    ---R1--      |---Rs--|     |---Rs--|
//    |        |  |       |     |       |     |       |
//     ----R1--|-----[A>--|--R-----[A>--|--R-----[A>--|
//             |          |             |             |
// vi -----R1--           |             |             |
//
//                       vhp           vbp           vlp
//
//
// vi  - input voltage
// vhp - highpass output
// vbp - bandpass output
// vlp - lowpass output
// [A> - op-amp
// R1  - summer resistor
// Rq  - resistor array controlling resonance (4 resistors)
// R   - NMOS FET voltage controlled resistor controlling cutoff frequency
// Rs  - shunt resitor
// C   - capacitor
//
//
//
// SID integrator
// --------------
//
//                                   V+
//
//                                   |
//                                   |
//                              -----|
//                             |     |
//                             | ||--
//                              -||
//                   ---C---     ||->
//                  |       |        |
//                  |---Rs-----------|---- vo
//                  |                |
//                  |            ||--
// vi ----     -----|------------||
//        |   ^     |            ||->
//        |___|     |                |
//        -----     |                |
//          |       |                |
//          |---R2--                 |
//          |
//          R1                       V-
//          |
//          |
//
//          Vw
//
// ----------------------------------------------------------------------------
class FilterFP
{
public:
  FilterFP();

  void enable_filter(bool enable);
  void set_chip_model(chip_model model);
  void set_distortion_properties(float, float, float);
  void set_type3_properties(float, float, float, float);
  void set_type4_properties(float, float);
  void set_clock_frequency(float);

  RESID_INLINE
  float clock(float voice1, float voice2, float voice3,
              float ext_in);
  void reset();

  // Write registers.
  void writeFC_LO(reg8);
  void writeFC_HI(reg8);
  void writeRES_FILT(reg8);
  void writeMODE_VOL(reg8);

private:
  void set_Q();
  void set_w0();
  float type3_w0(const float source, const float offset);
  float type4_w0();
  void calculate_helpers();
  void nuke_denormals();

  // Filter enabled.
  bool enabled;

  // 6581/8580 filter model (XXX: we should specialize in separate classes)
  chip_model model;

  // Filter cutoff frequency.
  reg12 fc;

  // Filter resonance.
  reg8 res;

  // Selects which inputs to route through filter.
  reg8 filt;

  // Switch voice 3 off.
  reg8 voice3off;

  // Highpass, bandpass, and lowpass filter modes.
  reg8 hp_bp_lp;

  // Output master volume.
  reg4 vol;
  float volf; /* avoid integer-to-float conversion at output */

  // clock
  float clock_frequency;

  /* Distortion params for Type3 */
  float distortion_rate, distortion_point, distortion_cf_threshold;

  /* Type3 params. */
  float type3_baseresistance, type3_offset, type3_steepness, type3_minimumfetresistance;

  /* Type4 params */
  float type4_k, type4_b;

  // State of filter.
  float Vhp, Vbp, Vlp;

  /* Resonance/Distortion/Type3/Type4 helpers. */
  float type4_w0_cache, _1_div_Q, type3_fc_kink_exp, distortion_CT,
        type3_fc_distortion_offset_bp, type3_fc_distortion_offset_hp;

friend class SIDFP;
};

// ----------------------------------------------------------------------------
// Inline functions.
// The following functions are defined inline because they are called every
// time a sample is calculated.
// ----------------------------------------------------------------------------

/* kinkiness of DAC:
 * some chips have more, some less. We should make this tunable. */
const float kinkiness = 0.966f;
const float sidcaps_6581 = 470e-12f;
const float outputleveldifference_lp_bp = 1.4f;
const float outputleveldifference_bp_hp = 1.2f;

RESID_INLINE
static float fastexp(float val) {
    typedef union {
        int i;
        float f;
    } conv;

    conv tmp;

    /* single precision fp has 1 + 8 + 23 bits, exponent bias is 127.
     * It therefore follows that we need to shift left by 23 bits, and to
     * calculate exp(x) instead of pow(2, x) we divide the power by ln(2). */
    const float a = (1 << 23) / M_LN2_f;
    /* The other factor corrects for the exponent bias so that 2^0 = 1. */
    const float b = (1 << 23) * 127;
    /* According to "A Fast, Compact Approximation of the Exponential Function"
     * by Nicol N. Schraudolph, 60801.48 yields the minimum RMS error for the
     * piecewise-linear approximation when using doubles (20 bits residual).
     * We have 23 bits, so we scale this value by 8. */
    const float c = 60801.48f * 8.f + 0.5f;

    /* Parenthesis are important: C standard disallows folding subtraction.
     * Unfortunately GCC appears to generate a write to memory rather than
     * handle this conversion entirely in registers. */
    tmp.i = (int)(a * val + (b - c));
    return tmp.f;
}

RESID_INLINE
float FilterFP::type3_w0(const float source, const float distoffset)
{
    /* The distortion appears to be the result of MOSFET entering saturation
     * mode. The conductance of a FET is proportional to:
     *
     * ohmic = 2 * (Vgs - Vt) * Vds - Vds^2
     * saturation = (Vgs - Vt)^2
     *
     * The FET switches to saturation mode when Vgs - Vt < Vds.
     *
     * In the circuit, the Vgs is mixed with the Vds signal, which gives
     * (Vgs + Vds) / 2 as the gate voltage. Doing the substitutions we get:
     *
     * ohmic = 2 * ((Vgs + Vds) / 2 - Vt) * Vds - Vds^2 = (Vgs - Vt) * Vds
     * saturation = ((Vgs + Vds) / 2 - Vt)^2
     *
     * Therefore: once the Vds crosses a threshold given by the gate and
     * threshold FET conductance begins to increase faster. The exact shape
     * for this effect is a parabola.
     *
     * The scaling term here tries to match the FC control level with
     * the signal level in simulation. On the chip, the FC control is
     * biased by forcing its highest DAC bit in the 1 position, thus
     * limiting the electrical range to half. Therefore one can guess that
     * the real FC range is half of the full voice range.
     *
     * On the simulation, FC goes to 2047 and the voices to 4095 * 255.
     * If the FC control was intact, then the scaling factor would be
     * 1/512. (Simulation voices are 512 times "louder" intrinsically.)
     * As the real chip's FC has reduced range, the scaling required to
     * match levels is 1/256. */

    float fetresistance = type3_fc_kink_exp;
    if (source > distoffset) {
        const float dist = source - distoffset;
        fetresistance *= fastexp(dist * type3_steepness * distortion_rate);
    }
    const float dynamic_resistance = type3_minimumfetresistance + fetresistance;

    /* 2 parallel resistors */
    const float _1_div_resistance = (type3_baseresistance + dynamic_resistance) / (type3_baseresistance * dynamic_resistance);
    /* 1.f / (clock * caps * resistance) */
    return distortion_CT * _1_div_resistance;
}

RESID_INLINE
float FilterFP::type4_w0()
{
    const float freq = type4_k * fc + type4_b;
    return 2.f * M_PI_f * freq / clock_frequency;
}

// ----------------------------------------------------------------------------
// SID clocking - 1 cycle.
// ----------------------------------------------------------------------------
RESID_INLINE
float FilterFP::clock(float voice1,
                   float voice2,
                   float voice3,
                   float ext_in)
{
    /* Avoid denormal numbers by using small offsets from 0 */
    float Vi = 0.f, Vnf = 0.f, Vf = 0.f;

    // Route voices into or around filter.
    ((filt & 1) ? Vi : Vnf) += voice1;
    ((filt & 2) ? Vi : Vnf) += voice2;
    // NB! Voice 3 is not silenced by voice3off if it is routed through
    // the filter.
    if (filt & 4)
        Vi += voice3;
    else if (! voice3off)
        Vnf += voice3;
    ((filt & 8) ? Vi : Vnf) += ext_in;

    if (! enabled)
        return (Vnf - Vi) * volf;

    if (hp_bp_lp & 1)
        Vf += Vlp;
    if (hp_bp_lp & 2)
        Vf += Vbp;
    if (hp_bp_lp & 4)
        Vf += Vhp;

    if (model == MOS6581FP) {
        float diff1, diff2;

        Vhp = Vbp * _1_div_Q * (1.f/outputleveldifference_bp_hp) - Vlp * (1.f/outputleveldifference_bp_hp) - Vi * 0.5f;

        /* the input summer mixing, or something like it... */
        diff1 = (Vlp - Vbp) * distortion_cf_threshold;
        diff2 = (Vhp - Vbp) * distortion_cf_threshold;
        Vlp -= diff1;
        Vbp += diff1;
        Vbp += diff2;
        Vhp -= diff2;

        /* Model output strip mixing. Doing it now that HP state
         * variable modifying still makes some difference.
         * (Phase error, though.) */
        if (hp_bp_lp & 1)
            Vlp += (Vf + Vnf - Vlp) * distortion_cf_threshold;
        if (hp_bp_lp & 2)
            Vbp += (Vf + Vnf - Vbp) * distortion_cf_threshold;
        if (hp_bp_lp & 4)
            Vhp += (Vf + Vnf - Vhp) * distortion_cf_threshold;

        /* Simulating the exponential VCR that the FET block is... */
        Vlp -= Vbp * type3_w0(Vbp, type3_fc_distortion_offset_bp);
        Vbp -= Vhp * type3_w0(Vhp, type3_fc_distortion_offset_hp) * outputleveldifference_bp_hp;

        /* Tuned based on Fred Gray's Break Thru. It is probably not a hard
         * discontinuity but a saturation effect... */
        if (Vnf > 3.2e6f)
            Vnf = 3.2e6f;

        Vf += Vnf + Vlp * (outputleveldifference_lp_bp - 1.f);
    } else {
        /* On the 8580, BP appears mixed in phase with the rest. */
        Vhp = -Vbp * _1_div_Q - Vlp - Vi;
        Vlp += Vbp * type4_w0_cache;
        Vbp += Vhp * type4_w0_cache;

        Vf += Vnf;
    }

    return Vf * volf;
}

RESID_INLINE
void FilterFP::nuke_denormals()
{
    /* We could use the flush-to-zero flag or denormals-are-zero on systems
     * where compiling with -msse and -mfpmath=sse is acceptable. Since this
     * doesn't include general VICE builds, we do this instead. */
    if (Vbp > -1e-12f && Vbp < 1e-12f)
        Vbp = 0;
    if (Vlp > -1e-12f && Vlp < 1e-12f)
        Vlp = 0;
}

#endif // not __FILTER_H__