Update CMake and tests.

Closes #1
Closes #2

.clang-format

@@ -0,0 +1,110 @@
+---
+Language:        Cpp
+# BasedOnStyle:  WebKit
+AccessModifierOffset: -4
+AlignAfterOpenBracket: DontAlign
+AlignConsecutiveAssignments: false
+AlignConsecutiveDeclarations: false
+AlignEscapedNewlines: Right
+AlignOperands:   false
+AlignTrailingComments: false
+AllowAllParametersOfDeclarationOnNextLine: true
+AllowShortBlocksOnASingleLine: true
+AllowShortCaseLabelsOnASingleLine: true
+AllowShortFunctionsOnASingleLine: All
+AllowShortIfStatementsOnASingleLine: false
+AllowShortLambdasOnASingleLine: All
+AllowShortLoopsOnASingleLine: false
+AlwaysBreakAfterDefinitionReturnType: None
+AlwaysBreakAfterReturnType: None
+AlwaysBreakBeforeMultilineStrings: false
+AlwaysBreakTemplateDeclarations: No
+BinPackArguments: true
+BinPackParameters: true
+BraceWrapping:   
+  AfterClass:      false
+  AfterControlStatement: false
+  AfterEnum:       false
+  AfterFunction:   true
+  AfterNamespace:  false
+  AfterObjCDeclaration: false
+  AfterStruct:     false
+  AfterUnion:      false
+  BeforeCatch:     false
+  BeforeElse:      false
+  IndentBraces:    false
+  SplitEmptyFunction: false
+  SplitEmptyRecord: false
+  SplitEmptyNamespace: true
+BreakBeforeBinaryOperators: All
+BreakBeforeBraces: Custom
+BreakBeforeInheritanceComma: false
+BreakBeforeTernaryOperators: true
+BreakConstructorInitializersBeforeComma: false
+BreakConstructorInitializers: BeforeComma
+BreakAfterJavaFieldAnnotations: false
+BreakStringLiterals: true
+ColumnLimit:     0
+CommentPragmas:  '^ IWYU pragma:'
+CompactNamespaces: false
+ConstructorInitializerAllOnOneLineOrOnePerLine: false
+ConstructorInitializerIndentWidth: 4
+ContinuationIndentWidth: 4
+Cpp11BracedListStyle: false
+DerivePointerAlignment: false
+DisableFormat:   false
+ExperimentalAutoDetectBinPacking: false
+FixNamespaceComments: true
+ForEachMacros:   
+  - foreach
+  - Q_FOREACH
+  - BOOST_FOREACH
+IncludeBlocks: Preserve
+IncludeCategories: 
+  - Regex:           '^"(llvm|llvm-c|clang|clang-c)/'
+    Priority:        2
+  - Regex:           '^(<|"(gtest|gmock|isl|json)/)'
+    Priority:        3
+  - Regex:           '.*'
+    Priority:        1
+IncludeIsMainRegex: '(Test)?$'
+IndentCaseLabels: false
+IndentWidth:     4
+IndentWrappedFunctionNames: false
+JavaScriptQuotes: Leave
+JavaScriptWrapImports: true
+KeepEmptyLinesAtTheStartOfBlocks: true
+MacroBlockBegin: ''
+MacroBlockEnd:   ''
+MaxEmptyLinesToKeep: 1
+NamespaceIndentation: None
+ObjCBlockIndentWidth: 4
+ObjCSpaceAfterProperty: true
+ObjCSpaceBeforeProtocolList: true
+PenaltyBreakAssignment: 2
+PenaltyBreakBeforeFirstCallParameter: 19
+PenaltyBreakComment: 300
+PenaltyBreakFirstLessLess: 120
+PenaltyBreakString: 1000
+PenaltyExcessCharacter: 1000000
+PenaltyReturnTypeOnItsOwnLine: 60
+PointerAlignment: Left
+ReflowComments:  true
+SortIncludes:    true
+SortUsingDeclarations: true
+SpaceAfterCStyleCast: false
+SpaceAfterTemplateKeyword: true
+SpaceBeforeAssignmentOperators: true
+SpaceBeforeParens: ControlStatements
+SpaceInEmptyParentheses: false
+SpacesBeforeTrailingComments: 1
+SpacesInAngles:  false
+SpacesInContainerLiterals: true
+SpacesInCStyleCastParentheses: false
+SpacesInParentheses: false
+SpacesInSquareBrackets: false
+Standard:        Cpp11
+TabWidth:        8
+UseTab:          Never
+...
+

