marvin_Interpolators.h

Go to the documentation of this file.

// ========================================================================================================
//  _______ _______ ______ ___ ___ _______ _______
// |   |   |   _   |   __ \   |   |_     _|    |  |
// |       |       |      <   |   |_|   |_|       |
// |__|_|__|___|___|___|__|\_____/|_______|__|____|
//
// This file is part of the Marvin open source library and is licensed under the terms of the MIT License.
//
// ========================================================================================================
 
#ifndef MARVIN_INTERPOLATORS_H
#define MARVIN_INTERPOLATORS_H
#include <marvin/library/marvin_Concepts.h>
#include <marvin/library/marvin_Literals.h>
#include <marvin/math/marvin_Math.h>
#include <marvin/math/marvin_Windows.h>
#include <span>
#include <cassert>
#include <vector>
#include <functional>
namespace marvin::math::interpolators {
    template <FloatType SampleType, size_t N, windows::WindowType WindowType>
    requires(N % 2 == 0)
    class WindowedSincInterpolator {
    public:
        WindowedSincInterpolator()
        requires(WindowType == windows::WindowType::Sine)
        {
            auto windowFunc = [](SampleType i) -> SampleType {
                const auto res = windows::sine(i, static_cast<SampleType>(N));
                return res;
            };
            fillLookupTable(std::move(windowFunc));
        }

        explicit WindowedSincInterpolator(SampleType alpha)
        requires(WindowType == windows::WindowType::Tukey)
        {
            auto windowFunc = [alpha](SampleType i) -> SampleType {
                const auto res = windows::tukey(i, static_cast<SampleType>(N), alpha);
                return res;
            };
            fillLookupTable(std::move(windowFunc));
        }

        [[nodiscard]] SampleType interpolate(std::span<SampleType> sampleContext, SampleType ratio) {
            assert(sampleContext.size() == N);
            auto sum = static_cast<SampleType>(0.0);
            const auto inverseRatio = static_cast<SampleType>(1.0) - ratio;
            constexpr static auto halfN = static_cast<int>(N) / 2;
            for (auto i = 0_sz; i < N; ++i) {
                auto n = static_cast<int>(i) - static_cast<int>(halfN);
                const auto windowedSinc = lookupSinc(static_cast<SampleType>(n), inverseRatio);
                const auto res = sampleContext[i] * windowedSinc;
                sum += res;
            }
            return sum;
        }
 
    private:
        void fillLookupTable(std::function<SampleType(SampleType)>&& windowFunction) {
            m_lut.reserve(N + 2);
            m_lut.resize(N);
            constexpr static auto halfWidth{ N / 2 };
            for (auto i = 0_sz; i < N; ++i) {
                const auto n{ static_cast<int>(i) - static_cast<int>(halfWidth) };
                const auto window = windowFunction(static_cast<SampleType>(i));
                const auto sincRes = sinc(static_cast<SampleType>(n));
                const auto windowedSinc = sincRes * window;
                m_lut[i] = windowedSinc;
            }
            const auto start = m_lut.front();
            const auto end = m_lut[m_lut.size() - 1];
            m_lut.emplace_back(end);
            m_lut.emplace_back(start);
            std::rotate(m_lut.begin(), m_lut.end() - 1, m_lut.end()); // [ x x x x e s ] -> [s x x x x e]
        }
 
 
        [[nodiscard]] SampleType lookupSinc(SampleType x, SampleType frac) {
            constexpr static auto halfN{ static_cast<int>(N / 2) };
            const auto nBelow = static_cast<int>(std::floor(x));
            const auto nAbove = static_cast<int>(nBelow) + 1;
            const auto ratio = frac;
            // const auto ratio = x - static_cast<SampleType>(nBelow);
            // So this is now in range -N to N
            // translating this to an actual index into the LUT means just adding (N / 2)..
            // LUT has an extra point at the start compared to what we actually expect....
            const auto idxBelow = static_cast<size_t>(nBelow + halfN) + 1;
            const auto idxAbove = static_cast<size_t>(nAbove + halfN) + 1;
            const auto start = m_lut[idxBelow];
            const auto end = m_lut[idxAbove];
            const auto interpolated = math::lerp<SampleType>(start, end, ratio);
            return interpolated;
        }
 
        std::vector<SampleType> m_lut;
    };
} // namespace marvin::math::interpolators
 
#endif