feature: replace OPL3 with NukedOpl3

src/Spice86.Core/Backend/Audio/Iir/BandPassTransform.cs

@@ -0,0 +1,120 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ * 
+ *  Copyright (c) 2009 by Vinnie Falco
+ *  Copyright (c) 2016 by Bernd Porr
+ */
+
+using System.Numerics;
+
+namespace Spice86.Core.Backend.Audio.Iir;
+
+/**
+ * Transforms from an analogue bandpass filter to a digital bandstop filter
+ */
+public class BandPassTransform {
+    private readonly double _wc2;
+    private readonly double _wc;
+    private readonly double _a, _b;
+    private readonly double _a2, _b2;
+    private readonly double _ab, _ab2;
+
+    public BandPassTransform(double fc, double fw, LayoutBase digital,
+            LayoutBase analog) {
+
+        digital.Reset();
+
+        if (fc < 0) {
+            throw new ArithmeticException("Cutoff frequency cannot be negative.");
+        }
+
+        if (!(fc < 0.5)) {
+            throw new ArithmeticException("Cutoff frequency must be less than the Nyquist frequency.");
+        }
+
+        double ww = 2 * Math.PI * fw;
+
+        // pre-calcs
+        _wc2 = (2 * Math.PI * fc) - ww / 2;
+        _wc = _wc2 + ww;
+
+        // what is this crap?
+        if (_wc2 < 1e-8) {
+            _wc2 = 1e-8;
+        }
+
+        if (_wc > Math.PI - 1e-8) {
+            _wc = Math.PI - 1e-8;
+        }
+
+        _a = Math.Cos((_wc + _wc2) * 0.5) / Math.Cos((_wc - _wc2) * 0.5);
+        _b = 1 / Math.Tan((_wc - _wc2) * 0.5);
+        _a2 = _a * _a;
+        _b2 = _b * _b;
+        _ab = _a * _b;
+        _ab2 = 2 * _ab;
+
+        int numPoles = analog.NumPoles;
+        int pairs = numPoles / 2;
+        for (int i = 0; i < pairs; ++i) {
+            PoleZeroPair pair = analog.GetPair(i);
+            ComplexPair p1 = Transform(pair.poles.First);
+            ComplexPair z1 = Transform(pair.zeros.First);
+
+            digital.AddPoleZeroConjugatePairs(p1.First, z1.First);
+            digital.AddPoleZeroConjugatePairs(p1.Second, z1.Second);
+        }
+
+        if ((numPoles & 1) == 1) {
+            ComplexPair poles = Transform(analog.GetPair(pairs).poles.First);
+            ComplexPair zeros = Transform(analog.GetPair(pairs).zeros.First);
+
+            digital.Add(poles, zeros);
+        }
+
+        double wn = analog.NormalW;
+        digital.SetNormal(
+                2 * Math.Atan(Math.Sqrt(Math.Tan((_wc + wn) * 0.5)
+                        * Math.Tan((_wc2 + wn) * 0.5))), analog.NormalGain);
+    }
+
+    private ComplexPair Transform(Complex c) {
+        if (Complex.IsInfinity(c)) {
+            return new ComplexPair(new Complex(-1, 0), new Complex(1, 0));
+        }
+
+        c = new Complex(1, 0).Add(c).Divide(new Complex(1, 0).Subtract(c)); // bilinear
+
+        var v = new Complex(0, 0);
+        v = MathSupplement.AddMul(v, 4 * (_b2 * (_a2 - 1) + 1), c);
+        v = v.Add(8 * (_b2 * (_a2 - 1) - 1));
+        v = v.Multiply(c);
+        v = v.Add(4 * (_b2 * (_a2 - 1) + 1));
+        v = v.Sqrt();
+
+        Complex u = v.Multiply(-1);
+        u = MathSupplement.AddMul(u, _ab2, c);
+        u = u.Add(_ab2);
+
+        v = MathSupplement.AddMul(v, _ab2, c);
+        v = v.Add(_ab2);
+
+        var d = new Complex(0, 0);
+        d = MathSupplement.AddMul(d, 2 * (_b - 1), c).Add(2 * (1 + _b));
+
+        return new ComplexPair(u.Divide(d), v.Divide(d));
+    }
+}

