SWIPE_PitchEstimation.cpp

Go to the documentation of this file.

//=====================================================================================
//=====================================================================================
namespace krotos
{
    SWIPE_PitchEstimation::SWIPE_PitchEstimation(const float sampleRate, const float timeStep,
                                                 const float log2PitchStep, const float ERBStep,
                                                 const float normOverlap, const float strenThresh)
    {
        // TODO: If they are not passed as arguments then resort to these values
        m_sampleRate = sampleRate;
        m_timeStep = timeStep;
        m_log2PitchStep = log2PitchStep;
        m_ERBStep = ERBStep;
        m_normOverlap = normOverlap;
        m_strenThresh = strenThresh;
 
        // initialize vectors with pitch candidates, winSizes and ERBs for pitch estimation calc
        initVectors();
    }
 
    void SWIPE_PitchEstimation::initVectors()
    {
        m_pitchLims << 50.0f, 500.0f;
 
        // Define pitch candidates
        Eigen::VectorXf log2PitchCands = calculateLog2PitchCands(m_log2PitchStep);
        m_pitchCands = calculatePitchCands(log2PitchCands);
 
        // Determine P2-Winsizes
        Eigen::Vector2f logWinSizeLimsEigen;
        logWinSizeLimsEigen[0] = std::round(std::log2(8.0f * m_sampleRate / m_pitchLims[0]));
        logWinSizeLimsEigen[1] = std::round(std::log2(8.0f * m_sampleRate / m_pitchLims[1]));
 
        size_t vectorSize = static_cast<size_t>(logWinSizeLimsEigen[0] - logWinSizeLimsEigen[1]) + 1;
        m_winSizes.resize(vectorSize);
        for (size_t i = 0; i < vectorSize; i++)
        {
            m_winSizes[i] = std::pow(2.0f, logWinSizeLimsEigen[0] - i);
        }
 
        // Calculate winSizeOptPitches
        Eigen::VectorXf winSizeOptPitchesEigen(m_winSizes.size());
        for (size_t i = 0; i < m_winSizes.size(); ++i)
        {
            winSizeOptPitchesEigen[i] = 8.0f * m_sampleRate / m_winSizes[i];
        }
 
        // Determine window sizes used by each pitch candidate
        m_log2DistWin1OptPitchAndPitchCands.resize(log2PitchCands.size());
        for (size_t i = 0; i < log2PitchCands.size(); ++i)
        {
            m_log2DistWin1OptPitchAndPitchCands[i] = log2PitchCands[i] - std::log2(winSizeOptPitchesEigen[0]) + 1.0f;
        }
 
        // Create ERB - scale uniformly - spaced frequencies(in Hertz)
        float start = krotos::ERB_FFTSpectrogram::hz2Erbs(m_pitchCands.minCoeff() / 4.0f);
        float end = krotos::ERB_FFTSpectrogram::hz2Erbs(m_sampleRate / 2.0f);
        size_t size = static_cast<size_t>((end - start) / m_ERBStep) + 1;
        Eigen::VectorXf erbsRangeEigen = Eigen::VectorXf::LinSpaced(size, start, end);
 
        // Calculate corresponding Hz values using erbs2Hz function
        m_ERBs = krotos::ERB_FFTSpectrogram::erbs2Hz(erbsRangeEigen);
    }
 
