Update CMake and tests.

Closes #1
Closes #2

.clang-format

@@ -13,11 +13,12 @@ AllowShortBlocksOnASingleLine: true
 AllowShortCaseLabelsOnASingleLine: true
 AllowShortFunctionsOnASingleLine: All
 AllowShortIfStatementsOnASingleLine: false
+AllowShortLambdasOnASingleLine: All
 AllowShortLoopsOnASingleLine: false
 AlwaysBreakAfterDefinitionReturnType: None
 AlwaysBreakAfterReturnType: None
 AlwaysBreakBeforeMultilineStrings: false
-AlwaysBreakTemplateDeclarations: MultiLine
+AlwaysBreakTemplateDeclarations: No
 BinPackArguments: true
 BinPackParameters: true
 BraceWrapping:   

CMakeLists.txt

@@ -26,19 +26,18 @@
 # either expressed or implied, of the FreeBSD Project.
 
 # Version minimum
-cmake_minimum_required(VERSION 3.1.3)
+cmake_minimum_required(VERSION 3.8.2)
 
-set(PROJECT_NAME diFi++)
+set(PROJECT_NAME difi++)
 set(PROJECT_DESCRIPTION "Filter using rational transfer function")
 set(PROJECT_URL "https://github.com/vsamy/DiFipp")
 SET(PROJECT_DEBUG_POSTFIX "_d")
 set(INSTALL_GENERATED_HEADERS OFF)
 
 set(DOXYGEN_USE_MATHJAX "NO")
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 
 option(BUILD_TESTING "Disable unit tests." ON)
-option(BUILD_TEST_STATIC_BOOST "Build unit tests with static boost libraries" OFF)
 option(THROW_ON_CONTRACT_VIOLATION "Throw an error when program fails." ON)
 option(TERMINATE_ON_CONTRACT_VIOLATION "Terminate program when an error occurs. (Default)" OFF)
 option(UNENFORCED_ON_CONTRACT_VIOLATION "Do not perform any check." OFF)
@@ -71,6 +70,3 @@ add_subdirectory(include)
 if(${BUILD_TESTING})
     add_subdirectory(tests)
 endif()
-
-setup_project_finalize()
-setup_project_package_finalize()

README.md

@@ -7,10 +7,16 @@ DiFi++ is a small c++ header-only library  for **DI**gital **FI**lters based on
 
 The implementation is based on well written article from Neil Robertson.
 Please check out the followings
- * https://www.dsprelated.com/showarticle/1119.php
- * https://www.dsprelated.com/showarticle/1135.php
- * https://www.dsprelated.com/showarticle/1128.php
- * https://www.dsprelated.com/showarticle/1131.php
+
+* [Butterworth filter](https://www.dsprelated.com/showarticle/1119.php)
+
+* [Highpass filters](https://www.dsprelated.com/showarticle/1135.php)
+
+* [Bandpass filters](https://www.dsprelated.com/showarticle/1128.php)
+
+* [Band-reject filters](https://www.dsprelated.com/showarticle/1131.php)
+
+and the differentiators from [Pavel Holoborodko](http://www.holoborodko.com/pavel/)
 
 The library has been tested against Matlab results.
 
@@ -32,4 +38,5 @@ make install
 
 Note
 -----
+
 The method used is close but somewhat different from Matlab methods and Butterworth band-reject has quite different results (precision of 1e-8).

cmake

@@ -1 +1 @@
-Subproject commit 192812cc2c5d8aaffbc8536d0b0a2a7c11dad53b
+Subproject commit f4997a81cebfa2dfb69733bc088c55688965dfe8

include/BaseFilter.h

@@ -0,0 +1,137 @@
+// Copyright (c) 2019, Vincent SAMY
+// All rights reserved.
+
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+
+// 1. Redistributions of source code must retain the above copyright notice,
+//    this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+// The views and conclusions contained in the software and documentation are those
+// of the authors and should not be interpreted as representing official policies,
+// either expressed or implied, of the FreeBSD Project.
+
+#pragma once
+
+#include "gsl/gsl_assert.h"
+#include "type_checks.h"
+#include "typedefs.h"
+#include <string>
+
+namespace difi {
+
+// TODO: noexcept(Function of gsl variable)
+// TODO: constructor with universal refs
+// TODO: setCoeffs with universal refs
+
+/*! \brief Low-level filter.
+ * 
+ * It creates the basic and common functions of all linear filter that can written as a digital filter.
+ * This class can not be instantiated directly.
+ * 
+ * \warning In Debug mode, all functions may throw if a filter is badly initialized.
+ * This not the case in Realese mode.
+ * 
+ * \tparam T Floating type.
+ */
+template <typename T, typename Derived>
+class BaseFilter {
+    static_assert(std::is_floating_point<T>::value && !std::is_const<T>::value, "Only accept non-complex floating point types.");
+    friend Derived;
+
+public:
+    /*! \brief Reset the data and filtered data. */
+    void resetFilter() noexcept { derived().resetFilter(); };
+
+    /*!< \brief Return the filter type */
+    FilterType type() const noexcept { return m_type; }
+    /*! \brief Get how far back is the filtered value.
+     * The filtered value is computed at time \f$tc=t - center()*\Delta T\f$ where \f$t\f$ is the current time
+     * \f$\Delta T\f$ is your sampling time. 
+     * \note For FilterType::Backward filter, the function returns 0.
+     */
+    Eigen::Index center() const noexcept { return (m_type == FilterType::Backward ? 0 : ((m_bCoeff.size() - 1) / 2)); }
+    /*! \brief Get digital filter coefficients.
+     * 
+     * It will automatically resize the given vectors.
+     * \param[out] aCoeff Denominator coefficients of the filter in decreasing order.
+     * \param[out] bCoeff Numerator coefficients of the filter in decreasing order.
+     */
+    void getCoeffs(vectX_t<T>& aCoeff, vectX_t<T>& bCoeff) const noexcept;
+    /*! \brief Return coefficients of the denominator polynome. */
+    const vectX_t<T>& aCoeff() const noexcept { return m_aCoeff; }
+    /*! \brief Return coefficients of the numerator polynome. */
+    const vectX_t<T>& bCoeff() const noexcept { return m_bCoeff; }
+    /*! \brief Return the order the denominator polynome order of the filter. */
+    Eigen::Index aOrder() const noexcept { return m_aCoeff.size(); }
+    /*! \brief Return the order the numerator polynome order of the filter. */
+    Eigen::Index bOrder() const noexcept { return m_bCoeff.size(); }
+    /*! \brief Return the initialization state of the filter0 */
+    bool isInitialized() const noexcept { return m_isInitialized; }
+    /*! \brief Set type of filter (one-sided or centered)
+     * 
+     * \param type The filter type.
+     * \warning bCoeff must be set before.
+     */
+    void setType(FilterType type);
+    /*! \brief Set the new coefficients of the filters.
+     *
+     * It awaits a universal reference.
+     * \param aCoeff Denominator coefficients of the filter in decreasing order.
+     * \param bCoeff Numerator coefficients of the filter in decreasing order.
+     */
+    void setCoeffs(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff);
+
+protected:
+    /*! \brief Normalized the filter coefficients such that aCoeff(0) = 1. */
+    void normalizeCoeffs();
+    /*! \brief Check for bad coefficients.
+     * 
+     * Set the filter status to ready is everything is fine.
+     * \param aCoeff Denominator coefficients of the filter.
+     * \param bCoeff Numerator coefficients of the filter.
+     * \return True if the filter status is set on READY.
+     */
+    bool checkCoeffs(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff);
+
+private:
+    /*! \brief Default uninitialized constructor. */
+    BaseFilter() = default;
+    /*! \brief Constructor.
+     * \param aCoeff Denominator coefficients of the filter in decreasing order.
+     * \param bCoeff Numerator coefficients of the filter in decreasing order.
+     * \param type Type of the filter.
+     */
+    BaseFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, FilterType type = FilterType::Backward);
+    /*! \brief Default destructor. */
+    virtual ~BaseFilter() = default;
+
+    Derived& derived() noexcept { return *static_cast<Derived*>(this); }
+    const Derived& derived() const noexcept { return *static_cast<const Derived*>(this); }
+
+private:
+    FilterType m_type; /*!< Type of filter. Default is FilterType::Backward. */
+    bool m_isInitialized = false; /*!< Initialization state of the filter. Default is false */
+    vectX_t<T> m_aCoeff; /*!< Denominator coefficients of the filter */
+    vectX_t<T> m_bCoeff; /*!< Numerator coefficients of the filter */
+    vectX_t<T> m_filteredData; /*!< Last set of filtered data */
+    vectX_t<T> m_rawData; /*!< Last set of non-filtered data */
+};
+
+} // namespace difi
+
+#include "BaseFilter.tpp"

include/BaseFilter.tpp

@@ -0,0 +1,93 @@
+// Copyright (c) 2019, Vincent SAMY
+// All rights reserved.
+
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+
+// 1. Redistributions of source code must retain the above copyright notice,
+//    this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+// The views and conclusions contained in the software and documentation are those
+// of the authors and should not be interpreted as representing official policies,
+// either expressed or implied, of the FreeBSD Project.
+
+#include <limits>
+
+namespace difi {
+
+// Public functions
+
+template <typename T, typename Derived>
+void BaseFilter<T, Derived>::setType(FilterType type)
+{
+    Expects(type == FilterType::Centered ? m_bCoeff.size() > 2 && m_bCoeff.size() % 2 == 1 : true);
+    m_type = type;
+}
+
+template <typename T, typename Derived>
+void BaseFilter<T, Derived>::setCoeffs(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff)
+{
+    Expects(checkCoeffs(aCoeff, bCoeff));
+    m_aCoeff = aCoeff;
+    m_bCoeff = bCoeff;
+    normalizeCoeffs();
+    resetFilter();
+    m_isInitialized = true;
+}
+
+template <typename T, typename Derived>
+void BaseFilter<T, Derived>::getCoeffs(vectX_t<T>& aCoeff, vectX_t<T>& bCoeff) const noexcept
+{
+    aCoeff = m_aCoeff;
+    bCoeff = m_bCoeff;
+}
+
+// Protected functions
+
+template <typename T, typename Derived>
+BaseFilter<T, Derived>::BaseFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, FilterType type)
+    : m_type(type)
+    , m_aCoeff(aCoeff)
+    , m_bCoeff(bCoeff)
+    , m_filteredData(aCoeff.size())
+    , m_rawData(bCoeff.size())
+{
+    Expects(checkCoeffs(aCoeff, bCoeff));
+    normalizeCoeffs();
+    resetFilter();
+    m_isInitialized = true;
+}
+
+template <typename T, typename Derived>
+void BaseFilter<T, Derived>::normalizeCoeffs()
+{
+    T a0 = m_aCoeff(0);
+    if (std::abs(a0 - T(1)) < std::numeric_limits<T>::epsilon())
+        return;
+
+    m_aCoeff /= a0;
+    m_bCoeff /= a0;
+}
+
+template <typename T, typename Derived>
+bool BaseFilter<T, Derived>::checkCoeffs(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff)
+{
+    bool centering = (m_type == FilterType::Centered ? (bCoeff.size() % 2 == 1) : true);
+    return aCoeff.size() > 0 && std::abs(aCoeff[0]) > std::numeric_limits<T>::epsilon() && bCoeff.size() > 0 && centering;
+}
+
+} // namespace difi

include/Butterworth.h

@@ -60,7 +60,7 @@ public:
     /*! \brief Function to help you design a Butterworth filter.
      * 
      * It finds optimal values of the order and cut-off frequency.
-     * \warning Works only for low-pass and high-pass filters.
+     * \note Works only for low-pass and high-pass filters.
      * \see http://www.matheonics.com/Tutorials/Butterworth.html#Paragraph_3.2
      * \see https://www.mathworks.com/help/signal/ref/buttord.html#d120e11079
      * \param wPass Pass band edge.

include/Butterworth.tpp

@@ -131,7 +131,7 @@ void Butterworth<T>::computeDigitalRep(T fc)
     }
 
     scaleAmplitude(aCoeff, bCoeff);
-    setCoeffs(std::move(aCoeff), std::move(bCoeff));
+    this->setCoeffs(std::move(aCoeff), std::move(bCoeff));
 }
 
 template <typename T>
@@ -164,14 +164,14 @@ void Butterworth<T>::computeBandDigitalRep(T fLower, T fUpper)
     else
         scaleAmplitude(aCoeff, bCoeff);
 
-    setCoeffs(std::move(aCoeff), std::move(bCoeff));
+    this->setCoeffs(std::move(aCoeff), std::move(bCoeff));
 }
 
 template <typename T>
 std::complex<T> Butterworth<T>::generateAnalogPole(int k, T fpw1)
 {
     auto thetaK = [pi = pi<T>, order = m_order](int k) -> T {
-        return (2 * k - 1) * pi / (2 * order);
+        return static_cast<float>(2 * k - 1) * pi / static_cast<float>(2 * order);
     };
 
     std::complex<T> analogPole(-std::sin(thetaK(k)), std::cos(thetaK(k)));
@@ -189,7 +189,7 @@ template <typename T>
 std::pair<std::complex<T>, std::complex<T>> Butterworth<T>::generateBandAnalogPole(int k, T fpw0, T bw)
 {
     auto thetaK = [pi = pi<T>, order = m_order](int k) -> T {
-        return (2 * k - 1) * pi / (2 * order);
+        return static_cast<float>(2 * k - 1) * pi / static_cast<float>(2 * order);
     };
 
     std::complex<T> analogPole(-std::sin(thetaK(k)), std::cos(thetaK(k)));

include/CMakeLists.txt

@@ -26,13 +26,17 @@
 # either expressed or implied, of the FreeBSD Project.
 
 set(HEADERS
+    BaseFilter.h
+    BaseFilter.tpp
     BilinearTransform.h
     Butterworth.h
     Butterworth.tpp
+    differentiators.h
     difi
     DigitalFilter.h
     GenericFilter.h
     GenericFilter.tpp
+    math_utils.h
     MovingAverage.h
     polynome_functions.h
     type_checks.h
@@ -43,7 +47,7 @@ set(GSL_HEADERS gsl/gsl_assert.h)
 
 add_library(${PROJECT_NAME} INTERFACE)
 target_include_directories(${PROJECT_NAME} INTERFACE $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}> $<INSTALL_INTERFACE:include>)
-target_include_directories(${PROJECT_NAME} SYSTEM INTERFACE "${EIGEN3_INCLUDE_DIR}") 
+target_include_directories(${PROJECT_NAME} SYSTEM INTERFACE "${EIGEN3_INCLUDE_DIR}")
 install(TARGETS ${PROJECT_NAME}
     EXPORT "${TARGETS_EXPORT_NAME}"
     RUNTIME DESTINATION bin

include/DigitalFilter.h

@@ -45,40 +45,12 @@ public:
     /*! \brief Constructor.
      * \param aCoeff Denominator coefficients of the filter in decreasing order.
      * \param bCoeff Numerator coefficients of the filter in decreasing order.
+     * \param type Type of the filter.
      */
-    DigitalFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff)
-        : GenericFilter<T>(aCoeff, bCoeff)
+    DigitalFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, FilterType type = FilterType::Backward)
+        : GenericFilter<T>(aCoeff, bCoeff, type)
     {
     }
-    void setCoefficients(vectX_t<T>&& aCoeff, vectX_t<T>&& bCoeff)
-    {
-        setCoeffs(std::forward(aCoeff), std::forward(bCoeff));
-    }
-};
-
-/*! \brief Basic centered digital filter.
- * 
- * This filter allows you to set any centered digital filter based on its coefficients.
- * \tparam T Floating type.
- */
-template <typename T>
-class CenteredDigitalFilter : public GenericFilter<T> {
-public:
-    /*! \brief Default uninitialized constructor. */
-    CenteredDigitalFilter() = default;
-    /*! \brief Constructor.
-     * \param aCoeff Denominator coefficients of the filter in decreasing order.
-     * \param bCoeff Numerator coefficients of the filter in decreasing order.
-     */
-    CenteredDigitalFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff)
-        : GenericFilter<T>(aCoeff, bCoeff, Type::Centered)
-    {
-    }
-    void setCoefficients(vectX_t<T>&& aCoeff, vectX_t<T>&& bCoeff)
-    {
-        setCoeffs(std::forward(aCoeff), std::forward(bCoeff));
-        setType(Type::Centered);
-    }
 };
 
 } // namespace difi

