kopia lustrzana https://github.com/vsamy/DiFipp
Almost updated all merging errors.
rodzic
ef7b6dfba1
commit
3a6bc791ff
|
@ -1,6 +1,7 @@
|
|||
#pragma once
|
||||
|
||||
#include <vector>
|
||||
#include "type_checks.h"
|
||||
#include "typedefs.h"
|
||||
|
||||
namespace fratio {
|
||||
|
||||
|
@ -10,9 +11,9 @@ struct BilinearTransform {
|
|||
static_assert(std::is_floating_point<SubType>::value, "This struct can only accept floating point types (real and complex).");
|
||||
|
||||
static void SToZ(SubType fs, const T& sPlanePole, T& zPlanePole);
|
||||
static void SToZ(SubType fs, const std::vector<T>& sPlanePoles, std::vector<T>& zPlanePoles);
|
||||
static void SToZ(SubType fs, const Eigen::VectorX<T>& sPlanePoles, Eigen::VectorX<T>& zPlanePoles); // Can be optimized
|
||||
static void ZToS(SubType fs, const T& zPlanePole, T& sPlanePole);
|
||||
static void ZToS(SubType fs, const std::vector<T>& zPlanePoles, std::vector<T>& sPlanePoles);
|
||||
static void ZToS(SubType fs, const Eigen::VectorX<T>& zPlanePoles, Eigen::VectorX<T>& sPlanePoles); // Can be optimized
|
||||
};
|
||||
|
||||
template <typename T>
|
||||
|
@ -23,7 +24,7 @@ void BilinearTransform<T>::SToZ(SubType fs, const T& sPlanePole, T& zPlanePole)
|
|||
}
|
||||
|
||||
template <typename T>
|
||||
void BilinearTransform<T>::SToZ(SubType fs, const std::vector<T>& sPlanePoles, std::vector<T>& zPlanePoles)
|
||||
void BilinearTransform<T>::SToZ(SubType fs, const Eigen::VectorX<T>& sPlanePoles, Eigen::VectorX<T>& zPlanePoles)
|
||||
{
|
||||
assert(sPlanePoles.size() == zPlanePoles.size());
|
||||
for (size_t k = 0; k < sPlanePoles.size(); ++k)
|
||||
|
@ -38,7 +39,7 @@ void BilinearTransform<T>::ZToS(SubType fs, const T& zPlanePole, T& sPlanePole)
|
|||
}
|
||||
|
||||
template <typename T>
|
||||
void BilinearTransform<T>::ZToS(SubType fs, const std::vector<T>& zPlanePoles, std::vector<T>& sPlanePoles)
|
||||
void BilinearTransform<T>::ZToS(SubType fs, const Eigen::VectorX<T>& zPlanePoles, Eigen::VectorX<T>& sPlanePoles)
|
||||
{
|
||||
assert(zPlanePoles.size() == sPlanePoles.size());
|
||||
for (size_t k = 0; k < sPlanePoles.size(); ++k)
|
||||
|
|
|
@ -2,6 +2,7 @@
|
|||
|
||||
#include "GenericFilter.h"
|
||||
#include "typedefs.h"
|
||||
#include <cmath>
|
||||
#include <complex>
|
||||
|
||||
namespace fratio {
|
||||
|
@ -31,7 +32,7 @@ private:
|
|||
void computeDigitalRep();
|
||||
void updateCoeffSize();
|
||||
std::complex<T> generateAnalogPole(T fpw, size_t k);
|
||||
std::vector<std::complex<T>> generateAnalogZeros();
|
||||
Eigen::VectorX<std::complex<T>> generateAnalogZeros();
|
||||
void scaleAmplitude();
|
||||
|
||||
private:
|
||||
|
|
|
@ -1,145 +0,0 @@
|
|||
#include "BilinearTransform.h"
|
||||
#include "polynome_functions.h"
|
||||
#include <cmath>
|
||||
#include <sstream>
|
||||
|
||||
namespace fratio {
|
||||
|
||||
template <typename T>
|
||||
Butterworth<T>::Butterworth(Type type)
|
||||
: m_type(type)
|
||||
{
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
Butterworth<T>::Butterworth(size_t order, T fc, T fs, Type type)
|
||||
: m_type(type)
|
||||
{
|
||||
initialize(order, fc, fs);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void Butterworth<T>::setFilterParameters(size_t order, T fc, T fs)
|
||||
{
|
||||
initialize(order, fc, fs);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void Butterworth<T>::initialize(size_t order, T fc, T fs)
|
||||
{
|
||||
if (m_fc > m_fs / 2.) {
|
||||
std::stringstream ss;
|
||||
ss << "The cut-off-frequency must be inferior to the sampling frequency"
|
||||
<< "\n Given cut-off-frequency is " << m_fc
|
||||
<< "\n Given sample frequency is " << m_fs;
|
||||
throw std::runtime_error(ss.str());
|
||||
}
|
||||
|
||||
m_order = order;
|
||||
m_fc = fc;
|
||||
m_fs = fs;
|
||||
updateCoeffSize();
|
||||
computeDigitalRep();
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void Butterworth<T>::computeDigitalRep()
|
||||
{
|
||||
// Continuous pre-warped frequency
|
||||
T fpw = (m_fs / PI) * std::tan(PI * m_fc / m_fs);
|
||||
|
||||
// Compute poles
|
||||
std::complex<T> analogPole;
|
||||
std::vector<std::complex<T>> poles(m_order);
|
||||
for (size_t k = 1; k <= m_order; ++k) {
|
||||
analogPole = generateAnalogPole(fpw, k);
|
||||
BilinearTransform<std::complex<T>>::SToZ(m_fs, analogPole, poles[k - 1]);
|
||||
}
|
||||
|
||||
std::vector<std::complex<T>> zeros = generateAnalogZeros();
|
||||
std::vector<std::complex<T>> a = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(poles);
|
||||
std::vector<std::complex<T>> b = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(zeros);
|
||||
for (size_t i = 0; i < m_order + 1; ++i) {
|
||||
m_aCoeff[i] = a[i].real();
|
||||
m_bCoeff[i] = b[i].real();
|
||||
}
|
||||
|
||||
scaleAmplitude();
|
||||
checkCoeff(m_aCoeff, m_bCoeff);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void Butterworth<T>::updateCoeffSize()
|
||||
{
|
||||
m_aCoeff.resize(m_order + 1);
|
||||
m_bCoeff.resize(m_order + 1);
|
||||
resetFilter();
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
std::complex<T> Butterworth<T>::generateAnalogPole(T fpw, size_t k)
|
||||
{
|
||||
T scaleFactor = 2 * PI * fpw;
|
||||
|
||||
auto thetaK = [pi = PI, order = m_order](size_t k) -> T {
|
||||
return (2 * k - 1) * pi / (2 * order);
|
||||
};
|
||||
|
||||
std::complex<T> analogPole(-std::sin(thetaK(k)), std::cos(thetaK(k)));
|
||||
switch (m_type) {
|
||||
case Type::HighPass:
|
||||
return scaleFactor / analogPole;
|
||||
|
||||
case Type::LowPass:
|
||||
default:
|
||||
return scaleFactor * analogPole;
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
std::vector<std::complex<T>> Butterworth<T>::generateAnalogZeros()
|
||||
{
|
||||
switch (m_type) {
|
||||
case Type::HighPass:
|
||||
return std::vector<std::complex<T>>(m_order, std::complex<T>(1));
|
||||
|
||||
case Type::LowPass:
|
||||
default:
|
||||
return std::vector<std::complex<T>>(m_order, std::complex<T>(-1));
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void Butterworth<T>::scaleAmplitude()
|
||||
{
|
||||
T scale = 0;
|
||||
T sumB = 0;
|
||||
|
||||
switch (m_type) {
|
||||
case Type::HighPass:
|
||||
for (size_t i = 0; i < m_order + 1; ++i) {
|
||||
if (i % 2 == 0) {
|
||||
scale += m_aCoeff[i];
|
||||
sumB += m_bCoeff[i];
|
||||
} else {
|
||||
scale -= m_aCoeff[i];
|
||||
sumB -= m_bCoeff[i];
|
||||
}
|
||||
}
|
||||
break;
|
||||
|
||||
case Type::LowPass:
|
||||
default:
|
||||
for (size_t i = 0; i < m_order + 1; ++i) {
|
||||
scale += m_aCoeff[i];
|
||||
sumB += m_bCoeff[i];
|
||||
}
|
||||
break;
|
||||
}
|
||||
|
||||
scale /= sumB;
|
||||
for (auto& b : m_bCoeff)
|
||||
b *= scale;
|
||||
}
|
||||
|
||||
} // namespace fratio
|
|
@ -1,6 +1,5 @@
|
|||
#include "BilinearTransform.h"
|
||||
#include "polynome_functions.h"
|
||||
#include <cmath>
|
||||
#include <sstream>
|
||||
|
||||
namespace fratio {
|
||||
|
||||
|
@ -34,7 +33,6 @@ void Butterworth<T>::initialize(size_t order, T fc, T fs)
|
|||
m_order = order;
|
||||
m_fc = fc;
|
||||
m_fs = fs;
|
||||
m_poles.resize(order);
|
||||
updateCoeffSize();
|
||||
computeDigitalRep();
|
||||
}
|
||||
|
@ -42,39 +40,27 @@ void Butterworth<T>::initialize(size_t order, T fc, T fs)
|
|||
template <typename T>
|
||||
void Butterworth<T>::computeDigitalRep()
|
||||
{
|
||||
T pi = static_cast<T>(M_PI);
|
||||
// Continuous pre-warped frequency
|
||||
T fpw = (m_fs / pi) * std::tan(pi * m_fc / m_fs);
|
||||
T scaleFactor = T(2) * pi * fpw;
|
||||
|
||||
auto thetaK = [pi, order = m_order](size_t k) -> T {
|
||||
return (T(2) * k - T(1)) * pi / (T(2) * order);
|
||||
};
|
||||
T fpw = (m_fs / PI) * std::tan(PI * m_fc / m_fs);
|
||||
|
||||
// Compute poles
|
||||
std::complex<T> scalePole;
|
||||
std::complex<T> analogPole;
|
||||
Eigen::VectorX<std::complex<T>> poles(m_order);
|
||||
for (size_t k = 1; k <= m_order; ++k) {
|
||||
scalePole = scaleFactor * std::complex<T>(-std::sin(thetaK(k)), std::cos(thetaK(k)));
|
||||
scalePole /= T(2) * m_fs;
|
||||
m_poles(k - 1) = (T(1) + scalePole) / (T(1) - scalePole);
|
||||
analogPole = generateAnalogPole(fpw, k);
|
||||
BilinearTransform<std::complex<T>>::SToZ(m_fs, analogPole, poles[k - 1]);
|
||||
}
|
||||
|
||||
Eigen::VectorX<std::complex<T>> numPoles = Eigen::VectorX::Constant(m_order, std::complex<T>(-1));
|
||||
Eigen::VectorX<std::complex<T>> a = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(m_poles);
|
||||
Eigen::VectorX<std::complex<T>> b = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(numPoles);
|
||||
T norm = 0;
|
||||
T sumB = 0;
|
||||
Eigen::VectorX<std::complex<T>> zeros = generateAnalogZeros();
|
||||
Eigen::VectorX<std::complex<T>> a = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(poles);
|
||||
Eigen::VectorX<std::complex<T>> b = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(zeros);
|
||||
for (size_t i = 0; i < m_order + 1; ++i) {
|
||||
m_aCoeff(i) = a(i).real();
|
||||
m_bCoeff(i) = b(i).real();
|
||||
norm += m_aCoeff(i);
|
||||
sumB += m_bCoeff(i);
|
||||
m_aCoeff[i] = a[i].real();
|
||||
m_bCoeff[i] = b[i].real();
|
||||
}
|
||||
|
||||
norm /= sumB;
|
||||
m_bCoeff *= norm;
|
||||
|
||||
checkCoeffs(m_aCoeff, m_bCoeff);
|
||||
scaleAmplitude();
|
||||
checkCoeff(m_aCoeff, m_bCoeff);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
|
@ -85,4 +71,66 @@ void Butterworth<T>::updateCoeffSize()
|
|||
resetFilter();
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
std::complex<T> Butterworth<T>::generateAnalogPole(T fpw, size_t k)
|
||||
{
|
||||
T scaleFactor = 2 * PI * fpw;
|
||||
|
||||
auto thetaK = [pi = PI, order = m_order](size_t k) -> T {
|
||||
return (2 * k - 1) * pi / (2 * order);
|
||||
};
|
||||
|
||||
std::complex<T> analogPole(-std::sin(thetaK(k)), std::cos(thetaK(k)));
|
||||
switch (m_type) {
|
||||
case Type::HighPass:
|
||||
return scaleFactor / analogPole;
|
||||
|
||||
case Type::LowPass:
|
||||
default:
|
||||
return scaleFactor * analogPole;
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
Eigen::VectorX<std::complex<T>> Butterworth<T>::generateAnalogZeros()
|
||||
{
|
||||
switch (m_type) {
|
||||
case Type::HighPass:
|
||||
return Eigen::VectorX<std::complex<T>>::Constant(m_order, std::complex<T>(1));
|
||||
|
||||
case Type::LowPass:
|
||||
default:
|
||||
return Eigen::VectorX<std::complex<T>>::Constant(m_order, std::complex<T>(-1));
|
||||
}
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void Butterworth<T>::scaleAmplitude()
|
||||
{
|
||||
T scale = 0;
|
||||
T sumB = 0;
|
||||
|
||||
switch (m_type) {
|
||||
case Type::HighPass:
|
||||
for (size_t i = 0; i < m_order + 1; ++i) {
|
||||
if (i % 2 == 0) {
|
||||
scale += m_aCoeff(i);
|
||||
sumB += m_bCoeff(i);
|
||||
} else {
|
||||
scale -= m_aCoeff(i);
|
||||
sumB -= m_bCoeff(i);
|
||||
}
|
||||
}
|
||||
break;
|
||||
|
||||
case Type::LowPass:
|
||||
default:
|
||||
scale = m_aCoeff.sum();
|
||||
sumB = m_bCoeff.sum();
|
||||
break;
|
||||
}
|
||||
|
||||
m_bCoeff *= scale / sumB;
|
||||
}
|
||||
|
||||
} // namespace fratio
|
|
@ -1,9 +1,11 @@
|
|||
set(HEADERS
|
||||
BilinearTransform.h
|
||||
Butterworth.h
|
||||
DigitalFilter.h
|
||||
GenericFilter.h
|
||||
MovingAverage.h
|
||||
polynome_functions.h
|
||||
typedefs.h
|
||||
)
|
||||
|
||||
add_library(${PROJECT_NAME} INTERFACE)
|
||||
|
|
|
@ -2,7 +2,6 @@
|
|||
|
||||
#include "GenericFilter.h"
|
||||
#include "typedefs.h"
|
||||
#include <vector>
|
||||
|
||||
namespace fratio {
|
||||
|
||||
|
@ -10,10 +9,6 @@ template <typename T>
|
|||
class DigitalFilter : public GenericFilter<T> {
|
||||
public:
|
||||
DigitalFilter() = default;
|
||||
DigitalFilter(const std::vector<T>& aCoeff, const std::vector<T>& bCoeff)
|
||||
{
|
||||
setCoeffs(aCoeff, bCoeff);
|
||||
}
|
||||
DigitalFilter(const Eigen::VectorX<T>& aCoeff, const Eigen::VectorX<T>& bCoeff)
|
||||
: GenericFilter<T>(aCoeff, bCoeff)
|
||||
{
|
||||
|
@ -23,6 +18,4 @@ public:
|
|||
size_t bOrder() const noexcept { return m_bCoeff.size(); }
|
||||
};
|
||||
|
||||
} // namespace fratio
|
||||
|
||||
#include "DigitalFilter.tpp"
|
||||
} // namespace fratio
|
|
@ -4,7 +4,6 @@
|
|||
#include "typedefs.h"
|
||||
#include <stddef.h>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
namespace fratio {
|
||||
|
||||
|
@ -22,9 +21,7 @@ public:
|
|||
bool getFilterResults(Eigen::Ref<Eigen::VectorX<T>> results, const Eigen::VectorX<T>& data);
|
||||
void resetFilter();
|
||||
|
||||
bool setCoeffs(const std::vector<T>& aCoeff, const std::vector<T>& bCoeff);
|
||||
bool setCoeffs(const Eigen::VectorX<T>& aCoeff, const Eigen::VectorX<T>& bCoeff);
|
||||
void getCoeffs(std::vector<T>& aCoeff, std::vector<T>& bCoeff) const;
|
||||
void getCoeffs(Eigen::Ref<Eigen::VectorX<T>> aCoeff, Eigen::Ref<Eigen::VectorX<T>> bCoeff) const;
|
||||
FilterStatus status() const noexcept { return m_status; }
|
||||
|
||||
|
@ -34,8 +31,7 @@ protected:
|
|||
virtual ~GenericFilter() = default;
|
||||
|
||||
void normalizeCoeffs();
|
||||
template <typename T2>
|
||||
bool checkCoeffs(const T2& aCoeff, const T2& bCoeff);
|
||||
bool checkCoeffs(const Eigen::VectorX<T>& aCoeff, const Eigen::VectorX<T>& bCoeff);
|
||||
|
||||
protected:
|
||||
FilterStatus m_status;
|
||||
|
|
|
@ -1,95 +0,0 @@
|
|||
namespace fratio {
|
||||
|
||||
// Public functions
|
||||
|
||||
template <typename T>
|
||||
T GenericFilter<T>::stepFilter(T data)
|
||||
{
|
||||
// Slide data
|
||||
for (auto rit1 = m_rawData.rbegin(), rit2 = m_rawData.rbegin() + 1; rit2 != m_rawData.rend(); ++rit1, ++rit2)
|
||||
*rit1 = *rit2;
|
||||
for (auto rit1 = m_filteredData.rbegin(), rit2 = m_filteredData.rbegin() + 1; rit2 != m_filteredData.rend(); ++rit1, ++rit2)
|
||||
*rit1 = *rit2;
|
||||
|
||||
T filtData = 0.;
|
||||
m_rawData[0] = data;
|
||||
|
||||
for (size_t k = 0; k < m_bCoeff.size(); ++k)
|
||||
filtData += m_bCoeff[k] * m_rawData[k];
|
||||
for (size_t k = 1; k < m_aCoeff.size(); ++k)
|
||||
filtData -= m_aCoeff[k] * m_filteredData[k];
|
||||
|
||||
m_filteredData[0] = filtData;
|
||||
|
||||
return filtData;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
std::vector<T> GenericFilter<T>::filter(const std::vector<T>& data)
|
||||
{
|
||||
std::vector<T> results;
|
||||
results.reserve(data.size());
|
||||
|
||||
for (T d : data)
|
||||
results.emplace_back(stepFilter(d));
|
||||
|
||||
return results;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void GenericFilter<T>::resetFilter()
|
||||
{
|
||||
m_filteredData.assign(m_aCoeff.size(), 0);
|
||||
m_rawData.assign(m_bCoeff.size(), 0);
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void GenericFilter<T>::getCoeff(std::vector<T>& aCoeff, std::vector<T>& bCoeff) const noexcept
|
||||
{
|
||||
aCoeff = m_aCoeff;
|
||||
bCoeff = m_bCoeff;
|
||||
}
|
||||
|
||||
// Protected functions
|
||||
|
||||
template <typename T>
|
||||
GenericFilter<T>::GenericFilter(const std::vector<T>& aCoeff, const std::vector<T>& bCoeff)
|
||||
: m_aCoeff(aCoeff)
|
||||
, m_bCoeff(bCoeff)
|
||||
, m_filteredData(aCoeff.size(), 0)
|
||||
, m_rawData(bCoeff.size(), 0)
|
||||
{
|
||||
checkCoeff(aCoeff, bCoeff);
|
||||
normalize();
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void GenericFilter<T>::normalize()
|
||||
{
|
||||
T a0 = m_aCoeff.front();
|
||||
if (std::abs(a0) < 1e-8) // Divide by zero
|
||||
throw std::invalid_argument("By filtering value for coefficient a0. Should be superior to 1e-8");
|
||||
|
||||
if (std::abs(a0 - 1) < 1e-8)
|
||||
return;
|
||||
|
||||
for (T& a : m_aCoeff)
|
||||
a /= a0;
|
||||
for (T& b : m_bCoeff)
|
||||
b /= a0;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void GenericFilter<T>::checkCoeff(const std::vector<T>& aCoeff, const std::vector<T>& bCoeff)
|
||||
{
|
||||
std::stringstream err;
|
||||
if (aCoeff.size() == 0)
|
||||
err << "The size of coefficient 'a' should greater than 0\n";
|
||||
if (bCoeff.size() == 0)
|
||||
err << "The size of coefficient 'b' should greater than 0\n";
|
||||
|
||||
if (err.str().size() > 0)
|
||||
throw std::runtime_error(err.str());
|
||||
}
|
||||
|
||||
} // namespace fratio
|
|
@ -32,10 +32,8 @@ T GenericFilter<T>::stepFilter(const T& data)
|
|||
assert(m_status == FilterStatus::READY);
|
||||
|
||||
// Slide data (can't use SIMD, but should be small)
|
||||
for (auto rit1 = m_rawData.rbegin(), rit2 = m_rawData.rbegin() + 1; rit2 != m_rawData.rend(); ++rit1, ++rit2)
|
||||
*rit1 = *rit2;
|
||||
for (auto rit1 = m_filteredData.rbegin(), rit2 = m_filteredData.rbegin() + 1; rit2 != m_filteredData.rend(); ++rit1, ++rit2)
|
||||
*rit1 = *rit2;
|
||||
m_rawData.tail(m_rawData.size() - 1) = m_rawData.head(m_rawData.size() - 1);
|
||||
m_filteredData.tail(m_rawData.size() - 1) = m_filteredData.head(m_rawData.size() - 1);
|
||||
|
||||
m_rawData[0] = data;
|
||||
m_filteredData[0] = m_bCoeff.dot(m_rawData) - m_aCoeff.dot(m_filteredData);
|
||||
|
@ -59,9 +57,8 @@ bool GenericFilter<T>::getFilterResults(Eigen::Ref<Eigen::VectorX<T>> results, c
|
|||
if (results.size() != data.size())
|
||||
return false;
|
||||
|
||||
T* res = results.data();
|
||||
for (T d : data)
|
||||
*(res++) = stepFilter(d);
|
||||
for (Eigen::Index i = 0; i < data.size(); ++i)
|
||||
results(i) = stepFilter(data(i));
|
||||
|
||||
return true;
|
||||
}
|
||||
|
@ -74,20 +71,7 @@ void GenericFilter<T>::resetFilter()
|
|||
}
|
||||
|
||||
template <typename T>
|
||||
bool GenericFilter<T>::setCoeffs(const std::vector<T>& aCoeff, const std::vector<T>& bCoeff)
|
||||
{
|
||||
if (!checkCoeffs(aCoeff, bCoeff))
|
||||
return false;
|
||||
|
||||
m_aCoeff = Eigen::Map<Eigen::VectorX<T>>(aCoeff.data(), aCoeff.size());
|
||||
m_bCoeff = Eigen::Map<Eigen::VectorX<T>>(bCoeff.data(), bCoeff.size());
|
||||
resetFilter();
|
||||
normalizeCoeffs();
|
||||
return true;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void GenericFilter<T>::setCoeffs(const Eigen::VectorX<T>& aCoeff, const Eigen::VectorX<T>& bCoeff)
|
||||
bool GenericFilter<T>::setCoeffs(const Eigen::VectorX<T>& aCoeff, const Eigen::VectorX<T>& bCoeff)
|
||||
{
|
||||
if (!checkCoeffs(aCoeff, bCoeff))
|
||||
return false;
|
||||
|
@ -99,13 +83,6 @@ void GenericFilter<T>::setCoeffs(const Eigen::VectorX<T>& aCoeff, const Eigen::V
|
|||
return true;
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void GenericFilter<T>::getCoeffs(std::vector<T>& aCoeff, std::vector<T>& bCoeff) const
|
||||
{
|
||||
aCoeff.assign(m_aCoeff.data(), m_aCoeff.data() + m_aCoeff.size());
|
||||
bCoeff.assign(m_bCoeff.data(), m_bCoeff.data() + m_bCoeff.size());
|
||||
}
|
||||
|
||||
template <typename T>
|
||||
void GenericFilter<T>::getCoeffs(Eigen::Ref<Eigen::VectorX<T>> aCoeff, Eigen::Ref<Eigen::VectorX<T>> bCoeff) const
|
||||
{
|
||||
|
@ -122,7 +99,7 @@ GenericFilter<T>::GenericFilter(const Eigen::VectorX<T>& aCoeff, const Eigen::Ve
|
|||
, m_filteredData(aCoeff.size())
|
||||
, m_rawData(bCoeff.size())
|
||||
{
|
||||
if(!checkCoeffs(aCoeff, bCoeff))
|
||||
if (!checkCoeffs(aCoeff, bCoeff))
|
||||
return;
|
||||
|
||||
resetFilter();
|
||||
|
@ -143,24 +120,21 @@ void GenericFilter<T>::normalizeCoeffs()
|
|||
}
|
||||
|
||||
template <typename T>
|
||||
template <typename T2>
|
||||
bool GenericFilter<T>::checkCoeffs(const T2& aCoeff, const T2& bCoeff)
|
||||
bool GenericFilter<T>::checkCoeffs(const Eigen::VectorX<T>& aCoeff, const Eigen::VectorX<T>& bCoeff)
|
||||
{
|
||||
using namespace FilterStatus;
|
||||
|
||||
m_status = NONE;
|
||||
m_status = FilterStatus::NONE;
|
||||
if (aCoeff.size() == 0)
|
||||
m_status = A_COEFF_MISSING;
|
||||
m_status = FilterStatus::A_COEFF_MISSING;
|
||||
else if (std::abs(aCoeff[0]) < std::numeric_limits<T>::epsilon())
|
||||
m_status = BAD_A_COEFF;
|
||||
m_status = FilterStatus::BAD_A_COEFF;
|
||||
|
||||
if (bCoeff.size() == 0)
|
||||
m_status = (m_status == A_COEFF_MISSING ? ALL_COEFF_MISSING : B_COEFF_MISSING);
|
||||
m_status = (m_status == FilterStatus::A_COEFF_MISSING ? FilterStatus::ALL_COEFF_MISSING : FilterStatus::B_COEFF_MISSING);
|
||||
|
||||
if (m_status == NONE)
|
||||
m_status = READY;
|
||||
if (m_status == FilterStatus::NONE)
|
||||
m_status = FilterStatus::READY;
|
||||
|
||||
return m_status == READY;
|
||||
return m_status == FilterStatus::READY;
|
||||
}
|
||||
|
||||
} // namespace fratio
|
|
@ -10,13 +10,7 @@ class MovingAverage : public DigitalFilter<T> {
|
|||
public:
|
||||
MovingAverage() = default;
|
||||
MovingAverage(size_t windowSize)
|
||||
: DigitalFilter<T>(Eigen::VectorX<T>::Constant(1, T(1)), Eigen::VectorX<T>::Constant(windowSize, T(1) / windowSize)))
|
||||
{
|
||||
}
|
||||
|
||||
void setWindowSize(size_t windowSize) { setCoeffs(Eigen::VectorX<T>::Constant(1, T(1)), Eigen::VectorX<T>::Constant(windowSize, T(1) / windowSize)); }
|
||||
size_t windowSize() const noexcept { return bOrder(); }
|
||||
: DigitalFilter<T>(Eigen::VectorX<T>::Constant(1, T(1)), Eigen::VectorX<T>::Constant(windowSize, T(1) / windowSize)))
|
||||
: DigitalFilter<T>(Eigen::VectorX<T>::Constant(1, T(1)), Eigen::VectorX<T>::Constant(windowSize, T(1) / windowSize))
|
||||
{
|
||||
}
|
||||
|
||||
|
|
|
@ -15,18 +15,18 @@ using MovingAveraged = MovingAverage<double>;
|
|||
using Butterworthf = Butterworth<float>;
|
||||
using Butterworthd = Butterworth<double>;
|
||||
|
||||
// Polynome helper functions
|
||||
using VietaAlgof = VietaAlgo<float>;
|
||||
using VietaAlgod = VietaAlgo<double>;
|
||||
using VietaAlgoi = VietaAlgo<int>;
|
||||
using VietaAlgocf = VietaAlgo<std::complex<float>>;
|
||||
using VietaAlgocd = VietaAlgo<std::complex<double>>;
|
||||
using VietaAlgoci = VietaAlgo<std::complex<int>>;
|
||||
// // Polynome helper functions
|
||||
// using VietaAlgof = VietaAlgo<float>;
|
||||
// using VietaAlgod = VietaAlgo<double>;
|
||||
// using VietaAlgoi = VietaAlgo<int>;
|
||||
// using VietaAlgocf = VietaAlgo<std::complex<float>>;
|
||||
// using VietaAlgocd = VietaAlgo<std::complex<double>>;
|
||||
// using VietaAlgoci = VietaAlgo<std::complex<int>>;
|
||||
|
||||
// Bilinear transformation functions
|
||||
using BilinearTransformf = BilinearTransform<float>;
|
||||
using BilinearTransformd = BilinearTransform<double>;
|
||||
using BilinearTransformcf = BilinearTransform<std::complex<float>>;
|
||||
using BilinearTransformcd = BilinearTransform<std::complex<double>>;
|
||||
// // Bilinear transformation functions
|
||||
// using BilinearTransformf = BilinearTransform<float>;
|
||||
// using BilinearTransformd = BilinearTransform<double>;
|
||||
// using BilinearTransformcf = BilinearTransform<std::complex<float>>;
|
||||
// using BilinearTransformcd = BilinearTransform<std::complex<double>>;
|
||||
|
||||
} // namespace fratio
|
|
@ -20,8 +20,8 @@ macro(addTest testName)
|
|||
GENERATE_MSVC_DOT_USER_FILE(${testName})
|
||||
endmacro(addTest)
|
||||
|
||||
addTest(ButterWorthFilterTests)
|
||||
addTest(GenericFilterTests)
|
||||
addTest(DigitalFilterTests)
|
||||
addTest(MovingAverageFilterTests)
|
||||
addTest(polynome_functions_tests)
|
||||
addTest(ButterWorthFilterTests)
|
||||
# addTest(polynome_functions_tests)
|
|
@ -1,6 +1,7 @@
|
|||
#define BOOST_TEST_MODULE DigitalFilterTests
|
||||
|
||||
#include "fratio.h"
|
||||
#include "typedefs.h"
|
||||
#include "warning_macro.h"
|
||||
#include <boost/test/unit_test.hpp>
|
||||
|
||||
|
@ -8,10 +9,10 @@ DISABLE_CONVERSION_WARNING_BEGIN
|
|||
|
||||
template <typename T>
|
||||
struct System {
|
||||
std::vector<T> data = { 1, 2, 3, 4 };
|
||||
std::vector<T> aCoeff = { 1, -0.99993717 };
|
||||
std::vector<T> bCoeff = { 0.99996859, -0.99996859 };
|
||||
std::vector<T> results = { 0.99996859, 1.999874351973491, 2.999717289867956, 3.999497407630634 };
|
||||
Eigen::VectorX<T> data = (Eigen::VectorX<T>(4) << 1, 2, 3, 4).finished();
|
||||
Eigen::VectorX<T> aCoeff = (Eigen::VectorX<T>(4) << 1, -0.99993717).finished();
|
||||
Eigen::VectorX<T> bCoeff = (Eigen::VectorX<T>(4) << 0.99996859, -0.99996859).finished();
|
||||
Eigen::VectorX<T> results = (Eigen::VectorX<T>(4) << 0.99996859, 1.999874351973491, 2.999717289867956, 3.999497407630634).finished();
|
||||
};
|
||||
|
||||
DISABLE_CONVERSION_WARNING_END
|
||||
|
@ -24,16 +25,16 @@ BOOST_FIXTURE_TEST_CASE(DIGITAL_FILTER_FLOAT, System<float>)
|
|||
|
||||
std::vector<float> filteredData;
|
||||
|
||||
for (float d : data)
|
||||
filteredData.push_back(df.stepFilter(d));
|
||||
for (Eigen::Index i = 0; i < data.size(); ++i)
|
||||
filteredData.push_back(df.stepFilter(data(i)));
|
||||
|
||||
for (size_t i = 0; i < filteredData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results[i]), 1e-6f);
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results(i)), 1e-6f);
|
||||
|
||||
df.resetFilter();
|
||||
filteredData = df.filter(data);
|
||||
for (size_t i = 0; i < filteredData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results[i]), 1e-6f);
|
||||
Eigen::VectorXf fData = df.filter(data);
|
||||
for (Eigen::Index i = 0; i < fData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(fData(i) - results(i)), 1e-6f);
|
||||
}
|
||||
|
||||
BOOST_FIXTURE_TEST_CASE(DIGITAL_FILTER_DOUBLE, System<double>)
|
||||
|
@ -44,14 +45,14 @@ BOOST_FIXTURE_TEST_CASE(DIGITAL_FILTER_DOUBLE, System<double>)
|
|||
|
||||
std::vector<double> filteredData;
|
||||
|
||||
for (double d : data)
|
||||
filteredData.push_back(df.stepFilter(d));
|
||||
for (Eigen::Index i = 0; i < data.size(); ++i)
|
||||
filteredData.push_back(df.stepFilter(data(i)));
|
||||
|
||||
for (size_t i = 0; i < filteredData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results[i]), 1e-14);
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results(i)), 1e-14);
|
||||
|
||||
df.resetFilter();
|
||||
filteredData = df.filter(data);
|
||||
for (size_t i = 0; i < filteredData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results[i]), 1e-14);
|
||||
Eigen::VectorXd fData = df.filter(data);
|
||||
for (Eigen::Index i = 0; i < fData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(fData(i) - results(i)), 1e-14);
|
||||
}
|
||||
|
|
|
@ -1,16 +1,22 @@
|
|||
#define BOOST_TEST_MODULE GenericFilterTests
|
||||
|
||||
#include "fratio.h"
|
||||
#include "typedefs.h"
|
||||
#include <boost/test/unit_test.hpp>
|
||||
#include <vector>
|
||||
|
||||
BOOST_AUTO_TEST_CASE(FilterThrows)
|
||||
{
|
||||
BOOST_REQUIRE_THROW(fratio::DigitalFilterd({}, { 1, 2 }), std::runtime_error);
|
||||
BOOST_REQUIRE_THROW(fratio::DigitalFilterd({ 1, 2 }, {}), std::runtime_error);
|
||||
BOOST_REQUIRE_THROW(fratio::DigitalFilterd({ 0 }, { 1 }), std::invalid_argument);
|
||||
|
||||
auto df = fratio::DigitalFilterd();
|
||||
BOOST_REQUIRE_THROW(df.setCoeff({}, { 1, 2 }), std::runtime_error);
|
||||
BOOST_REQUIRE_THROW(df.setCoeff({ 1, 2 }, {}), std::runtime_error);
|
||||
BOOST_REQUIRE_THROW(df.setCoeff({ 0 }, { 1 }), std::invalid_argument);
|
||||
auto dfd = fratio::DigitalFilterd(Eigen::VectorXd(), Eigen::VectorXd::Constant(2, 0));
|
||||
BOOST_REQUIRE(dfd.status() == fratio::FilterStatus::A_COEFF_MISSING);
|
||||
dfd = fratio::DigitalFilterd(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd());
|
||||
BOOST_REQUIRE(dfd.status() == fratio::FilterStatus::B_COEFF_MISSING);
|
||||
dfd = fratio::DigitalFilterd(Eigen::VectorXd(), Eigen::VectorXd());
|
||||
BOOST_REQUIRE(dfd.status() == fratio::FilterStatus::ALL_COEFF_MISSING);
|
||||
dfd = fratio::DigitalFilterd(Eigen::VectorXd::Constant(2, 0), Eigen::VectorXd::Constant(2, 0));
|
||||
BOOST_REQUIRE(dfd.status() == fratio::FilterStatus::BAD_A_COEFF);
|
||||
dfd = fratio::DigitalFilterd();
|
||||
BOOST_REQUIRE(dfd.status() == fratio::FilterStatus::NONE);
|
||||
dfd = fratio::DigitalFilterd(Eigen::VectorXd::Constant(2, 1), Eigen::VectorXd::Constant(2, 0));
|
||||
BOOST_REQUIRE(dfd.status() == fratio::FilterStatus::READY);
|
||||
}
|
||||
|
|
|
@ -5,39 +5,39 @@
|
|||
|
||||
template <typename T>
|
||||
struct System {
|
||||
std::vector<T> data = { 1, 2, 3, 4, 5, 6 };
|
||||
Eigen::VectorX<T> data = (Eigen::VectorX<T>(6) << 1, 2, 3, 4, 5, 6).finished();
|
||||
size_t windowSize = 4;
|
||||
std::vector<T> results = { 0.25, 0.75, 1.5, 2.5, 3.5, 4.5 };
|
||||
Eigen::VectorX<T> results = (Eigen::VectorX<T>(6) << 0.25, 0.75, 1.5, 2.5, 3.5, 4.5).finished();
|
||||
};
|
||||
|
||||
BOOST_FIXTURE_TEST_CASE(MOVING_AVERAGE_FLOAT, System<float>)
|
||||
{
|
||||
auto maf = fratio::MovingAveragef(windowSize);
|
||||
std::vector<float> filteredData;
|
||||
for (float d : data)
|
||||
filteredData.push_back(maf.stepFilter(d));
|
||||
for (Eigen::Index i = 0; i < data.size(); ++i)
|
||||
filteredData.push_back(maf.stepFilter(data(i)));
|
||||
|
||||
for (size_t i = 0; i < filteredData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results[i]), 1e-6f);
|
||||
|
||||
maf.resetFilter();
|
||||
filteredData = maf.filter(data);
|
||||
for (size_t i = 0; i < filteredData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results[i]), 1e-6f);
|
||||
Eigen::VectorXf fData = maf.filter(data);
|
||||
for (Eigen::Index i = 0; i < fData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(fData(i) - results(i)), 1e-6f);
|
||||
}
|
||||
|
||||
BOOST_FIXTURE_TEST_CASE(MOVING_AVERAGE_DOUBLE, System<double>)
|
||||
{
|
||||
auto maf = fratio::MovingAveraged(windowSize);
|
||||
std::vector<double> filteredData;
|
||||
for (double d : data)
|
||||
filteredData.push_back(maf.stepFilter(d));
|
||||
for (Eigen::Index i = 0; i < data.size(); ++i)
|
||||
filteredData.push_back(maf.stepFilter(data(i)));
|
||||
|
||||
for (size_t i = 0; i < filteredData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results[i]), 1e-14);
|
||||
|
||||
maf.resetFilter();
|
||||
filteredData = maf.filter(data);
|
||||
for (size_t i = 0; i < filteredData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(filteredData[i] - results[i]), 1e-14);
|
||||
Eigen::VectorXd fData = maf.filter(data);
|
||||
for (Eigen::Index i = 0; i < fData.size(); ++i)
|
||||
BOOST_CHECK_SMALL(std::abs(fData(i) - results(i)), 1e-14);
|
||||
}
|
||||
|
|
Ładowanie…
Reference in New Issue