src/Spice86.Core/Backend/Audio/Iir/BandStopTransform.cs

@@ -0,0 +1,122 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ * 
+ *  Copyright (c) 2009 by Vinnie Falco
+ *  Copyright (c) 2016 by Bernd Porr
+ */
+
+using System.Numerics;
+
+namespace Spice86.Core.Backend.Audio.Iir;
+
+/**
+ * Transforms from an analogue lowpass filter to a digital bandstop filter
+ */
+public class BandStopTransform {
+    private readonly double _wc;
+    private readonly double _wc2;
+    private readonly double _a;
+    private readonly double _b;
+    private readonly double _a2;
+    private readonly double _b2;
+
+
+    public BandStopTransform(double fc,
+        double fw,
+        LayoutBase digital,
+        LayoutBase analog) {
+        digital.Reset();
+
+        if (fc < 0) {
+            throw new ArithmeticException("Cutoff frequency cannot be negative.");
+        }
+
+        if (!(fc < 0.5)) {
+            throw new ArithmeticException("Cutoff frequency must be less than the Nyquist frequency.");
+        }
+
+        double ww = 2 * Math.PI * fw;
+
+        _wc2 = 2 * Math.PI * fc - ww / 2;
+        _wc = _wc2 + ww;
+
+        // this is crap
+        if (_wc2 < 1e-8) {
+            _wc2 = 1e-8;
+        }
+
+        if (_wc > Math.PI - 1e-8) {
+            _wc = Math.PI - 1e-8;
+        }
+
+        _a = Math.Cos((_wc + _wc2) * .5) /
+             Math.Cos((_wc - _wc2) * .5);
+        _b = Math.Tan((_wc - _wc2) * .5);
+        _a2 = _a * _a;
+        _b2 = _b * _b;
+
+        int numPoles = analog.NumPoles;
+        int pairs = numPoles / 2;
+        for (int i = 0; i < pairs; i++) {
+            PoleZeroPair pair = analog.GetPair(i);
+            ComplexPair p = Transform(pair.poles.First);
+            ComplexPair z = Transform(pair.zeros.First);
+            digital.AddPoleZeroConjugatePairs(p.First, z.First);
+            digital.AddPoleZeroConjugatePairs(p.Second, z.Second);
+        }
+
+        if ((numPoles & 1) == 1) {
+            ComplexPair poles = Transform(analog.GetPair(pairs).poles.First);
+            ComplexPair zeros = Transform(analog.GetPair(pairs).zeros.First);
+
+            digital.Add(poles, zeros);
+        }
+
+        if (fc < 0.25) {
+            digital.SetNormal(Math.PI, analog.NormalGain);
+        } else {
+            digital.SetNormal(0, analog.NormalGain);
+        }
+    }
+
+    private ComplexPair Transform(Complex c) {
+        if (c == Complex.Infinity) {
+            c = new Complex(-1, 0);
+        } else {
+            c = new Complex(1, 0).Add(c).Divide(new Complex(1, 0).Subtract(c)); // bilinear
+        }
+
+        var u = new Complex(0, 0);
+        u = MathSupplement.AddMul(u, 4 * (_b2 + _a2 - 1), c);
+        u = u.Add(8 * (_b2 - _a2 + 1));
+        u = u.Multiply(c);
+        u = u.Add(4 * (_a2 + _b2 - 1));
+        u = u.Sqrt();
+
+        Complex v = u.Multiply(-.5);
+        v = v.Add(_a);
+        v = MathSupplement.AddMul(v, -_a, c);
+
+        u = u.Multiply(.5);
+        u = u.Add(_a);
+        u = MathSupplement.AddMul(u, -_a, c);
+
+        var d = new Complex(_b + 1, 0);
+        d = MathSupplement.AddMul(d, _b - 1, c);
+
+        return new ComplexPair(u.Divide(d), v.Divide(d));
+    }
+}

src/Spice86.Core/Backend/Audio/Iir/Biquad.cs