include/GenericFilter.h

@@ -27,36 +27,18 @@
 
 #pragma once
 
-#include "gsl/gsl_assert.h"
-#include "type_checks.h"
-#include "typedefs.h"
-#include <stddef.h>
-#include <string>
+#include "BaseFilter.h"
 
 namespace difi {
 
-// TODO: noexcept(Function of gsl variable)
-// TODO: constructor with universal refs
-
-/*! \brief Low-level filter.
- * 
- * It creates the basic and common functions of all linear filter that can written as a digital filter.
- * This class can not be instantiated directly.
- * 
- * \warning In Debug mode, all functions may throw if a filter is badly initialized.
- * This not the case in Realese mode.
- * 
- * \tparam T Floating type.
- */
 template <typename T>
-class GenericFilter {
-    static_assert(std::is_floating_point<T>::value && !std::is_const<T>::value, "Only accept non-complex floating point types.");
-
-public:
-    enum class Type {
-        OneSided,
-        Centered
-    };
+class GenericFilter : public BaseFilter<T, GenericFilter<T>> {
+    using Base = BaseFilter<T, GenericFilter<T>>;
+    using Base::m_isInitialized;
+    using Base::m_aCoeff;
+    using Base::m_bCoeff;
+    using Base::m_rawData;
+    using Base::m_filteredData;
 
 public:
     /*! \brief Filter a new data.
@@ -73,75 +55,59 @@ public:
      * \return Filtered signal.
      */
     vectX_t<T> filter(const vectX_t<T>& data);
-    /*! \brief Reset the data and filtered data. */
+
     void resetFilter() noexcept;
 
-    /*!< \brief Return the filter type */
-    Type type() const noexcept { return (m_center == 0 ? Type::OneSided : Type::Centered); }
-    /*! \brief Get digital filter coefficients.
+protected:
+    GenericFilter() = default;
+    GenericFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, FilterType type = FilterType::Backward)
+        : Base(aCoeff, bCoeff, type)
+    {}
+};
+
+template <typename T>
+class TVGenericFilter : public BaseFilter<T, TVGenericFilter<T>> {
+    using Base = BaseFilter<T, TVGenericFilter<T>>;
+    using Base::m_isInitialized;
+    using Base::m_aCoeff;
+    using Base::m_bCoeff;
+    using Base::m_rawData;
+    using Base::m_filteredData;
+
+public:
+    /*! \brief Filter a new data.
      * 
-     * It will automatically resize the given vectors.
-     * \param[out] aCoeff Denominator coefficients of the filter in decreasing order.
-     * \param[out] bCoeff Numerator coefficients of the filter in decreasing order.
+     * This function is practical for online application that does not know the whole signal in advance.
+     * \param data New data to filter.
+     * \return Filtered data.
      */
-    void getCoeffs(vectX_t<T>& aCoeff, vectX_t<T>& bCoeff) const noexcept;
-    /*! \brief Return coefficients of the denominator polynome. */
-    const vectX_t<T>& aCoeff() const noexcept { return m_aCoeff; }
-    /*! \brief Return coefficients of the numerator polynome. */
-    const vectX_t<T>& bCoeff() const noexcept { return m_bCoeff; }
-    /*! \brief Return the order the denominator polynome order of the filter. */
-    Eigen::Index aOrder() const noexcept { return m_aCoeff.size(); }
-    /*! \brief Return the order the numerator polynome order of the filter. */
-    Eigen::Index bOrder() const noexcept { return m_bCoeff.size(); }
-    /*! \brief Return the initialization state of the filter0 */
-    bool isInitialized() const noexcept { return m_isInitialized; }
+    T stepFilter(const T& time, const T& data);
+    /*! \brief Filter a signal.
+     * 
+     * Filter all data given by the signal.
+     * \param data Signal.
+     * \return Filtered signal.
+     */
+    vectX_t<T> filter(const vectX_t<T>& data, const vectX_t<T>& time);
+
+    void resetFilter() noexcept;
 
 protected:
-    /*! \brief Default uninitialized constructor. */
-    GenericFilter() = default;
-    /*! \brief Constructor.
-     * \param aCoeff Denominator coefficients of the filter in decreasing order.
-     * \param bCoeff Numerator coefficients of the filter in decreasing order.
-     * \param center 
-     */
-    GenericFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, Type type = Type::OneSided);
-    /*! \brief Default destructor. */
-    virtual ~GenericFilter() = default;
-
-    /*! \brief Set type of filter (one-sided or centered)
-     * 
-     * \param type The filter type.
-     * \warning bCoeff must be set before.
-     */
-    void setType(Type type);
-    /*! \brief Set the new coefficients of the filters.
-     *
-     * It awaits a universal reference.
-     * \param aCoeff Denominator coefficients of the filter in decreasing order.
-     * \param bCoeff Numerator coefficients of the filter in decreasing order.
-     */
-    template <typename T2>
-    void setCoeffs(T2&& aCoeff, T2&& bCoeff);
-    /*! \brief Normalized the filter coefficients such that aCoeff(0) = 1. */
-    void normalizeCoeffs();
-    /*! \brief Check for bad coefficients.
-     * 
-     * Set the filter status to ready is everything is fine.
-     * \param aCoeff Denominator coefficients of the filter.
-     * \param bCoeff Numerator coefficients of the filter.
-     * \return True if the filter status is set on READY.
-     */
-    bool checkCoeffs(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, Type type);
+    TVGenericFilter() = default;
+    TVGenericFilter(size_t differentialOrder, const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, FilterType type = FilterType::Backward)
+        : Base()
+        , m_diffOrder(differentialOrder)
+    {
+        Expects(differentialOrder >= 1);
+        this->setCoeffs(aCoeff, bCoeff);
+        this->setType(type);
+    }
 
 private:
-    Eigen::Index m_center = 0; /*!< Center of the filter. 0 is a one-sided filter. Default is 0. */
-    bool m_isInitialized = false; /*!< Initialization state of the filter. Default is false */
-    vectX_t<T> m_aCoeff; /*!< Denominator coefficients of the filter */
-    vectX_t<T> m_bCoeff; /*!< Numerator coefficients of the filter */
-    vectX_t<T> m_filteredData; /*!< Last set of filtered data */
-    vectX_t<T> m_rawData; /*!< Last set of non-filtered data */
+    size_t m_diffOrder = 1;
+    vectX_t<T> m_timers;
 };
 
 } // namespace difi
 
-#include "GenericFilter.tpp"
+#include "GenericFilter.tpp"

include/GenericFilter.tpp

@@ -25,12 +25,8 @@
 // of the authors and should not be interpreted as representing official policies,
 // either expressed or implied, of the FreeBSD Project.
 
-#include <limits>
-
 namespace difi {
 
-// Public functions
-
 template <typename T>
 T GenericFilter<T>::stepFilter(const T& data)
 {
@@ -42,10 +38,10 @@ T GenericFilter<T>::stepFilter(const T& data)
     for (Eigen::Index i = m_filteredData.size() - 1; i > 0; --i)
         m_filteredData(i) = m_filteredData(i - 1);
 
-    m_rawData[0] = data;
-    m_filteredData[0] = 0;
-    m_filteredData[m_center] = m_bCoeff.dot(m_rawData) - m_aCoeff.dot(m_filteredData);
-    return m_filteredData[m_center];
+    m_rawData(0) = data;
+    m_filteredData(0) = 0;
+    m_filteredData(0) = m_bCoeff.dot(m_rawData) - m_aCoeff.dot(m_filteredData);
+    return m_filteredData(0);
 }
 
 template <typename T>
@@ -66,65 +62,50 @@ void GenericFilter<T>::resetFilter() noexcept
 }
 
 template <typename T>
-void GenericFilter<T>::setType(Type type)
+T TVGenericFilter<T>::stepFilter(const T& time, const T& data)
 {
-    Expects(type == Type::Centered ? m_bCoeff.size() > 2 && m_bCoeff.size() % 2 == 1 : true);
-    m_center = (type == Type::OneSided ? 0 : (m_bCoeff.size() - 1) / 2);
+    Expects(m_isInitialized);
+
+    // Slide data (can't use SIMD, but should be small)
+    for (Eigen::Index i = m_rawData.size() - 1; i > 0; --i) {
+        m_rawData(i) = m_rawData(i - 1);
+        m_timers(i) = m_timers(i - 1);
+    }
+    for (Eigen::Index i = m_filteredData.size() - 1; i > 0; --i)
+        m_filteredData(i) = m_filteredData(i - 1);
+
+    m_timers(0) = time;
+    const Eigen::Index M = (m_rawData.size() - 1) / 2;
+    auto bCoeff = m_bCoeff;
+    for (Eigen::Index i = 1; i < M + 1; ++i) {
+        const T diff = std::pow(m_timers(M - i) - m_timers(M + i), m_diffOrder);
+        bCoeff(M + i) /= diff;
+        bCoeff(M - i) /= diff;
+        bCoeff(M) -= (bCoeff(M - i) + bCoeff(M + i));
+    }
+    m_rawData(0) = data;
+    m_filteredData(0) = 0;
+    m_filteredData(0) = bCoeff.dot(m_rawData) - m_aCoeff.dot(m_filteredData);
+    return m_filteredData(0);
 }
 
 template <typename T>
-template <typename T2>
-void GenericFilter<T>::setCoeffs(T2&& aCoeff, T2&& bCoeff)
+vectX_t<T> TVGenericFilter<T>::filter(const vectX_t<T>& data, const vectX_t<T>& time)
 {
-    static_assert(std::is_convertible_v<T2, vectX_t<T>>, "The coefficients types should be convertible to vectX_t<T>");
-
-    Expects(checkCoeffs(aCoeff, bCoeff, (m_center == 0 ? Type::OneSided : Type::Centered)));
-    m_aCoeff = aCoeff;
-    m_bCoeff = bCoeff;
-    normalizeCoeffs();
-    resetFilter();
-    m_isInitialized = true;
+    Expects(m_isInitialized);
+    Expects(data.size() == time.size());
+    vectX_t<T> results(data.size());
+    for (Eigen::Index i = 0; i < data.size(); ++i)
+        results(i) = stepFilter(time(i), data(i));
+    return results;
 }
 
 template <typename T>