    std::pair<Eigen::VectorXf, Eigen::VectorXf> SWIPE_PitchEstimation::calculatePitchStrengthPairs(
            std::vector<float>& signalSTL)
    {
        // At this point everything needs to be intialized
        jassert(m_sampleRate != -1.0f && m_timeStep != -1.0f && m_log2PitchStep != -1.0f && m_ERBStep != -1.0f &&
                m_normOverlap != -1.0f && m_strenThresh != -1.0f);
 
        // TODO:: There is a copy here, need to go one level above and pass VectorXf when called!!!!
        Eigen::VectorXf signal = Eigen::VectorXf::Map(&signalSTL[0], signalSTL.size());
 
        // Calculate time vector
        m_times = Eigen::VectorXf::LinSpaced(signal.size() / m_sampleRate / m_timeStep + 1, 0.0f,
                                             signal.size() / m_sampleRate);
 
        // Initialize pitch strength matrix
        Eigen::MatrixXf globStrenMtrx(m_pitchCands.size(), m_times.size());
        globStrenMtrx.setZero();
 
        // Vectors holding the STFT results
        Eigen::MatrixXf distr;
        // Frequency bins axis and time frame axis, will be used for interpolation later
        Eigen::VectorXf fSupport, tSupport;
 
        // optPitchCandsIdcs and imperfectFitIdcs vectors
        std::pair<Eigen::VectorXf, Eigen::VectorXf> candidates;
 
        // Main loop
        for (size_t winSizeIdx = 0; winSizeIdx < m_winSizes.size(); winSizeIdx++)
        {
            // Hop size
            float hopSize = std::max(1.0f, (std::round((1.0f - m_normOverlap) * m_winSizes[winSizeIdx])));
 
            // Zero Pad signal
            // Calculate the zero-padding size
            size_t padSize = static_cast<size_t>(m_winSizes[winSizeIdx] / 2);
            // Append zeros to the beginning using Eigen - the end zeros are handled from the STFT algorithm
            Eigen::VectorXf paddedSignal(signal.size() + padSize);
            paddedSignal << Eigen::VectorXf::Zero(padSize), signal;
 
            // -- Compute spectrum --
            // Init stft again
            std::unique_ptr<ShortTimeFourierTransform> stftObj = std::make_unique<ShortTimeFourierTransform>(
                    static_cast<int>(m_winSizes[winSizeIdx]), static_cast<int>(m_winSizes[winSizeIdx]),
                    static_cast<int>(hopSize), WindowType::Hann, ShortTimeFourierTransform::fftMode::freqMagOnly,
                    static_cast<int>(m_sampleRate));
 
            // Perform stft
            distr = stftObj->processSignal(paddedSignal);
 
            // Get the frequency and time vectors
            fSupport = stftObj->getBins(ShortTimeFourierTransform::FrequencyRange::Half);
            // TODO: Check why diff from MATLAB -> Problem here!!
            tSupport = stftObj->getColumns();
 
            // Select candidates that use this window size
            candidates = selectCandidatesForWindowSize(winSizeIdx + 1);
 
            // Compute loudness at ERBs uniformly-spaced frequencies
            auto ERBsSubset = getERBsSubset(candidates.first);
 
            // Interpolate abs stft matrix at ERB frequency axis and square
            auto ERBsDistrValues = ShapePreservingCubicInterpolator::interpolate2DEigen(fSupport, distr, ERBsSubset);
            ERBsDistrValues = ERBsDistrValues.unaryExpr([](float val) { return std::sqrt(std::max(0.0f, val)); });
 
            // Compute pitch strength
            // Collect pitchCands at specified indices
            Eigen::VectorXf selectedPitchCands;
            for (int i = 0; i < candidates.first.size(); ++i)
            {
                int index = static_cast<int>(candidates.first[i]);
 
                // to get the right output shapes for all window sizes we need <=
                if (index > 0 && index <= m_pitchCands.size())
                {
                    selectedPitchCands.conservativeResize(selectedPitchCands.size() + 1);
                    selectedPitchCands(selectedPitchCands.size() - 1) = m_pitchCands[index - 1];
                }
            }
 
            Eigen::MatrixXf locStrenMtrxEigen =
                    pitchStrengthAllCandidates(ERBsSubset, ERBsDistrValues, selectedPitchCands);
 
            // Interpolate pitch strength at desired times
            Eigen::MatrixXf locStrenMtrxInterpEigen;
            if (locStrenMtrxEigen.size() > 1)
            {
                locStrenMtrxEigen.transposeInPlace();
                locStrenMtrxInterpEigen = LinearInterpolator::interpolate2DEigen(tSupport, locStrenMtrxEigen, m_times);
                locStrenMtrxInterpEigen.transposeInPlace();
            }
            else
            {
                locStrenMtrxInterpEigen.resize(candidates.first.size(), m_times.size());
                locStrenMtrxInterpEigen.setConstant(std::numeric_limits<float>::quiet_NaN());
            }
 
            // Add pitch strength to combination
            Eigen::VectorXf currWinRelStrenEigen(candidates.first.size());
            currWinRelStrenEigen.setOnes();
            Eigen::VectorXf log2DistCurrWinOptPitchAndOptPitchCandsEigen;
 
            auto optPitchCandsIdcs = candidates.first.cast<int>().array() - 1;
            auto imperfectFitIdcs = candidates.second.cast<int>().array() - 1;
            auto log2DistCurrWinOptPitchAndOptPitchCands =
                    m_log2DistWin1OptPitchAndPitchCands(optPitchCandsIdcs(imperfectFitIdcs)).array() - (winSizeIdx + 1);
 
            currWinRelStrenEigen(imperfectFitIdcs) = (1 - log2DistCurrWinOptPitchAndOptPitchCands.abs().array());
 
            // vectorise accumulating into globStrenMtrxEigen
            globStrenMtrx(optPitchCandsIdcs, Eigen::all) =
                    globStrenMtrx(optPitchCandsIdcs, Eigen::all).array() +
                    (currWinRelStrenEigen.replicate(1, locStrenMtrxInterpEigen.cols()).array() *
                     locStrenMtrxInterpEigen.array());
        }
 
        return fineTunePitch(globStrenMtrx);
    }
 
