Split Bilinear tranform and update type checks.

include/BilinearTransform.h

@@ -0,0 +1,48 @@
+#pragma once
+
+#include <vector>
+
+namespace fratio {
+
+template <typename T>
+struct BilinearTransform {
+    using SubType = internal::complex_sub_type_t<T>;
+    static_assert(std::is_floating_point<SubType>::value, "This struct can only accept floating point types (real and complex).");
+
+    static void SToZ(SubType fs, const T& sPlanePole, T& zPlanePole);
+    static void SToZ(SubType fs, const std::vector<T>& sPlanePoles, std::vector<T>& zPlanePoles);
+    static void ZToS(SubType fs, const T& zPlanePole, T& sPlanePole);
+    static void ZToS(SubType fs, const std::vector<T>& zPlanePoles, std::vector<T>& sPlanePoles);
+};
+
+template <typename T>
+void BilinearTransform<T>::SToZ(SubType fs, const T& sPlanePole, T& zPlanePole)
+{
+    T scalePole = sPlanePole / (2 * fs);
+    zPlanePole = (T(1) + scalePole) / (T(1) - scalePole);
+}
+
+template <typename T>
+void BilinearTransform<T>::SToZ(SubType fs, const std::vector<T>& sPlanePoles, std::vector<T>& zPlanePoles)
+{
+    assert(sPlanePoles.size() == zPlanePoles.size());
+    for (size_t k = 0; k < sPlanePoles.size(); ++k)
+        SToZ(fs, sPlanePoles[k], zPlanePoles[k]);
+}
+
+template <typename T>
+void BilinearTransform<T>::ZToS(SubType fs, const T& zPlanePole, T& sPlanePole)
+{
+    T invPole = T(1) / zPlanePole;
+    sPlanePole = 2 * fs * (T(1) - invPole) / (T(1) + invPole);
+}
+
+template <typename T>
+void BilinearTransform<T>::ZToS(SubType fs, const std::vector<T>& zPlanePoles, std::vector<T>& sPlanePoles)
+{
+    assert(zPlanePoles.size() == sPlanePoles.size());
+    for (size_t k = 0; k < sPlanePoles.size(); ++k)
+        ZToS(fs, zPlanePoles[k], sPlanePoles[k]);
+}
+
+} // namespace fratio

include/Butterworth.h

@@ -11,7 +11,8 @@ template <typename T>
 class Butterworth : public GenericFilter<T> {
 public:
     enum class Type {
-        LowPass
+        LowPass,
+        HighPass
     };
 
     // public:
@@ -26,6 +27,7 @@ private:
     void initialize(size_t order, T fc, T fs);
     void computeDigitalRep();
     void updateCoeffSize();
+    void transformFilter();
 
 private:
     Type m_type;

include/Butterworth.inl

@@ -1,3 +1,4 @@
+#include "BilinearTransform.h"
 #include "polynome_functions.h"
 #include <cmath>
 #include <sstream>
@@ -58,8 +59,7 @@ void Butterworth<T>::computeDigitalRep()
     std::complex<T> scalePole;
     for (size_t k = 1; k <= m_order; ++k) {
         scalePole = scaleFactor * std::complex<T>(-std::sin(thetaK(k)), std::cos(thetaK(k)));
-        scalePole /= 2 * m_fs;
-        m_poles[k - 1] = (T(1) + scalePole) / (T(1) - scalePole);
+        BilinearTransform<std::complex<T>>::SToZ(m_fs, scalePole, m_poles[k - 1]);
     }
 
     std::vector<std::complex<T>> numPoles(m_order, std::complex<T>(-1));
@@ -89,4 +89,17 @@ void Butterworth<T>::updateCoeffSize()
     resetFilter();
 }
 
+template <typename T>
+void Butterworth<T>::transformFilter()
+{
+    switch (m_type) {
+    case Type::HighPass:
+        break;
+
+    case Type::LowPass:
+    default:
+        break;
+    }
+}
+
 } // namespace fratio

include/GenericFilter.h

