Add minimum requirement functionvfor HP and LP.

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 <typename T>
 class Butterworth : public DigitalFilter<T> {
 public:
-    T PI = static_cast<T>(M_PI);
+    static T PI;
 
 public:
     enum class Type {
@@ -24,18 +25,21 @@ public:
         BandReject
     };
 
-    // public:
+public:
-    //     static double minimumRequiredFrequency(...);
+    // http://www.matheonics.com/Tutorials/Butterworth.html#Paragraph_3.2
+    // https://www.mathworks.com/help/signal/ref/buttord.html#d120e11079
+    static std::pair<int, T> 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<T> generateAnalogPole(int k, T fpw);

include/Butterworth.tpp

@@ -3,6 +3,33 @@
 
 namespace fratio {
 
+template <typename T>
+T Butterworth<T>::PI = static_cast<T>(M_PI);
+
+template <typename T>
+std::pair<int, T> Butterworth<T>::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<int>(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<int, T>(order, T(2) * std::atan(ctf) / PI);
+}
+
 template <typename T>
 Butterworth<T>::Butterworth(Type type)
     : m_type(type)
@@ -30,30 +57,30 @@ void Butterworth<T>::setFilterParameters(int order, T fc, T fs)
 }
 
 template <typename T>
-void Butterworth<T>::setFilterParameters(int order, T f0, T fLimit, T fs)
+void Butterworth<T>::setFilterParameters(int order, T fLower, T fUpper, T fs)
 {
-    initialize(order, f0, fLimit, fs);
+    initialize(order, fLower, fUpper, fs);
 }
 
 template <typename T>
-void Butterworth<T>::initialize(int order, T f0, T fLimit, T fs)
+void Butterworth<T>::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<T>::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();
 }

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<float>(0.081038494957764) - butterRequirement.second), std::numeric_limits<float>::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<float>(0.296655824107340) - butterRequirement.second), std::numeric_limits<float>::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<double>::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<double>::epsilon() * 10);
+}
+
 BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_LP_FILTER_FLOAT, System<float>)
 {
     auto bf = fratio::Butterworthf(order, fc, fs);

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);
 }