-void GenericFilter<T>::getCoeffs(vectX_t<T>& aCoeff, vectX_t<T>& bCoeff) const noexcept
+void TVGenericFilter<T>::resetFilter() noexcept
 {
-    aCoeff = m_aCoeff;
-    bCoeff = m_bCoeff;
-}
-
-// Protected functions
-
-template <typename T>
-GenericFilter<T>::GenericFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, Type type)
-    : m_aCoeff(aCoeff)
-    , m_bCoeff(bCoeff)
-    , m_filteredData(aCoeff.size())
-    , m_rawData(bCoeff.size())
-{
-    Expects(checkCoeffs(aCoeff, bCoeff, type));
-    m_center = (type == Type::OneSided ? 0 : (bCoeff.size() - 1) / 2);
-    normalizeCoeffs();
-    resetFilter();
-    m_isInitialized = true;
-}
-
-template <typename T>
-void GenericFilter<T>::normalizeCoeffs()
-{
-    T a0 = m_aCoeff(0);
-    if (std::abs(a0 - T(1)) < std::numeric_limits<T>::epsilon())
-        return;
-
-    m_aCoeff /= a0;
-    m_bCoeff /= a0;
-}
-
-template <typename T>
-bool GenericFilter<T>::checkCoeffs(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, Type type)
-{
-    bool centering = (type == Type::Centered ? (bCoeff.size() % 2 == 1) : true);
-    return aCoeff.size() > 0 && std::abs(aCoeff[0]) > std::numeric_limits<T>::epsilon() && bCoeff.size() > 0 && centering;
+    m_filteredData.setZero(m_aCoeff.size());
+    m_rawData.setZero(m_bCoeff.size());
+    m_timers.setZero(m_bCoeff.size());
 }
 
 } // namespace difi

include/MovingAverage.h

@@ -45,47 +45,21 @@ public:
     MovingAverage() = default;
     /*! \brief Constructor.
      * \param windowSize Size of the moving average window.
+     * \param type Type of the filter.
      */
-    MovingAverage(int windowSize)
+    MovingAverage(int windowSize, FilterType type = FilterType::Backward)
     {
         setWindowSize(windowSize);
+        this->setType(type);
     }
     /*! \brief Set the size of the moving average window. */
     void setWindowSize(int windowSize)
     {
         Expects(windowSize > 0);
-        setCoeffs(vectX_t<T>::Constant(1, T(1)), vectX_t<T>::Constant(windowSize, T(1) / windowSize));
+        this->setCoeffs(vectX_t<T>::Constant(1, T(1)), vectX_t<T>::Constant(windowSize, T(1) / static_cast<T>(windowSize)));
     }
     /*! \brief Get the size of the moving average window. */
-    int windowSize() const noexcept { return bOrder(); }
-};
-
-/*! \brief Centered moving average digital filter.
- * 
- * This is a specialization of a digital filter in order to use a centered moving average.
- * \tparam T Floating type.
- */
-template <typename T>
-class CenteredMovingAverage : public DigitalFilter<T> {
-public:
-    /*! \brief Default uninitialized constructor. */
-    CenteredMovingAverage() = default;
-    /*! \brief Constructor.
-     * \param windowSize Size of the moving average window.
-     */
-    CenteredMovingAverage(int windowSize)
-    {
-        setWindowSize(windowSize);
-    }
-    /*! \brief Set the size of the moving average window. */
-    void setWindowSize(int windowSize)
-    {
-        Expects(windowSize > 2 && windowSize % 2 == 1);
-        setCoeffs(vectX_t<T>::Constant(1, T(1)), vectX_t<T>::Constant(windowSize, T(1) / windowSize));
-        setType(Type::Centered);
-    }
-    /*! \brief Get the size of the moving average window. */
-    int windowSize() const noexcept { return bOrder(); }
+    int windowSize() const noexcept { return this->bOrder(); }
 };
 
 } // namespace difi

include/differentiators.h

@@ -0,0 +1,368 @@
+// Copyright (c) 2019, Vincent SAMY
+// All rights reserved.
+
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+
+// 1. Redistributions of source code must retain the above copyright notice,
+//    this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+// The views and conclusions contained in the software and documentation are those
+// of the authors and should not be interpreted as representing official policies,
+// either expressed or implied, of the FreeBSD Project.
+
+#pragma once
+#include "DigitalFilter.h"
+#include "math_utils.h"
+
+// Read this if you are an adulator of the math god: https://arxiv.org/pdf/1709.08321.pdf
+
+namespace difi {
+
+namespace details {
+
+// T: type
+// N: Number of points
+
+// Centered differentiators: http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/central-differences/
+template <typename T, int N> struct GetCDCoeffs {
+    vectN_t<T, N> operator()() const;
+};
+template <typename T> struct GetCDCoeffs<T, 3> {
+    vectN_t<T, 3> operator()() const { return (vectN_t<T, 3>() << T(1), T(0), T(-1)).finished() / T(2); }
+};
+template <typename T> struct GetCDCoeffs<T, 5> {
+    vectN_t<T, 5> operator()() const { return (vectN_t<T, 5>() << T(-1), T(8), T(0), T(8), T(1)).finished() / T(12); }
+};
+template <typename T> struct GetCDCoeffs<T, 7> {
+    vectN_t<T, 7> operator()() const { return (vectN_t<T, 7>() << T(1), T(-9), T(45), T(0), T(-45), T(9), T(-1)).finished() / T(60); }
+};
+template <typename T> struct GetCDCoeffs<T, 9> {
+    vectN_t<T, 9> operator()() const { return (vectN_t<T, 9>() << T(-3), T(32), T(-168), T(672), T(0), T(-672), T(168), T(-32), T(3)).finished() / T(840); }
+};
+
+// Low-noise Lanczos differentiators: http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/lanczos-low-noise-differentiators/
+template <typename T, int N>
+struct GetLNLCoeffs {
+    vectN_t<T, N> operator()() const
+    {
+        static_assert(N > 2 && N % 2 == 1, "'N' must be odd.");
+        constexpr const int M = (N - 1) / 2;
+        constexpr const int Den = M * (M + 1) * (2 * M + 1);
+
+        vectN_t<T, N> v{};
+        v(M) = T(0);
+        for (int k = 0; k < M; ++k) {
+            v(k) = T(3) * static_cast<T>(M - k) / static_cast<T>(Den);
+            v(N - k - 1) = -v(k);
+        }
+        return v;
+    }
+};
+
+// Super Low-noise Lanczos differentiators: http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/lanczos-low-noise-differentiators/
+template <typename T, int N> struct GetSLNLCoeffs {
+    vectN_t<T, N> operator()() const;
+};
+template <typename T> struct GetSLNLCoeffs<T, 7> {
+    vectN_t<T, 7> operator()() const { return (vectN_t<T, 7>() << T(-22), T(67), T(58), T(0), T(-58), T(-67), T(22)).finished() / T(252); }
+};
+template <typename T> struct GetSLNLCoeffs<T, 9> {
+    vectN_t<T, 9> operator()() const { return (vectN_t<T, 9>() << T(-86), T(142), T(193), T(126), T(0), T(-126), T(-193), T(-142), T(86)).finished() / T(1188); }
+};
+template <typename T> struct GetSLNLCoeffs<T, 11> {
+    vectN_t<T, 11> operator()() const { return (vectN_t<T, 11>() << T(-300), T(294), T(532), T(503), T(296), T(0), T(-296), T(-503), T(-532), T(-294), T(300)).finished() / T(5148); }
+};
+
+// Backward Noise-Robust differentiators; http://www.holoborodko.com/pavel/wp-content/uploads/OneSidedNoiseRobustDifferentiators.pdf
+template <typename T, int N>
+struct GetFNRCoeffs {
+    vectN_t<T, N> operator()() const
+    {
+        static_assert(N >= 2, "N should be greater than 2");
+        constexpr const int BinCoeff = N - 2;
+        constexpr const int M = N / 2;
+
+        vectN_t<T, N> v{};
+        v(0) = T(1);
+        v(M) = T(0);
+        v(N - 1) = T(-1);
+        for (int i = 1; i < M; ++i) {
+            v(i) = Binomial<T>(BinCoeff, i) - Binomial<T>(BinCoeff, i - 1);
+            v(N - i - 1) = -v(i);
+        }
+
+        v /= std::pow(T(2), T(BinCoeff));
+        return v;
+    }
+};
+
+// Backward Hybrid Noise-Robust differentiators; http://www.holoborodko.com/pavel/wp-content/uploads/OneSidedNoiseRobustDifferentiators.pdf
+template <typename T, int N> struct GetFHNRCoeffs {
+    vectN_t<T, N> operator()() const;
+};
+template <typename T> struct GetFHNRCoeffs<T, 4> {
+    vectN_t<T, 4> operator()() const { return (vectN_t<T, 4>() << T(2), T(-1), T(-2), T(1)).finished() / T(2); }
+};
+template <typename T> struct GetFHNRCoeffs<T, 5> {
+    vectN_t<T, 5> operator()() const { return (vectN_t<T, 5>() << T(7), T(1), T(-10), T(-1), T(3)).finished() / T(10); }
+};
+template <typename T> struct GetFHNRCoeffs<T, 6> {
+    vectN_t<T, 6> operator()() const { return (vectN_t<T, 6>() << T(16), T(1), T(-10), T(-10), T(-6), T(9)).finished() / T(28); }
+};
+template <typename T> struct GetFHNRCoeffs<T, 7> {
+    vectN_t<T, 7> operator()() const { return (vectN_t<T, 7>() << T(12), T(5), T(-8), T(-6), T(-10), T(1), T(6)).finished() / T(28); }
+};
+template <typename T> struct GetFHNRCoeffs<T, 8> {
+    vectN_t<T, 8> operator()() const { return (vectN_t<T, 8>() << T(22), T(7), T(-6), T(-11), T(-14), T(-9), T(-2), T(13)).finished() / T(60); }
+};
+template <typename T> struct GetFHNRCoeffs<T, 9> {
+    vectN_t<T, 9> operator()() const { return (vectN_t<T, 9>() << T(52), T(29), T(-14), T(-17), T(-40), T(-23), T(-26), T(11), T(28)).finished() / T(180); }
+};
+template <typename T> struct GetFHNRCoeffs<T, 10> {
+    vectN_t<T, 10> operator()() const { return (vectN_t<T, 10>() << T(56), T(26), T(-2), T(-17), T(-30), T(-30), T(-28), T(-13), T(4), T(34)).finished() / T(220); }
+};
+template <typename T> struct GetFHNRCoeffs<T, 11> {
+    vectN_t<T, 11> operator()() const { return (vectN_t<T, 11>() << T(320), T(206), T(-8), T(-47), T(-186), T(-150), T(-214), T(-103), T(-92), T(94), T(180)).finished() / T(1540); }
+};
+template <typename T> struct GetFHNRCoeffs<T, 16> {
+    vectN_t<T, 16> operator()() const { return (vectN_t<T, 16>() << T(322), T(217), T(110), T(35), T(-42), T(-87), T(-134), T(-149), T(-166), T(-151), T(-138), T(-93), T(-50), T(28), T(98), T(203)).finished() / T(2856); }
+};
+
+template <typename T, int N, typename BackwardCoeffs> vectN_t<T, N> GetBackwardISDCoeffs()
+{
+    vectN_t<T, N> v{};
+    const vectN_t<T, N> v0 = BackwardCoeffs{}();
+    for (Eigen::Index k = 0; k < N; ++k)
+        v(k) = k * v0(k);
+    return v;
+}
+// Backward Noise-Robust differentiators for irregular space data
+template <typename T, int N> struct GetFNRISDCoeffs {
+    vectN_t<T, N> operator()() const { return GetBackwardISDCoeffs<T, N, GetFNRCoeffs<T, N>>(); }
+};
+// Backward Hybrid Noise-Robust differentiators for irregular space data
+template <typename T, int N> struct GetFHNRISDCoeffs {
+    vectN_t<T, N> operator()() const { return GetBackwardISDCoeffs<T, N, GetFHNRCoeffs<T, N>>(); }
+};
+
+// Centered Noise-Robust differentiators (tangency at 2nd order): http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
+template <typename T, int N>
+struct GetCNR2Coeffs {
+    vectN_t<T, N> operator()() const
+    {
+        static_assert(N % 2 == 1., "'N' must be odd.");
+        return GetFNRCoeffs<T, N>{}(); // Same coefficients
+    }
+};
+
+// Centered Noise-Robust  differentiators (tangency at 4th order): http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
+template <typename T, int N> struct GetCNR4Coeffs {
+    vectN_t<T, N> operator()() const;
+};
+template <typename T> struct GetCNR4Coeffs<T, 7> {
+    vectN_t<T, 7> operator()() const { return (vectN_t<T, 7>() << T(-5), T(12), T(39), T(0), T(-39), T(-12), T(5)).finished() / T(96); }
+};
+template <typename T> struct GetCNR4Coeffs<T, 9> {
+    vectN_t<T, 9> operator()() const { return (vectN_t<T, 9>() << T(-2), T(-1), T(16), T(27), T(0), T(-27), T(-16), T(1), T(2)).finished() / T(96); }
+};
+template <typename T> struct GetCNR4Coeffs<T, 11> {
+    vectN_t<T, 11> operator()() const { return (vectN_t<T, 11>() << T(-11), T(-32), T(39), T(256), T(322), T(0), T(-322), T(-256), T(-39), T(32), T(11)).finished() / T(1536); }
+};
+
+// Centered Noise-Robust differentiators for irregular space data: http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
+template <typename T, int N, typename CNRCoeffs> vectN_t<T, N> GetCNRISDCoeffs()
+{
+    constexpr const int M = (N - 1) / 2;
+    vectN_t<T, N> v{};
+    const vectN_t<T, N> v0 = CNRCoeffs{}();
+    v(M) = 0;
+    for (int k = 1; k < M + 1; ++k) {
+        v(M - k) = T(2) * k * v0(M - k);
+        v(M + k) = T(2) * k * v0(M + k);
+    }
+
+    return v;
+}
+template <typename T, int N> struct GetCNR2ISDCoeffs {
+    vectN_t<T, N> operator()() const { return GetCNRISDCoeffs<T, N, GetCNR2Coeffs<T, N>>(); }
+};
+template <typename T, int N> struct GetCNR4ISDCoeffs {
+    vectN_t<T, N> operator()() const { return GetCNRISDCoeffs<T, N, GetCNR4Coeffs<T, N>>(); }
+};
+
+/*
+ * Second order differentiators
+ */
+
+template <typename T>
+constexpr T GetSONRBaseCoeff(int N, int M, int k)
+{
+    if (k > M)
+        return T(0);
+    else if (k == M)
+        return T(1);
+    else
+        return ((T(2) * N - T(10)) * GetSONRBaseCoeff<T>(N, M, k + 1) - (N + T(2) * k + T(3)) * GetSONRBaseCoeff<T>(N, M, k + 2)) / (N - T(2) * k - T(1));
+}
+
+template <typename T, int N> vectN_t<T, (N - 1) / 2 + 1> GetSONRBaseCoeffs()
+{
+    static_assert(N >= 5 && N % 2 == 1, "N must be a odd number >= 5");
+    constexpr const int M = (N - 1) / 2;
+    vectN_t<T, M + 1> s{};
+    for (int k = 0; k < M + 1; ++k)
+        s(k) = GetSONRBaseCoeff<T>(N, M, k);
+
+    return s;
+}
+
+// Second-Order Centered Noise-Robust differentiator: http://www.holoborodko.com/pavel/downloads/NoiseRobustSecondDerivative.pdf
+template <typename T, int N>
+struct GetSOCNRCoeffs {
+    vectN_t<T, N> operator()() const
+    {
+        constexpr const int M = (N - 1) / 2;
+        constexpr const T Den = pow(2, N - 3);
+        vectN_t<T, N> v{};
+        vectN_t<T, M + 1> s = GetSONRBaseCoeffs<T, N>();
+        v(M) = s(0);
+        for (int k = 1; k < M + 1; ++k) {
+            v(M + k) = s(k);
+            v(M - k) = s(k);
+        }
+
+        v /= Den;
+        return v;
+    }
+};
+// Second-Order Backward Noise-Robust differentiator: http://www.holoborodko.com/pavel/downloads/NoiseRobustSecondDerivative.pdf
+template <typename T, int N> struct GetSOFNRCoeffs {
+    vectN_t<T, N> operator()() const { return GetSOCNRCoeffs<T, N>(); }
+}; // Coefficients are the same.
+
+// Second-Order Centered Noise-Robust Irregular Space Data differentiator: http://www.holoborodko.com/pavel/downloads/NoiseRobustSecondDerivative.pdf
+template <typename T, int N>
+struct GetSOCNRISDCoeffs {
+    vectN_t<T, N> operator()() const
+    {
+        constexpr const int M = (N - 1) / 2;
+        constexpr const T Den = pow(2, N - 3);
+
+        vectN_t<T, N> v{};
+        const vectN_t<T, M + 1> s = GetSONRBaseCoeffs<T, N>();
+
+        const auto alpha = [&s](int k) -> T { return T(4) * k * k * s(k); };
+
+        v(M) = T(0);
+        for (int k = 1; k < M + 1; ++k) {
+            auto alph = alpha(k);
+            v(M) -= T(2) * alph;
+            v(M + k) = alph;
+            v(M - k) = alph;
+        }
+
+        v /= Den;
+        return v;
+    }
+};
+
+// Second-Order Backward Noise-Robust Irregular Space Data differentiator: http://www.holoborodko.com/pavel/downloads/NoiseRobustSecondDerivative.pdf
+template <typename T, int N> struct GetSOFNRISDCoeffs {
+    vectN_t<T, N> operator()() const { return GetSOCNRISDCoeffs<T, N>(); }
+}; // Same coefficients
+
+/*
+ * Differentiator Generator
+ */
+
+template <typename T, int N, int Order, typename CoeffGetter>
+class BackwardDifferentiator : public GenericFilter<T> {
+public:
+    BackwardDifferentiator()
+        : GenericFilter<T>(vectX_t<T>::Constant(1, T(1)), CoeffGetter{}())
+    {}
+    BackwardDifferentiator(T timestep)
+        : GenericFilter<T>(vectX_t<T>::Constant(1, T(1)), CoeffGetter{}() / std::pow(timestep, Order))
+    {}
+    void setTimestep(T timestep) { this->setCoeffs(vectX_t<T>::Constant(1, T(1)), CoeffGetter{}() / std::pow(timestep, Order)); }
+    T timestep() const noexcept { return std::pow(this->bCoeff()(0) / CoeffGetter{}()(0), T(1) / Order); }
+};
+
+template <typename T, int N, int Order, typename CoeffGetter>
+class CenteredDifferentiator : public GenericFilter<T> {
+public:
+    CenteredDifferentiator()
+        : GenericFilter<T>(vectX_t<T>::Constant(1, T(1)), CoeffGetter{}(), FilterType::Centered)
+    {}
+    CenteredDifferentiator(T timestep)
+        : GenericFilter<T>(vectX_t<T>::Constant(1, T(1)), CoeffGetter{}() / std::pow(timestep, Order))
+    {}
+    void setTimestep(T timestep) { this->setCoeffs(vectX_t<T>::Constant(1, T(1)), CoeffGetter{}() / std::pow(timestep, Order)); }
+    T timestep() const noexcept { return std::pow(this->bCoeff()(0) / CoeffGetter{}()(0), T(1) / Order); }
+};
+
+template <typename T, int N, int Order, typename CoeffGetter>
+class TVBackwardDifferentiator : public TVGenericFilter<T> {
+    static_assert(Order >= 1, "Order must be greater or equal to 1");
+
+public:
+    TVBackwardDifferentiator()
+        : TVGenericFilter<T>(Order, vectX_t<T>::Constant(1, T(1)), CoeffGetter{}())
+    {}
+};
+
+template <typename T, int N, int Order, typename CoeffGetter>
+class TVCenteredDifferentiator : public TVGenericFilter<T> {
+    static_assert(Order >= 1, "Order must be greater or equal to 1");
+
+public:
+    TVCenteredDifferentiator()
+        : TVGenericFilter<T>(Order, vectX_t<T>::Constant(1, T(1)), CoeffGetter{}(), FilterType::Centered)
+    {}
+};
+
+} // namespace details
+
+// Backward differentiators
+template <typename T, int N> using BackwardDiffNoiseRobust = details::BackwardDifferentiator<T, N, 1, details::GetFNRCoeffs<T, N>>;
+template <typename T, int N> using BackwardDiffHybridNoiseRobust = details::BackwardDifferentiator<T, N, 1, details::GetFHNRCoeffs<T, N>>;
+// Time-Varying backward differentiators
+template <typename T, int N> using TVBackwardDiffNoiseRobust = details::TVBackwardDifferentiator<T, N, 1, details::GetFNRISDCoeffs<T, N>>;
+template <typename T, int N> using TVBackwardDiffHybridNoiseRobust = details::TVBackwardDifferentiator<T, N, 1, details::GetFHNRISDCoeffs<T, N>>;
+
+// Centered differentiators
+template <typename T, int N> using CenteredDiffBasic = details::CenteredDifferentiator<T, N, 1, details::GetCDCoeffs<T, N>>;
+template <typename T, int N> using CenteredDiffLowNoiseLanczos = details::CenteredDifferentiator<T, N, 1, details::GetLNLCoeffs<T, N>>;
+template <typename T, int N> using CenteredDiffSuperLowNoiseLanczos = details::CenteredDifferentiator<T, N, 1, details::GetSLNLCoeffs<T, N>>;
+template <typename T, int N> using CenteredDiffNoiseRobust2 = details::CenteredDifferentiator<T, N, 1, details::GetCNR2Coeffs<T, N>>;
+template <typename T, int N> using CenteredDiffNoiseRobust4 = details::CenteredDifferentiator<T, N, 1, details::GetCNR4Coeffs<T, N>>;
+// Time-Varying centered differentiators
+template <typename T, int N> using TVCenteredDiffNoiseRobust2 = details::TVCenteredDifferentiator<T, N, 1, details::GetCNR2ISDCoeffs<T, N>>;
+template <typename T, int N> using TVCenteredDiffNoiseRobust4 = details::TVCenteredDifferentiator<T, N, 1, details::GetCNR4ISDCoeffs<T, N>>;
+
+// Second-order backward differentiators
+template <typename T, int N> using BackwardDiffSecondOrder = details::BackwardDifferentiator<T, N, 2, details::GetSOFNRCoeffs<T, N>>;
+// Second-order Time-Varying backward differentiators
+template <typename T, int N> using TVBackwardDiffSecondOrder = details::TVBackwardDifferentiator<T, N, 2, details::GetSOFNRISDCoeffs<T, N>>;
+
+// Second-order centered differentiators
+template <typename T, int N> using CenteredDiffSecondOrder = details::CenteredDifferentiator<T, N, 2, details::GetSOCNRCoeffs<T, N>>;
+// Second-order Time-Varying centered differentiators
+template <typename T, int N> using TVCenteredDiffSecondOrder = details::TVCenteredDifferentiator<T, N, 2, details::GetSOCNRISDCoeffs<T, N>>;
+
+} // namespace difi