    Eigen::VectorXf SWIPE_PitchEstimation::calculateLog2PitchCands(float log2PitchStep)
    {
        // Calculate the size of the resulting vector
        int vectorSize = static_cast<int>((std::log2(m_pitchLims[1]) - std::log2(m_pitchLims[0])) / log2PitchStep) + 1;
 
        // Initialize Eigen vector with the calculated size
        Eigen::VectorXf log2PitchCands(vectorSize);
 
        // Populate the vector using Eigen's sequence generation
        for (int i = 0; i < vectorSize; ++i)
        {
            log2PitchCands(i) = std::log2(m_pitchLims[0]) + (float)i * log2PitchStep;
        }
 
        return log2PitchCands;
    }
 
    Eigen::VectorXf SWIPE_PitchEstimation::calculatePitchCands(const Eigen::VectorXf& log2PitchCands)
    {
        Eigen::VectorXf pitchCands(log2PitchCands.size());
 
        // Perform an element-wise power of 2
        pitchCands = Eigen::pow(2.0f, log2PitchCands.array());
 
        return pitchCands;
    }
 
    std::pair<Eigen::VectorXf, Eigen::VectorXf> SWIPE_PitchEstimation::selectCandidatesForWindowSize(size_t winSizeIDx)
    {
        Eigen::VectorXf optPitchCandsIdcs;
        Eigen::VectorXf imperfectFitIdcs;
 
        if (m_winSizes.size() == 1)
        {
            optPitchCandsIdcs =
                    Eigen::VectorXf::LinSpaced(m_pitchCands.size(), 1.0f, static_cast<float>(m_pitchCands.size()));
        }
        else if (winSizeIDx == m_winSizes.size())
        {
            for (size_t j = 0; j < m_log2DistWin1OptPitchAndPitchCands.size(); ++j)
            {
                if (m_log2DistWin1OptPitchAndPitchCands[j] - static_cast<float>(winSizeIDx) > -1.0f)
                {
                    optPitchCandsIdcs.conservativeResize(optPitchCandsIdcs.size() + 1);
                    optPitchCandsIdcs(optPitchCandsIdcs.size() - 1) = j + 1.0f;
                }
            }
 
            for (size_t j = 0; j < optPitchCandsIdcs.size(); ++j)
            {
                if (m_log2DistWin1OptPitchAndPitchCands[static_cast<size_t>(optPitchCandsIdcs(j)) - 1.0f] -
                            static_cast<float>(winSizeIDx) <
                    0.0f)
                {
                    imperfectFitIdcs.conservativeResize(imperfectFitIdcs.size() + 1);
                    imperfectFitIdcs(imperfectFitIdcs.size() - 1) = j + 1.0f;
                }
            }
        }
        else if (winSizeIDx == 1)
        {
            for (size_t j = 0; j < m_log2DistWin1OptPitchAndPitchCands.size(); ++j)
            {
                if (m_log2DistWin1OptPitchAndPitchCands[j] - static_cast<float>(winSizeIDx) < 1.0f)
                {
                    optPitchCandsIdcs.conservativeResize(optPitchCandsIdcs.size() + 1);
                    optPitchCandsIdcs(optPitchCandsIdcs.size() - 1) = j + 1.0f;
                }
            }
 
            for (size_t j = 0; j < optPitchCandsIdcs.size(); ++j)
            {
                if (m_log2DistWin1OptPitchAndPitchCands[static_cast<size_t>(optPitchCandsIdcs(j)) - 1.0f] -
                            static_cast<float>(winSizeIDx) >
                    0.0f)
                {
                    imperfectFitIdcs.conservativeResize(imperfectFitIdcs.size() + 1);
                    imperfectFitIdcs(imperfectFitIdcs.size() - 1) = j + 1.0f;
                }
            }
        }
        else
        {
            for (Eigen::Index j = 0; j < m_log2DistWin1OptPitchAndPitchCands.size(); ++j)
            {
                if (std::abs(m_log2DistWin1OptPitchAndPitchCands[j] - static_cast<float>(winSizeIDx)) < 1.0f)
                {
                    optPitchCandsIdcs.conservativeResize(optPitchCandsIdcs.size() + 1);
                    optPitchCandsIdcs(optPitchCandsIdcs.size() - 1) = j + 1.0f;
                }
            }
            imperfectFitIdcs = Eigen::VectorXf::LinSpaced(optPitchCandsIdcs.size(), 1.0f,
                                                          static_cast<float>(optPitchCandsIdcs.size()));
        }
 
        return {optPitchCandsIdcs, imperfectFitIdcs};
    }
 
