From d55cecd568a53d58fccca11fd6393ba45aab81a2 Mon Sep 17 00:00:00 2001 From: StefanoPietrosanti Date: Wed, 13 Oct 2021 12:03:11 +0100 Subject: [PATCH 1/2] Added static class that generates OnlineIirFilter objects from Butterworth filters Butterworth filters were only generating coefficients, while now it is possible to directly obtain an OnlineIirFilter object designed using the Butterworth filter designer. Added tests to check that the online filters produce the expected output. --- .../Butterworth/OnlineIirTest.cs | 129 ++++++++++++++++++ .../Butterworth/OnlineIirButterworthFilter.cs | 114 ++++++++++++++++ src/Filtering/IIR/OnlineIirFilter.cs | 51 ++++++- .../Properties/Resources.Designer.cs | 21 ++- src/Filtering/Properties/Resources.resx | 10 +- 5 files changed, 319 insertions(+), 6 deletions(-) create mode 100644 src/Filtering.Tests/Butterworth/OnlineIirTest.cs create mode 100644 src/Filtering/Butterworth/OnlineIirButterworthFilter.cs diff --git a/src/Filtering.Tests/Butterworth/OnlineIirTest.cs b/src/Filtering.Tests/Butterworth/OnlineIirTest.cs new file mode 100644 index 00000000..fe459121 --- /dev/null +++ b/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; + } + } +} diff --git a/src/Filtering/Butterworth/OnlineIirButterworthFilter.cs b/src/Filtering/Butterworth/OnlineIirButterworthFilter.cs new file mode 100644 index 00000000..9d53a783 --- /dev/null +++ b/src/Filtering/Butterworth/OnlineIirButterworthFilter.cs @@ -0,0 +1,114 @@ +// +// 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. +// + +using MathNet.Filtering.IIR; + +namespace MathNet.Filtering.Butterworth +{ + /// + /// Butterworth filter design generating objects. + /// + /// + public static class OnlineIirButterworthFilter + { + /// + /// Generates an online IIR low-pass Butterworth filter. + /// + /// Passband corner frequency (in Hz). + /// Stopband corner frequency (in Hz). + /// Maximum allowed passband ripple. + /// Minimum required stopband attenuation. + /// Online IIR filter + 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); + } + + /// + /// Generates an online IIR high-pass Butterworth filter. + /// + /// Stopband corner frequency (in Hz). + /// Passband corner frequency (in Hz). + /// Maximum allowed passband ripple. + /// Minimum required stopband attenuation. + /// Online IIR filter + 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); + } + + /// + /// Generates an online IIR band-pass Butterworth filter. + /// + /// Lower stopband corner frequency (in Hz). + /// Lower passband corner frequency (in Hz). + /// Higher passband corner frequency (in Hz). + /// Higher stopband corner frequency (in Hz). + /// Maximum allowed passband ripple. + /// Minimum required stopband attenuation. + /// Online IIR filter + 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); + } + + /// + /// Generates an online IIR band-stop Butterworth filter. + /// + /// Lower passband corner frequency (in Hz). + /// Lower stopband corner frequency (in Hz). + /// Higher stopband corner frequency (in Hz). + /// Higher passband corner frequency (in Hz). + /// Maximum allowed passband ripple. + /// Minimum required stopband attenuation. + /// Online IIR filter + 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); + } + + /// + /// Generates an online IIR notch Butterworth filter. + /// + /// Filter central frequency. + /// Quality factor. + /// Maximum allowed passband ripple. + /// Minimum required stopband attenuation. + /// Online IIR filter + 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); + } + } +} diff --git a/src/Filtering/IIR/OnlineIirFilter.cs b/src/Filtering/IIR/OnlineIirFilter.cs index ae1d642b..dc30fc5d 100644 --- a/src/Filtering/IIR/OnlineIirFilter.cs +++ b/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. /// /// - /// 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) /// public class OnlineIirFilter : OnlineFilter { @@ -54,6 +54,9 @@ public class OnlineIirFilter : OnlineFilter /// /// Infinite Impulse Response (IIR) Filter. /// + /// Vector of IIR coefficients in the following form: [b0 b1 b2 ... a0 a1 a2 ...] + /// + /// public OnlineIirFilter(double[] coefficients) { if (null == coefficients) @@ -81,6 +84,50 @@ public OnlineIirFilter(double[] coefficients) _bufferY = new double[size]; } + /// + /// Infinite Impulse Response (IIR) Filter + /// + /// Numerator coefficients [b0 b1 ...] + /// Denominator coefficients [a0 a1 ...] + /// + /// + 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]; + } + /// /// Process a single sample. /// diff --git a/src/Filtering/Properties/Resources.Designer.cs b/src/Filtering/Properties/Resources.Designer.cs index 62a1ee58..4f908846 100644 --- a/src/Filtering/Properties/Resources.Designer.cs +++ b/src/Filtering/Properties/Resources.Designer.cs @@ -1,7 +1,7 @@ //------------------------------------------------------------------------------ // // This code was generated by a tool. -// Runtime Version:4.0.30319.34014 +// Runtime Version:4.0.30319.42000 // // Changes to this file may cause incorrect behavior and will be lost if // the code is regenerated. @@ -21,7 +21,7 @@ namespace MathNet.Filtering.Properties { // class via a tool like ResGen or Visual Studio. // To add or remove a member, edit your .ResX file then rerun ResGen // with the /str option, or rebuild your VS project. - [global::System.CodeDom.Compiler.GeneratedCodeAttribute("System.Resources.Tools.StronglyTypedResourceBuilder", "4.0.0.0")] + [global::System.CodeDom.Compiler.GeneratedCodeAttribute("System.Resources.Tools.StronglyTypedResourceBuilder", "17.0.0.0")] [global::System.Diagnostics.DebuggerNonUserCodeAttribute()] [global::System.Runtime.CompilerServices.CompilerGeneratedAttribute()] internal class Resources { @@ -71,6 +71,23 @@ internal Resources() { } } + /// + /// Looks up a localized string similar to Numerator and denominator vectors must contain the same number of elements.. + /// + internal static string ArgumentCoefficientLengthMismatch { + get { + return ResourceManager.GetString("ArgumentCoefficientLengthMismatch", resourceCulture); + } + } + /// + /// Looks up a localized string similar to Coefficient vector must not be empty.. + /// + internal static string ArgumentEmptyCoefficientVector { + get { + return ResourceManager.GetString("ArgumentEmptyCoefficientVector", resourceCulture); + } + } + /// /// Looks up a localized string similar to Even number of coefficients required.. /// diff --git a/src/Filtering/Properties/Resources.resx b/src/Filtering/Properties/Resources.resx index 569e1856..cee14718 100644 --- a/src/Filtering/Properties/Resources.resx +++ b/src/Filtering/Properties/Resources.resx @@ -112,10 +112,10 @@ 2.0 - System.Resources.ResXResourceReader, System.Windows.Forms, Version=2.0.0.0, Culture=neutral, PublicKeyToken=b77a5c561934e089 + System.Resources.ResXResourceReader, System.Windows.Forms, Version=4.0.0.0, Culture=neutral, PublicKeyToken=b77a5c561934e089 - System.Resources.ResXResourceWriter, System.Windows.Forms, Version=2.0.0.0, Culture=neutral, PublicKeyToken=b77a5c561934e089 + System.Resources.ResXResourceWriter, System.Windows.Forms, Version=4.0.0.0, Culture=neutral, PublicKeyToken=b77a5c561934e089 Even number of coefficients required. @@ -153,4 +153,10 @@ Measurement covariance should be a square matrix with rows/cols equal to number of measurements. + + Numerator and denominator vectors must contain the same number of elements. + + + Coefficient vector must not be empty. + \ No newline at end of file From 746489f3f00f76f7fde88ebc5be9df94371cfb78 Mon Sep 17 00:00:00 2001 From: StefanoPietrosanti Date: Wed, 13 Oct 2021 12:06:03 +0100 Subject: [PATCH 2/2] Fixed issue in Butterworth.IirCoefficients.cs Fixed issue in IirCoefficients.Notch() where wc2 frequency was not pre-warped and wc1 was pre-warped twice. Removed unneeded System.Linq inclusion in MathFunctions.cs --- src/Filtering/Butterworth/IirCoefficients.cs | 2 +- src/Filtering/Helpers/MathFunctions.cs | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Filtering/Butterworth/IirCoefficients.cs b/src/Filtering/Butterworth/IirCoefficients.cs index 45d18fa4..bf45bd08 100644 --- a/src/Filtering/Butterworth/IirCoefficients.cs +++ b/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); diff --git a/src/Filtering/Helpers/MathFunctions.cs b/src/Filtering/Helpers/MathFunctions.cs index 99f5afef..f4fb37e4 100644 --- a/src/Filtering/Helpers/MathFunctions.cs +++ b/src/Filtering/Helpers/MathFunctions.cs @@ -28,7 +28,6 @@ // using System; -using System.Linq; using System.Numerics; using MathNet.Numerics;