Butterworth BandPass ok.

topic/diffentiators
Vincent Samy 2018-12-18 16:16:47 +09:00
rodzic a26692666f
commit fcd4859401
11 zmienionych plików z 183 dodań i 66 usunięć

Wyświetl plik

@ -8,6 +8,8 @@
namespace fratio {
// https://www.dsprelated.com/showarticle/1119.php
// https://www.dsprelated.com/showarticle/1135.php
// https://www.dsprelated.com/showarticle/1128.php
template <typename T>
class Butterworth : public DigitalFilter<T> {
public:
@ -16,7 +18,9 @@ public:
public:
enum class Type {
LowPass,
HighPass
HighPass,
BandPass,
BandReject
};
// public:
@ -24,22 +28,27 @@ public:
public:
Butterworth(Type type = Type::LowPass);
Butterworth(int order, T fc, T fs, Type type = Type::LowPass);
Butterworth(int order, T bw, T fs, T fCenter, Type type = Type::BandPass);
void setFilterParameters(int order, T fc, T fs);
void setFilterParameters(int order, T fc, T fs, T fCenter = T(0));
private:
void initialize(int order, T fc, T fs);
void computeDigitalRep();
std::complex<T> generateAnalogPole(T fpw, int k);
vectX_t<std::complex<T>> generateAnalogZeros();
void initialize(int order, T fc, T fs, T fCenter = T(0)); // fc = bw for bandPass filter
void computeDigitalRep(T fc);
void computeBandDigitalRep(T bw, T fCenter);
std::complex<T> generateAnalogPole(int k, T fpw);
std::pair<std::complex<T>, std::complex<T>> generateBandAnalogPole(int k, T fpw1, T fpw2);
vectXc_t<T> generateAnalogZeros();
void scaleAmplitude(Eigen::Ref<vectX_t<T>> aCoeff, Eigen::Ref<vectX_t<T>> bCoeff);
private:
Type m_type;
int m_order;
T m_fc;
T m_bw;
T m_fs;
vectX_t<std::complex<T>> m_poles;
T m_fCenter;
vectXc_t<T> m_poles;
};
} // namespace fratio

Wyświetl plik