    Eigen::VectorXf SWIPE_PitchEstimation::getERBsSubset(const Eigen::VectorXf& pitchCandsSubset)
    {
        float threshold = m_pitchCands[static_cast<int>(pitchCandsSubset[0])] / 4.0f;
 
        // Find the index of the first element greater than the threshold
        Eigen::VectorXf::Index startIndex = -1;
        for (Eigen::VectorXf::Index i = 0; i < m_ERBs.size(); ++i)
        {
            if (m_ERBs(i) > threshold)
            {
                startIndex = i;
                break;
            }
        }
 
        // Extract the subset from the found position to the end
        Eigen::VectorXf ERBsSubset;
        if (startIndex != -1)
        {
            ERBsSubset = m_ERBs.segment(startIndex, m_ERBs.size() - startIndex);
        }
 
        return ERBsSubset;
    }
 
    Eigen::MatrixXf SWIPE_PitchEstimation::pitchStrengthAllCandidates(const Eigen::VectorXf& _ERBs,
                                                                      Eigen::MatrixXf& _ERBsDistrValues,
                                                                      const Eigen::VectorXf& _pitchCands)
    {
        // Create pitch strength matrix
        Eigen::MatrixXf locStrenMtrx(_pitchCands.size(), _ERBsDistrValues.cols());
        locStrenMtrx.setZero();
 
        // Define integration regions
        Eigen::VectorXf k(_pitchCands.size() + 1);
        k.setZero();
 
        for (int j = 0; j < k.size() - 1; ++j)
        {
            // Find index
            auto startIdx = static_cast<size_t>(k[j]);
            auto it = std::find_if(_ERBs.data() + startIdx, _ERBs.data() + _ERBs.size(),
                                   [&](float erb) { return erb > _pitchCands[j] / 4.0f; });
 
            // Update k based on the found index
            if (it != _ERBs.data() + _ERBs.size())
            {
                float dist = static_cast<float>(std::distance(_ERBs.data(), it));
                k(j + 1) = dist;
            }
            else
            {
                k(j + 1) = k(j);
            }
        }
 
        // Remove the first element of k
        Eigen::VectorXf kKeep = k.segment(1, k.size() - 1);
 
        // Create loudness normalization matrix
        auto N = computeN(_ERBsDistrValues);
 
        for (size_t j = 0; j < _pitchCands.size(); j++)
        {
            // Normalize loudness
            auto n = N.row(kKeep[j]).transpose(); // Extract row k(j). Do I transpose ??
            // Make zero-loudness equal zero after normalization
            std::replace_if(
                    n.begin(), n.end(), [](float val) { return val == 0; }, INFINITY);
            // NL = ERBsDistrValues(k(j):end, :) ./ repmat(n, size(ERBsDistrValues, 1)-k(j)+1, 1)
            Eigen::MatrixXf NL = _ERBsDistrValues.bottomRows(_ERBsDistrValues.rows() - kKeep[j]);
            NL.array().rowwise() /= n.transpose().array();
 
            // Compute pitch strength
            locStrenMtrx.row(j) = pitchStrengthOneCandidate(_ERBs.segment(kKeep[j], _ERBs.size() - kKeep[j]), NL,
                                                            _pitchCands[j]); // sliiiightly different;
        }
 
        return locStrenMtrx;
    }
 