@@ -6,11 +6,10 @@
 
 namespace fratio {
 
-template <typename T, typename = std::enable_if_t<std::is_floating_point<T>::value && !std::is_const<T>::value>>
-class GenericFilter;
-
 template <typename T>
-class GenericFilter<T> {
+class GenericFilter {
+    static_assert(std::is_floating_point<T>::value && !std::is_const<T>::value, "Only accept non-complex floating point types.");
+
 public:
     T stepFilter(T data);
     std::vector<T> filter(const std::vector<T>& data);

include/fratio.h

@@ -7,6 +7,7 @@
 
 namespace fratio {
 
+// Filters
 using DigitalFilterf = DigitalFilter<float>;
 using DigitalFilterd = DigitalFilter<double>;
 using MovingAveragef = MovingAverage<float>;
@@ -14,6 +15,7 @@ using MovingAveraged = MovingAverage<double>;
 using Butterworthf = Butterworth<float>;
 using Butterworthd = Butterworth<double>;
 
+// Polynome helper functions
 using VietaAlgof = VietaAlgo<float>;
 using VietaAlgod = VietaAlgo<double>;
 using VietaAlgoi = VietaAlgo<int>;
@@ -21,4 +23,10 @@ using VietaAlgocf = VietaAlgo<std::complex<float>>;
 using VietaAlgocd = VietaAlgo<std::complex<double>>;
 using VietaAlgoci = VietaAlgo<std::complex<int>>;
 
+// Bilinear transformation functions
+using BilinearTransformf = BilinearTransform<float>;
+using BilinearTransformd = BilinearTransform<double>;
+using BilinearTransformcf = BilinearTransform<std::complex<float>>;
+using BilinearTransformcd = BilinearTransform<std::complex<double>>;
+
 } // namespace fratio

include/polynome_functions.h

@@ -6,11 +6,10 @@
 
 namespace fratio {
 
-template <typename T, typename = std::enable_if_t<std::is_arithmetic<T>::value || is_complex_t<T>::value>>
-struct VietaAlgo;
-
 template <typename T>
-struct VietaAlgo<T> {
+struct VietaAlgo {
+    static_assert(std::is_arithmetic<internal::complex_sub_type_t<T>>::value, "This struct can only accept arithmetic types or complex.");
+
     // Vieta's computation: https://en.wikipedia.org/wiki/Vieta%27s_formulas
     static std::vector<T> polyCoeffFromRoot(const std::vector<T>& poles);
 };

include/type_checks.h

@@ -4,12 +4,29 @@
 
 namespace fratio {
 
-template <typename T>
-struct is_complex_t : public std::false_type {
-};
+namespace internal {
 
-template <typename T>
-struct is_complex_t<std::complex<T>> : public std::true_type {
-};
+    // Based on https://stackoverflow.com/a/41449096/9774052
+    template <typename T>
+    std::true_type is_tt_impl(std::complex<T>);
+    std::false_type is_tt_impl(...);
+
+    template <typename T>
+    using is_complex = decltype(is_tt_impl(std::declval<T>()));
+
+    template <typename T, bool B>
+    struct sub_type {
+        using type = T;
+    };
+
+    template <typename T>
+    struct sub_type<T, true> {
+        using type = typename T::value_type;
+    };
+
+    template <typename T>
+    using complex_sub_type_t = typename sub_type<T, is_complex<T>::value>::type;
+
+} // namespace internal
 
 } // namespace fratio

tests/ButterworthFilterTests.cpp

@@ -8,7 +8,7 @@ DISABLE_CONVERSION_WARNING_BEGIN
 
 template <typename T>
 struct System {
-    std::vector<T> data = { 1., 2., 3., 4., 5., 6., 7., 8. };
+    std::vector<T> data = { 1, 2, 3, 4, 5, 6, 7, 8 };
     size_t order = 5;
     T fc = 10;
     T fs = 100;
@@ -88,7 +88,7 @@ BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_FILTER_DOUBLE, System<double>)
 
     bf2.setFilterParameters(order, fc, fs);
     filteredData.resize(0);
-    for (float d : data)
+    for (double d : data)
         filteredData.push_back(bf2.stepFilter(d));
 
     for (size_t i = 0; i < filteredData.size(); ++i)

tests/DigitalFilterTests.cpp

@@ -8,8 +8,8 @@ DISABLE_CONVERSION_WARNING_BEGIN
 
 template <typename T>
 struct System {
-    std::vector<T> data = { 1., 2., 3., 4. };
-    std::vector<T> aCoeff = { 1., -0.99993717 };
+    std::vector<T> data = { 1, 2, 3, 4 };
+    std::vector<T> aCoeff = { 1, -0.99993717 };
     std::vector<T> bCoeff = { 0.99996859, -0.99996859 };
     std::vector<T> results = { 0.99996859, 1.999874351973491, 2.999717289867956, 3.999497407630634 };
 };

tests/GenericFilterTests.cpp

@@ -5,12 +5,12 @@
 
 BOOST_AUTO_TEST_CASE(FilterThrows)
 {
-    BOOST_REQUIRE_THROW(fratio::DigitalFilterd({}, { 1., 2. }), std::runtime_error);
-    BOOST_REQUIRE_THROW(fratio::DigitalFilterd({ 1., 2. }, {}), std::runtime_error);
-    BOOST_REQUIRE_THROW(fratio::DigitalFilterd({ 0. }, { 1. }), std::invalid_argument);
+    BOOST_REQUIRE_THROW(fratio::DigitalFilterd({}, { 1, 2 }), std::runtime_error);
+    BOOST_REQUIRE_THROW(fratio::DigitalFilterd({ 1, 2 }, {}), std::runtime_error);
+    BOOST_REQUIRE_THROW(fratio::DigitalFilterd({ 0 }, { 1 }), std::invalid_argument);
 
     auto df = fratio::DigitalFilterd();
-    BOOST_REQUIRE_THROW(df.setCoeff({}, { 1., 2. }), std::runtime_error);
-    BOOST_REQUIRE_THROW(df.setCoeff({ 1., 2. }, {}), std::runtime_error);
-    BOOST_REQUIRE_THROW(df.setCoeff({ 0. }, { 1. }), std::invalid_argument);
+    BOOST_REQUIRE_THROW(df.setCoeff({}, { 1, 2 }), std::runtime_error);
+    BOOST_REQUIRE_THROW(df.setCoeff({ 1, 2 }, {}), std::runtime_error);
+    BOOST_REQUIRE_THROW(df.setCoeff({ 0 }, { 1 }), std::invalid_argument);
 }

tests/MovingAverageFilterTests.cpp

@@ -5,7 +5,7 @@
 
 template <typename T>
 struct System {
-    std::vector<T> data = { 1., 2., 3., 4., 5., 6. };
+    std::vector<T> data = { 1, 2, 3, 4, 5, 6 };
     size_t windowSize = 4;
     std::vector<T> results = { 0.25, 0.75, 1.5, 2.5, 3.5, 4.5 };
 };

tests/polynome_functions_tests.cpp

@@ -19,13 +19,13 @@ DISABLE_CONVERSION_WARNING_BEGIN
 template <typename T>
 struct SystemFloat {
     std::vector<T> data = { 0.32, -0.0518, 41.4, 0.89 };
-    std::vector<T> results = { 1., -42.558199999999999, 48.171601999999993, -9.181098159999999, -0.610759296 };
+    std::vector<T> results = { 1, -42.558199999999999, 48.171601999999993, -9.181098159999999, -0.610759296 };
 };
 
 template <typename T>
 struct SystemCFloat {
     std::vector<std::complex<T>> data = { { 1.35, 0.2 }, { -1.5, 4.45 }, { 12.8, -3.36 }, { 5.156, 2.12 } };
-    std::vector<std::complex<T>> results = { { 1., 0. }, { -17.806, -3.41 }, { 73.2776, 99.20074 }, { 101.857496, -444.634694 }, { -269.1458768, 388.7308864 } };
+    std::vector<std::complex<T>> results = { { 1, 0 }, { -17.806, -3.41 }, { 73.2776, 99.20074 }, { 101.857496, -444.634694 }, { -269.1458768, 388.7308864 } };
 };
 
 DISABLE_CONVERSION_WARNING_END