@ -17,13 +17,20 @@ Butterworth<T>::Butterworth(int order, T fc, T fs, Type type)
}
template <typename T>
void Butterworth<T>::setFilterParameters(int order, T fc, T fs)
Butterworth<T>::Butterworth(int order, T fc, T fs, T fCenter, Type type)
: m_type(type)
{
initialize(order, fc, fs);
initialize(order, fc, fs, fCenter);
}
template <typename T>
void Butterworth<T>::initialize(int order, T fc, T fs)
void Butterworth<T>::setFilterParameters(int order, T fc, T fs, T fCenter)
{
initialize(order, fc, fs, fCenter);
}
template <typename T>
void Butterworth<T>::initialize(int order, T fc, T fs, T fCenter)
{
if (order <= 0) {
m_status = FilterStatus::BAD_ORDER_SIZE;
@ -35,35 +42,44 @@ void Butterworth<T>::initialize(int order, T fc, T fs)
return;
}
if ((m_type == Type::BandPass || m_type == Type::BandReject) && fCenter - fc / 2. <= 0) {
m_status = FilterStatus::BAD_BAND_FREQUENCY;
return;
}
if (m_fc > m_fs / 2.) {
m_status = FilterStatus::BAD_CUTOFF_FREQUENCY;
return;
}
m_order = order;
m_fc = fc;
m_fs = fs;
computeDigitalRep();
if (m_type == Type::LowPass || m_type == Type::HighPass)
computeDigitalRep(fc);
else
computeBandDigitalRep(fc, fCenter); // For band-like filters
resetFilter();
}
template <typename T>
void Butterworth<T>::computeDigitalRep()
void Butterworth<T>::computeDigitalRep(T fc)
{
m_fc = fc;
// Continuous pre-warped frequency
T fpw = (m_fs / PI) * std::tan(PI * m_fc / m_fs);
// Compute poles
vectXc_t<T> poles(m_order);
std::complex<T> analogPole;
vectX_t<std::complex<T>> poles(m_order);
for (int k = 1; k <= m_order; ++k) {
analogPole = generateAnalogPole(fpw, k);
BilinearTransform<std::complex<T>>::SToZ(m_fs, analogPole, poles(k - 1));
for (int k = 0; k < m_order; ++k) {
analogPole = generateAnalogPole(k + 1, fpw);
BilinearTransform<std::complex<T>>::SToZ(m_fs, analogPole, poles(k));
}
vectX_t<std::complex<T>> zeros = generateAnalogZeros();
vectX_t<std::complex<T>> a = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(poles);
vectX_t<std::complex<T>> b = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(zeros);
vectXc_t<T> zeros = generateAnalogZeros();
vectXc_t<T> a = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(poles);
vectXc_t<T> b = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(zeros);
vectX_t<T> aCoeff(m_order + 1);
vectX_t<T> bCoeff(m_order + 1);
for (int i = 0; i < m_order + 1; ++i) {
@ -76,10 +92,48 @@ void Butterworth<T>::computeDigitalRep()
}
template <typename T>
std::complex<T> Butterworth<T>::generateAnalogPole(T fpw, int k)
void Butterworth<T>::computeBandDigitalRep(T bw, T fCenter)
{
T scaleFactor = 2 * PI * fpw;
m_bw = bw;
m_fCenter = fCenter;
T f1 = m_fCenter - m_bw / T(2);
T f2 = m_fCenter + m_bw / T(2);
T fpw1 = (m_fs / PI) * std::tan(PI * f1 / m_fs);
T fpw2 = (m_fs / PI) * std::tan(PI * f2 / m_fs);
vectXc_t<T> poles(2 * m_order);
std::pair<std::complex<T>, std::complex<T>> analogPoles;
for (int k = 0; k < m_order; ++k) {
analogPoles = generateBandAnalogPole(k + 1, fpw1, fpw2);
BilinearTransform<std::complex<T>>::SToZ(m_fs, analogPoles.first, poles(2 * k));
BilinearTransform<std::complex<T>>::SToZ(m_fs, analogPoles.second, poles(2 * k + 1));
}
vectXc_t<T> zeros = generateAnalogZeros();
vectXc_t<T> a = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(poles);
vectXc_t<T> b = VietaAlgo<std::complex<T>>::polyCoeffFromRoot(zeros);
vectX_t<T> aCoeff(2 * m_order + 1);
vectX_t<T> bCoeff(2 * m_order + 1);
for (int i = 0; i < 2 * m_order + 1; ++i) {
aCoeff(i) = a(i).real();
bCoeff(i) = b(i).real();
}
std::complex<T> s = std::exp(std::complex<T>(T(0), T(2) * PI * std::sqrt(f1 * f2) / m_fs));
std::complex<T> num(b(0));
std::complex<T> denum(a(0));
for (int i = 1; i < 2 * m_order + 1; ++i) {
num = num * s + b(i);
denum = denum * s + a(i);
}
std::complex<T> hf0 = num / denum;
bCoeff *= T(1) / std::abs(hf0); // bCoeff *= 1 / abs(num / denum)
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);
};
@ -87,54 +141,76 @@ std::complex<T> Butterworth<T>::generateAnalogPole(T fpw, int k)
std::complex<T> analogPole(-std::sin(thetaK(k)), std::cos(thetaK(k)));
switch (m_type) {
case Type::HighPass:
return scaleFactor / analogPole;
return T(2) * PI * fpw1 / analogPole;
case Type::LowPass:
default:
return scaleFactor * analogPole;
return T(2) * PI * fpw1 * analogPole;
}
}
template <typename T>
vectX_t<std::complex<T>> Butterworth<T>::generateAnalogZeros()
std::pair<std::complex<T>, std::complex<T>> Butterworth<T>::generateBandAnalogPole(int k, T fpw1, T fpw2)
{
auto thetaK = [pi = PI, order = m_order](int k) -> T {
return (2 * k - 1) * pi / (2 * order);
};
std::complex<T> analogPole(-std::sin(thetaK(k)), std::cos(thetaK(k)));
std::pair<std::complex<T>, std::complex<T>> poles;
switch (m_type) {
case Type::BandReject:
return poles;
case Type::BandPass:
default: {
std::complex<T> fpw0 = std::sqrt(fpw1 * fpw2);
std::complex<T> s = T(0.5) * (fpw2 - fpw1) * analogPole / fpw0;
poles.first = T(2) * PI * fpw0 * (s + std::complex<T>(T(0), T(1)) * std::sqrt(std::complex<T>(T(1), T(0)) - s * s));
poles.second = T(2) * PI * fpw0 * (s - std::complex<T>(T(0), T(1)) * std::sqrt(std::complex<T>(T(1), T(0)) - s * s));
return poles;
}
}
}
template <typename T>
vectXc_t<T> Butterworth<T>::generateAnalogZeros()
{
switch (m_type) {
case Type::HighPass:
return vectX_t<std::complex<T>>::Constant(m_order, std::complex<T>(1));
return vectXc_t<T>::Constant(m_order, std::complex<T>(1));
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::LowPass:
default:
return vectX_t<std::complex<T>>::Constant(m_order, std::complex<T>(-1));
return vectXc_t<T>::Constant(m_order, std::complex<T>(-1));
}
}
template <typename T>
void Butterworth<T>::scaleAmplitude(Eigen::Ref<vectX_t<T>> aCoeff, Eigen::Ref<vectX_t<T>> bCoeff)
{
T scale = 0;
T sumB = 0;
T num = 0;
T denum = 0;
switch (m_type) {
case Type::HighPass:
for (int i = 0; i < m_order + 1; ++i) {
if (i % 2 == 0) {
scale += aCoeff(i);
sumB += bCoeff(i);
num += aCoeff(i);
denum += bCoeff(i);
} else {
scale -= aCoeff(i);
sumB -= bCoeff(i);
num -= aCoeff(i);
denum -= bCoeff(i);
}
}
break;
case Type::LowPass:
default:
scale = aCoeff.sum();
sumB = bCoeff.sum();
num = aCoeff.sum();
denum = bCoeff.sum();
break;
}
bCoeff *= scale / sumB;
bCoeff *= num / denum;
}
} // namespace fratio