include/difi

@@ -32,6 +32,7 @@
 #include "DigitalFilter.h"
 #include "GenericFilter.h"
 #include "MovingAverage.h"
+#include "differentiators.h"
 #include "polynome_functions.h"
 #include "typedefs.h"
 
@@ -59,4 +60,42 @@ using BilinearTransformd = BilinearTransform<double>;
 using BilinearTransformcf = BilinearTransform<std::complex<float>>;
 using BilinearTransformcd = BilinearTransform<std::complex<double>>;
 
+// 1st order centered differentiators
+template <int N> using CenteredDiffBasicf = CenteredDiffBasic<float, N>;
+template <int N> using CenteredDiffBasicd = CenteredDiffBasic<double, N>;
+template <int N> using CenteredDiffLowNoiseLanczosf = CenteredDiffLowNoiseLanczos<float, N>;
+template <int N> using CenteredDiffLowNoiseLanczosd = CenteredDiffLowNoiseLanczos<double, N>;
+template <int N> using CenteredDiffSuperLowNoiseLanczosf = CenteredDiffSuperLowNoiseLanczos<float, N>;
+template <int N> using CenteredDiffSuperLowNoiseLanczosd = CenteredDiffSuperLowNoiseLanczos<double, N>;
+template <int N> using CenteredDiffNoiseRobust2f = CenteredDiffNoiseRobust2<float, N>;
+template <int N> using CenteredDiffNoiseRobust2d = CenteredDiffNoiseRobust2<double, N>;
+template <int N> using CenteredDiffNoiseRobust4f = CenteredDiffNoiseRobust4<float, N>;
+template <int N> using CenteredDiffNoiseRobust4d = CenteredDiffNoiseRobust4<double, N>;
+template <int N> using TVCenteredDiffNoiseRobust2f = TVCenteredDiffNoiseRobust2<float, N>;
+template <int N> using TVCenteredDiffNoiseRobust2d = TVCenteredDiffNoiseRobust2<double, N>;
+template <int N> using TVCenteredDiffNoiseRobust4f = TVCenteredDiffNoiseRobust4<float, N>;
+template <int N> using TVCenteredDiffNoiseRobust4d = TVCenteredDiffNoiseRobust4<double, N>;
+
+// 2nd order centered differentiators
+template <int N> using CenteredDiffSecondOrderf = CenteredDiffSecondOrder<float, N>;
+template <int N> using CenteredDiffSecondOrderd = CenteredDiffSecondOrder<double, N>;
+template <int N> using TVCenteredDiffSecondOrderf = TVCenteredDiffSecondOrder<float, N>;
+template <int N> using TVCenteredDiffSecondOrderd = TVCenteredDiffSecondOrder<double, N>;
+
+// 1st order backward differentiators
+template <int N> using BackwardDiffNoiseRobustf = BackwardDiffNoiseRobust<float, N>;
+template <int N> using BackwardDiffNoiseRobustd = BackwardDiffNoiseRobust<double, N>;
+template <int N> using BackwardDiffHybridNoiseRobustf = BackwardDiffHybridNoiseRobust<float, N>;
+template <int N> using BackwardDiffHybridNoiseRobustd = BackwardDiffHybridNoiseRobust<double, N>;
+template <int N> using TVBackwardDiffNoiseRobustf = TVBackwardDiffNoiseRobust<float, N>;
+template <int N> using TVBackwardDiffNoiseRobustd = TVBackwardDiffNoiseRobust<double, N>;
+template <int N> using TVBackwardDiffHybridNoiseRobustf = TVBackwardDiffHybridNoiseRobust<float, N>;
+template <int N> using TVBackwardDiffHybridNoiseRobustd = TVBackwardDiffHybridNoiseRobust<double, N>;
+
+// 2nd order backward differentiators
+template <int N> using BackwardDiffSecondOrderf = BackwardDiffSecondOrder<float, N>;
+template <int N> using BackwardDiffSecondOrderd = BackwardDiffSecondOrder<double, N>;
+template <int N> using TVBackwardDiffSecondOrderf = TVBackwardDiffSecondOrder<float, N>;
+template <int N> using TVBackwardDiffSecondOrderd = TVBackwardDiffSecondOrder<double, N>;
+
 } // namespace difi

