diff --git a/include/algorithms/public/TransientExtraction.hpp b/include/algorithms/public/TransientExtraction.hpp index 304b2265f..e9cff04a6 100644 --- a/include/algorithms/public/TransientExtraction.hpp +++ b/include/algorithms/public/TransientExtraction.hpp @@ -42,7 +42,10 @@ class TransientExtraction mBackwardError(asUnsigned(maxBlockSize + maxOrder), alloc), mForwardWindowedError(asUnsigned(maxBlockSize), alloc), mBackwardWindowedError(asUnsigned(maxBlockSize), alloc) - {} + { + assert(maxBlockSize >= maxOrder && "Max block size needs to be >= max order"); + assert(maxPadSize >= maxOrder && "Max pad size needs to be >= max order"); + } void init(index order, index blockSize, index padSize) { diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index b15aa6084..3c4e4f0c3 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -10,8 +10,8 @@ under the European Union’s Horizon 2020 research and innovation programme #pragma once -#include "FluidEigenMappings.hpp" #include "FFT.hpp" +#include "FluidEigenMappings.hpp" #include "Toeplitz.hpp" #include "../public/WindowFuncs.hpp" #include "../../data/FluidIndex.hpp" @@ -35,7 +35,7 @@ class ARModel public: ARModel(index maxOrder, index maxBlockSize, Allocator& alloc) - : mParameters(maxOrder, alloc), + : mOrder{maxOrder}, mParameters(maxOrder, alloc), mWindow(maxBlockSize + maxOrder, alloc), mSpectralIn(nextPower2(2 * maxBlockSize), alloc), mSpectralOut(nextPower2(2 * maxBlockSize), alloc), @@ -43,15 +43,21 @@ class ARModel mIFFT(nextPower2(2 * maxBlockSize), alloc) {} - void init(index order) { mParameters.head(order).setZero(); } + void init(index order) + { + assert(order <= mParameters.size() && + "Order > maxOrder. Check your constructor arguments!"); + mOrder = order; + mParameters.head(order).setZero(); + } const double* getParameters() const { return mParameters.data(); } double variance() const { return mVariance; } - index order() const { return mParameters.size(); } + index order() const { return mOrder; } void estimate(FluidTensorView input, index nIterations, - double robustFactor, Allocator& alloc) + double robustFactor, Allocator& alloc) { if (nIterations > 0) robustEstimate(input, nIterations, robustFactor, alloc); @@ -63,7 +69,7 @@ class ARModel { double prediction; modelPredict(input, FluidTensorView(&prediction, 0, 1), - std::negate{}, Predict{}); + std::negate{}, Predict{}); return prediction; } @@ -71,7 +77,7 @@ class ARModel { double prediction; modelPredict(input, FluidTensorView(&prediction, 0, 1), - Identity{}, Predict{}); + Identity{}, Predict{}); return prediction; } @@ -89,14 +95,14 @@ class ARModel return error; } - void forwardErrorArray( - FluidTensorView input, FluidTensorView errors) + void forwardErrorArray(FluidTensorView input, + FluidTensorView errors) { modelPredict(input, errors, std::negate{}, Error{}); } - void backwardErrorArray( - FluidTensorView input, FluidTensorView errors) + void backwardErrorArray(FluidTensorView input, + FluidTensorView errors) { modelPredict(input, errors, Identity{}, Error{}); } @@ -124,25 +130,25 @@ class ARModel template void modelPredict(FluidTensorView input, - FluidTensorView output, Indexer fIdx, OutputFn fOut) + FluidTensorView output, Indexer fIdx, + OutputFn fOut) { index numPredictions = output.size(); assert(fIdx(1) >= -input.descriptor().start && "array bounds error in AR model prediction: input too short"); assert(fIdx(1) < input.size() && "array bounds error in AR model prediction: input too short"); - assert(fIdx(mParameters.size()) >= -input.descriptor().start && + assert(fIdx(mOrder) >= -input.descriptor().start && "array bounds error in AR model prediction: input too short"); - assert(fIdx(mParameters.size()) < input.size() && + assert(fIdx(mOrder) < input.size() && "array bounds error in AR model prediction: input too short"); assert(((numPredictions - 1) + fIdx(1)) >= -input.descriptor().start && "array bounds error in AR model prediction: input too short"); assert(((numPredictions - 1) + fIdx(1)) < input.size() && "array bounds error in AR model prediction: input too short"); - assert(((numPredictions - 1) + fIdx(mParameters.size())) >= - -input.descriptor().start && + assert(((numPredictions - 1) + fIdx(mOrder)) >= -input.descriptor().start && "array bounds error in AR model prediction: input too short"); - assert(((numPredictions - 1) + fIdx(mParameters.size())) < input.size() && + assert(((numPredictions - 1) + fIdx(mOrder)) < input.size() && "array bounds error in AR model prediction: input too short"); const double* input_ptr = input.data(); @@ -150,19 +156,18 @@ class ARModel for (index p = 0; p < numPredictions; p++) { double estimate = 0; - for (index i = 0; i < mParameters.size(); i++) + for (index i = 0; i < mOrder; i++) estimate += mParameters(i) * (input_ptr + p)[fIdx(i + 1)]; output[p] = fOut(input_ptr[p], estimate); } } void directEstimate(FluidTensorView input, - bool updateVariance, Allocator& alloc) + bool updateVariance, Allocator& alloc) { index size = input.size(); // copy input to a 32 byte aligned block (otherwise risk segfaults on Linux) - // FluidEigenMap frame = ScopedEigenMap frame(size, alloc); frame = _impl::asEigen(input); @@ -170,33 +175,25 @@ class ARModel { if (mWindow.size() != size) { - mWindow.setZero(); // = ArrayXd::Zero(size); - WindowFuncs::map()[WindowFuncs::WindowTypes::kHann]( - size, mWindow.head(size)); + mWindow.setZero(); + WindowFuncs::map()[WindowFuncs::WindowTypes::kHann](size, + mWindow.head(size)); } frame.array() *= mWindow.head(size); } - // VectorXd autocorrelation(size); + ScopedEigenMap autocorrelation(size, alloc); autocorrelate(frame, autocorrelation); - // algorithm::autocorrelateReal(autocorrelation.data(), frame.data(), - - // asUnsigned(size)); - // Resize to the desired order (only keep coefficients for up to the order - // we need) - double pN = mParameters.size() < size ? autocorrelation(mParameters.size()) - : autocorrelation(0); - // autocorrelation.conservativeResize(mParameters.size()); + double pN = mOrder < size ? autocorrelation(mOrder) : autocorrelation(0); // Form a toeplitz matrix - // MatrixXd mat = toeplitz(autocorrelation); - ScopedEigenMap mat(mParameters.size(), mParameters.size(), alloc); - mat = MatrixXd::NullaryExpr(mParameters.size(), mParameters.size(), - [&autocorrelation, n = mParameters.size()]( - Eigen::Index row, Eigen::Index col) { + ScopedEigenMap mat(mOrder, mOrder, alloc); + mat = MatrixXd::NullaryExpr( + mOrder, mOrder, + [&autocorrelation, n = mOrder](Eigen::Index row, Eigen::Index col) { index idx = row - col; if (idx < 0) idx += n; return autocorrelation(idx); @@ -205,54 +202,53 @@ class ARModel // Yule Walker autocorrelation(0) = pN; std::rotate(autocorrelation.data(), autocorrelation.data() + 1, - autocorrelation.data() + mParameters.size()); - mParameters = mat.llt().solve(autocorrelation.head(mParameters.size())); + autocorrelation.data() + mOrder); + mParameters.head(mOrder) = mat.llt().solve(autocorrelation.head(mOrder)); if (updateVariance) { // Calculate variance double variance = mat(0, 0); - for (index i = 0; i < mParameters.size() - 1; i++) + for (index i = 0; i < mOrder - 1; i++) variance -= mParameters(i) * mat(0, i + 1); - setVariance( - (variance - (mParameters(mParameters.size() - 1) * pN)) / size); + setVariance((variance - (mParameters(mOrder - 1) * pN)) / size); } } void robustEstimate(FluidTensorView input, index nIterations, - double robustFactor, Allocator& alloc) + double robustFactor, Allocator& alloc) { - FluidTensor estimates(input.size() + mParameters.size(), alloc); + FluidTensor estimates(input.size() + mOrder, alloc); // Calculate an initial estimate of parameters directEstimate(input, true, alloc); - assert(input.descriptor().start >= mParameters.size() && + assert(input.descriptor().start >= mOrder && "too little offset into input data"); // Initialise Estimates - for (index i = 0; i < input.size() + mParameters.size(); i++) - estimates[i] = input.data()[i - mParameters.size()]; + for (index i = 0; i < input.size() + mOrder; i++) + estimates[i] = input.data()[i - mOrder]; // Variance - robustVariance(estimates(Slice(mParameters.size())), input, robustFactor); + robustVariance(estimates(Slice(mOrder)), input, robustFactor); // Iterate for (index iterations = nIterations; iterations--;) - robustIteration( - estimates(Slice(mParameters.size())), input, robustFactor, alloc); + robustIteration(estimates(Slice(mOrder)), input, robustFactor, alloc); } double robustResidual(double input, double prediction, double cs) { - assert(cs > 0); + assert(cs > 0); return cs * psiFunction((input - prediction) / cs); } - void robustVariance(FluidTensorView estimates, - FluidTensorView input, double robustFactor) + void robustVariance(FluidTensorView estimates, + FluidTensorView input, + double robustFactor) { double residualSqSum = 0.0; @@ -261,23 +257,24 @@ class ARModel { const double residual = robustResidual(input[i], fowardPrediction(estimates(Slice(i))), - robustFactor * sqrt(mVariance)); + robustFactor * sqrt(mVariance)); residualSqSum += residual * residual; } setVariance(residualSqSum / input.size()); } - void robustIteration(FluidTensorView estimates, - FluidTensorView input, double robustFactor, - Allocator& alloc) + void robustIteration(FluidTensorView estimates, + FluidTensorView input, + double robustFactor, Allocator& alloc) { // Iterate to find new filtered input for (index i = 0; i < input.size(); i++) { const double prediction = fowardPrediction(estimates(Slice(i))); - estimates[i] = prediction + robustResidual(input[i], prediction, - robustFactor * sqrt(mVariance)); + estimates[i] = + prediction + + robustResidual(input[i], prediction, robustFactor * sqrt(mVariance)); } // New parameters @@ -289,7 +286,7 @@ class ARModel { variance = std::max(mMinVariance, variance); mVariance = variance; - assert(mVariance >= 0); + assert(mVariance >= 0); } // Huber PSI function @@ -299,8 +296,8 @@ class ARModel } template - void autocorrelate( - const Eigen::DenseBase& in, Eigen::DenseBase& out) + void autocorrelate(const Eigen::DenseBase& in, + Eigen::DenseBase& out) { mSpectralIn.head(in.size()) = in; @@ -324,6 +321,7 @@ class ARModel return static_cast(std::pow(2, std::ceil(std::log2(n)))); } + index mOrder; ScopedEigenMap mParameters; double mVariance{0.0}; ScopedEigenMap mWindow; diff --git a/include/clients/rt/TransientClient.hpp b/include/clients/rt/TransientClient.hpp index ccf5ce76e..56887efaa 100644 --- a/include/clients/rt/TransientClient.hpp +++ b/include/clients/rt/TransientClient.hpp @@ -37,11 +37,16 @@ enum TransientParamIndex { kDebounce }; +constexpr index transientMaxOrder = 200; +constexpr index transientMaxBlockSize = 4096; +constexpr index transientMaxPadding = 1024; + constexpr auto TransientParams = defineParameters( LongParam("order", "Order", 20, Min(10), LowerLimit(), - UpperLimit(), Max(200)), - LongParam("blockSize", "Block Size", 256, Min(100), LowerLimit(), Max(4096)), - LongParam("padSize", "Padding", 128, Min(0)), + UpperLimit(), Max(transientMaxOrder)), + LongParam("blockSize", "Block Size", 256, Min(100), LowerLimit(), + Max(transientMaxBlockSize)), + LongParam("padSize", "Padding", 128, Min(0), Max(transientMaxPadding)), FloatParam("skew", "Skew", 0, Min(-10), Max(10)), FloatParam("threshFwd", "Forward Threshold", 2, Min(0)), FloatParam("threshBack", "Backward Threshold", 1.1, Min(0)), @@ -67,10 +72,13 @@ class TransientClient : public FluidBaseClient, public AudioIn, public AudioOut static constexpr auto& getParameterDescriptors() { return TransientParams; } TransientClient(ParamSetViewType& p, FluidContext const& c) - : mParams(p), maxWindowIn{2 * get() + get()}, - maxWindowOut{get()}, - mBufferedProcess{maxWindowIn, maxWindowOut, 1, 2, c.hostVectorSize(), c.allocator()}, - mExtractor{get(), get(), get(), c.allocator()} + : mParams(p), + maxWindowIn{2 * transientMaxBlockSize + transientMaxPadding}, + maxWindowOut{transientMaxBlockSize}, + mExtractor{transientMaxOrder, transientMaxBlockSize, + transientMaxPadding, c.allocator()}, + mBufferedProcess{maxWindowIn, maxWindowOut, 1, 2, + c.hostVectorSize(), c.allocator()} { audioChannelsIn(1); audioChannelsOut(2); @@ -88,26 +96,20 @@ class TransientClient : public FluidBaseClient, public AudioIn, public AudioOut index blockSize = get(); index padding = get(); index hostVecSize = input[0].size(); -// index maxWinIn = 2 * blockSize + padding; -// index maxWinOut = blockSize - order; if (mTrackValues.changed(order, blockSize, padding, hostVecSize) || !mExtractor.initialized()) { mExtractor.init(order, blockSize, padding); -// mBufferedProcess.hostSize(hostVecSize); -// mBufferedProcess.maxSize(maxWinIn, maxWinOut, -// FluidBaseClient::audioChannelsIn(), -// FluidBaseClient::audioChannelsOut()); } double skew = pow(2, get()); double threshFwd = get(); - double thresBack = get(); + double threshBack = get(); index halfWindow = static_cast(round(get() / 2)); index debounce = get(); - mExtractor.setDetectionParameters(skew, threshFwd, thresBack, halfWindow, + mExtractor.setDetectionParameters(skew, threshFwd, threshBack, halfWindow, debounce); RealMatrix in(1, hostVecSize, c.allocator()); diff --git a/include/clients/rt/TransientSliceClient.hpp b/include/clients/rt/TransientSliceClient.hpp index 9469b550d..e343cc67c 100644 --- a/include/clients/rt/TransientSliceClient.hpp +++ b/include/clients/rt/TransientSliceClient.hpp @@ -36,11 +36,16 @@ enum TransientParamIndex { kMinSeg }; +constexpr index transientMaxOrder = 200; +constexpr index transientMaxBlockSize = 4096; +constexpr index transientMaxPadding = 1024; + constexpr auto TransientSliceParams = defineParameters( LongParam("order", "Order", 20, Min(10), LowerLimit(), - UpperLimit(), Max(200)), - LongParam("blockSize", "Block Size", 256, Min(100), LowerLimit(), Max(4096)), - LongParam("padSize", "Padding", 128, Min(0)), + UpperLimit(), Max(transientMaxOrder)), + LongParam("blockSize", "Block Size", 256, Min(100), LowerLimit(), + Max(transientMaxBlockSize)), + LongParam("padSize", "Padding", 128, Min(0), Max(transientMaxPadding)), FloatParam("skew", "Skew", 0, Min(-10), Max(10)), FloatParam("threshFwd", "Forward Threshold", 2, Min(0)), FloatParam("threshBack", "Backward Threshold", 1.1, Min(0)), @@ -72,10 +77,13 @@ class TransientSliceClient : public FluidBaseClient, } TransientSliceClient(ParamSetViewType& p, FluidContext& c) - : mParams{p},mMaxWindowIn{2 * get() + get()}, - mMaxWindowOut{get()}, - mExtractor{get(), get(), get(), c.allocator()}, - mBufferedProcess{mMaxWindowIn, mMaxWindowOut, 1, 1, c.hostVectorSize(), c.allocator()} + : mParams{p}, + mMaxWindowIn{2 * transientMaxBlockSize + transientMaxPadding}, + mMaxWindowOut{transientMaxBlockSize}, + mExtractor{transientMaxOrder, transientMaxBlockSize, + transientMaxPadding, c.allocator()}, + mBufferedProcess{mMaxWindowIn, mMaxWindowOut, 1, 1, + c.hostVectorSize(), c.allocator()} { audioChannelsIn(1); audioChannelsOut(1);