Skip to content

Commit

Permalink
Merge pull request #61 from petiaccja/interpolation-cleanup
Browse files Browse the repository at this point in the history
Interpolation cleanup
  • Loading branch information
petiaccja committed May 21, 2022
2 parents b3af8e2 + 9e5b7fb commit 2bcf452
Show file tree
Hide file tree
Showing 9 changed files with 891 additions and 251 deletions.
2 changes: 1 addition & 1 deletion include/dspbb/Filtering/FIR.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -155,7 +155,7 @@ namespace impl {

} // namespace impl

template <class SignalR, template <typename, typename...> class Desc, class... Params>
template <class SignalR, template <typename, typename...> class Desc, class... Params, std::enable_if_t<is_mutable_signal_v<SignalR>, int> = 0>
auto FirFilter(SignalR&& out, const Desc<impl::FirMethodLeastSquares, Params...>& desc)
-> decltype(void(impl::TranslateLeastSquares(desc))) {
const auto [response, weight] = impl::TranslateLeastSquares(desc);
Expand Down
332 changes: 242 additions & 90 deletions include/dspbb/Filtering/Interpolation.hpp

Large diffs are not rendered by default.

95 changes: 78 additions & 17 deletions include/dspbb/Filtering/Polyphase.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,41 +10,92 @@ namespace dspbb {


template <class T, eSignalDomain Domain>
struct PolyphaseDecomposition {
BasicSignalView<T, Domain> filterBank;
size_t numFilters;
class PolyphaseView {
public:
PolyphaseView() noexcept = default;
PolyphaseView(BasicSignalView<T, Domain> data, size_t numFilters) noexcept
: m_bufferView(data), m_filterCount(numFilters) {}

BasicSignalView<T, Domain> operator[](size_t index) {
auto loc = SubSignalLocation(index);
return filterBank.SubSignal(loc.first, loc.second);
return m_bufferView.SubSignal(loc.first, loc.second);
}
BasicSignalView<const T, Domain> operator[](size_t index) const {
auto loc = SubSignalLocation(index);
return filterBank.SubSignal(loc.first, loc.second);
return m_bufferView.SubSignal(loc.first, loc.second);
}
size_t PhaseSize() const noexcept {
return (m_bufferView.Size() + m_filterCount - 1) / m_filterCount;
}
size_t OriginalSize() const noexcept {
return m_bufferView.Size();
}
size_t Size() const {
return (filterBank.Size() + numFilters - 1) / numFilters;
size_t FilterCount() const noexcept {
return m_filterCount;
}

private:
std::pair<size_t, size_t> SubSignalLocation(size_t index) const {
assert(index < numFilters);
const size_t numExtended = filterBank.Size() % numFilters;
const size_t baseFilterSize = filterBank.Size() / numFilters;
assert(index < m_filterCount);
const size_t numExtended = m_bufferView.Size() % m_filterCount;
const size_t baseFilterSize = m_bufferView.Size() / m_filterCount;
const size_t thisFilterSize = baseFilterSize + size_t(index < numExtended);
const size_t offset = baseFilterSize * index + std::min(numExtended, index);
return { offset, thisFilterSize };
}

private:
BasicSignalView<T, Domain> m_bufferView;
size_t m_filterCount = 1;
};

template <class T, eSignalDomain Domain>
void PolyphaseNormalize(PolyphaseDecomposition<T, Domain>& polyphase) {
for (size_t i = 0; i < polyphase.numFilters; ++i) {
class PolyphaseFilter : public PolyphaseView<T, Domain> {
public:
PolyphaseFilter() = default;
PolyphaseFilter(size_t hrFilterSize, size_t numPhases) : m_buffer(hrFilterSize) {
PolyphaseView<T, Domain>::operator=({ AsView(m_buffer), numPhases });
}
PolyphaseFilter(const PolyphaseFilter& rhs) : m_buffer(rhs.m_buffer) {
PolyphaseView<T, Domain>::operator=({ AsView(m_buffer), rhs.FilterCount() });
}
PolyphaseFilter(PolyphaseFilter&& rhs) noexcept : m_buffer(std::move(rhs.m_buffer)) {
PolyphaseView<T, Domain>::operator=(rhs);
rhs.PolyphaseView<T, Domain>::operator=({});
}
PolyphaseFilter& operator=(const PolyphaseFilter& rhs) {
m_buffer = rhs.m_buffer;
PolyphaseView<T, Domain>::operator=({ AsView(m_buffer), rhs.FilterCount() });
return *this;
}
PolyphaseFilter& operator=(PolyphaseFilter&& rhs) noexcept {
m_buffer = std::move(rhs.m_buffer);
PolyphaseView<T, Domain>::operator=(rhs);
rhs.PolyphaseView<T, Domain>::operator=({});
return *this;
}
~PolyphaseFilter() = default;

BasicSignalView<T, Domain> Buffer() {
return AsView(m_buffer);
}
BasicSignalView<const T, Domain> Buffer() const {
return AsView(m_buffer);
}

private:
BasicSignal<T, Domain> m_buffer;
};

template <class T, eSignalDomain Domain>
void PolyphaseNormalize(PolyphaseView<T, Domain>& polyphase) {
for (size_t i = 0; i < polyphase.FilterCount(); ++i) {
polyphase[i] *= T(1) / Sum(polyphase[i]);
}
}

template <class T, eSignalDomain Domain>
auto PolyphaseNormalized(PolyphaseDecomposition<T, Domain>&& polyphase) {
auto PolyphaseNormalized(PolyphaseFilter<T, Domain> polyphase) {
PolyphaseNormalize(polyphase);
return polyphase;
}
Expand All @@ -54,10 +105,9 @@ auto PolyphaseDecompose(SignalR&& output, const SignalT& filter, size_t numFilte
assert(output.Size() == filter.Size());
assert(output.Data() != filter.Data());

PolyphaseDecomposition<typename signal_traits<std::decay_t<SignalR>>::type, signal_traits<std::decay_t<SignalR>>::domain> view{
AsView(output),
numFilters,
};
using R = typename signal_traits<std::decay_t<SignalR>>::type;
constexpr auto Domain = signal_traits<std::decay_t<SignalR>>::domain;
PolyphaseView<R, Domain> view{ AsView(output), numFilters };

for (size_t phaseIdx = 0; phaseIdx < numFilters; ++phaseIdx) {
auto filterPhase = view[phaseIdx];
Expand All @@ -72,5 +122,16 @@ auto PolyphaseDecompose(SignalR&& output, const SignalT& filter, size_t numFilte
return view;
}

template <class SignalT>
auto PolyphaseDecompose(const SignalT& filter, size_t numFilters) {
using R = typename signal_traits<std::decay_t<SignalT>>::type;
constexpr auto Domain = signal_traits<std::decay_t<SignalT>>::domain;
PolyphaseFilter<R, Domain> polyphase{ filter.Size(), numFilters };

PolyphaseDecompose(polyphase.Buffer(), filter, numFilters);

return polyphase;
}


} // namespace dspbb
4 changes: 2 additions & 2 deletions include/dspbb/Math/Convolution.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ using impl::CONV_FULL;
/// <summary> Calculates the length of the result of the convolution U*V. </summary>
/// <param name="lengthU"> Length of U. </param>
/// <param name="lengthV"> Length of V. </param>
inline size_t ConvolutionLength(size_t lengthU, size_t lengthV, impl::ConvCentral) {
constexpr size_t ConvolutionLength(size_t lengthU, size_t lengthV, impl::ConvCentral) {
if (lengthU == 0 || lengthV == 0) {
return 0;
}
Expand All @@ -37,7 +37,7 @@ inline size_t ConvolutionLength(size_t lengthU, size_t lengthV, impl::ConvCentra
/// <summary> Calculates the length of the result of the convolution U*V. </summary>
/// <param name="lengthU"> Length of U. </param>
/// <param name="lengthV"> Length of V. </param>
inline size_t ConvolutionLength(size_t lengthU, size_t lengthV, impl::ConvFull) {
constexpr size_t ConvolutionLength(size_t lengthU, size_t lengthV, impl::ConvFull) {
if (lengthU == 0 || lengthV == 0) {
return 0;
}
Expand Down
104 changes: 87 additions & 17 deletions include/dspbb/Math/Rational.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ class Rational {
Rational& operator/=(const Rational& rhs) noexcept;

template <class FloatT, std::enable_if_t<std::is_floating_point_v<FloatT>, int> = 0>
explicit operator FloatT() const;
constexpr explicit operator FloatT() const;

private:
int_t num;
Expand Down Expand Up @@ -197,7 +197,7 @@ Rational<T>& Rational<T>::operator/=(const Rational& rhs) noexcept {

template <class T>
template <class FloatT, std::enable_if_t<std::is_floating_point_v<FloatT>, int>>
Rational<T>::operator FloatT() const {
constexpr Rational<T>::operator FloatT() const {
return static_cast<FloatT>(num) / static_cast<FloatT>(den);
}

Expand All @@ -220,42 +220,112 @@ constexpr Rational<T> operator/(T lhs, const Rational<T>& rhs) noexcept {
}

namespace impl {
template <class T>
constexpr std::pair<T, T> CommonNumerators(const Rational<T>& lhs, const Rational<T>& rhs) noexcept {
template <class T1, class T2>
constexpr auto CommonNumerators(const Rational<T1>& lhs, const Rational<T2>& rhs) noexcept {
const auto common = std::lcm(lhs.Denominator(), rhs.Denominator());
const auto lhsNum = lhs.Numerator() * (common / lhs.Denominator());
const auto rhsNum = rhs.Numerator() * (common / rhs.Denominator());
return { lhsNum, rhsNum };
return std::pair{ lhsNum, rhsNum };
}
} // namespace impl

template <class T>
constexpr bool operator==(const Rational<T>& lhs, const Rational<T>& rhs) noexcept {

template <class T1, class T2>
constexpr bool operator==(const Rational<T1>& lhs, const Rational<T2>& rhs) noexcept {
const auto [lhsNum, rhsNum] = impl::CommonNumerators(lhs, rhs);
return lhsNum == rhsNum;
}
template <class T>
constexpr bool operator!=(const Rational<T>& lhs, const Rational<T>& rhs) noexcept {
template <class T1, class T2>
constexpr bool operator!=(const Rational<T1>& lhs, const Rational<T2>& rhs) noexcept {
return !(lhs == rhs);
}
template <class T>
constexpr bool operator<(const Rational<T>& lhs, const Rational<T>& rhs) noexcept {
template <class T1, class T2>
constexpr bool operator<(const Rational<T1>& lhs, const Rational<T2>& rhs) noexcept {
const auto [lhsNum, rhsNum] = impl::CommonNumerators(lhs, rhs);
return lhsNum < rhsNum;
}
template <class T>
constexpr bool operator>(const Rational<T>& lhs, const Rational<T>& rhs) noexcept {
template <class T1, class T2>
constexpr bool operator>(const Rational<T1>& lhs, const Rational<T2>& rhs) noexcept {
const auto [lhsNum, rhsNum] = impl::CommonNumerators(lhs, rhs);
return lhsNum > rhsNum;
}
template <class T>
constexpr bool operator<=(const Rational<T>& lhs, const Rational<T>& rhs) noexcept {
template <class T1, class T2>
constexpr bool operator<=(const Rational<T1>& lhs, const Rational<T2>& rhs) noexcept {
return !(lhs > rhs);
}
template <class T>
constexpr bool operator>=(const Rational<T>& lhs, const Rational<T>& rhs) noexcept {
template <class T1, class T2>
constexpr bool operator>=(const Rational<T1>& lhs, const Rational<T2>& rhs) noexcept {
return !(lhs < rhs);
}


template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator==(const Rational<T>& lhs, S rhs) noexcept {
return lhs == Rational{ rhs };
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator!=(const Rational<T>& lhs, S rhs) noexcept {
return !(lhs == rhs);
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator<(const Rational<T>& lhs, S rhs) noexcept {
return lhs < Rational{ rhs };
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator>(const Rational<T>& lhs, S rhs) noexcept {
return lhs > Rational{ rhs };
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator<=(const Rational<T>& lhs, S rhs) noexcept {
return !(lhs > rhs);
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator>=(const Rational<T>& lhs, S rhs) noexcept {
return !(lhs < rhs);
}


template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator==(S lhs, const Rational<T>& rhs) noexcept {
return Rational{ lhs } == rhs;
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator!=(S lhs, const Rational<T>& rhs) noexcept {
return !(lhs == rhs);
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator<(S lhs, const Rational<T>& rhs) noexcept {
return Rational{ lhs } < rhs;
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator>(S lhs, const Rational<T>& rhs) noexcept {
return Rational{ lhs } > rhs;
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator<=(S lhs, const Rational<T>& rhs) noexcept {
return !(lhs > rhs);
}
template <class S, class T, std::enable_if_t<std::is_constructible_v<Rational<T>, S>, int> = 0>
constexpr bool operator>=(S lhs, const Rational<T>& rhs) noexcept {
return !(lhs < rhs);
}



template <class T>
constexpr T ceil(const Rational<T>& rational) {
return (rational.Numerator() + rational.Denominator() - 1) / rational.Denominator();
}

template <class T>
constexpr T floor(const Rational<T>& rational) {
return rational.Numerator() / rational.Denominator();
}

template <class T>
constexpr Rational<T> frac(const Rational<T> rational) {
return rational - floor(rational);
}


} // namespace dspbb
7 changes: 6 additions & 1 deletion include/dspbb/Utility/TypeTraits.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,9 @@
#include <iterator>
#include <type_traits>

#pragma warning(push)
#pragma warning(disable : 4244)

namespace dspbb {

template <class T>
Expand Down Expand Up @@ -150,4 +153,6 @@ template <class Iter>
constexpr bool is_random_access_iterator_v = is_random_access_iterator<Iter>::value;


} // namespace dspbb
} // namespace dspbb

#pragma warning(pop)
Loading

0 comments on commit 2bcf452

Please sign in to comment.