    Eigen::VectorXf SWIPE_PitchEstimation::pitchStrengthOneCandidate(const Eigen::VectorXf& ERBS,
                                                                     const Eigen::MatrixXf& normERBsDistrValues,
                                                                     float pitchCand)
    {
        int n = static_cast<int>(std::floor(ERBS[ERBS.size() - 1] / pitchCand - 0.75)); // Number of harmonics
 
        if (n == 0)
        {
            return Eigen::VectorXf::Constant(normERBsDistrValues.cols(), std::numeric_limits<float>::quiet_NaN());
        }
 
        Eigen::VectorXf k = Eigen::VectorXf::Zero(ERBS.size()); // Kernel
 
        // Normalize frequency w.r.t. candidate
        Eigen::VectorXf q = ERBS.array() / pitchCand;
 
        auto primesVec = primes(n);
        // Create Kernel
        for (int i : primesVec) // primes up to 17
        {
            Eigen::VectorXf a = (q.array() - static_cast<float>(i)).cwiseAbs();
 
            // Peak's weight
            for (size_t p = 0; p < a.size(); ++p)
            {
                if (a[p] < 0.25f)
                {
                    k[p] = std::cos(2.0f * M_PI * q[p]);
                }
            }
 
            // Valleys' weights
            for (size_t v = 0; v < a.size(); ++v)
            {
                if (a[v] > 0.25f && a[v] < 0.75f)
                {
                    k[v] += std::cos(2.0f * M_PI * q[v]) / 2.0f;
                }
            }
        }
 
        // Apply envelope
        k.array() *= (1.0f / ERBS.array()).array().sqrt();
 
        // K+-normalize kernel
        Eigen::VectorXf nonZeroElements;
        for (size_t i = 0; i < k.size(); i++)
        {
            if (k[i] > 0.0f)
            {
                nonZeroElements.conservativeResize(nonZeroElements.size() + 1);
                nonZeroElements(nonZeroElements.size() - 1) = k[i];
            }
        }
        float normNonZero = nonZeroElements.norm();
        k.array() /= normNonZero;
 
        // Compute pitch strength
        Eigen::VectorXf locStrenLine = k.transpose() * normERBsDistrValues;
        return locStrenLine;
    }
 