Wyświetl plik

@ -22,7 +22,7 @@ public:
void resetFilter();
template <typename T2>
bool setCoeffs(T2&& aCoeff, T2&& bCoeff);
void setCoeffs(T2&& aCoeff, T2&& bCoeff);
void getCoeffs(vectX_t<T>& aCoeff, vectX_t<T>& bCoeff) const;
FilterStatus status() const noexcept { return m_status; }

Wyświetl plik

@ -23,6 +23,8 @@ std::string GenericFilter<T>::filterStatus(FilterStatus status)
return "Filter has a received a frequency that is negative or equal to zero";
case FilterStatus::BAD_CUTOFF_FREQUENCY:
return "Filter has a received a bad cut-off frequency. It must be inferior to the sampling frequency";
case FilterStatus::BAD_BAND_FREQUENCY:
return "You try to initialize the filter with a bad combination of the frequency and bandwith, you must have fCenter > bw/2";
default:
return "I forgot to implement this error documentation";
}
@ -79,18 +81,17 @@ void GenericFilter<T>::resetFilter()
template <typename T>
template <typename T2>
bool GenericFilter<T>::setCoeffs(T2&& aCoeff, T2&& bCoeff)
void GenericFilter<T>::setCoeffs(T2&& aCoeff, T2&& bCoeff)
{
static_assert(std::is_same_v<T2, vectX_t<T>>, "The coefficents should be of type vectX_t<T>");
static_assert(std::is_convertible_v<T2, vectX_t<T>>, "The coefficients types should be convertible to vectX_t<T>");
if (!checkCoeffs(aCoeff, bCoeff))
return false;
return;
m_aCoeff = aCoeff;
m_bCoeff = bCoeff;
resetFilter();
normalizeCoeffs();
return true;
}
template <typename T>

Wyświetl plik

