From 09a697bb5c359126e3b1ea7e28afda1c40d69ccc Mon Sep 17 00:00:00 2001 From: Vincent Samy Date: Tue, 29 Oct 2019 13:27:48 +0900 Subject: [PATCH] Change Base filter class. --- include/BaseFilter.h | 132 ++++++++++++++++++++++++++++++++++++++ include/BaseFilter.tpp | 96 +++++++++++++++++++++++++++ include/Butterworth.h | 2 +- include/DigitalFilter.h | 1 + include/GenericFilter.h | 116 +++++++++------------------------ include/GenericFilter.tpp | 101 +++++++++++++---------------- include/differentiators.h | 2 +- include/typedefs.h | 5 ++ 8 files changed, 311 insertions(+), 144 deletions(-) create mode 100644 include/BaseFilter.h create mode 100644 include/BaseFilter.tpp diff --git a/include/BaseFilter.h b/include/BaseFilter.h new file mode 100644 index 0000000..9d398ff --- /dev/null +++ b/include/BaseFilter.h @@ -0,0 +1,132 @@ +// 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 +#include + +namespace difi { + +// TODO: noexcept(Function of gsl variable) +// TODO: constructor with universal refs + +/*! \brief Low-level filter. + * + * It creates the basic and common functions of all linear filter that can written as a digital filter. + * This class can not be instantiated directly. + * + * \warning In Debug mode, all functions may throw if a filter is badly initialized. + * This not the case in Realese mode. + * + * \tparam T Floating type. + */ +template +class BaseFilter { + static_assert(std::is_floating_point::value && !std::is_const::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_center == 0 ? Type::Forward : Type::Centered); } + /*! \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& aCoeff, vectX_t& bCoeff) const noexcept; + /*! \brief Return coefficients of the denominator polynome. */ + const vectX_t& aCoeff() const noexcept { return m_aCoeff; } + /*! \brief Return coefficients of the numerator polynome. */ + const vectX_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; } + +protected: + /*! \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. + */ + template + void setCoeffs(T2&& aCoeff, T2&& bCoeff); + /*! \brief Normalized the filter coefficients such that aCoeff(0) = 1. */ + void normalizeCoeffs(); + /*! \brief Check for bad coefficients. + * + * Set the filter status to ready is everything is fine. + * \param aCoeff Denominator coefficients of the filter. + * \param bCoeff Numerator coefficients of the filter. + * \return True if the filter status is set on READY. + */ + bool checkCoeffs(const vectX_t& aCoeff, const vectX_t& bCoeff, Type type); + +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 center + */ + BaseFilter(const vectX_t& aCoeff, const vectX_t& bCoeff, FilterType type = FilterType::Forward); + /*! \brief Default destructor. */ + virtual ~BaseFilter() = default; + + Derived& derived() noexcept { return *static_cast(this); } + const Derived& derived() const noexcept { return *static_cast(this); } + +private: + Eigen::Index m_center = 0; /*!< Center of the filter. 0 is a one-sided filter. Default is 0. */ + bool m_isInitialized = false; /*!< Initialization state of the filter. Default is false */ + vectX_t m_aCoeff; /*!< Denominator coefficients of the filter */ + vectX_t m_bCoeff; /*!< Numerator coefficients of the filter */ + vectX_t m_filteredData; /*!< Last set of filtered data */ + vectX_t m_rawData; /*!< Last set of non-filtered data */ +}; + +} // namespace difi + +#include "BaseFilter.tpp" diff --git a/include/BaseFilter.tpp b/include/BaseFilter.tpp new file mode 100644 index 0000000..35d9b2a --- /dev/null +++ b/include/BaseFilter.tpp @@ -0,0 +1,96 @@ +// 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 + +namespace difi { + +// Public functions + +template +void BaseFilter::setType(Type type) +{ + Expects(type == Type::Centered ? m_bCoeff.size() > 2 && m_bCoeff.size() % 2 == 1 : true); + m_center = (type == Type::Forward ? 0 : (m_bCoeff.size() - 1) / 2); +} + +template +template +void BaseFilter::setCoeffs(T2&& aCoeff, T2&& bCoeff) +{ + static_assert(std::is_convertible_v>, "The coefficients types should be convertible to vectX_t"); + + Expects(checkCoeffs(aCoeff, bCoeff, (m_center == 0 ? Type::Forward : Type::Centered))); + m_aCoeff = aCoeff; + m_bCoeff = bCoeff; + normalizeCoeffs(); + resetFilter(); + m_isInitialized = true; +} + +template +void BaseFilter::getCoeffs(vectX_t& aCoeff, vectX_t& bCoeff) const noexcept +{ + aCoeff = m_aCoeff; + bCoeff = m_bCoeff; +} + +// Protected functions + +template +BaseFilter::BaseFilter(const vectX_t& aCoeff, const vectX_t& bCoeff, FilterType type) + : m_aCoeff(aCoeff) + , m_bCoeff(bCoeff) + , m_filteredData(aCoeff.size()) + , m_rawData(bCoeff.size()) +{ + Expects(checkCoeffs(aCoeff, bCoeff, type)); + m_center = (type == Type::Forward ? 0 : (bCoeff.size() - 1) / 2); + normalizeCoeffs(); + resetFilter(); + m_isInitialized = true; +} + +template +void BaseFilter::normalizeCoeffs() +{ + T a0 = m_aCoeff(0); + if (std::abs(a0 - T(1)) < std::numeric_limits::epsilon()) + return; + + m_aCoeff /= a0; + m_bCoeff /= a0; +} + +template +bool BaseFilter::checkCoeffs(const vectX_t& aCoeff, const vectX_t& bCoeff, FilterType type) +{ + bool centering = (type == Type::Centered ? (bCoeff.size() % 2 == 1) : true); + return aCoeff.size() > 0 && std::abs(aCoeff[0]) > std::numeric_limits::epsilon() && bCoeff.size() > 0 && centering; +} + +} // namespace difi \ No newline at end of file diff --git a/include/Butterworth.h b/include/Butterworth.h index 43dafb0..1b365fa 100644 --- a/include/Butterworth.h +++ b/include/Butterworth.h @@ -60,7 +60,7 @@ public: /*! \brief Function to help you design a Butterworth filter. * * It finds optimal values of the order and cut-off frequency. - * \warning Works only for low-pass and high-pass filters. + * \note Works only for low-pass and high-pass filters. * \see http://www.matheonics.com/Tutorials/Butterworth.html#Paragraph_3.2 * \see https://www.mathworks.com/help/signal/ref/buttord.html#d120e11079 * \param wPass Pass band edge. diff --git a/include/DigitalFilter.h b/include/DigitalFilter.h index f5dae97..c210408 100644 --- a/include/DigitalFilter.h +++ b/include/DigitalFilter.h @@ -50,6 +50,7 @@ public: : GenericFilter(aCoeff, bCoeff) { } + void setCoefficients(vectX_t&& aCoeff, vectX_t&& bCoeff) { setCoeffs(std::forward(aCoeff), std::forward(bCoeff)); diff --git a/include/GenericFilter.h b/include/GenericFilter.h index f8f85ce..9eb9566 100644 --- a/include/GenericFilter.h +++ b/include/GenericFilter.h @@ -27,37 +27,12 @@ #pragma once -#include "gsl/gsl_assert.h" -#include "type_checks.h" -#include "typedefs.h" -#include -#include +#include "BaseFilter.h" namespace difi { -// TODO: noexcept(Function of gsl variable) -// TODO: constructor with universal refs - -/*! \brief Low-level filter. - * - * It creates the basic and common functions of all linear filter that can written as a digital filter. - * This class can not be instantiated directly. - * - * \warning In Debug mode, all functions may throw if a filter is badly initialized. - * This not the case in Realese mode. - * - * \tparam T Floating type. - */ template -class GenericFilter { - static_assert(std::is_floating_point::value && !std::is_const::value, "Only accept non-complex floating point types."); - -public: - enum class Type { - OneSided, - Centered - }; - +class GenericFilter : BaseFilter> { public: /*! \brief Filter a new data. * @@ -73,75 +48,48 @@ public: * \return Filtered signal. */ vectX_t filter(const vectX_t& data); - /*! \brief Reset the data and filtered data. */ - void resetFilter() noexcept; - /*!< \brief Return the filter type */ - Type type() const noexcept { return (m_center == 0 ? Type::OneSided : Type::Centered); } - /*! \brief Get digital filter coefficients. - * - * 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& aCoeff, vectX_t& bCoeff) const noexcept; - /*! \brief Return coefficients of the denominator polynome. */ - const vectX_t& aCoeff() const noexcept { return m_aCoeff; } - /*! \brief Return coefficients of the numerator polynome. */ - const vectX_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; } + void resetFilter() noexcept protected: - /*! \brief Default uninitialized constructor. */ GenericFilter() = default; - /*! \brief Constructor. - * \param aCoeff Denominator coefficients of the filter in decreasing order. - * \param bCoeff Numerator coefficients of the filter in decreasing order. - * \param center - */ - GenericFilter(const vectX_t& aCoeff, const vectX_t& bCoeff, Type type = Type::OneSided); - /*! \brief Default destructor. */ - virtual ~GenericFilter() = default; + GenericFilter(const vectX_t& aCoeff, const vectX_t& bCoeff, Type type = Type::Forward) + : BaseFilter(aCoeff, bCoeff, type) + {} +}; - /*! \brief Set type of filter (one-sided or centered) +template +class TVGenericFilter : BaseFilter> { +public: + /*! \brief Filter a new data. * - * \param type The filter type. - * \warning bCoeff must be set before. + * 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 setType(Type type); - /*! \brief Set the new coefficients of the filters. - * - * It awaits a universal reference. - * \param aCoeff Denominator coefficients of the filter in decreasing order. - * \param bCoeff Numerator coefficients of the filter in decreasing order. - */ - template - void setCoeffs(T2&& aCoeff, T2&& bCoeff); - /*! \brief Normalized the filter coefficients such that aCoeff(0) = 1. */ - void normalizeCoeffs(); - /*! \brief Check for bad coefficients. + T stepFilter(const T& data, const T& time); + /*! \brief Filter a signal. * - * 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. + * Filter all data given by the signal. + * \param data Signal. + * \return Filtered signal. */ - bool checkCoeffs(const vectX_t& aCoeff, const vectX_t& bCoeff, Type type); + vectX_t filter(const vectX_t& data, const vectX_t& time); + +protected: + GenericFilter() = default; + GenericFilter(const vectX_t& aCoeff, const vectX_t& bCoeff, Type type = Type::Forward) + : BaseFilter(aCoeff, bCoeff, type) + { + m_timers.resize(m_bCoeffs.size()); + m_timeDiffs.resize(m_bCoeffs.size()); + } private: - Eigen::Index m_center = 0; /*!< Center of the filter. 0 is a one-sided filter. Default is 0. */ - bool m_isInitialized = false; /*!< Initialization state of the filter. Default is false */ - vectX_t m_aCoeff; /*!< Denominator coefficients of the filter */ - vectX_t m_bCoeff; /*!< Numerator coefficients of the filter */ - vectX_t m_filteredData; /*!< Last set of filtered data */ - vectX_t m_rawData; /*!< Last set of non-filtered data */ + vectX_t m_timers; + vectX_t m_timeDiffs; }; } // namespace difi -#include "GenericFilter.tpp" +#include "GenericFilter.tpp" \ No newline at end of file diff --git a/include/GenericFilter.tpp b/include/GenericFilter.tpp index 611a849..0e5b936 100644 --- a/include/GenericFilter.tpp +++ b/include/GenericFilter.tpp @@ -25,12 +25,8 @@ // of the authors and should not be interpreted as representing official policies, // either expressed or implied, of the FreeBSD Project. -#include - namespace difi { -// Public functions - template T GenericFilter::stepFilter(const T& data) { @@ -42,10 +38,10 @@ T GenericFilter::stepFilter(const T& data) for (Eigen::Index i = m_filteredData.size() - 1; i > 0; --i) m_filteredData(i) = m_filteredData(i - 1); - m_rawData[0] = data; - m_filteredData[0] = 0; - m_filteredData[m_center] = m_bCoeff.dot(m_rawData) - m_aCoeff.dot(m_filteredData); - return m_filteredData[m_center]; + m_rawData(0) = data; + m_filteredData(0) = 0; + m_filteredData(m_center) = m_bCoeff.dot(m_rawData) - m_aCoeff.dot(m_filteredData); + return m_filteredData(m_center); } template @@ -66,65 +62,54 @@ void GenericFilter::resetFilter() noexcept } template -void GenericFilter::setType(Type type) +T TVGenericFilter::stepFilter(const T& data, const T& time) { - Expects(type == Type::Centered ? m_bCoeff.size() > 2 && m_bCoeff.size() % 2 == 1 : true); - m_center = (type == Type::OneSided ? 0 : (m_bCoeff.size() - 1) / 2); + Expects(m_isInitialized); + + // Slide data (can't use SIMD, but should be small) + for (Eigen::Index i = m_rawData.size() - 1; i > 0; --i) { + m_rawData(i) = m_rawData(i - 1); + m_timers(i) = m_timers(i - 1); + } + for (Eigen::Index i = m_filteredData.size() - 1; i > 0; --i) + m_filteredData(i) = m_filteredData(i - 1); + + m_timers(0) = time; + if (m_center == 0) { + for (Eigen::Index i = m_rawData.size() - 1; i > 0; --i) + m_timeDiffs(i) = m_timeDiffs(i - 1); + m_timeDiffs(0) = time - m_timers(1); + } else { + const Eigen::Index S = data.size() - 1; + const Eigen::Index M = S / 2; + m_timeDiffs(M) = T(1); + for (Eigen::Index i = 0; i < M; ++i) + m_timeDiffs(i) = m_timers(i) - m_timers(S - i); + } + m_rawData(0) = data; + m_filteredData(0) = 0; + m_filteredData(m_center) = (m_bCoeff.cwiseQuotient(m_timeDiffs)).dot(m_rawData) - m_aCoeff.dot(m_filteredData); + return m_filteredData(m_center); } template -template -void GenericFilter::setCoeffs(T2&& aCoeff, T2&& bCoeff) +vectX_t TVGenericFilter::filter(const vectX_t& data, const vectX_t& time) { - static_assert(std::is_convertible_v>, "The coefficients types should be convertible to vectX_t"); - - Expects(checkCoeffs(aCoeff, bCoeff, (m_center == 0 ? Type::OneSided : Type::Centered))); - m_aCoeff = aCoeff; - m_bCoeff = bCoeff; - normalizeCoeffs(); - resetFilter(); - m_isInitialized = true; + Expects(m_isInitialized); + Expects(data.size() == time.size()); + vectX_t results(data.size()); + for (Eigen::Index i = 0; i < data.size(); ++i) + results(i) = stepFilter(data(i), time(i)); + return results; } template -void GenericFilter::getCoeffs(vectX_t& aCoeff, vectX_t& bCoeff) const noexcept +void GenericFilter::resetFilter() noexcept { - aCoeff = m_aCoeff; - bCoeff = m_bCoeff; -} - -// Protected functions - -template -GenericFilter::GenericFilter(const vectX_t& aCoeff, const vectX_t& bCoeff, Type type) - : m_aCoeff(aCoeff) - , m_bCoeff(bCoeff) - , m_filteredData(aCoeff.size()) - , m_rawData(bCoeff.size()) -{ - Expects(checkCoeffs(aCoeff, bCoeff, type)); - m_center = (type == Type::OneSided ? 0 : (bCoeff.size() - 1) / 2); - normalizeCoeffs(); - resetFilter(); - m_isInitialized = true; -} - -template -void GenericFilter::normalizeCoeffs() -{ - T a0 = m_aCoeff(0); - if (std::abs(a0 - T(1)) < std::numeric_limits::epsilon()) - return; - - m_aCoeff /= a0; - m_bCoeff /= a0; -} - -template -bool GenericFilter::checkCoeffs(const vectX_t& aCoeff, const vectX_t& bCoeff, Type type) -{ - bool centering = (type == Type::Centered ? (bCoeff.size() % 2 == 1) : true); - return aCoeff.size() > 0 && std::abs(aCoeff[0]) > std::numeric_limits::epsilon() && bCoeff.size() > 0 && centering; + m_filteredData.setZero(m_aCoeff.size()); + m_rawData.setZero(m_bCoeff.size()); + m_timers.setZero(m_bCoeff.size()); + m_timerDiffs.setZero(m_bCoeff.size()); } } // namespace difi \ No newline at end of file diff --git a/include/differentiators.h b/include/differentiators.h index cd5ac92..a4f4f91 100644 --- a/include/differentiators.h +++ b/include/differentiators.h @@ -226,7 +226,7 @@ public: T timestep() const noexcept { return std::pow(bCoeff()(0) / CoeffGetter()(0), T(1) / Order); } }; -template typename CoeffGetter> +template typename CoeffGetter> class CentralDifferentiator : public GenericFilter { public: CentralDifferentiator() diff --git a/include/typedefs.h b/include/typedefs.h index d67f479..f85ac96 100644 --- a/include/typedefs.h +++ b/include/typedefs.h @@ -43,4 +43,9 @@ using vectN_t = Eigen::Matrix; /*!< Fixed Eigen column-vector */ template using vectXc_t = vectX_t>; /*!< Eigen complex column-vector */ +enum class FilterType { + Forward, + Centered +}; + } // namespace difi \ No newline at end of file