From 655db4a39a2de3a08d539b284e42c761f68d5e6b Mon Sep 17 00:00:00 2001 From: Vincent Samy Date: Fri, 4 Jan 2019 15:13:08 +0900 Subject: [PATCH] Add minimum requirement functionvfor HP and LP. --- include/Butterworth.h | 14 +++++++---- include/Butterworth.tpp | 43 ++++++++++++++++++++++++++------ tests/ButterworthFilterTests.cpp | 26 +++++++++++++++++++ tests/GenericFilterTests.cpp | 2 +- 4 files changed, 71 insertions(+), 14 deletions(-) diff --git a/include/Butterworth.h b/include/Butterworth.h index 058be42..e646766 100644 --- a/include/Butterworth.h +++ b/include/Butterworth.h @@ -11,10 +11,11 @@ namespace fratio { // https://www.dsprelated.com/showarticle/1135.php // https://www.dsprelated.com/showarticle/1128.php // https://www.dsprelated.com/showarticle/1131.php +// https://www.mathworks.com/help/signal/ref/butter.html template class Butterworth : public DigitalFilter { public: - T PI = static_cast(M_PI); + static T PI; public: enum class Type { @@ -24,18 +25,21 @@ public: BandReject }; - // public: - // static double minimumRequiredFrequency(...); +public: + // http://www.matheonics.com/Tutorials/Butterworth.html#Paragraph_3.2 + // https://www.mathworks.com/help/signal/ref/buttord.html#d120e11079 + static std::pair findMinimumButter(T wPass, T wStop, T APass, T AStop); // Only for low and high pass + public: Butterworth(Type type = Type::LowPass); Butterworth(int order, T fc, T fs, Type type = Type::LowPass); Butterworth(int order, T fLower, T fUpper, T fs, Type type = Type::BandPass); void setFilterParameters(int order, T fc, T fs); - void setFilterParameters(int order, T f0, T fLimit, T fs); + void setFilterParameters(int order, T fLower, T fUpper, T fs); private: - void initialize(int order, T f0, T fLimit, T fs); // f0 = fc for LowPass/HighPass filter + void initialize(int order, T f1, T f2, T fs); // f1 = fc for LowPass/HighPass filter void computeDigitalRep(T fc); void computeBandDigitalRep(T fLower, T fUpper); std::complex generateAnalogPole(int k, T fpw); diff --git a/include/Butterworth.tpp b/include/Butterworth.tpp index ba9c1f6..f6fef49 100644 --- a/include/Butterworth.tpp +++ b/include/Butterworth.tpp @@ -3,6 +3,33 @@ namespace fratio { +template +T Butterworth::PI = static_cast(M_PI); + +template +std::pair Butterworth::findMinimumButter(T wPass, T wStop, T APass, T AStop) +{ + T num = std::log10((std::pow(T(10), T(0.1) * std::abs(AStop)) - 1) / (std::pow(T(10), T(0.1) * std::abs(APass)) - 1)); + // pre-warp + T fwPass = std::tan(0.5 * PI * wPass); + T fwStop = std::tan(0.5 * PI * wStop); + T w; + if (wPass < wStop) + w = std::abs(fwStop / fwPass); + else + w = std::abs(fwPass / fwStop); + T denum = T(2) * std::log10(w); + + int order = static_cast(std::ceil(num / denum)); + T ctf = w / std::pow(std::pow(T(10), T(0.1) * std::abs(AStop)) - 1, T(1) / T(2 * order)); + if (wPass < wStop) + ctf *= fwPass; + else + ctf = fwPass / ctf; + + return std::pair(order, T(2) * std::atan(ctf) / PI); +} + template Butterworth::Butterworth(Type type) : m_type(type) @@ -30,30 +57,30 @@ void Butterworth::setFilterParameters(int order, T fc, T fs) } template -void Butterworth::setFilterParameters(int order, T f0, T fLimit, T fs) +void Butterworth::setFilterParameters(int order, T fLower, T fUpper, T fs) { - initialize(order, f0, fLimit, fs); + initialize(order, fLower, fUpper, fs); } template -void Butterworth::initialize(int order, T f0, T fLimit, T fs) +void Butterworth::initialize(int order, T f1, T f2, T fs) { if (order <= 0) { m_status = FilterStatus::BAD_ORDER_SIZE; return; } - if (f0 <= 0 || fs <= 0) { + if (f1 <= 0 || fs <= 0) { m_status = FilterStatus::BAD_FREQUENCY_VALUE; return; } - if ((m_type == Type::BandPass || m_type == Type::BandReject) && f0 >= fLimit) { + if ((m_type == Type::BandPass || m_type == Type::BandReject) && f1 >= f2) { m_status = FilterStatus::BAD_BAND_FREQUENCY; return; } - if ((m_type == Type::LowPass || m_type == Type::HighPass) && f0 > fs / 2.) { + if ((m_type == Type::LowPass || m_type == Type::HighPass) && f1 > fs / 2.) { m_status = FilterStatus::BAD_CUTOFF_FREQUENCY; return; } @@ -61,9 +88,9 @@ void Butterworth::initialize(int order, T f0, T fLimit, T fs) m_order = order; m_fs = fs; if (m_type == Type::LowPass || m_type == Type::HighPass) - computeDigitalRep(f0); + computeDigitalRep(f1); else - computeBandDigitalRep(f0, fLimit); // For band-like filters + computeBandDigitalRep(f1, f2); // For band-like filters resetFilter(); } diff --git a/tests/ButterworthFilterTests.cpp b/tests/ButterworthFilterTests.cpp index c675101..630aa1a 100644 --- a/tests/ButterworthFilterTests.cpp +++ b/tests/ButterworthFilterTests.cpp @@ -37,6 +37,32 @@ struct System { DISABLE_CONVERSION_WARNING_END +BOOST_AUTO_TEST_CASE(FIND_BUTTERWORTH_LP_HP_FLOAT) +{ + // LP + auto butterRequirement = fratio::Butterworthf::findMinimumButter(40.f / 500.f, 150.f / 500.f, 3.f, 60.f); + BOOST_REQUIRE_EQUAL(5, butterRequirement.first); + BOOST_REQUIRE_SMALL(std::abs(static_cast(0.081038494957764) - butterRequirement.second), std::numeric_limits::epsilon() * 10); + + // HP + butterRequirement = fratio::Butterworthf::findMinimumButter(150.f / 500.f, 40.f / 500.f, 3.f, 60.f); + BOOST_REQUIRE_EQUAL(5, butterRequirement.first); + BOOST_REQUIRE_SMALL(std::abs(static_cast(0.296655824107340) - butterRequirement.second), std::numeric_limits::epsilon() * 10); +} + +BOOST_AUTO_TEST_CASE(FIND_BUTTERWORTH_LP_HP_DOUBLE) +{ + // LP + auto butterRequirement = fratio::Butterworthd::findMinimumButter(40. / 500., 150. / 500., 3., 60.); + BOOST_REQUIRE_EQUAL(5, butterRequirement.first); + BOOST_REQUIRE_SMALL(std::abs(0.081038494957764 - butterRequirement.second), std::numeric_limits::epsilon() * 10); + + // HP + butterRequirement = fratio::Butterworthd::findMinimumButter(150. / 500., 40. / 500., 3., 60.); + BOOST_REQUIRE_EQUAL(5, butterRequirement.first); + BOOST_REQUIRE_SMALL(std::abs(0.296655824107340 - butterRequirement.second), std::numeric_limits::epsilon() * 10); +} + BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_LP_FILTER_FLOAT, System) { auto bf = fratio::Butterworthf(order, fc, fs); diff --git a/tests/GenericFilterTests.cpp b/tests/GenericFilterTests.cpp index f14ccb6..4b94ca3 100644 --- a/tests/GenericFilterTests.cpp +++ b/tests/GenericFilterTests.cpp @@ -22,6 +22,6 @@ BOOST_AUTO_TEST_CASE(FILTER_FAILURES) BOOST_REQUIRE(mad.status() == fratio::FilterStatus::BAD_ORDER_SIZE); auto bfd = fratio::Butterworthd(0, 10, 100); BOOST_REQUIRE(bfd.status() == fratio::FilterStatus::BAD_ORDER_SIZE); - bfd = fratio::Butterworthd(5, 5, 6, 100); + bfd = fratio::Butterworthd(5, 6, 5, 100); BOOST_REQUIRE(bfd.status() == fratio::FilterStatus::BAD_BAND_FREQUENCY); }