@ -10,8 +10,8 @@ class MovingAverage : public DigitalFilter<T> {
public:
MovingAverage() = default;
MovingAverage(int windowSize)
: DigitalFilter<T>(vectX_t<T>::Constant(1, T(1)), vectX_t<T>::Constant(windowSize, T(1) / windowSize))
{
setWindowSize(windowSize);
}
void setWindowSize(int windowSize)

Wyświetl plik

@ -6,6 +6,8 @@ namespace fratio {
template <typename T>
using vectX_t = Eigen::Matrix<T, Eigen::Dynamic, 1>;
template <typename T>
using vectXc_t = vectX_t<std::complex<T>>;
enum class FilterStatus {
// Generic filter
@ -19,7 +21,8 @@ enum class FilterStatus {
// Butterworth filter
BAD_FREQUENCY_VALUE,
BAD_CUTOFF_FREQUENCY
BAD_CUTOFF_FREQUENCY,
BAD_BAND_FREQUENCY
};
} // namespace fratio

Wyświetl plik

@ -1,5 +1,7 @@
#define BOOST_TEST_MODULE ButterworthFilterTests
// Note: In term of precision, LP > HP > BP ~= BR
#include "fratio"
#include "test_functions.h"
#include "warning_macro.h"
@ -13,6 +15,8 @@ struct System {
int order = 5;
T fc = 10;
T fs = 100;
T bw = 10;
T fCenter = 10;
// LP
fratio::vectX_t<T> lpACoeffRes = (fratio::vectX_t<T>(6) << 1, -2.975422109745684, 3.806018119320413, -2.545252868330468, 0.881130075437837, -0.125430622155356).finished();
fratio::vectX_t<T> lpBCoeffRes = (fratio::vectX_t<T>(6) << 0.001282581078961, 0.006412905394803, 0.012825810789607, 0.012825810789607, 0.006412905394803, 0.001282581078961).finished();
@ -21,6 +25,10 @@ struct System {
fratio::vectX_t<T> hpACoeffRes = (fratio::vectX_t<T>(6) << 1, -2.975422109745683, 3.806018119320411, -2.545252868330467, 0.8811300754378368, -0.1254306221553557).finished();
fratio::vectX_t<T> hpBCoeffRes = (fratio::vectX_t<T>(6) << 0.3541641810934298, -1.770820905467149, 3.541641810934299, -3.541641810934299, 1.770820905467149, -0.3541641810934298).finished();
fratio::vectX_t<T> hpResults = (fratio::vectX_t<T>(8) << 0.3541641810934298, -0.008704608374924483, -0.3113626313910076, -0.3460321436983160, -0.1787600153274098, 0.04471440201428267, 0.2059279258827846, 0.2533941579793959).finished();
// BP
fratio::vectX_t<T> bpACoeffRes = (fratio::vectX_t<T>(11) << 1, -6.784299264603903, 21.577693329895588, -42.338550072279737, 56.729081385507655, -54.208087151300411, 37.399203252161037, -18.397491390111661, 6.180883710485754, -1.283022311577260, 0.125430622155356).finished();
fratio::vectX_t<T> bpBCoeffRes = (fratio::vectX_t<T>(11) << 0.001282581078963, 0, -0.006412905394817, 0, 0.012825810789633, 0, -0.012825810789633, 0, 0.006412905394817, 0, -0.001282581078963).finished();
fratio::vectX_t<T> bpResults = (fratio::vectX_t<T>(8) << 0.001282581078963, 0.011266576028733, 0.046195520115810, 0.116904647483408, 0.200574194600111, 0.232153315136604, 0.141350142008155, -0.086403129422609).finished();
};
DISABLE_CONVERSION_WARNING_END
@ -29,30 +37,46 @@ BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_LP_FILTER_FLOAT, System<float>)
{
auto bf = fratio::Butterworthf(order, fc, fs);
BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
test_coeffs(lpACoeffRes, lpBCoeffRes, bf);
test_results(lpResults, data, bf);
test_coeffs(lpACoeffRes, lpBCoeffRes, bf, std::numeric_limits<float>::epsilon() * 10);
test_results(lpResults, data, bf, std::numeric_limits<float>::epsilon() * 100);
}
BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_LP_FILTER_DOUBLE, System<double>)
{
auto bf = fratio::Butterworthd(order, fc, fs);
BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
test_coeffs(lpACoeffRes, lpBCoeffRes, bf);
test_results(lpResults, data, bf);
test_coeffs(lpACoeffRes, lpBCoeffRes, bf, std::numeric_limits<double>::epsilon() * 10);
test_results(lpResults, data, bf, std::numeric_limits<double>::epsilon() * 100);
}
BOOST_FIXTURE_TEST_CASE(BUTTERWORTH_HP_FILTER_FLOAT, System<float>)
{
auto bf = fratio::Butterworthf(order, fc, fs, fratio::Butterworthf::Type::HighPass);
BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
test_coeffs(hpACoeffRes, hpBCoeffRes, bf);
test_results(hpResults, data, bf);
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 = fratio::Butterworthd(order, fc, fs, fratio::Butterworthd::Type::HighPass);
BOOST_REQUIRE_EQUAL(bf.aOrder(), bf.bOrder());
test_coeffs(hpACoeffRes, hpBCoeffRes, bf);
test_results(hpResults, data, bf);
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 = fratio::Butterworthf(order, bw, fs, fCenter);
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 = fratio::Butterworthd(order, bw, fs, fCenter);
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);
}

