Add custom band limited saw generator and chamberlin two-pole low pas…

…s multi output filter.

packages/math/include/math/MathFunc.h

@@ -125,7 +125,6 @@ namespace l::math::functions {
 		}
 	}
 
-
 	template<class T>
 	T cos(T val) {
 		if constexpr (std::is_floating_point_v<T>) {
@@ -138,6 +137,54 @@ namespace l::math::functions {
 		}
 	}
 
+	template<class T>
+	T tan(T val) {
+		if constexpr (std::is_floating_point_v<T>) {
+			if constexpr (sizeof(T) == 4) {
+				return tanf(val);
+			}
+			else if constexpr (sizeof(T) == 8) {
+				return ::tan(val);
+			}
+		}
+	}
+
+	template<class T>
+	T asin(T val) {
+		if constexpr (std::is_floating_point_v<T>) {
+			if constexpr (sizeof(T) == 4) {
+				return asinf(val);
+			}
+			else if constexpr (sizeof(T) == 8) {
+				return ::asin(val);
+			}
+		}
+	}
+
+	template<class T>
+	T acos(T val) {
+		if constexpr (std::is_floating_point_v<T>) {
+			if constexpr (sizeof(T) == 4) {
+				return acosf(val);
+			}
+			else if constexpr (sizeof(T) == 8) {
+				return ::acos(val);
+			}
+		}
+	}
+
+	template<class T>
+	T atan(T val) {
+		if constexpr (std::is_floating_point_v<T>) {
+			if constexpr (sizeof(T) == 4) {
+				return atanf(val);
+			}
+			else if constexpr (sizeof(T) == 8) {
+				return ::atan(val);
+			}
+		}
+	}
+
 	template<class T>
 	T mod(T val, T mod) {
 		if constexpr (std::is_floating_point_v<T>) {

packages/nodegraph/include/nodegraph/NodeGraphSchema.h

@@ -59,6 +59,7 @@ namespace l::nodegraph {
             RegisterNodeType("Logic", 102, "Xor");
             RegisterNodeType("Filter", 150, "Lowpass");
             RegisterNodeType("Filter", 151, "Highpass");
+            RegisterNodeType("Filter", 152, "Chamberlin two-pole (4 mode)");
             RegisterNodeType("Output", 200, "Debug");
             RegisterNodeType("Output", 201, "Speaker");
             RegisterNodeType("Output", 202, "Plot");
@@ -84,6 +85,7 @@ namespace l::nodegraph {
             RegisterNodeType("Signal", 353, "Sine FM 3");
             RegisterNodeType("Signal", 354, "Saw");
             RegisterNodeType("Signal", 355, "Sine 2");
+            RegisterNodeType("Signal", 356, "Saw 2");
 
         }
 

packages/nodegraph/include/nodegraph/operations/NodeGraphOpFilter.h

@@ -86,5 +86,49 @@ namespace l::nodegraph {
         float mState1 = 0.0f;
     };
 
+    /*********************************************************************/
+    // source: https://www.musicdsp.org/en/latest/Filters/23-state-variable.html
+    class GraphFilterChamberlain2pole : public NodeGraphOp {
+    public:
+        std::string defaultInStrings[4] = { "In", "Cutoff", "Resonance", "Mode"};
+        std::string defaultOutStrings[1] = { "Out"};
+
+        GraphFilterChamberlain2pole(NodeGraphBase* node) :
+            NodeGraphOp(node, 4, 1)
+        {
+            mState.resize(4);
+        }
+
+        virtual ~GraphFilterChamberlain2pole() = default;
+        virtual void Reset() override;
+        virtual void Process(int32_t numSamples, std::vector<NodeGraphInput>& inputs, std::vector<NodeGraphOutput>& outputs) override;
+        virtual bool IsDataVisible(int8_t) override { return true; }
+        virtual bool IsDataEditable(int8_t channel) override { return channel > 0 ? true : false; }
+        virtual std::string_view GetInputName(int8_t inputChannel) override {
+            return defaultInStrings[inputChannel];
+        }
+
+        virtual std::string_view GetOutputName(int8_t outputChannel) override {
+            return defaultOutStrings[outputChannel];
+        }
+
+        virtual std::string_view GetName() override {
+            return "Chamberlin two-pole";
+        }
+    protected:
+        float mSamplesUntilUpdate = 0.0f;
+        float mUpdateSamples = 16.0f;
+
+        float mInputValuePrev = 0.0f;
+        float mCutoff = 0.0f;
+        float mResonance = 0.0f;
+
+        float mSampleRate = 44100.0f;
+        float mFreq = 0.0f;
+        float mScale = 0.0f;
+
+        std::vector<float> mState;
+    };
+
 }
 

packages/nodegraph/include/nodegraph/operations/NodeGraphOpSignal.h

@@ -25,15 +25,15 @@ namespace l::nodegraph {
     class GraphSignalBase : public NodeGraphOp {
     public:
 
-        static const int8_t mNumDefaultInputs = 7;
+        static const int8_t mNumDefaultInputs = 4;
         static const int8_t mNumDefaultOutputs = 1;
 
         GraphSignalBase(NodeGraphBase* node, std::string_view name, int32_t numInputs = 0, int32_t numOutputs = 0, int32_t numConstants = 0) :
             NodeGraphOp(node, mNumDefaultInputs + numInputs, mNumDefaultOutputs + numOutputs, numConstants),
             mName(name)
         {}
 
-        std::string defaultInStrings[mNumDefaultInputs] = { "Reset", "Freq", "Volume", "Smooth", "Cutoff", "Resonance", "Phase expansion"};
+        std::string defaultInStrings[mNumDefaultInputs] = { "Reset", "Freq", "Volume", "Smooth"};
         std::string defaultOutStrings[mNumDefaultOutputs] = { "Out" };
 
         virtual ~GraphSignalBase() = default;
@@ -73,12 +73,6 @@ namespace l::nodegraph {
         float mVolumeTarget = 0.0f;
         float mSamplesUntilUpdate = 0.0f;
         float mUpdateSamples = 256.0f;
-
-        // high pass
-        float mHPCutoff = 0.5f;
-        float mHPResonance = 0.0001f;
-        float mHPState0 = 0.0f;
-        float mHPState1 = 0.0f;
     };
 
     /*********************************************************************/
@@ -104,6 +98,139 @@ namespace l::nodegraph {
         float mPhaseFmod = 0.0f;
     };
 
+    /*********************************************************************/
+    struct LowpassType {
+        double x1 = 0.0;
+        double y1 = 0.0;
+        double a = 0.0;
+        double b = 0.0;
+    };
+
+    struct WaveformBlit {
+        double phase = 0.0;    /* phase accumulator */
+        double aNQ = 0.0;      /* attenuation at nyquist */
+        double curcps = 0.0;   /* current frequency, updated once per cycle */
+        double curper = 0.0;   /* current period, updated once per cycle */
+        LowpassType leaky; /* leaky integrator */
+        double N = 0.0;        /* # partials */
+        double a = 0.0;        /* dsf parameter which controls roll-off */
+        double aN = 0.0;       /* former to the N */
+    };
+
+    class GraphSignalSaw2 : public GraphSignalBase {
+    public:
+        GraphSignalSaw2(NodeGraphBase* node) :
+            GraphSignalBase(node, "Saw 2", 2)
+        {}
+        std::string extraString[2] = { "Attenuation", "Cutoff" };
+
+        virtual ~GraphSignalSaw2() = default;
+        virtual std::string_view GetInputNameExtra(int8_t extraInputChannel) override {
+            if (extraInputChannel < 2) return extraString[static_cast<uint8_t>(extraInputChannel)];
+            return "";
+        }
+        void ResetSignal() override;
+        void UpdateSignal(std::vector<NodeGraphInput>& inputs, std::vector<NodeGraphOutput>& outputs) override;
+        float GenerateSignal(float deltaTime, float freq, float deltaPhase) override;
+
+        void InitSaw(WaveformBlit* b, double aNQ, double cutoff)
+        {
+            b->phase = 0.0;
+            b->aNQ = aNQ;
+            b->curcps = 0.0;
+            b->curper = 0.0;
+            InitLowpass(&b->leaky, cutoff);
+        }
+
+        void UpdateSaw(WaveformBlit* b, double aNQ, double cutoff) {
+            b->aNQ = aNQ;
+            UpdateLowpass(&b->leaky, cutoff + 0.00001);
+        }
+
+        /* Returns a sawtooth computed from a leaky integration
+         * of a DSF bandlimited impulse train.
+         *
+         * cps (cycles per sample) is the fundamental
+         * frequency: 0 -> 0.5 == 0 -> nyquist
+         */
+
+        double ProcessSaw(WaveformBlit* b, double cps) {
+            double P2, beta, Nbeta, cosbeta, n, d, blit, saw;
+
+            if (b->phase >= 1.0 || b->curcps == 0.0)
+            {
+                /* New cycle, update frequency and everything
+                 * that depends on it
+                 */
+
+                if (b->phase >= 1.0)
+                    b->phase -= 1.0;
+                double cpsClamped = l::math::functions::clamp(cps, 2.0 / 44100, 0.5);
+                b->curcps = cpsClamped;        /* this cycle\'s frequency */
+                b->curper = 1.0 / cpsClamped;  /* this cycle\'s period */
+
+                P2 = b->curper / 2.0;
+                b->N = 1.0 + l::math::functions::floor(P2); /* # of partials incl. dc */
+
+                /* find the roll-off parameter which gives
+                 * the desired attenuation at nyquist
+                 */
+
+                b->a = l::math::functions::pow(b->aNQ, 1.0 / P2);
+                b->aN = l::math::functions::pow(b->a, b->N);
+            }
+
+            beta = 2.0 * l::math::constants::PI * b->phase;
+
+            Nbeta = b->N * beta;
+            cosbeta = l::math::functions::cos(beta);
+
+            /* The dsf blit is scaled by 1 / period to give approximately the same
+             * peak-to-peak over a wide range of frequencies.
+             */
+
+            n = 1.0 -
+                b->aN * l::math::functions::cos(Nbeta) -
+                b->a * (cosbeta - b->aN * l::math::functions::cos(Nbeta - beta));
+            d = b->curper * (1.0 + b->a * (-2.0 * cosbeta + b->a));
+
+            b->phase += b->curcps; /* update phase */
+
+            blit = n / d - b->curcps; /* This division can only fail if |a| == 1.0
+               * Subtracting the fundamental frq rids of DC
+               */
+
+            saw = ProcessLowpass(&b->leaky, blit); /* shape blit spectrum into a saw */
+
+            return 4.0 * saw;
+        }
+
+        void InitLowpass(LowpassType* lp, double cutoff) {
+            lp->x1 = lp->y1 = 0.0;
+            UpdateLowpass(lp, cutoff);
+        }
+
+        void UpdateLowpass(LowpassType* lp, double cutoff) {
+            double Omega = l::math::functions::atan(l::math::constants::PI * cutoff);
+            lp->a = -(1.0 - Omega) / (1.0 + Omega);
+            lp->b = (1.0 - lp->b) / 2.0;
+        }
+
+        double ProcessLowpass(LowpassType* lp, double x) {
+            double y;
+            y = lp->b * (x + lp->x1) - lp->a * lp->y1;
+            lp->x1 = x;
+            lp->y1 = y;
+            return y;
+        }
+
+    protected:
+        double mAttenuation = 0.0f;
+        double mCutoff = 0.0f;
+
+        WaveformBlit mSaw;
+    };
+
     /*********************************************************************/
     class GraphSignalSine : public NodeGraphOp {
     public:

packages/nodegraph/source/common/NodeGraphSchema.cpp

@@ -79,6 +79,9 @@ namespace l::nodegraph {
         case 151:
             node = mMainNodeGraph.NewNode<l::nodegraph::GraphFilterHighpass>(OutputType::Default);
             break;
+        case 152:
+            node = mMainNodeGraph.NewNode<l::nodegraph::GraphFilterChamberlain2pole>(OutputType::Default);
+            break;
         case 200:
             node = mMainNodeGraph.NewNode<l::nodegraph::GraphOutputDebug>(OutputType::ExternalOutput);
             break;
@@ -154,6 +157,9 @@ namespace l::nodegraph {
         case 355:
             node = mMainNodeGraph.NewNode<l::nodegraph::GraphSignalSine2>(OutputType::Default);
             break;
+        case 356:
+            node = mMainNodeGraph.NewNode<l::nodegraph::GraphSignalSaw2>(OutputType::Default);
+            break;
 
 
         default:

packages/nodegraph/source/common/operations/NodeGraphOpFilter.cpp

@@ -59,4 +59,52 @@ namespace l::nodegraph {
 
         outputs.at(0).mOutput = inputValue - mState1;
     }
+
+    /*********************************************************************/
+    void GraphFilterChamberlain2pole::Reset() {
+        for (int32_t i = 0; i < 4; i++) {
+            mState.at(i) = 0.0f;
+        }
+
+        mSamplesUntilUpdate = 0.0f;
+        mUpdateSamples = 4.0f;
+
+        mNode->SetInput(1, 0.99f);
+        mNode->SetInput(2, 0.01f);
+        mNode->SetInput(3, 0.0f);
+        mNode->SetInputBound(1, InputBound::INPUT_0_TO_1);
+        mNode->SetInputBound(2, InputBound::INPUT_0_TO_1);
+        mNode->SetInputBound(3, InputBound::INPUT_0_TO_1);
+    }
+
+    void GraphFilterChamberlain2pole::Process(int32_t numSamples, std::vector<NodeGraphInput>& inputs, std::vector<NodeGraphOutput>& outputs) {
+        auto input = &inputs.at(0).Get(numSamples);
+
+        auto mode = static_cast<int32_t>(3.0f * inputs.at(3).Get() + 0.5f);
+        auto output0 = &outputs.at(0).GetOutput(numSamples);
+
+        mSamplesUntilUpdate = l::audio::BatchUpdate(mUpdateSamples, mSamplesUntilUpdate, 0, numSamples,
+            [&]() {
+                mCutoff = inputs.at(1).Get();
+                mResonance = 1.0f - inputs.at(2).Get();
+
+                mFreq = l::math::functions::sin(l::math::constants::PI_f * mCutoff * mCutoff / 2.0f);
+                mScale = l::math::functions::sqrt(mResonance);
+            },
+            [&](int32_t start, int32_t end, bool) {
+                for (int32_t i = start; i < end; i++) {
+                    float inputValue = *input++;
+                    float inputValueInbetween = (mInputValuePrev + inputValue) * 0.5f;
+                    for (int32_t oversample = 0; oversample < 2; oversample++) {
+                        mState.at(0) = mState.at(0) + mFreq * mState.at(2);
+                        mState.at(1) = mScale * (oversample == 0 ? inputValueInbetween : inputValue) - mState.at(0) - mResonance * mState.at(2);
+                        mState.at(2) = mFreq * mState.at(1) + mState.at(2);
+                        mState.at(3) = mState.at(1) + mState.at(0);
+                    }
+                    *output0++ = mState.at(mode);
+                    mInputValuePrev = inputValue;
+                }
+            }
+        );
+    }
 }