Merge pull request #5030 from jriwanek-forks/ssi2001

Update resid-fp to 2.12.0 from libsidplayfp

src/sound/resid-fp/CMakeLists.txt

@@ -13,6 +13,8 @@
 #          Copyright 2020-2021 David Hrdlička.
 #
 
+set(CMAKE_CXX_STANDARD 17)
+
 add_library(resid-fp STATIC Dac.cpp EnvelopeGenerator.cpp ExternalFilter.cpp
     Filter.cpp Filter6581.cpp Filter8580.cpp FilterModelConfig.cpp
     FilterModelConfig6581.cpp FilterModelConfig8580.cpp

src/sound/resid-fp/Dac.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2016 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004,2010 Dag Lem <resid@nimrod.no>
  *
@@ -22,9 +22,14 @@
 
 #include "Dac.h"
 
+#include "sidcxx11.h"
+
 namespace reSIDfp
 {
 
+constexpr double MOSFET_LEAKAGE_6581 = 0.0075;
+constexpr double MOSFET_LEAKAGE_8580 = 0.0035;
+
 Dac::Dac(unsigned int bits) :
     dac(new double[bits]),
     dacLength(bits)
@@ -41,10 +46,8 @@ double Dac::getOutput(unsigned int input) const
 
     for (unsigned int i = 0; i < dacLength; i++)
     {
-        if ((input & (1 << i)) != 0)
+        const bool transistor_on = (input & (1 << i)) != 0;
-        {
+        dacValue += transistor_on ? dac[i] : dac[i] * leakage;
-            dacValue += dac[i];
-        }
     }
 
     return dacValue;
@@ -52,7 +55,7 @@ double Dac::getOutput(unsigned int input) const
 
 void Dac::kinkedDac(ChipModel chipModel)
 {
-    const double R_INFINITY = 1e6;
+    constexpr double R_INFINITY = 1e6;
 
     // Non-linearity parameter, 8580 DACs are perfectly linear
     const double _2R_div_R = chipModel == MOS6581 ? 2.20 : 2.00;
@@ -60,6 +63,10 @@ void Dac::kinkedDac(ChipModel chipModel)
     // 6581 DACs are not terminated by a 2R resistor
     const bool term = chipModel == MOS8580;
 
+    leakage = chipModel == MOS6581 ? MOSFET_LEAKAGE_6581 : MOSFET_LEAKAGE_8580;
+
+    double Vsum = 0.;
+
     // Calculate voltage contribution by each individual bit in the R-2R ladder.
     for (unsigned int set_bit = 0; set_bit < dacLength; set_bit++)
     {
@@ -102,18 +109,10 @@ void Dac::kinkedDac(ChipModel chipModel)
         }
 
         dac[set_bit] = Vn;
+        Vsum += Vn;
     }
 
     // Normalize to integerish behavior
-    double Vsum = 0.;
-
-    for (unsigned int i = 0; i < dacLength; i++)
-    {
-        Vsum += dac[i];
-    }
-
-    Vsum /= 1 << dacLength;
-
     for (unsigned int i = 0; i < dacLength; i++)
     {
         dac[i] /= Vsum;

src/sound/resid-fp/Dac.h

@@ -75,6 +75,15 @@ namespace reSIDfp
 class Dac
 {
 private:
+    /**
+     * DAC leakage
+     *
+     * "Even in standard transistors a small amount of current leaks even when they are technically switched off."
+     *
+     * https://en.wikipedia.org/wiki/Subthreshold_conduction
+     */
+    double leakage;
+
     /// analog values
     double * const dac;
 

src/sound/resid-fp/EnvelopeGenerator.cpp

@@ -79,7 +79,7 @@ void EnvelopeGenerator::reset()
     exponential_counter_period = 1;
     new_exponential_counter_period = 0;
 
-    state = RELEASE;
+    state = State::RELEASE;
     counter_enabled = true;
     rate = adsrtable[release];
 }
@@ -98,7 +98,7 @@ void EnvelopeGenerator::writeCONTROL_REG(unsigned char control)
         if (gate_next)
         {
             // Gate bit on:  Start attack, decay, sustain.
-            next_state = ATTACK;
+            next_state = State::ATTACK;
             state_pipeline = 2;
 
             if (resetLfsr || (exponential_pipeline == 2))
@@ -113,7 +113,7 @@ void EnvelopeGenerator::writeCONTROL_REG(unsigned char control)
         else
         {
             // Gate bit off: Start release.
-            next_state = RELEASE;
+            next_state = State::RELEASE;
             state_pipeline = envelope_pipeline > 0 ? 3 : 2;
         }
     }
@@ -124,11 +124,11 @@ void EnvelopeGenerator::writeATTACK_DECAY(unsigned char attack_decay)
     attack = (attack_decay >> 4) & 0x0f;
     decay = attack_decay & 0x0f;
 
-    if (state == ATTACK)
+    if (state == State::ATTACK)
     {
         rate = adsrtable[attack];
     }
-    else if (state == DECAY_SUSTAIN)
+    else if (state == State::DECAY_SUSTAIN)
     {
         rate = adsrtable[decay];
     }
@@ -146,7 +146,7 @@ void EnvelopeGenerator::writeSUSTAIN_RELEASE(unsigned char sustain_release)
 
     release = sustain_release & 0x0f;
 
-    if (state == RELEASE)
+    if (state == State::RELEASE)
     {
         rate = adsrtable[release];
     }

src/sound/resid-fp/EnvelopeGenerator.h

@@ -47,68 +47,68 @@ private:
      * The envelope state machine's distinct states. In addition to this,
      * envelope has a hold mode, which freezes envelope counter to zero.
      */
-    enum State
+    enum class State
     {
         ATTACK, DECAY_SUSTAIN, RELEASE
     };
 
 private:
     /// XOR shift register for ADSR prescaling.
-    unsigned int lfsr;
+    unsigned int lfsr = 0x7fff;
 
     /// Comparison value (period) of the rate counter before next event.
-    unsigned int rate;
+    unsigned int rate = 0;
 
     /**
      * During release mode, the SID approximates envelope decay via piecewise
      * linear decay rate.
      */
-    unsigned int exponential_counter;
+    unsigned int exponential_counter = 0;
 
     /**
      * Comparison value (period) of the exponential decay counter before next
      * decrement.
      */
-    unsigned int exponential_counter_period;
+    unsigned int exponential_counter_period = 1;
-    unsigned int new_exponential_counter_period;
+    unsigned int new_exponential_counter_period = 0;
 
-    unsigned int state_pipeline;
+    unsigned int state_pipeline = 0;
 
     ///
-    unsigned int envelope_pipeline;
+    unsigned int envelope_pipeline = 0;
 
-    unsigned int exponential_pipeline;
+    unsigned int exponential_pipeline = 0;
 
     /// Current envelope state
-    State state;
+    State state = State::RELEASE;
-    State next_state;
+    State next_state = State::RELEASE;
 
     /// Whether counter is enabled. Only switching to ATTACK can release envelope.
-    bool counter_enabled;
+    bool counter_enabled = true;
 
     /// Gate bit
-    bool gate;
+    bool gate = false;
 
     ///
-    bool resetLfsr;
+    bool resetLfsr = false;
 
     /// The current digital value of envelope output.
-    unsigned char envelope_counter;
+    unsigned char envelope_counter = 0xaa;
 
     /// Attack register
-    unsigned char attack;
+    unsigned char attack = 0;
 
     /// Decay register
-    unsigned char decay;
+    unsigned char decay = 0;
 
     /// Sustain register
-    unsigned char sustain;
+    unsigned char sustain = 0;
 
     /// Release register
-    unsigned char release;
+    unsigned char release = 0;
 
     /// The ENV3 value, sampled at the first phase of the clock
-    unsigned char env3;
+    unsigned char env3 = 0;
 
 private:
     static const unsigned int adsrtable[16];
@@ -129,31 +129,6 @@ public:
      */
     unsigned int output() const { return envelope_counter; }
 
-    /**
-     * Constructor.
-     */
-    EnvelopeGenerator() :
-        lfsr(0x7fff),
-        rate(0),
-        exponential_counter(0),
-        exponential_counter_period(1),
-        new_exponential_counter_period(0),
-        state_pipeline(0),
-        envelope_pipeline(0),
-        exponential_pipeline(0),
-        state(RELEASE),
-        next_state(RELEASE),
-        counter_enabled(true),
-        gate(false),
-        resetLfsr(false),
-        envelope_counter(0xaa),
-        attack(0),
-        decay(0),
-        sustain(0),
-        release(0),
-        env3(0)
-    {}
-
     /**
      * SID reset.
      */
@@ -218,15 +193,15 @@ void EnvelopeGenerator::clock()
     {
         if (likely(counter_enabled))
         {
-            if (state == ATTACK)
+            if (state == State::ATTACK)
             {
                 if (++envelope_counter==0xff)
                 {
-                    next_state = DECAY_SUSTAIN;
+                    next_state = State::DECAY_SUSTAIN;
                     state_pipeline = 3;
                 }
             }
-            else if ((state == DECAY_SUSTAIN) || (state == RELEASE))
+            else if ((state == State::DECAY_SUSTAIN) || (state == State::RELEASE))
             {
                 if (--envelope_counter==0x00)
                 {
@@ -241,8 +216,8 @@ void EnvelopeGenerator::clock()
     {
         exponential_counter = 0;
 
-        if (((state == DECAY_SUSTAIN) && (envelope_counter != sustain))
+        if (((state == State::DECAY_SUSTAIN) && (envelope_counter != sustain))
-            || (state == RELEASE))
+            || (state == State::RELEASE))
         {
             // The envelope counter can flip from 0x00 to 0xff by changing state to
             // attack, then to release. The envelope counter will then continue
@@ -257,7 +232,7 @@ void EnvelopeGenerator::clock()
         lfsr = 0x7fff;
         resetLfsr = false;
 
-        if (state == ATTACK)
+        if (state == State::ATTACK)
         {
             // The first envelope step in the attack state also resets the exponential
             // counter. This has been verified by sampling ENV3.
@@ -344,7 +319,7 @@ void EnvelopeGenerator::state_change()
 
     switch (next_state)
     {
-    case ATTACK:
+    case State::ATTACK:
         if (state_pipeline == 1)
         {
             // The decay rate is "accidentally" enabled during first cycle of attack phase
@@ -352,24 +327,24 @@ void EnvelopeGenerator::state_change()
         }
         else if (state_pipeline == 0)
         {
-            state = ATTACK;
+            state = State::ATTACK;
             // The attack rate is correctly enabled during second cycle of attack phase
             rate = adsrtable[attack];
             counter_enabled = true;
         }
         break;
-    case DECAY_SUSTAIN:
+    case State::DECAY_SUSTAIN:
         if (state_pipeline == 0)
         {
-            state = DECAY_SUSTAIN;
+            state = State::DECAY_SUSTAIN;
             rate = adsrtable[decay];
         }
         break;
-    case RELEASE:
+    case State::RELEASE:
-        if (((state == ATTACK) && (state_pipeline == 0))
+        if (((state == State::ATTACK) && (state_pipeline == 0))
-            || ((state == DECAY_SUSTAIN) && (state_pipeline == 1)))
+            || ((state == State::DECAY_SUSTAIN) && (state_pipeline == 1)))
         {
-            state = RELEASE;
+            state = State::RELEASE;
             rate = adsrtable[release];
         }
         break;

src/sound/resid-fp/ExternalFilter.cpp

@@ -38,9 +38,7 @@ inline double getRC(double res, double cap)
     return res * cap;
 }
 
-ExternalFilter::ExternalFilter() :
+ExternalFilter::ExternalFilter()
-    w0lp_1_s7(0),
-    w0hp_1_s17(0)
 {
     reset();
 }

src/sound/resid-fp/ExternalFilter.h

@@ -34,8 +34,6 @@ namespace reSIDfp
  * acts as a high-pass filter with a cutoff dependent on the attached audio
  * equipment impedance. Here we suppose an impedance of 10kOhm resulting
  * in a 3 dB attenuation at 1.6Hz.
- * To operate properly the 6581 audio output needs a pull-down resistor
- *(1KOhm recommended, not needed on 8580)
  *
  * ~~~
  *                                 9/12V
@@ -47,15 +45,18 @@ namespace reSIDfp
  *          |        |  pF    +-C----o-----C-----+ 10k
  *                             470   |           |
  *         GND      GND         pF   R 1K        | amp
- *          *                   *    |           +-----
+ *          *                   **   |           +-----
  *
  *                                  GND
  * ~~~
  *
  * The STC networks are connected with a [BJT] based [common collector]
  * used as a voltage follower (featuring a 2SC1815 NPN transistor).
- * * The C64c board additionally includes a [bootstrap] condenser to increase
+ *
- * the input impedance of the common collector.
+ * * To operate properly the 6581 audio output needs a pull-down resistor
+ *   (1KOhm recommended, not needed on 8580)
+ * ** The C64c board additionally includes a [bootstrap] condenser to increase
+ *    the input impedance of the common collector.
  *
  * [BJT]: https://en.wikipedia.org/wiki/Bipolar_junction_transistor
  * [common collector]: https://en.wikipedia.org/wiki/Common_collector
@@ -70,9 +71,9 @@ private:
     /// Highpass filter voltage
     int Vhp;
 
-    int w0lp_1_s7;
+    int w0lp_1_s7 = 0;
 
-    int w0hp_1_s17;
+    int w0hp_1_s17 = 0;
 
 public:
     /**
@@ -80,7 +81,7 @@ public:
      *
      * @param input
      */
-    int clock(unsigned short input);
+    int clock(int input);
 
     /**
      * Constructor.
@@ -108,9 +109,9 @@ namespace reSIDfp
 {
 
 RESID_INLINE
-int ExternalFilter::clock(unsigned short input)
+int ExternalFilter::clock(int input)
 {
-    const int Vi = (static_cast<unsigned int>(input)<<11) - (1 << (11+15));
+    const int Vi = (input<<11) - (1 << (11+15));
     const int dVlp = (w0lp_1_s7 * (Vi - Vlp) >> 7);
     const int dVhp = (w0hp_1_s17 * (Vlp - Vhp) >> 17);
     Vlp += dVlp;

src/sound/resid-fp/Filter.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2013 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004 Dag Lem <resid@nimrod.no>
  *
@@ -20,11 +20,87 @@
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  */
 
+#define FILTER_CPP
+
 #include "Filter.h"
 
 namespace reSIDfp
 {
 
+void Filter::updateMixing()
+{
+    currentVolume = volume[vol];
+
+    unsigned int Nsum = 0;
+    unsigned int Nmix = 0;
+
+    (filt1 ? Nsum : Nmix)++;
+    (filt2 ? Nsum : Nmix)++;
+
+    if (filt3) Nsum++;
+    else if (!voice3off) Nmix++;
+
+    (filtE ? Nsum : Nmix)++;
+
+    currentSummer = summer[Nsum];
+
+    if (lp) Nmix++;
+    if (bp) Nmix++;
+    if (hp) Nmix++;
+
+    currentMixer = mixer[Nmix];
+}
+
+void Filter::writeFC_LO(unsigned char fc_lo)
+{
+    fc = (fc & 0x7f8) | (fc_lo & 0x007);
+    updateCenterFrequency();
+}
+
+void Filter::writeFC_HI(unsigned char fc_hi)
+{
+    fc = (fc_hi << 3 & 0x7f8) | (fc & 0x007);
+    updateCenterFrequency();
+}
+
+void Filter::writeRES_FILT(unsigned char res_filt)
+{
+    filt = res_filt;
+
+    updateResonance((res_filt >> 4) & 0x0f);
+
+    if (enabled)
+    {
+        filt1 = (filt & 0x01) != 0;
+        filt2 = (filt & 0x02) != 0;
+        filt3 = (filt & 0x04) != 0;
+        filtE = (filt & 0x08) != 0;
+    }
+
+    updateMixing();
+}
+
+void Filter::writeMODE_VOL(unsigned char mode_vol)
+{
+    vol = mode_vol & 0x0f;
+    lp = (mode_vol & 0x10) != 0;
+    bp = (mode_vol & 0x20) != 0;
+    hp = (mode_vol & 0x40) != 0;
+    voice3off = (mode_vol & 0x80) != 0;
+
+    updateMixing();
+}
+
+Filter::Filter(FilterModelConfig& fmc) :
+    mixer(fmc.getMixer()),
+    summer(fmc.getSummer()),
+    resonance(fmc.getResonance()),
+    volume(fmc.getVolume()),
+    fmc(fmc)
+{
+    input(0);
+}
+
 void Filter::enable(bool enable)
 {
     enabled = enable;
@@ -47,44 +123,4 @@ void Filter::reset()
     writeRES_FILT(0);
 }
 
-void Filter::writeFC_LO(unsigned char fc_lo)
-{
-    fc = (fc & 0x7f8) | (fc_lo & 0x007);
-    updatedCenterFrequency();
-}
-
-void Filter::writeFC_HI(unsigned char fc_hi)
-{
-    fc = (fc_hi << 3 & 0x7f8) | (fc & 0x007);
-    updatedCenterFrequency();
-}
-
-void Filter::writeRES_FILT(unsigned char res_filt)
-{
-    filt = res_filt;
-
-    updateResonance((res_filt >> 4) & 0x0f);
-
-    if (enabled)
-    {
-        filt1 = (filt & 0x01) != 0;
-        filt2 = (filt & 0x02) != 0;
-        filt3 = (filt & 0x04) != 0;
-        filtE = (filt & 0x08) != 0;
-    }
-
-    updatedMixing();
-}
-
-void Filter::writeMODE_VOL(unsigned char mode_vol)
-{
-    vol = mode_vol & 0x0f;
-    lp = (mode_vol & 0x10) != 0;
-    bp = (mode_vol & 0x20) != 0;
-    hp = (mode_vol & 0x40) != 0;
-    voice3off = (mode_vol & 0x80) != 0;
-
-    updatedMixing();
-}
-
 } // namespace reSIDfp

src/sound/resid-fp/Filter.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2017 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004 Dag Lem <resid@nimrod.no>
  *
@@ -23,6 +23,10 @@
 #ifndef FILTER_H
 #define FILTER_H
 
+#include "FilterModelConfig.h"
+
+#include "siddefs-fp.h"
+
 namespace reSIDfp
 {
 
@@ -31,93 +35,97 @@ namespace reSIDfp
  */
 class Filter
 {
+private:
+    unsigned short** mixer;
+    unsigned short** summer;
+    unsigned short** resonance;
+    unsigned short** volume;
+
 protected:
-    /// Current volume amplifier setting.
+    FilterModelConfig& fmc;
-    unsigned short* currentGain;
 
     /// Current filter/voice mixer setting.
-    unsigned short* currentMixer;
+    unsigned short* currentMixer = nullptr;
 
     /// Filter input summer setting.
-    unsigned short* currentSummer;
+    unsigned short* currentSummer = nullptr;
 
     /// Filter resonance value.
-    unsigned short* currentResonance;
+    unsigned short* currentResonance = nullptr;
+
+    /// Current volume amplifier setting.
+    unsigned short* currentVolume = nullptr;
 
     /// Filter highpass state.
-    int Vhp;
+    int Vhp = 0;
 
     /// Filter bandpass state.
-    int Vbp;
+    int Vbp = 0;
 
     /// Filter lowpass state.
-    int Vlp;
+    int Vlp = 0;
 
     /// Filter external input.
-    int ve;
+    int Ve = 0;
 
     /// Filter cutoff frequency.
-    unsigned int fc;
+    unsigned int fc = 0;
 
     /// Routing to filter or outside filter
-    bool filt1, filt2, filt3, filtE;
+    //@{
+    bool filt1 = false;
+    bool filt2 = false;
+    bool filt3 = false;
+    bool filtE = false;
+    //@}
 
     /// Switch voice 3 off.
-    bool voice3off;
+    bool voice3off = false;
 
     /// Highpass, bandpass, and lowpass filter modes.
-    bool hp, bp, lp;
+    //@{
-
+    bool hp = false;
-    /// Current volume.
+    bool bp = false;
-    unsigned char vol;
+    bool lp = false;
+    //@}
 
 private:
+    /// Current volume.
+    unsigned char vol = 0;
+
     /// Filter enabled.
-    bool enabled;
+    bool enabled = true;
 
     /// Selects which inputs to route through filter.
-    unsigned char filt;
+    unsigned char filt = 0;
 
 protected:
     /**
-     * Set filter cutoff frequency.
+     * Update filter cutoff frequency.
      */
-    virtual void updatedCenterFrequency() = 0;
+    virtual void updateCenterFrequency() = 0;
 
     /**
-     * Set filter resonance.
+     * Update filter resonance.
+     *
+     * @param res the new resonance value
      */
-    virtual void updateResonance(unsigned char res) = 0;
+    void updateResonance(unsigned char res) { currentResonance = resonance[res]; }
 
     /**
      * Mixing configuration modified (offsets change)
      */
-    virtual void updatedMixing() = 0;
+    void updateMixing();
+
+    /**
+     * Get the filter cutoff register value
+     */
+    unsigned int getFC() const { return fc; }
 
 public:
-    Filter() :
+    Filter(FilterModelConfig& fmc);
-        currentGain(nullptr),
-        currentMixer(nullptr),
-        currentSummer(nullptr),
-        currentResonance(nullptr),
-        Vhp(0),
-        Vbp(0),
-        Vlp(0),
-        ve(0),
-        fc(0),
-        filt1(false),
-        filt2(false),
-        filt3(false),
-        filtE(false),
-        voice3off(false),
-        hp(false),
-        bp(false),
-        lp(false),
-        vol(0),
-        enabled(true),
-        filt(0) {}
 
-    virtual ~Filter() {}
+    virtual ~Filter() = default;
 
     /**
      * SID clocking - 1 cycle
@@ -169,7 +177,14 @@ public:
      */
     void writeMODE_VOL(unsigned char mode_vol);
 
-    virtual void input(int input) = 0;
+    /**
+     * Apply a signal to EXT-IN
+     *
+     * @param input a signed 16 bit sample
+     */
+    void input(short input) { Ve = fmc.getNormalizedVoice(input/32768.f, 0); }
+
+    inline int getNormalizedVoice(float value, unsigned int env) const { return fmc.getNormalizedVoice(value, env); }
 };
 
 } // namespace reSIDfp

src/sound/resid-fp/Filter6581.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2015 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004,2010 Dag Lem <resid@nimrod.no>
  *
@@ -20,8 +20,6 @@
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  */
 
-#define FILTER6581_CPP
-
 #include "Filter6581.h"
 
 #include "Integrator6581.h"
@@ -29,47 +27,60 @@
 namespace reSIDfp
 {
 
+unsigned short Filter6581::clock(int voice1, int voice2, int voice3)
+{
+    const int V1 = voice1;
+    const int V2 = voice2;
+    // Voice 3 is silenced by voice3off if it is not routed through the filter.
+    const int V3 = (filt3 || !voice3off) ? voice3 : 0;
+
+    int Vsum = 0;
+    int Vmix = 0;
+
+    (filt1 ? Vsum : Vmix) += V1;
+    (filt2 ? Vsum : Vmix) += V2;
+    (filt3 ? Vsum : Vmix) += V3;
+    (filtE ? Vsum : Vmix) += Ve;
+
+    Vhp = currentSummer[currentResonance[Vbp] + Vlp + Vsum];
+    Vbp = hpIntegrator.solve(Vhp);
+    Vlp = bpIntegrator.solve(Vbp);
+
+    int Vfilt = 0;
+    if (lp) Vfilt += Vlp;
+    if (bp) Vfilt += Vbp;
+    if (hp) Vfilt += Vhp;
+
+    // The filter input resistors are slightly bigger than the voice ones
+    // Scale the values accordingly
+    constexpr int filterGain = static_cast<int>(0.93 * (1 << 12));
+    Vfilt = (Vfilt * filterGain) >> 12;
+
+    return currentVolume[currentMixer[Vmix + Vfilt]];
+}
+
 Filter6581::~Filter6581()
 {
     delete [] f0_dac;
 }
 
-void Filter6581::updatedCenterFrequency()
+void Filter6581::updateCenterFrequency()
 {
-    const unsigned short Vw = f0_dac[fc];
+    const unsigned short Vw = f0_dac[getFC()];
-    hpIntegrator->setVw(Vw);
+    hpIntegrator.setVw(Vw);
-    bpIntegrator->setVw(Vw);
+    bpIntegrator.setVw(Vw);
-}
-
-void Filter6581::updatedMixing()
-{
-    currentGain = gain_vol[vol];
-
-    unsigned int ni = 0;
-    unsigned int no = 0;
-
-    (filt1 ? ni : no)++;
-    (filt2 ? ni : no)++;
-
-    if (filt3) ni++;
-    else if (!voice3off) no++;
-
-    (filtE ? ni : no)++;
-
-    currentSummer = summer[ni];
-
-    if (lp) no++;
-    if (bp) no++;
-    if (hp) no++;
-
-    currentMixer = mixer[no];
 }
 
 void Filter6581::setFilterCurve(double curvePosition)
 {
     delete [] f0_dac;
     f0_dac = FilterModelConfig6581::getInstance()->getDAC(curvePosition);
-    updatedCenterFrequency();
+    updateCenterFrequency();
+}
+
+void Filter6581::setFilterRange(double adjustment)
+{
+    FilterModelConfig6581::getInstance()->setFilterRange(adjustment);
 }
 
 } // namespace reSIDfp

src/sound/resid-fp/Filter6581.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2022 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004,2010 Dag Lem <resid@nimrod.no>
  *
@@ -23,12 +23,9 @@
 #ifndef FILTER6581_H
 #define FILTER6581_H
 
-#include "siddefs-fp.h"
-
-#include <memory>
-
 #include "Filter.h"
 #include "FilterModelConfig6581.h"
+#include "Integrator6581.h"
 
 #include "sidcxx11.h"
 
@@ -108,7 +105,7 @@ class Integrator6581;
  *     |   |   |   v1       |    |   |                               |
  * D0  |   |   |   \ ---R8--+    |   |   +---------------------------+
  *     |   |   |   |             |   |   |
- *     R6  R6  R6  R6            R6  R6  R6
+ *     R6  R6  R6  R6            R6* R6* R6*
  *     |   |   |   | $18         |   |   |  $18
  *     |    \  |   | D7: 1=open   \   \   \ D6 - D4: 0=open
  *     |   |   |   |             |   |   |
@@ -143,6 +140,7 @@ class Integrator6581;
  *
  *     R2  ~  2.0*R1
  *     R6  ~  6.0*R1
+ *     R6* ~  1.07*R6
  *     R8  ~  8.0*R1
  *     R24 ~ 24.0*R1
  *
@@ -322,104 +320,49 @@ class Integrator6581;
 class Filter6581 final : public Filter
 {
 private:
-    const unsigned short* f0_dac;
-
-    unsigned short** mixer;
-    unsigned short** summer;
-    unsigned short** gain_res;
-    unsigned short** gain_vol;
-
-    const int voiceScaleS11;
-    const int voiceDC;
-
     /// VCR + associated capacitor connected to highpass output.
-    std::unique_ptr<Integrator6581> const hpIntegrator;
+    Integrator6581 hpIntegrator;
 
     /// VCR + associated capacitor connected to bandpass output.
-    std::unique_ptr<Integrator6581> const bpIntegrator;
+    Integrator6581 bpIntegrator;
+
+    const unsigned short* f0_dac;
 
 protected:
     /**
      * Set filter cutoff frequency.
      */
-    void updatedCenterFrequency() override;
+    void updateCenterFrequency() override;
-
-    /**
-     * Set filter resonance.
-     *
-     * In the MOS 6581, 1/Q is controlled linearly by res.
-     */
-    void updateResonance(unsigned char res) override { currentResonance = gain_res[res]; }
-
-    void updatedMixing() override;
 
 public:
     Filter6581() :
-        f0_dac(FilterModelConfig6581::getInstance()->getDAC(0.5)),
+        Filter(*FilterModelConfig6581::getInstance()),
-        mixer(FilterModelConfig6581::getInstance()->getMixer()),
+        hpIntegrator(*FilterModelConfig6581::getInstance()),
-        summer(FilterModelConfig6581::getInstance()->getSummer()),
+        bpIntegrator(*FilterModelConfig6581::getInstance()),
-        gain_res(FilterModelConfig6581::getInstance()->getGainRes()),
+        f0_dac(FilterModelConfig6581::getInstance()->getDAC(0.5))
-        gain_vol(FilterModelConfig6581::getInstance()->getGainVol()),
+    {}
-        voiceScaleS11(FilterModelConfig6581::getInstance()->getVoiceScaleS11()),
-        voiceDC(FilterModelConfig6581::getInstance()->getNormalizedVoiceDC()),
-        hpIntegrator(FilterModelConfig6581::getInstance()->buildIntegrator()),
-        bpIntegrator(FilterModelConfig6581::getInstance()->buildIntegrator())
-    {
-        input(0);
-    }
 
-    ~Filter6581();
+    ~Filter6581() override;
 
-    unsigned short clock(int voice1, int voice2, int voice3) override;
+    unsigned short clock(int v1, int v2, int v3) override;
-
-    void input(int sample) override { ve = (sample * voiceScaleS11 * 3 >> 11) + mixer[0][0]; }
 
     /**
      * Set filter curve type based on single parameter.
      *
-     * @param curvePosition 0 .. 1, where 0 sets center frequency high ("light") and 1 sets it low ("dark"), default is 0.5
+     * @param curvePosition 0 .. 1, where 0 sets center frequency high ("bright") and 1 sets it low ("dark").
+     *                      Default is 0.5
      */
     void setFilterCurve(double curvePosition);
+
+    /**
+     * Set filter offset and range based on single parameter.
+     *
+     * @param adjustment 0 .. 1, where 0 sets center frequency low ("dark"), 1 sets it high ("bright").
+     *                   This also affects the range. Default is 0.5
+     */
+    void setFilterRange(double adjustment);
 };
 
 } // namespace reSIDfp
 
-#if RESID_INLINING || defined(FILTER6581_CPP)
-
-#include "Integrator6581.h"
-
-namespace reSIDfp
-{
-
-RESID_INLINE
-unsigned short Filter6581::clock(int voice1, int voice2, int voice3)
-{
-    voice1 = (voice1 * voiceScaleS11 >> 15) + voiceDC;
-    voice2 = (voice2 * voiceScaleS11 >> 15) + voiceDC;
-    // Voice 3 is silenced by voice3off if it is not routed through the filter.
-    voice3 = (filt3 || !voice3off) ? (voice3 * voiceScaleS11 >> 15) + voiceDC : 0;
-
-    int Vi = 0;
-    int Vo = 0;
-
-    (filt1 ? Vi : Vo) += voice1;
-    (filt2 ? Vi : Vo) += voice2;
-    (filt3 ? Vi : Vo) += voice3;
-    (filtE ? Vi : Vo) += ve;
-
-    Vhp = currentSummer[currentResonance[Vbp] + Vlp + Vi];
-    Vbp = hpIntegrator->solve(Vhp);
-    Vlp = bpIntegrator->solve(Vbp);
-
-    if (lp) Vo += Vlp;
-    if (bp) Vo += Vbp;
-    if (hp) Vo += Vhp;
-
-    return currentGain[currentMixer[Vo]];
-}
-
-} // namespace reSIDfp
-
-#endif
-
 #endif

src/sound/resid-fp/Filter8580.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2019 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004,2010 Dag Lem <resid@nimrod.no>
  *
@@ -20,8 +20,6 @@
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  */
 
-#define FILTER8580_CPP
-
 #include "Filter8580.h"
 
 #include "Integrator8580.h"
@@ -29,6 +27,32 @@
 namespace reSIDfp
 {
 
+unsigned short Filter8580::clock(int voice1, int voice2, int voice3)
+{
+    const int V1 = voice1;
+    const int V2 = voice2;
+    // Voice 3 is silenced by voice3off if it is not routed through the filter.
+    const int V3 = (filt3 || !voice3off) ? voice3 : 0;
+
+    int Vsum = 0;
+    int Vmix = 0;
+
+    (filt1 ? Vsum : Vmix) += V1;
+    (filt2 ? Vsum : Vmix) += V2;
+    (filt3 ? Vsum : Vmix) += V3;
+    (filtE ? Vsum : Vmix) += Ve;
+
+    Vhp = currentSummer[currentResonance[Vbp] + Vlp + Vsum];
+    Vbp = hpIntegrator.solve(Vhp);
+    Vlp = bpIntegrator.solve(Vbp);
+
+    if (lp) Vmix += Vlp;
+    if (bp) Vmix += Vbp;
+    if (hp) Vmix += Vhp;
+
+    return currentVolume[currentMixer[Vmix]];
+}
+
 /**
  * W/L ratio of frequency DAC bit 0,
  * other bit are proportional.
@@ -37,18 +61,18 @@ namespace reSIDfp
  */
 const double DAC_WL0 = 0.00615;
 
-Filter8580::~Filter8580() {}
+Filter8580::~Filter8580() = default;
 
-void Filter8580::updatedCenterFrequency()
+void Filter8580::updateCenterFrequency()
 {
     double wl;
     double dacWL = DAC_WL0;
-    if (fc)
+    if (getFC())
     {
         wl = 0.;
         for (unsigned int i = 0; i < 11; i++)
         {
-            if (fc & (1 << i))
+            if (getFC() & (1 << i))
             {
                 wl += dacWL;
             }
@@ -60,32 +84,8 @@ void Filter8580::updatedCenterFrequency()
         wl = dacWL/2.;
     }
 
-    hpIntegrator->setFc(wl);
+    hpIntegrator.setFc(wl);
-    bpIntegrator->setFc(wl);
+    bpIntegrator.setFc(wl);
-}
-
-void Filter8580::updatedMixing()
-{
-    currentGain = gain_vol[vol];
-
-    unsigned int ni = 0;
-    unsigned int no = 0;
-
-    (filt1 ? ni : no)++;
-    (filt2 ? ni : no)++;
-
-    if (filt3) ni++;
-    else if (!voice3off) no++;
-
-    (filtE ? ni : no)++;
-
-    currentSummer = summer[ni];
-
-    if (lp) no++;
-    if (bp) no++;
-    if (hp) no++;
-
-    currentMixer = mixer[no];
 }
 
 void Filter8580::setFilterCurve(double curvePosition)
@@ -94,8 +94,8 @@ void Filter8580::setFilterCurve(double curvePosition)
     // 1.2 <= cp <= 1.8
     cp = 1.8 - curvePosition * 3./5.;
 
-    hpIntegrator->setV(cp);
+    hpIntegrator.setV(cp);
-    bpIntegrator->setV(cp);
+    bpIntegrator.setV(cp);
 }
 
 } // namespace reSIDfp

src/sound/resid-fp/Filter8580.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2022 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004,2010 Dag Lem <resid@nimrod.no>
  *
@@ -23,10 +23,6 @@
 #ifndef FILTER8580_H
 #define FILTER8580_H
 
-#include "siddefs-fp.h"
-
-#include <memory>
-
 #include "Filter.h"
 #include "FilterModelConfig8580.h"
 #include "Integrator8580.h"
@@ -281,58 +277,32 @@ class Integrator8580;
 class Filter8580 final : public Filter
 {
 private:
-    unsigned short** mixer;
-    unsigned short** summer;
-    unsigned short** gain_res;
-    unsigned short** gain_vol;
-
-    const int voiceScaleS11;
-    const int voiceDC;
-
-    double cp;
-
     /// VCR + associated capacitor connected to highpass output.
-    std::unique_ptr<Integrator8580> const hpIntegrator;
+    Integrator8580 hpIntegrator;
 
     /// VCR + associated capacitor connected to bandpass output.
-    std::unique_ptr<Integrator8580> const bpIntegrator;
+    Integrator8580 bpIntegrator;
+
+    double cp;
 
 protected:
     /**
      * Set filter cutoff frequency.
      */
-    void updatedCenterFrequency() override;
+    void updateCenterFrequency() override;
-
-    /**
-     * Set filter resonance.
-	 *
-     * @param res the new resonance value
-     */
-    void updateResonance(unsigned char res) override { currentResonance = gain_res[res]; }
-
-    void updatedMixing() override;
 
 public:
     Filter8580() :
-        mixer(FilterModelConfig8580::getInstance()->getMixer()),
+        Filter(*FilterModelConfig8580::getInstance()),
-        summer(FilterModelConfig8580::getInstance()->getSummer()),
+        hpIntegrator(*FilterModelConfig8580::getInstance()),
-        gain_res(FilterModelConfig8580::getInstance()->getGainRes()),
+        bpIntegrator(*FilterModelConfig8580::getInstance())
-        gain_vol(FilterModelConfig8580::getInstance()->getGainVol()),
-        voiceScaleS11(FilterModelConfig8580::getInstance()->getVoiceScaleS11()),
-        voiceDC(FilterModelConfig8580::getInstance()->getNormalizedVoiceDC()),
-        cp(0.5),
-        hpIntegrator(FilterModelConfig8580::getInstance()->buildIntegrator()),
-        bpIntegrator(FilterModelConfig8580::getInstance()->buildIntegrator())
     {
-        setFilterCurve(cp);
+        setFilterCurve(0.5);
-        input(0);
     }
 
-    ~Filter8580();
+    ~Filter8580() override;
 
-    unsigned short clock(int voice1, int voice2, int voice3) override;
+    unsigned short clock(int v1, int v2, int v3) override;
-
-    void input(int sample) override { ve = (sample * voiceScaleS11 * 3 >> 11) + mixer[0][0]; }
 
     /**
      * Set filter curve type based on single parameter.
@@ -344,40 +314,4 @@ public:
 
 } // namespace reSIDfp
 
-#if RESID_INLINING || defined(FILTER8580_CPP)
-
-namespace reSIDfp
-{
-
-RESID_INLINE
-unsigned short Filter8580::clock(int voice1, int voice2, int voice3)
-{
-    voice1 = (voice1 * voiceScaleS11 >> 15) + voiceDC;
-    voice2 = (voice2 * voiceScaleS11 >> 15) + voiceDC;
-    // Voice 3 is silenced by voice3off if it is not routed through the filter.
-    voice3 = (filt3 || !voice3off) ? (voice3 * voiceScaleS11 >> 15) + voiceDC : 0;
-
-    int Vi = 0;
-    int Vo = 0;
-
-    (filt1 ? Vi : Vo) += voice1;
-    (filt2 ? Vi : Vo) += voice2;
-    (filt3 ? Vi : Vo) += voice3;
-    (filtE ? Vi : Vo) += ve;
-
-    Vhp = currentSummer[currentResonance[Vbp] + Vlp + Vi];
-    Vbp = hpIntegrator->solve(Vhp);
-    Vlp = bpIntegrator->solve(Vbp);
-
-    if (lp) Vo += Vlp;
-    if (bp) Vo += Vbp;
-    if (hp) Vo += Vhp;
-
-    return currentGain[currentMixer[Vo]];
-}
-
-} // namespace reSIDfp
-
-#endif
-
 #endif

src/sound/resid-fp/FilterModelConfig.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2022 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004,2010 Dag Lem
  *
@@ -29,7 +29,6 @@ namespace reSIDfp
 
 FilterModelConfig::FilterModelConfig(
     double vvr,
-    double vdv,
     double c,
     double vdd,
     double vth,
@@ -37,21 +36,19 @@ FilterModelConfig::FilterModelConfig(
     const Spline::Point *opamp_voltage,
     int opamp_size
 ) :
-    voice_voltage_range(vvr),
-    voice_DC_voltage(vdv),
     C(c),
     Vdd(vdd),
     Vth(vth),
-    Ut(26.0e-3),
-    uCox(ucox),
     Vddt(Vdd - Vth),
     vmin(opamp_voltage[0].x),
     vmax(std::max(Vddt, opamp_voltage[0].y)),
     denorm(vmax - vmin),
     norm(1.0 / denorm),
     N16(norm * ((1 << 16) - 1)),
-    currFactorCoeff(denorm * (uCox / 2. * 1.0e-6 / C))
+    voice_voltage_range(vvr)
 {
+    setUCox(ucox);
+
     // Convert op-amp voltage transfer to 16 bit values.
 
     std::vector<Spline::Point> scaled_voltage(opamp_size);
@@ -79,4 +76,29 @@ FilterModelConfig::FilterModelConfig(
     }
 }
 
+FilterModelConfig::~FilterModelConfig()
+{
+    for (int i = 0; i < 8; i++)
+    {
+        delete [] mixer[i];
+    }
+
+    for (int i = 0; i < 5; i++)
+    {
+        delete [] summer[i];
+    }
+
+    for (int i = 0; i < 16; i++)
+    {
+        delete [] volume[i];
+        delete [] resonance[i];
+    }
+}
+
+void FilterModelConfig::setUCox(double new_uCox)
+{
+    uCox = new_uCox;
+    currFactorCoeff = denorm * (uCox / 2. * 1.0e-6 / C);
+}
+
 } // namespace reSIDfp

src/sound/resid-fp/FilterModelConfig.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2023 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004,2010 Dag Lem
  *
@@ -24,8 +24,10 @@
 #define FILTERMODELCONFIG_H
 
 #include <algorithm>
+#include <random>
 #include <cassert>
 
+#include "OpAmp.h"
 #include "Spline.h"
 
 #include "sidcxx11.h"
@@ -35,20 +37,46 @@ namespace reSIDfp
 
 class FilterModelConfig
 {
-protected:
+private:
-    const double voice_voltage_range;
+    /*
-    const double voice_DC_voltage;
+     * Hack to add quick dither when converting values from float to int
+     * and avoid quantization noise.
+     * Hopefully this can be removed the day we move all the analog part
+     * processing to floats.
+     *
+     * Not sure about the effect of using such small buffer of numbers
+     * since the random sequence repeats every 1024 values but for
+     * now it seems to do the job.
+     */
+    class Randomnoise
+    {
+    private:
+        double buffer[1024];
+        mutable int index = 0;
+    public:
+        Randomnoise()
+        {
+            std::uniform_real_distribution<double> unif(0., 1.);
+            std::default_random_engine re;
+            for (int i=0; i<1024; i++)
+                buffer[i] = unif(re);
+        }
+        double getNoise() const { index = (index + 1) & 0x3ff; return buffer[index]; }
+    };
 
+protected:
     /// Capacitor value.
     const double C;
 
     /// Transistor parameters.
     //@{
-    const double Vdd;
+    /// Thermal voltage: Ut = kT/q = 8.61734315e-5*T ~ 26mV
+    static constexpr double Ut = 26.0e-3;
+
+    const double Vdd;           ///< Positive supply voltage
     const double Vth;           ///< Threshold voltage
-    const double Ut;            ///< Thermal voltage: Ut = kT/q = 8.61734315e-5*T ~ 26mV
-    const double uCox;          ///< Transconductance coefficient: u*Cox
     const double Vddt;          ///< Vdd - Vth
+    double uCox;                ///< Transconductance coefficient: u*Cox
     //@}
 
     // Derived stuff
@@ -58,38 +86,46 @@ protected:
     /// Fixed point scaling for 16 bit op-amp output.
     const double N16;
 
+    const double voice_voltage_range;
+
     /// Current factor coefficient for op-amp integrators.
-    const double currFactorCoeff;
+    double currFactorCoeff;
 
     /// Lookup tables for gain and summer op-amps in output stage / filter.
     //@{
     unsigned short* mixer[8];       //-V730_NOINIT this is initialized in the derived class constructor
     unsigned short* summer[5];      //-V730_NOINIT this is initialized in the derived class constructor
-    unsigned short* gain_vol[16];   //-V730_NOINIT this is initialized in the derived class constructor
+    unsigned short* volume[16];     //-V730_NOINIT this is initialized in the derived class constructor
-    unsigned short* gain_res[16];   //-V730_NOINIT this is initialized in the derived class constructor
+    unsigned short* resonance[16];  //-V730_NOINIT this is initialized in the derived class constructor
     //@}
 
     /// Reverse op-amp transfer function.
     unsigned short opamp_rev[1 << 16]; //-V730_NOINIT this is initialized in the derived class constructor
 
 private:
-    FilterModelConfig (const FilterModelConfig&) DELETE;
+    Randomnoise rnd;
-    FilterModelConfig& operator= (const FilterModelConfig&) DELETE;
+
+private:
+    FilterModelConfig(const FilterModelConfig&) = delete;
+    FilterModelConfig& operator= (const FilterModelConfig&) = delete;
+
+    inline double getVoiceVoltage(float value, unsigned int env) const
+    {
+        return value * voice_voltage_range + getVoiceDC(env);
+    }
 
 protected:
     /**
      * @param vvr voice voltage range
-     * @param vdv voice DC voltage
      * @param c   capacitor value
-     * @param vdd Vdd
+     * @param vdd Vdd supply voltage
      * @param vth threshold voltage
      * @param ucox u*Cox
-     * @param ominv opamp min voltage
+     * @param opamp_voltage opamp voltage array
-     * @param omaxv opamp max voltage
+     * @param opamp_size opamp voltage array size
      */
     FilterModelConfig(
         double vvr,
-        double vdv,
         double c,
         double vdd,
         double vth,
@@ -98,52 +134,139 @@ protected:
         int opamp_size
     );
 
-    ~FilterModelConfig()
+    ~FilterModelConfig();
+
+    void setUCox(double new_uCox);
+
+    virtual double getVoiceDC(unsigned int env) const = 0;
+
+    /**
+     * The filter summer operates at n ~ 1, and has 5 fundamentally different
+     * input configurations (2 - 6 input "resistors").
+     *
+     * Note that all "on" transistors are modeled as one. This is not
+     * entirely accurate, since the input for each transistor is different,
+     * and transistors are not linear components. However modeling all
+     * transistors separately would be extremely costly.
+     */
+    inline void buildSummerTable(const OpAmp& opampModel)
     {
-        for (int i = 0; i < 8; i++)
+        const double r_N16 = 1. / N16;
-        {
-            delete [] mixer[i];
-        }
 
         for (int i = 0; i < 5; i++)
         {
-            delete [] summer[i];
+            const int idiv = 2 + i;        // 2 - 6 input "resistors".
-        }
+            const int size = idiv << 16;
+            const double n = idiv;
+            const double r_idiv = 1. / idiv;
+            opampModel.reset();
+            summer[i] = new unsigned short[size];
 
-        for (int i = 0; i < 16; i++)
+            for (int vi = 0; vi < size; vi++)
+            {
+                const double vin = vmin + vi * r_N16 * r_idiv; /* vmin .. vmax */
+                summer[i][vi] = getNormalizedValue(opampModel.solve(n, vin));
+            }
+        }
+    }
+
+    /**
+     * The audio mixer operates at n ~ 8/6 (6581) or 8/5 (8580),
+     * and has 8 fundamentally different input configurations
+     * (0 - 7 input "resistors").
+     *
+     * All "on", transistors are modeled as one - see comments above for
+     * the filter summer.
+     */
+    inline void buildMixerTable(const OpAmp& opampModel, double nRatio)
+    {
+        const double r_N16 = 1. / N16;
+
+        for (int i = 0; i < 8; i++)
         {
-            delete [] gain_vol[i];
+            const int idiv = (i == 0) ? 1 : i;
-            delete [] gain_res[i];
+            const int size = (i == 0) ? 1 : i << 16;
+            const double n = i * nRatio;
+            const double r_idiv = 1. / idiv;
+            opampModel.reset();
+            mixer[i] = new unsigned short[size];
+
+            for (int vi = 0; vi < size; vi++)
+            {
+                const double vin = vmin + vi * r_N16 * r_idiv; /* vmin .. vmax */
+                mixer[i][vi] = getNormalizedValue(opampModel.solve(n, vin));
+            }
+        }
+    }
+
+    /**
+     * 4 bit "resistor" ladders in the audio output gain
+     * necessitate 16 gain tables.
+     * From die photographs of the volume "resistor" ladders
+     * it follows that gain ~ vol/12 (6581) or vol/16 (8580)
+     * (assuming ideal op-amps and ideal "resistors").
+     */
+    inline void buildVolumeTable(const OpAmp& opampModel, double nDivisor)
+    {
+        const double r_N16 = 1. / N16;
+
+        for (int n8 = 0; n8 < 16; n8++)
+        {
+            const int size = 1 << 16;
+            const double n = n8 / nDivisor;
+            opampModel.reset();
+            volume[n8] = new unsigned short[size];
+
+            for (int vi = 0; vi < size; vi++)
+            {
+                const double vin = vmin + vi * r_N16; /* vmin .. vmax */
+                volume[n8][vi] = getNormalizedValue(opampModel.solve(n, vin));
+            }
+        }
+    }
+
+    /**
+     * 4 bit "resistor" ladders in the bandpass resonance gain
+     * necessitate 16 gain tables.
+     * From die photographs of the bandpass "resistor" ladders
+     * it follows that 1/Q ~ ~res/8 (6581) or 2^((4 - res)/8) (8580)
+     * (assuming ideal op-amps and ideal "resistors").
+     */
+    inline void buildResonanceTable(const OpAmp& opampModel, const double resonance_n[16])
+    {
+        const double r_N16 = 1. / N16;
+
+        for (int n8 = 0; n8 < 16; n8++)
+        {
+            const int size = 1 << 16;
+            opampModel.reset();
+            resonance[n8] = new unsigned short[size];
+
+            for (int vi = 0; vi < size; vi++)
+            {
+                const double vin = vmin + vi * r_N16; /* vmin .. vmax */
+                resonance[n8][vi] = getNormalizedValue(opampModel.solve(resonance_n[n8], vin));
+            }
         }
     }
 
 public:
-    unsigned short** getGainVol() { return gain_vol; }
+    unsigned short** getVolume() { return volume; }
-    unsigned short** getGainRes() { return gain_res; }
+    unsigned short** getResonance() { return resonance; }
     unsigned short** getSummer() { return summer; }
     unsigned short** getMixer() { return mixer; }
 
-    /**
-     * The digital range of one voice is 20 bits; create a scaling term
-     * for multiplication which fits in 11 bits.
-     */
-    int getVoiceScaleS11() const { return static_cast<int>((norm * ((1 << 11) - 1)) * voice_voltage_range); }
-
-    /**
-     * The "zero" output level of the voices.
-     */
-    int getNormalizedVoiceDC() const { return static_cast<int>(N16 * (voice_DC_voltage - vmin)); }
-
     inline unsigned short getOpampRev(int i) const { return opamp_rev[i]; }
     inline double getVddt() const { return Vddt; }
     inline double getVth() const { return Vth; }
 
     // helper functions
+
     inline unsigned short getNormalizedValue(double value) const
     {
         const double tmp = N16 * (value - vmin);
-        assert(tmp > -0.5 && tmp < 65535.5);
+        assert(tmp >= 0. && tmp <= 65535.);
-        return static_cast<unsigned short>(tmp + 0.5);
+        return static_cast<unsigned short>(tmp + rnd.getNoise());
     }
 
     inline unsigned short getNormalizedCurrentFactor(double wl) const
@@ -153,11 +276,17 @@ public:
         return static_cast<unsigned short>(tmp + 0.5);
     }
 
-    inline unsigned short getNVmin() const {
+    inline unsigned short getNVmin() const
+    {
         const double tmp = N16 * vmin;
         assert(tmp > -0.5 && tmp < 65535.5);
         return static_cast<unsigned short>(tmp + 0.5);
     }
+
+    inline int getNormalizedVoice(float value, unsigned int env) const
+    {
+        return static_cast<int>(getNormalizedValue(getVoiceVoltage(value, env)));
+    }
 };
 
 } // namespace reSIDfp

src/sound/resid-fp/FilterModelConfig6581.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2023 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2010 Dag Lem
  *
@@ -22,28 +22,20 @@
 
 #include "FilterModelConfig6581.h"
 
-#include <cmath>
-
 #include "Integrator6581.h"
 #include "OpAmp.h"
 
+#include "sidcxx11.h"
+
+#include <algorithm>
+#include <mutex>
+#include <thread>
+#include <cmath>
+
 namespace reSIDfp
 {
 
-#ifndef HAVE_CXX11
+constexpr unsigned int OPAMP_SIZE = 33;
-/**
- * Compute log(1+x) without losing precision for small values of x
- *
- * @note when compiling with -ffastm-math the compiler will
- * optimize the expression away leaving a plain log(1. + x)
- */
-inline double log1p(double x)
-{
-    return log(1. + x) - (((1. + x) - 1.) - x) / (1. + x);
-}
-#endif
-
-const unsigned int OPAMP_SIZE = 33;
 
 /**
  * This is the SID 6581 op-amp voltage transfer function, measured on
@@ -51,7 +43,7 @@ const unsigned int OPAMP_SIZE = 33;
  * All measured chips have op-amps with output voltages (and thus input
  * voltages) within the range of 0.81V - 10.31V.
  */
-const Spline::Point opamp_voltage[OPAMP_SIZE] =
+constexpr Spline::Point opamp_voltage[OPAMP_SIZE] =
 {
   {  0.81, 10.31 },  // Approximate start of actual range
   {  2.40, 10.31 },
@@ -90,8 +82,12 @@ const Spline::Point opamp_voltage[OPAMP_SIZE] =
 
 std::unique_ptr<FilterModelConfig6581> FilterModelConfig6581::instance(nullptr);
 
+std::mutex Instance6581_Lock;
+
 FilterModelConfig6581* FilterModelConfig6581::getInstance()
 {
+    std::lock_guard<std::mutex> lock(Instance6581_Lock);
+
     if (!instance.get())
     {
         instance.reset(new FilterModelConfig6581());
@@ -100,14 +96,32 @@ FilterModelConfig6581* FilterModelConfig6581::getInstance()
     return instance.get();
 }
 
+void FilterModelConfig6581::setFilterRange(double adjustment)
+{
+    // clamp into allowed range
+#ifdef HAVE_CXX17
+     adjustment = std::clamp(adjustment, 0.0, 1.0);
+#else
+     adjustment = std::max(std::min(adjustment, 1.0), 0.);
+#endif
+
+     // Get the new uCox value, in the range [1,40]
+     const double new_uCox = (1. + 39. * adjustment) * 1e-6;
+
+    // Ignore small changes
+    if (std::abs(uCox - new_uCox) < 1e-12)
+        return;
+
+    setUCox(new_uCox);
+}
+
 FilterModelConfig6581::FilterModelConfig6581() :
     FilterModelConfig(
-        1.5,     // voice voltage range
+        1.5,                    // voice voltage range FIXME should theoretically be ~3,571V
-        5.075,   // voice DC voltage
+        470e-12,                // capacitor value
-        470e-12, // capacitor value
+        12. * VOLTAGE_SKEW,     // Vdd
-        12.18,   // Vdd
+        1.31,                   // Vth
-        1.31,    // Vth
+        20e-6,                  // uCox
-        20e-6,   // uCox
         opamp_voltage,
         OPAMP_SIZE
     ),
@@ -119,190 +133,144 @@ FilterModelConfig6581::FilterModelConfig6581() :
 {
     dac.kinkedDac(MOS6581);
 
-    // Create lookup tables for gains / summers.
-
-#ifndef _OPENMP
-    OpAmp opampModel(
-        std::vector<Spline::Point>(
-            std::begin(opamp_voltage),
-            std::end(opamp_voltage)),
-        Vddt,
-        vmin,
-        vmax);
-#endif
-
-//    #pragma omp parallel sections
     {
-//        #pragma omp section
+        Dac envDac(8);
+        envDac.kinkedDac(MOS6581);
+        for(int i=0; i<256; i++)
         {
-#ifdef _OPENMP
+            const double envI = envDac.getOutput(i);
-            OpAmp opampModel(
+            voiceDC[i] = 5. * VOLTAGE_SKEW + (0.2143 * envI);
-                std::vector<Spline::Point>(
-                    std::begin(opamp_voltage),
-                    std::end(opamp_voltage)),
-                Vddt,
-                vmin,
-                vmax);
-#endif
-            // The filter summer operates at n ~ 1, and has 5 fundamentally different
-            // input configurations (2 - 6 input "resistors").
-            //
-            // Note that all "on" transistors are modeled as one. This is not
-            // entirely accurate, since the input for each transistor is different,
-            // and transistors are not linear components. However modeling all
-            // transistors separately would be extremely costly.
-            for (int i = 0; i < 5; i++)
-            {
-                const int idiv = 2 + i;        // 2 - 6 input "resistors".
-                const int size = idiv << 16;
-                const double n = idiv;
-                opampModel.reset();
-                summer[i] = new unsigned short[size];
-
-                for (int vi = 0; vi < size; vi++)
-                {
-                    const double vin = vmin + vi / N16 / idiv; /* vmin .. vmax */
-                    summer[i][vi] = getNormalizedValue(opampModel.solve(n, vin));
-                }
-            }
-        }
-
-//        #pragma omp section
-        {
-#ifdef _OPENMP
-            OpAmp opampModel(
-                std::vector<Spline::Point>(
-                    std::begin(opamp_voltage),
-                    std::end(opamp_voltage)),
-                Vddt,
-                vmin,
-                vmax);
-#endif
-            // The audio mixer operates at n ~ 8/6, and has 8 fundamentally different
-            // input configurations (0 - 7 input "resistors").
-            //
-            // All "on", transistors are modeled as one - see comments above for
-            // the filter summer.
-            for (int i = 0; i < 8; i++)
-            {
-                const int idiv = (i == 0) ? 1 : i;
-                const int size = (i == 0) ? 1 : i << 16;
-                const double n = i * 8.0 / 6.0;
-                opampModel.reset();
-                mixer[i] = new unsigned short[size];
-
-                for (int vi = 0; vi < size; vi++)
-                {
-                    const double vin = vmin + vi / N16 / idiv; /* vmin .. vmax */
-                    mixer[i][vi] = getNormalizedValue(opampModel.solve(n, vin));
-                }
-            }
-        }
-
-//        #pragma omp section
-        {
-#ifdef _OPENMP
-            OpAmp opampModel(
-                std::vector<Spline::Point>(
-                    std::begin(opamp_voltage),
-                    std::end(opamp_voltage)),
-                Vddt,
-                vmin,
-                vmax);
-#endif
-            // 4 bit "resistor" ladders in the audio output gain
-            // necessitate 16 gain tables.
-            // From die photographs of the volume "resistor" ladders
-            // it follows that gain ~ vol/12 (assuming ideal
-            // op-amps and ideal "resistors").
-            for (int n8 = 0; n8 < 16; n8++)
-            {
-                const int size = 1 << 16;
-                const double n = n8 / 12.0;
-                opampModel.reset();
-                gain_vol[n8] = new unsigned short[size];
-
-                for (int vi = 0; vi < size; vi++)
-                {
-                    const double vin = vmin + vi / N16; /* vmin .. vmax */
-                    gain_vol[n8][vi] = getNormalizedValue(opampModel.solve(n, vin));
-                }
-            }
-        }
-
-//        #pragma omp section
-        {
-#ifdef _OPENMP
-            OpAmp opampModel(
-                std::vector<Spline::Point>(
-                    std::begin(opamp_voltage),
-                    std::end(opamp_voltage)),
-                Vddt,
-                vmin,
-                vmax);
-#endif
-            // 4 bit "resistor" ladders in the bandpass resonance gain
-            // necessitate 16 gain tables.
-            // From die photographs of the bandpass "resistor" ladders
-            // it follows that 1/Q ~ ~res/8 (assuming ideal
-            // op-amps and ideal "resistors").
-            for (int n8 = 0; n8 < 16; n8++)
-            {
-                const int size = 1 << 16;
-                const double n = (~n8 & 0xf) / 8.0;
-                opampModel.reset();
-                gain_res[n8] = new unsigned short[size];
-
-                for (int vi = 0; vi < size; vi++)
-                {
-                    const double vin = vmin + vi / N16; /* vmin .. vmax */
-                    gain_res[n8][vi] = getNormalizedValue(opampModel.solve(n, vin));
-                }
-            }
-        }
-
-//        #pragma omp section
-        {
-            const double nVddt = N16 * (Vddt - vmin);
-
-            for (unsigned int i = 0; i < (1 << 16); i++)
-            {
-                // The table index is right-shifted 16 times in order to fit in
-                // 16 bits; the argument to sqrt is thus multiplied by (1 << 16).
-                const double tmp = nVddt - sqrt(static_cast<double>(i << 16));
-                assert(tmp > -0.5 && tmp < 65535.5);
-                vcr_nVg[i] = static_cast<unsigned short>(tmp + 0.5);
-            }
-        }
-
-//        #pragma omp section
-        {
-            //  EKV model:
-            //
-            //  Ids = Is * (if - ir)
-            //  Is = (2 * u*Cox * Ut^2)/k * W/L
-            //  if = ln^2(1 + e^((k*(Vg - Vt) - Vs)/(2*Ut))
-            //  ir = ln^2(1 + e^((k*(Vg - Vt) - Vd)/(2*Ut))
-
-            // moderate inversion characteristic current
-            const double Is = (2. * uCox * Ut * Ut) * WL_vcr;
-
-            // Normalized current factor for 1 cycle at 1MHz.
-            const double N15 = norm * ((1 << 15) - 1);
-            const double n_Is = N15 * 1.0e-6 / C * Is;
-
-            // kVgt_Vx = k*(Vg - Vt) - Vx
-            // I.e. if k != 1.0, Vg must be scaled accordingly.
-            for (int kVgt_Vx = 0; kVgt_Vx < (1 << 16); kVgt_Vx++)
-            {
-                const double log_term = log1p(exp((kVgt_Vx / N16) / (2. * Ut)));
-                // Scaled by m*2^15
-                const double tmp = n_Is * log_term * log_term;
-                assert(tmp > -0.5 && tmp < 65535.5);
-                vcr_n_Ids_term[kVgt_Vx] = static_cast<unsigned short>(tmp + 0.5);
-            }
         }
     }
+
+    // Create lookup tables for gains / summers.
+
+    //
+    // We spawn six threads to calculate these tables in parallel
+    //
+    auto filterSummer = [this]
+    {
+        OpAmp opampModel(
+            std::vector<Spline::Point>(
+                std::begin(opamp_voltage),
+                std::end(opamp_voltage)),
+            Vddt,
+            vmin,
+            vmax);
+
+        buildSummerTable(opampModel);
+    };
+
+    auto filterMixer = [this]
+    {
+        OpAmp opampModel(
+            std::vector<Spline::Point>(
+                std::begin(opamp_voltage),
+                std::end(opamp_voltage)),
+            Vddt,
+            vmin,
+            vmax);
+
+        buildMixerTable(opampModel, 8.0 / 6.0);
+    };
+
+    auto filterGain = [this]
+    {
+        OpAmp opampModel(
+            std::vector<Spline::Point>(
+                std::begin(opamp_voltage),
+                std::end(opamp_voltage)),
+            Vddt,
+            vmin,
+            vmax);
+
+        buildVolumeTable(opampModel, 12.0);
+    };
+
+    auto filterResonance = [this]
+    {
+        OpAmp opampModel(
+            std::vector<Spline::Point>(
+                std::begin(opamp_voltage),
+                std::end(opamp_voltage)),
+            Vddt,
+            vmin,
+            vmax);
+
+        // build temp n table
+        double resonance_n[16];
+        for (int n8 = 0; n8 < 16; n8++)
+        {
+            resonance_n[n8] = (~n8 & 0xf) / 8.0;
+        }
+
+        buildResonanceTable(opampModel, resonance_n);
+    };
+
+    auto filterVcrVg = [this]
+    {
+        const double nVddt = N16 * (Vddt - vmin);
+
+        for (unsigned int i = 0; i < (1 << 16); i++)
+        {
+            // The table index is right-shifted 16 times in order to fit in
+            // 16 bits; the argument to sqrt is thus multiplied by (1 << 16).
+            const double tmp = nVddt - std::sqrt(static_cast<double>(i << 16));
+            assert(tmp > -0.5 && tmp < 65535.5);
+            vcr_nVg[i] = static_cast<unsigned short>(tmp + 0.5);
+        }
+    };
+
+    auto filterVcrIds = [this]
+    {
+        //  EKV model:
+        //
+        //  Ids = Is * (if - ir)
+        //  Is = (2 * u*Cox * Ut^2)/k * W/L
+        //  if = ln^2(1 + e^((k*(Vg - Vt) - Vs)/(2*Ut))
+        //  ir = ln^2(1 + e^((k*(Vg - Vt) - Vd)/(2*Ut))
+
+        // moderate inversion characteristic current
+        // will be multiplied by uCox later
+        const double Is = (2. * Ut * Ut) * WL_vcr;
+
+        // Normalized current factor for 1 cycle at 1MHz.
+        const double N15 = norm * ((1 << 15) - 1);
+        const double n_Is = N15 * 1.0e-6 / C * Is;
+
+        // kVgt_Vx = k*(Vg - Vt) - Vx
+        // I.e. if k != 1.0, Vg must be scaled accordingly.
+        const double  r_N16_2Ut = 1.0 / (N16 * 2.0 * Ut);
+        for (int i = 0; i < (1 << 16); i++)
+        {
+            const int kVgt_Vx = i - (1 << 15);
+            const double log_term = std::log1p(std::exp(kVgt_Vx * r_N16_2Ut));
+            // Scaled by m*2^15
+            vcr_n_Ids_term[i] = n_Is * log_term * log_term;
+        }
+    };
+
+#if defined(HAVE_CXX20) && defined(__cpp_lib_jthread)
+    using sidThread = std::jthread;
+#else
+    using sidThread = std::thread;
+#endif
+
+    sidThread thdSummer(filterSummer);
+    sidThread thdMixer(filterMixer);
+    sidThread thdGain(filterGain);
+    sidThread thdResonance(filterResonance);
+    sidThread thdVcrVg(filterVcrVg);
+    sidThread thdVcrIds(filterVcrIds);
+
+#if !defined(HAVE_CXX20) || !defined(__cpp_lib_jthread)
+    thdSummer.join();
+    thdMixer.join();
+    thdGain.join();
+    thdResonance.join();
+    thdVcrVg.join();
+    thdVcrIds.join();
+#endif
 }
 
 unsigned short* FilterModelConfig6581::getDAC(double adjustment) const
@@ -314,15 +282,10 @@ unsigned short* FilterModelConfig6581::getDAC(double adjustment) const
     for (unsigned int i = 0; i < (1 << DAC_BITS); i++)
     {
         const double fcd = dac.getOutput(i);
-        f0_dac[i] = getNormalizedValue(dac_zero + fcd * dac_scale / (1 << DAC_BITS));
+        f0_dac[i] = getNormalizedValue(dac_zero + fcd * dac_scale);
     }
 
     return f0_dac;
 }
 
-std::unique_ptr<Integrator6581> FilterModelConfig6581::buildIntegrator()
-{
-    return MAKE_UNIQUE(Integrator6581, this, WL_snake);
-}
-
 } // namespace reSIDfp

src/sound/resid-fp/FilterModelConfig6581.h

@@ -41,17 +41,18 @@ class Integrator6581;
  */
 class FilterModelConfig6581 final : public FilterModelConfig
 {
-private:
-    static const unsigned int DAC_BITS = 11;
-
 private:
     static std::unique_ptr<FilterModelConfig6581> instance;
     // This allows access to the private constructor
-#ifdef HAVE_CXX11
     friend std::unique_ptr<FilterModelConfig6581>::deleter_type;
-#else
+
-    friend class std::auto_ptr<FilterModelConfig6581>;
+private:
-#endif
+    static constexpr unsigned int DAC_BITS = 11;
+
+    /**
+     * Power bricks generate voltages slightly out of spec
+     */
+    static constexpr double VOLTAGE_SKEW = 1.015;
 
     /// Transistor parameters.
     //@{
@@ -68,21 +69,36 @@ private:
     /// DAC lookup table
     Dac dac;
 
-    /// VCR - 6581 only.
+    /// Voltage Controlled Resistors
     //@{
     unsigned short vcr_nVg[1 << 16];
-    unsigned short vcr_n_Ids_term[1 << 16];
+    double vcr_n_Ids_term[1 << 16];
     //@}
 
+    // Voice DC offset LUT
+    double voiceDC[256];
+
 private:
     double getDacZero(double adjustment) const { return dac_zero + (1. - adjustment); }
 
     FilterModelConfig6581();
-    ~FilterModelConfig6581() DEFAULT;
+    ~FilterModelConfig6581() = default;
+
+protected:
+    /**
+     * On 6581 the DC offset varies between ~5.0V and ~5.214V depending on
+     * the envelope value.
+     */
+    inline double getVoiceDC(unsigned int env) const override
+    {
+        return voiceDC[env];
+    }
 
 public:
     static FilterModelConfig6581* getInstance();
 
+    void setFilterRange(double adjustment);
+
     /**
      * Construct an 11 bit cutoff frequency DAC output voltage table.
      * Ownership is transferred to the requester which becomes responsible
@@ -93,17 +109,17 @@ public:
      */
     unsigned short* getDAC(double adjustment) const;
 
-    /**
+    inline double getWL_snake() const { return WL_snake; }
-     * Construct an integrator solver.
-     *
-     * @return the integrator
-     */
-    std::unique_ptr<Integrator6581> buildIntegrator();
 
     inline unsigned short getVcr_nVg(int i) const { return vcr_nVg[i]; }
-    inline unsigned short getVcr_n_Ids_term(int i) const { return vcr_n_Ids_term[i]; }
+    inline unsigned short getVcr_n_Ids_term(int i) const
+    {
+        const double tmp = vcr_n_Ids_term[i] * uCox;
+        assert(tmp > -0.5 && tmp < 65535.5);
+        return static_cast<unsigned short>(tmp + 0.5);
+    }
     // only used if SLOPE_FACTOR is defined
-    inline double getUt() const { return Ut; }
+    inline constexpr double getUt() const { return Ut; }
     inline double getN16() const { return N16; }
 };
 

src/sound/resid-fp/FilterModelConfig8580.cpp

@@ -25,6 +25,10 @@
 #include "Integrator8580.h"
 #include "OpAmp.h"
 
+#include "sidcxx11.h"
+
+#include <mutex>
+#include <thread>
 
 namespace reSIDfp
 {
@@ -57,7 +61,7 @@ namespace reSIDfp
  *  E   Rf|R2     RC
  *  F   Rf|R3     RC
  */
-const double resGain[16] =
+constexpr double resGain[16] =
 {
                         1.4/1.0, //      Rf/Ri   1.4
     ((1.4*15.3)/(1.4+15.3))/1.0, // (Rf|R1)/Ri   1.28263
@@ -77,13 +81,13 @@ const double resGain[16] =
       ((1.4*4.7)/(1.4+4.7))/2.8, // (Rf|R3)/RC   0.385246
 };
 
-const unsigned int OPAMP_SIZE = 21;
+constexpr unsigned int OPAMP_SIZE = 21;
 
 /**
  * This is the SID 8580 op-amp voltage transfer function, measured on
  * CAP1B/CAP1A on a chip marked CSG 8580R5 1690 25.
  */
-const Spline::Point opamp_voltage[OPAMP_SIZE] =
+constexpr Spline::Point opamp_voltage[OPAMP_SIZE] =
 {
     {  1.30,  8.91 },  // Approximate start of actual range
     {  4.76,  8.91 },
@@ -110,8 +114,12 @@ const Spline::Point opamp_voltage[OPAMP_SIZE] =
 
 std::unique_ptr<FilterModelConfig8580> FilterModelConfig8580::instance(nullptr);
 
+std::mutex Instance8580_Lock;
+
 FilterModelConfig8580* FilterModelConfig8580::getInstance()
 {
+    std::lock_guard<std::mutex> lock(Instance8580_Lock);
+
     if (!instance.get())
     {
         instance.reset(new FilterModelConfig8580());
@@ -122,161 +130,89 @@ FilterModelConfig8580* FilterModelConfig8580::getInstance()
 
 FilterModelConfig8580::FilterModelConfig8580() :
     FilterModelConfig(
-        0.30,   // voice voltage range FIXME measure
+        0.24,               // voice voltage range FIXME should theoretically be ~0,474V
-        4.84,   // voice DC voltage FIXME measure
+        22e-9,              // capacitor value
-        22e-9,  // capacitor value
+        9. * VOLTAGE_SKEW,  // Vdd
-        9.09,   // Vdd
+        0.80,               // Vth
-        0.80,   // Vth
+        100e-6,             // uCox
-        100e-6, // uCox
         opamp_voltage,
         OPAMP_SIZE
     )
 {
     // Create lookup tables for gains / summers.
-#ifndef _OPENMP
-    OpAmp opampModel(
-        std::vector<Spline::Point>(
-            std::begin(opamp_voltage),
-            std::end(opamp_voltage)),
-        Vddt,
-        vmin,
-        vmax);
-#endif
 
-//    #pragma omp parallel sections
+    //
+    // We spawn four threads to calculate these tables in parallel
+    //
+    auto filterSummer = [this]
     {
-//        #pragma omp section
+        OpAmp opampModel(
-        {
+            std::vector<Spline::Point>(
-#ifdef _OPENMP
+                std::begin(opamp_voltage),
-            OpAmp opampModel(
+                std::end(opamp_voltage)),
-                std::vector<Spline::Point>(
+            Vddt,
-                    std::begin(opamp_voltage),
+            vmin,
-                    std::end(opamp_voltage)),
+            vmax);
-                Vddt,
+
-                vmin,
+        buildSummerTable(opampModel);
-                vmax);
+    };
+
+    auto filterMixer = [this]
+    {
+        OpAmp opampModel(
+            std::vector<Spline::Point>(
+                std::begin(opamp_voltage),
+                std::end(opamp_voltage)),
+            Vddt,
+            vmin,
+            vmax);
+
+        buildMixerTable(opampModel, 8.0 / 5.0);
+    };
+
+    auto filterGain = [this]
+    {
+        OpAmp opampModel(
+            std::vector<Spline::Point>(
+                std::begin(opamp_voltage),
+                std::end(opamp_voltage)),
+            Vddt,
+            vmin,
+            vmax);
+
+        buildVolumeTable(opampModel, 16.0);
+    };
+
+    auto filterResonance = [this]
+    {
+        OpAmp opampModel(
+            std::vector<Spline::Point>(
+                std::begin(opamp_voltage),
+                std::end(opamp_voltage)),
+            Vddt,
+            vmin,
+            vmax);
+
+        buildResonanceTable(opampModel, resGain);
+    };
+
+#if defined(HAVE_CXX20) && defined(__cpp_lib_jthread)
+    using sidThread = std::jthread;
+#else
+    using sidThread = std::thread;
 #endif
-            // The filter summer operates at n ~ 1, and has 5 fundamentally different
-            // input configurations (2 - 6 input "resistors").
-            //
-            // Note that all "on" transistors are modeled as one. This is not
-            // entirely accurate, since the input for each transistor is different,
-            // and transistors are not linear components. However modeling all
-            // transistors separately would be extremely costly.
-            for (int i = 0; i < 5; i++)
-            {
-                const int idiv = 2 + i;        // 2 - 6 input "resistors".
-                const int size = idiv << 16;
-                const double n = idiv;
-                opampModel.reset();
-                summer[i] = new unsigned short[size];
 
-                for (int vi = 0; vi < size; vi++)
+    sidThread thdSummer(filterSummer);
-                {
+    sidThread thdMixer(filterMixer);
-                    const double vin = vmin + vi / N16 / idiv; /* vmin .. vmax */
+    sidThread thdGain(filterGain);
-                    summer[i][vi] = getNormalizedValue(opampModel.solve(n, vin));
+    sidThread thdResonance(filterResonance);
-                }
-            }
-        }
 
-//        #pragma omp section
+#if !defined(HAVE_CXX20) || !defined(__cpp_lib_jthread)
-        {
+    thdSummer.join();
-#ifdef _OPENMP
+    thdMixer.join();
-            OpAmp opampModel(
+    thdGain.join();
-                std::vector<Spline::Point>(
+    thdResonance.join();
-                    std::begin(opamp_voltage),
-                    std::end(opamp_voltage)),
-                Vddt,
-                vmin,
-                vmax);
 #endif
-            // The audio mixer operates at n ~ 8/5, and has 8 fundamentally different
-            // input configurations (0 - 7 input "resistors").
-            //
-            // All "on", transistors are modeled as one - see comments above for
-            // the filter summer.
-            for (int i = 0; i < 8; i++)
-            {
-                const int idiv = (i == 0) ? 1 : i;
-                const int size = (i == 0) ? 1 : i << 16;
-                const double n = i * 8.0 / 5.0;
-                opampModel.reset();
-                mixer[i] = new unsigned short[size];
-
-                for (int vi = 0; vi < size; vi++)
-                {
-                    const double vin = vmin + vi / N16 / idiv; /* vmin .. vmax */
-                    mixer[i][vi] = getNormalizedValue(opampModel.solve(n, vin));
-                }
-            }
-        }
-
-//        #pragma omp section
-        {
-#ifdef _OPENMP
-            OpAmp opampModel(
-                std::vector<Spline::Point>(
-                    std::begin(opamp_voltage),
-                    std::end(opamp_voltage)),
-                Vddt,
-                vmin,
-                vmax);
-#endif
-            // 4 bit "resistor" ladders in the audio output gain
-            // necessitate 16 gain tables.
-            // From die photographs of the volume "resistor" ladders
-            // it follows that gain ~ vol/16 (assuming ideal
-            // op-amps and ideal "resistors").
-            for (int n8 = 0; n8 < 16; n8++)
-            {
-                const int size = 1 << 16;
-                const double n = n8 / 16.0;
-                opampModel.reset();
-                gain_vol[n8] = new unsigned short[size];
-
-                for (int vi = 0; vi < size; vi++)
-                {
-                    const double vin = vmin + vi / N16; /* vmin .. vmax */
-                    gain_vol[n8][vi] = getNormalizedValue(opampModel.solve(n, vin));
-                }
-            }
-        }
-
-//        #pragma omp section
-        {
-#ifdef _OPENMP
-            OpAmp opampModel(
-                std::vector<Spline::Point>(
-                    std::begin(opamp_voltage),
-                    std::end(opamp_voltage)),
-                Vddt,
-                vmin,
-                vmax);
-#endif
-            // 4 bit "resistor" ladders in the bandpass resonance gain
-            // necessitate 16 gain tables.
-            // From die photographs of the bandpass "resistor" ladders
-            // it follows that 1/Q ~ 2^((4 - res)/8) (assuming ideal
-            // op-amps and ideal "resistors").
-            for (int n8 = 0; n8 < 16; n8++)
-            {
-                const int size = 1 << 16;
-                opampModel.reset();
-                gain_res[n8] = new unsigned short[size];
-
-                for (int vi = 0; vi < size; vi++)
-                {
-                    const double vin = vmin + vi / N16; /* vmin .. vmax */
-                    gain_res[n8][vi] = getNormalizedValue(opampModel.solve(resGain[n8], vin));
-                }
-            }
-        }
-    }
-}
-
-std::unique_ptr<Integrator8580> FilterModelConfig8580::buildIntegrator()
-{
-    return MAKE_UNIQUE(Integrator8580, this);
 }
 
 } // namespace reSIDfp

src/sound/resid-fp/FilterModelConfig8580.h

@@ -42,25 +42,30 @@ class FilterModelConfig8580 final : public FilterModelConfig
 private:
     static std::unique_ptr<FilterModelConfig8580> instance;
     // This allows access to the private constructor
-#ifdef HAVE_CXX11
     friend std::unique_ptr<FilterModelConfig8580>::deleter_type;
-#else
+
-    friend class std::auto_ptr<FilterModelConfig8580>;
+private:
-#endif
+    /**
+     * Reference voltage generated from Vcc by a voltage divider
+     */
+    static constexpr double Vref = 4.75;
+
+    /**
+     * Power bricks generate voltages slightly out of spec
+     */
+    static constexpr double VOLTAGE_SKEW = 1.01;
 
 private:
     FilterModelConfig8580();
-    ~FilterModelConfig8580() DEFAULT;
+    ~FilterModelConfig8580() = default;
+
+protected:
+    inline double getVoiceDC(unsigned int) const override { return getVref(); }
 
 public:
     static FilterModelConfig8580* getInstance();
 
-    /**
+    inline constexpr double getVref() const { return Vref * VOLTAGE_SKEW; }
-     * Construct an integrator solver.
-     *
-     * @return the integrator
-     */
-    std::unique_ptr<Integrator8580> buildIntegrator();
 };
 
 } // namespace reSIDfp

src/sound/resid-fp/Integrator.h

@@ -0,0 +1,47 @@
+/*
+ * This file is part of libsidplayfp, a SID player engine.
+ *
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2007-2010 Antti Lankila
+ * Copyright 2004, 2010 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
+ */
+
+#ifndef INTEGRATOR_H
+#define INTEGRATOR_H
+
+namespace reSIDfp
+{
+
+class Integrator
+{
+protected:
+    mutable int vx;
+    mutable int vc;
+
+    Integrator() :
+        vx(0),
+        vc(0) {}
+
+public:
+    virtual int solve(int vi) const = 0;
+
+    virtual ~Integrator() = default;
+};
+
+} // namespace reSIDfp
+
+#endif

src/sound/resid-fp/Integrator6581.cpp

@@ -18,8 +18,80 @@
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  */
 
-#define INTEGRATOR_CPP
-
 #include "Integrator6581.h"
 
-// This is needed when compiling with --disable-inline
+#ifdef SLOPE_FACTOR
+#  include <cmath>
+#  include "sidcxx11.h"
+#endif
+
+namespace reSIDfp
+{
+
+int Integrator6581::solve(int vi) const
+{
+    // Make sure Vgst>0 so we're not in subthreshold mode
+    assert(vx < nVddt);
+
+    // Check that transistor is actually in triode mode
+    // Vds < Vgs - Vth
+    assert(vi < nVddt);
+
+    // "Snake" voltages for triode mode calculation.
+    const unsigned int Vgst = nVddt - vx;
+    const unsigned int Vgdt = nVddt - vi;
+
+    const unsigned int Vgst_2 = Vgst * Vgst;
+    const unsigned int Vgdt_2 = Vgdt * Vgdt;
+
+    // "Snake" current, scaled by (1/m)*2^13*m*2^16*m*2^16*2^-15 = m*2^30
+    const int n_I_snake = fmc.getNormalizedCurrentFactor(wlSnake) * (static_cast<int>(Vgst_2 - Vgdt_2) >> 15);
+
+    // VCR gate voltage.       // Scaled by m*2^16
+    // Vg = Vddt - sqrt(((Vddt - Vw)^2 + Vgdt^2)/2)
+    const int nVg = static_cast<int>(fmc.getVcr_nVg((nVddt_Vw_2 + (Vgdt_2 >> 1)) >> 16));
+#ifdef SLOPE_FACTOR
+    const double nVp = static_cast<double>(nVg - nVt) / n; // Pinch-off voltage
+    const int kVgt = static_cast<int>(nVp + 0.5) - nVmin;
+#else
+    const int kVgt = (nVg - nVt) - nVmin;
+#endif
+
+    // VCR voltages for EKV model table lookup.
+    const int kVgt_Vs = (kVgt - vx) + (1 << 15);
+    assert((kVgt_Vs >= 0) && (kVgt_Vs < (1 << 16)));
+    const int kVgt_Vd = (kVgt - vi) + (1 << 15);
+    assert((kVgt_Vd >= 0) && (kVgt_Vd < (1 << 16)));
+
+    // VCR current, scaled by m*2^15*2^15 = m*2^30
+    const unsigned int If = static_cast<unsigned int>(fmc.getVcr_n_Ids_term(kVgt_Vs)) << 15;
+    const unsigned int Ir = static_cast<unsigned int>(fmc.getVcr_n_Ids_term(kVgt_Vd)) << 15;
+#ifdef SLOPE_FACTOR
+    const double iVcr = static_cast<double>(If - Ir);
+    const int n_I_vcr = static_cast<int>(iVcr * n);
+#else
+    const int n_I_vcr = If - Ir;
+#endif
+
+#ifdef SLOPE_FACTOR
+    // estimate new slope factor based on gate voltage
+    constexpr double gamma = 1.0;   // body effect factor
+    constexpr double phi = 0.8;     // bulk Fermi potential
+    const double Vp = nVp / fmc.getN16();
+    n = 1. + (gamma / (2. * std::sqrt(Vp + phi + 4. * fmc.getUt())));
+    assert((n > 1.2) && (n < 1.8));
+#endif
+
+    // Change in capacitor charge.
+    vc += n_I_snake + n_I_vcr;
+
+    // vx = g(vc)
+    const int tmp = (vc >> 15) + (1 << 15);
+    assert(tmp < (1 << 16));
+    vx = fmc.getOpampRev(tmp);
+
+    // Return vo.
+    return vx - (vc >> 14);
+}
+
+} // namespace reSIDfp

src/sound/resid-fp/Integrator6581.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2022 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2023 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004, 2010 Dag Lem <resid@nimrod.no>
  *
@@ -23,6 +23,7 @@
 #ifndef INTEGRATOR6581_H
 #define INTEGRATOR6581_H
 
+#include "Integrator.h"
 #include "FilterModelConfig6581.h"
 
 #include <stdint.h>
@@ -33,10 +34,6 @@
 // actually produces worse results, needs investigation
 //#define SLOPE_FACTOR
 
-#ifdef SLOPE_FACTOR
-#  include <cmath>
-#endif
-
 #include "siddefs-fp.h"
 
 namespace reSIDfp
@@ -164,12 +161,10 @@ namespace reSIDfp
  *
  *     Vg = nVddt - sqrt(((nVddt - vi)^2 + (nVddt - Vw)^2)/2)
  */
-class Integrator6581
+class Integrator6581 : public Integrator
 {
 private:
-    unsigned int nVddt_Vw_2;
+    const double wlSnake;
-    mutable int vx;
-    mutable int vc;
 
 #ifdef SLOPE_FACTOR
     // Slope factor n = 1/k
@@ -177,109 +172,32 @@ private:
     // k = Cox/(Cox+Cdep) ~ 0.7 (depends on gate voltage)
     mutable double n;
 #endif
+
+    unsigned int nVddt_Vw_2;
+
     const unsigned short nVddt;
     const unsigned short nVt;
     const unsigned short nVmin;
-    const unsigned short nSnake;
 
-    const FilterModelConfig6581* fmc;
+    FilterModelConfig6581& fmc;
 
 public:
-    Integrator6581(const FilterModelConfig6581* fmc,
+    Integrator6581(FilterModelConfig6581& fmc) :
-               double WL_snake) :
+        wlSnake(fmc.getWL_snake()),
-        nVddt_Vw_2(0),
-        vx(0),
-        vc(0),
 #ifdef SLOPE_FACTOR
         n(1.4),
 #endif
-        nVddt(fmc->getNormalizedValue(fmc->getVddt())),
+        nVddt_Vw_2(0),
-        nVt(fmc->getNormalizedValue(fmc->getVth())),
+        nVddt(fmc.getNormalizedValue(fmc.getVddt())),
-        nVmin(fmc->getNVmin()),
+        nVt(fmc.getNormalizedValue(fmc.getVth())),
-        nSnake(fmc->getNormalizedCurrentFactor(WL_snake)),
+        nVmin(fmc.getNVmin()),
         fmc(fmc) {}
 
     void setVw(unsigned short Vw) { nVddt_Vw_2 = ((nVddt - Vw) * (nVddt - Vw)) >> 1; }
 
-    int solve(int vi) const;
+    int solve(int vi) const override;
 };
 
 } // namespace reSIDfp
 
-#if RESID_INLINING || defined(INTEGRATOR_CPP)
-
-namespace reSIDfp
-{
-
-RESID_INLINE
-int Integrator6581::solve(int vi) const
-{
-    // Make sure Vgst>0 so we're not in subthreshold mode
-    assert(vx < nVddt);
-
-    // Check that transistor is actually in triode mode
-    // Vds < Vgs - Vth
-    assert(vi < nVddt);
-
-    // "Snake" voltages for triode mode calculation.
-    const unsigned int Vgst = nVddt - vx;
-    const unsigned int Vgdt = nVddt - vi;
-
-    const unsigned int Vgst_2 = Vgst * Vgst;
-    const unsigned int Vgdt_2 = Vgdt * Vgdt;
-
-    // "Snake" current, scaled by (1/m)*2^13*m*2^16*m*2^16*2^-15 = m*2^30
-    const int n_I_snake = nSnake * (static_cast<int>(Vgst_2 - Vgdt_2) >> 15);
-
-    // VCR gate voltage.       // Scaled by m*2^16
-    // Vg = Vddt - sqrt(((Vddt - Vw)^2 + Vgdt^2)/2)
-    const int nVg = static_cast<int>(fmc->getVcr_nVg((nVddt_Vw_2 + (Vgdt_2 >> 1)) >> 16));
-#ifdef SLOPE_FACTOR
-    const double nVp = static_cast<double>(nVg - nVt) / n; // Pinch-off voltage
-    const int kVgt = static_cast<int>(nVp + 0.5) - nVmin;
-#else
-    const int kVgt = (nVg - nVt) - nVmin;
-#endif
-
-    // VCR voltages for EKV model table lookup.
-    const int kVgt_Vs = (vx < kVgt) ? kVgt - vx : 0;
-    assert(kVgt_Vs < (1 << 16));
-    const int kVgt_Vd = (vi < kVgt) ? kVgt - vi : 0;
-    assert(kVgt_Vd < (1 << 16));
-
-    // VCR current, scaled by m*2^15*2^15 = m*2^30
-    const unsigned int If = static_cast<unsigned int>(fmc->getVcr_n_Ids_term(kVgt_Vs)) << 15;
-    const unsigned int Ir = static_cast<unsigned int>(fmc->getVcr_n_Ids_term(kVgt_Vd)) << 15;
-#ifdef SLOPE_FACTOR
-    const double iVcr = static_cast<double>(If - Ir);
-    const int n_I_vcr = static_cast<int>(iVcr * n);
-#else
-    const int n_I_vcr = If - Ir;
-#endif
-
-#ifdef SLOPE_FACTOR
-    // estimate new slope factor based on gate voltage
-    const double gamma = 1.0;   // body effect factor
-    const double phi = 0.8;     // bulk Fermi potential
-    const double Vp = nVp / fmc->getN16();
-    n = 1. + (gamma / (2. * sqrt(Vp + phi + 4. * fmc->getUt())));
-    assert((n > 1.2) && (n < 1.8));
-#endif
-
-    // Change in capacitor charge.
-    vc += n_I_snake + n_I_vcr;
-
-    // vx = g(vc)
-    const int tmp = (vc >> 15) + (1 << 15);
-    assert(tmp < (1 << 16));
-    vx = fmc->getOpampRev(tmp);
-
-    // Return vo.
-    return vx - (vc >> 14);
-}
-
-} // namespace reSIDfp
-
-#endif
-
 #endif

src/sound/resid-fp/Integrator8580.cpp

@@ -18,8 +18,36 @@
  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  */
 
-#define INTEGRATOR8580_CPP
-
 #include "Integrator8580.h"
 
-// This is needed when compiling with --disable-inline
+namespace reSIDfp
+{
+
+int Integrator8580::solve(int vi) const
+{
+    // Make sure we're not in subthreshold mode
+    assert(vx < nVgt);
+
+    // DAC voltages
+    const unsigned int Vgst = nVgt - vx;
+    const unsigned int Vgdt = (vi < nVgt) ? nVgt - vi : 0;  // triode/saturation mode
+
+    const unsigned int Vgst_2 = Vgst * Vgst;
+    const unsigned int Vgdt_2 = Vgdt * Vgdt;
+
+    // DAC current, scaled by (1/m)*2^13*m*2^16*m*2^16*2^-15 = m*2^30
+    const int n_I_dac = n_dac * (static_cast<int>(Vgst_2 - Vgdt_2) >> 15);
+
+    // Change in capacitor charge.
+    vc += n_I_dac;
+
+    // vx = g(vc)
+    const int tmp = (vc >> 15) + (1 << 15);
+    assert(tmp < (1 << 16));
+    vx = fmc.getOpampRev(tmp);
+
+    // Return vo.
+    return vx - (vc >> 14);
+}
+
+} // namespace reSIDfp

src/sound/resid-fp/Integrator8580.h

@@ -23,6 +23,7 @@
 #ifndef INTEGRATOR8580_H
 #define INTEGRATOR8580_H
 
+#include "Integrator.h"
 #include "FilterModelConfig8580.h"
 
 #include <stdint.h>
@@ -51,21 +52,16 @@ namespace reSIDfp
  *
  * Rfc gate voltage is generated by an OP Amp and depends on chip temperature.
  */
-class Integrator8580
+class Integrator8580 : public Integrator
 {
 private:
-    mutable int vx;
-    mutable int vc;
-
     unsigned short nVgt;
     unsigned short n_dac;
 
-    const FilterModelConfig8580* fmc;
+    FilterModelConfig8580& fmc;
 
 public:
-    Integrator8580(const FilterModelConfig8580* fmc) :
+    Integrator8580(FilterModelConfig8580& fmc) :
-        vx(0),
-        vc(0),
         fmc(fmc)
     {
         setV(1.5);
@@ -78,7 +74,7 @@ public:
     {
         // Normalized current factor, 1 cycle at 1MHz.
         // Fit in 5 bits.
-        n_dac = fmc->getNormalizedCurrentFactor(wl);
+        n_dac = fmc.getNormalizedCurrentFactor(wl);
     }
 
     /**
@@ -87,56 +83,19 @@ public:
     void setV(double v)
     {
         // Gate voltage is controlled by the switched capacitor voltage divider
-        // Ua = Ue * v = 4.76v  1<v<2
+        // Ua = Ue * v = 4.75v  1<v<2
         assert(v > 1.0 && v < 2.0);
-        const double Vg = 4.76 * v;
+        const double Vg = fmc.getVref() * v;
-        const double Vgt = Vg - fmc->getVth();
+        const double Vgt = Vg - fmc.getVth();
 
         // Vg - Vth, normalized so that translated values can be subtracted:
         // Vgt - x = (Vgt - t) - (x - t)
-        nVgt = fmc->getNormalizedValue(Vgt);
+        nVgt = fmc.getNormalizedValue(Vgt);
     }
 
-    int solve(int vi) const;
+    int solve(int vi) const override;
 };
 
 } // namespace reSIDfp
 
-#if RESID_INLINING || defined(INTEGRATOR8580_CPP)
-
-namespace reSIDfp
-{
-
-RESID_INLINE
-int Integrator8580::solve(int vi) const
-{
-    // Make sure we're not in subthreshold mode
-    assert(vx < nVgt);
-
-    // DAC voltages
-    const unsigned int Vgst = nVgt - vx;
-    const unsigned int Vgdt = (vi < nVgt) ? nVgt - vi : 0;  // triode/saturation mode
-
-    const unsigned int Vgst_2 = Vgst * Vgst;
-    const unsigned int Vgdt_2 = Vgdt * Vgdt;
-
-    // DAC current, scaled by (1/m)*2^13*m*2^16*m*2^16*2^-15 = m*2^30
-    const int n_I_dac = n_dac * (static_cast<int>(Vgst_2 - Vgdt_2) >> 15);
-
-    // Change in capacitor charge.
-    vc += n_I_dac;
-
-    // vx = g(vc)
-    const int tmp = (vc >> 15) + (1 << 15);
-    assert(tmp < (1 << 16));
-    vx = fmc->getOpampRev(tmp);
-
-    // Return vo.
-    return vx - (vc >> 14);
-}
-
-} // namespace reSIDfp
-
-#endif
-
 #endif

src/sound/resid-fp/OpAmp.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2015 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  *
  * This program is free software; you can redistribute it and/or modify
@@ -28,7 +28,7 @@
 namespace reSIDfp
 {
 
-const double EPSILON = 1e-8;
+constexpr double EPSILON = 1e-8;
 
 double OpAmp::solve(double n, double vi) const
 {
@@ -48,7 +48,7 @@ double OpAmp::solve(double n, double vi) const
 
         // Calculate f and df.
 
-        Spline::Point out = opamp->evaluate(x);
+        Spline::Point out = opamp.evaluate(x);
         const double vo = out.x;
         const double dvo = out.y;
 
@@ -64,9 +64,9 @@ double OpAmp::solve(double n, double vi) const
         // Newton-Raphson step: xk1 = xk - f(xk)/f'(xk)
         x -= f / df;
 
-        if (unlikely(fabs(x - xk) < EPSILON))
+        if (unlikely(std::fabs(x - xk) < EPSILON))
         {
-            out = opamp->evaluate(x);
+            out = opamp.evaluate(x);
             return out.x;
         }
 

src/sound/resid-fp/OpAmp.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2023 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004,2010 Dag Lem
  *
@@ -23,7 +23,6 @@
 #ifndef OPAMP_H
 #define OPAMP_H
 
-#include <memory>
 #include <vector>
 
 #include "Spline.h"
@@ -72,13 +71,13 @@ class OpAmp
 {
 private:
     /// Current root position (cached as guess to speed up next iteration)
-    mutable double x;
+    mutable double x = 0.;
 
     const double Vddt;
     const double vmin;
     const double vmax;
 
-    std::unique_ptr<Spline> const opamp;
+    Spline opamp;
 
 public:
     /**
@@ -89,14 +88,13 @@ public:
      * @param vmin
      * @param vmax
      */
-    OpAmp(const std::vector<Spline::Point> &opamp, double Vddt,
+    OpAmp(const std::vector<Spline::Point> &opamp_voltages, double Vddt,
             double vmin, double vmax
     ) :
-        x(0.),
         Vddt(Vddt),
         vmin(vmin),
         vmax(vmax),
-        opamp(new Spline(opamp)) {}
+        opamp(opamp_voltages) {}
 
     /**
      * Reset root position

src/sound/resid-fp/SID.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2016 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004 Dag Lem <resid@nimrod.no>
  *
@@ -26,11 +26,12 @@
 
 #include <limits>
 
+#include "sidcxx11.h"
+
 #include "array.h"
 #include "Dac.h"
 #include "Filter6581.h"
 #include "Filter8580.h"
-#include "Potentiometer.h"
 #include "WaveformCalculator.h"
 #include "resample/TwoPassSincResampler.h"
 #include "resample/ZeroOrderResampler.h"
@@ -38,8 +39,8 @@
 namespace reSIDfp
 {
 
-const unsigned int ENV_DAC_BITS = 8;
+constexpr unsigned int ENV_DAC_BITS = 8;
-const unsigned int OSC_DAC_BITS = 12;
+constexpr unsigned int OSC_DAC_BITS = 12;
 
 /**
  * The waveform D/A converter introduces a DC offset in the signal
@@ -106,8 +107,8 @@ const unsigned int OSC_DAC_BITS = 12;
  * On my 6581R4AR has 0x3A as the only value giving the same output level as 1.prg
  */
 //@{
-unsigned int constexpr OFFSET_6581 = 0x380;
+constexpr unsigned int OFFSET_6581 = 0x380;
-unsigned int constexpr OFFSET_8580 = 0x9c0;
+constexpr unsigned int OFFSET_8580 = 0x9c0;
 //@}
 
 /**
@@ -128,31 +129,24 @@ unsigned int constexpr OFFSET_8580 = 0x9c0;
  * [2]: http://noname.c64.org/csdb/forums/?roomid=11&topicid=29025&showallposts=1
  */
 //@{
-int constexpr BUS_TTL_6581 = 0x01d00;
+constexpr int BUS_TTL_6581 = 0x01d00;
-int constexpr BUS_TTL_8580 = 0xa2000;
+constexpr int BUS_TTL_8580 = 0xa2000;
 //@}
 
 SID::SID() :
     filter6581(new Filter6581()),
     filter8580(new Filter8580()),
-    externalFilter(new ExternalFilter()),
     resampler(nullptr),
-    potX(new Potentiometer()),
+    cws(AVERAGE)
-    potY(new Potentiometer())
 {
-    voice[0].reset(new Voice());
+    setChipModel(MOS6581);
-    voice[1].reset(new Voice());
-    voice[2].reset(new Voice());
-
-    muted[0] = muted[1] = muted[2] = false;
-
     reset();
-    setChipModel(MOS8580);
 }
 
 SID::~SID()
 {
-    // Needed to delete auto_ptr with complete type
+    delete filter6581;
+    delete filter8580;
 }
 
 void SID::setFilter6581Curve(double filterCurve)
@@ -160,6 +154,11 @@ void SID::setFilter6581Curve(double filterCurve)
     filter6581->setFilterCurve(filterCurve);
 }
 
+void SID::setFilter6581Range(double adjustment)
+{
+    filter6581->setFilterRange(adjustment);
+}
+
 void SID::setFilter8580Curve(double filterCurve)
 {
     filter8580->setFilterCurve(filterCurve);
@@ -178,7 +177,7 @@ void SID::voiceSync(bool sync)
         // Synchronize the 3 waveform generators.
         for (int i = 0; i < 3; i++)
         {
-            voice[i]->wave()->synchronize(voice[(i + 1) % 3]->wave(), voice[(i + 2) % 3]->wave());
+            voice[i].wave()->synchronize(voice[(i + 1) % 3].wave(), voice[(i + 2) % 3].wave());
         }
     }
 
@@ -187,10 +186,10 @@ void SID::voiceSync(bool sync)
 
     for (int i = 0; i < 3; i++)
     {
-        WaveformGenerator* const wave = voice[i]->wave();
+        WaveformGenerator* const wave = voice[i].wave();
         const unsigned int freq = wave->readFreq();
 
-        if (wave->readTest() || freq == 0 || !voice[(i + 1) % 3]->wave()->readSync())
+        if (wave->readTest() || freq == 0 || !voice[(i + 1) % 3].wave()->readSync())
         {
             continue;
         }
@@ -210,12 +209,14 @@ void SID::setChipModel(ChipModel model)
     switch (model)
     {
     case MOS6581:
-        filter = filter6581.get();
+        filter = filter6581;
+        scaleFactor = 3;
         modelTTL = BUS_TTL_6581;
         break;
 
     case MOS8580:
-        filter = filter8580.get();
+        filter = filter8580;
+        scaleFactor = 5;
         modelTTL = BUS_TTL_8580;
         break;
 
@@ -227,7 +228,7 @@ void SID::setChipModel(ChipModel model)
 
     // calculate waveform-related tables
     matrix_t* wavetables = WaveformCalculator::getInstance()->getWaveTable();
-    matrix_t* pulldowntables = WaveformCalculator::getInstance()->buildPulldownTable(model);
+    matrix_t* pulldowntables = WaveformCalculator::getInstance()->buildPulldownTable(model, cws);
 
     // calculate envelope DAC table
     {
@@ -247,7 +248,8 @@ void SID::setChipModel(ChipModel model)
         Dac dacBuilder(OSC_DAC_BITS);
         dacBuilder.kinkedDac(model);
 
-        const double offset = dacBuilder.getOutput(is6581 ? OFFSET_6581 : OFFSET_8580);
+        //const double offset = dacBuilder.getOutput(is6581 ? OFFSET_6581 : OFFSET_8580);
+        const double offset = dacBuilder.getOutput(0x7ff);
 
         for (unsigned int i = 0; i < (1 << OSC_DAC_BITS); i++)
         {
@@ -259,11 +261,35 @@ void SID::setChipModel(ChipModel model)
     // set voice tables
     for (int i = 0; i < 3; i++)
     {
-        voice[i]->setEnvDAC(envDAC);
+        voice[i].setEnvDAC(envDAC);
-        voice[i]->setWavDAC(oscDAC);
+        voice[i].setWavDAC(oscDAC);
-        voice[i]->wave()->setModel(is6581);
+        voice[i].wave()->setModel(is6581);
-        voice[i]->wave()->setWaveformModels(wavetables);
+        voice[i].wave()->setWaveformModels(wavetables);
-        voice[i]->wave()->setPulldownModels(pulldowntables);
+        voice[i].wave()->setPulldownModels(pulldowntables);
+    }
+}
+
+void SID::setCombinedWaveforms(CombinedWaveforms cws)
+{
+    switch (cws)
+    {
+    case AVERAGE:
+    case WEAK:
+    case STRONG:
+        break;
+
+    default:
+        throw SIDError("Unknown combined waveforms type");
+    }
+
+    this->cws = cws;
+
+    // rebuild waveform-related tables
+    matrix_t* pulldowntables = WaveformCalculator::getInstance()->buildPulldownTable(model, cws);
+
+    for (int i = 0; i < 3; i++)
+    {
+        voice[i].wave()->setPulldownModels(pulldowntables);
     }
 }
 
@@ -271,12 +297,12 @@ void SID::reset()
 {
     for (int i = 0; i < 3; i++)
     {
-        voice[i]->reset();
+        voice[i].reset();
     }
 
     filter6581->reset();
     filter8580->reset();
-    externalFilter->reset();
+    externalFilter.reset();
 
     if (resampler.get())
     {
@@ -299,22 +325,22 @@ unsigned char SID::read(int offset)
     switch (offset)
     {
     case 0x19: // X value of paddle
-        busValue = potX->readPOT();
+        busValue = potX.readPOT();
         busValueTtl = modelTTL;
         break;
 
     case 0x1a: // Y value of paddle
-        busValue = potY->readPOT();
+        busValue = potY.readPOT();
         busValueTtl = modelTTL;
         break;
 
     case 0x1b: // Voice #3 waveform output
-        busValue = voice[2]->wave()->readOSC();
+        busValue = voice[2].wave()->readOSC();
         busValueTtl = modelTTL;
         break;
 
     case 0x1c: // Voice #3 ADSR output
-        busValue = voice[2]->envelope()->readENV();
+        busValue = voice[2].envelope()->readENV();
         busValueTtl = modelTTL;
         break;
 
@@ -337,87 +363,87 @@ void SID::write(int offset, unsigned char value)
     switch (offset)
     {
     case 0x00: // Voice #1 frequency (Low-byte)
-        voice[0]->wave()->writeFREQ_LO(value);
+        voice[0].wave()->writeFREQ_LO(value);
         break;
 
     case 0x01: // Voice #1 frequency (High-byte)
-        voice[0]->wave()->writeFREQ_HI(value);
+        voice[0].wave()->writeFREQ_HI(value);
         break;
 
     case 0x02: // Voice #1 pulse width (Low-byte)
-        voice[0]->wave()->writePW_LO(value);
+        voice[0].wave()->writePW_LO(value);
         break;
 
     case 0x03: // Voice #1 pulse width (bits #8-#15)
-        voice[0]->wave()->writePW_HI(value);
+        voice[0].wave()->writePW_HI(value);
         break;
 
     case 0x04: // Voice #1 control register
-        voice[0]->writeCONTROL_REG(muted[0] ? 0 : value);
+        voice[0].writeCONTROL_REG(value);
         break;
 
     case 0x05: // Voice #1 Attack and Decay length
-        voice[0]->envelope()->writeATTACK_DECAY(value);
+        voice[0].envelope()->writeATTACK_DECAY(value);
         break;
 
     case 0x06: // Voice #1 Sustain volume and Release length
-        voice[0]->envelope()->writeSUSTAIN_RELEASE(value);
+        voice[0].envelope()->writeSUSTAIN_RELEASE(value);
         break;
 
     case 0x07: // Voice #2 frequency (Low-byte)
-        voice[1]->wave()->writeFREQ_LO(value);
+        voice[1].wave()->writeFREQ_LO(value);
         break;
 
     case 0x08: // Voice #2 frequency (High-byte)
-        voice[1]->wave()->writeFREQ_HI(value);
+        voice[1].wave()->writeFREQ_HI(value);
         break;
 
     case 0x09: // Voice #2 pulse width (Low-byte)
-        voice[1]->wave()->writePW_LO(value);
+        voice[1].wave()->writePW_LO(value);
         break;
 
     case 0x0a: // Voice #2 pulse width (bits #8-#15)
-        voice[1]->wave()->writePW_HI(value);
+        voice[1].wave()->writePW_HI(value);
         break;
 
     case 0x0b: // Voice #2 control register
-        voice[1]->writeCONTROL_REG(muted[1] ? 0 : value);
+        voice[1].writeCONTROL_REG(value);
         break;
 
     case 0x0c: // Voice #2 Attack and Decay length
-        voice[1]->envelope()->writeATTACK_DECAY(value);
+        voice[1].envelope()->writeATTACK_DECAY(value);
         break;
 
     case 0x0d: // Voice #2 Sustain volume and Release length
-        voice[1]->envelope()->writeSUSTAIN_RELEASE(value);
+        voice[1].envelope()->writeSUSTAIN_RELEASE(value);
         break;
 
     case 0x0e: // Voice #3 frequency (Low-byte)
-        voice[2]->wave()->writeFREQ_LO(value);
+        voice[2].wave()->writeFREQ_LO(value);
         break;
 
     case 0x0f: // Voice #3 frequency (High-byte)
-        voice[2]->wave()->writeFREQ_HI(value);
+        voice[2].wave()->writeFREQ_HI(value);
         break;
 
     case 0x10: // Voice #3 pulse width (Low-byte)
-        voice[2]->wave()->writePW_LO(value);
+        voice[2].wave()->writePW_LO(value);
         break;
 
     case 0x11: // Voice #3 pulse width (bits #8-#15)
-        voice[2]->wave()->writePW_HI(value);
+        voice[2].wave()->writePW_HI(value);
         break;
 
     case 0x12: // Voice #3 control register
-        voice[2]->writeCONTROL_REG(muted[2] ? 0 : value);
+        voice[2].writeCONTROL_REG(value);
         break;
 
     case 0x13: // Voice #3 Attack and Decay length
-        voice[2]->envelope()->writeATTACK_DECAY(value);
+        voice[2].envelope()->writeATTACK_DECAY(value);
         break;
 
     case 0x14: // Voice #3 Sustain volume and Release length
-        voice[2]->envelope()->writeSUSTAIN_RELEASE(value);
+        voice[2].envelope()->writeSUSTAIN_RELEASE(value);
         break;
 
     case 0x15: // Filter cut off frequency (bits #0-#2)
@@ -448,9 +474,9 @@ void SID::write(int offset, unsigned char value)
     voiceSync(false);
 }
 
-void SID::setSamplingParameters(double clockFrequency, SamplingMethod method, double samplingFrequency, double highestAccurateFrequency)
+void SID::setSamplingParameters(double clockFrequency, SamplingMethod method, double samplingFrequency)
 {
-    externalFilter->setClockFrequency(clockFrequency);
+    externalFilter.setClockFrequency(clockFrequency);
 
     switch (method)
     {
@@ -459,7 +485,7 @@ void SID::setSamplingParameters(double clockFrequency, SamplingMethod method, do
         break;
 
     case RESAMPLE:
-        resampler.reset(TwoPassSincResampler::create(clockFrequency, samplingFrequency, highestAccurateFrequency));
+        resampler.reset(TwoPassSincResampler::create(clockFrequency, samplingFrequency));
         break;
 
     default:
@@ -480,16 +506,16 @@ void SID::clockSilent(unsigned int cycles)
             for (int i = 0; i < delta_t; i++)
             {
                 // clock waveform generators (can affect OSC3)
-                voice[0]->wave()->clock();
+                voice[0].wave()->clock();
-                voice[1]->wave()->clock();
+                voice[1].wave()->clock();
-                voice[2]->wave()->clock();
+                voice[2].wave()->clock();
 
-                voice[0]->wave()->output(voice[2]->wave());
+                voice[0].wave()->output(voice[2].wave());
-                voice[1]->wave()->output(voice[0]->wave());
+                voice[1].wave()->output(voice[0].wave());
-                voice[2]->wave()->output(voice[1]->wave());
+                voice[2].wave()->output(voice[1].wave());
 
                 // clock ENV3 only
-                voice[2]->envelope()->clock();
+                voice[2].envelope()->clock();
             }
 
             cycles -= delta_t;

src/sound/resid-fp/Spline.cpp

@@ -92,11 +92,11 @@ Spline::Point Spline::evaluate(double x) const
 {
     if ((x < c->x1) || (x > c->x2))
     {
-        for (size_t i = 0; i < params.size(); i++)
+        for (const auto & param : params)
         {
-            if (x <= params[i].x2)
+            if (x <= param.x2)
             {
-                c = &params[i];
+                c = &param;
                 break;
             }
         }

src/sound/resid-fp/Spline.h

@@ -38,14 +38,14 @@ namespace reSIDfp
 class Spline
 {
 public:
-    typedef struct
+    using Point = struct
     {
         double x;
         double y;
-    } Point;
+    };
 
 private:
-    typedef struct
+    using Param = struct
     {
         double x1;
         double x2;
@@ -53,9 +53,9 @@ private:
         double b;
         double c;
         double d;
-    } Param;
+    };
 
-    typedef std::vector<Param> ParamVector;
+    using ParamVector = std::vector<Param>;
 
 private:
     /// Interpolation parameters

src/sound/resid-fp/Voice.h

@@ -23,14 +23,10 @@
 #ifndef VOICE_H
 #define VOICE_H
 
-#include <memory>
-
 #include "siddefs-fp.h"
 #include "WaveformGenerator.h"
 #include "EnvelopeGenerator.h"
 
-#include "sidcxx11.h"
-
 namespace reSIDfp
 {
 
@@ -40,9 +36,9 @@ namespace reSIDfp
 class Voice
 {
 private:
-    std::unique_ptr<WaveformGenerator> const waveformGenerator;
+    WaveformGenerator waveformGenerator;
 
-    std::unique_ptr<EnvelopeGenerator> const envelopeGenerator;
+    EnvelopeGenerator envelopeGenerator;
 
     /// The DAC LUT for analog waveform output
     float* wavDAC; //-V730_NOINIT this is initialized in the SID constructor
@@ -67,23 +63,16 @@ public:
      * @return the voice analog output
      */
     RESID_INLINE
-    int output(const WaveformGenerator* ringModulator) const
+    float output(const WaveformGenerator* ringModulator)
     {
-        unsigned int const wav = waveformGenerator->output(ringModulator);
+        unsigned int const wav = waveformGenerator.output(ringModulator);
-        unsigned int const env = envelopeGenerator->output();
+        unsigned int const env = envelopeGenerator.output();
 
         // DAC imperfections are emulated by using the digital output
         // as an index into a DAC lookup table.
-        return static_cast<int>(wavDAC[wav] * envDAC[env]);
+        return wavDAC[wav] * envDAC[env];
     }
 
-    /**
-     * Constructor.
-     */
-    Voice() :
-        waveformGenerator(new WaveformGenerator()),
-        envelopeGenerator(new EnvelopeGenerator()) {}
-
     /**
      * Set the analog DAC emulation for waveform generator.
      * Must be called before any operation.
@@ -100,9 +89,9 @@ public:
      */
     void setEnvDAC(float* dac) { envDAC = dac; }
 
-    WaveformGenerator* wave() const { return waveformGenerator.get(); }
+    WaveformGenerator* wave() { return &waveformGenerator; }
 
-    EnvelopeGenerator* envelope() const { return envelopeGenerator.get(); }
+    EnvelopeGenerator* envelope() { return &envelopeGenerator; }
 
     /**
      * Write control register.
@@ -111,8 +100,8 @@ public:
      */
     void writeCONTROL_REG(unsigned char control)
     {
-        waveformGenerator->writeCONTROL_REG(control);
+        waveformGenerator.writeCONTROL_REG(control);
-        envelopeGenerator->writeCONTROL_REG(control);
+        envelopeGenerator.writeCONTROL_REG(control);
     }
 
     /**
@@ -120,8 +109,8 @@ public:
      */
     void reset()
     {
-        waveformGenerator->reset();
+        waveformGenerator.reset();
-        envelopeGenerator->reset();
+        envelopeGenerator.reset();
     }
 };
 

src/sound/resid-fp/WaveformCalculator.cpp

@@ -21,66 +21,154 @@
 
 #include "WaveformCalculator.h"
 
+#include "sidcxx11.h"
+
+#include <map>
+#include <mutex>
 #include <cmath>
 
 namespace reSIDfp
 {
 
+/**
+ * Combined waveform model parameters.
+ */
+using distance_t = float (*)(float, int);
+
+using CombinedWaveformConfig = struct
+{
+    distance_t distFunc;
+    float threshold;
+    float topbit;
+    float pulsestrength;
+    float distance1;
+    float distance2;
+};
+
+using cw_cache_t = std::map<const CombinedWaveformConfig*, matrix_t>;
+
+cw_cache_t PULLDOWN_CACHE;
+
+std::mutex PULLDOWN_CACHE_Lock;
+
 WaveformCalculator* WaveformCalculator::getInstance()
 {
     static WaveformCalculator instance;
     return &instance;
 }
 
-/**
- * Parameters derived with the Monte Carlo method based on
- * samplings by kevtris. Code and data available in the project repository [1].
- *
- * The score here reported is the acoustic error
- * calculated XORing the estimated and the sampled values.
- * In parentheses the number of mispredicted bits.
- *
- * [1] https://github.com/libsidplayfp/combined-waveforms
- */
-const CombinedWaveformConfig config[2][5] =
-{
-    { /* kevtris chip G (6581 R2) */
-        {0.862147212f, 0.f,          10.8962431f,    2.50848103f }, // TS  error  1941 (327/28672)
-        {0.932746708f, 2.07508397f,   1.03668225f,   1.14876997f }, // PT  error  5992 (126/32768)
-        {0.860927045f, 2.43506575f,   0.908603609f,  1.07907593f }, // PS  error  3693 (521/28672)
-        {0.741343081f, 0.0452554375f, 1.1439606f,    1.05711341f }, // PTS error   338 ( 29/28672)
-        {0.96f,        2.5f,          1.1f,          1.2f        }, // NP  guessed
-    },
-    { /* kevtris chip V (8580 R5) */
-        {0.715788841f, 0.f,           1.32999945f,   2.2172699f  }, // TS  error   928 (135/32768)
-        {0.93500334f,  1.05977178f,   1.08629429f,   1.43518543f }, // PT  error  7991 (212/32768)
-        {0.920648575f, 0.943601072f,  1.13034654f,   1.41881108f }, // PS  error 12566 (394/32768)
-        {0.90921098f,  0.979807794f,  0.942194462f,  1.40958893f }, // PTS error  2092 ( 60/32768)
-        {0.95f,        1.15f,         1.f,           1.45f       }, // NP  guessed
-    },
-};
-
-typedef float (*distance_t)(float, int);
-
 // Distance functions
 static float exponentialDistance(float distance, int i)
 {
     return pow(distance, -i);
 }
 
-#if 0
 MAYBE_UNUSED static float linearDistance(float distance, int i)
 {
     return 1.f / (1.f + i * distance);
 }
-#endif
 
-#if 0
+static float quadraticDistance(float distance, int i)
-MAYBE_UNUSED static float quadraticDistance(float distance, int i)
 {
     return 1.f / (1.f + (i*i) * distance);
 }
-#endif
+
+/**
+ * Parameters derived with the Monte Carlo method based on
+ * samplings from real machines.
+ * Code and data available in the project repository [1].
+ * Sampling program made by Dag Lem [2].
+ *
+ * The score here reported is the acoustic error
+ * calculated XORing the estimated and the sampled values.
+ * In parentheses the number of mispredicted bits.
+ *
+ * [1] https://github.com/libsidplayfp/combined-waveforms
+ * [2] https://github.com/daglem/reDIP-SID/blob/master/research/combsample.d64
+ */
+const CombinedWaveformConfig configAverage[2][5] =
+{
+    { /* 6581 R3 0486S sampled by Trurl */
+        // TS  error  3555 (324/32768) [RMS: 73.98]
+        { exponentialDistance,  0.877322257f, 1.11349654f, 0.f, 2.14537621f, 9.08618164f },
+        // PT  error  4590 (124/32768) [RMS: 68.90]
+        { linearDistance, 0.941692829f, 1.f, 1.80072665f, 0.033124879f, 0.232303441f },
+        // PS  error 19352 (763/32768) [RMS: 96.91]
+        { linearDistance, 1.66494179f, 1.03760982f, 5.62705326f, 0.291590303f, 0.283631504f },
+        // PTS error  5068 ( 94/32768) [RMS: 41.69]
+        { linearDistance, 1.09762526f, 0.975265801f, 1.52196741f, 0.151528224f, 0.841949463f },
+        // NP  guessed
+        { exponentialDistance, 0.96f, 1.f, 2.5f, 1.1f, 1.2f },
+    },
+    { /* 8580 R5 1088 sampled by reFX-Mike */
+        // TS  error 10660 (353/32768) [RMS: 58.34]
+        { exponentialDistance, 0.853578329f, 1.09615636f, 0.f, 1.8819375f, 6.80794907f },
+        // PT  error 10635 (289/32768) [RMS: 108.81]
+        { exponentialDistance, 0.929835618f, 1.f, 1.12836814f, 1.10453653f, 1.48065746f },
+        // PS  error 12255 (554/32768) [RMS: 102.27]
+        { quadraticDistance, 0.911938608f, 0.996440411f, 1.2278074f, 0.000117214302f, 0.18948476f },
+        // PTS error  6913 (127/32768) [RMS: 55.80]
+        { exponentialDistance, 0.938004673f, 1.04827631f, 1.21178246f, 0.915959001f, 1.42698038f },
+        // NP  guessed
+        { exponentialDistance, 0.95f, 1.f, 1.15f, 1.f, 1.45f },
+    },
+};
+
+const CombinedWaveformConfig configWeak[2][5] =
+{
+    { /* 6581 R2 4383 sampled by ltx128 */
+        // TS  error 1474 (198/32768) [RMS: 62.81]
+        { exponentialDistance, 0.892563999f, 1.11905622f, 0.f, 2.21876144f, 9.63837719f },
+        // PT  error  612 (102/32768) [RMS: 43.71]
+        { linearDistance, 1.01262534f, 1.f, 2.46070528f, 0.0537485816f, 0.0986242667f },
+        // PS  error 8135 (575/32768) [RMS: 75.10]
+        { linearDistance, 2.14896345f, 1.0216713f, 10.5400085f, 0.244498149f, 0.126134038f },
+        // PTS error 2489 (60/32768) [RMS: 24.41]
+        { linearDistance, 1.22330308f, 0.933797896f, 2.83245254f, 0.0615176819f, 0.323831677f },
+        // NP  guessed
+        { exponentialDistance, 0.96f, 1.f, 2.5f, 1.1f, 1.2f },
+    },
+    { /* 8580 R5 4887 sampled by reFX-Mike */
+        // TS  error  741 (76/32768) [RMS: 53.74]
+        { exponentialDistance, 0.812351167f, 1.1727736f, 0.f, 1.87459648f, 2.31578159f },
+        // PT  error 7199 (192/32768) [RMS: 88.43]
+        { exponentialDistance, 0.917997837f, 1.f, 1.01248944f, 1.05761552f, 1.37529826f },
+        // PS  error 9856 (332/32768) [RMS: 86.29]
+        { quadraticDistance, 0.968754232f, 1.00669801f, 1.29909098f, 0.00962483883f, 0.146850556f },
+        // PTS error 4809 (60/32768) [RMS: 45.37]
+        { exponentialDistance, 0.941834152f, 1.06401193f, 0.991132736f, 0.995310068f, 1.41105855f },
+        // NP  guessed
+        { exponentialDistance, 0.95f, 1.f, 1.15f, 1.f, 1.45f },
+    },
+};
+
+const CombinedWaveformConfig configStrong[2][5] =
+{
+    { /* 6581 R2 0384 sampled by Trurl */
+        // TS  error 20337 (1579/32768) [RMS: 88.57]
+        { exponentialDistance, 0.000637792516f, 1.56725872f, 0.f, 0.00036806846f, 1.51800942f },
+        // PT  error  5190 (238/32768) [RMS: 83.54]
+        { linearDistance, 0.924780309f, 1.f, 1.96809769f, 0.0888123438f, 0.234606609f },
+        // PS  error 31015 (2181/32768) [RMS: 114.99]
+        { linearDistance, 1.2328074f, 0.73079139f, 3.9719491f, 0.00156516861f, 0.314677745f },
+        // PTS error  9874 (201/32768) [RMS: 52.30]
+        { linearDistance, 1.08558261f, 0.857638359f, 1.52781796f, 0.152927235f, 1.02657032f },
+        // NP  guessed
+        { exponentialDistance, 0.96f, 1.f, 2.5f, 1.1f, 1.2f },
+    },
+    { /* 8580 R5 1489 sampled by reFX-Mike */
+        // TS  error  4837 (388/32768) [RMS: 76.07]
+        { exponentialDistance, 0.89762634f, 56.7594185f, 0.f, 7.68995237f, 12.0754194f },
+        // PT  error  9266 (508/32768) [RMS: 127.83]
+        { exponentialDistance,  0.87147671f, 1.f, 1.44887495f, 1.05899632f, 1.43786001f },
+        // PS  error 13168 (718/32768) [RMS: 123.35]
+        { quadraticDistance, 0.89255774f, 1.2253896f, 1.75615835f, 0.0245045591f, 0.12982437f },
+        // PTS error  6702 (300/32768) [RMS: 71.01]
+        { linearDistance, 0.91124934f, 0.963609755f, 0.909965038f, 1.07445884f, 1.82399702f },
+        // NP  guessed
+        { exponentialDistance, 0.95f, 1.f, 1.15f, 1.f, 1.45f },
+    },
+};
 
 /// Calculate triangle waveform
 static unsigned int triXor(unsigned int val)
@@ -96,15 +184,17 @@ static unsigned int triXor(unsigned int val)
  * @param threshold
  * @param accumulator the high bits of the accumulator value
  */
-short calculatePulldown(float distancetable[], float pulsestrength, float threshold, unsigned int accumulator)
+short calculatePulldown(float distancetable[], float topbit, float pulsestrength, float threshold, unsigned int accumulator)
 {
-    unsigned char bit[12];
+    float bit[12];
 
     for (unsigned int i = 0; i < 12; i++)
     {
-        bit[i] = (accumulator & (1u << i)) != 0 ? 1 : 0;
+        bit[i] = (accumulator & (1u << i)) != 0 ? 1.f : 0.f;
     }
 
+    bit[11] *= topbit;
+
     float pulldown[12];
 
     for (int sb = 0; sb < 12; sb++)
@@ -117,7 +207,7 @@ short calculatePulldown(float distancetable[], float pulsestrength, float thresh
             if (cb == sb)
                 continue;
             const float weight = distancetable[sb - cb + 12];
-            avg += static_cast<float>(1 - bit[cb]) * weight;
+            avg += (1.f - bit[cb]) * weight;
             n += weight;
         }
 
@@ -131,7 +221,7 @@ short calculatePulldown(float distancetable[], float pulsestrength, float thresh
 
     for (unsigned int i = 0; i < 12; i++)
     {
-        const float bitValue = bit[i] != 0 ? 1.f - pulldown[i] : 0.f;
+        const float bitValue = bit[i] > 0.f ? 1.f - pulldown[i] : 0.f;
         if (bitValue > threshold)
         {
             value |= 1u << i;
@@ -157,9 +247,26 @@ WaveformCalculator::WaveformCalculator() :
     }
 }
 
-matrix_t* WaveformCalculator::buildPulldownTable(ChipModel model)
+matrix_t* WaveformCalculator::buildPulldownTable(ChipModel model, CombinedWaveforms cws)
 {
-    const CombinedWaveformConfig* cfgArray = config[model == MOS6581 ? 0 : 1];
+    std::lock_guard<std::mutex> lock(PULLDOWN_CACHE_Lock);
+
+    const int modelIdx = model == MOS6581 ? 0 : 1;
+    const CombinedWaveformConfig* cfgArray;
+
+    switch (cws)
+    {
+    default:
+    case AVERAGE:
+        cfgArray = configAverage[modelIdx];
+        break;
+    case WEAK:
+        cfgArray = configWeak[modelIdx];
+        break;
+    case STRONG:
+        cfgArray = configStrong[modelIdx];
+        break;
+    }
 
     cw_cache_t::iterator lb = PULLDOWN_CACHE.lower_bound(cfgArray);
 
@@ -174,7 +281,7 @@ matrix_t* WaveformCalculator::buildPulldownTable(ChipModel model)
     {
         const CombinedWaveformConfig& cfg = cfgArray[wav];
 
-        const distance_t distFunc = exponentialDistance;
+        const distance_t distFunc = cfg.distFunc;
 
         float distancetable[12 * 2 + 1];
         distancetable[12] = 1.f;
@@ -186,14 +293,11 @@ matrix_t* WaveformCalculator::buildPulldownTable(ChipModel model)
 
         for (unsigned int idx = 0; idx < (1u << 12); idx++)
         {
-            pdTable[wav][idx] = calculatePulldown(distancetable, cfg.pulsestrength, cfg.threshold, idx);
+            pdTable[wav][idx] = calculatePulldown(distancetable, cfg.topbit, cfg.pulsestrength, cfg.threshold, idx);
         }
     }
-#ifdef HAVE_CXX11
+
     return &(PULLDOWN_CACHE.emplace_hint(lb, cw_cache_t::value_type(cfgArray, pdTable))->second);
-#else
-    return &(PULLDOWN_CACHE.insert(lb, cw_cache_t::value_type(cfgArray, pdTable))->second);
-#endif
 }
 
 } // namespace reSIDfp

src/sound/resid-fp/WaveformCalculator.h

@@ -22,46 +22,33 @@
 #ifndef WAVEFORMCALCULATOR_h
 #define WAVEFORMCALCULATOR_h
 
-#include <map>
-
 #include "array.h"
-#include "sidcxx11.h"
+
 #include "siddefs-fp.h"
 
 
 namespace reSIDfp
 {
 
-/**
- * Combined waveform model parameters.
- */
-typedef struct
-{
-    float threshold;
-    float pulsestrength;
-    float distance1;
-    float distance2;
-} CombinedWaveformConfig;
-
 /**
  * Combined waveform calculator for WaveformGenerator.
  * By combining waveforms, the bits of each waveform are effectively short
- * circuited. A zero bit in one waveform will result in a zero output bit
+ * circuited, a zero bit in one waveform will result in a zero output bit,
- * (thus the infamous claim that the waveforms are AND'ed).
+ * thus the claim that the waveforms are AND'ed.
  * However, a zero bit in one waveform may also affect the neighboring bits
  * in the output.
  *
  * Example:
- * 
+ *
  *                 1 1
  *     Bit #       1 0 9 8 7 6 5 4 3 2 1 0
  *                 -----------------------
  *     Sawtooth    0 0 0 1 1 1 1 1 1 0 0 0
- *     
+ *
  *     Triangle    0 0 1 1 1 1 1 1 0 0 0 0
- *     
+ *
  *     AND         0 0 0 1 1 1 1 1 0 0 0 0
- *     
+ *
  *     Output      0 0 0 0 1 1 1 0 0 0 0 0
  *
  *
@@ -98,14 +85,9 @@ typedef struct
  */
 class WaveformCalculator
 {
-private:
-    typedef std::map<const CombinedWaveformConfig*, matrix_t> cw_cache_t;
-
 private:
     matrix_t wftable;
 
-    cw_cache_t PULLDOWN_CACHE;
-
 private:
     WaveformCalculator();
 
@@ -126,9 +108,10 @@ public:
      * Build pulldown table for use by WaveformGenerator.
      *
      * @param model Chip model to use
+     * @param cws strength of combined waveforms
      * @return Pulldown table
      */
-    matrix_t* buildPulldownTable(ChipModel model);
+    matrix_t* buildPulldownTable(ChipModel model, CombinedWaveforms cws);
 };
 
 } // namespace reSIDfp

src/sound/resid-fp/WaveformGenerator.cpp

@@ -40,13 +40,13 @@ namespace reSIDfp
  * and [VICE Bug #1128](http://sourceforge.net/p/vice-emu/bugs/1128/)
  */
 // ~95ms
-const unsigned int FLOATING_OUTPUT_TTL_6581R3  =   54000;
+constexpr unsigned int FLOATING_OUTPUT_TTL_6581R3  =   54000;
-const unsigned int FLOATING_OUTPUT_FADE_6581R3 =    1400;
+constexpr unsigned int FLOATING_OUTPUT_FADE_6581R3 =    1400;
 // ~1s
-//const unsigned int FLOATING_OUTPUT_TTL_6581R4  = 1000000;
+constexpr unsigned int FLOATING_OUTPUT_TTL_6581R4  = 1000000;
 // ~1s
-const unsigned int FLOATING_OUTPUT_TTL_8580R5  =  800000;
+constexpr unsigned int FLOATING_OUTPUT_TTL_8580R5  =  800000;
-const unsigned int FLOATING_OUTPUT_FADE_8580R5 =   50000;
+constexpr unsigned int FLOATING_OUTPUT_FADE_8580R5 =   50000;
 
 /**
  * Number of cycles after which the shift register is reset
@@ -58,15 +58,15 @@ const unsigned int FLOATING_OUTPUT_FADE_8580R5 =   50000;
  * only the big difference between the old and new models.
  */
 // ~210ms
-const unsigned int SHIFT_REGISTER_RESET_6581R3 =   50000;
+constexpr unsigned int SHIFT_REGISTER_RESET_6581R3 =   50000;
-const unsigned int SHIFT_REGISTER_FADE_6581R3  =   15000;
+constexpr unsigned int SHIFT_REGISTER_FADE_6581R3  =   15000;
 // ~2.15s
-//const unsigned int SHIFT_REGISTER_RESET_6581R4 = 2150000;
+constexpr unsigned int SHIFT_REGISTER_RESET_6581R4 = 2150000;
 // ~2.8s
-const unsigned int SHIFT_REGISTER_RESET_8580R5 =  986000;
+constexpr unsigned int SHIFT_REGISTER_RESET_8580R5 =  986000;
-const unsigned int SHIFT_REGISTER_FADE_8580R5  =  314300;
+constexpr unsigned int SHIFT_REGISTER_FADE_8580R5  =  314300;
 
-const unsigned int shift_mask =
+constexpr unsigned int shift_mask =
     ~(
         (1u <<  2) |  // Bit 20
         (1u <<  4) |  // Bit 18
@@ -107,15 +107,100 @@ const unsigned int shift_mask =
  * -----+-------+--------------+--------------
  * phi1 |   1   |   X --> X    |   A --> A      <- shift phase 2
  * phi2 |   1   |   X <-> X    |   A <-> A
+ *
+ *
+ * Normal cycles
+ * -------------
+ * Normally, when noise is selected along with another waveform,
+ * c1 and c2 are closed and the output bits pull down the corresponding
+ * shift register bits.
+ * 
+ *        noi_out_x             noi_out_x+1
+ *          ^                     ^
+ *          |                     |
+ *          +-------------+       +-------------+
+ *          |             |       |             |
+ *          +---o<|---+   |       +---o<|---+   |
+ *          |         |   |       |         |   |
+ *       c2 |      c1 |   |    c2 |      c1 |   |
+ *          |         |   |       |         |   |
+ *  >---/---+---|>o---+   +---/---+---|>o---+   +---/--->
+ *      LC                    LC                    LC
+ *
+ *
+ * Shift phase 1
+ * -------------
+ * During shift phase 1 c1 and c2 are open, the SR bits are floating
+ * and will be driven by the output of combined waveforms,
+ * or slowly turn high.
+ *
+ *        noi_out_x             noi_out_x+1
+ *          ^                     ^
+ *          |                     |
+ *          +-------------+       +-------------+
+ *          |             |       |             |
+ *          +---o<|---+   |       +---o<|---+   |
+ *          |         |   |       |         |   |
+ *       c2 /      c1 /   |    c2 /      c1 /   |
+ *          |         |   |       |         |   |
+ *  >-------+---|>o---+   +-------+---|>o---+   +------->
+ *      LC                    LC                    LC
+ *
+ *
+ * Shift phase 2 (phi1)
+ * --------------------
+ * During the first half cycle of shift phase 2 c1 is closed
+ * so the value from of noi_out_x-1 enters the bit.
+ *
+ *        noi_out_x             noi_out_x+1
+ *          ^                     ^
+ *          |                     |
+ *          +-------------+       +-------------+
+ *          |             |       |             |
+ *          +---o<|---+   |       +---o<|---+   |
+ *          |         |   |       |         |   |
+ *       c2 /      c1 |   |    c2 /      c1 |   |
+ *          |         |   |       |         |   |
+ *  >---/---+---|>o---+   +---/---+---|>o---+   +---/--->
+ *      LC                    LC                    LC
+ *
+ *
+ * Shift phase 2 (phi2)
+ * --------------------
+ * On the second half of shift phase 2 c2 closes and
+ * we're back to normal cycles.
  */
 
 inline bool do_writeback(unsigned int waveform_old, unsigned int waveform_new, bool is6581)
 {
     // no writeback without combined waveforms
-    if (waveform_new <= 8)
+
-        return false;
     if (waveform_old <= 8)
-        return false; // fixes SID/noisewriteback/noise_writeback_test2-{old,new}
+        // fixes SID/noisewriteback/noise_writeback_test2-{old,new}
+        return false;
+
+    if (waveform_new < 8)
+        return false;
+
+    if ((waveform_new == 8)
+        // breaks noise_writeback_check_F_to_8_old
+        // but fixes simple and scan
+        && (waveform_old != 0xf))
+    {
+        // fixes
+        // noise_writeback_check_9_to_8_old
+        // noise_writeback_check_A_to_8_old
+        // noise_writeback_check_B_to_8_old
+        // noise_writeback_check_D_to_8_old
+        // noise_writeback_check_E_to_8_old
+        // noise_writeback_check_F_to_8_old
+        // noise_writeback_check_9_to_8_new
+        // noise_writeback_check_A_to_8_new
+        // noise_writeback_check_D_to_8_new
+        // noise_writeback_check_E_to_8_new
+        // noise_writeback_test1-{old,new}
+        return false;
+    }
 
     // What's happening here?
     if (is6581 &&
@@ -190,8 +275,16 @@ void WaveformGenerator::write_shift_register()
 {
     if (unlikely(waveform > 0x8))
     {
+#if 0
+        // FIXME this breaks SID/wf12nsr/wf12nsr
         if (waveform == 0xc)
-            return; // breaks SID/wf12nsr/wf12nsr
+            // fixes
+            // noise_writeback_check_8_to_C_old
+            // noise_writeback_check_9_to_C_old
+            // noise_writeback_check_A_to_C_old
+            // noise_writeback_check_C_to_C_old
+            return;
+#endif
 
         // Write changes to the shift register output caused by combined waveforms
         // back into the shift register.

src/sound/resid-fp/WaveformGenerator.h

@@ -93,64 +93,64 @@ namespace reSIDfp
 class WaveformGenerator
 {
 private:
-    matrix_t* model_wave;
+    matrix_t* model_wave = nullptr;
-    matrix_t* model_pulldown;
+    matrix_t* model_pulldown = nullptr;
 
-    short* wave;
+    short* wave = nullptr;
-    short* pulldown;
+    short* pulldown = nullptr;
 
     // PWout = (PWn/40.95)%
-    unsigned int pw;
+    unsigned int pw = 0;
 
-    unsigned int shift_register;
+    unsigned int shift_register = 0;
 
     /// Shift register is latched when transitioning to shift phase 1.
-    unsigned int shift_latch;
+    unsigned int shift_latch = 0;
 
     /// Emulation of pipeline causing bit 19 to clock the shift register.
-    int shift_pipeline;
+    int shift_pipeline = 0;
 
-    unsigned int ring_msb_mask;
+    unsigned int ring_msb_mask = 0;
-    unsigned int no_noise;
+    unsigned int no_noise = 0;
-    unsigned int noise_output;
+    unsigned int noise_output = 0;
-    unsigned int no_noise_or_noise_output;
+    unsigned int no_noise_or_noise_output = 0;
-    unsigned int no_pulse;
+    unsigned int no_pulse = 0;
-    unsigned int pulse_output;
+    unsigned int pulse_output = 0;
 
     /// The control register right-shifted 4 bits; used for output function table lookup.
-    unsigned int waveform;
+    unsigned int waveform = 0;
 
-    unsigned int waveform_output;
+    unsigned int waveform_output = 0;
 
     /// Current accumulator value.
-    unsigned int accumulator;
+    unsigned int accumulator = 0x555555; // Accumulator's even bits are high on powerup
 
     // Fout = (Fn*Fclk/16777216)Hz
-    unsigned int freq;
+    unsigned int freq = 0;
 
     /// 8580 tri/saw pipeline
-    unsigned int tri_saw_pipeline;
+    unsigned int tri_saw_pipeline = 0x555;
 
     /// The OSC3 value
-    unsigned int osc3;
+    unsigned int osc3 = 0;
 
     /// Remaining time to fully reset shift register.
-    unsigned int shift_register_reset;
+    unsigned int shift_register_reset = 0;
 
     // The wave signal TTL when no waveform is selected.
-    unsigned int floating_output_ttl;
+    unsigned int floating_output_ttl = 0;
 
     /// The control register bits. Gate is handled by EnvelopeGenerator.
     //@{
-    bool test;
+    bool test = false;
-    bool sync;
+    bool sync = false;
     //@}
 
     /// Test bit is latched at phi2 for the noise XOR.
     bool test_or_reset;
 
     /// Tell whether the accumulator MSB was set high on this cycle.
-    bool msb_rising;
+    bool msb_rising = false;
 
     bool is6581; //-V730_NOINIT this is initialized in the SID constructor
 
@@ -160,7 +160,7 @@ private:
     void write_shift_register();
 
     void set_noise_output();
-    
+
     void set_no_noise_or_noise_output();
 
     void waveBitfade();
@@ -194,35 +194,6 @@ public:
      */
     void synchronize(WaveformGenerator* syncDest, const WaveformGenerator* syncSource) const;
 
-    /**
-     * Constructor.
-     */
-    WaveformGenerator() :
-        model_wave(nullptr),
-        model_pulldown(nullptr),
-        wave(nullptr),
-        pulldown(nullptr),
-        pw(0),
-        shift_register(0),
-        shift_pipeline(0),
-        ring_msb_mask(0),
-        no_noise(0),
-        noise_output(0),
-        no_noise_or_noise_output(0),
-        no_pulse(0),
-        pulse_output(0),
-        waveform(0),
-        waveform_output(0),
-        accumulator(0x555555),          // Accumulator's even bits are high on powerup
-        freq(0),
-        tri_saw_pipeline(0x555),
-        osc3(0),
-        shift_register_reset(0),
-        floating_output_ttl(0),
-        test(false),
-        sync(false),
-        msb_rising(false) {}
-
     /**
      * Write FREQ LO register.
      *
@@ -397,13 +368,13 @@ unsigned int WaveformGenerator::output(const WaveformGenerator* ringModulator)
         {
             osc3 = waveform_output;
         }
-
         // In the 6581 the top bit of the accumulator may be driven low by combined waveforms
         // when the sawtooth is selected
-        if (is6581
+        if (is6581 && (waveform & 0x2) && ((waveform_output & 0x800) == 0))
-                && (waveform & 0x2)
+        {
-                && ((waveform_output & 0x800) == 0))
+            msb_rising = 0;
             accumulator &= 0x7fffff;
+        }
 
         write_shift_register();
     }

src/sound/resid-fp/array.h

@@ -26,9 +26,7 @@
 #  include "config.h"
 #endif
 
-#ifdef HAVE_CXX11
+#include <atomic>
-#  include <atomic>
-#endif
 
 /**
  * Counter.
@@ -36,11 +34,7 @@
 class counter
 {
 private:
-#ifndef HAVE_CXX11
-    volatile unsigned int c;
-#else
     std::atomic<unsigned int> c;
-#endif
 
 public:
     counter() : c(1) {}
@@ -81,6 +75,6 @@ public:
     T const* operator[](unsigned int a) const { return &data[a * y]; }
 };
 
-typedef matrix<short> matrix_t;
+using matrix_t = matrix<short>;
 
 #endif

src/sound/resid-fp/config.h

@@ -1 +1 @@
-#define HAVE_CXX14
+#define HAVE_CXX17

src/sound/resid-fp/resample/Resampler.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2020 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  *
  * This program is free software; you can redistribute it and/or modify
@@ -23,6 +23,7 @@
 #define RESAMPLER_H
 
 #include <cmath>
+#include <cassert>
 
 #include "../sidcxx11.h"
 
@@ -37,28 +38,45 @@ namespace reSIDfp
  */
 class Resampler
 {
-protected:
+private:
-    inline short softClip(int x) const
+    template<int m>
+    static inline int clipper(int x)
     {
+        assert(x >= 0);
         constexpr int threshold = 28000;
         if (likely(x < threshold))
             return x;
 
-        constexpr double t = threshold / 32768.;
+        constexpr double max_val = static_cast<double>(m);
+        constexpr double t = threshold / max_val;
         constexpr double a = 1. - t;
         constexpr double b = 1. / a;
 
-        double value = static_cast<double>(x - threshold) / 32768.;
+        double value = static_cast<double>(x - threshold) / max_val;
-        value = t + a * tanh(b * value);
+        value = t + a * std::tanh(b * value);
-        return static_cast<short>(value * 32768.);
+        return static_cast<int>(value * max_val);
     }
 
+    /*
+     * Soft Clipping implementation, splitted for test.
+     */
+    static inline int softClipImpl(int x)
+    {
+        return x < 0 ? -clipper<32768>(-x) : clipper<32767>(x);
+    }
+
+protected:
+    /*
+     * Soft Clipping into 16 bit range [-32768,32767]
+     */
+    static inline short softClip(int x) { return static_cast<short>(softClipImpl(x)); }
+
     virtual int output() const = 0;
 
     Resampler() {}
 
 public:
-    virtual ~Resampler() {}
+    virtual ~Resampler() = default;
 
     /**
      * Input a sample into resampler. Output "true" when resampler is ready with new sample.
@@ -73,9 +91,10 @@ public:
      *
      * @return resampled sample
      */
-    short getOutput() const
+    inline short getOutput(int scaleFactor) const
     {
-        return softClip(output());
+        const int out = (scaleFactor * output()) / 2;
+        return softClip(out);
     }
 
     virtual void reset() = 0;

src/sound/resid-fp/resample/SincResampler.cpp

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2020 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004 Dag Lem <resid@nimrod.no>
  *
@@ -22,11 +22,16 @@
 
 #include "SincResampler.h"
 
+#ifdef HAVE_CXX20
+#  include <numbers>
+#endif
+
+#include <algorithm>
+#include <iterator>
 #include <cassert>
 #include <cstring>
 #include <cmath>
-#include <iostream>
+#include <cstdint>
-#include <sstream>
 
 #include "../siddefs-fp.h"
 
@@ -34,10 +39,8 @@
 #  include "config.h"
 #endif
 
-#ifdef HAVE_EMMINTRIN_H
+#ifdef HAVE_SMMINTRIN_H
-#  include <emmintrin.h>
+#  include <smmintrin.h>
-#elif defined HAVE_MMINTRIN_H
-#  include <mmintrin.h>
 #elif defined(HAVE_ARM_NEON_H)
 #  include <arm_neon.h>
 #endif
@@ -45,15 +48,10 @@
 namespace reSIDfp
 {
 
-typedef std::map<std::string, matrix_t> fir_cache_t;
-
-/// Cache for the expensive FIR table computation results.
-fir_cache_t FIR_CACHE;
-
 /// Maximum error acceptable in I0 is 1e-6, or ~96 dB.
-const double I0E = 1e-6;
+constexpr double I0E = 1e-6;
 
-const int BITS = 16;
+constexpr int BITS = 16;
 
 /**
  * Compute the 0th order modified Bessel function of the first kind.
@@ -90,7 +88,7 @@ double I0(double x)
  * @param bLength length of the sinc buffer
  * @return convolved result
  */
-int convolve(const short* a, const short* b, int bLength)
+int convolve(const int* a, const short* b, int bLength)
 {
 #ifdef HAVE_EMMINTRIN_H
     int out = 0;
@@ -102,7 +100,7 @@ int convolve(const short* a, const short* b, int bLength)
     {
         if (offset)
         {
-            const int l = (0x10 - offset)/2;
+            const int l = (0x10 - offset) / 2;
 
             for (int i = 0; i < l; i++)
             {
@@ -208,9 +206,9 @@ int convolve(const short* a, const short* b, int bLength)
     bLength &= 3;
 #else
     int32x4_t acc = vdupq_n_s32(0);
-    
+
     const int n = bLength / 4;
-    
+
     for (int i = 0; i < n; i++)
     {
         const int16x4_t h_vec = vld1_s16(a);
@@ -219,12 +217,12 @@ int convolve(const short* a, const short* b, int bLength)
         a += 4;
         b += 4;
     }
-    
+
     int out = vgetq_lane_s32(acc, 0) +
               vgetq_lane_s32(acc, 1) +
               vgetq_lane_s32(acc, 2) +
               vgetq_lane_s32(acc, 3);
-    
+
     bLength &= 3;
 #endif
 #else
@@ -233,7 +231,7 @@ int convolve(const short* a, const short* b, int bLength)
 
     for (int i = 0; i < bLength; i++)
     {
-        out += *a++ * *b++;
+        out += a[i] * static_cast<int>(b[i]);
     }
 
     return (out + (1 << 14)) >> 15;
@@ -265,17 +263,27 @@ int SincResampler::fir(int subcycle)
     return v1 + (firTableOffset * (v2 - v1) >> 10);
 }
 
-SincResampler::SincResampler(double clockFrequency, double samplingFrequency, double highestAccurateFrequency) :
+SincResampler::SincResampler(
-    sampleIndex(0),
+        double clockFrequency,
-    cyclesPerSample(static_cast<int>(clockFrequency / samplingFrequency * 1024.)),
+        double samplingFrequency,
-    sampleOffset(0),
+        double highestAccurateFrequency) :
-    outputValue(0)
+    cyclesPerSample(static_cast<int>(clockFrequency / samplingFrequency * 1024.))
 {
+#if defined(HAVE_CXX20) && defined(__cpp_lib_constexpr_cmath)
+    constexpr double PI = std::numbers::pi;
+#else
+#  ifdef M_PI
+        constexpr double PI = M_PI;
+#else
+        constexpr double PI = 3.14159265358979323846;
+#  endif
+#endif
+
     // 16 bits -> -96dB stopband attenuation.
-    const double A = -20. * log10(1.0 / (1 << BITS));
+    const double A = -20. * std::log10(1.0 / (1 << BITS));
     // A fraction of the bandwidth is allocated to the transition band, which we double
     // because we design the filter to transition halfway at nyquist.
-    const double dw = (1. - 2.*highestAccurateFrequency / samplingFrequency) * M_PI * 2.;
+    const double dw = (1. - 2.*highestAccurateFrequency / samplingFrequency) * PI * 2.;
 
     // For calculation of beta and N see the reference for the kaiserord
     // function in the MATLAB Signal Processing Toolbox:
@@ -283,6 +291,7 @@ SincResampler::SincResampler(double clockFrequency, double samplingFrequency, do
     const double beta = 0.1102 * (A - 8.7);
     const double I0beta = I0(beta);
     const double cyclesPerSampleD = clockFrequency / samplingFrequency;
+    const double inv_cyclesPerSampleD = samplingFrequency / clockFrequency;
 
     {
         // The filter order will maximally be 124 with the current constraints.
@@ -302,40 +311,22 @@ SincResampler::SincResampler(double clockFrequency, double samplingFrequency, do
         assert(firN < RINGSIZE);
 
         // Error is bounded by err < 1.234 / L^2, so L = sqrt(1.234 / (2^-16)) = sqrt(1.234 * 2^16).
-        firRES = static_cast<int>(ceil(sqrt(1.234 * (1 << BITS)) / cyclesPerSampleD));
+        firRES = static_cast<int>(std::ceil(std::sqrt(1.234 * (1 << BITS)) * inv_cyclesPerSampleD));
 
         // firN*firRES represent the total resolution of the sinc sampling. JOS
         // recommends a length of 2^BITS, but we don't quite use that good a filter.
         // The filter test program indicates that the filter performs well, though.
     }
 
-    // Create the map key
-    std::ostringstream o;
-    o << firN << "," << firRES << "," << cyclesPerSampleD;
-    const std::string firKey = o.str();
-    fir_cache_t::iterator lb = FIR_CACHE.lower_bound(firKey);
-
-    // The FIR computation is expensive and we set sampling parameters often, but
-    // from a very small set of choices. Thus, caching is used to speed initialization.
-    if (lb != FIR_CACHE.end() && !(FIR_CACHE.key_comp()(firKey, lb->first)))
-    {
-        firTable = &(lb->second);
-    }
-    else
     {
         // Allocate memory for FIR tables.
-        matrix_t tempTable(firRES, firN);
+        firTable = new matrix_t(firRES, firN);
-#ifdef HAVE_CXX11
-        firTable = &(FIR_CACHE.emplace_hint(lb, fir_cache_t::value_type(firKey, tempTable))->second);
-#else
-        firTable = &(FIR_CACHE.insert(lb, fir_cache_t::value_type(firKey, tempTable))->second);
-#endif
 
         // The cutoff frequency is midway through the transition band, in effect the same as nyquist.
-        const double wc = M_PI;
+        const double wc = PI;
 
         // Calculate the sinc tables.
-        const double scale = 32768.0 * wc / cyclesPerSampleD / M_PI;
+        const double scale = 32768.0 * wc * inv_cyclesPerSampleD / PI;
 
         // we're not interested in the fractional part
         // so use int division before converting to double
@@ -351,10 +342,10 @@ SincResampler::SincResampler(double clockFrequency, double samplingFrequency, do
                 const double x = j - jPhase;
 
                 const double xt = x / firN_2;
-                const double kaiserXt = fabs(xt) < 1. ? I0(beta * sqrt(1. - xt * xt)) / I0beta : 0.;
+                const double kaiserXt = std::fabs(xt) < 1. ? I0(beta * std::sqrt(1. - xt * xt)) / I0beta : 0.;
 
-                const double wt = wc * x / cyclesPerSampleD;
+                const double wt = wc * x * inv_cyclesPerSampleD;
-                const double sincWt = fabs(wt) >= 1e-8 ? sin(wt) / wt : 1.;
+                const double sincWt = std::fabs(wt) >= 1e-8 ? std::sin(wt) / wt : 1.;
 
                 (*firTable)[i][j] = static_cast<short>(scale * sincWt * kaiserXt);
             }
@@ -362,18 +353,16 @@ SincResampler::SincResampler(double clockFrequency, double samplingFrequency, do
     }
 }
 
+SincResampler::~SincResampler()
+{
+    delete firTable;
+}
+
 bool SincResampler::input(int input)
 {
     bool ready = false;
 
-    /*
+    sample[sampleIndex] = sample[sampleIndex + RINGSIZE] = input;
-     * Clip the input as it may overflow the 16 bit range.
-     *
-     * Approximate measured input ranges:
-     * 6581: [-24262,+25080]  (Kawasaki_Synthesizer_Demo)
-     * 8580: [-21514,+35232]  (64_Forever, Drum_Fool)
-     */
-    sample[sampleIndex] = sample[sampleIndex + RINGSIZE] = softClip(input);
     sampleIndex = (sampleIndex + 1) & (RINGSIZE - 1);
 
     if (sampleOffset < 1024)
@@ -390,7 +379,7 @@ bool SincResampler::input(int input)
 
 void SincResampler::reset()
 {
-    memset(sample, 0, sizeof(sample));
+    std::fill(std::begin(sample), std::end(sample), 0);
     sampleOffset = 0;
 }
 

src/sound/resid-fp/resample/SincResampler.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2013 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004 Dag Lem <resid@nimrod.no>
  *
@@ -25,13 +25,8 @@
 
 #include "Resampler.h"
 
-#include <string>
-#include <map>
-
 #include "../array.h"
 
-#include "../sidcxx11.h"
-
 namespace reSIDfp
 {
 
@@ -54,13 +49,13 @@ class SincResampler final : public Resampler
 {
 private:
     /// Size of the ring buffer, must be a power of 2
-    static const int RINGSIZE = 2048;
+    static constexpr int RINGSIZE = 2048;
 
 private:
     /// Table of the fir filter coefficients
     matrix_t* firTable;
 
-    int sampleIndex;
+    int sampleIndex = 0;
 
     /// Filter resolution
     int firRES;
@@ -70,11 +65,11 @@ private:
 
     const int cyclesPerSample;
 
-    int sampleOffset;
+    int sampleOffset = 0;
 
-    int outputValue;
+    int outputValue = 0;
 
-    short sample[RINGSIZE * 2];
+    int sample[RINGSIZE * 2];
 
 private:
     int fir(int subcycle);
@@ -82,25 +77,25 @@ private:
 public:
     /**
      * Use a clock freqency of 985248Hz for PAL C64, 1022730Hz for NTSC C64.
-     * The default end of passband frequency is pass_freq = 0.9*sample_freq/2
-     * for sample frequencies up to ~ 44.1kHz, and 20kHz for higher sample frequencies.
      *
-     * For resampling, the ratio between the clock frequency and the sample frequency
+     * For resampling, the ratio between the clock frequency
-     * is limited as follows: 125*clock_freq/sample_freq < 16384
+     * and the sample frequency is limited as follows:
+     *   125*clock_freq/sample_freq < 16384
+     *
      * E.g. provided a clock frequency of ~ 1MHz, the sample frequency
      * can not be set lower than ~ 8kHz.
-     * A lower sample frequency would make the resampling code overfill its 16k sample ring buffer.
+     * A lower sample frequency would make the resampling code overfill
-     *
+     * its 16k sample ring buffer.
-     * The end of passband frequency is also limited: pass_freq <= 0.9*sample_freq/2
-     *
-     * E.g. for a 44.1kHz sampling rate the end of passband frequency is limited
-     * to slightly below 20kHz. This constraint ensures that the FIR table is not overfilled.
      *
      * @param clockFrequency System clock frequency at Hz
      * @param samplingFrequency Desired output sampling rate
-     * @param highestAccurateFrequency
+     * @param highestAccurateFrequency passband frequency limit
      */
-    SincResampler(double clockFrequency, double samplingFrequency, double highestAccurateFrequency);
+    SincResampler(
+        double clockFrequency,
+        double samplingFrequency,
+        double highestAccurateFrequency);
+    ~SincResampler() override;
 
     bool input(int input) override;
 

src/sound/resid-fp/resample/TwoPassSincResampler.h

@@ -51,14 +51,25 @@ private:
 
 public:
     // Named constructor
-    static TwoPassSincResampler* create(double clockFrequency, double samplingFrequency, double highestAccurateFrequency)
+    static TwoPassSincResampler* create(double clockFrequency, double samplingFrequency)
     {
-        // Calculation according to Laurent Ganier. It evaluates to about 120 kHz at typical settings.
+        // Set the passband frequency slightly below half sampling frequency
+        //   pass_freq <= 0.9*sample_freq/2
+        //
+        // This constraint ensures that the FIR table is not overfilled.
+        // For higher sampling frequencies we're fine with 20KHz
+        const double halfFreq = (samplingFrequency > 44000.)
+            ? 20000. : samplingFrequency * 0.45;
+
+        // Calculation according to Laurent Ganier.
+        // It evaluates to about 120 kHz at typical settings.
         // Some testing around the chosen value seems to confirm that this does work.
-        double const intermediateFrequency = 2. * highestAccurateFrequency
+        double const intermediateFrequency = 2. * halfFreq
-            + sqrt(2. * highestAccurateFrequency * clockFrequency
+            + std::sqrt(2. * halfFreq * clockFrequency
-                * (samplingFrequency - 2. * highestAccurateFrequency) / samplingFrequency);
+                * (samplingFrequency - 2. * halfFreq) / samplingFrequency);
-        return new TwoPassSincResampler(clockFrequency, samplingFrequency, highestAccurateFrequency, intermediateFrequency);
+
+        return new TwoPassSincResampler(
+            clockFrequency, samplingFrequency, halfFreq, intermediateFrequency);
     }
 
     bool input(int sample) override

src/sound/resid-fp/resample/test.cpp

@@ -35,6 +35,10 @@
 #  define unique_ptr auto_ptr
 #endif
 
+#ifndef M_PI
+#  define M_PI    3.14159265358979323846
+#endif
+
 /**
  * Simple sin waveform in, power output measurement function.
  * It would be far better to use FFT.
@@ -57,7 +61,7 @@ int main(int, const char*[])
 
         for (int j = 0; j < RINGSIZE; j ++)
         {
-            int signal = static_cast<int>(32768.0 * sin(k++ * omega) * sqrt(2));
+            int signal = static_cast<int>(32768.0 * std::sin(k++ * omega) * sqrt(2));
             r->input(signal);
         }
 
@@ -67,7 +71,7 @@ int main(int, const char*[])
         /* Now, during measurement stage, put 100 cycles of waveform through filter. */
         for (int j = 0; j < 100000; j ++)
         {
-            int signal = static_cast<int>(32768.0 * sin(k++ * omega) * sqrt(2));
+            int signal = static_cast<int>(32768.0 * std::sin(k++ * omega) * std::sqrt(2));
 
             if (r->input(signal))
             {
@@ -77,7 +81,7 @@ int main(int, const char*[])
             }
         }
 
-        results.insert(std::make_pair(freq, 10 * log10(pwr / n)));
+        results.insert(std::make_pair(freq, 10 * std::log10(pwr / n)));
     }
 
     clock_t end = clock();

src/sound/resid-fp/sid.h

@@ -1,7 +1,7 @@
 /*
  * This file is part of libsidplayfp, a SID player engine.
  *
- * Copyright 2011-2016 Leandro Nini <drfiemost@users.sourceforge.net>
+ * Copyright 2011-2024 Leandro Nini <drfiemost@users.sourceforge.net>
  * Copyright 2007-2010 Antti Lankila
  * Copyright 2004 Dag Lem <resid@nimrod.no>
  *
@@ -26,6 +26,9 @@
 #include <memory>
 
 #include "siddefs-fp.h"
+#include "ExternalFilter.h"
+#include "Potentiometer.h"
+#include "Voice.h"
 
 #include "sidcxx11.h"
 
@@ -35,9 +38,6 @@ namespace reSIDfp
 class Filter;
 class Filter6581;
 class Filter8580;
-class ExternalFilter;
-class Potentiometer;
-class Voice;
 class Resampler;
 
 /**
@@ -64,28 +64,31 @@ private:
     Filter* filter;
 
     /// Filter used, if model is set to 6581
-    std::unique_ptr<Filter6581> const filter6581;
+    Filter6581* const filter6581;
 
     /// Filter used, if model is set to 8580
-    std::unique_ptr<Filter8580> const filter8580;
+    Filter8580* const filter8580;
+
+    /// Resampler used by audio generation code.
+    std::unique_ptr<Resampler> resampler;
 
     /**
      * External filter that provides high-pass and low-pass filtering
      * to adjust sound tone slightly.
      */
-    std::unique_ptr<ExternalFilter> const externalFilter;
+    ExternalFilter externalFilter;
-
-    /// Resampler used by audio generation code.
-    std::unique_ptr<Resampler> resampler;
 
     /// Paddle X register support
-    std::unique_ptr<Potentiometer> const potX;
+    Potentiometer potX;
 
     /// Paddle Y register support
-    std::unique_ptr<Potentiometer> const potY;
+    Potentiometer potY;
 
     /// SID voices
-    std::unique_ptr<Voice> voice[3];
+    Voice voice[3];
+
+    /// Used to amplify the output by x/2 to get an adequate playback volume
+    int scaleFactor;
 
     /// Time to live for the last written value
     int busValueTtl;
@@ -99,12 +102,12 @@ private:
     /// Currently active chip model.
     ChipModel model;
 
+    /// Currently selected combined waveforms strength.
+    CombinedWaveforms cws;
+
     /// Last written value
     unsigned char busValue;
 
-    /// Flags for muted channels
-    bool muted[3];
-
     /**
      * Emulated nonlinearity of the envelope DAC.
      *
@@ -132,7 +135,7 @@ private:
      *
      * @return the output sample
      */
-    int output() const;
+    int output();
 
     /**
      * Calculate the numebr of cycles according to current parameters
@@ -159,6 +162,14 @@ public:
      */
     ChipModel getChipModel() const { return model; }
 
+    /**
+     * Set combined waveforms strength.
+     *
+     * @param cws strength of combined waveforms
+     * @throw SIDError
+     */
+    void setCombinedWaveforms(CombinedWaveforms cws);
+
     /**
      * SID reset.
      */
@@ -204,14 +215,6 @@ public:
      */
     void write(int offset, unsigned char value);
 
-    /**
-     * SID voice muting.
-     *
-     * @param channel channel to modify
-     * @param enable is muted?
-     */
-    void mute(int channel, bool enable) { muted[channel] = enable; }
-
     /**
      * Setting of SID sampling parameters.
      *
@@ -237,7 +240,11 @@ public:
      * @param highestAccurateFrequency
      * @throw SIDError
      */
-    void setSamplingParameters(double clockFrequency, SamplingMethod method, double samplingFrequency, double highestAccurateFrequency);
+    void setSamplingParameters(
+        double clockFrequency,
+        SamplingMethod method,
+        double samplingFrequency
+    );
 
     /**
      * Clock SID forward using chosen output sampling algorithm.
@@ -267,6 +274,13 @@ public:
      */
     void setFilter6581Curve(double filterCurve);
 
+    /**
+    * Set filter range parameter for 6581 model
+    *
+    * @see Filter6581::setFilterRange(double)
+    */
+    void setFilter6581Range ( double adjustment );
+
     /**
      * Set filter curve parameter for 8580 model.
      *
@@ -312,13 +326,22 @@ void SID::ageBusValue(unsigned int n)
 }
 
 RESID_INLINE
-int SID::output() const
+int SID::output()
 {
-    const int v1 = voice[0]->output(voice[2]->wave());
+    const float o1 = voice[0].output(voice[2].wave());
-    const int v2 = voice[1]->output(voice[0]->wave());
+    const float o2 = voice[1].output(voice[0].wave());
-    const int v3 = voice[2]->output(voice[1]->wave());
+    const float o3 = voice[2].output(voice[1].wave());
 
-    return externalFilter->clock(filter->clock(v1, v2, v3));
+    const unsigned int env1 = voice[0].envelope()->output();
+    const unsigned int env2 = voice[1].envelope()->output();
+    const unsigned int env3 = voice[2].envelope()->output();
+
+    const int v1 = filter->getNormalizedVoice(o1, env1);
+    const int v2 = filter->getNormalizedVoice(o2, env2);
+    const int v3 = filter->getNormalizedVoice(o3, env3);
+
+    const int input = static_cast<int>(filter->clock(v1, v2, v3));
+    return externalFilter.clock(input);
 }
 
 
@@ -337,18 +360,18 @@ int SID::clock(unsigned int cycles, short* buf)
             for (unsigned int i = 0; i < delta_t; i++)
             {
                 // clock waveform generators
-                voice[0]->wave()->clock();
+                voice[0].wave()->clock();
-                voice[1]->wave()->clock();
+                voice[1].wave()->clock();
-                voice[2]->wave()->clock();
+                voice[2].wave()->clock();
 
                 // clock envelope generators
-                voice[0]->envelope()->clock();
+                voice[0].envelope()->clock();
-                voice[1]->envelope()->clock();
+                voice[1].envelope()->clock();
-                voice[2]->envelope()->clock();
+                voice[2].envelope()->clock();
 
                 if (unlikely(resampler->input(output())))
                 {
-                    buf[s++] = resampler->getOutput();
+                    buf[s++] = resampler->getOutput(scaleFactor);
                 }
             }
 

src/sound/resid-fp/siddefs-fp.h

@@ -26,10 +26,6 @@
 // Compiler specifics.
 #define HAVE_BUILTIN_EXPECT true
 
-#ifndef M_PI
-#  define M_PI    3.14159265358979323846
-#endif
-
 // Branch prediction macros, lifted off the Linux kernel.
 #if RESID_BRANCH_HINTS && HAVE_BUILTIN_EXPECT
 #  define likely(x)      __builtin_expect(!!(x), 1)
@@ -43,6 +39,8 @@ namespace reSIDfp {
 
 typedef enum { MOS6581=1, MOS8580 } ChipModel;
 
+typedef enum { AVERAGE=1, WEAK, STRONG } CombinedWaveforms;
+
 typedef enum { DECIMATE=1, RESAMPLE } SamplingMethod;
 }
 

src/sound/resid-fp/siddefs-fp.h.in

@@ -26,10 +26,6 @@
 // Compiler specifics.
 #define HAVE_BUILTIN_EXPECT @HAVE_BUILTIN_EXPECT@
 
-#ifndef M_PI
-#  define M_PI    3.14159265358979323846
-#endif
-
 // Branch prediction macros, lifted off the Linux kernel.
 #if RESID_BRANCH_HINTS && HAVE_BUILTIN_EXPECT
 #  define likely(x)      __builtin_expect(!!(x), 1)
@@ -43,6 +39,8 @@ namespace reSIDfp {
 
 typedef enum { MOS6581=1, MOS8580 } ChipModel;
 
+typedef enum { AVERAGE=1, WEAK, STRONG } CombinedWaveforms;
+
 typedef enum { DECIMATE=1, RESAMPLE } SamplingMethod;
 }
 

src/sound/snd_resid.cpp

@@ -25,7 +25,7 @@ sid_init(void)
 #if 0
     psid_t *psid;
 #endif
-    reSIDfp::SamplingMethod method         = reSIDfp::DECIMATE;
+    reSIDfp::SamplingMethod method         = reSIDfp::RESAMPLE;
     float                   cycles_per_sec = 14318180.0 / 16.0;
 
     psid = new psid_t;
@@ -34,8 +34,7 @@ sid_init(void)
 #endif
     psid->sid = new SID;
 
-    psid->sid->setChipModel(reSIDfp::MOS8580);
+    psid->sid->setChipModel(reSIDfp::MOS6581);
-    psid->sid->enableFilter(true);
 
     psid->sid->reset();
 
@@ -43,14 +42,13 @@ sid_init(void)
         psid->sid->write(c, 0);
 
     try {
-        psid->sid->setSamplingParameters(cycles_per_sec, method, (float) RESID_FREQ, 0.9 * (float) RESID_FREQ / 2.0);
+        psid->sid->setSamplingParameters(cycles_per_sec, method, (float) RESID_FREQ);
     } catch (reSIDfp::SIDError) {
 #if 0
         printf("reSID failed!\n");
 #endif
     }
 
-    psid->sid->setChipModel(reSIDfp::MOS6581);
     psid->sid->input(0);
 
     return (void *) psid;

src/sound/snd_ssi2001.c

@@ -125,8 +125,7 @@ static const device_config_t ssi2001_config[] = {
 // clang-format off
 };
 
-const device_t ssi2001_device =
+const device_t ssi2001_device = {
-{
     .name = "Innovation SSI-2001",
     .internal_name = "ssi2001",
     .flags = DEVICE_ISA,