Wyświetl plik

@ -20,13 +20,13 @@ DISABLE_CONVERSION_WARNING_END
BOOST_FIXTURE_TEST_CASE(DIGITAL_FILTER_FLOAT, System<float>)
{
auto df = fratio::DigitalFilterf(aCoeff, bCoeff);
test_coeffs(aCoeff, bCoeff, df);
test_results(results, data, df);
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 = fratio::DigitalFilterd(aCoeff, bCoeff);
test_coeffs(aCoeff, bCoeff, df);
test_results(results, data, df);
test_coeffs(aCoeff, bCoeff, df, std::numeric_limits<double>::epsilon() * 10);
test_results(results, data, df, std::numeric_limits<double>::epsilon() * 10);
}

Wyświetl plik

@ -18,4 +18,8 @@ BOOST_AUTO_TEST_CASE(FILTER_FAILURES)
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);
auto mad = fratio::MovingAveraged(0);
BOOST_REQUIRE(mad.status() == fratio::FilterStatus::BAD_ORDER_SIZE);
auto bfd = fratio::Butterworthd(0, 10, 100);
BOOST_REQUIRE(bfd.status() == fratio::FilterStatus::BAD_ORDER_SIZE);
}

Wyświetl plik

@ -14,11 +14,11 @@ struct System {
BOOST_FIXTURE_TEST_CASE(MOVING_AVERAGE_FLOAT, System<float>)
{
auto maf = fratio::MovingAveragef(windowSize);
test_results(results, data, maf);
test_results(results, data, maf, std::numeric_limits<float>::epsilon() * 10);
}
BOOST_FIXTURE_TEST_CASE(MOVING_AVERAGE_DOUBLE, System<double>)
{
auto maf = fratio::MovingAveraged(windowSize);
test_results(results, data, maf);
test_results(results, data, maf, std::numeric_limits<double>::epsilon() * 10);
}

Wyświetl plik

@ -5,20 +5,20 @@
#include <limits>
template <typename T>
void test_coeffs(const fratio::vectX_t<T>& aCoeff, const fratio::vectX_t<T>& bCoeff, const fratio::GenericFilter<T>& filter)
void test_coeffs(const fratio::vectX_t<T>& aCoeff, const fratio::vectX_t<T>& bCoeff, const fratio::GenericFilter<T>& filter, T prec)
{
BOOST_REQUIRE_EQUAL(aCoeff.size(), filter.aOrder());
BOOST_REQUIRE_EQUAL(bCoeff.size(), filter.bOrder());
fratio::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)), std::numeric_limits<T>::epsilon() * 10);
BOOST_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)), std::numeric_limits<T>::epsilon() * 10);
BOOST_REQUIRE_SMALL(std::abs(bCoeff(i) - fbCoeff(i)), prec);
}
template <typename T>
void test_results(const fratio::vectX_t<T>& results, const fratio::vectX_t<T>& data, fratio::GenericFilter<T>& filter)
void test_results(const fratio::vectX_t<T>& results, const fratio::vectX_t<T>& data, fratio::GenericFilter<T>& filter, T prec)
{
fratio::vectX_t<T> filteredData(results.size());
@ -26,10 +26,10 @@ void test_results(const fratio::vectX_t<T>& results, const fratio::vectX_t<T>& d
filteredData(i) = filter.stepFilter(data(i));
for (Eigen::Index i = 0; i < filteredData.size(); ++i)
BOOST_REQUIRE_SMALL(std::abs(filteredData(i) - results(i)), std::numeric_limits<T>::epsilon() * 1000);
BOOST_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)), std::numeric_limits<T>::epsilon() * 1000);
BOOST_REQUIRE_SMALL(std::abs(filteredData(i) - results(i)), prec);
}