Butterworth Online IIR filter

src/Filtering.Tests/Butterworth/OnlineIirTest.cs

@@ -0,0 +1,129 @@
+using MathNet.Filtering.Butterworth;
+using MathNet.Filtering.IIR;
+using NUnit.Framework;
+
+namespace MathNet.Filtering.Tests.Butterworth
+{
+    [TestFixture]
+    public class OnlineIirTest
+    {
+        private const double Tolerance = 10e-14;
+
+        private const double SamplingFrequency = 20000;
+
+        private int InputDataLength => InputData.Length;
+        private readonly double[] InputData = { 0.228897927516298, -2.619954349660920,
+                                                -17.502123684467897, -2.856509715953298,
+                                                -8.313665115676244, -9.792063051673022,
+                                                -11.564016556640022, -5.335571093159874,
+                                                -20.026357358830605, 9.642294226316274 };
+
+
+        [Test]
+        public void LowPass()
+        {
+            const double passbandFreq = 2000 / SamplingFrequency;
+            const double stopbandFreq = 2500 / SamplingFrequency;
+            const double passbandRipple = 5;
+            const double stopbandAttenuation = 6;
+            double[] ExpectedOutputData = {0.022257300742183576,
+                                          -0.21456981840243389,
+                                          -2.1294480226855419 ,
+                                          -3.6949348755798441 ,
+                                          -4.06251962388313   ,
+                                          -5.0330098214473757 ,
+                                          -6.1308190623358048 ,
+                                          -6.58179840726892   ,
+                                          -7.7679250740922559 ,
+                                          -7.2669818186786257 };
+
+            OnlineIirFilter filter = OnlineIirButterworthFilter.LowPass(passbandFreq, stopbandFreq, passbandRipple, stopbandAttenuation);
+
+            double[] output = RunFilter(filter);
+            Assert.That(output, Is.EqualTo(ExpectedOutputData).Within(Tolerance));
+        }
+
+        [Test]
+        public void HighPass()
+        {
+            const double passbandFreq = 2500 / SamplingFrequency;
+            const double stopbandFreq = 2000 / SamplingFrequency;
+            const double passbandRipple = 5;
+            const double stopbandAttenuation = 6;
+            double[] ExpectedOutputData = { 0.17709781344898412, -2.1072080720443367,
+                -12.667777308342028, 4.3969941901531158,
+                -1.8152939925238485, -2.1375169628929394,
+                -2.5410236580725973, 3.4279893350358845,
+                -9.4897624101777485, 17.759915828621352 };
+
+            OnlineIirFilter filter = OnlineIirButterworthFilter.HighPass(stopbandFreq, passbandFreq, passbandRipple, stopbandAttenuation);
+
+            double[] output = RunFilter(filter);
+            Assert.That(output, Is.EqualTo(ExpectedOutputData).Within(Tolerance));
+        }
+
+        [Test]
+        public void BandPass()
+        {
+            const double lowPassbandFreq = 2500 / SamplingFrequency;
+            const double lowStopbandFreq = 2000 / SamplingFrequency;
+            const double highPassbandFreq = 3000 / SamplingFrequency;
+            const double highStopbandFreq = 3500 / SamplingFrequency;
+            const double passbandRipple = 5;
+            const double stopbandAttenuation = 6;
+            double[] ExpectedOutputData = {0.0059568621907906414,
+                                          -0.057636016095946464,
+                                          -0.56911855227151587,
+                                          -0.95907837367120408,
+                                          -0.9193187331454562 ,
+                                          -0.89887865697280822,
+                                          -0.80447803861636846,
+                                          -0.45616524702168815,
+                                          -0.2652060530769399 ,
+                                          0.35269154712711837};
+
+            OnlineIirFilter filter = OnlineIirButterworthFilter.BandPass(lowStopbandFreq, lowPassbandFreq, highPassbandFreq, highStopbandFreq, passbandRipple, stopbandAttenuation);
+            double[] output = RunFilter(filter);
+            Assert.That(output, Is.EqualTo(ExpectedOutputData).Within(Tolerance));
+
+        }
+
+        [Test]
+        public void BandStop()
+        {
+            const double lowPassbandFreq = 2000 / SamplingFrequency;
+            const double lowStopbandFreq = 2500 / SamplingFrequency;
+            const double highPassbandFreq = 3500 / SamplingFrequency;
+            const double highStopbandFreq = 3000 / SamplingFrequency;
+            const double passbandRipple = 5;
+            const double stopbandAttenuation = 6;
+            double[] ExpectedOutputData = { 0.19763688241819988,
+                                            -2.3112043872011712,
+                                            -14.573275217116462,
+                                            1.5480395909989819,
+                                            -4.7847154775654985,
+                                            -6.5078504178710119,
+                                            -8.5307650341308729,
+                                            -3.5708389643430767,
+                                            -18.305730196535904,
+                                            9.0144565837006585};
+
+            OnlineIirFilter filter = OnlineIirButterworthFilter.BandStop(lowPassbandFreq, lowStopbandFreq, highStopbandFreq, highPassbandFreq, passbandRipple, stopbandAttenuation);
+            double[] output = RunFilter(filter);
+            Assert.That(output, Is.EqualTo(ExpectedOutputData).Within(Tolerance));
+
+        }
+
+        protected double[] RunFilter(OnlineIirFilter filter)
+        {
+            double[] output = new double[InputDataLength];
+            for (int i = 0; i < InputDataLength; i++)
+            {
+                output[i] = filter.ProcessSample(InputData[i]);
+                //double outputSample = filter.ProcessSample(InputData[i]);
+                //Assert.AreEqual(outputSample, ExpectedOutputData[i], Tolerance);
+            }
+            return output;
+        }
+    }
+}