    std::pair<Eigen::VectorXf, Eigen::VectorXf> SWIPE_PitchEstimation::fineTunePitch(Eigen::MatrixXf globStrMatrix)
    {
        size_t numTimes = globStrMatrix.cols();
        size_t numPitches = globStrMatrix.rows();
 
        Eigen::VectorXf pitches(numTimes);
        pitches.setConstant(std::numeric_limits<float>::quiet_NaN());
 
        Eigen::VectorXf strengths(numTimes);
        strengths.setConstant(std::numeric_limits<float>::quiet_NaN());
 
        for (size_t timeIdx = 0; timeIdx < numTimes; ++timeIdx)
        {
            // Find the pitch with maximum strength
            // Extract the column
            Eigen::VectorXf column = globStrMatrix.col(timeIdx);
            // Find the maximum and it's index
            // Eigen::Index maxStrenPitchIdx;
            // auto maxStrength = column.maxCoeff(&maxStrenPitchIdx);
 
            float maxStrength = std::numeric_limits<float>::quiet_NaN();
            Eigen::Index maxStrenPitchIdx = 0;
            for (Eigen::Index i = 0; i < column.size(); ++i)
            {
                if (!std::isnan(column[i]))
                {
                    if (std::isnan(maxStrength) || column[i] > maxStrength)
                    {
                        maxStrength = column[i];
                        maxStrenPitchIdx = i;
                    }
                }
            }
 
            strengths[timeIdx] = maxStrength;
 
            // Check strength threshold
            if (strengths[timeIdx] < m_strenThresh)
                continue;
 
            // Check conditions for pitch assignment
            if (maxStrenPitchIdx == 0 || maxStrenPitchIdx == numPitches - 1)
            {
                pitches[timeIdx] = m_pitchCands[maxStrenPitchIdx];
            }
            else
            {
                // Perform polyfit for pitch assignment
                Eigen::VectorXf maxStrenPitchLocIdcs(3);
                maxStrenPitchLocIdcs << maxStrenPitchIdx - 1.0f, maxStrenPitchIdx, maxStrenPitchIdx + 1.0f;
 
                Eigen::VectorXf tc(maxStrenPitchLocIdcs.size());
                Eigen::VectorXf ntc(tc.size());
 
                for (size_t i = 0; i < tc.size(); ++i)
                {
                    tc[i] = 1.0f / m_pitchCands[maxStrenPitchLocIdcs[i]];
                }
 
                for (size_t i = 0; i < tc.size(); ++i)
                {
                    ntc[i] = ((tc[i] / tc[1]) - 1.0f) * (2.0f * static_cast<float>(M_PI));
                }
 
                // Perform polyfit
                Eigen::VectorXf y;
                for (size_t i = 0; i < maxStrenPitchLocIdcs.size(); i++)
                {
                    y.conservativeResize(y.size() + 1);
                    y(y.size() - 1) = globStrMatrix.row(maxStrenPitchLocIdcs[i])[timeIdx];
                }
                Eigen::VectorXf c = polyfit(ntc, y, 2);
 
                // Calculate ftcs
                float startVal = std::log2(m_pitchCands[maxStrenPitchLocIdcs[0]]);
                float endVal = std::log2(m_pitchCands[maxStrenPitchLocIdcs[2]]);
                int numSteps = static_cast<int>((endVal - startVal) / (1.0f / 12.0f / 100.0f) + 1);
                Eigen::VectorXf ftc = Eigen::VectorXf::LinSpaced(numSteps, startVal, endVal);
                ftc = Eigen::pow(2.0f, ftc.array());
                ftc = 1.0f / ftc.array();
 
                // Calculate ftcs
                Eigen::VectorXf nftc(ftc.size());
                nftc = (ftc.array() / tc(1) - 1.0f) * 2.0f * M_PI;
 
                // Evaluate polynomial
                Eigen::VectorXf polyValRes(nftc.size());
                for (int i = 0; i < nftc.size(); i++)
                {
                    polyValRes[i] = Eigen::poly_eval(c, nftc[i]);
                }
 
                // Find the maximum value and its index
                auto maxIt = std::max_element(polyValRes.begin(), polyValRes.end());
                strengths[timeIdx] = *maxIt;
                size_t maxPolyFitIdcs = std::distance(polyValRes.begin(), maxIt);
 
                float pitchValue = m_pitchCands[static_cast<size_t>(maxStrenPitchLocIdcs[0])];
                float maxPolyFitValue = static_cast<float>(maxPolyFitIdcs) / 12.0f / 100.0f;
 
                pitches[timeIdx] = std::pow(2.0f, (std::log2(pitchValue) + maxPolyFitValue));
            }
        }
 
        return {pitches, strengths};
    }
 