include/math_utils.h

@@ -0,0 +1,61 @@
+// Copyright (c) 2019, Vincent SAMY
+// All rights reserved.
+
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+
+// 1. Redistributions of source code must retain the above copyright notice,
+//    this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+// The views and conclusions contained in the software and documentation are those
+// of the authors and should not be interpreted as representing official policies,
+// either expressed or implied, of the FreeBSD Project.
+
+#pragma once
+#include <type_traits>
+
+namespace difi {
+
+// https://stackoverflow.com/questions/44718971/calculate-binomial-coffeficient-very-reliably
+template <typename T>
+constexpr T Binomial(int n, int k)
+{
+    if (k > n)
+        return T(0);
+    else if (k == 0 || k == n)
+        return T(1);
+    else if (k == 1 || k == n - 1)
+        return T(n);
+    else if (2 * k < n)
+        return Binomial<T>(n - 1, k - 1) * n / k;
+    else
+        return Binomial<T>(n - 1, k) * n / (n - k);
+}
+
+template <typename T>
+constexpr T pow(T n, T k)
+{
+    static_assert(std::is_integral<T>::value, "For integral values only, eitherwise, use std::pow");
+    if (n == 0)
+        return (k == 0) ? 1 : 0;
+    else if (k == 0)
+        return T(1);
+    else
+        return n * pow(n, k - 1);
+}
+
+} // namespace difi

include/typedefs.h

@@ -37,7 +37,15 @@ constexpr T pi = T(3.14159265358979323846264338327950288419716939937510582097494
 template <typename T>
 using vectX_t = Eigen::Matrix<T, Eigen::Dynamic, 1>; /*!< Eigen column-vector */
 
+template <typename T, int N>
+using vectN_t = Eigen::Matrix<T, N, 1>; /*!< Fixed Eigen column-vector */
+
 template <typename T>
 using vectXc_t = vectX_t<std::complex<T>>; /*!< Eigen complex column-vector */
 
+enum class FilterType {
+    Backward,
+    Centered
+};
+
 } // namespace difi

tests/ButterworthFilterTests.cpp

@@ -25,14 +25,13 @@
 // of the authors and should not be interpreted as representing official policies,
 // either expressed or implied, of the FreeBSD Project.
 
-#define BOOST_TEST_MODULE ButterworthFilterTests
-
 // Note: In term of precision, LP > HP > BP ~= BR
 
 #include "difi"
+#include "doctest/doctest.h"
+#include "doctest_helper.h"
 #include "test_functions.h"
 #include "warning_macro.h"
-#include <boost/test/unit_test.hpp>
 
 DISABLE_CONVERSION_WARNING_BEGIN
 
@@ -64,92 +63,51 @@ struct System {
 
 DISABLE_CONVERSION_WARNING_END
 
-BOOST_AUTO_TEST_CASE(FIND_BUTTERWORTH_LP_HP_FLOAT)
+TEST_CASE_TEMPLATE("Find butterworth Low pass and High pass", T, float, double)
 {
     // LP
-    auto butterRequirement = difi::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);
+    auto butterRequirement = difi::Butterworth<T>::findMinimumButter(T(40. / 500.), T(150. / 500.), T(3.), T(60.));
+    REQUIRE_EQUAL(5, butterRequirement.first);
+    REQUIRE_SMALL(std::abs(T(0.081038494957764) - butterRequirement.second), std::numeric_limits<T>::epsilon() * 10);
 
     // HP
-    butterRequirement = difi::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);
+    butterRequirement = difi::Butterworth<T>::findMinimumButter(T(150. / 500.), T(40. / 500.), T(3.), T(60.));
+    REQUIRE_EQUAL(5, butterRequirement.first);
+    REQUIRE_SMALL(std::abs(T(0.296655824107340) - butterRequirement.second), std::numeric_limits<T>::epsilon() * 10);
 }
 
-BOOST_AUTO_TEST_CASE(FIND_BUTTERWORTH_LP_HP_DOUBLE)
+TEST_CASE_TEMPLATE("Butterworth low pass filter", T, float, double)
 {
-    // LP
-    auto butterRequirement = difi::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 = difi::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);
+    System<T> s;
+    auto bf = difi::Butterworth<T>(s.order, s.fc, s.fs);
+    REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
+    test_coeffs(s.lpACoeffRes, s.lpBCoeffRes, bf, std::numeric_limits<T>::epsilon() * 10);
+    test_results(s.lpResults, s.data, bf, std::numeric_limits<T>::epsilon() * 100);
 }
 
-BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_LP_FILTER_FLOAT, System<float>)
+TEST_CASE_TEMPLATE("Butterworth high pass filter", T, float, double)
 {
-    auto bf = difi::Butterworthf(order, fc, fs);
-    BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
-    test_coeffs(lpACoeffRes, lpBCoeffRes, bf, std::numeric_limits<float>::epsilon() * 10);
-    test_results(lpResults, data, bf, std::numeric_limits<float>::epsilon() * 100);
+    System<T> s;
+    auto bf = difi::Butterworth<T>(s.order, s.fc, s.fs, difi::Butterworth<T>::Type::HighPass);
+    REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
+    test_coeffs(s.hpACoeffRes, s.hpBCoeffRes, bf, std::numeric_limits<T>::epsilon() * 10);
+    test_results(s.hpResults, s.data, bf, std::numeric_limits<T>::epsilon() * 1000);
 }
 
-BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_LP_FILTER_DOUBLE, System<double>)
+TEST_CASE_TEMPLATE("Butterworth band pass filter", T, float, double)
 {
-    auto bf = difi::Butterworthd(order, fc, fs);
-    BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
-    test_coeffs(lpACoeffRes, lpBCoeffRes, bf, std::numeric_limits<double>::epsilon() * 10);
-    test_results(lpResults, data, bf, std::numeric_limits<double>::epsilon() * 100);
+    System<T> s;
+    auto bf = difi::Butterworth<T>(s.order, s.fLower, s.fUpper, s.fs);
+    REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
+    test_coeffs(s.bpACoeffRes, s.bpBCoeffRes, bf, std::numeric_limits<T>::epsilon() * 1000);
+    test_results(s.bpResults, s.data, bf, std::numeric_limits<T>::epsilon() * 10000);
 }
 
-BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_HP_FILTER_FLOAT, System<float>)
+TEST_CASE_TEMPLATE("Butterworth band-reject filter", T, float, double)
 {
-    auto bf = difi::Butterworthf(order, fc, fs, difi::Butterworthf::Type::HighPass);
-    BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
-    test_coeffs(hpACoeffRes, hpBCoeffRes, bf, std::numeric_limits<float>::epsilon() * 10);
-    test_results(hpResults, data, bf, std::numeric_limits<float>::epsilon() * 1000);
-}
-
-BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_HP_FILTER_DOUBLE, System<double>)
-{
-    auto bf = difi::Butterworthd(order, fc, fs, difi::Butterworthd::Type::HighPass);
-    BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
-    test_coeffs(hpACoeffRes, hpBCoeffRes, bf, std::numeric_limits<double>::epsilon() * 10);
-    test_results(hpResults, data, bf, std::numeric_limits<double>::epsilon() * 100);
-}
-
-BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_BP_FILTER_FLOAT, System<float>)
-{
-    auto bf = difi::Butterworthf(order, fLower, fUpper, fs);
-    BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
-    test_coeffs(bpACoeffRes, bpBCoeffRes, bf, std::numeric_limits<float>::epsilon() * 1000);
-    test_results(bpResults, data, bf, std::numeric_limits<float>::epsilon() * 10000);
-}
-
-BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_BP_FILTER_DOUBLE, System<double>)
-{
-    auto bf = difi::Butterworthd(order, fLower, fUpper, fs);
-    BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
-    test_coeffs(bpACoeffRes, bpBCoeffRes, bf, std::numeric_limits<double>::epsilon() * 1000);
-    test_results(bpResults, data, bf, std::numeric_limits<double>::epsilon() * 10000);
-}
-
-BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_BR_FILTER_FLOAT, System<float>)
-{
-    auto bf = difi::Butterworthf(order, fLower, fUpper, fs, difi::Butterworthf::Type::BandReject);
-    BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
-    test_coeffs(brACoeffRes, brBCoeffRes, bf, 1.f);
-    test_results(brResults, data, bf, 1.f);
-}
-
-BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_BR_FILTER_DOUBLE, System<double>)
-{
-    auto bf = difi::Butterworthd(order, fLower, fUpper, fs, difi::Butterworthd::Type::BandReject);
-    BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
-    test_coeffs(brACoeffRes, brBCoeffRes, bf, 1e-8);
-    test_results(brResults, data, bf, 1e-8);
+    System<T> s;
+    auto bf = difi::Butterworth<T>(s.order, s.fLower, s.fUpper, s.fs, difi::Butterworth<T>::Type::BandReject);
+    REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
+    test_coeffs(s.brACoeffRes, s.brBCoeffRes, bf, std::numeric_limits<T>::epsilon() * T(1e8));
+    test_results(s.brResults, s.data, bf, std::numeric_limits<T>::epsilon() * T(1e8));
 }

tests/CMakeLists.txt

@@ -27,20 +27,13 @@
 
 enable_testing()
 
-if(${BUILD_TEST_STATIC_BOOST})
-    set(Boost_USE_STATIC_LIBS ON)
-    set(BUILD_SHARED_LIBS OFF)
-    set(BOOST_DEFS "")
-else()
-    set(Boost_USE_STATIC_LIBS OFF)
-    set(BUILD_SHARED_LIBS ON)
-    set(BOOST_DEFS Boost::dynamic_linking)
-endif()
-find_package(Boost REQUIRED COMPONENTS unit_test_framework)
-
 macro(addTest testName)
     add_executable(${testName} ${testName}.cpp)
-    target_link_libraries(${testName} PRIVATE Boost::unit_test_framework Boost::disable_autolinking ${BOOST_DEFS} ${PROJECT_NAME})
+    if (MSVC)
+        target_compile_definitions(${testName} PUBLIC _SILENCE_ALL_CXX17_DEPRECATION_WARNINGS)
+    endif()
+    target_compile_definitions(${testName} PUBLIC DOCTEST_CONFIG_IMPLEMENT_WITH_MAIN)
+    target_link_libraries(${testName} PUBLIC Eigen3::Eigen)
     # Adding a project configuration file (for MSVC only)
     generate_msvc_dot_user_file(${testName})
 
@@ -51,4 +44,7 @@ addTest(GenericFilterTests)
 addTest(polynome_functions_tests)
 addTest(DigitalFilterTests)
 addTest(MovingAverageFilterTests)
-addTest(ButterWorthFilterTests)
+addTest(ButterworthFilterTests)
+
+# Differentiators
+addTest(differentiator_tests)

tests/DigitalFilterTests.cpp

@@ -25,12 +25,10 @@
 // of the authors and should not be interpreted as representing official policies,
 // either expressed or implied, of the FreeBSD Project.
 
-#define BOOST_TEST_MODULE DigitalFilterTests
-
 #include "difi"
+#include "doctest/doctest.h"
 #include "test_functions.h"
 #include "warning_macro.h"
-#include <boost/test/unit_test.hpp>
 
 DISABLE_CONVERSION_WARNING_BEGIN
 