src/Filtering/Butterworth/IirCoefficients.cs

@@ -147,7 +147,7 @@ public static (double[] numerator, double[] denominator) Notch(double centralFre
             var (gain, zeros, poles) = TransferFunction(n);
 
             wc1 = Helpers.MathFunctions.WarpFrequency(wc1, T);
-            wc1 = Helpers.MathFunctions.WarpFrequency(wc2, T);
+            wc2 = Helpers.MathFunctions.WarpFrequency(wc2, T);
             (gain, zeros, poles) = TransferFunctionTransformer.BandStop(gain, zeros, poles, wc1, wc2);
 
             return Coefficients(gain, zeros, poles, T);

src/Filtering/Butterworth/OnlineIirButterworthFilter.cs

@@ -0,0 +1,114 @@
+// <copyright file="IirCoefficients.cs" company="Math.NET">
+// Math.NET Filtering, part of the Math.NET Project
+// http://filtering.mathdotnet.com
+// http://github.com/mathnet/mathnet-filtering
+//
+// Copyright (c) 2009-2021 Math.NET
+//
+// Permission is hereby granted, free of charge, to any person
+// obtaining a copy of this software and associated documentation
+// files (the "Software"), to deal in the Software without
+// restriction, including without limitation the rights to use,
+// copy, modify, merge, publish, distribute, sublicense, and/or sell
+// copies of the Software, and to permit persons to whom the
+// Software is furnished to do so, subject to the following
+// conditions:
+//
+// The above copyright notice and this permission notice shall be
+// included in all copies or substantial portions of the Software.
+//
+// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+// EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES
+// OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
+// NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
+// HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+// WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
+// OTHER DEALINGS IN THE SOFTWARE.
+// </copyright>
+
+using MathNet.Filtering.IIR;
+
+namespace MathNet.Filtering.Butterworth
+{
+    /// <summary>
+    /// Butterworth filter design generating <see cref="OnlineIirFilter"/> objects.
+    /// </summary>
+    /// <seealso cref="IirCoefficients"/>
+    public static class OnlineIirButterworthFilter
+    {
+        /// <summary>
+        /// Generates an online IIR low-pass Butterworth filter.
+        /// </summary>
+        /// <param name="passbandFreq">Passband corner frequency (in Hz).</param>
+        /// <param name="stopbandFreq">Stopband corner frequency (in Hz).</param>
+        /// <param name="passbandRipple">Maximum allowed passband ripple.</param>
+        /// <param name="stopbandAttenuation">Minimum required stopband attenuation.</param>
+        /// <returns>Online IIR filter</returns>
+        public static OnlineIirFilter LowPass(double passbandFreq, double stopbandFreq, double passbandRipple, double stopbandAttenuation)
+        {
+            (double[] numerator, double[] denominator) = IirCoefficients.LowPass(passbandFreq, stopbandFreq, passbandRipple, stopbandAttenuation);
+            return new OnlineIirFilter(numerator, denominator);
+        }
+
+        /// <summary>
+        /// Generates an online IIR high-pass Butterworth filter.
+        /// </summary>
+        /// <param name="stopbandFreq">Stopband corner frequency (in Hz).</param>
+        /// <param name="passbandFreq">Passband corner frequency (in Hz).</param>
+        /// <param name="passbandRipple">Maximum allowed passband ripple.</param>
+        /// <param name="stopbandAttenuation">Minimum required stopband attenuation.</param>
+        /// <returns>Online IIR filter</returns>
+        public static OnlineIirFilter HighPass(double stopbandFreq, double passbandFreq, double passbandRipple, double stopbandAttenuation)
+        {
+            (double[] numerator, double[] denominator) = IirCoefficients.HighPass(stopbandFreq, passbandFreq, passbandRipple, stopbandAttenuation);
+            return new OnlineIirFilter(numerator, denominator);
+        }
+
+        /// <summary>
+        /// Generates an online IIR band-pass Butterworth filter.
+        /// </summary>
+        /// <param name="lowStopbandFreq">Lower stopband corner frequency (in Hz).</param>
+        /// <param name="lowPassbandFreq">Lower passband corner frequency (in Hz).</param>
+        /// <param name="highPassbandFreq">Higher passband corner frequency (in Hz).</param>
+        /// <param name="highStopbandFreq">Higher stopband corner frequency (in Hz).</param>
+        /// <param name="passbandRipple">Maximum allowed passband ripple.</param>
+        /// <param name="stopbandAttenuation">Minimum required stopband attenuation.</param>
+        /// <returns>Online IIR filter</returns>
+        public static OnlineIirFilter BandPass(double lowStopbandFreq, double lowPassbandFreq, double highPassbandFreq, double highStopbandFreq, double passbandRipple, double stopbandAttenuation)
+        {
+            (double[] numerator, double[] denominator) = IirCoefficients.BandPass(lowStopbandFreq, lowPassbandFreq, highPassbandFreq, highStopbandFreq, passbandRipple, stopbandAttenuation);
+            return new OnlineIirFilter(numerator, denominator);
+        }
+
+        /// <summary>
+        /// Generates an online IIR band-stop Butterworth filter.
+        /// </summary>
+        /// <param name="lowPassbandFreq">Lower passband corner frequency (in Hz).</param>
+        /// <param name="lowStopbandFreq">Lower stopband corner frequency (in Hz).</param>
+        /// <param name="highStopbandFreq">Higher stopband corner frequency (in Hz).</param>
+        /// <param name="highPassbandFreq">Higher passband corner frequency (in Hz).</param>
+        /// <param name="passbandRipple">Maximum allowed passband ripple.</param>
+        /// <param name="stopbandAttenuation">Minimum required stopband attenuation.</param>
+        /// <returns>Online IIR filter</returns>
+        public static OnlineIirFilter BandStop(double lowPassbandFreq, double lowStopbandFreq, double highStopbandFreq, double highPassbandFreq, double passbandRipple, double stopbandAttenuation)
+        {
+            (double[] numerator, double[] denominator) = IirCoefficients.BandStop(lowPassbandFreq, lowStopbandFreq, highStopbandFreq, highPassbandFreq, passbandRipple, stopbandAttenuation);
+            return new OnlineIirFilter(numerator, denominator);
+        }
+
+        /// <summary>
+        /// Generates an online IIR notch Butterworth filter.
+        /// </summary>
+        /// <param name="centralFreq">Filter central frequency.</param>
+        /// <param name="Q">Quality factor.</param>
+        /// <param name="passbandRipple">Maximum allowed passband ripple.</param>
+        /// <param name="stopbandAttenuation">Minimum required stopband attenuation.</param>
+        /// <returns>Online IIR filter</returns>
+        public static OnlineIirFilter Notch(double centralFreq, double Q, double passbandRipple, double stopbandAttenuation)
+        {
+            (double[] numerator, double[] denominator) = IirCoefficients.Notch(centralFreq, Q, passbandRipple, stopbandAttenuation);
+            return new OnlineIirFilter(numerator, denominator);
+        }
+    }
+}

src/Filtering/Helpers/MathFunctions.cs

@@ -28,7 +28,6 @@
 // </copyright>
 
 using System;
-using System.Linq;
 using System.Numerics;
 using MathNet.Numerics;
 

src/Filtering/IIR/OnlineIirFilter.cs

@@ -3,7 +3,7 @@
 // http://filtering.mathdotnet.com
 // http://github.com/mathnet/mathnet-filtering
 //
-// Copyright (c) 2009-2014 Math.NET
+// Copyright (c) 2009-2021 Math.NET
 //
 // Permission is hereby granted, free of charge, to any person
 // obtaining a copy of this software and associated documentation
@@ -40,7 +40,7 @@ namespace MathNet.Filtering.IIR
     /// Filter implements the canonic Direct Form II structure.
     /// </summary>
     /// <remarks>
-    /// System Description: H(z) = (b0 + b1*z^-1 + b2*z^-2) / (1 + a1*z^-1 + a2*z^-2)
+    /// System Description: H(z) = (b0 + b1*z^-1 + b2*z^-2) / (a0 + a1*z^-1 + a2*z^-2)
     /// </remarks>
     public class OnlineIirFilter : OnlineFilter
     {
@@ -54,6 +54,9 @@ public class OnlineIirFilter : OnlineFilter
         /// <summary>
         /// Infinite Impulse Response (IIR) Filter.
         /// </summary>
+        /// <param name="coefficients">Vector of IIR coefficients in the following form: [b0 b1 b2 ... a0 a1 a2 ...]</param>
+        /// <exception cref="ArgumentNullException"></exception>
+        /// <exception cref="ArgumentException"></exception>
         public OnlineIirFilter(double[] coefficients)
         {
             if (null == coefficients)
@@ -81,6 +84,50 @@ public OnlineIirFilter(double[] coefficients)
             _bufferY = new double[size];
         }
 
+        /// <summary>
+        /// Infinite Impulse Response (IIR) Filter
+        /// </summary>
+        /// <param name="numerator">Numerator coefficients [b0 b1 ...]</param>
+        /// <param name="denominator">Denominator coefficients [a0 a1 ...]</param>
+        /// <exception cref="ArgumentNullException"></exception>
+        /// <exception cref="ArgumentException"></exception>
+        public OnlineIirFilter(double[] numerator, double[] denominator)
+        {
+            if (null == numerator)
+            {
+                throw new ArgumentNullException(nameof(numerator));
+            }
+            if (null == denominator)
+            {
+                throw new ArgumentNullException(nameof(denominator));
+            }
+            if (denominator.Length != numerator.Length)
+            {
+                throw new ArgumentException(Resources.ArgumentCoefficientLengthMismatch);
+            }
+            if (denominator.Length == 0)
+            {
+                throw new ArgumentException(Resources.ArgumentEmptyCoefficientVector, nameof(denominator));
+            }
+            if (numerator.Length == 0)
+            {
+                throw new ArgumentException(Resources.ArgumentEmptyCoefficientVector, nameof(numerator));
+            }
+
+            _halfSize = denominator.Length;
+            int size = denominator.Length * 2;
+            _b = new double[size];
+            _a = new double[size];
+            for (int i = 0; i < _halfSize; i++)
+            {
+                _b[i] = _b[_halfSize + i] = numerator[i];
+                _a[i] = _a[_halfSize + i] = denominator[i];
+            }
+
+            _bufferX = new double[size];
+            _bufferY = new double[size];
+        }
+
         /// <summary>
         /// Process a single sample.
         /// </summary>