diff --git a/include/dspbb/Filtering/FIR.hpp b/include/dspbb/Filtering/FIR.hpp index 603e8f5..0485906 100644 --- a/include/dspbb/Filtering/FIR.hpp +++ b/include/dspbb/Filtering/FIR.hpp @@ -155,7 +155,7 @@ namespace impl { } // namespace impl -template class Desc, class... Params> +template class Desc, class... Params, std::enable_if_t, int> = 0> auto FirFilter(SignalR&& out, const Desc& desc) -> decltype(void(impl::TranslateLeastSquares(desc))) { const auto [response, weight] = impl::TranslateLeastSquares(desc); diff --git a/include/dspbb/Filtering/Interpolation.hpp b/include/dspbb/Filtering/Interpolation.hpp index 4db8bcd..b0d5e63 100644 --- a/include/dspbb/Filtering/Interpolation.hpp +++ b/include/dspbb/Filtering/Interpolation.hpp @@ -1,18 +1,167 @@ #pragma once +#include "../Generators/Spaces.hpp" #include "../Math/DotProduct.hpp" +#include "../Math/Rational.hpp" #include "../Primitives/Signal.hpp" #include "../Primitives/SignalTraits.hpp" #include "../Primitives/SignalView.hpp" #include "Polyphase.hpp" - namespace dspbb { +//------------------------------------------------------------------------------ +// Public utilities +//------------------------------------------------------------------------------ + +struct InterpSuspensionPoint { + size_t firstInputSample; + size_t startPoint; +}; + + +struct ResamplingSuspensionPoint { + size_t firstInputSample; + Rational startPoint; +}; + + +template +constexpr size_t InterpLength(size_t inputSize, + size_t filterSize, + size_t numPhases, + const ConvType&) { + static_assert(std::is_same_v || std::is_same_v); + + const ptrdiff_t hrInputSize = inputSize * numPhases; + return ConvolutionLength(hrInputSize, filterSize, ConvType{}); +} + + +constexpr double InterpFilterCutoff(size_t numPhases) { + return 1.0 / double(numPhases); +} + + +template +constexpr Rational ResamplingLength(size_t inputSize, + size_t filterSize, + size_t numPhases, + Rational sampleRates, + const ConvType&) { + static_assert(std::is_same_v || std::is_same_v); + const int64_t interpolatedSize = int64_t(numPhases) * inputSize; + const int64_t filteredInterpolatedSize = ConvolutionLength(interpolatedSize, filterSize, ConvType{}); + + return filteredInterpolatedSize / sampleRates / int64_t(numPhases); +} + + +constexpr double ResamplingFilterCutoff(Rational sampleRates, size_t numPhases) { + const double base = 1.0 / double(numPhases); + const double rate = std::min(1.0, 1.0 / double(sampleRates)); + return base * rate; +} + + +constexpr Rational ResamplingDelay(size_t filterSize, + size_t numPhases, + Rational sampleRates) { + return Rational{ int64_t(filterSize) - 1, 2 * int64_t(numPhases) } / sampleRates; +} + + +//------------------------------------------------------------------------------ +// Internal utilities +//------------------------------------------------------------------------------ + +namespace impl { + inline InterpSuspensionPoint FindInterpSuspensionPoint(size_t nextOutputSample, size_t filterSize, size_t numPhases) { + const ptrdiff_t firstOutputSample = ptrdiff_t(nextOutputSample) - ptrdiff_t(filterSize - 1); + if (firstOutputSample < 0) { + return { 0, nextOutputSample }; + } + + const size_t firstInputSample = firstOutputSample / numPhases; + const size_t startPoint = firstOutputSample - (numPhases * firstInputSample) + (filterSize - 1); + + return { firstInputSample, startPoint }; + } + + constexpr Rational ChangeSampleRate(int64_t sourceRate, + int64_t targetRate, + Rational sample) { + return sample * Rational{ targetRate, sourceRate }; + } + + constexpr ResamplingSuspensionPoint FindResamplingSuspensionPoint(Rational nextOutputSample, + size_t filterSize, + size_t numPhases, + Rational sampleRates) { + const auto nextInputSample = ChangeSampleRate(sampleRates.Denominator(), sampleRates.Numerator(), nextOutputSample); + const auto convolutionOffset = Rational{ int64_t(filterSize) - 1, int64_t(numPhases) }; + const auto firstInputSample = nextInputSample - convolutionOffset; + + if (firstInputSample <= 0) { + return { 0, nextOutputSample }; + } + else { + const size_t firstInputSampleWhole = floor(firstInputSample); + const auto inputStartPoint = frac(firstInputSample) + convolutionOffset; + const auto outputStartPoint = ChangeSampleRate(sampleRates.Numerator(), sampleRates.Denominator(), inputStartPoint); + return { firstInputSampleWhole, outputStartPoint }; + } + } + + struct PhaseSample { + size_t inputIndex; + size_t phaseIndex; + uint64_t weight; + }; + + constexpr std::pair InputIndex2Sample(Rational inputIndex, size_t numPhases) { + const Rational indexFrac = frac(inputIndex); + + const size_t firstPhase = floor(indexFrac * int64_t(numPhases)); + const size_t secondPhase = (firstPhase + 1) % numPhases; + + const Rational t = frac(indexFrac * int64_t(numPhases)); + const size_t secondWeight = t.Numerator(); + const size_t firstWeight = t.Denominator() - t.Numerator(); + + const size_t firstIndex = floor(inputIndex); + const size_t secondIndex = secondPhase == 0 ? firstIndex + 1 : firstIndex; + + return { + PhaseSample{ firstIndex, firstPhase, firstWeight }, + PhaseSample{ secondIndex, secondPhase, secondWeight } + }; + } + + template + auto DotProductSample(const SignalT& input, const SignalU& filter, size_t inputReverseFirst) { + const ptrdiff_t desiredFirst = ptrdiff_t(inputReverseFirst) - filter.Size() + 1; + const ptrdiff_t desiredLast = ptrdiff_t(inputReverseFirst) + 1; + const ptrdiff_t possibleFirst = std::max(ptrdiff_t(0), desiredFirst); + const ptrdiff_t possibleLast = std::min(ptrdiff_t(input.Size()), desiredLast); + const ptrdiff_t count = possibleLast - possibleFirst; + assert(count >= 0); + const ptrdiff_t offset = possibleFirst - desiredFirst; + + const auto inputView = AsConstView(input).SubSignal(possibleFirst, count); + const auto filterView = AsConstView(filter).SubSignal(offset, count); + return DotProduct(inputView, filterView); + } + + +} // namespace impl + + +//------------------------------------------------------------------------------ +// Expansion & Interpolation & Resampling +//------------------------------------------------------------------------------ + -/// -/// Erases all but every th sample. -/// template && is_mutable_signal_v, int> = 0> @@ -27,6 +176,7 @@ void Decimate(SignalR&& output, } } + template , int> = 0> auto Decimate(const SignalT& input, size_t rate) { using T = std::remove_const_t::type>; @@ -37,10 +187,6 @@ auto Decimate(const SignalT& input, size_t rate) { } -/// -/// Inserts zeros between samples to increase sample rate by a factor of . -/// -/// Follow expansion by a low-pass filter to interpolate a signal. template && is_mutable_signal_v, int> = 0> @@ -69,106 +215,112 @@ auto Expand(const SignalT& input, size_t rate) { } -/// -/// Inserts meaningful samples to increase sample rate by a factor of .numFilters. -/// -/// A polyphase decomposition of an appropriate low-pass filter. -/// The number of phases defines the interpolation ratio. -/// No need to follow up with low-pass filtering. -/// The polyphase filter must have the appropriate cutoff-frequency of -/// (input sample rate / 2), and the polyphase filter must operate at -/// the output sample rate. template > && is_mutable_signal_v, int> = 0> -void Interpolate(SignalR&& output, - const SignalT& input, - const PolyphaseDecomposition& polyphase, - intptr_t offset) { - const auto rate = polyphase.numFilters; - //assert(output.Size() == input.Size() * rate); - size_t count = output.Size() / rate; +InterpSuspensionPoint Interpolate(SignalR&& hrOutput, + const SignalT& lrInput, + const PolyphaseView& polyphase, + size_t hrOffset) { + const ptrdiff_t rate = polyphase.FilterCount(); + const ptrdiff_t hrFilterSize = polyphase.OriginalSize(); + const ptrdiff_t lrPhaseSize = polyphase.PhaseSize(); + const ptrdiff_t hrOutputSize = hrOutput.Size(); - auto writeIt = output.begin(); - for (size_t inputIdx = 0; inputIdx < count; ++inputIdx) { - for (size_t i = 0; i < rate; ++i, ++writeIt) { - const auto& filter = polyphase[i]; - intptr_t offsetInputIdx = intptr_t(inputIdx) - intptr_t(filter.Size()) + 1 + intptr_t(offset); - const intptr_t inputIndexClamped = std::max(offsetInputIdx, intptr_t(0)); - const auto subInput = AsConstView(input).SubSignal(inputIndexClamped); - const auto subFilter = filter.SubSignal(-std::min(intptr_t(0), offsetInputIdx)); - const size_t dotLength = std::min(subFilter.Size(), subInput.Size()); - *writeIt = DotProduct(subFilter.SubSignal(0, dotLength), subInput.SubSignal(0, dotLength)); + const ptrdiff_t hrOutputMaxSize = InterpLength(lrInput.Size(), hrFilterSize, rate, CONV_FULL); + assert(ptrdiff_t(hrOffset) + hrOutputSize <= hrOutputMaxSize); + + size_t hrOutputIdx = hrOffset; + for (; hrOutputIdx < hrOffset + hrOutputSize; ++hrOutputIdx) { + const ptrdiff_t hrInputIdx = 1 - hrFilterSize + hrOutputIdx; + const ptrdiff_t lrInputIdx = (hrInputIdx + hrFilterSize - 1) / rate - lrPhaseSize + 1; + const ptrdiff_t polyphaseIdx = (hrInputIdx + hrFilterSize - 1) % rate; + + const auto& phase = polyphase[polyphaseIdx]; + + const Interval inputSpan = { ptrdiff_t(0), ptrdiff_t(lrInput.Size()) }; + const Interval lrInputInterval = { lrInputIdx, lrInputIdx + lrPhaseSize }; + const Interval lrPhaseInterval = { lrInputInterval.last - ptrdiff_t(phase.Size()), lrInputInterval.last }; + const Interval lrInputProductInterval = Intersection(inputSpan, Intersection(lrInputInterval, lrPhaseInterval)); + const Interval lrPhaseProductInterval = lrInputProductInterval - lrInputIdx; + + if (lrInputProductInterval.Size() > 0) { + const auto lrInputView = AsView(lrInput).SubSignal(lrInputProductInterval.first, + lrInputProductInterval.last - lrInputProductInterval.first); + const auto lrPhaseView = phase.SubSignal(lrPhaseProductInterval.first - lrPhaseSize + ptrdiff_t(phase.Size()), + lrPhaseProductInterval.last - lrPhaseProductInterval.first); + const auto value = DotProduct(lrInputView, lrPhaseView); + hrOutput[hrOutputIdx - hrOffset] = value; } } + + return impl::FindInterpSuspensionPoint(hrOutputIdx, polyphase.OriginalSize(), polyphase.FilterCount()); +} + + +template >, int> = 0> +auto Interpolate(const SignalT& lrInput, + const PolyphaseView& polyphase, + size_t hrOffset, + size_t hrLength) { + using T = typename signal_traits>::type; + using R = multiplies_result_t; + + BasicSignal out(hrLength, R(0)); + Interpolate(out, lrInput, polyphase, hrOffset); + return out; } -/// -/// Arbitrary rational resampling of a signal using approximate polyphase interpolation. -/// -/// A polyphase decomposition of an appropriate low-pass filter. -/// The number of phases is arbitrary (see remarks), but must be at least 2. -/// A pair of the {input, output} sample rates. -/// For example, use {16000, 44100} if you want to increase the sample rate of a signal by 2.75625 times. -/// A rational number that tells where to take the first output sample. For example, specify {1, 2} -/// so that the first output sample will be taken from halfway between the first and second input samples. Specify {3, 77} -/// to take the first output sample from just after the first input sample. The denominator is arbitrary, the numerator -/// must be larger than zero, and float(num)/float(den) must be smaller or equal to length()-1. -/// A rational number that tells from where the next output sample would have been taken. -/// Analogues in meaning to , but its denominator is not necessarily the same. -/// When resampling streams in batches, you can use this as a base for the to a subsequent Resample call to ensure -/// the output stream is continuous. template > && is_mutable_signal_v, int> = 0> -std::pair Resample(SignalR&& output, - const SignalT& input, - const PolyphaseDecomposition& polyphase, - std::pair sampleRates, - std::pair startPoint = { 0, 1 }) { - using R = typename signal_traits>::type; - assert(sampleRates.first > 0); - assert(sampleRates.second > 0); - assert(startPoint.second > 0); - assert(polyphase.numFilters > 0); - const size_t commonRate = std::lcm(std::lcm(std::lcm(sampleRates.first, sampleRates.second), startPoint.second), polyphase.numFilters); - - // All these are in common rate. - const int64_t commonStartPoint = startPoint.first * commonRate / startPoint.second; - const uint64_t commonStep = sampleRates.first * commonRate / sampleRates.second; - int64_t commonSample = commonStartPoint; - - for (auto& o : output) { - const int64_t commonFraction = commonSample % commonRate; +ResamplingSuspensionPoint Resample(SignalR&& output, + const SignalT& input, + const PolyphaseView& polyphase, + Rational sampleRates, + Rational startPoint = { 0, 1 }) { + assert(sampleRates >= 0ll); + assert(startPoint >= 0ll); + assert(polyphase.FilterCount() > 0); - const int64_t polyphaseIndex = commonFraction * polyphase.numFilters / commonRate; - const int64_t polyphaseIndexNext = (polyphaseIndex + 1) % polyphase.numFilters; - const int64_t polyphaseFraction = commonFraction * polyphase.numFilters % commonRate; + [[maybe_unused]] const auto maxLength = ResamplingLength(input.Size(), polyphase.OriginalSize(), polyphase.FilterCount(), sampleRates, CONV_FULL); + assert(startPoint + int64_t(output.Size()) < maxLength); - const int64_t inputIndex = commonSample / commonRate - polyphase[polyphaseIndex].Size() + 1; - const int64_t inputIndexNext = commonSample / commonRate - int64_t(polyphase[polyphaseIndexNext].Size()) + 1 + int64_t(polyphaseIndexNext == 0); + auto outputIndex = startPoint; + for (auto outputIt = output.begin(); outputIt != output.end(); ++outputIt, outputIndex += 1) { + const auto inputIndex = impl::ChangeSampleRate(sampleRates.Denominator(), sampleRates.Numerator(), outputIndex); + const auto [firstSampleLoc, secondSampleLoc] = impl::InputIndex2Sample(inputIndex, polyphase.FilterCount()); + const auto firstSampleVal = impl::DotProductSample(input, polyphase[firstSampleLoc.phaseIndex], firstSampleLoc.inputIndex); + const auto secondSampleVal = impl::DotProductSample(input, polyphase[secondSampleLoc.phaseIndex], secondSampleLoc.inputIndex); + using CommonType = decltype(firstSampleVal); + *outputIt = (firstSampleVal * CommonType(firstSampleLoc.weight) + secondSampleVal * CommonType(secondSampleLoc.weight)) + / (CommonType(firstSampleLoc.weight) + CommonType(secondSampleLoc.weight)); + } - int64_t inputIndexClamped = std::max(inputIndex, int64_t(0)); - int64_t inputIndexClampedNext = std::max(inputIndexNext, int64_t(0)); + return impl::FindResamplingSuspensionPoint(outputIndex, polyphase.OriginalSize(), polyphase.FilterCount(), sampleRates); +} - const auto filter = polyphase[polyphaseIndex].SubSignal(inputIndexClamped - inputIndex); - const auto filterNext = polyphase[polyphaseIndexNext].SubSignal(inputIndexClampedNext - inputIndexNext); - const auto inputSection = AsView(input).SubSignal(inputIndexClamped); - const auto inputSectionNext = AsView(input).SubSignal(inputIndexClampedNext); - const auto dotLength = std::min(inputSection.Size(), filter.Size()); - const auto sample = DotProduct(inputSection.SubSignal(0, dotLength), filter.SubSignal(0, dotLength)); - const auto dotLengthNext = std::min(inputSectionNext.Size(), filterNext.Size()); - const auto sampleNext = DotProduct(inputSectionNext.SubSignal(0, dotLengthNext), filterNext.SubSignal(0, dotLengthNext)); - o = (sample * R(commonRate - polyphaseFraction) + sampleNext * R(polyphaseFraction)) / R(commonRate); +template >, int> = 0> +auto Resample(const SignalT& input, + const PolyphaseView& polyphase, + Rational sampleRates, + Rational startPoint, + size_t outputLength) { + using T = typename signal_traits>::type; + using R = multiplies_result_t; - // Next sample - commonSample += commonStep; - } - return { commonSample, commonRate }; + BasicSignal out(outputLength, R(0)); + Resample(out, input, polyphase, sampleRates, startPoint); + return out; } diff --git a/include/dspbb/Filtering/Polyphase.hpp b/include/dspbb/Filtering/Polyphase.hpp index d247e85..62a20a9 100644 --- a/include/dspbb/Filtering/Polyphase.hpp +++ b/include/dspbb/Filtering/Polyphase.hpp @@ -10,41 +10,92 @@ namespace dspbb { template -struct PolyphaseDecomposition { - BasicSignalView filterBank; - size_t numFilters; +class PolyphaseView { +public: + PolyphaseView() noexcept = default; + PolyphaseView(BasicSignalView data, size_t numFilters) noexcept + : m_bufferView(data), m_filterCount(numFilters) {} + BasicSignalView operator[](size_t index) { auto loc = SubSignalLocation(index); - return filterBank.SubSignal(loc.first, loc.second); + return m_bufferView.SubSignal(loc.first, loc.second); } BasicSignalView 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 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 m_bufferView; + size_t m_filterCount = 1; }; template -void PolyphaseNormalize(PolyphaseDecomposition& polyphase) { - for (size_t i = 0; i < polyphase.numFilters; ++i) { +class PolyphaseFilter : public PolyphaseView { +public: + PolyphaseFilter() = default; + PolyphaseFilter(size_t hrFilterSize, size_t numPhases) : m_buffer(hrFilterSize) { + PolyphaseView::operator=({ AsView(m_buffer), numPhases }); + } + PolyphaseFilter(const PolyphaseFilter& rhs) : m_buffer(rhs.m_buffer) { + PolyphaseView::operator=({ AsView(m_buffer), rhs.FilterCount() }); + } + PolyphaseFilter(PolyphaseFilter&& rhs) noexcept : m_buffer(std::move(rhs.m_buffer)) { + PolyphaseView::operator=(rhs); + rhs.PolyphaseView::operator=({}); + } + PolyphaseFilter& operator=(const PolyphaseFilter& rhs) { + m_buffer = rhs.m_buffer; + PolyphaseView::operator=({ AsView(m_buffer), rhs.FilterCount() }); + return *this; + } + PolyphaseFilter& operator=(PolyphaseFilter&& rhs) noexcept { + m_buffer = std::move(rhs.m_buffer); + PolyphaseView::operator=(rhs); + rhs.PolyphaseView::operator=({}); + return *this; + } + ~PolyphaseFilter() = default; + + BasicSignalView Buffer() { + return AsView(m_buffer); + } + BasicSignalView Buffer() const { + return AsView(m_buffer); + } + +private: + BasicSignal m_buffer; +}; + +template +void PolyphaseNormalize(PolyphaseView& polyphase) { + for (size_t i = 0; i < polyphase.FilterCount(); ++i) { polyphase[i] *= T(1) / Sum(polyphase[i]); } } template -auto PolyphaseNormalized(PolyphaseDecomposition&& polyphase) { +auto PolyphaseNormalized(PolyphaseFilter polyphase) { PolyphaseNormalize(polyphase); return polyphase; } @@ -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>::type, signal_traits>::domain> view{ - AsView(output), - numFilters, - }; + using R = typename signal_traits>::type; + constexpr auto Domain = signal_traits>::domain; + PolyphaseView view{ AsView(output), numFilters }; for (size_t phaseIdx = 0; phaseIdx < numFilters; ++phaseIdx) { auto filterPhase = view[phaseIdx]; @@ -72,5 +122,16 @@ auto PolyphaseDecompose(SignalR&& output, const SignalT& filter, size_t numFilte return view; } +template +auto PolyphaseDecompose(const SignalT& filter, size_t numFilters) { + using R = typename signal_traits>::type; + constexpr auto Domain = signal_traits>::domain; + PolyphaseFilter polyphase{ filter.Size(), numFilters }; + + PolyphaseDecompose(polyphase.Buffer(), filter, numFilters); + + return polyphase; +} + } // namespace dspbb \ No newline at end of file diff --git a/include/dspbb/Math/Convolution.hpp b/include/dspbb/Math/Convolution.hpp index b07531f..c1d9156 100644 --- a/include/dspbb/Math/Convolution.hpp +++ b/include/dspbb/Math/Convolution.hpp @@ -24,7 +24,7 @@ using impl::CONV_FULL; /// Calculates the length of the result of the convolution U*V. /// Length of U. /// Length of V. -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; } @@ -37,7 +37,7 @@ inline size_t ConvolutionLength(size_t lengthU, size_t lengthV, impl::ConvCentra /// Calculates the length of the result of the convolution U*V. /// Length of U. /// Length of V. -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; } diff --git a/include/dspbb/Math/Rational.hpp b/include/dspbb/Math/Rational.hpp index dd90181..b87a102 100644 --- a/include/dspbb/Math/Rational.hpp +++ b/include/dspbb/Math/Rational.hpp @@ -40,7 +40,7 @@ class Rational { Rational& operator/=(const Rational& rhs) noexcept; template , int> = 0> - explicit operator FloatT() const; + constexpr explicit operator FloatT() const; private: int_t num; @@ -197,7 +197,7 @@ Rational& Rational::operator/=(const Rational& rhs) noexcept { template template , int>> -Rational::operator FloatT() const { +constexpr Rational::operator FloatT() const { return static_cast(num) / static_cast(den); } @@ -220,42 +220,112 @@ constexpr Rational operator/(T lhs, const Rational& rhs) noexcept { } namespace impl { - template - constexpr std::pair CommonNumerators(const Rational& lhs, const Rational& rhs) noexcept { + template + constexpr auto CommonNumerators(const Rational& lhs, const Rational& 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 -constexpr bool operator==(const Rational& lhs, const Rational& rhs) noexcept { + +template +constexpr bool operator==(const Rational& lhs, const Rational& rhs) noexcept { const auto [lhsNum, rhsNum] = impl::CommonNumerators(lhs, rhs); return lhsNum == rhsNum; } -template -constexpr bool operator!=(const Rational& lhs, const Rational& rhs) noexcept { +template +constexpr bool operator!=(const Rational& lhs, const Rational& rhs) noexcept { return !(lhs == rhs); } -template -constexpr bool operator<(const Rational& lhs, const Rational& rhs) noexcept { +template +constexpr bool operator<(const Rational& lhs, const Rational& rhs) noexcept { const auto [lhsNum, rhsNum] = impl::CommonNumerators(lhs, rhs); return lhsNum < rhsNum; } -template -constexpr bool operator>(const Rational& lhs, const Rational& rhs) noexcept { +template +constexpr bool operator>(const Rational& lhs, const Rational& rhs) noexcept { const auto [lhsNum, rhsNum] = impl::CommonNumerators(lhs, rhs); return lhsNum > rhsNum; } -template -constexpr bool operator<=(const Rational& lhs, const Rational& rhs) noexcept { +template +constexpr bool operator<=(const Rational& lhs, const Rational& rhs) noexcept { return !(lhs > rhs); } -template -constexpr bool operator>=(const Rational& lhs, const Rational& rhs) noexcept { +template +constexpr bool operator>=(const Rational& lhs, const Rational& rhs) noexcept { + return !(lhs < rhs); +} + + +template , S>, int> = 0> +constexpr bool operator==(const Rational& lhs, S rhs) noexcept { + return lhs == Rational{ rhs }; +} +template , S>, int> = 0> +constexpr bool operator!=(const Rational& lhs, S rhs) noexcept { + return !(lhs == rhs); +} +template , S>, int> = 0> +constexpr bool operator<(const Rational& lhs, S rhs) noexcept { + return lhs < Rational{ rhs }; +} +template , S>, int> = 0> +constexpr bool operator>(const Rational& lhs, S rhs) noexcept { + return lhs > Rational{ rhs }; +} +template , S>, int> = 0> +constexpr bool operator<=(const Rational& lhs, S rhs) noexcept { + return !(lhs > rhs); +} +template , S>, int> = 0> +constexpr bool operator>=(const Rational& lhs, S rhs) noexcept { + return !(lhs < rhs); +} + + +template , S>, int> = 0> +constexpr bool operator==(S lhs, const Rational& rhs) noexcept { + return Rational{ lhs } == rhs; +} +template , S>, int> = 0> +constexpr bool operator!=(S lhs, const Rational& rhs) noexcept { + return !(lhs == rhs); +} +template , S>, int> = 0> +constexpr bool operator<(S lhs, const Rational& rhs) noexcept { + return Rational{ lhs } < rhs; +} +template , S>, int> = 0> +constexpr bool operator>(S lhs, const Rational& rhs) noexcept { + return Rational{ lhs } > rhs; +} +template , S>, int> = 0> +constexpr bool operator<=(S lhs, const Rational& rhs) noexcept { + return !(lhs > rhs); +} +template , S>, int> = 0> +constexpr bool operator>=(S lhs, const Rational& rhs) noexcept { return !(lhs < rhs); } + +template +constexpr T ceil(const Rational& rational) { + return (rational.Numerator() + rational.Denominator() - 1) / rational.Denominator(); +} + +template +constexpr T floor(const Rational& rational) { + return rational.Numerator() / rational.Denominator(); +} + +template +constexpr Rational frac(const Rational rational) { + return rational - floor(rational); +} + + } // namespace dspbb \ No newline at end of file diff --git a/include/dspbb/Utility/TypeTraits.hpp b/include/dspbb/Utility/TypeTraits.hpp index 13d3e36..afad6c3 100644 --- a/include/dspbb/Utility/TypeTraits.hpp +++ b/include/dspbb/Utility/TypeTraits.hpp @@ -6,6 +6,9 @@ #include #include +#pragma warning(push) +#pragma warning(disable : 4244) + namespace dspbb { template @@ -150,4 +153,6 @@ template constexpr bool is_random_access_iterator_v = is_random_access_iterator::value; -} // namespace dspbb \ No newline at end of file +} // namespace dspbb + +#pragma warning(pop) \ No newline at end of file diff --git a/test/Filtering/Test_Interpolation.cpp b/test/Filtering/Test_Interpolation.cpp index 6a8416e..4aaeafb 100644 --- a/test/Filtering/Test_Interpolation.cpp +++ b/test/Filtering/Test_Interpolation.cpp @@ -1,3 +1,5 @@ +#include "../TestUtils.hpp" + #include #include @@ -6,14 +8,13 @@ using namespace dspbb; -auto MakeRamp(size_t size) { - Signal signal; - for (size_t i = 0; i < size; ++i) { - signal.PushBack(float(i)); - } - return signal; + +template +SignalT InterpolateRefImpl(const SignalT& signal, const SignalT& filter, size_t rate, size_t offset, size_t length) { + return Convolution(Expand(signal, rate), filter, offset, length) * rate; } + TEST_CASE("Decimate", "[Interpolation]") { const Signal s = { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 }; const Signal d = Decimate(s, 3); @@ -34,159 +35,514 @@ TEST_CASE("Expand", "[Interpolation]") { REQUIRE(Max(Abs(e - exp)) == Approx(0.0f)); } +TEST_CASE("Interpolation full", "[Interpolation]") { + constexpr int interpRate = 5; + constexpr int signalSize = 1024; + + for (const int filterSize : { 31, 33, 2047 }) { + const auto signal = RandomSignal(signalSize); + const auto filter = FirFilter(filterSize, Lowpass(WINDOWED).Cutoff(1.0f / interpRate)); + const auto polyphase = PolyphaseDecompose(filter, interpRate); + + const auto length = ConvolutionLength(signal.Size() * interpRate, filter.Size(), CONV_FULL); + const auto reference = InterpolateRefImpl(signal, filter, interpRate, 0, length); + const auto answer = Interpolate(signal, polyphase, 0, length); + + INFO("filterSize=" << filterSize) + REQUIRE(reference.Size() == answer.Size()); + REQUIRE(Max(Abs(reference - answer)) < 1e-6f); + } +} + +TEST_CASE("Interpolation central", "[Interpolation]") { + constexpr int interpRate = 5; + constexpr int signalSize = 1024; -TEST_CASE("Polyphase interpolation", "[Interpolation]") { - constexpr int numFilters = 3; - const Signal filter = FirFilter(31, Lowpass(WINDOWED).Cutoff(1.0f / numFilters)); - Signal scratch(filter.Size()); - const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(scratch, filter, numFilters)); + for (const int filterSize : { 31, 33, 2047 }) { + const auto signal = RandomSignal(signalSize); + const auto filter = FirFilter(filterSize, Lowpass(WINDOWED).Cutoff(1.0f / interpRate)); + const auto polyphase = PolyphaseDecompose(filter, interpRate); - const Signal signal = MakeRamp(150); - Signal output(numFilters * ConvolutionLength(signal.Size(), polyphase.Size(), CONV_FULL)); - Interpolate(output, signal, polyphase, 0); + const auto length = ConvolutionLength(signal.Size() * interpRate, filter.Size(), CONV_CENTRAL); + const auto reference = InterpolateRefImpl(signal, filter, interpRate, filterSize - 1, length); + const auto answer = Interpolate(signal, polyphase, filterSize - 1, length); - REQUIRE(output[0] == Approx(0).margin(0.001f)); + INFO("filterSize=" << filterSize) + REQUIRE(reference.Size() == answer.Size()); + REQUIRE(Max(Abs(reference - answer)) < 1e-6f); + } +} - const SignalView frontView{ output.begin() + 100, output.begin() + 350 }; - const SignalView backView{ output.begin() + 101, output.begin() + 351 }; - const auto diff = backView - frontView; - REQUIRE(Mean(diff) == Approx(1.0f / numFilters).epsilon(0.001f)); - REQUIRE(Min(diff) == Approx(1.0f / numFilters).epsilon(0.02f)); - REQUIRE(Max(diff) == Approx(1.0f / numFilters).epsilon(0.02f)); + +TEST_CASE("Resampling length full", "[Interpolation]") { + SECTION("Upsample exact") { + constexpr Rational sampleRates = { 2, 3 }; + constexpr size_t signalSize = 2000; + constexpr size_t filterSize = 1001; + constexpr size_t numPhases = 5; + + constexpr auto size = ResamplingLength(signalSize, filterSize, numPhases, sampleRates, CONV_FULL); + REQUIRE(double(size) == Approx(16500.0 / 5).margin(0.01)); + } + SECTION("Upsample inexact") { + constexpr Rational sampleRates = { 3, 5 }; + constexpr size_t signalSize = 2000; + constexpr size_t filterSize = 1001; + constexpr size_t numPhases = 5; + + constexpr auto size = ResamplingLength(signalSize, filterSize, numPhases, sampleRates, CONV_FULL); + REQUIRE(double(size) == Approx(18333.333 / 5).margin(0.01)); + } + SECTION("Downsample exact") { + constexpr Rational sampleRates = { 11000, 3500 }; + constexpr size_t signalSize = 2000; + constexpr size_t filterSize = 1001; + constexpr size_t numPhases = 5; + + const auto size = ResamplingLength(signalSize, filterSize, numPhases, sampleRates, CONV_FULL); + REQUIRE(double(size) == Approx(3500.0 / 5).margin(0.01)); + } + SECTION("Downsample inexact") { + constexpr Rational sampleRates = { 22000, 7001 }; + constexpr size_t signalSize = 2000; + constexpr size_t filterSize = 1001; + constexpr size_t numPhases = 5; + + constexpr auto size = ResamplingLength(signalSize, filterSize, numPhases, sampleRates, CONV_FULL); + REQUIRE(double(size) == Approx(3500.5 / 5).margin(0.01)); + } +} + +TEST_CASE("Resampling length central", "[Interpolation]") { + SECTION("Upsample exact") { + constexpr Rational sampleRates = { 9000, 14000 }; + constexpr size_t signalSize = 2000; + constexpr size_t filterSize = 1001; + constexpr size_t numPhases = 5; + + constexpr auto size = ResamplingLength(signalSize, filterSize, numPhases, sampleRates, CONV_CENTRAL); + REQUIRE(double(size) == Approx(14000.0 / 5).margin(0.01)); + } + SECTION("Upsample inexact") { + constexpr Rational sampleRates = { 27000, 14000 }; + constexpr size_t signalSize = 2000; + constexpr size_t filterSize = 1001; + constexpr size_t numPhases = 5; + + constexpr auto size = ResamplingLength(signalSize, filterSize, numPhases, sampleRates, CONV_CENTRAL); + REQUIRE(double(size) == Approx(4666.667 / 5).margin(0.01)); + } + SECTION("Downsample exact") { + constexpr Rational sampleRates = { 9000, 3500 }; + constexpr size_t signalSize = 2000; + constexpr size_t filterSize = 1001; + constexpr size_t numPhases = 5; + + const auto size = ResamplingLength(signalSize, filterSize, numPhases, sampleRates, CONV_CENTRAL); + REQUIRE(double(size) == Approx(3500.0 / 5).margin(0.01)); + } + SECTION("Downsample inexact") { + constexpr Rational sampleRates = { 18000, 7001 }; + constexpr size_t signalSize = 2000; + constexpr size_t filterSize = 1001; + constexpr size_t numPhases = 5; + + constexpr auto size = ResamplingLength(signalSize, filterSize, numPhases, sampleRates, CONV_CENTRAL); + REQUIRE(double(size) == Approx(3500.5 / 5).margin(0.01)); + } } +TEST_CASE("Resampling change sample rate", "[Interpolation]") { + constexpr int inputRate = 7; + constexpr int outputRate = 17; + + constexpr Rational originalSample = { 28, 42 }; -TEST_CASE("Polyphase resampling replicate convolution full", "[Interpolation]") { - constexpr int numFilters = 4; - const Signal filter = FirFilter(31, Lowpass(WINDOWED).Cutoff(1.0f / numFilters)); - Signal scratch(filter.Size()); - const auto polyphase = PolyphaseDecompose(scratch, filter, numFilters); + SECTION("Regular") { + constexpr auto newSample = impl::ChangeSampleRate(inputRate, outputRate, originalSample); - const Signal signal = MakeRamp(150); - Signal padded(signal.Size() * numFilters); - Expand(padded, signal, numFilters); + const double inputIndexRealExpected = double(originalSample) / double(inputRate) * double(outputRate); - Signal outputConv = Convolution(padded, filter * numFilters, CONV_FULL); - Signal output(outputConv.Size()); + REQUIRE(double(newSample) == Approx(inputIndexRealExpected)); + } + SECTION("Simplify") { + constexpr auto newSample = impl::ChangeSampleRate(inputRate, outputRate, originalSample); - Resample(output, signal, polyphase, { 1, numFilters }); + const double inputIndexRealExpected = double(originalSample) / double(inputRate) * double(outputRate); - REQUIRE(Max(Abs(output - outputConv)) < 0.0001f); + REQUIRE(double(newSample) == Approx(inputIndexRealExpected)); + } } +TEST_CASE("Resampling input index 2 samples", "[Interpolation]") { + SECTION("Zero weight") { + const auto [firstSample, secondSample] = impl::InputIndex2Sample({ 43, 7 }, 7); + REQUIRE(firstSample.inputIndex == 6); + REQUIRE(firstSample.phaseIndex == 1); + REQUIRE(firstSample.weight == 1); + + REQUIRE(secondSample.inputIndex == 6); + REQUIRE(secondSample.phaseIndex == 2); + REQUIRE(secondSample.weight == 0); + } + SECTION("Split weight") { + const auto [firstSample, secondSample] = impl::InputIndex2Sample({ 87, 14 }, 5); + REQUIRE(firstSample.inputIndex == 6); + REQUIRE(firstSample.phaseIndex == 1); + REQUIRE(firstSample.weight == 13); + + REQUIRE(secondSample.inputIndex == 6); + REQUIRE(secondSample.phaseIndex == 2); + REQUIRE(secondSample.weight == 1); + } + SECTION("Rollover") { + const auto [firstSample, secondSample] = impl::InputIndex2Sample({ 27, 14 }, 5); + REQUIRE(firstSample.inputIndex == 1); + REQUIRE(firstSample.phaseIndex == 4); + REQUIRE(firstSample.weight == 5); + + REQUIRE(secondSample.inputIndex == 2); + REQUIRE(secondSample.phaseIndex == 0); + REQUIRE(secondSample.weight == 9); + } +} -TEST_CASE("Polyphase resampling upsample constant", "[Interpolation]") { - constexpr int numFilters = 4; - const Signal filter = FirFilter(31, Lowpass(WINDOWED).Cutoff(1.0f / numFilters)); - Signal scratch(filter.Size()); - const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(scratch, filter, numFilters)); +TEST_CASE("Resampling dot product sample", "[Interpolation]") { + const Signal signal = { 1, 2, 3, 6, 5, 7 }; + const Signal filter = { -1, 3, -2 }; // Convolution: -2, 3, -1 + REQUIRE(-2 == impl::DotProductSample(signal, filter, 0)); + REQUIRE(-1 == impl::DotProductSample(signal, filter, 2)); + REQUIRE(-5 == impl::DotProductSample(signal, filter, 5)); + REQUIRE(-7 == impl::DotProductSample(signal, filter, 7)); +} - Signal signal(150, 1.0f); - Signal output(signal.Size() - polyphase[0].Size() - 1); +TEST_CASE("Resampling filter cutoff", "[Interpolation]") { + REQUIRE(ResamplingFilterCutoff({ 4, 6 }, 5) == Approx(0.2)); + REQUIRE(ResamplingFilterCutoff({ 6, 4 }, 5) == Approx(0.1333333333)); + REQUIRE(ResamplingFilterCutoff({ 4, 71 }, 12) == Approx(0.0833333333)); + REQUIRE(ResamplingFilterCutoff({ 40, 6 }, 12) == Approx(0.0125)); +} - Resample(output, signal, polyphase, { 7, 11 }, { filter.Size() * 100 + 62 * numFilters, numFilters * 100 }); - REQUIRE(Min(output) == Approx(1.f)); - REQUIRE(Max(output) == Approx(1.f)); +template +auto ResampledSimilarity(std::pair sampleRates, SignalT original, SignalU resampled) { + const size_t rescale = std::max(original.Size() / sampleRates.first, resampled.Size() / sampleRates.second) + 1; + original.Resize(rescale * sampleRates.first); + resampled.Resize(rescale * sampleRates.second); + + const auto fftSignal = Abs(Fft(original, FFT_HALF)); + const auto fftResampled = Abs(Fft(resampled, FFT_HALF)); + + const size_t fftCompareSize = std::min(fftSignal.Size(), fftResampled.Size()); + const auto fftSignalCompare = AsView(fftSignal).SubSignal(0, fftCompareSize); + const auto fftResampledCompare = AsView(fftResampled).SubSignal(0, fftCompareSize); + + const auto similarity = DotProduct(fftSignalCompare, fftResampledCompare) / Norm(fftSignalCompare) / Norm(fftResampledCompare); + + return similarity; } -TEST_CASE("Polyphase resampling upsample ramp", "[Interpolation]") { - constexpr int numFilters = 4; - const Signal filter = FirFilter(31, Lowpass(WINDOWED).Cutoff(1.0f / numFilters)); - Signal scratch(filter.Size()); - const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(scratch, filter, numFilters)); +TEST_CASE("Resampling spectrum invariance - upsample mild", "[Interpolation]") { + constexpr int inputRate = 7; + constexpr int outputRate = 11; + constexpr int supersamplingRate = 16; + constexpr int signalSize = 1024; + constexpr auto filterCutoff = ResamplingFilterCutoff({ inputRate, outputRate }, supersamplingRate); - const Signal signal = MakeRamp(150); - Signal output(signal.Size() - polyphase[0].Size() - 1); + for (const int filterSize : { 513, 2047 }) { + const auto signal = RandomSignal(signalSize); + const auto filter = FirFilter(filterSize, Lowpass(WINDOWED).Cutoff(filterCutoff)); + const auto polyphase = PolyphaseDecompose(filter, supersamplingRate); - Resample(output, signal, polyphase, { 7, 11 }, { filter.Size(), numFilters }); + const auto length = ResamplingLength(signalSize, filterSize, supersamplingRate, { inputRate, outputRate }, CONV_FULL); + const auto resampled = Resample(signal, polyphase, { inputRate, outputRate }, { 0, 1 }, floor(length)); + const auto similarity = ResampledSimilarity({ inputRate, outputRate }, signal, resampled); - const SignalView frontView{ output.begin(), output.Size() - 1 }; - const SignalView backView{ output.begin() + 1, output.Size() - 1 }; - const auto diff = backView - frontView; - REQUIRE(Mean(diff) == Approx(7.f / 11.f).epsilon(0.001f)); - REQUIRE(Min(diff) == Approx(7.f / 11.f).epsilon(0.02f)); - REQUIRE(Max(diff) == Approx(7.f / 11.f).epsilon(0.02f)); + INFO("filterSize=" << filterSize) + REQUIRE(similarity > 0.98f); + } } +TEST_CASE("Resampling spectrum invariance - upsample strong", "[Interpolation]") { + constexpr int inputRate = 9; + constexpr int outputRate = 210; + constexpr int supersamplingRate = 32; + constexpr int signalSize = 2048; + constexpr auto filterCutoff = ResamplingFilterCutoff({ inputRate, outputRate }, supersamplingRate); -TEST_CASE("Polyphase resampling downsample ramp mild", "[Interpolation]") { - constexpr int numFilters = 4; - const std::pair ratio = { 11, 7 }; - const float ratioReal = float(ratio.first) / float(ratio.second); - const Signal filter = FirFilter(31, Lowpass(WINDOWED).Cutoff(1.0f / ratioReal / numFilters)); + for (const int filterSize : { 1023, 4047 }) { + const auto signal = RandomSignal(signalSize); + const auto filter = FirFilter(filterSize, Lowpass(WINDOWED).Cutoff(filterCutoff)); + const auto polyphase = PolyphaseDecompose(filter, supersamplingRate); - Signal scratch(filter.Size()); - const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(scratch, filter, numFilters)); + const auto length = ResamplingLength(signalSize, filterSize, supersamplingRate, { inputRate, outputRate }, CONV_FULL); + const auto resampled = Resample(signal, polyphase, { inputRate, outputRate }, { 0, 1 }, floor(length)); + const auto similarity = ResampledSimilarity({ inputRate, outputRate }, signal, resampled); - const Signal signal = MakeRamp(150); - Signal output(signal.Size() * ratio.second / ratio.first - polyphase[0].Size() - 1); + INFO("filterSize=" << filterSize) + REQUIRE(similarity > 0.98f); + } +} + +TEST_CASE("Resampling spectrum invariance - downsample mild", "[Interpolation]") { + constexpr int inputRate = 11; + constexpr int outputRate = 7; + constexpr int supersamplingRate = 16; + constexpr int signalSize = 16384; + constexpr auto filterCutoff = ResamplingFilterCutoff({ inputRate, outputRate }, supersamplingRate); + + for (const int filterSize : { 4095, 20001 }) { + const auto signal = RandomSignal(signalSize); + const auto filter = FirFilter(filterSize, Lowpass(WINDOWED).Cutoff(filterCutoff)); + const auto polyphase = PolyphaseDecompose(filter, supersamplingRate); - Resample(output, signal, polyphase, ratio, { filter.Size(), numFilters }); + const auto length = ResamplingLength(signalSize, filterSize, supersamplingRate, { inputRate, outputRate }, CONV_FULL); + const auto resampled = Resample(signal, polyphase, { inputRate, outputRate }, { 0, 1 }, floor(length)); + const auto similarity = ResampledSimilarity({ inputRate, outputRate }, signal, resampled); - const SignalView frontView{ output.begin(), output.Size() - 1 }; - const SignalView backView{ output.begin() + 1, output.Size() - 1 }; - const auto diff = backView - frontView; - REQUIRE(Mean(diff) == Approx(ratioReal).epsilon(0.001f)); - REQUIRE(Min(diff) == Approx(ratioReal).epsilon(0.02f)); - REQUIRE(Max(diff) == Approx(ratioReal).epsilon(0.02f)); + INFO("filterSize=" << filterSize) + REQUIRE(similarity > 0.98f); + } } +TEST_CASE("Resampling spectrum invariance - downsample strong", "[Interpolation]") { + constexpr int inputRate = 210; + constexpr int outputRate = 9; + constexpr int supersamplingRate = 16; + constexpr int signalSize = 16384; + constexpr auto filterCutoff = ResamplingFilterCutoff({ inputRate, outputRate }, supersamplingRate); -TEST_CASE("Polyphase resampling downsample ramp strong", "[Interpolation]") { - constexpr int numFilters = 4; - const std::pair ratio = { 39, 7 }; - const float ratioReal = float(ratio.first) / float(ratio.second); - const Signal filter = FirFilter(31, Lowpass(WINDOWED).Cutoff(1.0f / ratioReal / numFilters)); + for (const int filterSize : { 4095, 20001 }) { + const auto signal = RandomSignal(signalSize); + const auto filter = FirFilter(filterSize, Lowpass(WINDOWED).Cutoff(filterCutoff)); + const auto polyphase = PolyphaseDecompose(filter, supersamplingRate); - Signal scratch(filter.Size()); - const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(scratch, filter, numFilters)); + const auto length = ResamplingLength(signalSize, filterSize, supersamplingRate, { inputRate, outputRate }, CONV_FULL); + const auto resampled = Resample(signal, polyphase, { inputRate, outputRate }, { 0, 1 }, floor(length)); + const auto similarity = ResampledSimilarity({ inputRate, outputRate }, signal, resampled); - const Signal signal = MakeRamp(150); - Signal output(signal.Size() * ratio.second / ratio.first - polyphase[0].Size() - 1); + INFO("filterSize=" << filterSize) + REQUIRE(similarity > 0.98f); + } +} - Resample(output, signal, polyphase, ratio, { filter.Size(), numFilters }); - const SignalView frontView{ output.begin(), output.Size() - 1 }; - const SignalView backView{ output.begin() + 1, output.Size() - 1 }; - const auto diff = backView - frontView; - REQUIRE(Mean(diff) == Approx(ratioReal).epsilon(0.001f)); - REQUIRE(Min(diff) == Approx(ratioReal).epsilon(0.02f)); - REQUIRE(Max(diff) == Approx(ratioReal).epsilon(0.02f)); +template +double FindCrossing(const SignalT& signal, double value) { + const auto it = std::adjacent_find(signal.begin(), signal.end(), [&value](auto left, auto right) { + return left <= value && value < right; + }); + if (it != signal.end()) { + const size_t firstIndex = it - signal.begin(); + const auto difference = it[1] - it[0]; + const auto t = (value - it[0]) / difference; + return double(firstIndex) + double(t); + } + return -1.0; } -TEST_CASE("Polyphase resampling shift ramp", "[Interpolation]") { - constexpr int numFilters = 2; - const Signal filter = FirFilter(63, Lowpass(WINDOWED).Cutoff(1.0f / numFilters)); - Signal scratch(filter.Size()); - const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(scratch, filter, numFilters)); +TEST_CASE("Resampling delay - upsample mild", "[Interpolation]") { + // Resample a ramp function. + // The exact crossing (i.e. f(x) = 10, x = ?) can be easily found by linear interpolation. + // The exact crossing can be used to correlate delays on the input and output signals. + + constexpr int inputRate = 7; + constexpr int outputRate = 11; + constexpr int supersamplingRate = 16; + constexpr int signalSize = 1024; + constexpr auto filterCutoff = ResamplingFilterCutoff({ inputRate, outputRate }, supersamplingRate); - const Signal signal = MakeRamp(150); - Signal output(signal.Size() - polyphase[0].Size() - 1); + for (const int filterSize : { 513, 2047 }) { + auto signal = Signal(signalSize); + std::iota(signal.begin(), signal.end(), 0.0f); + const auto filter = FirFilter(filterSize, Lowpass(WINDOWED).Cutoff(filterCutoff)); + const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(filter, supersamplingRate)); - const std::pair offset = { filter.Size() * 100 + 42 * numFilters, numFilters * 100 }; - const float offsetReal = float(offset.first) / float(offset.second) - float(filter.Size() / numFilters) / 2; - Resample(output, signal, polyphase, { 1, 1 }, offset); + const auto length = ResamplingLength(signalSize, filterSize, supersamplingRate, { inputRate, outputRate }, CONV_FULL); + const auto resampled = Resample(signal, polyphase, { inputRate, outputRate }, { 0, 1 }, floor(length)); - REQUIRE(output[0] == Approx(offsetReal).epsilon(0.02f)); + const double crossingSignal = FindCrossing(signal, 500.0); + const double crossingResampled = FindCrossing(resampled, 500.0); + const auto resamplingDelay = ResamplingDelay(filterSize, supersamplingRate, { inputRate, outputRate }); + const double crossingExpected = double(resamplingDelay) + crossingSignal * outputRate / inputRate; + + INFO("filterSize=" << filterSize) + REQUIRE(crossingExpected == Approx(crossingResampled)); + } } -TEST_CASE("Polyphase resampling returned offset", "[Interpolation]") { - constexpr int numFilters = 5; - const Signal filter = FirFilter(63, Lowpass(WINDOWED).Cutoff(1.0f / numFilters)); - Signal scratch(filter.Size()); - const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(scratch, filter, numFilters)); +TEST_CASE("Interplation continuation calculation", "[Interpolation]") { + constexpr size_t numPhases = 6; + constexpr size_t filterSize = 31; + + SECTION("Initial point") { + constexpr size_t nextOutputSample = 0; + const auto [inputIndex, startPoint] = impl::FindInterpSuspensionPoint(nextOutputSample, filterSize, numPhases); - const Signal signal = MakeRamp(150); - Signal output(17); + REQUIRE(inputIndex == 0); + REQUIRE(startPoint == 0); + } + SECTION("One off") { + constexpr size_t nextOutputSample = 2; + const auto [inputIndex, startPoint] = impl::FindInterpSuspensionPoint(nextOutputSample, filterSize, numPhases); + + REQUIRE(inputIndex == 0); + REQUIRE(startPoint == 2); + } + SECTION("Middle point") { + constexpr size_t nextOutputSample = 36; + const auto [inputIndex, startPoint] = impl::FindInterpSuspensionPoint(nextOutputSample, filterSize, numPhases); - const std::pair offset = { 173, 982 }; - const std::pair ratio = { 7743, 9235 }; - const auto last = Resample(output, signal, polyphase, ratio, offset); + REQUIRE(inputIndex == 1); + REQUIRE(startPoint == 30); + } + SECTION("Far point") { + constexpr size_t nextOutputSample = 158; + const auto [inputIndex, startPoint] = impl::FindInterpSuspensionPoint(nextOutputSample, filterSize, numPhases); - REQUIRE(last.first / last.second == 17 * ratio.first / ratio.second); - REQUIRE((last.first % last.second) / double(last.second) == Approx(0.4296632)); + REQUIRE(inputIndex == 21); + REQUIRE(startPoint == 32); + } } + + +TEST_CASE("Resampling continuation calculation", "[Interpolation]") { + constexpr size_t numPhases = 6; + constexpr size_t filterSize = 31; + constexpr Rational sampleRates = { 4, 7 }; + + SECTION("Initial point") { + constexpr Rational nextOutputSample = { 0, 1 }; + const auto [inputIndex, startPoint] = impl::FindResamplingSuspensionPoint(nextOutputSample, filterSize, numPhases, sampleRates); + + REQUIRE(inputIndex == 0); + REQUIRE(double(startPoint) == Approx(0)); + } + SECTION("One off") { + constexpr Rational nextOutputSample = { 7, 7 }; + const auto [inputIndex, startPoint] = impl::FindResamplingSuspensionPoint(nextOutputSample, filterSize, numPhases, sampleRates); + + REQUIRE(inputIndex == 0); + REQUIRE(double(startPoint) == Approx(1)); + } + SECTION("Middle point") { + constexpr Rational nextOutputSample = { 6 * 7, 4 }; + const auto [inputIndex, startPoint] = impl::FindResamplingSuspensionPoint(nextOutputSample, filterSize, numPhases, sampleRates); + + REQUIRE(inputIndex == 1); + const double expectedTotalOffset = double(nextOutputSample); + const double actualTotalOffset = double(inputIndex) / double(sampleRates) + double(startPoint); + REQUIRE(expectedTotalOffset == Approx(actualTotalOffset)); + } + SECTION("Far point") { + constexpr Rational nextOutputSample = { 156, 1 }; + const auto [inputIndex, startPoint] = impl::FindResamplingSuspensionPoint(nextOutputSample, filterSize, numPhases, sampleRates); + + REQUIRE(inputIndex == 84); + const double expectedTotalOffset = double(nextOutputSample); + const double actualTotalOffset = double(inputIndex) / double(sampleRates) + double(startPoint); + REQUIRE(expectedTotalOffset == Approx(actualTotalOffset)); + } +} + + +TEST_CASE("Interpolation continuation output", "[Interpolation]") { + constexpr size_t numPhases = 6; + constexpr size_t filterSize = 511; + constexpr float filterCutoff = float(InterpFilterCutoff(numPhases)); + + const auto filter = FirFilter(filterSize, Lowpass(LEAST_SQUARES).Cutoff(0.90f * filterCutoff, filterCutoff)); + const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(filter, numPhases)); + + // This creates a linearly increasing ramp-like function + const auto signal = LinSpace(0.0f, 100.f, 2500); + + const size_t maxLength = InterpLength(signal.Size(), filterSize, numPhases, CONV_FULL); + + auto output = Signal(maxLength, 0.0f); + + size_t chunkSize = 1; + size_t outputWritten = 0; + size_t firstInputSample = 0; + size_t startPoint{ 0 }; + while (outputWritten < output.Size() / 2) { + const auto [newFirstInputSample, newStartPoint] = Interpolate(AsView(output).SubSignal(outputWritten, chunkSize), + AsView(signal).SubSignal(firstInputSample), + polyphase, + startPoint); + + startPoint = newStartPoint; + firstInputSample += newFirstInputSample; + outputWritten += chunkSize; + chunkSize *= 2; + } + + // Find the linear part of the output + const auto first = std::find_if(output.begin(), output.end(), [](float v) { return v >= 3.0f; }); + const auto last = std::max_element(output.begin(), output.end()); + + REQUIRE(first != last); + REQUIRE(first != output.end()); + REQUIRE(last != output.end()); + REQUIRE(first - output.begin() < ptrdiff_t(output.Size()) / 30 + filterSize - 1); + REQUIRE(last - output.begin() >= ptrdiff_t(output.Size()) / 2); + + // Check if increments between adjacent elements of the output ramp are roughly equal + SignalView left{ first, last - 1 }; + SignalView right{ first + 1, last }; + REQUIRE(Max(right - left) == Approx(Min(right - left)).epsilon(0.02)); +} + + + +TEST_CASE("Resampling continuation output", "[Interpolation]") { + constexpr size_t numPhases = 6; + constexpr size_t filterSize = 511; + constexpr Rational sampleRates = { 4, 7 }; + constexpr float filterCutoff = float(ResamplingFilterCutoff(sampleRates, numPhases)); + + const auto filter = FirFilter(filterSize, Lowpass(LEAST_SQUARES).Cutoff(0.90f * filterCutoff, filterCutoff)); + const auto polyphase = PolyphaseNormalized(PolyphaseDecompose(filter, numPhases)); + + // This creates a linearly increasing ramp-like function + const auto signal = LinSpace(0.0f, 100.f, 2500); + + const Rational maxLength = ResamplingLength(signal.Size(), filterSize, numPhases, sampleRates, CONV_FULL); + + auto output = Signal(floor(maxLength), 0.0f); + + size_t chunkSize = 1; + size_t outputWritten = 0; + size_t firstInputSample = 0; + Rational startPoint{ 0 }; + while (outputWritten < output.Size() / 2) { + const auto [newFirstInputSample, newStartPoint] = Resample(AsView(output).SubSignal(outputWritten, chunkSize), + AsView(signal).SubSignal(firstInputSample), + polyphase, + sampleRates, + startPoint); + + startPoint = newStartPoint; + firstInputSample += newFirstInputSample; + outputWritten += chunkSize; + chunkSize *= 2; + } + + // Find the linear part of the output + const auto first = std::find_if(output.begin(), output.end(), [](float v) { return v >= 3.0f; }); + const auto last = std::max_element(output.begin(), output.end()); + + REQUIRE(first != last); + REQUIRE(first != output.end()); + REQUIRE(last != output.end()); + REQUIRE(first - output.begin() < ptrdiff_t(output.Size()) / 30 + ceil(ResamplingDelay(filterSize, numPhases, sampleRates))); + REQUIRE(last - output.begin() >= ptrdiff_t(output.Size()) / 2); + + // Check if increments between adjacent elements of the output ramp are roughly equal + SignalView left{ first, last - 1 }; + SignalView right{ first + 1, last }; + REQUIRE(Max(right - left) == Approx(Min(right - left)).epsilon(0.02)); +} \ No newline at end of file diff --git a/test/Filtering/Test_Polyphase.cpp b/test/Filtering/Test_Polyphase.cpp index 3c73247..ca53da7 100644 --- a/test/Filtering/Test_Polyphase.cpp +++ b/test/Filtering/Test_Polyphase.cpp @@ -9,9 +9,8 @@ using namespace dspbb; TEST_CASE("Polyphase view filter non-uniform", "[Polyphase]") { const Signal filter = { 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2 }; const std::array filterSizes = { 3, 3, 3, 2 }; - Signal scratch(filter.Size(), std::numeric_limits::quiet_NaN()); - const auto view = PolyphaseDecompose(scratch, filter, 4); - REQUIRE(view.numFilters == 4); + const auto view = PolyphaseDecompose(filter, 4); + REQUIRE(view.FilterCount() == 4); for (size_t i = 0; i < 4; ++i) { REQUIRE(std::all_of(view[i].begin(), view[i].end(), [i](float c) { return c == float(4 * i); })); REQUIRE(view[i].Size() == filterSizes[i]); @@ -21,9 +20,8 @@ TEST_CASE("Polyphase view filter non-uniform", "[Polyphase]") { TEST_CASE("Polyphase view filter uniform", "[Polyphase]") { const Signal filter = { 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3 }; const std::array filterSizes = { 3, 3, 3, 3 }; - Signal scratch(filter.Size(), std::numeric_limits::quiet_NaN()); - const auto view = PolyphaseDecompose(scratch, filter, 4); - REQUIRE(view.numFilters == 4); + const auto view = PolyphaseDecompose(filter, 4); + REQUIRE(view.FilterCount() == 4); for (size_t i = 0; i < 4; ++i) { REQUIRE(std::all_of(view[i].begin(), view[i].end(), [i](float c) { return c == float(4 * i); })); REQUIRE(view[i].Size() == filterSizes[i]); @@ -32,8 +30,7 @@ TEST_CASE("Polyphase view filter uniform", "[Polyphase]") { TEST_CASE("Polyphase normalize", "[Polyphase]") { const Signal filter = { 1, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3 }; - Signal scratch(filter.Size(), std::numeric_limits::quiet_NaN()); - const auto view = PolyphaseNormalized(PolyphaseDecompose(scratch, filter, 4)); + const auto view = PolyphaseNormalized(PolyphaseDecompose(filter, 4)); for (size_t i = 0; i < 4; ++i) { REQUIRE(Sum(view[i]) == Approx(1.0f)); } @@ -42,8 +39,7 @@ TEST_CASE("Polyphase normalize", "[Polyphase]") { TEST_CASE("Polyphase reverse", "[Polyphase]") { const Signal filter = { 0, 1, 2, 3 }; - Signal scratch(filter.Size(), std::numeric_limits::quiet_NaN()); - const auto view = PolyphaseDecompose(scratch, filter, 2); + const auto view = PolyphaseDecompose(filter, 2); REQUIRE(view[0][0] == 2 * 2); REQUIRE(view[0][1] == 2 * 0); REQUIRE(view[1][0] == 2 * 3); diff --git a/test/TestUtils.hpp b/test/TestUtils.hpp index f5854d4..613e93f 100644 --- a/test/TestUtils.hpp +++ b/test/TestUtils.hpp @@ -121,7 +121,7 @@ dspbb::Signal RandomPositiveSignal(size_t size) { template dspbb::BasicSignal RandomSignal(size_t length) { thread_local std::mt19937_64 rne(723574); - thread_local std::uniform_real_distribution rng; + thread_local std::uniform_real_distribution rng(-1, 1); dspbb::BasicSignal s(length); for (auto& v : s) { if constexpr (dspbb::is_complex_v) {