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/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/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/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; 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