Skip to content

Commit

Permalink
added step- and pulse-response analysis functions
Browse files Browse the repository at this point in the history 
 * added 2nd-order resonance filter
 * added apply[Symmetric]Filter(..) math function
 * added enum to declare if math operation is performed `ProcessMode::InPlace` or on a `ProcessMode::Copy`
 * added enum to declare if estimates are added as meta-info (`ProcessMode::None`)
 * improved recurring `gr::dataset::detail::checkIndexRange(..)` helper function
 * fixed some missing `(ds, minIndex, maxIndex, ..)` method signatures.
 * added `gr::cast<T>(..)` to suppress float-double compiler warnings in templating code that cannot be resolved programmatically
 * added `filter::applyMedian(..)`
 * added `filter::applyRms(..)`
 * added `filter::applyPeakToPeak(..)`

Signed-off-by: Ralph J. Steinhagen <[email protected]>
  • Loading branch information
@RalphSteinhagen
RalphSteinhagen committed Jan 11, 2025
1 parent 0f7e5a1 commit 83bcce7
Show file tree
Hide file tree
Showing 10 changed files with 1,208 additions and 191 deletions.
515 changes: 363 additions & 152 deletions algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetEstimators.hpp
View file

Large diffs are not rendered by default.

41 changes: 33 additions & 8 deletions algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetHelper.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,11 +10,31 @@ inline static constexpr std::size_t Y = 1UZ; /// Y-axis index
inline static constexpr std::size_t Z = 2UZ; /// Z-axis index
} // namespace dim

enum class ProcessMode : std::uint8_t {
InPlace = 0U, /// in-place processing
Copy /// copy first, then process
};

enum class MetaInfo : std::uint8_t {
None = 0U, /// do nothing
Apply /// add meta-info to DataSet<T>::meta_information[] or ::timing_events[] where applicable
};