CMakeLists.txt

@@ -2,52 +2,51 @@
 # All rights reserved.
 
 # Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions are met:
-#     * Redistributions of source code must retain the above copyright
-#       notice, this list of conditions and the following disclaimer.
-#     * 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.
-#     * Neither the name of the AIST nor the
-#       names of its contributors may be used to endorse or promote products
-#       derived from this software without specific prior written permission.
+# 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# 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.
+
 # Version minimum
-cmake_minimum_required(VERSION 3.0.2)
+cmake_minimum_required(VERSION 3.8.2)
 
-include(cmake/base.cmake)
-include(cmake/eigen.cmake)
-include(cmake/msvc-specific.cmake)
-
-set(PROJECT_NAME DiFi++)
+set(PROJECT_NAME difi++)
 set(PROJECT_DESCRIPTION "Filter using rational transfer function")
-set(PROJECT_URL "...")
+set(PROJECT_URL "https://github.com/vsamy/DiFipp")
 SET(PROJECT_DEBUG_POSTFIX "_d")
+set(INSTALL_GENERATED_HEADERS OFF)
 
-#SET(CXX_DISABLE_WERROR True)
 set(DOXYGEN_USE_MATHJAX "NO")
+set(CMAKE_CXX_STANDARD 17)
 
-project(${PROJECT_NAME} CXX)
-set(CMAKE_CXX_STANDARD 14)
-
-setup_project()
-
-option(DISABLE_TESTS "Disable unit tests." OFF)
+option(BUILD_TESTING "Disable unit tests." ON)
 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)
 
+include(cmake/base.cmake)
+include(cmake/msvc-specific.cmake)
+
+project(${PROJECT_NAME} CXX)
+
 # Handle contracts specifications
 if(${THROW_ON_CONTRACT_VIOLATION})
     add_compile_options(-DGSL_THROW_ON_CONTRACT_VIOLATION)
@@ -59,21 +58,15 @@ endif()
 
 # for MSVC
 if(MSVC)
-    set(CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} /W1 /MP")
     set(CMAKE_MSVCIDE_RUN_PATH "\$(SolutionDir)/src/\$(Configuration)")
 endif()
 
-add_compile_options("-D_USE_MATH_DEFINES")
-
 # Eigen
 set(Eigen_REQUIRED "eigen3 >= 3.3")
-search_for_eigen()
+add_project_dependency(Eigen3 REQUIRED)
 
 add_subdirectory(include)
 
-if(NOT ${DISABLE_TESTS})
+if(${BUILD_TESTING})
     add_subdirectory(tests)
 endif()
-
-pkg_config_append_libs(${PROJECT_NAME})
-setup_project_finalize()

LICENSE

@@ -2,27 +2,28 @@ 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:
-    * Redistributions of source code must retain the above copyright
-      notice, this list of conditions and the following disclaimer.
-    * 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.
-    * Neither the name of the AIST nor the
-      names of its contributors may be used to endorse or promote products
-      derived from this software without specific prior written permission.
+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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+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.
 
 
 Note: GSL has a different license

README.md

@@ -1,16 +1,22 @@
 DiFi++
 ======
 
