diff --git a/CMakeLists.txt b/CMakeLists.txt index 40e667098..8f1c2544f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -12,6 +12,9 @@ set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) +if (CMAKE_VERSION VERSION_GREATER_EQUAL "3.24.0") + cmake_policy(SET CMP0135 NEW) +endif() list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/script") @@ -53,7 +56,7 @@ FetchContent_Declare( HISSTools GIT_REPOSITORY https://github.com/AlexHarker/HISSTools_Library GIT_PROGRESS TRUE - GIT_TAG 5dd8530 + GIT_TAG 559b5c8 ) FetchContent_Declare( @@ -67,10 +70,11 @@ FetchContent_Declare( FetchContent_Declare( Spectra + GIT_SHALLOW TRUE GIT_REPOSITORY https://github.com/yixuan/spectra GIT_PROGRESS TRUE GIT_BRANCH "master" - GIT_TAG "v0.9.0" + GIT_TAG "v1.0.1" ) FetchContent_Declare( @@ -80,16 +84,25 @@ FetchContent_Declare( GIT_PROGRESS TRUE ) +#see https://json.nlohmann.me/integration/cmake/#fetchcontent +FetchContent_Declare(json URL https://github.com/nlohmann/json/releases/download/v3.11.2/json.tar.xz) +set(JSON_SystemInclude ON CACHE BOOL "") + FetchContent_Declare( - json - GIT_SHALLOW TRUE - # GIT_REPOSITORY https://github.com/nlohmann/json.git - # advice on nlohmann repo is to use this mirror unless we really want ~150 meg of test data as well as headers: - GIT_REPOSITORY https://github.com/ArthurSonzogni/nlohmann_json_cmake_fetchcontent + Memory + GIT_SHALLOW TRUE + GIT_REPOSITORY https://github.com/foonathan/memory.git GIT_PROGRESS TRUE - GIT_TAG v3.10.5 + GIT_TAG main +) + +FetchContent_Declare( + fmt + GIT_SHALLOW TRUE + GIT_REPOSITORY https://github.com/fmtlib/fmt + GIT_PROGRESS TRUE + GIT_TAG master ) -set(JSON_SystemInclude ON CACHE BOOL "") if(HISS_PATH) #if hiss path is set, this will stop it downloading get_filename_component(FETCHCONTENT_SOURCE_DIR_HISSTOOLS ${HISS_PATH} ABSOLUTE) @@ -138,6 +151,18 @@ endif() add_subdirectory(${json_SOURCE_DIR} ${json_BINARY_DIR} EXCLUDE_FROM_ALL) +set(FOONATHAN_MEMORY_BUILD_TOOLS OFF CACHE BOOL "" FORCE) +set(FOONATHAN_MEMORY_BUILD_EXAMPLES OFF CACHE BOOL "" FORCE) +set(FOONATHAN_MEMORY_BUILD_TESTS OFF CACHE BOOL "" FORCE) + +FetchContent_GetProperties(Memory) +if(NOT memory_POPULATED) + FetchContent_Populate(Memory) +endif() +add_subdirectory(${memory_SOURCE_DIR} ${memory_BINARY_DIR} EXCLUDE_FROM_ALL) + +FetchContent_MakeAvailable(fmt) + # HISSTools FFT target add_library( HISSTools_FFT STATIC "${hisstools_SOURCE_DIR}/HISSTools_FFT/HISSTools_FFT.cpp" @@ -183,6 +208,8 @@ target_include_directories( "${eigen_SOURCE_DIR}" "${spectra_SOURCE_DIR}/include" "${hisstools_SOURCE_DIR}" + "${memory_SOURCE_DIR}/include/foonathan" + "${fmt_SOURCE_DIR}/include" ) target_link_libraries( FLUID_DECOMPOSITION INTERFACE @@ -190,6 +217,8 @@ target_link_libraries( flucoma_VERSION_LIB tl::optional nlohmann_json::nlohmann_json + foonathan_memory + fmt::fmt ) target_sources( FLUID_DECOMPOSITION INTERFACE ${HEADERS} @@ -224,9 +253,14 @@ if(DEFINED FLUID_ARCH) endif() #Examples + +option(BUILD_EXAMPLES "Build C++ example code (off by default)" OFF) + +if(BUILD_EXAMPLES) add_subdirectory( "${CMAKE_CURRENT_SOURCE_DIR}/examples" ) +endif() enable_testing() diff --git a/FlucomaClients.cmake b/FlucomaClients.cmake index 3aaede947..94ccc38be 100644 --- a/FlucomaClients.cmake +++ b/FlucomaClients.cmake @@ -110,6 +110,7 @@ add_client(BufSTFT clients/nrt/BufSTFTClient.hpp CLASS NRTThreadedBufferSTFTClie add_client(BufScale clients/nrt/BufScaleClient.hpp CLASS NRTThreadedBufferScaleClient ) add_client(BufSelect clients/nrt/BufSelectClient.hpp CLASS NRTThreadingSelectClient ) add_client(BufSelectEvery clients/nrt/BufSelectEveryClient.hpp CLASS NRTThreadingSelectEveryClient ) +add_client(BufSineFeature clients/rt/SineFeatureClient.hpp CLASS NRTThreadedSineFeatureClient ) add_client(BufSines clients/rt/SinesClient.hpp CLASS NRTThreadedSinesClient ) add_client(BufSpectralShape clients/rt/SpectralShapeClient.hpp CLASS NRTThreadedSpectralShapeClient ) add_client(BufStats clients/nrt/BufStatsClient.hpp CLASS NRTThreadedBufferStatsClient ) @@ -132,6 +133,7 @@ add_client(OnsetFeature clients/rt/OnsetFeatureClient.hpp CLASS RTOnsetFeatureCl add_client(OnsetSlice clients/rt/OnsetSliceClient.hpp CLASS RTOnsetSliceClient ) add_client(Pitch clients/rt/PitchClient.hpp CLASS RTPitchClient ) add_client(STFTPass clients/rt/BaseSTFTClient.hpp CLASS RTSTFTPassClient NOINSTALL) +add_client(SineFeature clients/rt/SineFeatureClient.hpp CLASS RTSineFeatureClient ) add_client(Sines clients/rt/SinesClient.hpp CLASS RTSinesClient ) add_client(SpectralShape clients/rt/SpectralShapeClient.hpp CLASS RTSpectralShapeClient ) add_kr_in_client(Stats clients/rt/RunningStatsClient.hpp CLASS RunningStatsClient ) diff --git a/FlucomaVersion.cmake b/FlucomaVersion.cmake index fd6a1866f..b6437a6fe 100644 --- a/FlucomaVersion.cmake +++ b/FlucomaVersion.cmake @@ -14,7 +14,7 @@ find_package(Git REQUIRED) set(flucoma_VERSION_MAJOR 1) set(flucoma_VERSION_MINOR 0) -set(flucoma_VERSION_PATCH 4) +set(flucoma_VERSION_PATCH 5) set(flucoma_VERSION_SUFFIX "") function(make_flucoma_version_string output_variable) diff --git a/include/algorithms/public/AudioTransport.hpp b/include/algorithms/public/AudioTransport.hpp index a86563ce6..b8cc082b8 100644 --- a/include/algorithms/public/AudioTransport.hpp +++ b/include/algorithms/public/AudioTransport.hpp @@ -19,6 +19,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/FFT.hpp" #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -41,77 +42,88 @@ class AudioTransport using ArrayXcd = Eigen::ArrayXcd; template using Ref = Eigen::Ref; - using TransportMatrix = std::vector>; + using TransportMatrix = rt::vector>; template - using vector = std::vector; + using vector = rt::vector; public: - AudioTransport(index maxFFTSize) + AudioTransport(index maxFFTSize, Allocator& alloc) : mWindowSize(maxFFTSize), mFFTSize(maxFFTSize), mBins(maxFFTSize / 2 + 1), mFFT(maxFFTSize), mSTFT(maxFFTSize, maxFFTSize, maxFFTSize / 2), mISTFT(maxFFTSize, maxFFTSize, maxFFTSize / 2), - mReassignSTFT(maxFFTSize, maxFFTSize, maxFFTSize / 2) + mReassignSTFT(maxFFTSize, maxFFTSize, maxFFTSize / 2, + static_cast(WindowFuncs::WindowTypes::kHannD)), + mBinFreqs(maxFFTSize / 2 + 1, alloc), mWindow(maxFFTSize, alloc), + mWindowSquared(maxFFTSize, alloc), mPhase(maxFFTSize / 2 + 1, alloc), + mPhaseDiff(maxFFTSize / 2 + 1, alloc), + mChanged(maxFFTSize / 2 + 1, alloc) {} void init(index windowSize, index fftSize, index hopSize) { mWindowSize = windowSize; - mWindow = ArrayXd::Zero(mWindowSize); - WindowFuncs::map()[WindowFuncs::WindowTypes::kHann](mWindowSize, mWindow); + mWindow.head(mWindowSize).setZero(); + WindowFuncs::map()[WindowFuncs::WindowTypes::kHann]( + mWindowSize, mWindow.head(mWindowSize)); mWindowSquared = mWindow * mWindow; mFFTSize = fftSize; mHopSize = hopSize; mBins = fftSize / 2 + 1; - mPhase = ArrayXd::Zero(mBins); - mChanged = ArrayXi::Zero(mBins); - mBinFreqs = ArrayXd::LinSpaced(mBins, 0, mBins - 1) * (2 * pi) / mFFTSize; - mPhaseDiff = mBinFreqs * mHopSize; - mSTFT = STFT(windowSize, fftSize, hopSize); - mISTFT = ISTFT(windowSize, fftSize, hopSize); - mReassignSTFT = STFT(windowSize, fftSize, hopSize, - static_cast(WindowFuncs::WindowTypes::kHannD)); + mPhase.setZero(); + mChanged.setZero(); + mBinFreqs.head(mBins) = + ArrayXd::LinSpaced(mBins, 0, mBins - 1) * (2 * pi) / mFFTSize; + mPhaseDiff.head(mBins) = mBinFreqs * mHopSize; + mSTFT.resize(windowSize, fftSize, hopSize); + mISTFT.resize(windowSize, fftSize, hopSize); + mReassignSTFT.resize(windowSize, fftSize, hopSize); + mFFT.resize(fftSize); mInitialized = true; } bool initialized() const { return mInitialized; } void processFrame(RealVectorView in1, RealVectorView in2, double weight, - RealMatrixView out) + RealMatrixView out, Allocator& alloc) { using namespace _impl; using namespace Eigen; assert(mInitialized); - ArrayXd frame1 = asEigen(in1); - ArrayXd frame2 = asEigen(in2); - ArrayXcd spectrum1(mBins); - ArrayXcd spectrum1Dh(mBins); - ArrayXcd spectrum2(mBins); - ArrayXcd spectrum2Dh(mBins); - ArrayXd output(frame1.size()); + ScopedEigenMap frame1(in1.size(), alloc); + frame1 = asEigen(in1); + ScopedEigenMap frame2(in2.size(), alloc); + frame2 = asEigen(in2); + ScopedEigenMap spectrum1(mBins, alloc); + ScopedEigenMap spectrum1Dh(mBins, alloc); + ScopedEigenMap spectrum2(mBins, alloc); + ScopedEigenMap spectrum2Dh(mBins, alloc); + ScopedEigenMap output(frame1.size(), alloc); mSTFT.processFrame(frame1, spectrum1); mReassignSTFT.processFrame(frame1, spectrum1Dh); mSTFT.processFrame(frame2, spectrum2); mReassignSTFT.processFrame(frame2, spectrum2Dh); - ArrayXcd result = - interpolate(spectrum1, spectrum1Dh, spectrum2, spectrum2Dh, weight); + ScopedEigenMap result = interpolate( + spectrum1, spectrum1Dh, spectrum2, spectrum2Dh, weight, alloc); mISTFT.processFrame(result, output); - out.row(0) <<= asFluid(output); - out.row(1) <<= asFluid(mWindowSquared); + _impl::asEigen(out.row(0)) = output; + _impl::asEigen(out.row(1)) = mWindowSquared.head(mWindowSize); } vector segmentSpectrum(const Ref mag, - const Ref reasignedFreq) + const Ref reasignedFreq, + Allocator& alloc) { - vector masses; - double totalMass = mag.sum() + epsilon; - ArrayXi sign = (reasignedFreq > mBinFreqs).cast(); + vector masses(alloc); + double totalMass = mag.sum() + epsilon; + ScopedEigenMap sign(reasignedFreq.size(), alloc); + sign = (reasignedFreq > mBinFreqs.head(reasignedFreq.size())).cast(); mChanged.setZero(); mChanged.segment(1, mBins - 1) = sign.segment(1, mBins - 1) - sign.segment(0, mBins - 1); SpetralMass currentMass{0, 0, 0, 0}; - for (index i = 1; i < mChanged.size(); i++) + for (index i = 1; i < mBins; i++) { if (mChanged(i) == -1) { @@ -137,10 +149,11 @@ class AudioTransport return masses; } - TransportMatrix computeTransportMatrix(std::vector m1, - std::vector m2) + TransportMatrix computeTransportMatrix(rt::vector& m1, + rt::vector& m2, + Allocator& alloc) { - TransportMatrix matrix; + TransportMatrix matrix(alloc); index index1 = 0, index2 = 0; double mass1 = m1[0].mass; double mass2 = m2[0].mass; @@ -186,41 +199,52 @@ class AudioTransport } } - ArrayXcd interpolate(Ref in1, Ref in1Dh, - Ref in2, Ref in2Dh, - double interpolation) + ScopedEigenMap interpolate(Ref in1, Ref in1Dh, + Ref in2, Ref in2Dh, + double interpolation, Allocator& alloc) { - ArrayXd mag1 = in1.abs().real(); - ArrayXd mag2 = in2.abs().real(); - ArrayXcd result = ArrayXcd::Zero(mBins); - double mag1Sum = mag1.sum(); - double mag2Sum = mag2.sum(); + ScopedEigenMap mag1(in1.size(), alloc); + mag1 = in1.abs().real(); + ScopedEigenMap mag2(in2.size(), alloc); + mag2 = in2.abs().real(); + ScopedEigenMap result(mBins, alloc); + result.setZero(); + double mag1Sum = mag1.sum(); + double mag2Sum = mag2.sum(); if (mag1Sum <= 0 && mag2Sum <= 0) { return result; } else if (mag1Sum > 0 && mag2Sum <= 0) { - return in1; + result = in1; + return result; } else if (mag1Sum <= 0 && mag2Sum > 0) { - return in2; + result = in2; + return result; } - ArrayXd phase1 = in1.arg().real(); - ArrayXd phase2 = in2.arg().real(); - ArrayXd reasignedW1 = mBinFreqs - (in1Dh / in1).imag(); - ArrayXd reasignedW2 = mBinFreqs - (in2Dh / in2).imag(); - ArrayXd newAmplitudes = ArrayXd::Zero(mBins); - ArrayXd newPhases = ArrayXd::Zero(mBins); - std::vector s1 = segmentSpectrum(mag1, reasignedW1); - std::vector s2 = segmentSpectrum(mag2, reasignedW2); + ScopedEigenMap phase1(in1.size(), alloc); + phase1 = in1.arg().real(); + ScopedEigenMap phase2(in2.size(), alloc); + phase2 = in2.arg().real(); + ScopedEigenMap reasignedW1(in1.size(), alloc); + reasignedW1 = mBinFreqs.head(in1.size()) - (in1Dh / in1).imag(); + ScopedEigenMap reasignedW2(in2.size(), alloc); + reasignedW2 = mBinFreqs.head(in2.size()) - (in2Dh / in2).imag(); + ScopedEigenMap newAmplitudes(mBins, alloc); + newAmplitudes.setZero(); + ScopedEigenMap newPhases(mBins, alloc); + newPhases.setZero(); + rt::vector s1 = segmentSpectrum(mag1, reasignedW1, alloc); + rt::vector s2 = segmentSpectrum(mag2, reasignedW2, alloc); if (s1.size() == 0 || s2.size() == 0) { return result; } - TransportMatrix matrix = computeTransportMatrix(s1, s2); + TransportMatrix matrix = computeTransportMatrix(s1, s2, alloc); for (auto t : matrix) { SpetralMass m1 = s1[asUnsigned(std::get<0>(t))]; SpetralMass m2 = s2[asUnsigned(std::get<1>(t))]; index interpolatedBin = std::lrint((1 - interpolation) * m1.centerBin + - interpolation * m2.centerBin); + interpolation * m2.centerBin); double interpolationFactor = interpolation; if (m1.centerBin != m2.centerBin) { @@ -243,21 +267,21 @@ class AudioTransport return result; } - index mWindowSize{1024}; - index mHopSize{512}; - ArrayXd mBinFreqs; - ArrayXd mWindow; - ArrayXd mWindowSquared; - index mFFTSize{1024}; - index mBins{513}; - FFT mFFT; - bool mInitialized{false}; - ArrayXd mPhase; - ArrayXd mPhaseDiff; - ArrayXi mChanged; - STFT mSTFT; - ISTFT mISTFT; - STFT mReassignSTFT; + index mWindowSize{1024}; + index mHopSize{512}; + index mFFTSize{1024}; + index mBins{513}; + FFT mFFT; + STFT mSTFT; + ISTFT mISTFT; + STFT mReassignSTFT; + bool mInitialized{false}; + ScopedEigenMap mBinFreqs; + ScopedEigenMap mWindow; + ScopedEigenMap mWindowSquared; + ScopedEigenMap mPhase; + ScopedEigenMap mPhaseDiff; + ScopedEigenMap mChanged; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/CepstrumF0.hpp b/include/algorithms/public/CepstrumF0.hpp index 4b0df471f..b8504ef9b 100644 --- a/include/algorithms/public/CepstrumF0.hpp +++ b/include/algorithms/public/CepstrumF0.hpp @@ -16,6 +16,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/PeakDetection.hpp" #include "../../data/FluidIndex.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include #include @@ -28,31 +29,33 @@ class CepstrumF0 public: using ArrayXd = Eigen::ArrayXd; - CepstrumF0(index maxSize) : mCepstrumStorage(maxSize) {} + CepstrumF0(index maxSize, Allocator& alloc) + : mDCT{maxSize, maxSize, alloc}, mCepstrum(maxSize, alloc) {} - void init(index size) + void init(index size, Allocator& alloc) { // avoid allocation of maxSize^2 at constructor - mDCT = DCT(size, size); - mDCT.init(size, size); - - mCepstrum = mCepstrumStorage.segment(0, size); +// mDCT = DCT(size, size); + mDCT.init(size, size, alloc); +// mCepstrum = mCepstrumStorage.segment(0, size); mCepstrum.setZero(); } void processFrame(const RealVectorView& input, RealVectorView output, - double minFreq, double maxFreq, double sampleRate) + double minFreq, double maxFreq, double sampleRate, Allocator& alloc) { using namespace Eigen; using namespace std; PeakDetection pd; - ArrayXd mag = _impl::asEigen(input); - ArrayXd logMag = mag.max(epsilon).log(); +// ArrayXd mag = _impl::asEigen(input); + ScopedEigenMap logMag(input.size(), alloc); + logMag = _impl::asEigen(input).max(epsilon).log(); + double pitch = 0; double confidence = 0; - index minBin = min(lrint(sampleRate / maxFreq), mag.size()); - index maxBin = min(lrint(sampleRate / minFreq), mag.size()); + index minBin = min(lrint(sampleRate / maxFreq), logMag.size()); + index maxBin = min(lrint(sampleRate / minFreq), logMag.size()); mDCT.processFrame(logMag, mCepstrum); @@ -60,7 +63,7 @@ class CepstrumF0 { auto seg = mCepstrum.segment(minBin, maxBin - minBin); auto vec = pd.process(mCepstrum.segment(minBin, maxBin - minBin), 1, - seg.minCoeff()); + seg.minCoeff(), true, true, alloc); if (vec.size() > 0) { pitch = sampleRate / (vec[0].first + minBin); @@ -72,9 +75,8 @@ class CepstrumF0 } private: - DCT mDCT{0, 0}; - ArrayXd mCepstrumStorage; - ArrayXd mCepstrum; + DCT mDCT; + ScopedEigenMap mCepstrum; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/ChromaFilterBank.hpp b/include/algorithms/public/ChromaFilterBank.hpp index 34af5abd8..7bac353a7 100644 --- a/include/algorithms/public/ChromaFilterBank.hpp +++ b/include/algorithms/public/ChromaFilterBank.hpp @@ -16,6 +16,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/AlgorithmUtils.hpp" #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include #include #include @@ -26,42 +27,53 @@ namespace algorithm { class ChromaFilterBank { public: - ChromaFilterBank(index maxBins, index maxFFT) - : mFiltersStorage(maxBins, maxFFT / 2 + 1) + ChromaFilterBank(index maxBins, index maxFFT, Allocator& alloc) + : mMaxChroma(maxBins), mMaxFFT(maxFFT), mFilters(maxBins, maxFFT/2 + 1, alloc) {} - void init(index nChroma, index nBins, double ref, double sampleRate) + void init(index nChroma, index nBins, double ref, double sampleRate, Allocator& alloc) { using namespace Eigen; - const double log2E = 1.44269504088896340736; + constexpr double log2E = 1.44269504088896340736; index fftSize = 2 * (nBins - 1); - ArrayXd freqs = ArrayXd::LinSpaced(fftSize, 0, sampleRate); + assert(fftSize <= mMaxFFT); + assert(nChroma <= mMaxChroma); + + ScopedEigenMap freqs(fftSize, alloc); + freqs = ArrayXd::LinSpaced(fftSize, 0, sampleRate); freqs = nChroma * (freqs / (ref / 16)).log() * log2E; freqs[0] = freqs[1] - 1.5 * nChroma; - ArrayXd widths = ArrayXd::Ones(fftSize); + ScopedEigenMap widths(fftSize, alloc); + widths = ArrayXd::Ones(fftSize); widths.segment(0, fftSize - 1) = freqs.segment(1, fftSize - 1) - freqs.segment(0, fftSize - 1); widths = widths.max(1.0); - ArrayXXd diffs = + + ScopedEigenMap diffs(nChroma, fftSize, alloc); + diffs = freqs.replicate(1, nChroma).transpose() - ArrayXd::LinSpaced(nChroma, 0, nChroma - 1).transpose().replicate(fftSize, 1).transpose(); index halfChroma = std::lrint(nChroma / 2); - ArrayXXd remainder = diffs.unaryExpr([&](const double x){ + + ScopedEigenMap remainder(nChroma, fftSize, alloc); + remainder = diffs.unaryExpr([&](const double x){ return std::fmod(x + 10* nChroma + halfChroma, nChroma) - halfChroma; }); - MatrixXd filters = (-0.5 * (2 * remainder / widths.replicate(1, nChroma) + + mFilters.topLeftCorner(nChroma,nBins) = (-0.5 * (2 * remainder / widths.replicate(1, nChroma) .transpose()) .square()) .exp() .block(0, 0, nChroma, nBins); - - filters.colwise().normalize(); - - mFiltersStorage.setZero(); - mFiltersStorage.block(0, 0, nChroma, nBins) = filters; + + ScopedEigenMap colNorms(nBins, alloc); + colNorms = mFilters.topLeftCorner(nChroma, nBins).colwise().norm(); + mFilters.topLeftCorner(nChroma, nBins).array().rowwise() /= + colNorms.transpose(); + mNChroma = nChroma; mNBins = nBins; mScale = 2.0 / (fftSize * mNChroma); @@ -73,9 +85,9 @@ class ChromaFilterBank { using namespace Eigen; using namespace std; - ArrayXd frame = _impl::asEigen(in); - Eigen::Ref filters = mFiltersStorage.block(0, 0, mNChroma, mNBins); - + + FluidEigenMap frame = _impl::asEigen(in); + if(minFreq != 0 || maxFreq != -1){ maxFreq = (maxFreq == -1) ? (mSampleRate / 2) : min(maxFreq, mSampleRate / 2); double binHz = mSampleRate / ((mNBins - 1) * 2.); @@ -83,24 +95,32 @@ class ChromaFilterBank minFreq == 0 ? 0 : static_cast(ceil(minFreq / binHz)); index maxBin = min(static_cast(floorl(maxFreq / binHz)), (mNBins - 1)); - frame.segment(0, minBin).setZero(); - frame.segment(maxBin, frame.size() - maxBin).setZero(); + + using Eigen::seq; + using Eigen::seqN; + frame(seq(0,minBin),seqN(0,1)).setZero(); + frame(seqN(maxBin, frame.size() - maxBin), seqN(0,1)).setZero(); } + frame = frame.square(); - ArrayXd result = mScale * (filters * frame.square().matrix()).array(); + FluidEigenMap result = _impl::asEigen(out); + result.matrix().noalias() = + mFilters.topLeftCorner(mNChroma, mNBins) * frame.matrix(); + result *= mScale; if (normalize > 0) { double norm = normalize == 1? result.sum() : result.maxCoeff(); result = result / std::max(norm, epsilon); } - out <<= _impl::asFluid(result); } index mNChroma; index mNBins; double mScale; double mSampleRate; - Eigen::MatrixXd mFiltersStorage; + index mMaxChroma; + index mMaxFFT; + ScopedEigenMap mFilters; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/DCT.hpp b/include/algorithms/public/DCT.hpp index a7a0236c3..b3ad21757 100644 --- a/include/algorithms/public/DCT.hpp +++ b/include/algorithms/public/DCT.hpp @@ -13,6 +13,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/AlgorithmUtils.hpp" #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -27,44 +28,58 @@ class DCT using ArrayXd = Eigen::ArrayXd; using MatrixXd = Eigen::MatrixXd; - DCT(index maxInputSize, index maxOutputSize) - { - mTableStorage = MatrixXd::Zero(maxOutputSize, maxInputSize); - } + DCT(index maxInputSize, index maxOutputSize, + Allocator& alloc = FluidDefaultAllocator()) + : mTable{maxOutputSize, maxInputSize, alloc} + {} - void init(index inputSize, index outputSize) + void init(index inputSize, index outputSize, + Allocator& alloc = FluidDefaultAllocator()) { using namespace std; assert(inputSize >= outputSize); + assert(inputSize <= mTable.cols()); + assert(outputSize <= mTable.rows()); + mInputSize = inputSize; mOutputSize = outputSize; - mTable = mTableStorage.block(0, 0, mOutputSize, mInputSize); mTable.setZero(); + + ScopedEigenMap freqs(inputSize, alloc); for (index i = 0; i < mOutputSize; i++) { - double scale = i == 0 ? 1.0 / sqrt(inputSize) : sqrt(2.0 / inputSize); - ArrayXd freqs = ((pi / inputSize) * i) * - ArrayXd::LinSpaced(inputSize, 0.5, inputSize - 0.5); - mTable.row(i) = freqs.cos() * scale; + double scale = i == 0 ? 1.0 / sqrt(inputSize) : sqrt(2.0 / inputSize); + freqs = ((pi / inputSize) * i) * + ArrayXd::LinSpaced(inputSize, 0.5, inputSize - 0.5); + mTable.topLeftCorner(outputSize, inputSize).row(i) = freqs.cos() * scale; } + mInitialized = true; } void processFrame(const RealVectorView in, RealVectorView out) { - assert(in.size() == mInputSize); - ArrayXd frame = _impl::asEigen(in); - ArrayXd result = (mTable * frame.matrix()).array(); - out <<= _impl::asFluid(result); + assert(mInitialized && "DCT: processFrame() called before init()"); + assert(in.size() == mInputSize && + "DCT: actual input size doesn't match expected size"); + assert(out.size() == mOutputSize && + "DCT: actual output size doesn't maatch expected size"); + + FluidEigenMap frame = _impl::asEigen(in); + _impl::asEigen(out).noalias() = + (mTable.topLeftCorner(mOutputSize, mInputSize) * frame); } void processFrame(Eigen::Ref input, Eigen::Ref output) { - output = (mTable * input.matrix()).array(); + output.matrix().noalias() = + (mTable.topLeftCorner(mOutputSize, mInputSize) * input.matrix()); } - index mInputSize{40}; - index mOutputSize{13}; - MatrixXd mTable; - MatrixXd mTableStorage; + +private: + index mInputSize{40}; + index mOutputSize{13}; + bool mInitialized{false}; + ScopedEigenMap mTable; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/DataSetQuery.hpp b/include/algorithms/public/DataSetQuery.hpp index 9852c39e1..18bebe411 100644 --- a/include/algorithms/public/DataSetQuery.hpp +++ b/include/algorithms/public/DataSetQuery.hpp @@ -70,8 +70,8 @@ class DataSetQuery index maxColumn() { return mColumns.empty() ? 0 : *mColumns.rbegin(); } - bool addCondition(index column, string comparison, double value, - bool conjunction) + bool addCondition( + index column, string comparison, double value, bool conjunction) { auto pos = std::find(mComparisons.begin(), mComparisons.end(), comparison); if (pos == mComparisons.end()) return false; @@ -87,7 +87,7 @@ class DataSetQuery { auto data = input.getData(); auto ids = input.getIds(); - mTmpPoint = RealVector(asUnsigned(current.pointSize()) + mColumns.size()); + mTmpPoint = RealVector(current.pointSize() + asSigned(mColumns.size())); index limit = mLimit == 0 ? input.size() : mLimit; index count = 0; for (index i = 0; i < input.size() && count < limit; i++) @@ -128,8 +128,8 @@ class DataSetQuery } private: - void addRow(string id, RealVectorView point, const DataSet& current, - DataSet& out) + void addRow( + string id, RealVectorView point, const DataSet& current, DataSet& out) { mTmpPoint.fill(0); index currentSize = current.pointSize(); diff --git a/include/algorithms/public/Envelope.hpp b/include/algorithms/public/Envelope.hpp index c625ef15f..ee6e98ecb 100644 --- a/include/algorithms/public/Envelope.hpp +++ b/include/algorithms/public/Envelope.hpp @@ -23,9 +23,6 @@ namespace algorithm { class Envelope { - - using ArrayXd = Eigen::ArrayXd; - public: void init(double floor, double hiPassFreq) { diff --git a/include/algorithms/public/EnvelopeGate.hpp b/include/algorithms/public/EnvelopeGate.hpp index b19d4b607..21206c021 100644 --- a/include/algorithms/public/EnvelopeGate.hpp +++ b/include/algorithms/public/EnvelopeGate.hpp @@ -14,6 +14,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/FluidEigenMappings.hpp" #include "../util/SlideUDFilter.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -27,15 +28,14 @@ class EnvelopeGate using ArrayXd = Eigen::ArrayXd; public: - EnvelopeGate(index maxSize) - { - mInputStorage = ArrayXd(maxSize); - mOutputStorage = ArrayXd(maxSize); - } + EnvelopeGate(index maxSize, Allocator& alloc = FluidDefaultAllocator()) + : mInputBuffer(maxSize, alloc), + mOutputBuffer(maxSize, alloc) + {} void init(double onThreshold, double offThreshold, double hiPassFreq, - index minTimeAboveThreshold, index upwardLookupTime, - index minTimeBelowThreshold, index downwardLookupTime) + index minTimeAboveThreshold, index upwardLookupTime, + index minTimeBelowThreshold, index downwardLookupTime) { using namespace std; @@ -44,9 +44,10 @@ class EnvelopeGate mMinTimeBelowThreshold = minTimeBelowThreshold, mDownwardLookupTime = downwardLookupTime; mDownwardLatency = max(minTimeBelowThreshold, mDownwardLookupTime); - mLatency = max(mMinTimeAboveThreshold + mUpwardLookupTime, - mDownwardLatency); + mLatency = max( + mMinTimeAboveThreshold + mUpwardLookupTime, mDownwardLatency); if (mLatency < 0) mLatency = 1; + assert(mLatency <= mInputBuffer.size()); mHiPassFreq = hiPassFreq; initFilters(mHiPassFreq); double initVal = min(onThreshold, offThreshold) - 1; @@ -62,8 +63,8 @@ class EnvelopeGate } double processSample(const double in, double onThreshold, double offThreshold, - index rampUpTime, index rampDownTime, double hiPassFreq, - index minEventDuration, index minSilenceDuration) + index rampUpTime, index rampDownTime, double hiPassFreq, + index minEventDuration, index minSilenceDuration) { using namespace std; assert(mInitialized); @@ -121,7 +122,7 @@ class EnvelopeGate { index onsetIndex = refineStart(mWriteHead - mMinTimeAboveThreshold - mUpwardLookupTime, - mUpwardLookupTime); + mUpwardLookupTime); index blockSize = mWriteHead > onsetIndex ? mWriteHead - onsetIndex @@ -158,20 +159,20 @@ class EnvelopeGate } mOutputBuffer(mWriteHead) = mOutputState ? 1 : 0; - + mInputState = nextState; } mInputBuffer(mWriteHead) = smoothed; - + if (mFillCount < mLatency) mFillCount++; - double result = mOutputBuffer(mReadHead); - mWriteHead++; - mWriteHead = mWriteHead % mOutputBuffer.size(); - mReadHead++; - mReadHead = mReadHead % mOutputBuffer.size(); - return result; -} + double result = mOutputBuffer(mReadHead); + + if (++mWriteHead >= max(mLatency, 1)) mWriteHead = 0; + if (++mReadHead >= max(mLatency, 1)) mReadHead = 0; + + return result; + } index getLatency() { return mLatency; } bool initialized() { return mInitialized; } @@ -180,14 +181,12 @@ class EnvelopeGate void initBuffers(double initialValue) { using namespace std; - mInputBuffer = mInputStorage.segment(0, max(mLatency, 1)) - .setConstant(initialValue); - mOutputBuffer = - mOutputStorage.segment(0, max(mLatency, 1)).setZero(); + mInputBuffer.segment(0, max(mLatency, 1)).setConstant(initialValue); + mOutputBuffer.segment(0, max(mLatency, 1)).setZero(); mInputState = false; mOutputState = false; mFillCount = max(mLatency, 1); - mWriteHead = mOutputBuffer.size() - 1; + mWriteHead = max(mLatency, 1) - 1; mReadHead = 0; } @@ -269,12 +268,11 @@ class EnvelopeGate index mMinTimeBelowThreshold{10}; index mUpwardLookupTime{24}; - ArrayXd mInputBuffer; - ArrayXd mOutputBuffer; - ArrayXd mInputStorage; - ArrayXd mOutputStorage; - index mWriteHead; - index mReadHead; + ScopedEigenMap mInputBuffer; + ScopedEigenMap mOutputBuffer; + + index mWriteHead; + index mReadHead; bool mInputState{false}; bool mOutputState{false}; diff --git a/include/algorithms/public/EnvelopeSegmentation.hpp b/include/algorithms/public/EnvelopeSegmentation.hpp index 296c6bae1..644a301a9 100644 --- a/include/algorithms/public/EnvelopeSegmentation.hpp +++ b/include/algorithms/public/EnvelopeSegmentation.hpp @@ -24,21 +24,13 @@ namespace algorithm { class EnvelopeSegmentation { - - using ArrayXd = Eigen::ArrayXd; - public: void init(double floor, double hiPassFreq) { - // mFastSlide.init(floor); - // mSlowSlide.init(floor); - mEnvelope.init(floor,hiPassFreq); + mEnvelope.init(floor,hiPassFreq); mDebounceCount = 1; - // initFilters(hiPassFreq); - // mHiPassFreq = hiPassFreq; mPrevValue = 0; mState = false; - // mInitialized = true; } double processSample(const double in, double onThreshold, double offThreshold, @@ -46,27 +38,6 @@ class EnvelopeSegmentation index fastRampDownTime, index slowRampDownTime, double hiPassFreq, index debounce) { - // using namespace std; - // assert(mInitialized); - // mFastSlide.updateCoeffs(fastRampUpTime, fastRampDownTime); - // mSlowSlide.updateCoeffs(slowRampUpTime, slowRampDownTime); - // double filtered = in; - // if (hiPassFreq != mHiPassFreq) - // { - // initFilters(hiPassFreq); - // mHiPassFreq = hiPassFreq; - // } - // if (mHiPassFreq > 0){ - // filtered = mHiPass2.processSample(mHiPass1.processSample(in)); - // } - // double rectified = abs(filtered); - // double dB = 20 * log10(rectified); - // double clipped = max(dB, floor); - // double fast = mFastSlide.processSample(clipped); - // double slow = mSlowSlide.processSample(clipped); - // double value = fast - slow; - - double value = mEnvelope.processSample(in, floor, fastRampUpTime, slowRampUpTime, fastRampDownTime, slowRampDownTime, hiPassFreq); @@ -91,23 +62,10 @@ class EnvelopeSegmentation bool initialized() { return mEnvelope.initialized(); } private: -// void initFilters(double cutoff) -// { -// mHiPass1.init(cutoff); -// mHiPass2.init(cutoff); -// } - - Envelope mEnvelope; - // double mHiPassFreq{0}; + Envelope mEnvelope; index mDebounceCount{1}; double mPrevValue{0}; - // bool mInitialized{false}; bool mState{false}; - - // ButterworthHPFilter mHiPass1; - // ButterworthHPFilter mHiPass2; - // SlideUDFilter mFastSlide; - // SlideUDFilter mSlowSlide; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/HPS.hpp b/include/algorithms/public/HPS.hpp index 693f2eacd..3389a3d6d 100644 --- a/include/algorithms/public/HPS.hpp +++ b/include/algorithms/public/HPS.hpp @@ -24,15 +24,18 @@ class HPS public: void processFrame(const RealVectorView& input, RealVectorView output, index nHarmonics, double minFreq, double maxFreq, - double sampleRate) + double sampleRate, Allocator& alloc) { using namespace Eigen; using namespace std; ArrayXd::Index maxIndex; - ArrayXd mag = _impl::asEigen(input); - ArrayXd hps = mag; + ScopedEigenMap mag(input.size(), alloc); + mag = _impl::asEigen(input); + + ScopedEigenMap hps(mag.size(), alloc); + hps = mag; index nBins = mag.size(); double binHz = sampleRate / ((nBins - 1) * 2); index minBin = lrint(minFreq / binHz); @@ -41,21 +44,21 @@ class HPS double confidence = 0; double hpsSum = 0; + ScopedEigenMap h(nBins, alloc); + ScopedEigenMap hp(nBins, alloc); for (index i = 2; i < nHarmonics; i++) { index hBins = nBins / i; - ArrayXd h = ArrayXd::Zero(hBins); - for (index j = 0; j < hBins; j++) h(j) = mag(j * i); - ArrayXd hp = ArrayXd::Zero(nBins); - hp.segment(0, hBins) = h; + for (index j = 0; j < hBins; j++) h.head(hBins)(j) = mag(j * i); + hp.setZero(); + hp.head(hBins) = h.head(hBins); hps = hps * hp; } hpsSum = hps.sum(); if (maxBin > minBin && hpsSum > 0) { - hps = hps.segment(minBin, maxBin - minBin); - double maxVal = hps.maxCoeff(&maxIndex); + double maxVal = hps.segment(minBin, maxBin - minBin).maxCoeff(&maxIndex); confidence = maxVal / hpsSum; f0 = (minBin + maxIndex) * binHz; } diff --git a/include/algorithms/public/HPSS.hpp b/include/algorithms/public/HPSS.hpp index 8511b5641..3373e7b5f 100644 --- a/include/algorithms/public/HPSS.hpp +++ b/include/algorithms/public/HPSS.hpp @@ -15,6 +15,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/MedianFilter.hpp" #include "../../data/FluidIndex.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include namespace fluid { @@ -23,37 +24,41 @@ namespace algorithm { class HPSS { public: - using ArrayXXd = Eigen::ArrayXXd; - using ArrayXXcd = Eigen::ArrayXXcd; - using ArrayXcd = Eigen::ArrayXcd; - + + template + using Container = rt::vector; + enum HPSSMode { kClassic, kCoupled, kAdvanced }; - HPSS(index maxFFTSize, index maxHSize) - : mMaxH(maxFFTSize / 2 + 1, maxHSize), - mMaxV(maxFFTSize / 2 + 1, maxHSize), - mMaxBuf(maxFFTSize / 2 + 1, maxHSize) - { - mMaxH.setZero(); - mMaxV.setZero(); - mMaxBuf.setZero(); - } + HPSS(index maxFFTSize, index maxHSize,Allocator& alloc) + : mMaxBins(maxFFTSize / 2 + 1),mMaxHSize(maxHSize), + mHBuf(asUnsigned(mMaxBins * maxHSize), 0, alloc), + mVBuf(asUnsigned(mMaxBins * maxHSize), 0, alloc), + mFrameBuf(asUnsigned(mMaxBins * maxHSize), 0, alloc), + mPaddedBuf(asUnsigned(mMaxBins * 3), 0, alloc), + mHFilters(asUnsigned(mMaxBins),MedianFilter(mMaxHSize, alloc), alloc), + mVFilter(mMaxBins, alloc), + mHarmMaskBuf(asUnsigned(mMaxBins),alloc), + mPercMaskBuf(asUnsigned(mMaxBins),alloc), + mResMaskBuf(asUnsigned(mMaxBins),alloc), + mMaskNormBuf(asUnsigned(mMaxBins),alloc), + mMaskThreshBuf(asUnsigned(mMaxBins),alloc) + {} void init(index nBins, index hSize) { using namespace Eigen; assert(hSize % 2); - assert(nBins <= mMaxBuf.rows()); - assert(hSize <= mMaxBuf.cols()); + assert(nBins <= mMaxBins); + assert(hSize <= mMaxHSize); - mH = mMaxH.block(0, 0, nBins, hSize); - mV = mMaxV.block(0, 0, nBins, hSize); - mBuf = mMaxBuf.block(0, 0, nBins, hSize); - mH.setZero(); - mV.setZero(); - mBuf.setZero(); + ArrayXXMap v(mVBuf.data(),nBins, hSize); + ArrayXXMap h(mHBuf.data(),nBins,hSize); + ArrayXXcMap buf(mFrameBuf.data(),nBins,hSize); + h.setZero(); + v.setZero(); + buf.setZero(); - mHFilters = std::vector(asUnsigned(nBins)); for (index i = 0; i < nBins; i++) { mHFilters[asUnsigned(i)].init(hSize); } mInitialized = true; } @@ -67,46 +72,53 @@ class HPSS { using namespace Eigen; assert(mInitialized); + assert(vSize <= in.size()); + assert(in.size() <= mMaxBins); index h2 = (hSize - 1) / 2; index v2 = (vSize - 1) / 2; index nBins = in.size(); - ArrayXcd frame = _impl::asEigen(in); - ArrayXd mag = frame.abs().real(); - - mV.block(0, 0, nBins, hSize - 1) = mV.block(0, 1, nBins, hSize - 1); - mH.block(0, 0, nBins, hSize - 1) = mH.block(0, 1, nBins, hSize - 1); - mBuf.block(0, 0, nBins, hSize - 1) = mBuf.block(0, 1, nBins, hSize - 1); - - ArrayXd padded = ArrayXd::Zero(2 * vSize + nBins); - ArrayXd resultV = ArrayXd::Zero(padded.size()); - ArrayXd tmp = ArrayXd::Zero(padded.size()); - padded.segment(v2, nBins) = mag; + ArrayXXMap v(mVBuf.data(),nBins, hSize); + ArrayXXMap h(mHBuf.data(),nBins,hSize); + ArrayXXcMap buf(mFrameBuf.data(),nBins,hSize); + + v.block(0, 0, nBins, hSize - 1) = v.block(0, 1, nBins, hSize - 1); + h.block(0, 0, nBins, hSize - 1) = h.block(0, 1, nBins, hSize - 1); + buf.block(0, 0, nBins, hSize - 1) = buf.block(0, 1, nBins, hSize - 1); + + ArrayXMap padded(mPaddedBuf.data(),2 * vSize + nBins); + padded.setZero(); + + padded.segment(v2, nBins) = _impl::asEigen(in).abs().real(); mVFilter.init(vSize); for (index i = 0; i < padded.size(); i++) - { tmp(i) = mVFilter.processSample(padded(i)); } - mV.block(0, hSize - 1, nBins, 1) = tmp.segment(v2 * 3, nBins); - mBuf.block(0, hSize - 1, nBins, 1) = frame; - ArrayXd tmpRow = ArrayXd::Zero(2 * hSize); - for (index i = 0; i < nBins; i++) - { mH(i, h2 + 1) = mHFilters[asUnsigned(i)].processSample(mag(i)); } - ArrayXXcd result(nBins, 3); - ArrayXd harmonicMask = ArrayXd::Ones(nBins); - ArrayXd percussiveMask = ArrayXd::Ones(nBins); - ArrayXd residualMask = + { + padded(i) = mVFilter.processSample(padded(i)); + } + + v.block(0, hSize - 1, nBins, 1) = padded.segment(v2 * 3, nBins); + buf.block(0, hSize - 1, nBins, 1) = _impl::asEigen(in); + h.block(0, h2 + 1, nBins, 1) = ArrayXd::NullaryExpr(nBins,[&in,this](int i){ + return mHFilters[asUnsigned(i)].processSample(std::abs(in(i))); + }); + + ArrayXMap harmonicMask(mHarmMaskBuf.data(),nBins); + ArrayXMap percussiveMask(mPercMaskBuf.data(),nBins); + ArrayXMap residualMask(mResMaskBuf.data(),nBins); + residualMask = mode == kAdvanced ? ArrayXd::Ones(nBins) : ArrayXd::Zero(nBins); switch (mode) { case kClassic: { - ArrayXd HV = mH.col(0) + mV.col(0); - ArrayXd mult = (1.0 / HV.max(epsilon)); - harmonicMask = (mH.col(0) * mult); - percussiveMask = (mV.col(0) * mult); + ArrayXMap mult(mMaskNormBuf.data(),nBins); + mult = (1.0 / (h.col(0) + v.col(0)).max(epsilon)); + harmonicMask = (h.col(0) * mult); + percussiveMask = (v.col(0) * mult); break; } case kCoupled: { - harmonicMask = ((mH.col(0) / mV.col(0)) > + harmonicMask = ((h.col(0) / v.col(0)) > makeThreshold(nBins, hThresholdX1, hThresholdY1, hThresholdX2, hThresholdY2)) .cast(); @@ -114,17 +126,19 @@ class HPSS break; } case kAdvanced: { - harmonicMask = ((mH.col(0) / mV.col(0)) > + harmonicMask = ((h.col(0) / v.col(0)) > makeThreshold(nBins, hThresholdX1, hThresholdY1, hThresholdX2, hThresholdY2)) .cast(); - percussiveMask = ((mV.col(0) / mH.col(0)) > + percussiveMask = ((v.col(0) / h.col(0)) > makeThreshold(nBins, pThresholdX1, pThresholdY1, pThresholdX2, pThresholdY2)) .cast(); + residualMask = residualMask * (1 - harmonicMask); residualMask = residualMask * (1 - percussiveMask); - ArrayXd maskNorm = + ArrayXMap maskNorm(mMaskNormBuf.data(),nBins); + maskNorm = (1. / (harmonicMask + percussiveMask + residualMask)).max(epsilon); harmonicMask = harmonicMask * maskNorm; percussiveMask = percussiveMask * maskNorm; @@ -132,19 +146,20 @@ class HPSS break; } } - result.col(0) = mBuf.col(0) * harmonicMask.min(1.0); - result.col(1) = mBuf.col(0) * percussiveMask.min(1.0); - result.col(2) = mBuf.col(0) * residualMask.min(1.0); - out <<= _impl::asFluid(result); + _impl::asEigen(out).col(0) = buf.col(0) * harmonicMask.min(1.0); + _impl::asEigen(out).col(1) = buf.col(0) * percussiveMask.min(1.0); + _impl::asEigen(out).col(2) = buf.col(0) * residualMask.min(1.0); } + bool initialized() { return mInitialized; } private: - Eigen::ArrayXd makeThreshold(index nBins, double x1, double y1, double x2, + ArrayXMap makeThreshold(index nBins, double x1, double y1, double x2, double y2) { using namespace Eigen; - ArrayXd threshold = ArrayXd::Ones(nBins); + ArrayXMap threshold(mMaskThreshBuf.data(),nBins); + threshold = ArrayXd::Ones(nBins); index kneeStart = static_cast(std::floor(x1 * nBins)); index kneeEnd = static_cast(std::floor(x2 * nBins)); index kneeLength = kneeEnd - kneeStart; @@ -158,16 +173,20 @@ class HPSS return threshold; } - std::vector mHFilters; - MedianFilter mVFilter; - - ArrayXXd mMaxH; - ArrayXXd mMaxV; - ArrayXXcd mMaxBuf; - ArrayXXd mV; - ArrayXXd mH; - ArrayXXcd mBuf; - bool mInitialized{false}; + index mMaxBins; + index mMaxHSize; + Container mHBuf; + Container mVBuf; + Container> mFrameBuf; + Container mPaddedBuf; + Container mHFilters; + MedianFilter mVFilter; + Container mHarmMaskBuf; + Container mPercMaskBuf; + Container mResMaskBuf; + Container mMaskNormBuf; + Container mMaskThreshBuf; + bool mInitialized{false}; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/KDTree.hpp b/include/algorithms/public/KDTree.hpp index 6cbf71ab7..4d8e0bb9e 100644 --- a/include/algorithms/public/KDTree.hpp +++ b/include/algorithms/public/KDTree.hpp @@ -15,6 +15,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../../data/FluidIndex.hpp" #include "../../data/FluidTensor.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include #include #include @@ -34,8 +35,9 @@ class KDTree struct Node; using NodePtr = std::shared_ptr; using knnCandidate = std::pair; - using knnQueue = std::priority_queue, - std::less>; + using knnQueue = rt::vector; + using KNNResult = + std::pair, rt::vector>; using iterator = const std::vector::iterator; struct Node @@ -76,26 +78,25 @@ class KDTree mNPoints++; } - DataSet kNearest(ConstRealVectorView data, index k = 1, - double radius = 0) const + KNNResult kNearest(ConstRealVectorView data, index k = 1, double radius = 0, + Allocator& alloc = FluidDefaultAllocator()) const { assert(data.size() == mDims); - knnQueue queue; - auto result = DataSet(1); + rt::vector queue(alloc); + if (k > 0) queue.reserve(asUnsigned(k)); + kNearest(mRoot.get(), data, queue, k, radius, 0); - index numFound = asSigned(queue.size()); - std::vector sorted(asUnsigned(numFound)); - for (index i = numFound - 1; i >= 0; i--) - { - sorted[asUnsigned(i)] = queue.top(); - queue.pop(); - } - for (index i = 0; i < numFound; i++) - { - auto dist = FluidTensor{sorted[asUnsigned(i)].first}; - auto id = sorted[asUnsigned(i)].second->id; - result.add(id, dist); - } + std::sort_heap(queue.begin(), queue.end()); + + KNNResult result = + std::make_pair(rt::vector(queue.size(), alloc), + rt::vector(queue.size(), alloc)); + + std::for_each(queue.begin(), queue.end(), + [&result, i = 0u](knnCandidate const& x) mutable { + result.first[i] = x.first; + result.second[i++] = &(x.second->id); + }); return result; } @@ -205,11 +206,15 @@ class KDTree const double currentDist = distance(current->data, data); bool withinRadius = radius > 0 ? currentDist < radius : true; if (withinRadius && (knn.size() < asUnsigned(k) || k == 0)) - { knn.push(std::make_pair(currentDist, current)); } - else if (withinRadius && currentDist < knn.top().first) { - knn.pop(); - knn.push(std::make_pair(currentDist, current)); + knn.push_back(std::make_pair(currentDist, current)); + std::push_heap(knn.begin(), knn.end()); + } + else if (withinRadius && currentDist < knn.front().first) + { + std::pop_heap(knn.begin(), knn.end()); + knn.back() = std::make_pair(currentDist, current); + std::push_heap(knn.begin(), knn.end()); } const index d = depth % mDims; const double dimDif = current->data(d) - data(d); @@ -222,10 +227,12 @@ class KDTree } kNearest(firstBranch, data, knn, k, radius, depth + 1); if (k == 0 || knn.size() < asUnsigned(k) || - dimDif < knn.top().first) // ball centered at query with diametre - // kthDist intersects with current partition - // (or need to get more neighbors) - { kNearest(secondBranch, data, knn, k, radius, depth + 1); } + dimDif < knn.front().first) // ball centered at query with diametre + // kthDist intersects with current partition + // (or need to get more neighbors) + { + kNearest(secondBranch, data, knn, k, radius, depth + 1); + } } index flatten(index nodeId, const Node* current, FlatData& store) const diff --git a/include/algorithms/public/KMeans.hpp b/include/algorithms/public/KMeans.hpp index 6029e1402..04c53a1a9 100644 --- a/include/algorithms/public/KMeans.hpp +++ b/include/algorithms/public/KMeans.hpp @@ -16,6 +16,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../../data/FluidIndex.hpp" #include "../../data/FluidTensor.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include #include #include @@ -82,7 +83,8 @@ class KMeans index vq(RealVectorView point) const { assert(point.size() == mDims); - return assignPoint(_impl::asEigen(point)); + // transpose() allows us to avoid a temporary further down the call stack + return assignPoint(_impl::asEigen(point).transpose()); } void getMeans(RealMatrixView out) const @@ -119,12 +121,15 @@ class KMeans } protected: - double distance(const Eigen::ArrayXd& v1, const Eigen::ArrayXd& v2) const + template + double distance(const Eigen::ArrayBase& v1, + const Eigen::ArrayBase& v2) const { return (v1 - v2).matrix().norm(); } - index assignPoint(Eigen::ArrayXd point) const + template + index assignPoint(const Eigen::ArrayBase& point) const { double minDistance = std::numeric_limits::infinity(); index minK; @@ -140,15 +145,17 @@ class KMeans return minK; } - Eigen::VectorXi assignClusters(Eigen::ArrayXXd dataPoints) const + Eigen::VectorXi assignClusters(const Eigen::ArrayXXd& dataPoints) const { Eigen::VectorXi assignments = Eigen::VectorXi::Zero(dataPoints.rows()); for (index i = 0; i < dataPoints.rows(); i++) - { assignments(i) = static_cast(assignPoint(dataPoints.row(i))); } + { + assignments(i) = static_cast(assignPoint(dataPoints.row(i))); + } return assignments; } - void computeMeans(Eigen::ArrayXXd dataPoints) + void computeMeans(const Eigen::ArrayXXd& dataPoints) { using namespace Eigen; for (index k = 0; k < mK; k++) @@ -168,13 +175,15 @@ class KMeans ArrayXXd clusterPoints = ArrayXXd::Zero(asSigned(kAssignment.size()), mDims); for (index i = 0; asUnsigned(i) < kAssignment.size(); i++) - { clusterPoints.row(i) = dataPoints.row(kAssignment[asUnsigned(i)]); } + { + clusterPoints.row(i) = dataPoints.row(kAssignment[asUnsigned(i)]); + } ArrayXd mean = clusterPoints.colwise().mean(); mMeans.row(k) = mean; } } - bool changed(Eigen::VectorXi newAssignments) const + bool changed(const Eigen::VectorXi& newAssignments) const { auto dif = (newAssignments - mAssignments).cwiseAbs().sum(); return dif > 0; diff --git a/include/algorithms/public/KNNClassifier.hpp b/include/algorithms/public/KNNClassifier.hpp index 016b2ce8e..d68b41d10 100644 --- a/include/algorithms/public/KNNClassifier.hpp +++ b/include/algorithms/public/KNNClassifier.hpp @@ -17,6 +17,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../../data/FluidIndex.hpp" #include "../../data/FluidTensor.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include namespace fluid { @@ -28,62 +29,59 @@ class KNNClassifier public: using LabelSet = FluidDataSet; - std::string predict(KDTree tree, RealVectorView point, LabelSet labels, - index k, bool weighted) const + std::string const& predict(KDTree const& tree, RealVectorView point, + LabelSet const& labels, index k, bool weighted, + Allocator& alloc = FluidDefaultAllocator()) const { using namespace std; - unordered_map labelsMap; - auto nearest = tree.kNearest(point, k); - auto ids = nearest.getIds(); - auto distances = nearest.getData(); + unordered_map labelsMap; + auto [distances, ids] = tree.kNearest(point, k, 0, alloc); - double uniformWeight = 1.0 / k; - std::vector weights; - double sum = 0; + double uniformWeight = 1.0 / k; + rt::vector weights(asUnsigned(k), weighted ? 0 : uniformWeight, + alloc); + + double sum = 0; if (weighted) { - weights = std::vector(asUnsigned(k), 0); bool binaryWeights = false; - for (index i = 0; i < k; i++) + for (size_t i = 0; i < asUnsigned(k); i++) { - if (distances(i, 0) < epsilon) + if (distances[i] < epsilon) { binaryWeights = true; - weights[asUnsigned(i)] = 1; + weights[i] = 1; } else - sum += (1.0 / distances(i, 0)); + sum += (1.0 / distances[i]); } if (!binaryWeights) { - for (index i = 0; i < k; i++) - { weights[asUnsigned(i)] = (1.0 / distances(i, 0)) / sum; } + for (size_t i = 0; i < asUnsigned(k); i++) + { + weights[i] = (1.0 / distances[i]) / sum; + } } } - else - { - weights = std::vector(asUnsigned(k), uniformWeight); - } - string prediction; - FluidTensor tmp(1); - double maxWeight = 0; - for (index i = 0; i < k; i++) + string* prediction; + double maxWeight = 0; + for (size_t i = 0; i < asUnsigned(k); i++) { - labels.get(ids(i), tmp); - string label = tmp(0); - auto pos = labelsMap.find(label); + string* label = labels.get(*ids[i]).data(); + assert(label && "KNNClassifier: ID not mapped to label"); + auto pos = labelsMap.find(label); if (pos == labelsMap.end()) - labelsMap[label] = weights[asUnsigned(i)]; + labelsMap[label] = weights[i]; else - labelsMap[label] += weights[asUnsigned(i)]; + labelsMap[label] += weights[i]; if (labelsMap[label] > maxWeight) { maxWeight = labelsMap[label]; prediction = label; } } - return prediction; + return *prediction; } }; } // namespace algorithm diff --git a/include/algorithms/public/KNNRegressor.hpp b/include/algorithms/public/KNNRegressor.hpp index facb71a25..82dcf1251 100644 --- a/include/algorithms/public/KNNRegressor.hpp +++ b/include/algorithms/public/KNNRegressor.hpp @@ -17,6 +17,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../../data/FluidIndex.hpp" #include "../../data/FluidTensor.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include namespace fluid { @@ -28,46 +29,43 @@ class KNNRegressor public: using DataSet = FluidDataSet; - double predict(KDTree tree, DataSet targets, RealVectorView point, index k, - bool weighted) const + double predict(KDTree const& tree, DataSet const& targets, + RealVectorView point, index k, bool weighted, + Allocator& alloc = FluidDefaultAllocator()) const { using namespace std; - auto nearest = tree.kNearest(point, k); - double prediction = 0; - auto ids = nearest.getIds(); - auto distances = nearest.getData(); - double uniformWeight = 1.0 / k; - std::vector weights; - double sum = 0; + double prediction = 0; + auto [distances, ids] = tree.kNearest(point, k, 0, alloc); + double uniformWeight = 1.0 / k; + rt::vector weights(asUnsigned(k), weighted ? 0 : uniformWeight, + alloc); + double sum = 0; if (weighted) { - weights = std::vector(asUnsigned(k), 0); bool binaryWeights = false; - for (index i = 0; i < k; i++) + for (size_t i = 0; i < asUnsigned(k); i++) { - if (distances(i, 0) < epsilon) + if (distances[i] < epsilon) { binaryWeights = true; - weights[asUnsigned(i)] = 1; + weights[i] = 1; } else - sum += (1.0 / distances(i, 0)); + sum += (1.0 / distances[i]); } if (!binaryWeights) { - for (index i = 0; i < k; i++) - { weights[asUnsigned(i)] = (1.0 / distances(i, 0)) / sum; } + for (size_t i = 0; i < asUnsigned(k); i++) + { + weights[i] = (1.0 / distances[i]) / sum; + } } } - else - { - weights = std::vector(asUnsigned(k), uniformWeight); - } - for (index i = 0; i < k; i++) + + for (size_t i = 0; i < asUnsigned(k); i++) { - auto point = FluidTensor(1); - targets.get(ids(i), point); - prediction += (weights[asUnsigned(i)] * point(0)); + auto point = targets.get(*ids[i]); + prediction += (weights[i] * point(0)); } return prediction; } diff --git a/include/algorithms/public/LabelSetEncoder.hpp b/include/algorithms/public/LabelSetEncoder.hpp index b01ae448c..5f35ebfc1 100644 --- a/include/algorithms/public/LabelSetEncoder.hpp +++ b/include/algorithms/public/LabelSetEncoder.hpp @@ -44,7 +44,7 @@ class LabelSetEncoder mInitialized = true; } - index encodeIndex(string label) const + index encodeIndex(string const& label) const { auto pos = mLabelsMap.find(label); if (pos != mLabelsMap.end()) @@ -53,6 +53,7 @@ class LabelSetEncoder return -1; } + //todo: is this ever used? Will it be? std::string decodeIndex(index in) const { if (in < mLabels.size()) @@ -61,7 +62,7 @@ class LabelSetEncoder return ""; } - void encodeOneHot(string label, RealVectorView out) const + void encodeOneHot(string const& label, RealVectorView out) const { assert(out.size() == mNumLabels); RealVector result(mNumLabels); @@ -71,7 +72,7 @@ class LabelSetEncoder if (pos != mLabelsMap.end()) out(pos->second) = 1.0; } - std::string decodeOneHot(RealVectorView in) const + std::string const& decodeOneHot(RealVectorView in) const { double maxVal = 0; index maxIndex = 0; @@ -83,15 +84,16 @@ class LabelSetEncoder maxVal = in(i); } } + assert(maxIndex < mLabels.size()); return mLabels(maxIndex); } index numLabels() const { return mNumLabels; } // from JSON: expecting unique labels, as opposed to fit - void init(FluidTensor labels) + void init(FluidTensorView labels) { mLabelsMap.clear(); - mLabels = labels; + mLabels <<= labels; for (index i = 0; i < mLabels.size(); i++) { mLabelsMap[mLabels(i)] = i; } mNumLabels = mLabels.size(); mInitialized = true; diff --git a/include/algorithms/public/Loudness.hpp b/include/algorithms/public/Loudness.hpp index b9cf9367c..84a8fce1e 100644 --- a/include/algorithms/public/Loudness.hpp +++ b/include/algorithms/public/Loudness.hpp @@ -15,6 +15,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/KWeightingFilter.hpp" #include "../util/TruePeak.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -26,30 +27,34 @@ class Loudness { public: - Loudness(index maxSize) : mTP(maxSize) {} + Loudness(index maxSize, Allocator& alloc = FluidDefaultAllocator()) + : mTP(maxSize, alloc) + {} - void init(index size, double sampleRate) + void init( + index size, double sampleRate, Allocator& alloc = FluidDefaultAllocator()) { mFilter.init(sampleRate); - mTP.init(size, sampleRate); + mTP.init(size, sampleRate, alloc); mSize = size; mInitialized = true; } void processFrame(const RealVectorView& input, RealVectorView output, - bool weighting, bool truePeak) + bool weighting, bool truePeak, Allocator& alloc = FluidDefaultAllocator()) { using namespace Eigen; using namespace std; assert(mInitialized); assert(output.size() == 2); assert(input.size() == mSize); - ArrayXd in = _impl::asEigen(input); - ArrayXd filtered(mSize); + FluidEigenMap in = _impl::asEigen(input); + ScopedEigenMap filtered(mSize, alloc); for (index i = 0; i < mSize; i++) - filtered(i) = weighting ? mFilter.processSample(in(i)) : in(i); + filtered(i) = weighting ? mFilter.processSample(input(i)) : input(i); double loudness = -0.691 + 10 * log10(filtered.square().mean() + epsilon); - double peak = truePeak ? mTP.processFrame(input) : in.abs().maxCoeff(); + double peak = + truePeak ? mTP.processFrame(input, alloc) : in.abs().maxCoeff(); peak = 20 * log10(peak + epsilon); output(0) = loudness; output(1) = peak; diff --git a/include/algorithms/public/MLP.hpp b/include/algorithms/public/MLP.hpp index 652c8e353..640f95194 100644 --- a/include/algorithms/public/MLP.hpp +++ b/include/algorithms/public/MLP.hpp @@ -17,6 +17,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../../data/FluidIndex.hpp" #include "../../data/FluidTensor.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include #include @@ -44,6 +45,7 @@ class MLP activations.push_back(hiddenAct); } sizes.push_back(outputSize); + mMaxLayerSize = *std::max_element(sizes.begin(), sizes.end()); activations.push_back(outputAct); for (index i = 0; i < asSigned(sizes.size() - 1); i++) { @@ -100,17 +102,16 @@ class MLP } void processFrame(RealVectorView in, RealVectorView out, index startLayer, - index endLayer) const + index endLayer, + Allocator& alloc = FluidDefaultAllocator()) const { using namespace _impl; using namespace Eigen; - ArrayXd tmpIn = asEigen(in); - ArrayXXd input(1, tmpIn.size()); - input.row(0) = tmpIn; - ArrayXXd output = ArrayXXd::Zero(1, out.size()); - forward(input, output, startLayer, endLayer); - ArrayXd tmpOut = output.row(0); - out <<= asFluid(tmpOut); + ScopedEigenMap input(in.size(), alloc); + input = asEigen(in); + ScopedEigenMap output(out.size(), alloc); + forwardFrame(input, output, startLayer, endLayer, alloc); + asEigen(out) = output; } void forward(Eigen::Ref in, Eigen::Ref out) const @@ -137,6 +138,31 @@ class MLP out = output; } + void forwardFrame(Eigen::Ref in, Eigen::Ref out, + index startLayer, index endLayer, + Allocator& alloc = FluidDefaultAllocator()) const + { + if (startLayer >= asSigned(mLayers.size()) || + endLayer > asSigned(mLayers.size())) + return; + if (startLayer < 0 || endLayer <= 0) return; + ScopedEigenMap input(mMaxLayerSize, alloc); + input.head(in.size()) = in; + ScopedEigenMap output(mMaxLayerSize, alloc); + index inSize = in.size(); + for (index i = startLayer; i < endLayer; i++) + { + auto& l = mLayers[asUnsigned(i)]; + auto inputBlock = input.head(inSize); + auto outputBlock = output.head(l.outputSize()); + outputBlock.setZero(); + l.forwardFrame(inputBlock, outputBlock, alloc); + input.head(l.outputSize()) = outputBlock; + inSize = l.outputSize(); + } + out = output.head(out.size()); + } + void backward(Eigen::Ref out) { index nRows = out.rows(); @@ -192,6 +218,7 @@ class MLP std::vector mLayers; bool mInitialized{false}; bool mTrained{false}; + index mMaxLayerSize; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/MelBands.hpp b/include/algorithms/public/MelBands.hpp index e81617c34..c8ffa084c 100644 --- a/include/algorithms/public/MelBands.hpp +++ b/include/algorithms/public/MelBands.hpp @@ -13,6 +13,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/AlgorithmUtils.hpp" #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include #include #include @@ -23,8 +24,9 @@ namespace algorithm { class MelBands { public: - MelBands(index maxBands, index maxFFT) - : mFiltersStorage(maxBands, maxFFT / 2 + 1) + MelBands( + index maxBands, index maxFFT, Allocator& alloc = FluidDefaultAllocator()) + : mFilters(maxBands, maxFFT / 2 + 1, alloc) {} /*static inline double mel2hz(double x) { @@ -37,61 +39,70 @@ class MelBands } void init(double lo, double hi, index nBands, index nBins, double sampleRate, - index windowSize) + index windowSize, Allocator& alloc = FluidDefaultAllocator()) { using namespace Eigen; assert(hi > lo); assert(nBands > 1); + assert(nBins <= mFilters.cols()); mScale1 = 1.0 / (windowSize / 4.0); // scale to original amplitude index fftSize = 2 * (nBins - 1); + mScale2 = 1.0 / (2.0 * double(fftSize) / windowSize); - ArrayXd melFreqs = ArrayXd::LinSpaced(nBands + 2, hz2mel(lo), hz2mel(hi)); + ScopedEigenMap melFreqs(nBands + 2, alloc); + melFreqs = ArrayXd::LinSpaced(nBands + 2, hz2mel(lo), hz2mel(hi)); melFreqs = 700.0 * ((melFreqs / 1127.01048).exp() - 1.0); - mFilters = mFiltersStorage.block(0, 0, nBands, nBins); - mFilters.setZero(); - ArrayXd fftFreqs = ArrayXd::LinSpaced(nBins, 0, sampleRate / 2.0); - ArrayXd melD = - (melFreqs.segment(0, nBands + 1) - melFreqs.segment(1, nBands + 1)) - .abs(); - ArrayXXd ramps = melFreqs.replicate(1, nBins); + // mFilters = mFiltersStorage.block(0, 0, nBands, nBins); + mFilters.topLeftCorner(nBands, nBins).setZero(); + ScopedEigenMap fftFreqs(nBins, alloc); + fftFreqs = ArrayXd::LinSpaced(nBins, 0, sampleRate / 2.0); + ScopedEigenMap melD(nBands + 1, alloc); + melD = (melFreqs.segment(0, nBands + 1) - melFreqs.segment(1, nBands + 1)) + .abs(); + ScopedEigenMap ramps(melFreqs.rows(), nBins, alloc); + ramps = melFreqs.replicate(1, nBins); ramps.rowwise() -= fftFreqs.transpose(); + + ScopedEigenMap lower(nBins, alloc); + ScopedEigenMap upper(nBins, alloc); for (index i = 0; i < nBands; i++) { - ArrayXd lower = -ramps.row(i) / melD(i); - ArrayXd upper = ramps.row(i + 2) / melD(i + 1); - mFilters.row(i) = lower.min(upper).max(0); + lower = -ramps.row(i) / melD(i); + upper = ramps.row(i + 2) / melD(i + 1); + mFilters.row(i).head(nBins) = lower.min(upper).max(0); } + mNBands = nBands; + mNBins = nBins; } void processFrame(const RealVectorView in, RealVectorView out, bool magNorm, - bool usePower, bool logOutput) + bool usePower, bool logOutput, Allocator&) { using namespace Eigen; - ArrayXd frame = _impl::asEigen(in); + FluidEigenMap frame = _impl::asEigen(in); + FluidEigenMap result = _impl::asEigen(out); + if (magNorm) frame = frame * mScale1; - ArrayXd result; - if (usePower) { result = (mFilters * frame.square().matrix()).array(); } - else - { - result = (mFilters * frame.matrix()).array(); - } - if (magNorm) - { - double energy = frame.sum() * mScale2; - result = result * energy / std::max(epsilon, result.sum()); - } + double energy = frame.sum() * mScale2; + if (usePower) frame = frame.square(); + + result.matrix().noalias() = + (mFilters.topLeftCorner(mNBands, mNBins) * frame.matrix()); + + if (magNorm) { result = result * energy / std::max(epsilon, result.sum()); } if (logOutput) result = 20 * result.max(epsilon).log10(); - out <<= _impl::asFluid(result); } double mScale1{1.0}; double mScale2{1.0}; - Eigen::MatrixXd mFilters; - Eigen::MatrixXd mFiltersStorage; +private: + ScopedEigenMap mFilters; + index mNBands; + index mNBins; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/NMF.hpp b/include/algorithms/public/NMF.hpp index 98df683e3..e439b56b4 100644 --- a/include/algorithms/public/NMF.hpp +++ b/include/algorithms/public/NMF.hpp @@ -14,6 +14,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include #include @@ -34,47 +35,57 @@ class NMF using namespace Eigen; using namespace _impl; - MatrixXd W1 = asEigen(W).transpose(); - MatrixXd H1 = asEigen(H).transpose(); - MatrixXd result = (W1.col(idx) * H1.row(idx)).transpose(); - V <<= asFluid(result); + auto W1 = asEigen(W).transpose(); + auto H1 = asEigen(H).transpose(); + asEigen(V).transpose().noalias() = (W1.col(idx) * H1.row(idx)); } // processFrame computes activations of a dictionary W in a given frame void processFrame(const RealVectorView x, const RealMatrixView W0, - RealVectorView out, index nIterations = 10, - RealVectorView v = RealVectorView(nullptr, 0, 0)) + RealVectorView out, index nIterations, + RealVectorView v, Allocator& alloc) { using namespace Eigen; using namespace _impl; index rank = W0.extent(0); - MatrixXd W = asEigen(W0).transpose(); - VectorXd h = - MatrixXd::Random(rank, 1) * 0.5 + MatrixXd::Constant(rank, 1, 0.5); - VectorXd v0 = asEigen(x); + FluidEigenMap W = asEigen(W0); + + ScopedEigenMap h(rank, alloc); + h = VectorXd::Random(rank) * 0.5 + VectorXd::Constant(rank, 0.5); + + ScopedEigenMap v0(x.size(), alloc); + v0 = asEigen(x); W = W.array().max(epsilon).matrix(); h = h.array().max(epsilon).matrix(); v0 = v0.array().max(epsilon).matrix(); - MatrixXd WT = W.transpose(); - W.colwise().normalize(); - VectorXd ones = VectorXd::Ones(x.extent(0)); +// MatrixXd WT = W.transpose(); + ScopedEigenMap norm(W.rowwise().norm(), alloc); + W.array().colwise() /= norm.array(); + index nBins = x.extent(0); + + ScopedEigenMap v1{nBins, alloc}; + ScopedEigenMap vRatio{nBins, alloc}; + ScopedEigenMap hNum{rank, alloc}; + ScopedEigenMap hDen{rank, alloc}; + auto ones = VectorXd::Ones(nBins); while (nIterations--) { - ArrayXd v1 = (W * h).array().max(epsilon); - ArrayXXd hNum = (WT * (v0.array() / v1).matrix()).array(); - ArrayXXd hDen = (WT * ones).array(); + v1.matrix().noalias() = (W.transpose() * h); + v1 = v1.max(epsilon); + vRatio = v0.array() / v1; + hNum.matrix().noalias() = W * vRatio.matrix(); + hDen.matrix().noalias() = W * ones; h = (h.array() * hNum / hDen.max(epsilon)).matrix(); // VectorXd r = W * h; // double divergence = (v.cwiseProduct(v.cwiseQuotient(r)) - v + r).sum(); // std::cout<<"Divergence "< 0) - { - ArrayXd v2 = (W * h).array(); - v <<= asFluid(v2); - } + + if(out.data()) + _impl::asEigen(out) = h; + + if (v.data()) _impl::asEigen(v).noalias() = (W.transpose() * h); } void process(const RealMatrixView X, RealMatrixView W1, RealMatrixView H1, diff --git a/include/algorithms/public/NMFMorph.hpp b/include/algorithms/public/NMFMorph.hpp index 59309436c..b5ea740cc 100644 --- a/include/algorithms/public/NMFMorph.hpp +++ b/include/algorithms/public/NMFMorph.hpp @@ -16,6 +16,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/Munkres.hpp" #include "../util/OptimalTransport.hpp" #include "../util/RTPGHI.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -30,39 +31,56 @@ class NMFMorph public: using MatrixXd = Eigen::MatrixXd; + NMFMorph(index maxFFTSize, Allocator& alloc) + : mW1(0, 0, alloc), mW2(0, 0, alloc), mH(0, 0, alloc), + mRTPGHI(maxFFTSize, alloc), mOT(alloc) + {} + void init(RealMatrixView W1, RealMatrixView W2, RealMatrixView H, - index winSize, index fftSize, index hopSize, bool assign) + index winSize, index fftSize, index hopSize, bool assign, + Allocator& alloc) { using namespace Eigen; using namespace _impl; using namespace std; mInitialized = false; + index maxSize = std::max(W1.cols(), W2.cols()); + mW1 = ScopedEigenMap(W1.cols(), W1.rows(), alloc); mW1 = asEigen(W1).transpose(); + + mH = ScopedEigenMap(H.rows(), H.cols(), alloc); mH = asEigen(H); - MatrixXd tmpW2 = asEigen(W2).transpose(); - ArrayXXd cost = ArrayXXd::Zero(mW1.cols(), tmpW2.cols()); + if (assign) { + ScopedEigenMap tmpW2(W2.cols(), W2.rows(), alloc); + tmpW2 = asEigen(W2).transpose(); + ScopedEigenMap cost(mW1.cols(), tmpW2.cols(), alloc); + cost.setZero(); + OptimalTransport tmpOT(maxSize, alloc); for (index i = 0; i < mW1.cols(); i++) { for (index j = 0; j < tmpW2.cols(); j++) { - OptimalTransport tmpOT; - tmpOT.init(mW1.col(i), tmpW2.col(j)); - if(!tmpOT.initialized()) return; + tmpOT.init(mW1.col(i), tmpW2.col(j), alloc); + if (!tmpOT.initialized()) return; cost(i, j) = tmpOT.mDistance; } } - Munkres munk; - munk.init(mW1.cols(), tmpW2.cols()); - ArrayXi result = ArrayXi::Zero(mW1.cols()); - munk.process(cost, result); - mW2 = MatrixXd::Zero(tmpW2.rows(), tmpW2.cols()); + Munkres munk(mW1.cols(), tmpW2.cols(), alloc); + ScopedEigenMap result(mW1.cols(), alloc); + result.setZero(); + munk.process(cost, result, alloc); + mW2 = ScopedEigenMap(tmpW2.rows(), tmpW2.cols(), alloc); + mW2.setZero(); for (index i = 0; i < result.size(); i++) - { mW2.col(i) = tmpW2.col(result(i)); } + { + mW2.col(i) = tmpW2.col(result(i)); + } } else { + mW2 = ScopedEigenMap(W2.cols(), W2.rows(), alloc); mW2 = asEigen(W2).transpose(); } mWindowSize = winSize; @@ -71,44 +89,52 @@ class NMFMorph mRTPGHI.init(fftSize); index rank = mW1.cols(); - mOT = std::vector(asUnsigned(rank)); - for (index i = 0; i < rank; i++) { mOT[asUnsigned(i)].init(mW1.col(i), mW2.col(i)); } + mOT = rt::vector(alloc); + mOT.reserve(asUnsigned(rank)); + for (index i = 0; i < rank; i++) + { + mOT.emplace_back(maxSize, alloc); + mOT.back().init(mW1.col(i), mW2.col(i), alloc); + } mPos = 0; mInitialized = true; } - + bool initialized() { return mInitialized; } - void processFrame(ComplexVectorView v, double interpolation) + void processFrame(ComplexVectorView v, double interpolation, Allocator& alloc) { using namespace Eigen; using namespace _impl; - MatrixXd W = MatrixXd::Zero(mW1.rows(), mW1.cols()); + ScopedEigenMap W(mW1.rows(), mW1.cols(), alloc); + W.setZero(); + ScopedEigenMap out(mW2.rows(), alloc); for (int i = 0; i < W.cols(); i++) { - ArrayXd out = ArrayXd::Zero(mW2.rows()); + out.setZero(); mOT[asUnsigned(i)].interpolate(interpolation, out); W.col(i) = out; } - - VectorXd hFrame = mH.col(mPos); - VectorXd frame = W * hFrame; + ScopedEigenMap hFrame(mH.rows(), alloc); + hFrame = mH.col(mPos); + ScopedEigenMap frame(W.rows(), alloc); + frame = W * hFrame; RealVectorView mag1 = asFluid(frame); - mRTPGHI.processFrame(mag1, v, mWindowSize, mFFTSize, mHopSize, 1e-6); + mRTPGHI.processFrame(mag1, v, mWindowSize, mFFTSize, mHopSize, 1e-6, alloc); mPos = (mPos + 1) % mH.cols(); } private: - MatrixXd mW1; - MatrixXd mW2; - MatrixXd mH; - index mWindowSize; - index mHopSize; - index mFFTSize; - RTPGHI mRTPGHI; - std::vector mOT; - int mPos{0}; - bool mInitialized; + ScopedEigenMap mW1; + ScopedEigenMap mW2; + ScopedEigenMap mH; + index mWindowSize; + index mHopSize; + index mFFTSize; + RTPGHI mRTPGHI; + rt::vector mOT; + int mPos{0}; + bool mInitialized; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/Normalization.hpp b/include/algorithms/public/Normalization.hpp index f265db8bc..e5d17a24a 100644 --- a/include/algorithms/public/Normalization.hpp +++ b/include/algorithms/public/Normalization.hpp @@ -61,8 +61,8 @@ class Normalization { using namespace Eigen; using namespace _impl; - ArrayXd input = asEigen(in); - ArrayXd result; + FluidEigenMap input = asEigen(in); + FluidEigenMap result = asEigen(out); if (!inverse) { result = (input - mDataMin) / mDataRange; @@ -73,7 +73,6 @@ class Normalization result = (input - mMin) / mRange; result = mDataMin + (result * mDataRange); } - out <<= asFluid(result); } void process(const RealMatrixView in, RealMatrixView out, diff --git a/include/algorithms/public/NoveltyFeature.hpp b/include/algorithms/public/NoveltyFeature.hpp index 6590c9631..3ba121572 100644 --- a/include/algorithms/public/NoveltyFeature.hpp +++ b/include/algorithms/public/NoveltyFeature.hpp @@ -13,6 +13,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/FluidEigenMappings.hpp" #include "../util/Novelty.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include @@ -25,38 +26,43 @@ class NoveltyFeature public: using ArrayXd = Eigen::ArrayXd; - NoveltyFeature(index maxKernelSize, index maxFilterSize) - : mFilterBufferStorage(maxFilterSize), mNovelty(maxKernelSize) + NoveltyFeature(index maxKernelSize, index maxDims, index maxFilterSize, + Allocator& alloc) + : mFilterBuffer(maxFilterSize, alloc), + mNovelty(maxKernelSize, maxDims, alloc) {} - void init(index kernelSize, index filterSize, index nDims) + void init(index kernelSize, index filterSize, index nDims, Allocator& alloc) { assert(kernelSize % 2); - mNovelty.init(kernelSize, nDims); - mFilterBuffer = mFilterBufferStorage.segment(0, filterSize); - mFilterBuffer.setZero(); + mNovelty.init(kernelSize, nDims, alloc); + mFilterBuffer.head(filterSize).setZero(); + mFilterSize = filterSize; + mInitialized = true; } - double processFrame(const RealVectorView input) + double processFrame(const RealVectorView input, Allocator& alloc) { - double novelty = mNovelty.processFrame(_impl::asEigen(input)); - index filterSize = mFilterBuffer.size(); + assert(mInitialized); + double novelty = + mNovelty.processFrame(_impl::asEigen(input), alloc); - if (filterSize > 1) + if (mFilterSize > 1) { - mFilterBuffer.segment(0, filterSize - 1) = - mFilterBuffer.segment(1, filterSize - 1); + mFilterBuffer.segment(0, mFilterSize - 1) = + mFilterBuffer.segment(1, mFilterSize - 1); } - mFilterBuffer(filterSize - 1) = novelty; + mFilterBuffer(mFilterSize - 1) = novelty; - return mFilterBuffer.mean(); + return mFilterBuffer.head(mFilterSize).mean(); } private: - ArrayXd mFilterBuffer; - ArrayXd mFilterBufferStorage; - Novelty mNovelty; + bool mInitialized; + ScopedEigenMap mFilterBuffer; + Novelty mNovelty; + index mFilterSize; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/NoveltySegmentation.hpp b/include/algorithms/public/NoveltySegmentation.hpp index abb6480de..c63dba497 100644 --- a/include/algorithms/public/NoveltySegmentation.hpp +++ b/include/algorithms/public/NoveltySegmentation.hpp @@ -13,6 +13,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "NoveltyFeature.hpp" #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include @@ -25,24 +26,27 @@ class NoveltySegmentation public: using ArrayXd = Eigen::ArrayXd; - NoveltySegmentation(index maxKernelSize, index maxFilterSize) - : mNovelty(maxKernelSize, maxFilterSize) + NoveltySegmentation(index maxKernelSize, index maxDims, index maxFilterSize, + Allocator& alloc = FluidDefaultAllocator()) + : mNovelty(maxKernelSize, maxDims, maxFilterSize, alloc), + mPeakBuffer(3, alloc) {} - void init(index kernelSize, index filterSize, index nDims) + void init(index kernelSize, index filterSize, index nDims, + Allocator& alloc = FluidDefaultAllocator()) { - mNovelty.init(kernelSize, filterSize, nDims); + mNovelty.init(kernelSize, filterSize, nDims, alloc); mDebounceCount = 1; - mPeakBuffer.setZero(); + mPeakBuffer.setZero(); } double processFrame(const RealVectorView input, double threshold, - index minSliceLength) + index minSliceLength, Allocator& alloc = FluidDefaultAllocator()) { double detected = 0.; mPeakBuffer.segment(0, 2) = mPeakBuffer.segment(1, 2); - mPeakBuffer(2) = mNovelty.processFrame(input); + mPeakBuffer(2) = mNovelty.processFrame(input, alloc); if (mPeakBuffer(1) > mPeakBuffer(0) && mPeakBuffer(1) > mPeakBuffer(2) && mPeakBuffer(1) > threshold && mDebounceCount == 0) @@ -58,9 +62,9 @@ class NoveltySegmentation } private: - NoveltyFeature mNovelty; - ArrayXd mPeakBuffer{3}; - index mDebounceCount{1}; + NoveltyFeature mNovelty; + ScopedEigenMap mPeakBuffer; + index mDebounceCount{1}; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/OnsetDetectionFunctions.hpp b/include/algorithms/public/OnsetDetectionFunctions.hpp index 29359e2b8..9f6258ed2 100644 --- a/include/algorithms/public/OnsetDetectionFunctions.hpp +++ b/include/algorithms/public/OnsetDetectionFunctions.hpp @@ -16,6 +16,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/MedianFilter.hpp" #include "../util/OnsetDetectionFuncs.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -31,16 +32,20 @@ class OnsetDetectionFunctions using WindowTypes = WindowFuncs::WindowTypes; public: - OnsetDetectionFunctions(index maxSize) - : mFFT(maxSize), mWindowStorage(maxSize) + OnsetDetectionFunctions(index maxSize, index maxFilterSize, Allocator& alloc) + : mFFT(maxSize, alloc), mWindow(maxSize, alloc), + mPrevFrame(maxSize / 2 + 1, alloc), + mPrevPrevFrame(maxSize / 2 + 1, alloc), mFilter(maxFilterSize, alloc) {} void init(index windowSize, index fftSize, index filterSize) { + assert(windowSize <= mWindow.size()); makeWindow(windowSize); - prevFrame = ArrayXcd::Zero(fftSize / 2 + 1); - prevPrevFrame = ArrayXcd::Zero(fftSize / 2 + 1); - mFilter.init(filterSize); + + mPrevFrame.setZero(); + mPrevPrevFrame.setZero(); + mFilter.init(std::max(filterSize, 3)); mFFT.resize(fftSize); mDebounceCount = 1; mPrevFuncVal = 0; @@ -49,60 +54,71 @@ class OnsetDetectionFunctions void makeWindow(index windowSize) { - mWindowStorage.setZero(); - WindowFuncs::map()[mWindowType](windowSize, mWindowStorage); - mWindow = mWindowStorage.segment(0, windowSize); + mWindow.setZero(); + WindowFuncs::map()[mWindowType](windowSize, mWindow.head(windowSize)); mWindowSize = windowSize; } + index nextPower2(index n) + { + return static_cast(std::pow(2, std::ceil(std::log2(n)))); + } + + /// input window isn't necessarily a single framre because it should encompass + /// `frameDelta`'s worth of history double processFrame(RealVectorView input, index function, index filterSize, - index frameDelta = 0) + index frameDelta, Allocator& alloc) { assert(mInitialized); - ArrayXd in = _impl::asEigen(input); - double funcVal = 0; - double filteredFuncVal = 0; - + double funcVal = 0; + double filteredFuncVal = 0; + index frameSize = nextPower2(mWindowSize) / 2 + 1; + if (filterSize >= 3 && (!mFilter.initialized() || filterSize != mFilter.size())) mFilter.init(filterSize); - ArrayXcd frame = mFFT.process(in.segment(0, mWindowSize) * mWindow); - auto odf = static_cast(function); + ScopedEigenMap in(mWindowSize, alloc); + in = _impl::asEigen(input).col(0).head(mWindowSize) * + mWindow.head(mWindowSize); + + ScopedEigenMap frame(frameSize, alloc); + frame = mFFT.process(in); + if (function > 1 && function < 5 && frameDelta != 0) { - ArrayXcd frame2 = - mFFT.process(in.segment(frameDelta, mWindowSize) * mWindow); - funcVal = OnsetDetectionFuncs::map()[odf](frame2, frame, frame); + ScopedEigenMap frame2(frameSize, alloc); + frame2 = mFFT.process(in.segment(frameDelta, mWindowSize) * + mWindow.head(mWindowSize)); + funcVal = OnsetDetectionFuncs::map(function)(frame2, frame, frame, alloc); } else { - funcVal = - OnsetDetectionFuncs::map()[odf](frame, prevFrame, prevPrevFrame); + funcVal = OnsetDetectionFuncs::map(function)( + frame, mPrevFrame.head(frameSize), mPrevPrevFrame.head(frameSize), + alloc); } if (filterSize >= 3) filteredFuncVal = funcVal - mFilter.processSample(funcVal); else filteredFuncVal = funcVal - mPrevFuncVal; - prevPrevFrame = prevFrame; - prevFrame = frame; - + mPrevPrevFrame.head(frameSize) = mPrevFrame.head(frameSize); + mPrevFrame.head(frameSize) = frame; return filteredFuncVal; } private: - FFT mFFT{1024}; - ArrayXd mWindowStorage; - ArrayXd mWindow; - index mWindowSize{1024}; - index mDebounceCount{1}; - ArrayXcd prevFrame; - ArrayXcd prevPrevFrame; - double mPrevFuncVal{0.0}; - WindowTypes mWindowType{WindowTypes::kHann}; - MedianFilter mFilter; - bool mInitialized{false}; + FFT mFFT; + ScopedEigenMap mWindow; + index mWindowSize{1024}; + index mDebounceCount{1}; + ScopedEigenMap mPrevFrame; + ScopedEigenMap mPrevPrevFrame; + double mPrevFuncVal{0.0}; + WindowTypes mWindowType{WindowTypes::kHann}; + MedianFilter mFilter; + bool mInitialized{false}; }; } // namespace algorithm diff --git a/include/algorithms/public/OnsetSegmentation.hpp b/include/algorithms/public/OnsetSegmentation.hpp index 8492ba1da..fe1d1cf12 100644 --- a/include/algorithms/public/OnsetSegmentation.hpp +++ b/include/algorithms/public/OnsetSegmentation.hpp @@ -31,7 +31,9 @@ class OnsetSegmentation using ArrayXd = Eigen::ArrayXd; using ArrayXcd = Eigen::ArrayXcd; - OnsetSegmentation(index maxSize) : mODF{maxSize} {} + OnsetSegmentation(index maxSize, index maxFilterSize, Allocator& alloc = FluidDefaultAllocator()) + : mODF{maxSize, maxFilterSize, alloc} + {} void init(index windowSize, index fftSize, index filterSize) { @@ -39,12 +41,14 @@ class OnsetSegmentation mDebounceCount = 1; } + /// input window isn't necessarily a single framre because it should encompass + /// `frameDelta`'s worth of history double processFrame(RealVectorView input, index function, index filterSize, - double threshold, index debounce = 0, - index frameDelta = 0) + double threshold, index debounce, index frameDelta, + Allocator& alloc = fluid::FluidDefaultAllocator()) { double filteredFuncVal = - mODF.processFrame(input, function, filterSize, frameDelta); + mODF.processFrame(input, function, filterSize, frameDelta, alloc); double detected{0}; if (filteredFuncVal > threshold && mPrevFuncVal < threshold && diff --git a/include/algorithms/public/PCA.hpp b/include/algorithms/public/PCA.hpp index 28362f551..7bb849567 100644 --- a/include/algorithms/public/PCA.hpp +++ b/include/algorithms/public/PCA.hpp @@ -59,16 +59,18 @@ class PCA using namespace Eigen; using namespace _impl; if (k > mBases.cols()) return; - VectorXd input = asEigen(in); + FluidEigenMap input = asEigen(in); input = input - mMean; - VectorXd result = input.transpose() * mBases.block(0, 0, mBases.rows(), k); + FluidEigenMap result = asEigen(out); + // this avoids Eigen making a temporary (sorry): + result.transpose().noalias() = + input.transpose() * mBases.block(0, 0, mBases.rows(), k); if (whiten) { - ArrayXd norm = mExplainedVariance.segment(0, k).max(epsilon).rsqrt(); + auto norm = mExplainedVariance.segment(0, k).max(epsilon).rsqrt(); result.array() *= norm; } - out <<= _impl::asFluid(result); } void inverseProcessFrame(RealVectorView in, RealVectorView out, bool whiten = false) const diff --git a/include/algorithms/public/RatioMask.hpp b/include/algorithms/public/RatioMask.hpp index 95996e0ed..cd6bf95eb 100644 --- a/include/algorithms/public/RatioMask.hpp +++ b/include/algorithms/public/RatioMask.hpp @@ -14,6 +14,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include namespace fluid { @@ -25,29 +26,41 @@ class RatioMask using ArrayXXd = Eigen::ArrayXXd; public: + RatioMask(index maxRows, index maxCols, Allocator& alloc) + : mMultiplier(maxRows, maxCols, alloc) + {} + void init(RealMatrixView denominator) { using namespace _impl; using namespace Eigen; - mMultiplier = (1 / asEigen(denominator).max(epsilon)); + mRows = denominator.rows(); + mCols = denominator.cols(); + mMultiplier.topLeftCorner(mRows, mCols) = + (1 / asEigen(denominator).max(epsilon)); + mInitialized = true; } void process(const ComplexMatrixView& mixture, RealMatrixView targetMag, - index exponent, ComplexMatrixView result) + index exponent, ComplexMatrixView out) { using namespace _impl; using namespace Eigen; + assert(mInitialized); assert(mixture.cols() == targetMag.cols()); assert(mixture.rows() == targetMag.rows()); - ArrayXXcd tmp = + asEigen(out) = asEigen(mixture) * - (asEigen(targetMag).pow(exponent) * mMultiplier.pow(exponent)) + (asEigen(targetMag).pow(exponent) * + mMultiplier.topLeftCorner(mRows, mCols).pow(exponent)) .min(1.0); - result <<= asFluid(tmp); } private: - ArrayXXd mMultiplier; + ScopedEigenMap mMultiplier; + bool mInitialized{false}; + index mRows; + index mCols; }; } // namespace algorithm diff --git a/include/algorithms/public/RobustScaling.hpp b/include/algorithms/public/RobustScaling.hpp index bb52c800b..c4e8ab3c9 100644 --- a/include/algorithms/public/RobustScaling.hpp +++ b/include/algorithms/public/RobustScaling.hpp @@ -73,14 +73,13 @@ class RobustScaling { using namespace Eigen; using namespace _impl; - ArrayXd input = asEigen(in); - ArrayXd result; + FluidEigenMap input = asEigen(in); + FluidEigenMap result = asEigen(out); if (!inverse) { result = (input - mMedian) / mRange; } else { result = (input * mRange) + mMedian; } - out <<= asFluid(result); } void process(const RealMatrixView in, RealMatrixView out, diff --git a/include/algorithms/public/STFT.hpp b/include/algorithms/public/STFT.hpp index f6d68c13a..daf09e5d3 100644 --- a/include/algorithms/public/STFT.hpp +++ b/include/algorithms/public/STFT.hpp @@ -15,6 +15,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/FFT.hpp" #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/FluidTensor.hpp" #include "../../data/TensorTypes.hpp" #include @@ -29,36 +30,53 @@ class STFT using ArrayXXd = Eigen::ArrayXXd; using ArrayXcd = Eigen::ArrayXcd; using ArrayXXcd = Eigen::ArrayXXcd; + using ArrayXdMap = Eigen::Map; public: - STFT(index windowSize, index fftSize, index hopSize, index windowType = 0) + STFT(index windowSize, index fftSize, index hopSize, index windowType = 0, + Allocator& alloc = FluidDefaultAllocator()) : mWindowSize(windowSize), mHopSize(hopSize), mFrameSize(fftSize / 2 + 1), - mFFT(fftSize) + mMaxWindowSize(windowSize), + mWindowType(static_cast(windowType)), + mWindowBuffer(asUnsigned(mMaxWindowSize), alloc), + mWindowedFrameBuffer(asUnsigned(mMaxWindowSize), alloc), + mFFT(fftSize, alloc) { - mWindow = ArrayXd::Zero(mWindowSize); - auto windowTypeIndex = static_cast(windowType); - WindowFuncs::map()[windowTypeIndex](mWindowSize, mWindow); + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + WindowFuncs::map()[mWindowType](mWindowSize, window); + } + + void resize(index windowSize, index fftSize, index hopSize) + { + assert(windowSize <= mMaxWindowSize && + "STFT: Window Size greater than Max"); + mWindowSize = windowSize; + mHopSize = hopSize; + mFrameSize = fftSize / 2 + 1; + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + WindowFuncs::map()[mWindowType](mWindowSize, window); + mFFT.resize(fftSize); } static void magnitude(const FluidTensorView, 2> in, FluidTensorView out) { - ArrayXXd mag = _impl::asEigen(in).abs().real(); - out <<= _impl::asFluid(mag); + _impl::asEigen(out) = + _impl::asEigen(in).abs().real(); } static void magnitude(const FluidTensorView, 1> in, FluidTensorView out) { - ArrayXd mag = _impl::asEigen(in).abs().real(); - out <<= _impl::asFluid(mag); + _impl::asEigen(out) = + _impl::asEigen(in).abs().real(); } static void phase(const FluidTensorView, 2> in, FluidTensorView out) { - ArrayXXd phase = _impl::asEigen(in).arg().real(); - out <<= _impl::asFluid(phase); + _impl::asEigen(out) = + _impl::asEigen(in).arg().real(); } static void phase(const FluidTensorView, 1> in, @@ -71,8 +89,9 @@ class STFT void process(const RealVectorView audio, ComplexMatrixView spectrogram) { - index halfWindow = mWindowSize / 2; - ArrayXd padded(audio.size() + mWindowSize + mHopSize); + index halfWindow = mWindowSize / 2; + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + ArrayXd padded(audio.size() + mWindowSize + mHopSize); padded.fill(0); padded.segment(halfWindow, audio.size()) = Eigen::Map(audio.data(), audio.size()); @@ -83,7 +102,7 @@ class STFT for (index i = 0; i < nFrames; i++) { result.row(i) = - mFFT.process(padded.segment(i * mHopSize, mWindowSize) * mWindow); + mFFT.process(padded.segment(i * mHopSize, mWindowSize) * window); } spectrogram <<= _impl::asFluid(result); } @@ -91,29 +110,37 @@ class STFT void processFrame(const RealVectorView frame, ComplexVectorView out) { assert(frame.size() == mWindowSize); - ArrayXcd spectrum = - mFFT.process(_impl::asEigen(frame) * mWindow); - out <<= _impl::asFluid(spectrum); + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + ArrayXdMap windowedFrame(mWindowedFrameBuffer.data(), mWindowSize); + windowedFrame = _impl::asEigen(frame); + windowedFrame *= window; + _impl::asEigen(out) = mFFT.process(windowedFrame); } void processFrame(Eigen::Ref frame, Eigen::Ref out) { assert(frame.size() == mWindowSize); - out = mFFT.process(frame * mWindow); + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + ArrayXdMap windowedFrame(mWindowedFrameBuffer.data(), mWindowSize); + windowedFrame = frame; + windowedFrame *= window; + out = mFFT.process(windowedFrame); } - RealVectorView window() { - return RealVectorView(mWindow.data(), 0, mWindowSize); + return RealVectorView(mWindowBuffer.data(), 0, mWindowSize); } private: - index mWindowSize; - index mHopSize; - index mFrameSize; - ArrayXd mWindow; - FFT mFFT; + index mWindowSize; + index mHopSize; + index mFrameSize; + index mMaxWindowSize; + WindowFuncs::WindowTypes mWindowType; + rt::vector mWindowBuffer; + rt::vector mWindowedFrameBuffer; + FFT mFFT; }; class ISTFT @@ -121,67 +148,85 @@ class ISTFT using ArrayXd = Eigen::ArrayXd; using ArrayXcd = Eigen::ArrayXcd; using ArrayXXcd = Eigen::ArrayXXcd; + using ArrayXdMap = Eigen::Map; public: - ISTFT(index windowSize, index fftSize, index hopSize, index windowType = 0) - : mWindowSize(windowSize), mHopSize(hopSize), mScale(1 / double(fftSize)), - mIFFT(fftSize), mBuffer(mWindowSize) + ISTFT(index windowSize, index fftSize, index hopSize, index windowType = 0, + Allocator& alloc = FluidDefaultAllocator()) + : mWindowSize(windowSize), mMaxWindowSize(windowSize), mHopSize(hopSize), + mScale(1 / double(fftSize)), + mWindowType(static_cast(windowType)), + mIFFT(fftSize, alloc), mBuffer(asUnsigned(mMaxWindowSize), alloc), + mWindowBuffer(asUnsigned(mMaxWindowSize), alloc) { - mWindow = ArrayXd::Zero(mWindowSize); - auto windowTypeIndex = static_cast(windowType); - WindowFuncs::map()[windowTypeIndex](mWindowSize, mWindow); - mWindowSquared = mWindow * mWindow; + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + WindowFuncs::map()[mWindowType](mWindowSize, window); + } + + void resize(index windowSize, index fftSize, index hopSize) + { + assert(windowSize <= mMaxWindowSize && + "STFT: Window Size greater than Max"); + mWindowSize = windowSize; + mHopSize = hopSize; + mScale = 1 / double(fftSize); + mIFFT.resize(fftSize); + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + WindowFuncs::map()[mWindowType](mWindowSize, window); } void process(const ComplexMatrixView spectrogram, RealVectorView audio) { const auto& epsilon = std::numeric_limits::epsilon; - - index halfWindow = mWindowSize / 2; - index nFrames = spectrogram.rows(); - index outputSize = mWindowSize + (nFrames - 1) * mHopSize; + index halfWindow = mWindowSize / 2; + index nFrames = spectrogram.rows(); + index outputSize = mWindowSize + (nFrames - 1) * mHopSize; outputSize += mWindowSize + mHopSize; - ArrayXXcd specData = _impl::asEigen(spectrogram); - ArrayXd outputPadded = ArrayXd::Zero(outputSize); - ArrayXd norm = ArrayXd::Zero(outputSize); + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + ArrayXXcd specData = _impl::asEigen(spectrogram); + ArrayXd outputPadded = ArrayXd::Zero(outputSize); + ArrayXd norm = ArrayXd::Zero(outputSize); for (index i = 0; i < nFrames; i++) { ArrayXd frame = mIFFT.process(specData.row(i)).segment(0, mWindowSize); outputPadded.segment(i * mHopSize, mWindowSize) += - frame * mScale * mWindow; - norm.segment(i * mHopSize, mWindowSize) += mWindow * mWindow; + frame * mScale * window; + norm.segment(i * mHopSize, mWindowSize) += window * window; } outputPadded = outputPadded / norm.max(epsilon()); ArrayXd trimmed = outputPadded.segment(halfWindow, audio.size()); audio <<= _impl::asFluid(trimmed); } - void processFrame(const ComplexVectorView frame, RealVectorView audio) + void processFrame(ComplexVectorView frame, RealVectorView audio) { - mBuffer = mIFFT.process(_impl::asEigen(frame)) - .segment(0, mWindowSize) * - mWindow * mScale; - audio <<= _impl::asFluid(mBuffer); + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + Eigen::Map frameMap(mBuffer.data(), frame.size()); + frameMap = _impl::asEigen(frame); + _impl::asEigen(audio) = + mIFFT.process(frameMap).head(window.size()) * window * mScale; } void processFrame(Eigen::Ref frame, Eigen::Ref audio) { - audio = mIFFT.process(frame).segment(0, mWindowSize) * mWindow * mScale; + ArrayXdMap window(mWindowBuffer.data(), mWindowSize); + audio = mIFFT.process(frame).segment(0, mWindowSize) * window * mScale; } RealVectorView window() { - return RealVectorView(mWindow.data(), 0, mWindowSize); + return RealVectorView(mWindowBuffer.data(), 0, mWindowSize); } private: - index mWindowSize{1024}; - index mHopSize{512}; - ArrayXd mWindow; - ArrayXd mWindowSquared; - double mScale{1}; - IFFT mIFFT; - ArrayXd mBuffer; + index mWindowSize{1024}; + index mMaxWindowSize; + index mHopSize{512}; + double mScale{1}; + WindowFuncs::WindowTypes mWindowType; + IFFT mIFFT; + rt::vector> mBuffer; + rt::vector mWindowBuffer; }; } // namespace algorithm diff --git a/include/algorithms/public/SineExtraction.hpp b/include/algorithms/public/SineExtraction.hpp index abdad3ba9..72ec0a322 100644 --- a/include/algorithms/public/SineExtraction.hpp +++ b/include/algorithms/public/SineExtraction.hpp @@ -12,12 +12,12 @@ under the European Union’s Horizon 2020 research and innovation programme #include "WindowFuncs.hpp" #include "../util/AlgorithmUtils.hpp" -#include "../util/ConvolutionTools.hpp" #include "../util/FFT.hpp" #include "../util/FluidEigenMappings.hpp" #include "../util/PartialTracking.hpp" #include "../util/PeakDetection.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -33,16 +33,30 @@ class SineExtraction using VectorXd = Eigen::VectorXd; using ArrayXcd = Eigen::ArrayXcd; template - using vector = std::vector; + using vector = rt::vector; + + using Queue = rt::queue>; + + Queue makeEmptyQueue(Allocator& alloc) + { + return Queue{rt::deque>(alloc)}; + } + public: - void init(index windowSize, index fftSize, index transformSize) + SineExtraction(index maxFFT, Allocator& alloc) + : mTracking(alloc), mBuf{makeEmptyQueue(alloc)}, + mWindowTransform(maxFFT, alloc) + {} + + void init(index windowSize, index fftSize, index transformSize, + Allocator& alloc) { mBins = fftSize / 2 + 1; mCurrentFrame = 0; - mBuf = std::queue(); + // mBuf = std::queue(); mScale = 1.0 / (windowSize / 4.0); // scale to original amplitude - computeWindowTransform(windowSize, transformSize); + computeWindowTransform(windowSize, transformSize, alloc); mTracking.init(); mWindowBinIncr = mWindowTransform.size() / (mBins - 1) / 2; mInvWindowBinIncr = 1.0 / mWindowBinIncr; @@ -53,19 +67,28 @@ class SineExtraction double sampleRate, double detectionThreshold, index minTrackLength, double birthLowThreshold, double birthHighThreshold, index trackMethod, double zetaA, - double zetaF, double delta, index bandwidth) + double zetaF, double delta, index bandwidth, + Allocator& alloc) { assert(mInitialized); using namespace Eigen; - index fftSize = 2 * (mBins - 1); - ArrayXcd frame = _impl::asEigen(in); + index fftSize = 2 * (mBins - 1); + ScopedEigenMap frame(in.size(), alloc); + frame = _impl::asEigen(in); + if (minTrackLength != mTracking.minTrackLength()) - { mBuf = std::queue(); } + { + mBuf = makeEmptyQueue(alloc); + } + mBuf.push(frame); - ArrayXd mag = frame.abs().real(); + ScopedEigenMap mag(in.size(), alloc); + mag = frame.abs().real(); mag = mag * mScale; - ArrayXd logMag = 20 * mag.max(epsilon).log10(); - vector peaks; + ScopedEigenMap logMag(in.size(), alloc); + logMag = 20 * mag.max(epsilon).log10(); + + vector peaks(0, alloc); auto tmpPeaks = mPeakDetection.process(logMag, 0, -infinity, true, false); for (auto p : tmpPeaks) { @@ -75,24 +98,29 @@ class SineExtraction peaks.push_back({hz, p.second, false}); } } + double maxAmp = 20 * std::log10(mag.maxCoeff()); mTracking.processFrame(peaks, maxAmp, minTrackLength, birthLowThreshold, - birthHighThreshold, trackMethod, zetaA, zetaF, - delta); - vector sinePeaks = mTracking.getActivePeaks(); - ArrayXd frameSines = ArrayXd::Zero(mBins); + birthHighThreshold, trackMethod, zetaA, zetaF, delta, + alloc); + vector sinePeaks = mTracking.getActivePeaks(alloc); + ScopedEigenMap frameSines(mBins, alloc); + frameSines.setZero(); for (auto& p : sinePeaks) - { frameSines += synthesizePeak(p, sampleRate, bandwidth); } - ArrayXXcd result(mBins, 2); + { + frameSines += synthesizePeak(p, sampleRate, bandwidth, alloc); + } + ScopedEigenMap result(mBins, 2, alloc); if (asSigned(mBuf.size()) <= mTracking.minTrackLength()) { - result.col(0) = ArrayXd::Zero(mBins); - result.col(1) = ArrayXd::Zero(mBins); + result.col(0).setZero(); + result.col(1).setZero(); } else { - ArrayXcd resultFrame = mBuf.front(); - ArrayXd resultMag = resultFrame.abs().real(); + ScopedEigenMap& resultFrame = mBuf.front(); + ScopedEigenMap resultMag(resultFrame.size(), alloc); + resultMag = resultFrame.abs().real(); for (index i = 0; i < mBins; i++) { if (frameSines(i) >= resultMag(i)) @@ -110,25 +138,25 @@ class SineExtraction mBuf.pop(); } mTracking.prune(); - out <<= _impl::asFluid(result); + _impl::asEigen(out) = result; mCurrentFrame++; } - void reset() { mCurrentFrame = 0; } bool initialized() { return mInitialized; } private: - void computeWindowTransform(index windowSize, index transformSize) + void computeWindowTransform(index windowSize, index transformSize, + Allocator& alloc) { index halfBW = transformSize / 2; - mWindowTransform = ArrayXd::Zero(transformSize); - ArrayXd window = ArrayXd::Zero(windowSize); - FFT fft(transformSize); + mWindowTransform.head(transformSize) = ArrayXd::Zero(transformSize); + ScopedEigenMap window(ArrayXd::Zero(windowSize), alloc); + FFT fft(transformSize); WindowFuncs::map()[WindowFuncs::WindowTypes::kHann](windowSize, window); - ArrayXcd transform = - fft.process(Eigen::Map(window.data(), windowSize)); + ScopedEigenMap transform(transformSize / 2 + 1, alloc); + transform = fft.process(Eigen::Map(window.data(), windowSize)); for (index i = 0; i < halfBW; i++) { mWindowTransform(halfBW + i) = mWindowTransform(halfBW - i) = @@ -144,12 +172,13 @@ class SineExtraction return mWindowTransform(floor) + frac * mInvWindowBinIncr * dY; } - ArrayXd synthesizePeak(SinePeak p, double sampleRate, index bandwidth) + ScopedEigenMap synthesizePeak(SinePeak p, double sampleRate, + index bandwidth, Allocator& alloc) { using namespace std; - index halfBW = bandwidth / 2; - ArrayXd sine = ArrayXd::Zero(mBins); - double freqBin = p.freq * 2 * (mBins - 1) / sampleRate; + index halfBW = bandwidth / 2; + ScopedEigenMap sine(mBins, alloc); + double freqBin = p.freq * 2 * (mBins - 1) / sampleRate; if (freqBin >= mBins - 1) freqBin = mBins - 1; if (freqBin < 0) freqBin = 0; index freqBinFloor = lrint(floor(freqBin)); @@ -160,26 +189,30 @@ class SineExtraction for (index i = freqBinCeil; pos < mWindowTransform.size() - 2 && i < min(freqBinCeil + halfBW, mBins - 1); i++, pos += mWindowBinIncr) - { sine[i] = amp * interpolateWindow(pos); } + { + sine[i] = amp * interpolateWindow(pos); + } pos = (mWindowTransform.size() / 2) - ((freqBin - freqBinFloor) * mWindowBinIncr); for (index i = freqBinFloor; pos > 1 && i > max(freqBinFloor - halfBW, asSigned(0)); i--, pos -= mWindowBinIncr) - { sine[i] = amp * interpolateWindow(pos); } + { + sine[i] = amp * interpolateWindow(pos); + } return sine; } - PeakDetection mPeakDetection; - PartialTracking mTracking; - index mBins{513}; - index mCurrentFrame{0}; - std::queue mBuf; - ArrayXd mWindowTransform; - double mScale{1.0}; - bool mInitialized{false}; - double mWindowBinIncr; - double mInvWindowBinIncr; + PeakDetection mPeakDetection; + PartialTracking mTracking; + index mBins{513}; + index mCurrentFrame{0}; + Queue mBuf; + ScopedEigenMap mWindowTransform; + double mScale{1.0}; + bool mInitialized{false}; + double mWindowBinIncr; + double mInvWindowBinIncr; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/public/SineFeature.hpp b/include/algorithms/public/SineFeature.hpp new file mode 100644 index 000000000..634298351 --- /dev/null +++ b/include/algorithms/public/SineFeature.hpp @@ -0,0 +1,89 @@ +/* +Part of the Fluid Corpus Manipulation Project (http://www.flucoma.org/) +Copyright 2017-2019 University of Huddersfield. +Licensed under the BSD-3 License. +See license.md file in the project root for full license information. +This project has received funding from the European Research Council (ERC) +under the European Union’s Horizon 2020 research and innovation programme +(grant agreement No 725899). +*/ + +#pragma once + +#include "WindowFuncs.hpp" +#include "../util/AlgorithmUtils.hpp" +#include "../util/FFT.hpp" +#include "../util/FluidEigenMappings.hpp" +#include "../util/PartialTracking.hpp" +#include "../util/PeakDetection.hpp" +#include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" +#include "../../data/TensorTypes.hpp" +#include +#include +#include + +namespace fluid { +namespace algorithm { + + +class SineFeature +{ + using ArrayXd = Eigen::ArrayXd; + using VectorXd = Eigen::VectorXd; + using ArrayXcd = Eigen::ArrayXcd; + template + using vector = rt::vector; + +public: + SineFeature(Allocator&) {} + + void init(index windowSize, index fftSize) + { + mBins = fftSize / 2 + 1; + mScale = 1.0 / (windowSize / 4.0); // scale to original amplitude + mInitialized = true; + } + + index processFrame(const ComplexVectorView in, RealVectorView freqOut, + RealVectorView magOut, double sampleRate, + double detectionThreshold, index sortBy, Allocator& alloc) + { + assert(mInitialized); + using namespace Eigen; + index fftSize = 2 * (mBins - 1); + ScopedEigenMap frame(in.size(), alloc); + frame = _impl::asEigen(in); + + ScopedEigenMap mag(in.size(), alloc); + mag = frame.abs().real(); + mag = mag * mScale; + ScopedEigenMap logMagIn(in.size(), alloc); + logMagIn = 20 * mag.max(epsilon).log10(); + + auto tmpPeaks = + mPeakDetection.process(logMagIn, 0, detectionThreshold, true, sortBy); + + index maxNumOut = std::min(freqOut.size(), asSigned(tmpPeaks.size())); + + double ratio = sampleRate / fftSize; + std::transform(tmpPeaks.begin(), tmpPeaks.begin() + maxNumOut, + freqOut.begin(), + [ratio](auto peak) { return peak.first * ratio; }); + + std::transform(tmpPeaks.begin(), tmpPeaks.begin() + maxNumOut, + magOut.begin(), [](auto peak) { return peak.second; }); + + return maxNumOut; + } + + bool initialized() { return mInitialized; } + +private: + PeakDetection mPeakDetection; + index mBins{513}; + double mScale{1.0}; + bool mInitialized{false}; +}; +} // namespace algorithm +} // namespace fluid diff --git a/include/algorithms/public/SpectralShape.hpp b/include/algorithms/public/SpectralShape.hpp index c0ecff2fa..b87668db1 100644 --- a/include/algorithms/public/SpectralShape.hpp +++ b/include/algorithms/public/SpectralShape.hpp @@ -13,6 +13,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/AlgorithmUtils.hpp" #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -26,19 +27,20 @@ class SpectralShape using ArrayXd = Eigen::ArrayXd; public: - SpectralShape() {} + SpectralShape(Allocator& alloc) : mOutputBuffer(7, alloc) {} void processFrame(Eigen::Ref in, double sampleRate, double minFreq, double maxFreq, double rolloffTarget, bool logFreq, - bool usePower) + bool usePower, Allocator& alloc) { using namespace std; maxFreq = (maxFreq == -1) ? (sampleRate / 2) : min(maxFreq, sampleRate / 2); - ArrayXd mag = in.max(epsilon); - index nBins = mag.size(); - double binHz = sampleRate / ((nBins - 1) * 2.); - index minBin = static_cast(ceil(minFreq / binHz)); - index maxBin = + ScopedEigenMap mag(in.size(), alloc); + mag = in.max(epsilon); + index nBins = mag.size(); + double binHz = sampleRate / ((nBins - 1) * 2.); + index minBin = static_cast(ceil(minFreq / binHz)); + index maxBin = min(static_cast(floorl(maxFreq / binHz)), (nBins - 1)); if (maxBin <= minBin) { @@ -51,14 +53,20 @@ class SpectralShape minBin = 1; mag(1) += mag(0); } - ArrayXd amp = usePower ? mag.square() : mag; - amp = amp.segment(minBin, maxBin - minBin).eval(); - - double ampSum = amp.sum(); - ArrayXd freqs = - ArrayXd::LinSpaced(maxBin - minBin, minBin * binHz, maxBin * binHz); + + index size = maxBin - minBin; + + ScopedEigenMap amp(size, alloc); + amp = mag.segment(minBin, size); + if (usePower) amp = amp.square(); + + double ampSum = amp.sum(); + ScopedEigenMap freqs(size, alloc); + freqs = ArrayXd::LinSpaced(size, minBin * binHz, maxBin * binHz); if (logFreq) - { freqs = 69 + (12 * (freqs / 440).log() * log2E); } // MIDI cents + { + freqs = 69 + (12 * (freqs / 440).log() * log2E); + } // MIDI cents double centroid = (amp * freqs).sum() / ampSum; double spread = (amp * (freqs - centroid).square()).sum() / ampSum; @@ -93,20 +101,21 @@ class SpectralShape mOutputBuffer(6) = 20 * log10(max(crest, epsilon)); } - void processFrame(const RealVector& input, RealVectorView output, - double sampleRate, double minFreq = 0, double maxFreq = -1, - double rolloffTarget = 0.95, bool logFreq = false, - bool usePower = false) + void processFrame(const RealVectorView input, RealVectorView output, + double sampleRate, double minFreq, double maxFreq, + double rolloffTarget, bool logFreq, bool usePower, + Allocator& alloc) { assert(output.size() == 7); - ArrayXd in = _impl::asEigen(input); + ScopedEigenMap in(input.size(), alloc); + in = _impl::asEigen(input); processFrame(in, sampleRate, minFreq, maxFreq, rolloffTarget, logFreq, - usePower); - output <<= _impl::asFluid(mOutputBuffer); + usePower, alloc); + _impl::asEigen(output) = mOutputBuffer; } private: - ArrayXd mOutputBuffer{7}; + ScopedEigenMap mOutputBuffer; }; } // namespace algorithm diff --git a/include/algorithms/public/Standardization.hpp b/include/algorithms/public/Standardization.hpp index 78ad13da7..ff61996a3 100644 --- a/include/algorithms/public/Standardization.hpp +++ b/include/algorithms/public/Standardization.hpp @@ -53,14 +53,13 @@ class Standardization { using namespace Eigen; using namespace _impl; - ArrayXd input = asEigen(in); - ArrayXd result; + FluidEigenMap input = asEigen(in); + FluidEigenMap result = asEigen(out); if (!inverse) { result = (input - mMean) / mStd; } else { result = (input * mStd) + mMean; } - out <<= asFluid(result); } void process(const RealMatrixView in, RealMatrixView out, diff --git a/include/algorithms/public/TransientExtraction.hpp b/include/algorithms/public/TransientExtraction.hpp index 8bb39fb02..690b97c2a 100644 --- a/include/algorithms/public/TransientExtraction.hpp +++ b/include/algorithms/public/TransientExtraction.hpp @@ -13,6 +13,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/ARModel.hpp" #include "../util/FluidEigenMappings.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include #include @@ -31,15 +32,27 @@ class TransientExtraction using VectorXd = Eigen::VectorXd; public: + + TransientExtraction(index maxOrder, index maxBlockSize, index maxPadSize, + Allocator& alloc = FluidDefaultAllocator()) + : mModel(maxOrder, maxBlockSize + (2 * maxPadSize), alloc), + mInput(asUnsigned(maxBlockSize + (2 * maxPadSize) + maxOrder), alloc), + mDetect(asUnsigned(maxBlockSize), alloc), + mForwardError(asUnsigned(maxBlockSize + maxOrder), alloc), + mBackwardError(asUnsigned(maxBlockSize + maxOrder), alloc), + mForwardWindowedError(asUnsigned(maxBlockSize), alloc), + mBackwardWindowedError(asUnsigned(maxBlockSize), alloc) + {} + void init(index order, index blockSize, index padSize) { - mModel = ARModel(order); + mModel.init(order); prepareStream(blockSize, padSize); mInitialized = true; } void setDetectionParameters(double power, double threshHi, double threshLo, - index halfWindow = 7, index hold = 25) + index halfWindow = 7, index hold = 25) { mDetectPowerFactor = power; mDetectThreshHi = threshHi; @@ -80,36 +93,46 @@ class TransientExtraction return mBackwardWindowedError.data(); } - index detect(const double* input, index inSize) + index detect(const double* input, index inSize, + Allocator& alloc = FluidDefaultAllocator()) { + assert(mInitialized && "Call init() before processing"); + assert(modelOrder() >= mDetectHalfWindow && + "model order needs to be >= half filter size"); frame(input, inSize); - analyze(); + analyze(alloc); detection(); return mCount; } void process(const RealVectorView input, RealVectorView transients, - RealVectorView residual) + RealVectorView residual, Allocator& alloc = FluidDefaultAllocator()) { - assert(mInitialized); - index inSize = input.extent(0); + assert(mInitialized && "Call init() before processing"); + assert(modelOrder() >= mDetectHalfWindow && + "model order needs to be >= half filter size"); + index inSize = input.size(); frame(input.data(), inSize); - analyze(); + analyze(alloc); detection(); - interpolate(transients.data(), residual.data()); + interpolate(transients.data(), residual.data(), alloc); } void process(const RealVectorView input, const RealVectorView unknowns, - RealVectorView transients, RealVectorView residual) + RealVectorView transients, RealVectorView residual, + Allocator& alloc = FluidDefaultAllocator()) { - index inSize = input.extent(0); + assert(mInitialized && "Call init() before processing"); + assert(modelOrder() >= mDetectHalfWindow && + "model order needs to be >= half filter size"); + index inSize = input.size(); std::copy(unknowns.data(), unknowns.data() + hopSize(), mDetect.data()); mCount = 0; for (index i = 0, size = hopSize(); i < size; i++) if (mDetect[asUnsigned(i)] != 0) mCount++; frame(input.data(), inSize); - if (mCount) analyze(); - interpolate(transients.data(), residual.data()); + if (mCount) analyze(alloc); + interpolate(transients.data(), residual.data(), alloc); } bool initialized() { return mInitialized; } @@ -120,37 +143,49 @@ class TransientExtraction using namespace std; inSize = std::min(inSize, inputSize()); copy(mInput.data() + hopSize(), - mInput.data() + modelOrder() + padSize() + blockSize(), mInput.data()); + mInput.data() + modelOrder() + padSize() + blockSize(), mInput.data()); copy(input, input + inSize, - mInput.data() + modelOrder() + padSize() + modelOrder()); + mInput.data() + modelOrder() + padSize() + modelOrder()); fill(mInput.data() + modelOrder() + padSize() + modelOrder() + inSize, - mInput.data() + modelOrder() + analysisSize(), 0.0); + mInput.data() + modelOrder() + analysisSize(), 0.0); } - void analyze() + void analyze(Allocator& alloc) { mModel.setMinVariance(0.0000001); - mModel.estimate(mInput.data() + modelOrder(), analysisSize()); + mModel.estimate(FluidTensorView( + mInput.data(), modelOrder(), analysisSize()), + 3, 3.0, alloc); } void detection() { - const double* input = mInput.data() + modelOrder() + padSize(); + // const double* input = mInput.data() + modelOrder() + padSize(); // Forward and backward error const double normFactor = 1.0 / sqrt(mModel.variance()); + + auto inputView = FluidTensorView(mInput.data(), + modelOrder() + padSize(), + blockSize() + modelOrder() + mDetectHalfWindow); + auto fwdError = FluidTensorView( + mForwardError.data(), 0, blockSize() + mDetectHalfWindow); + auto backError = FluidTensorView( + mBackwardError.data(), 0, blockSize() + mDetectHalfWindow); + errorCalculation<&ARModel::forwardErrorArray>( - mForwardError.data(), input, blockSize() + mDetectHalfWindow + 1, - normFactor); + inputView, fwdError, normFactor); errorCalculation<&ARModel::backwardErrorArray>( - mBackwardError.data(), input, blockSize() + mDetectHalfWindow + 1, - normFactor); + inputView, backError, normFactor); // Window error functions (brute force convolution) - windowError(mForwardWindowedError.data(), - mForwardError.data() + modelOrder(), hopSize()); - windowError(mBackwardWindowedError.data(), - mBackwardError.data() + modelOrder(), hopSize()); + auto fwdWindowedError = + FluidTensorView(mForwardWindowedError.data(), 0, hopSize()); + auto backWindowedError = + FluidTensorView(mBackwardWindowedError.data(), 0, hopSize()); + + windowError(fwdError(Slice(modelOrder())), fwdWindowedError); + windowError(backError(Slice(modelOrder())), backWindowedError); // Detection index count = 0; @@ -164,7 +199,9 @@ class TransientExtraction { if (!click && (mBackwardWindowedError[asUnsigned(i)] > loThresh) && (mForwardWindowedError[asUnsigned(i)] > hiThresh)) - { click = true; } + { + click = true; + } else if (click && (mBackwardWindowedError[asUnsigned(i)] < loThresh)) { click = false; @@ -230,7 +267,7 @@ class TransientExtraction return Method(view); } - void interpolate(double* transients, double* residual) + void interpolate(double* transients, double* residual, Allocator& alloc) { const double* input = mInput.data() + padSize() + modelOrder(); const double* parameters = mModel.getParameters(); @@ -245,10 +282,15 @@ class TransientExtraction } // Declare matrices - MatrixXd A = MatrixXd::Zero(size - order, size); - MatrixXd U = MatrixXd::Zero(size, mCount); - MatrixXd K = MatrixXd::Zero(size, size - mCount); - VectorXd xK(size - mCount); + ScopedEigenMap A(size - order, size, alloc); + ScopedEigenMap U(size, mCount, alloc); + ScopedEigenMap K(size, size - mCount, alloc); + ScopedEigenMap xK(size - mCount, alloc); + + A.setZero(); + U.setZero(); + K.setZero(); + xK.setZero(); // Form data for (index i = 0; i < size - order; i++) @@ -271,9 +313,14 @@ class TransientExtraction } // Solve - MatrixXd Au = A * U; - MatrixXd M = -(Au.transpose() * Au); - MatrixXd u = M.fullPivLu().solve(Au.transpose() * (A * K) * xK); + ScopedEigenMap Au(size - order, mCount, alloc); + Au = A * U; + + ScopedEigenMap M(mCount, mCount, alloc); + M = -(Au.transpose() * Au); + + ScopedEigenMap u(mCount, alloc); + u = M.fullPivLu().solve(Au.transpose() * (A * K) * xK); // Write the output for (index i = 0, uCount = 0; i < (size - order); i++) @@ -284,27 +331,30 @@ class TransientExtraction residual[i] = input[i + order]; } - if (mRefine) refine(residual, size, Au, u); + // if (mRefine) refine(FluidTensorView(residual, 0, size), Au, + // u); for (index i = 0; i < (size - order); i++) transients[i] = input[i + order] - residual[i]; // Copy the residual indexo the correct place std::copy(residual, residual + (size - order), - mInput.data() + padSize() + order + order); + mInput.data() + padSize() + order + order); } - void refine(double* io, index size, Eigen::MatrixXd& Au, Eigen::MatrixXd& ls) + void refine( + FluidTensorView io, Eigen::MatrixXd& Au, Eigen::MatrixXd& ls) { const double energy = mModel.variance() * mCount; double energyLS = 0.0; index order = modelOrder(); + index size = io.size(); for (index i = 0; i < (size - order); i++) { if (mDetect[asUnsigned(i)] != 0) { - const double error = mModel.forwardError(io + i); + const double error = mModel.forwardError(io(Slice(i))); energyLS += error * error; } } @@ -312,8 +362,8 @@ class TransientExtraction if (energyLS < energy) { // Create the square matrix and solve - Eigen::LLT M(Au.transpose() * - Au); // Cholesky decomposition + Eigen::LLT M( + Au.transpose() * Au); // Cholesky decomposition Eigen::VectorXd u(mCount); @@ -322,7 +372,7 @@ class TransientExtraction // Write the output for (index i = 0, uCount = 0; i < (size - order); i++) { - if (mDetect[asUnsigned(i)] != 0) io[asUnsigned(i)] = u(uCount++); + if (mDetect[asUnsigned(i)] != 0) io[i] = u(uCount++); } } } @@ -341,44 +391,50 @@ class TransientExtraction return sum; } - template - void errorCalculation(double* error, const double* input, index size, - double normFactor) + template , FluidTensorView)> + void errorCalculation(FluidTensorView input, + FluidTensorView error, double normFactor) { - (mModel.*Method)(error, input, size); + (mModel.*Method)(input, error); // Take absolutes and normalise - for (index i = 0; i < size; i++) + for (index i = 0; i < error.size(); i++) error[i] = std::fabs(error[i]) * normFactor; } // Triangle window double calcWindow(double norm) { return std::min(norm, 1.0 - norm); } - void windowError(double* errorWindowed, const double* error, index size) + void windowError(FluidTensorView error, + FluidTensorView errorWindowed) { - const index windowSize = mDetectHalfWindow * 2 + 1; - const index windowOffset = mDetectHalfWindow; + assert(error.descriptor().start >= mDetectHalfWindow && + "insufficient offset for filter size"); + assert(error.size() > errorWindowed.size() - 1 + (mDetectHalfWindow) && + "insufficient input for filter size"); + + const index windowSize = mDetectHalfWindow * 2 + 3; + const index windowOffset = mDetectHalfWindow + 1; const double powFactor = mDetectPowerFactor; // Calculate window normalisation factor double windowNormFactor = 0.0; for (index j = 0; j < windowSize; j++) - windowNormFactor += calcWindow((double) j / windowSize); + windowNormFactor += calcWindow((double) j / (windowSize - 1)); windowNormFactor = 1.0 / windowNormFactor; // Do window processing - for (index i = 0; i < size; i++) + for (index i = 0; i < errorWindowed.size(); i++) { double windowed = 0.0; - for (index j = 1; j < windowSize; j++) + for (index j = 1; j < windowSize - 1; j++) { const double value = pow(fabs(error[i - windowOffset + j]), powFactor); - windowed += value * calcWindow((double) j / windowSize); - ; + windowed += value * calcWindow((double) j / (windowSize - 1)); } errorWindowed[i] = pow((windowed * windowNormFactor), 1.0 / powFactor); @@ -395,27 +451,27 @@ class TransientExtraction mBackwardWindowedError.resize(asUnsigned(hopSize()), 0.0); } - ARModel mModel{20}; + ARModel mModel; std::mt19937_64 mRandomGenerator{std::random_device()()}; - index mBlockSize{0}; - index mPadSize{0}; - index mCount{0}; - bool mRefine{false}; + index mBlockSize{0}; + index mPadSize{0}; + index mCount{0}; + // bool mRefine{false}; index mDetectHalfWindow{1}; index mDetectHold{25}; double mDetectPowerFactor{1.4}; double mDetectThreshHi{1.5}; double mDetectThreshLo{3.0}; - std::vector mInput; - std::vector mDetect; - std::vector mForwardError; - std::vector mBackwardError; - std::vector mForwardWindowedError; - std::vector mBackwardWindowedError; - bool mInitialized{false}; + rt::vector mInput; + rt::vector mDetect; + rt::vector mForwardError; + rt::vector mBackwardError; + rt::vector mForwardWindowedError; + rt::vector mBackwardWindowedError; + bool mInitialized{false}; }; } // namespace algorithm diff --git a/include/algorithms/public/TransientSegmentation.hpp b/include/algorithms/public/TransientSegmentation.hpp index e17561685..a15e8f9e6 100644 --- a/include/algorithms/public/TransientSegmentation.hpp +++ b/include/algorithms/public/TransientSegmentation.hpp @@ -12,6 +12,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "TransientExtraction.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include @@ -22,12 +23,13 @@ class TransientSegmentation : public algorithm::TransientExtraction { public: + using TransientExtraction::TransientExtraction; + void setDetectionParameters(double power, double threshHi, double threshLo, - index halfWindow = 7, index hold = 25, - index minSegment = 50) + index halfWindow = 7, index hold = 25, index minSegment = 50) { - TransientExtraction::setDetectionParameters(power, threshHi, threshLo, - halfWindow, hold); + TransientExtraction::setDetectionParameters( + power, threshHi, threshLo, halfWindow, hold); mMinSegment = minSegment; } @@ -38,9 +40,10 @@ class TransientSegmentation : public algorithm::TransientExtraction mDebounce = 0; } - void process(const RealVectorView input, RealVectorView output) + void process(const RealVectorView input, RealVectorView output, + Allocator& alloc = FluidDefaultAllocator()) { - detect(input.data(), input.extent(0)); + detect(input.data(), input.extent(0), alloc); const double* transientDetection = getDetect(); for (index i = 0; i < std::min(hopSize(), output.size()); i++) { diff --git a/include/algorithms/public/UMAP.hpp b/include/algorithms/public/UMAP.hpp index c9baa3aca..04d5feb68 100644 --- a/include/algorithms/public/UMAP.hpp +++ b/include/algorithms/public/UMAP.hpp @@ -4,8 +4,9 @@ #include "../util/FluidEigenMappings.hpp" #include "../util/SpectralEmbedding.hpp" #include "../../data/TensorTypes.hpp" +#include "../../data/FluidMemory.hpp" #include -#include +#include #include #include #include @@ -158,38 +159,54 @@ class UMAP } - void transformPoint(RealVectorView in, RealVectorView out) const + void transformPoint(RealVectorView in, RealVectorView out, + Allocator& alloc = FluidDefaultAllocator()) const { if (!mInitialized) return; - SparseMatrixXd knnGraph(1, mEmbedding.rows()); - ArrayXXd dists = ArrayXXd::Zero(1, mK); - knnGraph.reserve(mK); - auto nearest = mTree.kNearest(in, mK); - auto nearestIds = nearest.getIds(); - auto distances = nearest.getData().col(0); - for (index j = 0; j < mK; j++) + + auto [distances, nearestIds] = mTree.kNearest(in, mK, 0, alloc); + using SparseMap = Eigen::Map>; + + rt::vector data(alloc); + rt::vector inner(alloc); + data.reserve(asUnsigned(mK)); + inner.reserve(asUnsigned(mK)); + + rt::vector outer(2, 0, alloc); + outer[1] = static_cast(mK); + + ScopedEigenMap dists(1, mK, alloc); + for (size_t j = 0; j < asUnsigned(mK); j++) { - index neighborIndex = stoi(nearestIds(j)); - dists(0, j) = distances(j); - knnGraph.insert(0, neighborIndex) = distances(j); + int neighborIndex = stoi(*nearestIds[j]); + dists(0, asSigned(j)) = distances[j]; + data.push_back(distances[j]); + inner.push_back(neighborIndex); } - knnGraph.makeCompressed(); - ArrayXd sigma = findSigma(mK, dists); + int maxIndex = *std::max_element(inner.cbegin(), inner.cend()); + SparseMap knnGraph(1, maxIndex, mK, outer.data(), inner.data(), + data.data()); + ScopedEigenMap sigma = findSigma(mK, dists, 64, 1e-5, alloc); computeHighDimProb(dists, sigma, knnGraph); - normalizeRows(knnGraph); - ArrayXXd embedding = initTransformEmbedding(knnGraph, mEmbedding, 1); - ArrayXd result = embedding.row(0); - out <<= _impl::asFluid(result); + normalizeRows(knnGraph, alloc); + ScopedEigenMap embedding = + initTransformEmbedding(knnGraph, mEmbedding, 1, alloc); + _impl::asEigen(out) = embedding.row(0).transpose(); } private: - template - void traverseGraph(const SparseMatrixXd& graph, F func) const + template + void traverseGraph(Eigen::SparseCompressedBase& graph, F func) const { for (index i = 0; i < graph.outerSize(); i++) { - for (SparseMatrixXd::InnerIterator it(graph, i); it; ++it) { func(it); } + for (typename Eigen::SparseCompressedBase::InnerIterator it( + graph, i); + it; ++it) + { + func(it); + } } } @@ -203,12 +220,15 @@ class UMAP return CE.sum(); } - ArrayXd findSigma(index k, Ref dists, index maxIter = 64, - double tolerance = 1e-5) const + ScopedEigenMap + findSigma(index k, Ref dists, index maxIter = 64, + double tolerance = 1e-5, + Allocator& alloc = FluidDefaultAllocator()) const { using namespace std; - double target = log2(k); - ArrayXd result = ArrayXd::Zero(dists.rows()); + double target = log2(k); + ScopedEigenMap result(dists.rows(), alloc); + result.setZero(); for (index i = 0; i < dists.rows(); i++) { index iter = maxIter; @@ -241,8 +261,9 @@ class UMAP return result; } + template void computeHighDimProb(const Ref& dists, const Ref& sigma, - SparseMatrixXd& graph) const + Eigen::SparseCompressedBase& graph) const { traverseGraph(graph, [&](auto it) { it.valueRef() = @@ -269,15 +290,15 @@ class UMAP auto data = in.getData(); for (index i = 0; i < in.size(); i++) { - auto nearest = mTree.kNearest(data.row(i), discardFirst ? k + 1 : k); - auto nearestIds = nearest.getIds(); - auto distances = nearest.getData().col(0); - for (index j = 0; j < k; j++) + auto [distances, nearestIds] = + mTree.kNearest(data.row(i), discardFirst ? k + 1 : k); + + for (size_t j = 0; j < asUnsigned(k); j++) { - index pos = discardFirst ? j + 1 : j; - index neighborIndex = stoi(nearestIds(pos)); - dists(i, j) = distances(pos); - graph.insert(i, neighborIndex) = distances(pos); + size_t pos = discardFirst ? j + 1 : j; + index neighborIndex = stoi(*nearestIds[pos]); + dists(i, asSigned(j)) = distances[pos]; + graph.insert(i, neighborIndex) = distances[pos]; } } } @@ -297,8 +318,9 @@ class UMAP return 10.0 * result; } - void getGraphIndices(const SparseMatrixXd& graph, Ref rowIndices, - Ref colIndices) const + template + void getGraphIndices(Eigen::SparseCompressedBase& graph, + Ref rowIndices, Ref colIndices) const { index p = 0; traverseGraph(graph, [&](auto it) { @@ -308,8 +330,9 @@ class UMAP }); } - void computeEpochsPerSample(const SparseMatrixXd& graph, - Ref epochsPerSample) const + template + void computeEpochsPerSample(Eigen::SparseCompressedBase& graph, + Ref epochsPerSample) const { index p = 0; double maxVal = graph.coeffs().maxCoeff(); @@ -384,19 +407,25 @@ class UMAP } } - ArrayXXd initTransformEmbedding(const SparseMatrixXd& graph, - Ref reference, index N) const + template + ScopedEigenMap + initTransformEmbedding(Eigen::SparseCompressedBase& graph, + Ref reference, index N, + Allocator& alloc = FluidDefaultAllocator()) const { - ArrayXXd embedding = ArrayXXd::Zero(N, reference.cols()); + ScopedEigenMap embedding(N, reference.cols(), alloc); + embedding.setZero(); // todo: sort out 2D expression constructor? traverseGraph(graph, [&](auto it) { embedding.row(it.row()) += (reference.row(it.col()) * it.value()); }); return embedding; } - void normalizeRows(const SparseMatrixXd& graph) const + template + void normalizeRows(Eigen::SparseCompressedBase& graph, + Allocator& alloc = FluidDefaultAllocator()) const { - ArrayXd sums = ArrayXd::Zero(graph.innerSize()); + ScopedEigenMap sums(ArrayXd::Zero(graph.innerSize()), alloc); traverseGraph(graph, [&](auto it) { sums(it.row()) += it.value(); }); traverseGraph( graph, [&](auto it) { it.valueRef() = it.value() / sums(it.row()); }); diff --git a/include/algorithms/public/YINFFT.hpp b/include/algorithms/public/YINFFT.hpp index a241ceb1f..d12014ac9 100644 --- a/include/algorithms/public/YINFFT.hpp +++ b/include/algorithms/public/YINFFT.hpp @@ -14,6 +14,7 @@ under the European Union’s Horizon 2020 research and innovation programme #include "../util/FluidEigenMappings.hpp" #include "../util/PeakDetection.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include "../../data/TensorTypes.hpp" #include @@ -25,20 +26,27 @@ class YINFFT public: void processFrame(const RealVectorView& input, RealVectorView output, - double minFreq, double maxFreq, double sampleRate) + double minFreq, double maxFreq, double sampleRate, + Allocator& alloc = FluidDefaultAllocator()) { using namespace Eigen; PeakDetection pd; - ArrayXd mag = _impl::asEigen(input); - ArrayXd squareMag = mag.square(); - index nBins = mag.size(); - FFT fft(2 * (mag.size() - 1)); - double squareMagSum = 2 * squareMag.sum(); - ArrayXd squareMagSym(2 * (nBins - 1)); + // ScopedEigenMap mag = _; + ScopedEigenMap squareMag(input.size(), alloc); + squareMag = _impl::asEigen(input).square(); + + index nBins = input.size(); + FFT fft(2 * (input.size() - 1)); + double squareMagSum = 2 * squareMag.sum(); + + ScopedEigenMap squareMagSym(2 * (nBins - 1), alloc); squareMagSym << squareMag[0], squareMag.segment(1, nBins - 1), squareMag.segment(1, nBins - 2).reverse(); - ArrayXcd squareMagFFT = fft.process(squareMagSym); - ArrayXd yin = squareMagSum - squareMagFFT.real(); + + Eigen::Map squareMagFFT = fft.process(squareMagSym); + ScopedEigenMap yin(squareMagFFT.size(), alloc); + yin = squareMagSum - squareMagFFT.real(); + if (maxFreq == 0) maxFreq = 1; if (minFreq == 0) minFreq = 1; yin(0) = 1; @@ -52,7 +60,8 @@ class YINFFT double pitchConfidence = 0; if (tmpSum > 0) { - ArrayXd yinFlip = -yin; + ScopedEigenMap yinFlip(yin.size(), alloc); + yinFlip = -yin; // segment from max to min freq index minBin = std::lrint(sampleRate / maxFreq); index maxBin = std::lrint(sampleRate / minFreq); @@ -61,9 +70,9 @@ class YINFFT maxBin = yinFlip.size() - minBin - 1; if (maxBin > minBin) { - yinFlip = yinFlip.segment(minBin, maxBin - minBin).eval(); - - auto vec = pd.process(yinFlip, 1, yinFlip.minCoeff()); + // yinFlip = yinFlip.segment(minBin, maxBin - minBin).eval(); + auto yinSeg = yinFlip.segment(minBin, maxBin - minBin); + auto vec = pd.process(yinSeg, 1, yinSeg.minCoeff(), true, true, alloc); if (vec.size() > 0) { pitch = sampleRate / (minBin + vec[0].first); diff --git a/include/algorithms/util/ARModel.hpp b/include/algorithms/util/ARModel.hpp index 76296b350..7db4840fe 100644 --- a/include/algorithms/util/ARModel.hpp +++ b/include/algorithms/util/ARModel.hpp @@ -10,10 +10,13 @@ under the European Union’s Horizon 2020 research and innovation programme #pragma once -#include "ConvolutionTools.hpp" +#include "FluidEigenMappings.hpp" +#include "FFT.hpp" #include "Toeplitz.hpp" #include "../public/WindowFuncs.hpp" #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" +#include "../../data/FluidTensor.hpp" #include #include #include @@ -28,118 +31,182 @@ class ARModel using MatrixXd = Eigen::MatrixXd; using ArrayXd = Eigen::ArrayXd; using VectorXd = Eigen::VectorXd; + using ArrayXcd = Eigen::ArrayXcd; public: - ARModel(index order) : mParameters(VectorXd::Zero(order)) {} + ARModel(index maxOrder, index maxBlockSize, Allocator& alloc) + : mParameters(maxOrder, alloc), + mWindow(maxBlockSize + maxOrder, alloc), + mSpectralIn(nextPower2(2 * maxBlockSize), alloc), + mSpectralOut(nextPower2(2 * maxBlockSize), alloc), + mFFT(nextPower2(2 * maxBlockSize), alloc), + mIFFT(nextPower2(2 * maxBlockSize), alloc) + {} + + void init(index order) { mParameters.head(order).setZero(); } + const double* getParameters() const { return mParameters.data(); } double variance() const { return mVariance; } index order() const { return mParameters.size(); } - void estimate(const double* input, index size, index nIterations = 3, - double robustFactor = 3.0) + void estimate(FluidTensorView input, index nIterations, + double robustFactor, Allocator& alloc) { if (nIterations > 0) - robustEstimate(input, size, nIterations, robustFactor); + robustEstimate(input, nIterations, robustFactor, alloc); else - directEstimate(input, size, true); + directEstimate(input, true, alloc); } - double fowardPrediction(const double* input) + double fowardPrediction(FluidTensorView input) { - return modelPredict>(input); + double prediction; + modelPredict(input, FluidTensorView(&prediction, 0, 1), + std::negate{}, Predict{}); + return prediction; } - double backwardPrediction(const double* input) + double backwardPrediction(FluidTensorView input) { - struct Identity - { - index operator()(index a) { return a; } - }; - return modelPredict(input); + double prediction; + modelPredict(input, FluidTensorView(&prediction, 0, 1), + Identity{}, Predict{}); + return prediction; } - double forwardError(const double* input) + double forwardError(FluidTensorView input) { - return modelError<&ARModel::fowardPrediction>(input); + double error; + forwardErrorArray(input, FluidTensorView(&error, 0, 1)); + return error; } - double backwardError(const double* input) + double backwardError(FluidTensorView input) { - return modelError<&ARModel::backwardPrediction>(input); + double error; + backwardErrorArray(input, FluidTensorView(&error, 0, 1)); + return error; } - void forwardErrorArray(double* errors, const double* input, index size) + void forwardErrorArray( + FluidTensorView input, FluidTensorView errors) { - modelErrorArray<&ARModel::forwardError>(errors, input, size); + modelPredict(input, errors, std::negate{}, Error{}); } - void backwardErrorArray(double* errors, const double* input, index size) + void backwardErrorArray( + FluidTensorView input, FluidTensorView errors) { - modelErrorArray<&ARModel::backwardError>(errors, input, size); + modelPredict(input, errors, Identity{}, Error{}); } void setMinVariance(double variance) { mMinVariance = variance; } private: - template - double modelPredict(const double* input) + struct Identity { - double estimate = 0.0; + index operator()(index a) { return a; } + }; - for (index i = 0; i < mParameters.size(); i++) - estimate += mParameters(i) * input[Op()(i + 1)]; - - return estimate; - } + struct Predict + { + double operator()(double, double estimate) { return estimate; } + }; - template - double modelError(const double* input) + struct Error { - return input[0] - (this->*Method)(input); - } + double operator()(double input, double estimate) + { + return input - estimate; + } + }; - template - void modelErrorArray(double* errors, const double* input, index size) + template + void modelPredict(FluidTensorView input, + FluidTensorView output, Indexer fIdx, OutputFn fOut) { - for (index i = 0; i < size; i++) errors[i] = (this->*Method)(input + i); + 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 && + "array bounds error in AR model prediction: input too short"); + assert(fIdx(mParameters.size()) < 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 && + "array bounds error in AR model prediction: input too short"); + assert(((numPredictions - 1) + fIdx(mParameters.size())) < input.size() && + "array bounds error in AR model prediction: input too short"); + + const double* input_ptr = input.data(); + + for (index p = 0; p < numPredictions; p++) + { + double estimate = 0; + for (index i = 0; i < mParameters.size(); i++) + estimate += mParameters(i) * (input_ptr + p)[fIdx(i + 1)]; + output[p] = fOut(input_ptr[p], estimate); + } } - void directEstimate(const double* input, index size, bool updateVariance) + void directEstimate(FluidTensorView input, + bool updateVariance, Allocator& alloc) { + index size = input.size(); + // copy input to a 32 byte aligned block (otherwise risk segfaults on Linux) - VectorXd frame = Eigen::Map(input, size); + // FluidEigenMap frame = + ScopedEigenMap frame(size, alloc); + frame = _impl::asEigen(input); if (mUseWindow) { if (mWindow.size() != size) { - mWindow = ArrayXd::Zero(size); - WindowFuncs::map()[WindowFuncs::WindowTypes::kHann](size, mWindow); + mWindow.setZero(); // = ArrayXd::Zero(size); + WindowFuncs::map()[WindowFuncs::WindowTypes::kHann]( + size, mWindow.head(size)); } - frame.array() *= mWindow; + frame.array() *= mWindow.head(size); } + // VectorXd autocorrelation(size); + ScopedEigenMap autocorrelation(size, alloc); + autocorrelate(frame, autocorrelation); + // algorithm::autocorrelateReal(autocorrelation.data(), frame.data(), - VectorXd autocorrelation(size); - algorithm::autocorrelateReal(autocorrelation.data(), frame.data(), - asUnsigned(size)); + // 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()); + // autocorrelation.conservativeResize(mParameters.size()); // Form a toeplitz matrix - MatrixXd mat = toeplitz(autocorrelation); + // 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) { + index idx = row - col; + if (idx < 0) idx += n; + return autocorrelation(idx); + }); // Yule Walker autocorrelation(0) = pN; std::rotate(autocorrelation.data(), autocorrelation.data() + 1, - autocorrelation.data() + mParameters.size()); - mParameters = mat.llt().solve(autocorrelation); + autocorrelation.data() + mParameters.size()); + mParameters = mat.llt().solve(autocorrelation.head(mParameters.size())); if (updateVariance) { @@ -149,76 +216,80 @@ class ARModel for (index i = 0; i < mParameters.size() - 1; i++) variance -= mParameters(i) * mat(0, i + 1); - setVariance((variance - (mParameters(mParameters.size() - 1) * pN)) / - size); + setVariance( + (variance - (mParameters(mParameters.size() - 1) * pN)) / size); } } - void robustEstimate(const double* input, index size, index nIterations, - double robustFactor) + void robustEstimate(FluidTensorView input, index nIterations, + double robustFactor, Allocator& alloc) { - std::vector estimates(asUnsigned(size + mParameters.size())); + FluidTensor estimates(input.size() + mParameters.size(), alloc); // Calculate an initial estimate of parameters - directEstimate(input, size, true); + directEstimate(input, true, alloc); + + assert(input.descriptor().start >= mParameters.size() && + "too little offset into input data"); // Initialise Estimates - for (index i = mParameters.size(); i < mParameters.size() + size; i++) - estimates[asUnsigned(i)] = input[i - mParameters.size()]; + for (index i = 0; i < input.size() + mParameters.size(); i++) + estimates[i] = input.data()[i - mParameters.size()]; // Variance - robustVariance(estimates.data() + mParameters.size(), input, size, - robustFactor); + robustVariance(estimates(Slice(mParameters.size())), input, robustFactor); // Iterate for (index iterations = nIterations; iterations--;) - robustIteration(estimates.data() + mParameters.size(), input, size, - robustFactor); + robustIteration( + estimates(Slice(mParameters.size())), input, robustFactor, alloc); } double robustResidual(double input, double prediction, double cs) { + assert(cs > 0); return cs * psiFunction((input - prediction) / cs); } - void robustVariance(double* estimates, const double* input, index size, - double robustFactor) + void robustVariance(FluidTensorView estimates, + FluidTensorView input, double robustFactor) { double residualSqSum = 0.0; // Iterate to find new filtered input - for (index i = 0; i < size; i++) + for (index i = 0; i < input.size(); i++) { const double residual = - robustResidual(input[i], fowardPrediction(estimates + i), - robustFactor * sqrt(mVariance)); + robustResidual(input[i], fowardPrediction(estimates(Slice(i))), + robustFactor * sqrt(mVariance)); residualSqSum += residual * residual; } - setVariance(residualSqSum / size); + setVariance(residualSqSum / input.size()); } - void robustIteration(double* estimates, const double* input, index size, - double robustFactor) + void robustIteration(FluidTensorView estimates, + FluidTensorView input, double robustFactor, + Allocator& alloc) { // Iterate to find new filtered input - for (index i = 0; i < size; i++) + for (index i = 0; i < input.size(); i++) { - const double prediction = fowardPrediction(estimates); - estimates[0] = - prediction + - robustResidual(input[i], prediction, robustFactor * sqrt(mVariance)); + const double prediction = fowardPrediction(estimates(Slice(i))); + estimates[i] = prediction + robustResidual(input[i], prediction, + robustFactor * sqrt(mVariance)); } // New parameters - directEstimate(estimates, size, false); - robustVariance(estimates, input, size, robustFactor); + directEstimate(estimates(Slice(0, input.size())), false, alloc); + robustVariance(estimates(Slice(0, input.size())), input, robustFactor); } void setVariance(double variance) { - if (variance > 0) variance = std::max(mMinVariance, variance); + variance = std::max(mMinVariance, variance); mVariance = variance; + assert(mVariance >= 0); } // Huber PSI function @@ -227,11 +298,41 @@ class ARModel return fabs(x) > 1 ? std::copysign(1.0, x) : x; } - VectorXd mParameters; - double mVariance{0.0}; - ArrayXd mWindow; - bool mUseWindow{true}; - double mMinVariance{0.0}; + template + void autocorrelate( + const Eigen::DenseBase& in, Eigen::DenseBase& out) + { + mSpectralIn.head(in.size()) = in; + + index fftSize = nextPower2(in.size() * 2); + mFFT.resize(fftSize); + mIFFT.resize(fftSize); + + mSpectralIn.head(in.size()) = in; + + auto spec = mFFT.process(mSpectralIn.head(in.size())); + + mSpectralOut.head(spec.size()) = spec * spec.conjugate(); + + out.head(in.size()) = + mIFFT.process(mSpectralOut.head(spec.size())).head(in.size()); + out.head(in.size()) /= (fftSize * in.size()); + } + + index nextPower2(index n) + { + return static_cast(std::pow(2, std::ceil(std::log2(n)))); + } + + ScopedEigenMap mParameters; + double mVariance{0.0}; + ScopedEigenMap mWindow; + ScopedEigenMap mSpectralIn; + ScopedEigenMap mSpectralOut; + FFT mFFT; + IFFT mIFFT; + bool mUseWindow{true}; + double mMinVariance{1e-5}; }; } // namespace algorithm diff --git a/include/algorithms/util/ConvolutionTools.hpp b/include/algorithms/util/ConvolutionTools.hpp deleted file mode 100644 index 9d9568915..000000000 --- a/include/algorithms/util/ConvolutionTools.hpp +++ /dev/null @@ -1,590 +0,0 @@ -/* -Part of the Fluid Corpus Manipulation Project (http://www.flucoma.org/) -Copyright 2017-2019 University of Huddersfield. -Licensed under the BSD-3 License. -See license.md file in the project root for full license information. -This project has received funding from the European Research Council (ERC) -under the European Union’s Horizon 2020 research and innovation programme -(grant agreement No 725899). -*/ - -#pragma once - -#include -#include -#include - -namespace fluid { -namespace algorithm { - -// The edge mode determines wraparound etc. - -enum EdgeMode { kEdgeLinear, kEdgeWrap, kEdgeWrapCentre, kEdgeFold }; - -namespace impl { // Here is the underlying implementation - // See the bottom of the file for interface - // Note that the function calls will currently allocate and - // dellocate on the heap - -// Setups - -struct FFTComplexSetup -{ - FFTComplexSetup(size_t maxFFTLog2) - { - hisstools_create_setup(&mSetup, maxFFTLog2); - } - ~FFTComplexSetup() { hisstools_destroy_setup(mSetup); } - - FFTComplexSetup(const FFTComplexSetup&) = delete; - FFTComplexSetup operator=(const FFTComplexSetup&) = delete; - - FFT_SETUP_D mSetup; -}; - -struct FFTRealSetup : public FFTComplexSetup -{ - FFTRealSetup(size_t maxFFTLog2) : FFTComplexSetup(maxFFTLog2){}; -}; - -// Temporary Memory - -struct TempSpectra -{ - TempSpectra(size_t dataSize) - { - mData = allocate_aligned(dataSize * 2); - mSpectra.realp = mData; - mSpectra.imagp = mSpectra.realp + dataSize; - } - - ~TempSpectra() { deallocate_aligned(mData); } - - FFT_SPLIT_COMPLEX_D mSpectra; - double* mData; -}; - -struct ConvolveOp -{ - template - void operator()(T& outR, T& outI, const T a, const T b, const T c, const T d, - T scale) - { - outR = scale * (a * c - b * d); - outI = scale * (a * d + b * c); - } -}; - -struct CorrelateOp -{ - template - void operator()(T& outR, T& outI, const T a, const T b, const T c, const T d, - T scale) - { - outR = scale * (a * c + b * d); - outI = scale * (b * c - a * d); - } -}; - -// Complex Forward Transform - -void transformForward(FFTComplexSetup& setup, FFT_SPLIT_COMPLEX_D& io, - size_t fftSizelog2) -{ - hisstools_fft(setup.mSetup, &io, fftSizelog2); -} - -// Complex Inverse Transform - -void transformInverse(FFTComplexSetup& setup, FFT_SPLIT_COMPLEX_D& io, - size_t fftSizelog2) -{ - hisstools_ifft(setup.mSetup, &io, fftSizelog2); -} - -// Real Forward Transform (in-place) - -void transformForwardReal(FFTRealSetup& setup, FFT_SPLIT_COMPLEX_D& io, - size_t fftSizelog2) -{ - hisstools_rfft(setup.mSetup, &io, fftSizelog2); -} - -// Real Forward Tansform (with unzipping) - -void transformForwardReal(FFTRealSetup& setup, FFT_SPLIT_COMPLEX_D& output, - const double* input, size_t size, size_t fftSizelog2) -{ - hisstools_rfft(setup.mSetup, input, &output, size, fftSizelog2); -} - -// Real Inverse Transform (in-place) - -void transformInverseReal(FFTRealSetup& setup, FFT_SPLIT_COMPLEX_D& io, - size_t fftSizelog2) -{ - hisstools_rifft(setup.mSetup, &io, fftSizelog2); -} - -// Real Inverse Transform (with zipping) - -void transformInverseReal(FFTRealSetup& setup, double* output, - FFT_SPLIT_COMPLEX_D& input, size_t fftSizelog2) -{ - hisstools_rifft(setup.mSetup, &input, output, fftSizelog2); -} - -// Size calculations - -size_t ilog2(size_t value) -{ - size_t bitShift = value; - size_t bitCount = 0; - - while (bitShift) - { - bitShift >>= 1U; - bitCount++; - } - - if (bitCount && value == 1U << (bitCount - 1U)) - return bitCount - 1U; - else - return bitCount; -} - -size_t calcLinearSize(size_t size1, size_t size2) -{ - return (size1 && size2) ? size1 + size2 - 1 : 0; -} - -size_t calcSize(size_t size1, size_t size2, EdgeMode mode) -{ - size_t linearSize = size1 + size2 - 1; - size_t sizeOut = mode != kEdgeLinear ? std::max(size1, size2) : linearSize; - - if (mode == kEdgeFold && !(std::min(size1, size2) & 1U)) sizeOut++; - - return (!size1 || !size2) ? 0 : sizeOut; -} - -// Arrangement - -// Memory manipulation (complex) - -void copyZero(double* output, const double* input, size_t inSize, - size_t outSize) -{ - std::copy(input, input + inSize, output); - std::fill_n(output + inSize, outSize - inSize, 0.0); -} - -void copy(FFT_SPLIT_COMPLEX_D output, const FFT_SPLIT_COMPLEX_D spectrum, - size_t outOffset, size_t offset, size_t size) -{ - std::copy(spectrum.realp + offset, spectrum.realp + size + offset, - output.realp + outOffset); - std::copy(spectrum.imagp + offset, spectrum.imagp + size + offset, - output.imagp + outOffset); -} - -void wrap(FFT_SPLIT_COMPLEX_D output, const FFT_SPLIT_COMPLEX_D spectrum, - size_t outOffset, size_t offset, size_t size) -{ - for (size_t i = 0; i < size; i++) - { - output.realp[i + outOffset] += spectrum.realp[i + offset]; - output.imagp[i + outOffset] += spectrum.imagp[i + offset]; - } -} - -void fold(FFT_SPLIT_COMPLEX_D output, const FFT_SPLIT_COMPLEX_D spectrum, - size_t outOffset, size_t endOffset, size_t size) -{ - for (size_t i = 1; i <= size; i++) - { - output.realp[i + outOffset - 1] += spectrum.realp[endOffset - i]; - output.imagp[i + outOffset - 1] += spectrum.imagp[endOffset - i]; - } -} - -// Memory manipulation (real) - -struct Assign -{ - void operator()(double* output, const double* ptr) { *output = *ptr; } -}; -struct Accum -{ - void operator()(double* output, const double* ptr) { *output += *ptr; } -}; -struct Increment -{ - template - T* operator()(T*& ptr) - { - return ptr++; - } -}; -struct Decrement -{ - template - T* operator()(T*& ptr) - { - return ptr--; - } -}; - -template -void zip(double* output, const double* p1, const double* p2, size_t size, Op op, - PtrOp pOp) -{ - for (size_t i = 0; i < (size >> 1); i++) - { - op(pOp(output), pOp(p1)); - op(pOp(output), pOp(p2)); - } - - if (size & 1U) op(pOp(output), pOp(p1)); -} - -void copy(double* output, const FFT_SPLIT_COMPLEX_D spectrum, size_t outOffset, - size_t offset, size_t size) -{ - const double* p1 = (offset & 1U) ? spectrum.imagp + (offset >> 1) - : spectrum.realp + (offset >> 1); - const double* p2 = (offset & 1U) ? spectrum.realp + (offset >> 1) + 1 - : spectrum.imagp + (offset >> 1); - - zip(output + outOffset, p1, p2, size, Assign(), Increment()); -} - -void wrap(double* output, const FFT_SPLIT_COMPLEX_D spectrum, size_t outOffset, - size_t offset, size_t size) -{ - const double* p1 = (offset & 1U) ? spectrum.imagp + (offset >> 1) - : spectrum.realp + (offset >> 1); - const double* p2 = (offset & 1U) ? spectrum.realp + (offset >> 1) + 1 - : spectrum.imagp + (offset >> 1); - - zip(output + outOffset, p1, p2, size, Accum(), Increment()); -} - -void fold(double* output, const FFT_SPLIT_COMPLEX_D spectrum, size_t outOffset, - size_t endOffset, size_t size) -{ - const double* p1 = (endOffset & 1U) ? spectrum.realp + (endOffset >> 1) - : spectrum.imagp + (endOffset >> 1) - 1; - const double* p2 = (endOffset & 1U) ? spectrum.imagp + (endOffset >> 1) - 1 - : spectrum.realp + (endOffset >> 1) - 1; - - zip(output + outOffset, p1, p2, size, Accum(), Decrement()); -} - -template -void arrangeOutput(T output, FFT_SPLIT_COMPLEX_D spectrum, size_t minSize, - size_t sizeOut, size_t linearSize, size_t /*fftSize*/, - EdgeMode mode, ConvolveOp) -{ - size_t offset = - (mode == kEdgeFold || mode == kEdgeWrapCentre) ? (minSize - 1) / 2 : 0; - - copy(output, spectrum, 0, offset, sizeOut); - - if (mode == kEdgeWrap) - wrap(output, spectrum, 0, sizeOut, linearSize - sizeOut); - - if (mode == kEdgeWrapCentre) - { - size_t endWrap = (minSize - 1) - offset; - - wrap(output, spectrum, 0, linearSize - endWrap, endWrap); - wrap(output, spectrum, sizeOut - offset, 0, offset); - } - - if (mode == kEdgeFold) - { - size_t foldSize = minSize / 2; - size_t foldOffset = sizeOut - foldSize; - - fold(output, spectrum, 0, foldSize, foldSize); - fold(output, spectrum, foldOffset, linearSize, foldSize); - } -} - -template -void arrangeOutput(T output, FFT_SPLIT_COMPLEX_D spectrum, size_t minSize, - size_t sizeOut, size_t linearSize, size_t fftSize, - EdgeMode mode, CorrelateOp) -{ - size_t maxSize = (linearSize - minSize) + 1; - - if (mode == kEdgeLinear || mode == kEdgeWrap) - { - size_t extraSize = minSize - 1; - - copy(output, spectrum, 0, 0, maxSize); - - if (mode == kEdgeLinear) - copy(output, spectrum, maxSize, fftSize - extraSize, extraSize); - else - wrap(output, spectrum, (linearSize - (2 * (minSize - 1))), - fftSize - extraSize, extraSize); - } - else - { - size_t offset = minSize / 2; - - copy(output, spectrum, 0, fftSize - offset, offset); - copy(output, spectrum, offset, 0, sizeOut - offset); - - if (mode == kEdgeWrapCentre) - { - size_t endWrap = minSize - offset - 1; - - wrap(output, spectrum, 0, maxSize - offset, offset); - wrap(output, spectrum, sizeOut - endWrap, fftSize - (minSize - 1), - endWrap); - } - else - { - fold(output, spectrum, 0, fftSize - (offset - 1), offset); - fold(output, spectrum, sizeOut - offset, maxSize, offset); - } - } -} - -// Calculations - -template -void binaryOp(FFT_SPLIT_COMPLEX_D& io1, FFT_SPLIT_COMPLEX_D& in2, - size_t dataLength, double scale, Op op) -{ - const int vecSize = SIMDLimits::max_size; - - if (dataLength == 1 || dataLength < (vecSize / 2)) - { - op(io1.realp[0], io1.imagp[0], io1.realp[0], io1.imagp[0], in2.realp[0], - in2.imagp[0], scale); - } - else if (dataLength < vecSize) - { - const int currentVecSize = SIMDLimits::max_size / 2; - - SIMDType* real1 = - reinterpret_cast*>(io1.realp); - SIMDType* imag1 = - reinterpret_cast*>(io1.imagp); - SIMDType* real2 = - reinterpret_cast*>(in2.realp); - SIMDType* imag2 = - reinterpret_cast*>(in2.imagp); - - SIMDType scaleVec(scale); - - for (size_t i = 0; i < (dataLength / currentVecSize); i++) - op(real1[i], imag1[i], real1[i], imag1[i], real2[i], imag2[i], scaleVec); - } - else - { - SIMDType* real1 = - reinterpret_cast*>(io1.realp); - SIMDType* imag1 = - reinterpret_cast*>(io1.imagp); - SIMDType* real2 = - reinterpret_cast*>(in2.realp); - SIMDType* imag2 = - reinterpret_cast*>(in2.imagp); - - SIMDType scaleVec(scale); - - for (size_t i = 0; i < (dataLength / vecSize); i++) - op(real1[i], imag1[i], real1[i], imag1[i], real2[i], imag2[i], scaleVec); - } -} - -template -void binaryOpReal(FFT_SPLIT_COMPLEX_D& spectrum1, - FFT_SPLIT_COMPLEX_D& spectrum2, size_t dataLength, - double scale, Op op) -{ - // Store DC and Nyquist Results - - const double DC = spectrum1.realp[0] * spectrum2.realp[0] * scale; - const double nyquist = spectrum1.imagp[0] * spectrum2.imagp[0] * scale; - - binaryOp(spectrum1, spectrum2, dataLength, scale, op); - - // Set DC and Nyquist bins - - spectrum1.realp[0] = DC; - spectrum1.imagp[0] = nyquist; -} - -template -void binarySpectralOperation(double* rOut, double* iOut, const double* rIn1, - size_t sizeR1, const double* iIn1, size_t sizeI1, - const double* rIn2, size_t sizeR2, - const double* iIn2, size_t sizeI2, EdgeMode mode, - Op op) -{ - size_t size1 = std::max(sizeR1, sizeI1); - size_t size2 = std::max(sizeR2, sizeI2); - - size_t linearSize = calcLinearSize(size1, size2); - size_t sizeOut = calcSize(size1, size2, mode); - size_t fftSizelog2 = ilog2(linearSize); - size_t fftSize = 1 << fftSizelog2; - - FFTComplexSetup setup(fftSizelog2); - - // Special cases for short inputs - - if (!sizeOut) return; - - if (sizeOut == 1) - { - ConvolveOp()(rOut[0], iOut[0], rIn1[0], iIn1[0], rIn2[0], iIn2[0], 1.0); - return; - } - - // Assign temporary memory - - TempSpectra spectrum1(fftSize); - TempSpectra spectrum2(fftSize); - - FFT_SPLIT_COMPLEX_D output; - output.realp = rOut; - output.imagp = iOut; - - // Copy to the inputs - - copyZero(spectrum1.mSpectra.realp, rIn1, sizeR1, fftSize); - copyZero(spectrum1.mSpectra.imagp, iIn1, sizeI1, fftSize); - copyZero(spectrum2.mSpectra.realp, rIn2, sizeR2, fftSize); - copyZero(spectrum2.mSpectra.imagp, iIn2, sizeI2, fftSize); - - // Take the Forward FFTs - - transformForward(setup, spectrum1.mSpectra, fftSizelog2); - transformForward(setup, spectrum2.mSpectra, fftSizelog2); - - // Operate - - double scale = 1.0 / (double) fftSize; - binaryOp(spectrum1.mSpectra, spectrum2.mSpectra, fftSize, scale, Op()); - - // Inverse iFFT - - transformInverse(setup, spectrum1.mSpectra, fftSizelog2); - arrangeOutput(output, spectrum1.mSpectra, std::min(size1, size2), sizeOut, - linearSize, fftSize, mode, op); -} - -template -void binarySpectralOperationReal(double* output, const double* in1, - size_t size1, const double* in2, size_t size2, - EdgeMode mode, Op op) -{ - size_t linearSize = calcLinearSize(size1, size2); - size_t sizeOut = calcSize(size1, size2, mode); - size_t fftSizelog2 = ilog2(linearSize); - size_t fftSize = 1 << fftSizelog2; - - FFTRealSetup setup(fftSizelog2); - - // Special cases for short inputs - - if (!sizeOut) return; - - if (sizeOut == 1) - { - output[0] = in1[0] * in2[0]; - return; - } - - // Assign temporary memory - - TempSpectra spectrum1(fftSize >> 1); - TempSpectra spectrum2(fftSize >> 1); - - // Take the Forward Real FFTs - - transformForwardReal(setup, spectrum1.mSpectra, in1, size1, fftSizelog2); - transformForwardReal(setup, spectrum2.mSpectra, in2, size2, fftSizelog2); - - // Operate - - double scale = 0.25 / (double) fftSize; - binaryOpReal(spectrum1.mSpectra, spectrum2.mSpectra, fftSize >> 1, scale, - Op()); - - // Inverse iFFT - - transformInverseReal(setup, spectrum1.mSpectra, fftSizelog2); - arrangeOutput(output, spectrum1.mSpectra, std::min(size1, size2), sizeOut, - linearSize, fftSize, mode, op); -} -} // namespace impl - -// Interface - -// Convolution (Real) - -void convolveReal(double* output, const double* in1, size_t size1, - const double* in2, size_t size2, EdgeMode mode = kEdgeWrap) -{ - impl::binarySpectralOperationReal(output, in1, size1, in2, size2, mode, - impl::ConvolveOp()); -} - -// Correlation (Real) - -void correlateReal(double* output, const double* in1, size_t size1, - const double* in2, size_t size2, EdgeMode mode = kEdgeWrap) -{ - impl::binarySpectralOperationReal(output, in1, size1, in2, size2, mode, - impl::CorrelateOp()); -} - -// Autocorrelation (Real) (inefficient for now) - -void autocorrelateReal(double* output, const double* in, size_t size, - EdgeMode mode = kEdgeWrap) -{ - correlateReal(output, in, size, in, size, mode); -} - -// Convolution (Complex) - -void convolve(double* rOut, double* iOut, const double* rIn1, size_t sizeR1, - const double* iIn1, size_t sizeI1, const double* rIn2, - size_t sizeR2, const double* iIn2, size_t sizeI2, - EdgeMode mode = kEdgeWrap) -{ - impl::binarySpectralOperation(rOut, iOut, rIn1, sizeR1, iIn1, sizeI1, rIn2, - sizeR2, iIn2, sizeI2, mode, impl::ConvolveOp()); -} - -// Correlation (Complex) - -void correlate(double* rOut, double* iOut, const double* rIn1, size_t sizeR1, - const double* iIn1, size_t sizeI1, const double* rIn2, - size_t sizeR2, const double* iIn2, size_t sizeI2, - EdgeMode mode = kEdgeWrap) -{ - impl::binarySpectralOperation(rOut, iOut, rIn1, sizeR1, iIn1, sizeI1, rIn2, - sizeR2, iIn2, sizeI2, mode, - impl::CorrelateOp()); -} - -// Autocorrelation (Complex) (inefficient for now) - -void autocorrelate(double* rOut, double* iOut, const double* rIn, size_t sizeR, - const double* iIn, size_t sizeI, EdgeMode mode = kEdgeWrap) -{ - correlate(rOut, iOut, rIn, sizeR, iIn, sizeI, rIn, sizeR, iIn, sizeI, mode); -} - -} // namespace algorithm -} // namespace fluid diff --git a/include/algorithms/util/FFT.hpp b/include/algorithms/util/FFT.hpp index 3995d8a47..68dabef12 100644 --- a/include/algorithms/util/FFT.hpp +++ b/include/algorithms/util/FFT.hpp @@ -11,59 +11,77 @@ under the European Union’s Horizon 2020 research and innovation programme #pragma once #include "../../data/FluidIndex.hpp" +#include "../../data/FluidMemory.hpp" #include #include namespace fluid { namespace algorithm { -class FFT +namespace impl { +class FFTSetup { - public: - using ArrayXcd = Eigen::ArrayXcd; - using ArrayXcdRef = Eigen::Ref; - using ArrayXd = Eigen::ArrayXd; - using ArrayXdRef = Eigen::Ref; - - FFT() = delete; - - FFT(index size) - : mMaxSize(size), mSize(size), mFrameSize(size / 2 + 1), - mLog2Size(static_cast(std::log2(size))), - mOutputBuffer(mFrameSize), mRealBuffer(mFrameSize), - mImagBuffer(mFrameSize) + FFTSetup(index maxSize) : mMaxSize{maxSize} { - hisstools_create_setup(&mSetup, asUnsigned(mLog2Size)); - mSplit.realp = mRealBuffer.data(); - mSplit.imagp = mImagBuffer.data(); + assert(maxSize > 0 && "FFT Max Size must be > 0!"); + hisstools_create_setup(&mSetup, + asUnsigned(static_cast(std::log2(maxSize)))); } - ~FFT() { hisstools_destroy_setup(mSetup); } - - FFT(const FFT& other) = delete; - - FFT(FFT&& other) { *this = std::move(other); } + ~FFTSetup() + { + hisstools_destroy_setup(mSetup); + mSetup = 0; + } - FFT& operator=(const FFT&) = delete; + FFTSetup(FFTSetup const&) = delete; + FFTSetup& operator=(FFTSetup const&) = delete; - FFT& operator=(FFT&& other) + FFTSetup(FFTSetup&& other) { *this = std::move(other); }; + FFTSetup& operator=(FFTSetup&& other) { using std::swap; - mMaxSize = other.mMaxSize; - mSize = other.mSize; - mFrameSize = other.mFrameSize; - mLog2Size = other.mLog2Size; - swap(mOutputBuffer, other.mOutputBuffer); - swap(mRealBuffer, other.mRealBuffer); - swap(mImagBuffer, other.mImagBuffer); - swap(mSplit, other.mSplit); + swap(mMaxSize, other.mMaxSize); swap(mSetup, other.mSetup); - other.mSetup = nullptr; return *this; } - void resize(index newSize) + FFT_SETUP_D operator()() const noexcept { return mSetup; } + index maxSize() const noexcept { return mMaxSize; } + +private: + FFT_SETUP_D mSetup{nullptr}; + index mMaxSize; +}; +} // namespace impl + +class FFT +{ + +public: + using MapXcd = Eigen::Map; + + static void setup() { getFFTSetup(); } + + FFT() = delete; + + FFT(index size, Allocator& alloc = FluidDefaultAllocator()) + noexcept + : mMaxSize(size), mSize(size), mFrameSize(size / 2 + 1), + mLog2Size(static_cast(std::log2(size))), mSetup(getFFTSetup()), + mRealBuffer(asUnsigned(mFrameSize), alloc), + mImagBuffer(asUnsigned(mFrameSize), alloc), + mOutputBuffer(asUnsigned(mFrameSize), alloc) + {} + + FFT(const FFT& other) = delete; + FFT(FFT&& other) noexcept = default; + + FFT& operator=(const FFT&) = delete; + FFT& operator=(FFT&& other) noexcept = default; + + void resize(index newSize) noexcept { assert(newSize <= mMaxSize); mFrameSize = newSize / 2 + 1; @@ -71,22 +89,31 @@ class FFT mSize = newSize; } - Eigen::Ref process(const ArrayXdRef& input) + MapXcd process(const Eigen::Ref& input) { - hisstools_rfft(mSetup, input.data(), &mSplit, asUnsigned(input.size()), - asUnsigned(mLog2Size)); + + mSplit.realp = mRealBuffer.data(); + mSplit.imagp = mImagBuffer.data(); + hisstools_rfft(mSetup, input.derived().data(), &mSplit, + asUnsigned(input.size()), asUnsigned(mLog2Size)); mSplit.realp[mFrameSize - 1] = mSplit.imagp[0]; mSplit.imagp[mFrameSize - 1] = 0; mSplit.imagp[0] = 0; for (index i = 0; i < mFrameSize; i++) { - mOutputBuffer(i) = + mOutputBuffer[asUnsigned(i)] = 0.5 * std::complex(mSplit.realp[i], mSplit.imagp[i]); } - return mOutputBuffer.segment(0, mFrameSize); + return {mOutputBuffer.data(), mFrameSize}; } protected: + static FFT_SETUP_D getFFTSetup() + { + static const impl::FFTSetup static_setup(65536); + return static_setup(); + } + index mMaxSize{16384}; index mSize{1024}; index mFrameSize{513}; @@ -94,37 +121,42 @@ class FFT FFT_SETUP_D mSetup; FFT_SPLIT_COMPLEX_D mSplit; + rt::vector mRealBuffer; + rt::vector mImagBuffer; private: - ArrayXcd mOutputBuffer; - ArrayXd mRealBuffer; - ArrayXd mImagBuffer; + rt::vector> mOutputBuffer; }; class IFFT : public FFT { public: - IFFT(index size) : FFT(size), mOutputBuffer(size) {} + IFFT(index size, Allocator& alloc = FluidDefaultAllocator()) + : FFT(size, alloc), mOutputBuffer(asUnsigned(size), alloc) + {} - using ArrayXcdRef = Eigen::Ref; - using ArrayXdRef = Eigen::Ref; + using MapXd = Eigen::Map; - Eigen::Ref process(const Eigen::Ref& input) + MapXd process(const Eigen::Ref& input) { + assert(input.size() == mFrameSize); + + mSplit.realp = mRealBuffer.data(); + mSplit.imagp = mImagBuffer.data(); for (index i = 0; i < input.size(); i++) { - mSplit.realp[i] = input[i].real(); - mSplit.imagp[i] = input[i].imag(); + mSplit.realp[i] = input(i).real(); + mSplit.imagp[i] = input(i).imag(); } mSplit.imagp[0] = mSplit.realp[mFrameSize - 1]; hisstools_rifft(mSetup, &mSplit, mOutputBuffer.data(), asUnsigned(mLog2Size)); - return mOutputBuffer.segment(0, mSize); + return {mOutputBuffer.data(), mSize}; } private: - ArrayXd mOutputBuffer; + rt::vector mOutputBuffer; }; } // namespace algorithm } // namespace fluid diff --git a/include/algorithms/util/FluidEigenMappings.hpp b/include/algorithms/util/FluidEigenMappings.hpp index 62a67c70f..af65e0ca5 100644 --- a/include/algorithms/util/FluidEigenMappings.hpp +++ b/include/algorithms/util/FluidEigenMappings.hpp @@ -10,6 +10,7 @@ under the European Union’s Horizon 2020 research and innovation programme #pragma once +#include "../../data/FluidMemory.hpp" #include "../../data/FluidTensor.hpp" #include #include @@ -18,7 +19,6 @@ under the European Union’s Horizon 2020 research and innovation programme namespace fluid { namespace algorithm { -namespace _impl { using Eigen::Map; using Eigen::PlainObjectBase; @@ -29,6 +29,29 @@ using Eigen::ColMajor; using Eigen::Dynamic; using Eigen::RowMajor; +namespace _impl { + + +template +auto asFluid(ScopedEigenMap& a) + -> FluidTensorView +{ + constexpr size_t N = EigenType::IsVectorAtCompileTime ? 1 : 2; + + if constexpr (N == 2) + { + if (static_cast(a.Options) == static_cast(ColMajor)) + { + // Respect the colmajorness of an eigen type + auto slice = FluidTensorSlice(0, {a.cols(), a.rows()}); + return {slice.transpose(), a.data()}; + } + return {a.data(), 0, a.rows(), a.cols()}; + } + else + return {a.data(), 0, a.rows()}; +} /// Eigen Matrix/Array -> FluidTensorView template @@ -39,7 +62,7 @@ auto asFluid(PlainObjectBase& a) { constexpr size_t N = PlainObjectBase::IsVectorAtCompileTime ? 1 : 2; - if (N == 2) + if constexpr (N == 2) { if (static_cast(a.Options) == static_cast(ColMajor)) { @@ -62,7 +85,7 @@ auto asFluid(const PlainObjectBase& a) { constexpr size_t N = PlainObjectBase::IsVectorAtCompileTime ? 1 : 2; - if (N == 2) + if constexpr (N == 2) { if (static_cast(a.Options) == static_cast(ColMajor)) { @@ -201,7 +224,12 @@ auto asEigen(const FluidTensorView&& return asEigen(a); } - } // namespace _impl + +template