Skip to content

Commit

Permalink
added BasicFilter<T>
Browse files Browse the repository at this point in the history 
implementing configurable
 * IIR or FIR
 * low-, band-, high-pass or bandstop filter
 * decimating and non-decimating version
based on basic user-configurable parameter.

The version implements both fundamental and error propagating (i.e. `UncertainValue<T>`) floating-point types.

Minor fixes:
 * added default execution policy from `std::execution::seq` to `std::execution::unseq` (should be faster ... eventually
 * added missing default [ErrorPropagating]Filter<> constructor in filter tool
 * fixed erroneous 'std::forward<...>' in ErrorPropagatingFilter

Signed-off-by: rstein <[email protected]>
  • Loading branch information
@RalphSteinhagen
RalphSteinhagen committed Nov 20, 2024
1 parent 6b1966a commit 3b9d262
Show file tree
Hide file tree
Showing 4 changed files with 357 additions and 33 deletions.
21 changes: 14 additions & 7 deletions algorithm/include/gnuradio-4.0/algorithm/filter/FilterTool.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,7 @@ namespace detail {
template<typename T, std::size_t bufferSize = std::dynamic_extent, typename TBaseType = meta::fundamental_base_value_type_t<T>>
struct Section;

template<typename T, std::size_t bufferSize, Form form = std::is_floating_point_v<T> ? Form::DF_II : Form::DF_I, auto execPolicy = std::execution::seq>
template<typename T, std::size_t bufferSize, Form form = std::is_floating_point_v<T> ? Form::DF_II : Form::DF_I, auto execPolicy = std::execution::unseq>
[[nodiscard]] inline constexpr T computeFilter(const T& input, Section<T, bufferSize>& section) noexcept {
const auto& a = section.a;
const auto& b = section.b;
Expand Down Expand Up @@ -158,7 +158,7 @@ template<typename T, std::size_t bufferSize, Form form = std::is_floating_point_
}
}

template<typename T, std::size_t bufferSize, Form form = std::is_floating_point_v<T> ? Form::DF_II : Form::DF_I, auto execPolicy = std::execution::seq>
template<typename T, std::size_t bufferSize, Form form = std::is_floating_point_v<T> ? Form::DF_II : Form::DF_I, auto execPolicy = std::execution::unseq>
inline constexpr std::vector<T> computeImpulseResponse(Section<T, bufferSize>& section, std::size_t length) {
std::vector<T> impulseResponse(length, T(0));
T input = T(1); // impulse response: first input is 1
Expand Down Expand Up @@ -218,7 +218,7 @@ struct Section : public FilterCoefficients<TBaseType> {
* Filter<double> myFilter(filterSections);
* double outputSample = myFilter.processOne(inputSample);
*/
template<typename T, std::size_t bufferSize = std::dynamic_extent, Form form = std::is_floating_point_v<meta::fundamental_base_value_type_t<T>> ? Form::DF_II : Form::DF_I, auto execPolicy = std::execution::seq>
template<typename T, std::size_t bufferSize = std::dynamic_extent, Form form = std::is_floating_point_v<meta::fundamental_base_value_type_t<T>> ? Form::DF_II : Form::DF_I, auto execPolicy = std::execution::unseq>
struct Filter;

template<typename T, std::size_t bufferSize, Form form, auto execPolicy>
Expand All @@ -227,7 +227,8 @@ struct Filter<T, bufferSize, form, execPolicy> {
using TBaseType = meta::fundamental_base_value_type_t<T>;
alignas(64UZ) std::vector<detail::Section<T, bufferSize>> _sectionsMeanValue;

public:
constexpr Filter() noexcept { _sectionsMeanValue.emplace_back(FilterCoefficients<TBaseType>{.b = {1}, .a = {1}}); }

template<typename... TFilterCoefficients>
explicit Filter(TFilterCoefficients&&... filterSections) noexcept {
std::vector<FilterCoefficients<TBaseType>> filterSections_{std::forward<TFilterCoefficients>(filterSections)...};
Expand Down Expand Up @@ -298,7 +299,8 @@ struct Filter<UncertainValue<T>, bufferSize, form, execPolicy> {
return totalUncertainty;
}

public:
constexpr Filter() noexcept { _sectionsMeanValue.emplace_back(FilterCoefficients<TBaseType>{.b = {1}, .a = {1}}); }

template<typename... TFilterCoefficients>
explicit Filter(TFilterCoefficients&&... filterSections) noexcept {
std::vector<FilterCoefficients<TBaseType>> filterSections_{std::forward<TFilterCoefficients>(filterSections)...};
Expand All @@ -322,15 +324,20 @@ struct Filter<UncertainValue<T>, bufferSize, form, execPolicy> {
}
};

template<typename T, std::size_t bufferSize = std::dynamic_extent, Form form = std::is_floating_point_v<T> ? Form::DF_II : Form::DF_I, auto execPolicy = std::execution::seq>
template<typename T, std::size_t bufferSize = std::dynamic_extent, Form form = std::is_floating_point_v<T> ? Form::DF_II : Form::DF_I, auto execPolicy = std::execution::unseq>
struct ErrorPropagatingFilter {
using TBaseType = meta::fundamental_base_value_type_t<T>;

Filter<T, bufferSize, form, execPolicy> filterMean;
Filter<TBaseType, bufferSize, form, execPolicy> filterSquared;

constexpr ErrorPropagatingFilter() noexcept {
filterMean = Filter<T, bufferSize, form, execPolicy>(FilterCoefficients<TBaseType>{.b = {1}, .a = {1}});
filterSquared = Filter<TBaseType, bufferSize, form, execPolicy>(FilterCoefficients<TBaseType>{.b = {1}, .a = {1}});
}

template<typename... TFilterCoefficients>
explicit ErrorPropagatingFilter(TFilterCoefficients&&... filterSections) : filterMean(std::forward<TFilterCoefficients>(filterSections)...), filterSquared(std::forward<TFilterCoefficients>(filterSections)...) {}
explicit ErrorPropagatingFilter(TFilterCoefficients&&... filterSections) : filterMean(filterSections...), filterSquared(filterSections...) {}

constexpr void reset(T defaultValue = T()) {
filterMean.reset(defaultValue);
Expand Down
181 changes: 166 additions & 15 deletions blocks/filter/include/gnuradio-4.0/filter/time_domain_filter.hpp
View file
Original file line number Diff line number Diff line change
@@ -1,10 +1,18 @@
#ifndef GNURADIO_TIME_DOMAIN_FILTER_HPP
#define GNURADIO_TIME_DOMAIN_FILTER_HPP
#include <algorithm>
#include <execution>
#include <functional>
#include <numeric>
#include <vector>

#include <gnuradio-4.0/Block.hpp>
#include <gnuradio-4.0/BlockRegistry.hpp>
#include <gnuradio-4.0/HistoryBuffer.hpp>
#include <gnuradio-4.0/algorithm/filter/FilterTool.hpp>
#include <gnuradio-4.0/meta/UncertainValue.hpp>

#include <magic_enum.hpp>

namespace gr::filter {

Expand All @@ -21,21 +29,21 @@ H(z) = b[0] + b[1]*z^-1 + b[2]*z^-2 + ... + b[N]*z^-N
)"">;
PortIn<T> in;
PortOut<T> out;
std::vector<T> b{}; // feedforward coefficients
std::vector<T> b{T{1}}; // feedforward coefficients

GR_MAKE_REFLECTABLE(fir_filter, in, out, b);

HistoryBuffer<T> inputHistory{32};

void settingsChanged(const property_map& /*old_settings*/, const property_map& new_settings) noexcept {
if (new_settings.contains("b") && b.size() >= inputHistory.capacity()) {
if (new_settings.contains("b") && b.size() > inputHistory.capacity()) {
inputHistory = HistoryBuffer<T>(std::bit_ceil(b.size()));
}
}

constexpr T processOne(T input) noexcept {
inputHistory.push_back(input);
return std::inner_product(b.cbegin(), b.cend(), inputHistory.cbegin(), static_cast<T>(0));
return std::transform_reduce(std::execution::unseq, b.cbegin(), b.cend(), inputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{});
}
};

Expand Down Expand Up @@ -78,27 +86,28 @@ a are the feedback coefficients
// y[n] = b[0] * x[n] + b[1] * x[n-1] + ... + b[N] * x[n-N]
// - a[1] * y[n-1] - a[2] * y[n-2] - ... - a[M] * y[n-M]
inputHistory.push_back(input);
const T output = std::inner_product(b.cbegin(), b.cend(), inputHistory.cbegin(), static_cast<T>(0)) // feed-forward path
- std::inner_product(a.cbegin() + 1, a.cend(), outputHistory.cbegin(), static_cast<T>(0)); // feedback path
const T feedforward = std::transform_reduce(std::execution::unseq, b.cbegin(), b.cend(), inputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{});
const T feedback = std::transform_reduce(std::execution::unseq, a.cbegin() + 1, a.cend(), outputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{});
const T output = feedforward - feedback;
outputHistory.push_back(output);
return output;
} else if constexpr (form == IIRForm::DF_II) {
// w[n] = x[n] - a[1] * w[n-1] - a[2] * w[n-2] - ... - a[M] * w[n-M]
// y[n] = b[0] * w[n] + b[1] * w[n-1] + ... + b[N] * w[n-N]
const T w = input - std::inner_product(a.cbegin() + 1, a.cend(), inputHistory.cbegin(), T{0});
const T w = input - std::transform_reduce(std::execution::unseq, a.cbegin() + 1, a.cend(), inputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{});
inputHistory.push_back(w);
return std::inner_product(b.cbegin(), b.cend(), inputHistory.cbegin(), T{0});

return std::transform_reduce(std::execution::unseq, b.cbegin(), b.cend(), inputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{});
} else if constexpr (form == IIRForm::DF_I_TRANSPOSED) {
// w_1[n] = x[n] - a[1] * w_2[n-1] - a[2] * w_2[n-2] - ... - a[M] * w_2[n-M]
// y[n] = b[0] * w_2[n] + b[1] * w_2[n-1] + ... + b[N] * w_2[n-N]
T v0 = input - std::inner_product(a.cbegin() + 1, a.cend(), outputHistory.cbegin(), static_cast<T>(0));
const T v0 = input - std::transform_reduce(std::execution::unseq, a.cbegin() + 1, a.cend(), outputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{});
outputHistory.push_back(v0);
return std::inner_product(b.cbegin(), b.cend(), outputHistory.cbegin(), static_cast<T>(0));

return std::transform_reduce(std::execution::unseq, b.cbegin(), b.cend(), outputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{});
} else if constexpr (form == IIRForm::DF_II_TRANSPOSED) {
// y[n] = b_0*f[n] + \sum_(k=1)^N(b_k*f[n−k] − a_k*y[n−k])
const T output = b[0] * input //
+ std::inner_product(b.cbegin() + 1, b.cend(), inputHistory.cbegin(), static_cast<T>(0)) //
- std::inner_product(a.cbegin() + 1, a.cend(), outputHistory.cbegin(), static_cast<T>(0));
// y[n] = b_0 * f[n] + Σ (b_k * f[n−k] − a_k * y[n−k]) for k = 1 to N
const T output = b[0] * input + std::transform_reduce(std::execution::unseq, b.cbegin() + 1, b.cend(), inputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{}) - std::transform_reduce(std::execution::unseq, a.cbegin() + 1, a.cend(), outputHistory.cbegin(), T{0}, std::plus<>{}, std::multiplies<>{});

inputHistory.push_back(input);
outputHistory.push_back(output);
Expand All @@ -107,9 +116,151 @@ a are the feedback coefficients
}
};

template<typename T, typename... Args>
requires(std::floating_point<T> or std::is_arithmetic_v<meta::fundamental_base_value_type_t<T>>)
struct BasicFilterProto : Block<BasicFilterProto<T, Args...>, Args...> {
using TParent = Block<BasicFilterProto<T, Args...>, Args...>;
using Description = Doc<R""(@brief Basic Digital Filter class supporting FIR and IIR filters
This block implements a digital filter which can be configured as either FIR or IIR,
with selectable filter type (low-pass, high-pass, band-pass, band-stop), and supports resampling.
)"">;
using ValueType = meta::fundamental_base_value_type_t<T>;

PortIn<T> in;
PortOut<T> out;

// private internal state variables
enum class FilterType { FIR, IIR };

FilterType _filter_type{FilterType::IIR};
filter::Type _filter_response{filter::Type::LOWPASS};
filter::iir::Design _iir_design_method{filter::iir::Design::BUTTERWORTH};
algorithm::window::Type _fir_design_method{algorithm::window::Type::Kaiser};

using FilterImpl = std::conditional_t<UncertainValueLike<T>, filter::ErrorPropagatingFilter<T>, filter::Filter<T>>;

FilterImpl _filter;

// Public settings
Annotated<std::string, "filter_type", Doc<"Filter type ('FIR' or 'IIR')">, Visible> filter_type = std::string(magic_enum::enum_name(_filter_type));
Annotated<std::string, "filter_response", Doc<"Filter response ('LOWPASS', 'HIGHPASS', 'BANDPASS', 'BANDSTOP')">, Visible> filter_response = std::string(magic_enum::enum_name(_filter_response));
Annotated<gr::Size_t, "filter_order", Doc<"Filter order">> filter_order{3};
Annotated<double, "f_low", Doc<"Low cutoff frequency in Hz">, Visible> f_low{0.1};
Annotated<double, "f_high", Doc<"High cutoff frequency in Hz (only for BANDPASS/BANDSTOP)">, Visible> f_high{0.2};
Annotated<double, "sample rate", Doc<"Sample rate in Hz">, Visible> sample_rate{1.0};
Annotated<gr::Size_t, "decimation factor", Doc<"1: none, i.e. preserving the relationship: N_out = N_in/decimate">> decimate{1U};
Annotated<std::string, "iir_design_method", Doc<"IIR Filter design method ('BUTTERWORTH', 'BESSEL', 'CHEBYSHEV1', 'CHEBYSHEV2')">> iir_design_method = std::string(magic_enum::enum_name(_iir_design_method));
Annotated<std::string, "fir_design_method", Doc<"FIR Filter design method ('None', 'Rectangular', 'Hamming', 'Hann', 'HannExp', 'Blackman', 'Nuttall', 'BlackmanHarris', 'BlackmanNuttall', 'FlatTop', 'Exponential', 'Kaiser')">> //
fir_design_method = std::string(magic_enum::enum_name(_fir_design_method));

GR_MAKE_REFLECTABLE(BasicFilterProto, in, out, filter_type, filter_response, filter_order, f_low, f_high, sample_rate, decimate, iir_design_method, fir_design_method);

void settingsChanged(const property_map& /*oldSettings*/, const property_map& /*newSettings*/) { designFilter(); }

void designFilter() {
using namespace gr::filter;

if constexpr (not TParent::ResamplingControl::kIsConst) {
this->input_chunk_size = decimate;
}

auto cast_enum = []<typename V>(const V&, const std::string& str) {
auto result = magic_enum::enum_cast<V>(str);
if (!result.has_value()) {
throw gr::exception(fmt::format("invalid value for {}: {}", magic_enum::enum_type_name<V>(), str));
}
return result.value();
};

_filter_type = cast_enum(_filter_type, filter_type);
_filter_response = cast_enum(_filter_response, filter_response);
_iir_design_method = cast_enum(_iir_design_method, iir_design_method);
_fir_design_method = cast_enum(_fir_design_method, fir_design_method);

// Set up filter parameters
FilterParameters params;
params.order = filter_order;
params.fLow = f_low;
params.fHigh = f_high;
params.fs = sample_rate;

if (_filter_type == FilterType::FIR) { // design FIR filter
_filter = FilterImpl(fir::designFilter<ValueType>(_filter_response, params, _fir_design_method));
} else if (_filter_type == FilterType::IIR) { // design IIR filter
_filter = FilterImpl(iir::designFilter<ValueType>(_filter_response, params, _iir_design_method));
}
}

[[nodiscard]] T processOne(T input) noexcept
requires(TParent::ResamplingControl::kIsConst)
{
return _filter.processOne(input);
}

[[nodiscard]] work::Status processBulk(std::span<const T> input, std::span<T> output) noexcept
requires(not TParent::ResamplingControl::kIsConst)
{
assert(output.size() >= input.size() / decimate);

std::size_t out_sample_idx = 0;
for (std::size_t i = 0; i < input.size(); ++i) {
T output_sample = _filter.processOne(input[i]);

if (i % decimate == 0) {
output[out_sample_idx++] = output_sample;
}
}
return work::Status::OK;
}
};

template<typename T>
using BasicFilter = BasicFilterProto<T>;

template<typename T>
using BasicDecimatingFilter = BasicFilterProto<T, Resampling<1UZ, 1UZ, false>>;

template<typename T>
struct Decimator : Block<Decimator<T>, Resampling<1UZ, 1UZ, false>> {
using TParent = Block<Decimator<T>, Resampling<1UZ, 1UZ, false>>;
using Description = Doc<R""(@brief Basic Decimator Block
This block implements a decimator for downsampling (dropping) input data by a
configurable factor. Filtering is not included in this implementation so expect
aliasing and sub-sampling related effects.
)"">;

PortIn<T> in;
PortOut<T> out;

Annotated<gr::Size_t, "decimation factor", Doc<"Factor by which to downsample/drop input data">, Visible> decim{1};

GR_MAKE_REFLECTABLE(Decimator, in, out, decim);

void settingsChanged(const property_map& /*oldSettings*/, const property_map& /*newSettings*/) { this->input_chunk_size = decim; }

[[nodiscard]] work::Status processBulk(std::span<const T> input, std::span<T> output) noexcept {
assert(output.size() >= input.size() / decim);

std::size_t out_sample_idx = 0;
for (std::size_t i = 0; i < input.size(); ++i) {
if (i % decim == 0) {
output[out_sample_idx++] = input[i];
}
}
return work::Status::OK;
}
};

} // namespace gr::filter

auto registerFirFilter = gr::registerBlock<gr::filter::fir_filter, double, float>(gr::globalBlockRegistry());
auto registerIirFilter = gr::registerBlock<gr::filter::iir_filter, gr::filter::IIRForm::DF_I, double, float>(gr::globalBlockRegistry()) | gr::registerBlock<gr::filter::iir_filter, gr::filter::IIRForm::DF_II, double, float>(gr::globalBlockRegistry()) | gr::registerBlock<gr::filter::iir_filter, gr::filter::IIRForm::DF_I_TRANSPOSED, double, float>(gr::globalBlockRegistry()) | gr::registerBlock<gr::filter::iir_filter, gr::filter::IIRForm::DF_II_TRANSPOSED, double, float>(gr::globalBlockRegistry());
inline static auto registerFilter = gr::registerBlock<gr::filter::fir_filter, double, float>(gr::globalBlockRegistry()) //
+ gr::registerBlock<gr::filter::iir_filter, gr::filter::IIRForm::DF_I, double, float>(gr::globalBlockRegistry()) //
+ gr::registerBlock<gr::filter::iir_filter, gr::filter::IIRForm::DF_II, double, float>(gr::globalBlockRegistry()) //
+ gr::registerBlock<gr::filter::iir_filter, gr::filter::IIRForm::DF_I_TRANSPOSED, double, float>(gr::globalBlockRegistry()) + gr::registerBlock<gr::filter::iir_filter, gr::filter::IIRForm::DF_II_TRANSPOSED, double, float>(gr::globalBlockRegistry()) //
+ gr::registerBlock<gr::filter::BasicFilter, double, float, gr::UncertainValue<float>, gr::UncertainValue<double>>(gr::globalBlockRegistry()) //
+ gr::registerBlock<gr::filter::BasicDecimatingFilter, double, float, gr::UncertainValue<float>, gr::UncertainValue<double>>(gr::globalBlockRegistry()) //
+ gr::registerBlock<gr::filter::Decimator, uint8_t, int8_t, uint16_t, int16_t, uint32_t, int32_t, uint64_t, int64_t, float, double, std::complex<float>, std::complex<double>, gr::UncertainValue<float>, gr::UncertainValue<double>>(gr::globalBlockRegistry());

#endif // GNURADIO_TIME_DOMAIN_FILTER_HPP
Loading

0 comments on commit 3b9d262

Please sign in to comment.