-[License BSD 3-Clause](https://tldrlegal.com/license/bsd-3-clause-license-(revised)#fulltext)
+[License BSD 2-Clause](https://tldrlegal.com/license/bsd-2-clause-license-(freebsd)#fulltext)
 
 DiFi++ is a small c++ header-only library  for **DI**gital **FI**lters based on rational transfer functions such as the butterworth filter and the moving average filter. DiFi++ is using the Eigen library for math computations.
 
 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 6ccc9f9b2a2ff510ae7754c515c72f8e38410447
+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/BilinearTransform.h

@@ -3,26 +3,28 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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/Butterworth.h

@@ -3,31 +3,32 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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 "typedefs.h"
-#include <cmath>
 #include <complex>
 
 namespace difi {
@@ -44,9 +45,6 @@ namespace difi {
  */
 template <typename T>
 class Butterworth : public DigitalFilter<T> {
-public:
-    static T PI; /*!< pi depending on the type T */
-
 public:
     /*! \brief Type of butterworth filter0 */
     enum class Type {
@@ -62,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

@@ -3,34 +3,33 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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 "BilinearTransform.h"
 #include "polynome_functions.h"
 
 namespace difi {
 
-template <typename T>
-T Butterworth<T>::PI = static_cast<T>(M_PI);
-
 template <typename T>
 std::pair<int, T> Butterworth<T>::findMinimumButter(T wPass, T wStop, T APass, T AStop)
 {
@@ -38,8 +37,8 @@ std::pair<int, T> Butterworth<T>::findMinimumButter(T wPass, T wStop, T APass, T
     Expects(wStop > T(0) && wPass < T(1));
     T num = std::log10((std::pow(T(10), T(0.1) * std::abs(AStop)) - 1) / (std::pow(T(10), T(0.1) * std::abs(APass)) - 1));
     // pre-warp
-    T fwPass = std::tan(T(0.5) * PI * wPass);
-    T fwStop = std::tan(T(0.5) * PI * wStop);
+    T fwPass = std::tan(T(0.5) * pi<T> * wPass);
+    T fwStop = std::tan(T(0.5) * pi<T> * wStop);
     T w;
     if (wPass < wStop)
         w = std::abs(fwStop / fwPass);
@@ -54,7 +53,7 @@ std::pair<int, T> Butterworth<T>::findMinimumButter(T wPass, T wStop, T APass, T
     else
         ctf = fwPass / ctf;
 
-    return std::pair<int, T>(order, T(2) * std::atan(ctf) / PI);
+    return std::pair<int, T>(order, T(2) * std::atan(ctf) / pi<T>);
 }
 
 template <typename T>
@@ -111,7 +110,7 @@ template <typename T>
 void Butterworth<T>::computeDigitalRep(T fc)
 {
     // Continuous pre-warped frequency
-    T fpw = (m_fs / PI) * std::tan(PI * fc / m_fs);
+    T fpw = (m_fs / pi<T>)*std::tan(pi<T> * fc / m_fs);
 
     // Compute poles
     vectXc_t<T> poles(m_order);
@@ -132,14 +131,14 @@ 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>
 void Butterworth<T>::computeBandDigitalRep(T fLower, T fUpper)
 {
-    T fpw1 = (m_fs / PI) * std::tan(PI * fLower / m_fs);
-    T fpw2 = (m_fs / PI) * std::tan(PI * fUpper / m_fs);
+    T fpw1 = (m_fs / pi<T>)*std::tan(pi<T> * fLower / m_fs);
+    T fpw2 = (m_fs / pi<T>)*std::tan(pi<T> * fUpper / m_fs);
     T fpw0 = std::sqrt(fpw1 * fpw2);
 
     vectXc_t<T> poles(2 * m_order);
@@ -161,26 +160,26 @@ void Butterworth<T>::computeBandDigitalRep(T fLower, T fUpper)
     }
 
     if (m_type == Type::BandPass)
-        scaleAmplitude(aCoeff, bCoeff, std::exp(std::complex<T>(T(0), T(2) * PI * std::sqrt(fLower * fUpper) / m_fs)));
+        scaleAmplitude(aCoeff, bCoeff, std::exp(std::complex<T>(T(0), T(2) * pi<T> * std::sqrt(fLower * fUpper) / m_fs)));
     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, order = m_order](int k) -> T {
-        return (2 * k - 1) * pi / (2 * order);
+    auto thetaK = [pi = pi<T>, order = m_order](int k) -> T {
+        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)));
     switch (m_type) {
     case Type::HighPass:
-        return T(2) * PI * fpw1 / analogPole;
+        return T(2) * pi<T> * fpw1 / analogPole;
     case Type::LowPass:
-        return T(2) * PI * fpw1 * analogPole;
+        return T(2) * pi<T> * fpw1 * analogPole;
     default:
         GSL_ASSUME(0);
     }
@@ -189,13 +188,13 @@ std::complex<T> Butterworth<T>::generateAnalogPole(int k, T fpw1)
 template <typename T>
 std::pair<std::complex<T>, std::complex<T>> Butterworth<T>::generateBandAnalogPole(int k, T fpw0, T bw)
 {
-    auto thetaK = [pi = PI, order = m_order](int k) -> T {
-        return (2 * k - 1) * pi / (2 * order);
+    auto thetaK = [pi = pi<T>, order = m_order](int k) -> T {
+        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)));
     std::pair<std::complex<T>, std::complex<T>> poles;
-    std::complex<T> s0 = T(2) * PI * fpw0;
+    std::complex<T> s0 = T(2) * pi<T> * fpw0;
     std::complex<T> s = T(0.5) * bw / fpw0;
     switch (m_type) {
     case Type::BandReject:
@@ -222,7 +221,7 @@ vectXc_t<T> Butterworth<T>::generateAnalogZeros(T fpw0)
     case Type::BandPass:
         return (vectXc_t<T>(2 * m_order) << vectXc_t<T>::Constant(m_order, std::complex<T>(-1)), vectXc_t<T>::Constant(m_order, std::complex<T>(1))).finished();
     case Type::BandReject: {
-        T w0 = T(2) * std::atan(PI * fpw0 / m_fs);
+        T w0 = T(2) * std::atan(pi<T> * fpw0 / m_fs);
         return (vectXc_t<T>(2 * m_order) << vectXc_t<T>::Constant(m_order, std::exp(std::complex<T>(0, w0))), vectXc_t<T>::Constant(m_order, std::exp(std::complex<T>(0, -w0)))).finished();
     }
     case Type::LowPass:

include/CMakeLists.txt

@@ -2,46 +2,56 @@
 # All rights reserved.
 
 # Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions are met:
-#     * Redistributions of source code must retain the above copyright
-#       notice, this list of conditions and the following disclaimer.
-#     * 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.
-#     * Neither the name of the AIST nor the
-#       names of its contributors may be used to endorse or promote products
-#       derived from this software without specific prior written permission.
+# 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# 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.
+
 set(HEADERS
-    gsl/gsl_assert.h
+    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
     typedefs.h
 )
 
+set(GSL_HEADERS gsl/gsl_assert.h)
+
 add_library(${PROJECT_NAME} INTERFACE)
-target_include_directories(${PROJECT_NAME} INTERFACE ${HEADERS})
+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}")
 install(TARGETS ${PROJECT_NAME}
+    EXPORT "${TARGETS_EXPORT_NAME}"
     RUNTIME DESTINATION bin
     LIBRARY DESTINATION lib
     ARCHIVE DESTINATION lib)
-install(FILES ${HEADERS} DESTINATION ${INCLUDE_INSTALL_DESTINATION})
+install(FILES ${HEADERS} DESTINATION ${INCLUDE_INSTALL_DESTINATION})
+install(FILES ${GSL_HEADERS} DESTINATION ${INCLUDE_INSTALL_DESTINATION}/gsl)

include/DigitalFilter.h

@@ -3,26 +3,28 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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 "GenericFilter.h"
@@ -43,9 +45,10 @@ 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)
     {
     }
 };

include/GenericFilter.h

@@ -3,48 +3,42 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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_checks.h"
-#include "typedefs.h"
-#include <stddef.h>
-#include <string>
+#include "BaseFilter.h"
 
 namespace difi {
 
-/*! \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.");
+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.
@@ -61,68 +55,59 @@ public:
      * \return Filtered signal.
      */
     vectX_t<T> filter(const vectX_t<T>& data);
-    /*! \brief Filter a signal and store in a user-defined Eigen vector.
-     * 
-     * Useful if the length of the signal is known in advance.
-     * \param[out] results Filtered signal.
-     * \param data Signal.
-     * \return False if vector's lengths do not match.
-     */
-    void getFilterResults(Eigen::Ref<vectX_t<T>> results, const vectX_t<T>& data);
-    /*! \brief Reset the data and filtered data. */
+
     void resetFilter() noexcept;
-    /*! \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 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; }
 
 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.
-     */
-    GenericFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff);
-    /*! \brief Default destructor. */
-    virtual ~GenericFilter() = default;
+    GenericFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff, FilterType type = FilterType::Backward)
+        : Base(aCoeff, bCoeff, 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.
+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.
      * 
-     * 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.
+     * 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 checkCoeffs(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff);
+    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:
+    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:
-    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

@@ -3,32 +3,30 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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.
 
-#include <limits>
+// 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.
 
 namespace difi {
 
-// Public functions
-
 template <typename T>
 T GenericFilter<T>::stepFilter(const T& data)
 {
@@ -40,27 +38,20 @@ 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[0] = m_bCoeff.dot(m_rawData) - m_aCoeff.dot(m_filteredData);
-    return m_filteredData[0];
+    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>
 vectX_t<T> GenericFilter<T>::filter(const vectX_t<T>& data)
-{
-    vectX_t<T> results(data.size());
-    getFilterResults(results, data);
-    return results;
-}
-
-template <typename T>
-void GenericFilter<T>::getFilterResults(Eigen::Ref<vectX_t<T>> results, const vectX_t<T>& data)
 {
     Expects(m_isInitialized);
-    Expects(results.size() == data.size());
+    vectX_t<T> results(data.size());
     for (Eigen::Index i = 0; i < data.size(); ++i)
         results(i) = stepFilter(data(i));
+    return results;
 }
 
 template <typename T>
@@ -71,58 +62,50 @@ void GenericFilter<T>::resetFilter() noexcept
 }
 
 template <typename T>
-template <typename T2>
-void GenericFilter<T>::setCoeffs(T2&& aCoeff, T2&& bCoeff)
+T TVGenericFilter<T>::stepFilter(const T& time, const T& data)
 {
-    static_assert(std::is_convertible_v<T2, vectX_t<T>>, "The coefficients types should be convertible to vectX_t<T>");
+    Expects(m_isInitialized);
 
-    checkCoeffs(aCoeff, bCoeff);
-    m_aCoeff = aCoeff;
-    m_bCoeff = bCoeff;
-    normalizeCoeffs();
-    resetFilter();
-    m_isInitialized = true;
+    // 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>
-void GenericFilter<T>::getCoeffs(vectX_t<T>& aCoeff, vectX_t<T>& bCoeff) const noexcept
+vectX_t<T> TVGenericFilter<T>::filter(const vectX_t<T>& data, const vectX_t<T>& time)
 {
-    aCoeff = m_aCoeff;
-    bCoeff = m_bCoeff;
-}
-
-// Protected functions
-
-template <typename T>
-GenericFilter<T>::GenericFilter(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff)
-    : m_aCoeff(aCoeff)
-    , m_bCoeff(bCoeff)
-    , m_filteredData(aCoeff.size())
-    , m_rawData(bCoeff.size())
-{
-    checkCoeffs(aCoeff, 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>::normalizeCoeffs()
+void TVGenericFilter<T>::resetFilter() noexcept
 {
-    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>
-void GenericFilter<T>::checkCoeffs(const vectX_t<T>& aCoeff, const vectX_t<T>& bCoeff)
-{
-    Expects(aCoeff.size() > 0);
-    Expects(std::abs(aCoeff[0]) > std::numeric_limits<T>::epsilon());
-    Expects(bCoeff.size() > 0);
+    m_filteredData.setZero(m_aCoeff.size());
+    m_rawData.setZero(m_bCoeff.size());
+    m_timers.setZero(m_bCoeff.size());
 }
 
 } // namespace difi

include/MovingAverage.h

@@ -3,29 +3,32 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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 "gsl/gsl_assert.h"
 #include "typedefs.h"
 
 namespace difi {
@@ -42,19 +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(); }
+    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

@@ -3,26 +3,28 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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 "BilinearTransform.h"
@@ -30,6 +32,7 @@
 #include "DigitalFilter.h"
 #include "GenericFilter.h"
 #include "MovingAverage.h"
+#include "differentiators.h"
 #include "polynome_functions.h"
 #include "typedefs.h"
 
@@ -57,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/polynome_functions.h

@@ -3,26 +3,28 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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_checks.h"

include/type_checks.h

@@ -3,26 +3,28 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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 <complex>

include/typedefs.h

@@ -2,37 +2,50 @@
 // All rights reserved.
 
 // Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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 <Eigen/Core>
 
 namespace difi {
 
+template <typename T>
+constexpr T pi = T(3.14159265358979323846264338327950288419716939937510582097494459230781L);
+
 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

@@ -3,34 +3,35 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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.
 
-#define BOOST_TEST_MODULE ButterworthFilterTests
+// 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.
 
 // 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
 
@@ -62,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

@@ -2,51 +2,49 @@
 # All rights reserved.
 
 # Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions are met:
-#     * Redistributions of source code must retain the above copyright
-#       notice, this list of conditions and the following disclaimer.
-#     * 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.
-#     * Neither the name of the AIST nor the
-#       names of its contributors may be used to endorse or promote products
-#       derived from this software without specific prior written permission.
+# 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# 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.
 
-include(${CMAKE_SOURCE_DIR}/cmake/boost.cmake)
+# 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.
 
 enable_testing()
 
-set(BOOST_COMPONENTS unit_test_framework)
-search_for_boost()
-
-add_definitions(-DBOOST_TEST_DYN_LINK)
-include_directories(${Boost_INCLUDE_DIRS})
-include_directories(${CMAKE_SOURCE_DIR}/include)
-link_directories(${Boost_LIBRARY_DIR} ${Boost_LIBRARY_DIRS})
-
 macro(addTest testName)
     add_executable(${testName} ${testName}.cpp)
-    if(NOT MSVC)
-        target_link_libraries(${testName} PUBLIC ${Boost_LIBRARIES})
+    if (MSVC)
+        target_compile_definitions(${testName} PUBLIC _SILENCE_ALL_CXX17_DEPRECATION_WARNINGS)
     endif()
-    add_test(${testName}Unit ${testName})
+    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})
+    generate_msvc_dot_user_file(${testName})
+
+    add_test(${testName}Unit ${testName})
 endmacro(addTest)
 
 addTest(GenericFilterTests)
 addTest(polynome_functions_tests)
 addTest(DigitalFilterTests)
 addTest(MovingAverageFilterTests)
-addTest(ButterWorthFilterTests)
+addTest(ButterworthFilterTests)
+
+# Differentiators
+addTest(differentiator_tests)

tests/DigitalFilterTests.cpp

@@ -3,32 +3,32 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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.
 
-#define BOOST_TEST_MODULE DigitalFilterTests
+// 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 "difi"
+#include "doctest/doctest.h"
 #include "test_functions.h"
 #include "warning_macro.h"
-#include <boost/test/unit_test.hpp>
 
 DISABLE_CONVERSION_WARNING_BEGIN
 
@@ -42,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

@@ -3,55 +3,74 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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.
 
-#define BOOST_TEST_MODULE GenericFilterTests
+// 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 "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

@@ -3,31 +3,31 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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.
 
-#define BOOST_TEST_MODULE MovingAverageFilterTests
+// 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 "difi"
+#include "doctest/doctest.h"
 #include "test_functions.h"
-#include <boost/test/unit_test.hpp>
 
 template <typename T>
 struct System {
@@ -36,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

@@ -3,31 +3,32 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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.
 
-#define BOOST_TEST_MODULE polynome_functions_tests
+// 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 "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>;
@@ -60,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

@@ -3,43 +3,46 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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 "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>
@@ -51,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

@@ -3,40 +3,41 @@
 
 // Redistribution and use in source and binary forms, with or without
 // modification, are permitted provided that the following conditions are met:
-//     * Redistributions of source code must retain the above copyright
-//       notice, this list of conditions and the following disclaimer.
-//     * 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.
-//     * Neither the name of the AIST nor the
-//       names of its contributors may be used to endorse or promote products
-//       derived from this software without specific prior written permission.
+
+// 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 <COPYRIGHT HOLDER> BE LIABLE FOR ANY
-// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+// 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
 
-#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