    Eigen::MatrixXf SWIPE_PitchEstimation::computeN(const Eigen::MatrixXf& ERBsDistrValues)
    {
        // Step 1: Element-wise square
        Eigen::MatrixXf squaredValues = ERBsDistrValues.array().square();
 
        // Step 2: Flip rows upside down
        Eigen::MatrixXf flipSquaredValues = squaredValues.colwise().reverse().eval();
 
        // Step 3: Cumulative sum across columns
        Eigen::MatrixXf cumSumFlipSquaredValues = cumsum(flipSquaredValues);
 
        // Step 4: Flip again
        Eigen::MatrixXf flipCumSumFlipSquaredValues = cumSumFlipSquaredValues.colwise().reverse().eval();
 
        // Step 5: Calculate square root
        Eigen::MatrixXf rootValues = flipCumSumFlipSquaredValues.array().sqrt();
 
        return rootValues;
    }
 
    Eigen::VectorXi SWIPE_PitchEstimation::primes(int n)
    {
        Eigen::VectorXi primesList(2);
        primesList << 1, 2;
 
        for (int i = 3; i <= n; ++i)
        {
            Eigen::VectorXi range = Eigen::VectorXi::LinSpaced(i - 2, 2, i - 1);
            Eigen::VectorXi remainder = range.unaryExpr([&](int x) { return i % x; });
 
            bool isPrime = remainder.all();
 
            if (isPrime)
            {
                primesList.conservativeResize(primesList.size() + 1);
                primesList(primesList.size() - 1) = i;
            }
        }
 
        return primesList;
    }
 
    Eigen::VectorXf SWIPE_PitchEstimation::polyfit(const Eigen::VectorXf& x, const Eigen::VectorXf& y, int n)
    {
        jassert(x.size() == y.size());
 
        Eigen::MatrixXf xs(x.size(), n + 1);
        xs.col(0).setOnes();
        for (int i = 1; i <= n; ++i)
        {
            xs.col(i).array() = xs.col(i - 1).array() * x.array();
        }
 
        Eigen::VectorXf result(n + 1);
        // Use QR decomposition via Householder reflections. List of other decompositions in:
        // https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html and
        // https://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html
        auto decomposition = xs.householderQr();
        result = decomposition.solve(y);
 
        // Careful: Returns the vector of coef in the opposite order of MATLABs polyfit, but
        // combines correctly with Eigen::poly_eval to produce the same value!!
        return result;
    }
 
    Eigen::MatrixXf SWIPE_PitchEstimation::cumsum(const Eigen::MatrixXf& input)
    {
        Eigen::MatrixXf result(input.rows(), input.cols());
 
        for (int col = 0; col < input.cols(); ++col)
        {
            result(0, col) = input(0, col);
            for (int row = 1; row < input.rows(); ++row)
            {
                result(row, col) = result(row - 1, col) + input(row, col);
            }
        }
 
        return result;
    }
} // namespace krotos