namespace detail {
inline void checkRangeIndex(std::size_t minIndex, std::size_t maxIndex, std::size_t dataCount, std::source_location location = std::source_location::current()) {
if (minIndex > maxIndex || maxIndex > dataCount) {
throw gr::exception(fmt::format("{} invalid range: [{},{}), dataCount={}", location, minIndex, maxIndex, dataCount));
}

template<typename T>
std::size_t checkIndexRange(const DataSet<T>& ds, std::size_t minIndex = 0UZ, std::size_t maxIndex = 0UZ, std::size_t signalIndex = 0UZ, std::source_location location = std::source_location::current()) {
const std::size_t maxDataSetIndex = ds.axisValues(dim::X).size();
if (maxIndex == max_size_t) { // renormalise default range
maxIndex = maxDataSetIndex;
}
if (minIndex > maxIndex || minIndex >= maxDataSetIndex || maxIndex > maxDataSetIndex || signalIndex >= ds.size()) {
throw gr::exception(fmt::format("DataSet<{}> ({}/{}: \"{}\") indices [{}, {}] out of range [0, {}]", //
gr::meta::type_name<T>(), signalIndex, ds.size(), signalIndex < ds.size() ? ds.signalName(signalIndex) : "???", //
minIndex, maxIndex, maxDataSetIndex),
location);
}
return maxIndex;
}

template<typename T, typename U>
Expand Down Expand Up @@ -56,11 +76,11 @@ template<typename T, typename TValue = gr::meta::fundamental_base_value_type_t<T
}

template<typename T>
T getDistance(const DataSet<T>& dataSet, std::size_t dimIndex, std::size_t indexMin = 0UZ, std::size_t indexMax = 0UZ, std::size_t signalIndex = 0UZ) {
if (indexMax == 0UZ) { // renormalise default range
T getDistance(const DataSet<T>& dataSet, std::size_t dimIndex, std::size_t indexMin = 0UZ, std::size_t indexMax = max_size_t, std::size_t signalIndex = 0UZ) {
if (indexMax == max_size_t) { // renormalise default range
indexMax = dataSet.axisValues(dim::X).size() - 1UZ;
}
detail::checkRangeIndex(indexMin, indexMax, dataSet.axisValues(dim::X).size());
indexMax = detail::checkIndexRange(dataSet, indexMin, indexMax, signalIndex);

if (dimIndex == dim::X || dimIndex == dim::Y) {
return getIndexValue(dataSet, dimIndex, indexMax, signalIndex) - getIndexValue(dataSet, dimIndex, indexMin, signalIndex);
Expand Down Expand Up @@ -131,7 +151,7 @@ std::vector<T> getSubArrayCopy(const DataSet<T>& ds, std::size_t indexMin, std::
if (indexMax <= indexMin) {
return {};
}
detail::checkRangeIndex(indexMin, indexMax, ds.axisValues(0).size());
indexMax = detail::checkIndexRange(ds, indexMin, indexMax, signalIndex);
std::span<const T> values = ds.signalValues(signalIndex);
return std::vector<T>(values.begin() + static_cast<std::ptrdiff_t>(indexMin), values.begin() + static_cast<std::ptrdiff_t>(indexMax));
}
Expand All @@ -152,6 +172,11 @@ template<typename T, typename TValue = gr::meta::fundamental_base_value_type_t<T
return T(20) * gr::math::log10(x); // 20 * log10(x)
}

template<typename T>
[[nodiscard]] constexpr T inverseDecibel(T x) noexcept {
return gr::math::pow(T(10), x / T(20)); // Inverse decibel => 10^(value / 20)
}

template<bool ThrowOnFailure = true, typename T>
[[maybe_unused]] constexpr bool verify(const DataSet<T>& dataSet, std::source_location location = std::source_location::current()) {
auto handleFailure = [&](const std::string& message) -> bool {
Expand Down
256 changes: 255 additions & 1 deletion algorithm/include/gnuradio-4.0/algorithm/dataset/DataSetMath.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#define DATASETMATH_HPP

#include <fmt/format.h>
#include <random>

#include <gnuradio-4.0/DataSet.hpp>
#include <gnuradio-4.0/Message.hpp>
Expand Down Expand Up @@ -43,7 +44,7 @@ template<typename T, typename TValue = gr::meta::fundamental_base_value_type_t<T
case MathOp::SQRT: return (y1 + y2) > TValue(0) ? gr::math::sqrt(y1 + y2) : std::numeric_limits<TValue>::quiet_NaN();
case MathOp::LOG10: return tenLog10((y1 + y2));
case MathOp::DB: return decibel((y1 + y2));
case MathOp::INV_DB: return inverseDecibel(y1);
case MathOp::INV_DB: return inverseDecibel<T>(y1);
case MathOp::IDENTITY:
default: return (y1 + y2);
}
Expand Down Expand Up @@ -125,6 +126,259 @@ template<typename T>

// WIP complete other convenience and missing math functions

template<typename T>
std::vector<T> computeDerivative(const DataSet<T>& ds, std::size_t signalIndex = 0UZ) {
auto signal = ds.signalValues(signalIndex);
if (signal.size() < 2) {
throw gr::exception("signal must contain at least two samples to compute derivative.");
}

std::vector<T> derivative;
derivative.reserve(signal.size() - 1UZ);
for (std::size_t i = 1UZ; i < signal.size(); ++i) {
derivative.push_back(signal[i] - signal[i - 1UZ]);
}
return derivative;
}

template<ProcessMode mode = ProcessMode::Copy, typename T, typename TValue = gr::meta::fundamental_base_value_type_t<T>>
DataSet<T> addNoise(const DataSet<T>& ds, TValue noiseLevel, std::size_t signalIndex = 0UZ, std::uint64_t seed = 0U) {
if (noiseLevel < TValue(0)) {
throw gr::exception(fmt::format("noiseLevel {} must be a positive number.", noiseLevel));
}

DataSet<T> noisy;
if constexpr (mode == ProcessMode::Copy) { // copy
noisy = ds;
} else { // or move (in-place)
noisy = std::move(ds);
}
const auto signal = ds.signalValues(signalIndex);
auto noisySignal = noisy.signalValues(signalIndex);
std::random_device rd;
std::mt19937_64 rng(seed == 0 ? rd() : seed);
using Distribution = std::conditional_t<std::is_integral_v<TValue>, std::uniform_int_distribution<TValue>, std::uniform_real_distribution<TValue>>;

Distribution dist(-noiseLevel, +noiseLevel);

for (std::size_t i = 0UZ; i < signal.size(); ++i) {
noisySignal[i] += dist(rng);
}
return noisy;
}

namespace filter {

template<ProcessMode mode = ProcessMode::Copy, typename T, typename TValue = gr::meta::fundamental_base_value_type_t<T>>
DataSet<T> applyMovingAverage(const DataSet<T>& ds, std::size_t windowSize, std::size_t signalIndex = 0UZ) {
if (windowSize == 0 || !(windowSize & 1)) {
throw gr::exception("windowSize must be a positive odd number.");
}

DataSet<T> smoothed;
if constexpr (mode == ProcessMode::Copy) { // copy
smoothed = ds;
} else { // or move (in-place)
smoothed = std::move(ds);
}
const auto signal = ds.signalValues(signalIndex);
auto smoothedSignal = smoothed.signalValues(signalIndex);

const std::size_t halfWindow = windowSize / 2UZ;
for (std::size_t i = 0UZ; i < signal.size(); ++i) {
std::size_t start = (i >= halfWindow) ? i - halfWindow : 0UZ;
std::size_t end = std::min(i + halfWindow + 1UZ, signal.size());

T sum = std::accumulate(signal.begin() + static_cast<std::ptrdiff_t>(start), signal.begin() + static_cast<std::ptrdiff_t>(end), T(0));
smoothedSignal[i] = sum / TValue(end - start);
}
return smoothed;
}

template<ProcessMode mode = ProcessMode::Copy, typename T, typename TValue = gr::meta::fundamental_base_value_type_t<T>>
constexpr DataSet<T> applyMedian(const DataSet<T>& ds, std::size_t windowSize, std::size_t signalIndex = 0UZ, std::source_location location = std::source_location::current()) {
if (windowSize == 0) {
throw gr::exception(fmt::format("windowSize: {} must be a positive number.", windowSize), location);
}

DataSet<T> filtered;
if constexpr (mode == ProcessMode::Copy) {
filtered = ds;
} else {
filtered = std::move(ds);
}

const std::vector<T> signal{filtered.signalValues(signalIndex).begin(), filtered.signalValues(signalIndex).end()};
auto filteredSignal = filtered.signalValues(signalIndex);
const std::size_t N = filteredSignal.size();

std::vector<T> medianWindow(windowSize); // temporary mutable copy for in-place partitioning
const std::size_t halfWindow = windowSize / 2UZ;
for (std::size_t i = 0UZ; i < N; ++i) {
const auto start = static_cast<std::ptrdiff_t>(i > halfWindow ? i - halfWindow : 0UZ);
const auto end = static_cast<std::ptrdiff_t>(std::min(i + halfWindow + 1UZ, N));
const std::size_t size = static_cast<std::size_t>(end - start);

std::copy(signal.begin() + start, signal.begin() + end, medianWindow.begin());

auto medianWindowView = std::span(medianWindow.data(), size);
auto midIter = medianWindowView.begin() + (size / 2UZ);
std::ranges::nth_element(medianWindowView, midIter);

if ((size & 1UZ) == 0UZ) { // even-sized window -> take average around mid-point
auto midPrev = std::ranges::max_element(medianWindowView.begin(), midIter);
filteredSignal[i] = T(0.5) * (*midPrev + *midIter);
} else { // odd-sized window -> use exact mid-point
filteredSignal[i] = *midIter;
}
}

return filtered;
}

template<ProcessMode mode = ProcessMode::Copy, typename T, typename TValue = gr::meta::fundamental_base_value_type_t<T>>
constexpr DataSet<T> applyRms(const DataSet<T>& ds, std::size_t windowSize, std::size_t signalIndex = 0UZ, std::source_location location = std::source_location::current()) {
if (windowSize == 0UZ) {
throw gr::exception(fmt::format("windowSize: {} must be a positive number.", windowSize), location);
}

DataSet<T> filtered;
if constexpr (mode == ProcessMode::Copy) {
filtered = ds;
} else {
filtered = std::move(ds);
}

const std::vector<T> signal{filtered.signalValues(signalIndex).begin(), filtered.signalValues(signalIndex).end()};
auto filteredSignal = filtered.signalValues(signalIndex);
const std::size_t N = filteredSignal.size();

for (std::size_t i = 0; i < N; ++i) {
std::size_t start = (i > (windowSize / 2UZ)) ? i - (windowSize / 2UZ) : 0UZ;
std::size_t end = std::min(i + (windowSize / 2UZ) + 1UZ, N);
std::size_t size = (end - start);

T sum = T(0);
T sum2 = T(0);
for (std::size_t j = start; j < end; ++j) {
sum += signal[j];
sum2 += signal[j] * signal[j];
}
if (size > 1UZ) {
T mean = sum / gr::cast<TValue>(size);
filteredSignal[i] = gr::math::sqrt(gr::math::abs(sum2 / gr::cast<TValue>(size) - mean * mean));
} else {
filteredSignal[i] = gr::cast<TValue>(0);
}
}

return filtered;
}

template<ProcessMode mode = ProcessMode::Copy, typename T, typename TValue = gr::meta::fundamental_base_value_type_t<T>>
constexpr DataSet<T> applyPeakToPeak(const DataSet<T>& ds, std::size_t windowSize, std::size_t signalIndex = 0UZ, std::source_location location = std::source_location::current()) {
if (windowSize == 0UZ) {
throw gr::exception(fmt::format("windowSize: {} must be a positive number.", windowSize), location);
}

DataSet<T> filtered;
if constexpr (mode == ProcessMode::Copy) {
filtered = ds;
} else {
filtered = std::move(ds);
}

const std::vector<T> signal{filtered.signalValues(signalIndex).begin(), filtered.signalValues(signalIndex).end()};
auto filteredSignal = filtered.signalValues(signalIndex);
const std::size_t N = filteredSignal.size();

const std::size_t halfWindow = windowSize / 2UZ;

for (std::size_t i = 0; i < N; ++i) {
std::size_t start = (i > halfWindow) ? i - halfWindow : 0UZ;
std::size_t end = std::min(i + halfWindow + 1UZ, N);

auto minVal = std::numeric_limits<T>::max();
auto maxVal = std::numeric_limits<T>::lowest();
for (std::size_t j = start; j < end; ++j) {
if (signal[j] < minVal) {
minVal = signal[j];
}
if (signal[j] > maxVal) {
maxVal = signal[j];
}
}
filteredSignal[i] = maxVal - minVal;
}

return filtered;
}

template<ProcessMode mode = ProcessMode::Copy, bool symmetric = false, DataSetLike D, typename T = typename std::remove_cvref_t<D>::value_type, typename U>
DataSet<T> applyFilter(D&& dataSet, const gr::filter::FilterCoefficients<U>& coeffs, std::size_t signalIndex = max_size_t) {
using TValue = gr::meta::fundamental_base_value_type_t<T>;
constexpr bool isConstDataSet = std::is_const_v<std::remove_reference_t<D>>;

static_assert(!(isConstDataSet && mode == ProcessMode::InPlace), "cannot perform in-place computation on const DataSet<T>");

DataSet<T> smoothed;
if constexpr (mode == ProcessMode::Copy) { // copy
smoothed = dataSet;
} else { // or move (in-place)
smoothed = std::move(dataSet);
}

auto forwardPass = [&](std::span<T> signal) -> void {
auto filter = gr::filter::ErrorPropagatingFilter<T>(coeffs);
if constexpr (UncertainValueLike<TValue>) {
std::ranges::transform(signal, signal.begin(), [&](auto& val) { return filter.processOne(val); });
} else {
std::ranges::transform(signal, signal.begin(), [&](auto& val) { return gr::value(filter.processOne(val)); });
}
};

auto backwardPass = [&](std::span<T> signal) -> void {
auto filter = gr::filter::ErrorPropagatingFilter<T>(coeffs);
if constexpr (UncertainValueLike<TValue>) {
std::ranges::transform(signal | std::views::reverse, signal.rbegin(), [&](auto& val) { return filter.processOne(val); });
} else {
std::ranges::transform(signal | std::views::reverse, signal.rbegin(), [&](auto& val) { return gr::value(filter.processOne(val)); });
}
};

auto processSignal = [&](std::size_t sigIndex) {
auto signal = smoothed.signalValues(sigIndex);
if constexpr (!symmetric) {
forwardPass(signal); // forward pass only
} else {
std::vector<T> copy(signal.begin(), signal.end());
forwardPass(signal); // forward pass
backwardPass(copy); // reverse pass
// combine both passes
for (std::size_t i = 0UZ; i < signal.size(); i++) {
signal[i] = T{0.5} * (signal[i] + copy[i]);
}
}
};

if (signalIndex != max_size_t) {
processSignal(signalIndex);
} else {
for (std::size_t dsIndex = 0UZ; dsIndex < smoothed.size(); dsIndex++) {
processSignal(dsIndex);
}
}

return smoothed;
}

template<ProcessMode Mode = ProcessMode::Copy, typename T, typename... TFilterCoefficients>
DataSet<T> applySymmetricFilter(const DataSet<T>& ds, TFilterCoefficients&&... coeffs) {
return applyFilter<Mode, true>(ds, 0UZ, std::forward<TFilterCoefficients>(coeffs)...);
}

} // namespace filter

} // namespace gr::dataset

#endif // DATASETMATH_HPP
Loading

0 comments on commit 83bcce7

Please sign in to comment.