@@ -0,0 +1,148 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ * 
+ *  Copyright (c) 2009 by Vinnie Falco
+ *  Copyright (c) 2016 by Bernd Porr
+ */
+
+using System.Numerics;
+
+namespace Spice86.Core.Backend.Audio.Iir;
+
+/**
+ * Contains the coefficients of a 2nd order digital filter with two poles and two zeros
+ */
+public class Biquad {
+    public double A0 { get; set; }
+    public double A1 { get; set; }
+    public double A2 { get; set; }
+    public double B1 { get; set; }
+    public double B2 { get; set; }
+    public double B0 { get; set; }
+
+    public double GetA0 => A0;
+
+    public double GetA1 => A1 * A0;
+
+    public double GetA2 => A2 * A0;
+
+    public double GetB0 => B0 * A0;
+
+    public double GetB1 => B1 * A0;
+
+    public double GetB2 => B2 * A0;
+
+    public Complex Response(double normalizedFrequency) {
+        double a0 = GetA0;
+        double a1 = GetA1;
+        double a2 = GetA2;
+        double b0 = GetB0;
+        double b1 = GetB1;
+        double b2 = GetB2;
+
+        double w = 2 * Math.PI * normalizedFrequency;
+        Complex czn1 = ComplexUtils.PolarToComplex(1.0, -w);
+        Complex czn2 = ComplexUtils.PolarToComplex(1.0, -2 * w);
+        var ch = new Complex(1, 0);
+        var cbot = new Complex(1, 0);
+
+        var ct = new Complex(b0 / a0, 0);
+        var cb = new Complex(1, 0);
+        ct = MathSupplement.AddMul(ct, b1 / a0, czn1);
+        ct = MathSupplement.AddMul(ct, b2 / a0, czn2);
+        cb = MathSupplement.AddMul(cb, a1 / a0, czn1);
+        cb = MathSupplement.AddMul(cb, a2 / a0, czn2);
+        ch = Complex.Multiply(ch, ct);
+        cbot = Complex.Multiply(cbot, cb);
+
+        return Complex.Divide(ch, cbot);
+    }
+
+    public void SetCoefficients(double a0, double a1, double a2,
+                double b0, double b1, double b2) {
+        A0 = a0;
+        A1 = a1 / a0;
+        A2 = a2 / a0;
+        B0 = b0 / a0;
+        B1 = b1 / a0;
+        B2 = b2 / a0;
+    }
+
+    public void SetOnePole(Complex pole, Complex zero) {
+        const double a0 = 1;
+        double a1 = -pole.Real;
+        const double a2 = 0;
+        double b0 = -zero.Real;
+        const double b1 = 1;
+        const double b2 = 0;
+        SetCoefficients(a0, a1, a2, b0, b1, b2);
+    }
+
+    public void SetTwoPole(
+        Complex pole1, Complex zero1,
+        Complex pole2, Complex zero2) {
+        const double a0 = 1;
+        double a1;
+        double a2;
+
+        if (pole1.Imaginary != 0) {
+            a1 = -2 * pole1.Real;
+            a2 = Complex.Abs(pole1) * Complex.Abs(pole1);
+        } else {
+            a1 = -(pole1.Real + pole2.Real);
+            a2 = pole1.Real * pole2.Real;
+        }
+
+        const double b0 = 1;
+        double b1;
+        double b2;
+
+        if (zero1.Imaginary != 0) {
+            b1 = -2 * zero1.Real;
+            b2 = Complex.Abs(zero1) * Complex.Abs(zero1);
+        } else {
+            b1 = -(zero1.Real + zero2.Real);
+            b2 = zero1.Real * zero2.Real;
+        }
+
+        SetCoefficients(a0, a1, a2, b0, b1, b2);
+    }
+
+    public void SetPoleZeroForm(BiquadPoleState bps) {
+        SetPoleZeroPair(bps);
+        ApplyScale(bps.Gain);
+    }
+
+    public void SetIdentity() {
+        SetCoefficients(1, 0, 0, 1, 0, 0);
+    }
+
+    public void ApplyScale(double scale) {
+        B0 *= scale;
+        B1 *= scale;
+        B2 *= scale;
+    }
+
+
+    public void SetPoleZeroPair(PoleZeroPair pair) {
+        if (pair.IsSinglePole()) {
+            SetOnePole(pair.poles.First, pair.zeros.First);
+        } else {
+            SetTwoPole(pair.poles.First, pair.zeros.First,
+                    pair.poles.Second, pair.zeros.Second);
+        }
+    }
+}