@@ -44,16 +42,10 @@ struct System {
 
 DISABLE_CONVERSION_WARNING_END
 
-BOOST_FIXTURE_TEST_CASE(DIGITAL_FILTER_FLOAT, System<float>)
+TEST_CASE_TEMPLATE("Digital filter", T, float, double)
 {
-    auto df = difi::DigitalFilterf(aCoeff, bCoeff);
-    test_coeffs(aCoeff, bCoeff, df, std::numeric_limits<float>::epsilon() * 10);
-    test_results(results, data, df, std::numeric_limits<float>::epsilon() * 10);
-}
-
-BOOST_FIXTURE_TEST_CASE(DIGITAL_FILTER_DOUBLE, System<double>)
-{
-    auto df = difi::DigitalFilterd(aCoeff, bCoeff);
-    test_coeffs(aCoeff, bCoeff, df, std::numeric_limits<double>::epsilon() * 10);
-    test_results(results, data, df, std::numeric_limits<double>::epsilon() * 10);
+    System<T> s;
+    auto df = difi::DigitalFilter<T>(s.aCoeff, s.bCoeff);
+    test_coeffs(s.aCoeff, s.bCoeff, df, std::numeric_limits<T>::epsilon() * 10);
+    test_results(s.results, s.data, df, std::numeric_limits<T>::epsilon() * 10);
 }

tests/GenericFilterTests.cpp

@@ -25,35 +25,52 @@
 // of the authors and should not be interpreted as representing official policies,
 // either expressed or implied, of the FreeBSD Project.
 
-#define BOOST_TEST_MODULE GenericFilterTests
-
 #include "difi"
-#include <boost/test/unit_test.hpp>
+#include "doctest/doctest.h"
 #include <exception>
 #include <vector>
 
-BOOST_AUTO_TEST_CASE(FILTER_FAILURES)
+TEST_CASE("Filter failures")
 {
     // A coeff are missing
-    BOOST_REQUIRE_THROW(difi::DigitalFilterd(Eigen::VectorXd(), Eigen::VectorXd::Constant(2, 0)), std::logic_error);
+    REQUIRE_THROWS_AS(difi::DigitalFilterd(Eigen::VectorXd(), Eigen::VectorXd::Constant(2, 0)), std::logic_error);
     // B coeff are missing
-    BOOST_REQUIRE_THROW(difi::DigitalFilterd(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd()), std::logic_error);
+    REQUIRE_THROWS_AS(difi::DigitalFilterd(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd()), std::logic_error);
     // aCoeff(0) = 0
-    BOOST_REQUIRE_THROW(difi::DigitalFilterd(Eigen::VectorXd::Constant(2, 0), Eigen::VectorXd::Constant(2, 0)), std::logic_error);
+    REQUIRE_THROWS_AS(difi::DigitalFilterd(Eigen::VectorXd::Constant(2, 0), Eigen::VectorXd::Constant(2, 0)), std::logic_error);
     // Filter left uninitialized
-    BOOST_REQUIRE_NO_THROW(difi::DigitalFilterd());
+    REQUIRE_NOTHROW(difi::DigitalFilterd());
     auto df = difi::DigitalFilterd();
     // Filter data with uninitialized filter
-    BOOST_REQUIRE_THROW(df.stepFilter(10.), std::logic_error);
+    REQUIRE_THROWS_AS(df.stepFilter(10.), std::logic_error);
+    // Change type before settings coeffs (The difi::FilterType::Centered needs to have odd number of bCoeffs)
+    REQUIRE_THROWS_AS(df.setType(difi::FilterType::Centered), std::logic_error);
+    // Set type to difi::FilterType::Backward is always ok.
+    REQUIRE_NOTHROW(df.setType(difi::FilterType::Backward));
+    // Set even number of bCoeffs
+    REQUIRE_NOTHROW(df.setCoeffs(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd::Constant(2, 0)));
+    // Check again
+    REQUIRE_THROWS_AS(df.setType(difi::FilterType::Centered), std::logic_error);
+    // Set odd number of bCoeffs
+    REQUIRE_NOTHROW(df.setCoeffs(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd::Constant(3, 0)));
+    // Check again
+    REQUIRE_NOTHROW(df.setType(difi::FilterType::Centered));
+    // Put even number of bCoeffs on a Centered filter
+    REQUIRE_THROWS_AS(df.setCoeffs(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd::Constant(2, 0)), std::logic_error);
+    // Put odd number of bCoeffs
+    REQUIRE_NOTHROW(df.setCoeffs(Eigen::VectorXd::Constant(5, 2), Eigen::VectorXd::Constant(7, 3)));
+
     // window <= 0
-    BOOST_REQUIRE_THROW(difi::MovingAveraged(0), std::logic_error);
+    REQUIRE_THROWS_AS(difi::MovingAveraged(0), std::logic_error);
     // order <= 0
-    BOOST_REQUIRE_THROW(difi::Butterworthd(0, 10, 100), std::logic_error);
+    REQUIRE_THROWS_AS(difi::Butterworthd(0, 10, 100), std::logic_error);
     // fc > 2*fs
-    BOOST_REQUIRE_THROW(difi::Butterworthd(2, 60, 100), std::logic_error);
+    REQUIRE_THROWS_AS(difi::Butterworthd(2, 60, 100), std::logic_error);
     // Upper frequency < lower frequency
-    BOOST_REQUIRE_THROW(difi::Butterworthd(2, 6, 5, 100), std::logic_error);
+    REQUIRE_THROWS_AS(difi::Butterworthd(2, 6, 5, 100), std::logic_error);
 
     // Ok
-    BOOST_REQUIRE_NO_THROW(difi::DigitalFilterd(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd::Constant(2, 0)));
+    REQUIRE_NOTHROW(difi::DigitalFilterd(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd::Constant(2, 0)));
+    // Bad type. Need odd number of bCoeffs
+    REQUIRE_THROWS_AS(difi::DigitalFilterd(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd::Constant(2, 0), difi::FilterType::Centered), std::logic_error);
 }

tests/MovingAverageFilterTests.cpp

@@ -25,11 +25,9 @@
 // of the authors and should not be interpreted as representing official policies,
 // either expressed or implied, of the FreeBSD Project.
 
-#define BOOST_TEST_MODULE MovingAverageFilterTests
-
 #include "difi"
+#include "doctest/doctest.h"
 #include "test_functions.h"
-#include <boost/test/unit_test.hpp>
 
 template <typename T>
 struct System {
@@ -38,14 +36,9 @@ struct System {
     difi::vectX_t<T> results = (difi::vectX_t<T>(6) << 0.25, 0.75, 1.5, 2.5, 3.5, 4.5).finished();
 };
 
-BOOST_FIXTURE_TEST_CASE(MOVING_AVERAGE_FLOAT, System<float>)
+TEST_CASE_TEMPLATE("Moving average filter", T, float, double)
 {
-    auto maf = difi::MovingAveragef(windowSize);
-    test_results(results, data, maf, std::numeric_limits<float>::epsilon() * 10);
-}
-
-BOOST_FIXTURE_TEST_CASE(MOVING_AVERAGE_DOUBLE, System<double>)
-{
-    auto maf = difi::MovingAveraged(windowSize);
-    test_results(results, data, maf, std::numeric_limits<double>::epsilon() * 10);
-}
+    System<T> s;
+    auto maf = difi::MovingAverage<T>(s.windowSize);
+    test_results(s.results, s.data, maf, std::numeric_limits<T>::epsilon() * 10);
+}

tests/diffTesters.h

@@ -0,0 +1,277 @@
+// Copyright (c) 2019, Vincent SAMY
+// All rights reserved.
+
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+
+// 1. Redistributions of source code must retain the above copyright notice,
+//    this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+// The views and conclusions contained in the software and documentation are those
+// of the authors and should not be interpreted as representing official policies,
+// either expressed or implied, of the FreeBSD Project.
+
+#pragma once
+#include "doctest_helper.h"
+#include "difi"
+#include <array>
+#include "doctest/doctest.h"
+
+using namespace difi;
+
+namespace tester {
+
+namespace details {
+
+template <typename Tuple, size_t Index>
+struct set_time_step {
+    static void impl(Tuple& filters, double dt)
+    {
+        std::get<Index>(filters).setTimestep(dt);
+        set_time_step<Tuple, Index - 1>::impl(filters, dt);
+    }
+};
+
+template <typename Tuple>
+struct set_time_step<Tuple, 0> {
+    static void impl(Tuple& filters, double dt)
+    {
+        std::get<0>(filters).setTimestep(dt);
+    }
+};
+
+// For TF tests
+template <typename Filter>
+void test_filter(Filter& filt, const vectX_t<double>& f, const vectX_t<double>& df, double eps)
+{
+    for (int i = 0; i < 50; ++i) {
+        filt.stepFilter(f(i)); // First initialize, some steps
+        // auto value = filt.stepFilter(f(i)); // First initialize, some steps
+        // value = 2.;
+    }
+
+    for (int i = 50; i < f.size(); ++i) {
+        auto value = filt.stepFilter(f(i));
+        REQUIRE_SMALL(std::abs(value - df(i - filt.center())), eps);
+    }
+}
+
+// For TV tests
+template <typename Filter>
+void test_filter(Filter& filt, const vectX_t<double>& t, const vectX_t<double>& f, const vectX_t<double>& df, double eps)
+{
+    for (int i = 0; i < 50; ++i) {
+        filt.stepFilter(t(i), f(i)); // First initialize, some steps
+        // auto value = filt.stepFilter(f(i)); // First initialize, some steps
+        // value = 2.;
+    }
+
+    for (int i = 50; i < f.size(); ++i) {
+        auto value = filt.stepFilter(t(i), f(i));
+        REQUIRE_SMALL(std::abs(value - df(i - filt.center())), eps);
+    }
+}
+
+template <typename Tuple, size_t Index>
+struct run_test {
+    static constexpr size_t ind = Index;
+
+    // For TF tests
+    static void impl(Tuple& filters, const vectX_t<double>& f, const vectX_t<double>& df, const vectX_t<double>& eps)
+    {
+        test_filter(std::get<Index>(filters), f, df, eps(ind));
+        run_test<Tuple, Index - 1>::impl(filters, f, df, eps);
+    }
+
+    // For TV tests
+    static void impl(Tuple& filters, const vectX_t<double>& t, const vectX_t<double>& f, const vectX_t<double>& df, const vectX_t<double>& eps)
+    {
+        test_filter(std::get<Index>(filters), t, f, df, eps(ind));
+        run_test<Tuple, Index - 1>::impl(filters, t, f, df, eps);
+    }
+};
+
+template <typename Tuple>
+struct run_test<Tuple, 0> {
+    // For TF tests
+    static void impl(Tuple& filters, const vectX_t<double>& f, const vectX_t<double>& df, const vectX_t<double>& eps)
+    {
+        test_filter(std::get<0>(filters), f, df, eps(0));
+    }
+
+    // For TV tests
+    static void impl(Tuple& filters, const vectX_t<double>& t, const vectX_t<double>& f, const vectX_t<double>& df, const vectX_t<double>& eps)
+    {
+        test_filter(std::get<0>(filters), t, f, df, eps(0));
+    }
+};
+
+} // namespace details
+
+
+template <typename Tuple>
+void set_time_steps(Tuple& filters, double dt)
+{
+    details::set_time_step<Tuple, std::tuple_size<Tuple>::value - 1>::impl(filters, dt);
+}
+
+template <typename Tuple>
+void run_tests(Tuple& filters, const vectX_t<double>& f, const vectX_t<double>& df, const vectX_t<double>& eps)
+{
+    details::run_test<Tuple, std::tuple_size<Tuple>::value - 1>::impl(filters, f, df, eps);
+}
+
+template <typename Tuple>
+void run_tests(Tuple& filters, const vectX_t<double>& t, const vectX_t<double>& f, const vectX_t<double>& df, const vectX_t<double>& eps)
+{
+    details::run_test<Tuple, std::tuple_size<Tuple>::value - 1>::impl(filters, t, f, df, eps);
+}
+
+// template <typename T, typename Tuple, typename Derived>
+// class DiffTester {
+//     friend Derived;
+//     static constexpr const size_t TSize = std::tuple_size_v<Tuple>;
+
+// public:
+//     DiffTester(const Tuple& filters, const vectX_t<T>& f, const vectX_t<T>& df)
+//         : m_filters(filters)
+//         , m_f(f)
+//         , m_df(df)
+//     {}
+
+//     void run(const std::array<T, TSize>& eps) { testRunner<TSize - 1>(eps); }
+
+// private:
+//     template <typename Filter>
+//     void test(Filter& filt, T eps)
+//     {
+//         static_cast<Derived&>(*this).testFilter(filt, eps);
+//     }
+
+//     template <size_t Index>
+//     void testRunner(const std::array<T, TSize>& eps)
+//     {
+//         test(std::get<Index>(m_filters), eps[Index]);
+//         testRunner<Index - 1>(eps);
+//     }
+
+//     template <>
+//     void testRunner<0>(const std::array<T, TSize>& eps) { test(std::get<0>(m_filters), eps[0]); }
+
+// private:
+//     vectX_t<T> m_f;
+//     vectX_t<T> m_df;
+//     Tuple m_filters;
+// };
+
+// template <typename T, typename Tuple>
+// class TFDiffTester : public DiffTester<T, Tuple, TFDiffTester<T, Tuple>> {
+//     friend DiffTester<T, Tuple, TFDiffTester<T, Tuple>>;
+
+// public:
+//     TFDiffTester(const Tuple& filters, const vectX_t<T>& f, const vectX_t<T>& df)
+//         : DiffTester(filters, f, df)
+//     {}
+
+//     void setFilterTimestep(T step) { setStep<TSize - 1>(step); }
+
+// private:
+//     template <typename Filter>
+//     void testFilter(Filter& filt, T eps)
+//     {
+//         for (int i = 0; i < 50; ++i) {
+//             auto value = filt.stepFilter(m_f(i)); // First initialize, some steps
+//             value = 2.;
+//         }
+
+//         for (int i = 50; i < STEPS; ++i) {
+//             auto value = filt.stepFilter(m_f(i));
+//             REQUIRE_SMALL(std::abs(value - m_df(i - filt.center())), eps);
+//         }
+//     }
+
+//     template <size_t Index>
+//     void setStep(T step)
+//     {
+//         std::get<Index>(m_filters).setTimestep(step);
+//         setStep<Index - 1>(step);
+//     }
+
+//     template <>
+//     void setStep<0>(T step) { std::get<0>(m_filters).setTimestep(step); }
+// };
+
+// template <typename T, typename Tuple>
+// class TVDiffTester : public DiffTester<T, Tuple, TVDiffTester<T, Tuple>> {
+//     friend DiffTester<T, Tuple, TVDiffTester<T, Tuple>>;
+
+// public:
+//     TVDiffTester(const Tuple& filters, const vectX_t<T>& t, const vectX_t<T>& f, const vectX_t<T>& df)
+//         : DiffTester(filters, f, df)
+//         , m_time(t)
+//     {}
+
+// private:
+//     template <typename Filter>
+//     void testFilter(Filter& filt, T eps)
+//     {
+//         for (int i = 0; i < 50; ++i) {
+//             auto value = filt.stepFilter(m_time(i), m_f(i)); // First initialize, some steps
+//             value = 2.;
+//         }
+
+//         for (int i = 50; i < STEPS; ++i) {
+//             auto value = filt.stepFilter(m_time(i), m_f(i));
+//             REQUIRE_SMALL(std::abs(value - m_df(i - filt.center())), eps);
+//         }
+//     }
+
+// private:
+//     vectX_t<T> m_time;
+// };
+
+template <int N>
+using central_list = std::tuple<CenteredDiffBasicd<N>, CenteredDiffLowNoiseLanczosd<N>, CenteredDiffSuperLowNoiseLanczosd<N>, CenteredDiffNoiseRobust2d<N>, CenteredDiffNoiseRobust4d<N>>;
+
+// template <size_t N>
+// TFDiffTester<double, central_list<N>> generateCTester(const vectX_t<double>& f, const vectX_t<double>& df)
+// {
+//     return { central_list<N>{}, f, df };
+// }
+
+// template <size_t N>
+// TFDiffTester<double, std::tuple<CenteredDiffSecondOrderd<N>>> generateC2OTester(const vectX_t<double>& f, const vectX_t<double>& df)
+// {
+//     return { std::tuple<CenteredDiffSecondOrderd<N>>{}, f, df };
+// }
+
+template <int N>
+using tv_central_list = std::tuple<TVCenteredDiffNoiseRobust2d<N>, TVCenteredDiffNoiseRobust4d<N>>;
+
+// template <size_t N>
+// TVDiffTester<double, tv_central_list<N>> generateTVCTester(const vectX_t<double>& t, const vectX_t<double>& f, const vectX_t<double>& df)
+// {
+//     return { tv_central_list<N>{}, t, f, df };
+// }
+
+// template <size_t N>
+// TVDiffTester<double, std::tuple<TVCenteredDiffSecondOrderd<N>>> generateTVC2OTester(const vectX_t<double>& t, const vectX_t<double>& f, const vectX_t<double>& df)
+// {
+//     return { std::tuple<TVCenteredSecondOrderDiff<N>>{}, t, f, df };
+// }
+
+} // namespace tester

tests/differentiator_tests.cpp

@@ -0,0 +1,205 @@
+// Copyright (c) 2019, Vincent SAMY
+// All rights reserved.
+
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+
+// 1. Redistributions of source code must retain the above copyright notice,
+//    this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+// The views and conclusions contained in the software and documentation are those
+// of the authors and should not be interpreted as representing official policies,
+// either expressed or implied, of the FreeBSD Project.
+
+#include "diffTesters.h"
+#include "doctest/doctest.h"
+#include "doctest_helper.h"
+#include "noisy_function_generator.h"
+#include <limits>
+
+constexpr const int STEPS = 500;
+constexpr const int SIN_AMPLITUDE = 100;
+constexpr const int SIN_FREQUENCY = 1;
+template <typename T> constexpr const std::array<T, 5> POLY_4 = { 2, -3, 3, 2, -5 };
+template <typename T> constexpr const std::array<T, 8> POLY_7 = { 10, -12, 5, 6, -9, -3, -2, 4 };
+
+template <int N>
+vectN_t<double, N> generateLNLCoeffs()
+{
+    static_assert(N > 2 && N % 2 == 1, "N must be odd");
+    vectN_t<double, N> lnl;
+    const int M = (N - 1) / 2;
+    for (int i = 0; i < M + 1; ++i) {
+        lnl(i) = static_cast<double>(M - i);
+        lnl(M + i) = -static_cast<double>(i);
+    }
+
+    return lnl;
+}
+
+template <int N>
+void checkCoeffs(const vectN_t<double, N>& coeffs, const vectN_t<double, N>& goodCoeffs)
+{
+    for (int i = 0; i < N; ++i)
+        REQUIRE_SMALL(std::abs(coeffs(i) - goodCoeffs(i)), std::numeric_limits<double>::epsilon() * 5);
+}
+
+TEST_CASE("Coefficient calculation") // Some coeffs are computed, the rest are given
+{
+    // LNL coeffs
+    checkCoeffs<5>(details::GetLNLCoeffs<double, 5>{}(), generateLNLCoeffs<5>() / 10.);
+    checkCoeffs<7>(details::GetLNLCoeffs<double, 7>{}(), generateLNLCoeffs<7>() / 28.);
+    checkCoeffs<9>(details::GetLNLCoeffs<double, 9>{}(), generateLNLCoeffs<9>() / 60.);
+    checkCoeffs<11>(details::GetLNLCoeffs<double, 11>{}(), generateLNLCoeffs<11>() / 110.);
+
+    // FNR coeffs
+    checkCoeffs<3>(details::GetFNRCoeffs<double, 3>{}(), (vectN_t<double, 3>() << 1., 0., -1.).finished() / 2.);
+    checkCoeffs<4>(details::GetFNRCoeffs<double, 4>{}(), (vectN_t<double, 4>() << 1., 1., -1., -1.).finished() / 4.);
+    checkCoeffs<5>(details::GetFNRCoeffs<double, 5>{}(), (vectN_t<double, 5>() << 1., 2., 0., -2., -1.).finished() / 8.);
+    checkCoeffs<6>(details::GetFNRCoeffs<double, 6>{}(), (vectN_t<double, 6>() << 1., 3., 2., -2., -3., -1.).finished() / 16.);
+    checkCoeffs<7>(details::GetFNRCoeffs<double, 7>{}(), (vectN_t<double, 7>() << 1., 4., 5., 0., -5., -4., -1.).finished() / 32.);
+    checkCoeffs<8>(details::GetFNRCoeffs<double, 8>{}(), (vectN_t<double, 8>() << 1., 5., 9., 5., -5., -9., -5., -1.).finished() / 64.);
+    checkCoeffs<9>(details::GetFNRCoeffs<double, 9>{}(), (vectN_t<double, 9>() << 1., 6., 14., 14., 0., -14., -14., -6., -1.).finished() / 128.);
+    checkCoeffs<10>(details::GetFNRCoeffs<double, 10>{}(), (vectN_t<double, 10>() << 1., 7., 20., 28., 14., -14., -28., -20., -7., -1.).finished() / 256.);
+    checkCoeffs<11>(details::GetFNRCoeffs<double, 11>{}(), (vectN_t<double, 11>() << 1., 8., 27., 48., 42., 0., -42., -48., -27., -8., -1.).finished() / 512.);
+}
+
+TEST_CASE("Sinus time-fixed central derivative")
+{
+    double dt = 0.001;
+    auto sg = sinGenerator<double>(STEPS, SIN_AMPLITUDE, SIN_FREQUENCY, dt);
+    auto ct7 = tester::central_list<7>{};
+    auto ct9 = tester::central_list<9>{};
+
+    tester::set_time_steps(ct7, dt);
+    tester::set_time_steps(ct9, dt);
+
+    difi::vectX_t<double> eps{ 5 };
+    eps << 1e-10, 1e-1, 1e-6, 1e-1, 1e-6; // Value checked with MATLAB
+
+    tester::run_tests(ct7, std::get<0>(sg), std::get<1>(sg), eps);
+    tester::run_tests(ct9, std::get<0>(sg), std::get<1>(sg), eps);
+}
+
+TEST_CASE("Polynome time-fixed central derivative")
+{
+    double dt = 0.001;
+    {
+        auto pg = polyGenerator<double>(STEPS, POLY_4<double>, dt);
+        auto ct7 = tester::central_list<7>{};
+        auto ct9 = tester::central_list<9>{};
+
+        tester::set_time_steps(ct7, dt);
+        tester::set_time_steps(ct9, dt);
+
+        difi::vectX_t<double> eps{ 5 };
+        eps << 1e-12, 1e-4, 1e-12, 1e-4, 1e-12; // Value checked with MATLAB
+
+        tester::run_tests(ct7, std::get<0>(pg), std::get<1>(pg), eps);
+        tester::run_tests(ct9, std::get<0>(pg), std::get<1>(pg), eps);
+    }
+    {
+        auto pg = polyGenerator<double>(STEPS, POLY_7<double>, dt);
+        auto ct7 = tester::central_list<7>{};
+        auto ct9 = tester::central_list<9>{};
+
+        tester::set_time_steps(ct7, dt);
+        tester::set_time_steps(ct9, dt);
+
+        difi::vectX_t<double> eps{ 5 };
+        eps << 1e-11, 1e-3, 1e-9, 1e-4, 1e-9; // Value checked with MATLAB
+
+        tester::run_tests(ct7, std::get<0>(pg), std::get<1>(pg), eps);
+        tester::run_tests(ct9, std::get<0>(pg), std::get<1>(pg), eps);
+    }
+}
+
+TEST_CASE("2nd order sinus time-fixed center derivative")
+{
+    double dt = 0.001;
+    auto sg = sinGenerator<double>(STEPS, SIN_AMPLITUDE, SIN_FREQUENCY, dt);
+    auto ct7 = std::tuple<CenteredDiffSecondOrderd<7>>{};
+    auto ct9 = std::tuple<CenteredDiffSecondOrderd<9>>{};
+    auto ct11 = std::tuple<CenteredDiffSecondOrderd<11>>{};
+
+    tester::set_time_steps(ct7, dt);
+    tester::set_time_steps(ct9, dt);
+    tester::set_time_steps(ct11, dt);
+
+    difi::vectX_t<double> eps{ 1 };
+    eps << 2e-1;
+
+    tester::run_tests(ct7, std::get<0>(sg), std::get<2>(sg), eps);
+    tester::run_tests(ct9, std::get<0>(sg), std::get<2>(sg), eps);
+    tester::run_tests(ct11, std::get<0>(sg), std::get<2>(sg), eps);
+}
+
+TEST_CASE("Sinus time-varying central derivative")
+{
+    double dt = 0.001;
+    auto sg = tvSinGenerator<double>(STEPS, SIN_AMPLITUDE, SIN_FREQUENCY, dt);
+    auto ct7 = tester::tv_central_list<7>{};
+    auto ct9 = tester::tv_central_list<9>{};
+
+    difi::vectX_t<double> eps{ 2 };
+    eps << 1., 1.;
+
+    tester::run_tests(ct7, std::get<0>(sg), std::get<1>(sg), std::get<2>(sg), eps);
+    tester::run_tests(ct9, std::get<0>(sg), std::get<1>(sg), std::get<2>(sg), eps);
+}
+
+TEST_CASE("Polynome time-varying central derivative")
+{
+    double dt = 0.001;
+    {
+        auto pg = tvPolyGenerator<double>(STEPS, POLY_4<double>, dt);
+        auto ct7 = tester::tv_central_list<7>{};
+        auto ct9 = tester::tv_central_list<9>{};
+
+        difi::vectX_t<double> eps{ 2 };
+        eps << 1e-3, 1e-3;
+
+        tester::run_tests(ct7, std::get<0>(pg), std::get<1>(pg), std::get<2>(pg), eps);
+        tester::run_tests(ct9, std::get<0>(pg), std::get<1>(pg), std::get<2>(pg), eps);
+    }
+    {
+        auto pg = tvPolyGenerator<double>(STEPS, POLY_7<double>, dt);
+        auto ct7 = tester::tv_central_list<7>{};
+        auto ct9 = tester::tv_central_list<9>{};
+
+        difi::vectX_t<double> eps{ 2 };
+        eps << 1e-2, 1e-2;
+
+        tester::run_tests(ct7, std::get<0>(pg), std::get<1>(pg), std::get<2>(pg), eps);
+        tester::run_tests(ct9, std::get<0>(pg), std::get<1>(pg), std::get<2>(pg), eps);
+    }
+}
+
+// TEST_CASE("2nd order sinus time-varying center derivative", "[tv][sin][center][2nd]")
+// {
+//     // Test not passing.
+//     double dt = 0.001;
+//     auto sg = tvSinGenerator<double>(STEPS, SIN_AMPLITUDE, SIN_FREQUENCY, dt);
+//     auto ct7 = generateTVC2OTester<double, 7>(std::get<0>(sg), std::get<1>(sg), std::get<3>(sg));
+//     auto ct9 = generateTVC2OTester<double, 9>(std::get<0>(sg), std::get<1>(sg), std::get<3>(sg));
+//     auto ct11 = generateTVC2OTester<double, 11>(std::get<0>(sg), std::get<1>(sg), std::get<3>(sg));
+
+//     std::array<double, 1> eps = { 2e-1 };
+//     ct7.run(eps);
+//     ct9.run(eps);
+//     ct11.run(eps);
+// }

tests/doctest_helper.h

@@ -0,0 +1,4 @@
+#pragma once
+
+#define REQUIRE_EQUAL(left, right) REQUIRE((left) == (right))
+#define REQUIRE_SMALL(value, eps) REQUIRE((value) < (eps))

tests/noisy_function_generator.h

@@ -0,0 +1,197 @@
+// Copyright (c) 2019, Vincent SAMY
+// All rights reserved.
+
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+
+// 1. Redistributions of source code must retain the above copyright notice,
+//    this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright notice,
+//    this list of conditions and the following disclaimer in the documentation
+//    and/or other materials provided with the distribution.
+
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+// ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+// WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+// ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+// The views and conclusions contained in the software and documentation are those
+// of the authors and should not be interpreted as representing official policies,
+// either expressed or implied, of the FreeBSD Project.
+
+#pragma once
+
+#include "typedefs.h"
+#include <array>
+#include <cmath>
+#include <tuple>
+
+template <typename T>
+using FunctionGenerator = std::tuple<difi::vectX_t<T>, difi::vectX_t<T>, difi::vectX_t<T>>;
+
+template <typename T>
+FunctionGenerator<T> sinGenerator(int nrSteps, T amplitude, T frequency, T dt)
+{
+    using namespace difi;
+
+    vectX_t<T> f(nrSteps);
+    vectX_t<T> df(nrSteps);
+    vectX_t<T> ddf(nrSteps);
+
+    for (int i = 0; i < nrSteps; ++i) {
+        // function
+        f(i) = amplitude * std::sin(2 * pi<T> * frequency * i * dt);
+
+        // derivative 1
+        df(i) = 2 * pi<T> * frequency * amplitude * std::cos(2 * pi<T> * frequency * i * dt);
+
+        // derivative 2
+        ddf(i) = -4 * pi<T> * pi<T> * frequency * frequency * amplitude * std::sin(2 * pi<T> * frequency * i * dt);
+    }
+
+    return { f, df, ddf };
+}
+
+template <typename T, size_t N>
+FunctionGenerator<T> polyGenerator(int nrSteps, std::array<T, N> coeffs, T dt)
+{
+    using namespace difi;
+    static_assert(N >= 2, "N must be superior to 20");
+
+    auto computePoly = [](const auto& coeffs, T time) {
+        auto recursiveComputation = [time, &coeffs](size_t i, T result, const auto& self) -> T {
+            if (i > 0)
+                return self(i - 1, time * result + coeffs[i - 1], self);
+            else
+                return result;
+        };
+
+        return recursiveComputation(coeffs.size(), 0, recursiveComputation);
+    };
+
+    std::array<T, N - 1> dCoeffs;
+    for (size_t i = 1; i < N; ++i)
+        dCoeffs[i - 1] = static_cast<T>(i) * coeffs[i];
+
+    std::array<T, N - 2> ddCoeffs;
+    for (size_t i = 1; i < N - 1; ++i)
+        ddCoeffs[i - 1] = static_cast<T>(i) * dCoeffs[i];
+
+    vectX_t<T> f(nrSteps);
+    vectX_t<T> df(nrSteps);
+    vectX_t<T> ddf(nrSteps);
+    for (int i = 0; i < nrSteps; ++i) {
+        // function
+        f(i) = computePoly(coeffs, i * dt);
+
+        // derivative 1
+        df(i) = computePoly(dCoeffs, i * dt);
+
+        // derivative 2
+        ddf(i) = computePoly(ddCoeffs, i * dt);
+    }
+
+    return { f, df, ddf };
+}
+
+// TV sin generator
+
+template <typename T>
+using TVFunctionGenerator = std::tuple<difi::vectX_t<T>, difi::vectX_t<T>, difi::vectX_t<T>, difi::vectX_t<T>>;
+
+template <typename T>
+TVFunctionGenerator<T> tvSinGenerator(int nrSteps, T amplitude, T frequency, T meanDt)
+{
+    using namespace difi;
+
+    vectX_t<T> t(nrSteps);
+    vectX_t<T> f(nrSteps);
+    vectX_t<T> df(nrSteps);
+    vectX_t<T> ddf(nrSteps);
+
+    t(0) = T(0);
+    f(0) = T(0);
+    df(0) = 2 * pi<T> * frequency * amplitude;
+    ddf(0) = T(0);
+    for (int i = 1; i < nrSteps; ++i) {
+        // time
+        if (i % 3 == 0)
+            t(i) = t(i - 1) + 0.9 * meanDt;
+        else if (i % 3 == 1)
+            t(i) = t(i - 1) + meanDt;
+        else
+            t(i) = t(i - 1) + 1.1 * meanDt;
+
+        // function
+        f(i) = amplitude * std::sin(2 * pi<T> * frequency * t(i));
+
+        // derivative 1
+        df(i) = 2 * pi<T> * frequency * amplitude * std::cos(2 * pi<T> * frequency * t(i));
+
+        // derivative 2
+        ddf(i) = -4 * pi<T> * pi<T> * frequency * frequency * amplitude * std::sin(2 * pi<T> * frequency * t(i));
+    }
+
+    return { t, f, df, ddf };
+}
+
+template <typename T, size_t N>
+TVFunctionGenerator<T> tvPolyGenerator(int nrSteps, std::array<T, N> coeffs, T meanDt)
+{
+    using namespace difi;
+    static_assert(N >= 2, "N must be superior to 20");
+
+    auto computePoly = [](const auto& coeffs, T time) {
+        auto recursiveComputation = [time, &coeffs](size_t i, T result, const auto& self) -> T {
+            if (i > 0)
+                return self(i - 1, time * result + coeffs[i - 1], self);
+            else
+                return result;
+        };
+
+        return recursiveComputation(coeffs.size(), 0, recursiveComputation);
+    };
+
+    std::array<T, N - 1> dCoeffs;
+    for (size_t i = 1; i < N; ++i)
+        dCoeffs[i - 1] = static_cast<T>(i) * coeffs[i];
+
+    std::array<T, N - 2> ddCoeffs;
+    for (size_t i = 1; i < N - 1; ++i)
+        ddCoeffs[i - 1] = static_cast<T>(i) * dCoeffs[i];
+
+    vectX_t<T> t(nrSteps);
+    vectX_t<T> f(nrSteps);
+    vectX_t<T> df(nrSteps);
+    vectX_t<T> ddf(nrSteps);
+    t(0) = T(0);
+    f(0) = computePoly(coeffs, t(0));
+    df(0) = computePoly(dCoeffs, t(0));
+    ddf(0) = computePoly(ddCoeffs, t(0));
+    for (int i = 1; i < nrSteps; ++i) {
+        // time
+        if (i % 3 == 0)
+            t(i) = t(i - 1) + 0.9 * meanDt;
+        else if (i % 3 == 1)
+            t(i) = t(i - 1) + meanDt;
+        else
+            t(i) = t(i - 1) + 1.1 * meanDt;
+
+        // function
+        f(i) = computePoly(coeffs, t(i));
+
+        // derivative 1
+        df(i) = computePoly(dCoeffs, t(i));
+
+        // derivative 2
+        ddf(i) = computePoly(ddCoeffs, t(i));
+    }
+
+    return { t, f, df, ddf };
+}

tests/polynome_functions_tests.cpp

@@ -25,11 +25,10 @@
 // of the authors and should not be interpreted as representing official policies,
 // either expressed or implied, of the FreeBSD Project.
 
-#define BOOST_TEST_MODULE polynome_functions_tests
-
 #include "difi"
+#include "doctest/doctest.h"
+#include "doctest_helper.h"
 #include "warning_macro.h"
-#include <boost/test/unit_test.hpp>
 #include <limits>
 
 using c_int_t = std::complex<int>;
@@ -62,50 +61,36 @@ struct SystemCFloat {
 
 DISABLE_CONVERSION_WARNING_END
 
-BOOST_FIXTURE_TEST_CASE(POLYNOME_FUNCTION_INT, SystemInt)
+TEST_CASE("Polynome function for int")
 {
-    auto res = difi::VietaAlgoi::polyCoeffFromRoot(data);
-
+    SystemInt s;
+    auto res = difi::VietaAlgoi::polyCoeffFromRoot(s.data);
     for (Eigen::Index i = 0; i < res.size(); ++i)
-        BOOST_REQUIRE_EQUAL(res(i), results(i));
+        REQUIRE_EQUAL(res(i), s.results(i));
 }
 
-BOOST_FIXTURE_TEST_CASE(POLYNOME_FUNCTION_FLOAT, SystemFloat<float>)
+TEST_CASE_TEMPLATE("Polynome function for floating point", T, float, double)
 {
-    auto res = difi::VietaAlgof::polyCoeffFromRoot(data);
+    SystemFloat<T> s;
+    auto res = difi::VietaAlgo<T>::polyCoeffFromRoot(s.data);
 
     for (Eigen::Index i = 0; i < res.size(); ++i)
-        BOOST_REQUIRE_SMALL(std::abs(res(i) - results(i)), std::numeric_limits<float>::epsilon() * 1000);
+        REQUIRE_SMALL(std::abs(res(i) - s.results(i)), std::numeric_limits<T>::epsilon() * 1000);
 }
 
-BOOST_FIXTURE_TEST_CASE(POLYNOME_FUNCTION_DOUBLE, SystemFloat<double>)
+TEST_CASE("Polynome function for complex int")
 {
-    auto res = difi::VietaAlgod::polyCoeffFromRoot(data);
+    SystemCInt s;
+    auto res = difi::VietaAlgoci::polyCoeffFromRoot(s.data);
+    for (Eigen::Index i = 0; i < res.size(); ++i)
+        REQUIRE_EQUAL(res(i), s.results(i));
+}
+
+TEST_CASE_TEMPLATE("Polynome function for complex floating point", T, float, double)
+{
+    SystemCFloat<T> s;
+    auto res = difi::VietaAlgo<std::complex<T>>::polyCoeffFromRoot(s.data);
 
     for (Eigen::Index i = 0; i < res.size(); ++i)
-        BOOST_REQUIRE_SMALL(std::abs(res(i) - results(i)), std::numeric_limits<double>::epsilon() * 1000);
-}
-
-BOOST_FIXTURE_TEST_CASE(POLYNOME_FUNCTION_CINT, SystemCInt)
-{
-    auto res = difi::VietaAlgoci::polyCoeffFromRoot(data);
-
-    for (Eigen::Index i = 0; i < res.size(); ++i)
-        BOOST_REQUIRE_EQUAL(res(i), results(i));
-}
-
-BOOST_FIXTURE_TEST_CASE(POLYNOME_FUNCTION_CFLOAT, SystemCFloat<float>)
-{
-    auto res = difi::VietaAlgocf::polyCoeffFromRoot(data);
-
-    for (Eigen::Index i = 0; i < res.size(); ++i)
-        BOOST_REQUIRE_SMALL(std::abs(res(i) - results(i)), std::numeric_limits<float>::epsilon() * 1000);
-}
-
-BOOST_FIXTURE_TEST_CASE(POLYNOME_FUNCTION_CDOUBLE, SystemCFloat<double>)
-{
-    auto res = difi::VietaAlgocd::polyCoeffFromRoot(data);
-
-    for (Eigen::Index i = 0; i < res.size(); ++i)
-        BOOST_REQUIRE_SMALL(std::abs(res(i) - results(i)), std::numeric_limits<double>::epsilon() * 1000);
+        REQUIRE_SMALL(std::abs(res(i) - s.results(i)), std::numeric_limits<T>::epsilon() * 1000);
 }

tests/test_functions.h

@@ -28,20 +28,21 @@
 #pragma once
 
 #include "difi"
-#include <boost/test/unit_test.hpp>
+#include "doctest/doctest.h"
+#include "doctest_helper.h"
 #include <limits>
 
 template <typename T>
 void test_coeffs(const difi::vectX_t<T>& aCoeff, const difi::vectX_t<T>& bCoeff, const difi::GenericFilter<T>& filter, T prec)
 {
-    BOOST_REQUIRE_EQUAL(aCoeff.size(), filter.aOrder());
-    BOOST_REQUIRE_EQUAL(bCoeff.size(), filter.bOrder());
+    REQUIRE_EQUAL(aCoeff.size(), filter.aOrder());
+    REQUIRE_EQUAL(bCoeff.size(), filter.bOrder());
     difi::vectX_t<T> faCoeff, fbCoeff;
     filter.getCoeffs(faCoeff, fbCoeff);
     for (Eigen::Index i = 0; i < faCoeff.size(); ++i)
-        BOOST_REQUIRE_SMALL(std::abs(aCoeff(i) - faCoeff(i)), prec);
+        REQUIRE_SMALL(std::abs(aCoeff(i) - faCoeff(i)), prec);
     for (Eigen::Index i = 0; i < fbCoeff.size(); ++i)
-        BOOST_REQUIRE_SMALL(std::abs(bCoeff(i) - fbCoeff(i)), prec);
+        REQUIRE_SMALL(std::abs(bCoeff(i) - fbCoeff(i)), prec);
 }
 
 template <typename T>
@@ -53,10 +54,10 @@ void test_results(const difi::vectX_t<T>& results, const difi::vectX_t<T>& data,
         filteredData(i) = filter.stepFilter(data(i));
 
     for (Eigen::Index i = 0; i < filteredData.size(); ++i)
-        BOOST_REQUIRE_SMALL(std::abs(filteredData(i) - results(i)), prec);
+        REQUIRE_SMALL(std::abs(filteredData(i) - results(i)), prec);
 
     filter.resetFilter();
     filteredData = filter.filter(data);
     for (Eigen::Index i = 0; i < filteredData.size(); ++i)
-        BOOST_REQUIRE_SMALL(std::abs(filteredData(i) - results(i)), prec);
+        REQUIRE_SMALL(std::abs(filteredData(i) - results(i)), prec);
 }

tests/warning_macro.h

@@ -27,18 +27,17 @@
 
 #pragma once
 
-#ifdef _MSC_VER
+#if defined(_MSC_VER)
 #define DCW_BEGIN           \
     __pragma(warning(push)) \
-        __pragma(warning(disable : 4305))
+    __pragma(warning(disable : 4305))
 #define DCW_END __pragma(warning(pop))
-#else
-#ifdef __GNUC__ || __clang__
+#elif defined(__GNUC__) || defined(__clang__)
 #define DCW_BEGIN               \
-    _Pragma("GCC warning push") \
-        _Pragma("GCC diagnostic ignored \"-Wconversion\"")
-#define DCW_END _Pragma("GCC warning pop")
-#endif
+    _Pragma("GCC diagnostic push") \
+    _Pragma("GCC diagnostic ignored \"-Wfloat-conversion\"") \
+    _Pragma("GCC diagnostic ignored \"-Wconversion\"") // Only needed for clang
+#define DCW_END _Pragma("GCC diagnostic pop")
 #endif
 
 #ifdef DCW_BEGIN