HarmonicRepresentation.cpp

Go to the documentation of this file.

//=====================================================================================
//=====================================================================================
 
#include "HarmonicRepresentation.h"
 
//=====================================================================================
namespace krotos
{
    HarmonicRepresentation::HarmonicRepresentation(float sampleRate) : m_sampleRate(sampleRate)
    {
        float innerExpression = std::round(std::log2(8.0f * sampleRate / 50.0f)) - 1.0f;
        // calculate the maximum of hopSize and innerExpression divided by hopSize
        m_maxSwipepWinSize = config.hopSize * std::ceil(std::pow(2.0f, innerExpression) / config.hopSize);
    }
 
    Eigen::MatrixXf HarmonicRepresentation::processSignal(std::vector<float>& inputSignal)
    {
        // TODO:: There is a copy here, need to go one level above and pass VectorXf when called!!!!
        Eigen::VectorXf signal = Eigen::VectorXf::Map(&inputSignal[0], inputSignal.size());
 
        // -- Compute spectrum -- TODO::Debug and rewrite STFT as part of
        // https://krotosltd.atlassian.net/jira/software/c/projects/RD/boards/43?selectedIssue=RD-271 Init stft
        auto stftObj = std::make_unique<ShortTimeFourierTransform>(
                config.winSize, config.fftSize, config.hopSize, WindowType::BlackmanHarris,
                ShortTimeFourierTransform::fftMode::freqMagOnly, static_cast<int>(m_sampleRate));
 
        // If the window is centred at t, this is the starting index at which to
        // look up the sound.value which you want to multiply by the window.It is a
        // negative number because(almost) half of the window will be before time t
        // and half after.In fact, if the length of the window N is an even number,
        // it is set up so this number equals(N / 2 - 1).If the length of the window
        // is odd, this number equals(N - 1) / 2.
        int winFirstRelIdx = std::floor((config.winSize - 1) / 2.0f);
        // This is the last index at which to look up signal values and is equal to
        // (N - 1) / 2 if the length N of the window is odd and N / 2 if the length of the
        // window is even.This means that in the even case, the window has an
        // unequal number of past and future values, i.e., time t is not the centre
        // of the window, but slightly to the left of the centre of the window
        // (before it).
        int winLastRelIdx = std::ceil((config.winSize - 1) / 2.0f);
 
        // Perform stft
        auto distr = stftObj->processSignal(signal);
        // Only return frequencies below Nyquist rate
        distr = distr.block(0, 0, config.fftSize / 2, distr.cols()).eval();
        // Get the power distribution
        distr = distr.array().square();
        // Remove window energy
        float windowEnergy = stftObj->dftParams.window.sum();
        windowEnergy = windowEnergy * windowEnergy;
        distr = distr / windowEnergy;
 
        // Access stft column vector
        auto tSupport = stftObj->getColumns();
        tSupport = tSupport.array() - (static_cast<float>(winFirstRelIdx) / m_sampleRate);
 
        // Access stft frequency vector
        auto fSupport = stftObj->getBins(ShortTimeFourierTransform::FrequencyRange::Half);
        // Only return frequencies below Nyquist
        fSupport = fSupport.head(config.fftSize / 2).eval();
 
        constexpr int nFreqCorrs = static_cast<int>(10.0f / 0.1f) + 1;
        Eigen::VectorXf freqCorrs = Eigen::VectorXf::LinSpaced(nFreqCorrs, -5.0f, 5.0f);
        int nInharmCoeffs = static_cast<int>((0.001f - 0.0f) / 0.00005f) + 1;
        Eigen::VectorXf inharmCoeffs = Eigen::VectorXf::LinSpaced(nInharmCoeffs, 0.0f, 0.001f);
 
        // Adjusting time support vector
        int startIdx = 0;
        int endIdx = tSupport.size() - 1;
        tSupport.segment(startIdx, endIdx - startIdx + 1).array() -= tSupport(startIdx);
        int tSize = endIdx - startIdx + 1;
        Eigen::MatrixXf distrExtr = distr.middleCols(startIdx, endIdx - startIdx + 1); // Not the same size
 
        // Init Swipe -> TODO:: Make it a member var to avoid reinstantiation of vectors
        auto swipeObj = std::make_unique<SWIPE_PitchEstimation>(m_sampleRate, config.hopSize / m_sampleRate,
                                                                swipe.log2PitchStep, swipe.ERBStep, swipe.normOverlap,
                                                                swipe.strenThresh);
 
        auto pitchStrengthPairs = swipeObj->calculatePitchStrengthPairs(inputSignal); // 1. estPitches, 2. estStrengths
        Eigen::VectorXf estTimes = swipeObj->getTimes();
 
        // replace NaNs in pitches with median
        replaceNanWithMedian(pitchStrengthPairs.first);
        // replace NaN values with -Inf using array operations in pitch strengths
        pitchStrengthPairs.second = pitchStrengthPairs.second.unaryExpr(
                [](float x) { return std::isnan(x) ? -std::numeric_limits<float>::infinity() : x; });
 
        // This matrix is the return value of the function and contains the harmonic information of the signal
        Eigen::MatrixXf value;
        Eigen::VectorXf fundamentalFreqs;
 
        // determine if sample is harmonic or inharmonic with respect to threshold
        if (pitchStrengthPairs.second.maxCoeff() > harmRep.threshold)
        {
            Eigen::MatrixXf estTimePitchPairs(estTimes.size(), 2);
 
            jassert(m_maxSwipepWinSize != 0.0f && m_sampleRate != 0.0f); // Something went wrong
            estTimePitchPairs.col(0) = estTimes.array();
            estTimePitchPairs.col(1) = pitchStrengthPairs.first;
 
            // Estimate f0 at times at which spectrogram was evaluated by interpolating
            // between times when f0 was evaluated
            fundamentalFreqs = fevalbp(estTimePitchPairs, tSupport);
            fundamentalFreqs.transpose();
        }
        else
        {
            // If sound not harmonic, we just fill in with 0s as fundamental estimate,
            // otherwise this will cause there to be misalignment if the sound
            // analysed is a chunk of a whole sound.
            value.resize(2 * harmRep.nHarms, tSize);
            // Fill harmRep.value with zeros
            value.setZero();
 
            return value;
        }
 
        // Corrected Frequencies indexed by Time and Frequency(nb_frame, nb_FrqCorrs)
        // This is a range of frequencies around f0 the harmonics of which are
        // searched for spectral peaks.
        Eigen::MatrixXf corrFreqsTF =
                (fundamentalFreqs.replicate(1, nFreqCorrs).array() + freqCorrs.transpose().replicate(tSize, 1).array())
                        .matrix();
 
        Eigen::MatrixXf inharmFactorsHI(harmRep.nHarms, nInharmCoeffs);
        Eigen::VectorXf columnVectorHarms = Eigen::VectorXf::LinSpaced(harmRep.nHarms, 1, harmRep.nHarms);
        Eigen::MatrixXf repeatedColumnVectorHarms = columnVectorHarms.replicate(1, nInharmCoeffs);
        Eigen::VectorXf squares = columnVectorHarms.array().square();
        Eigen::MatrixXf secondTerm = ((squares * inharmCoeffs.transpose()).array() + 1).array().sqrt();
        inharmFactorsHI = repeatedColumnVectorHarms.array() * secondTerm.array();
 
        // Reshape corrFreqsTF and inharmFactorsHI
        Eigen::MatrixXf reshapedCorrFreqsTF = corrFreqsTF.reshaped(tSize * nFreqCorrs, 1);
        Eigen::MatrixXf reshapedInharmFactorsHI = inharmFactorsHI.reshaped(1, harmRep.nHarms * nInharmCoeffs);
        Eigen::MatrixXf multiplied = reshapedCorrFreqsTF * reshapedInharmFactorsHI;
 
        // Reshape multiplied into multiple little 2-D matrices (4-D MATLAB simulation)
        std::vector<Eigen::MatrixXf> reshapedMultiplied =
                calculateReshaped2D(multiplied, tSize, nFreqCorrs, nInharmCoeffs);
        // Calculate fSupIdcsTFHI by applying transformation and clipping on submatrices of reshapedMultiplied
        std::vector<Eigen::MatrixXi> fSupIdcsTFHI =
                applyTransformation(reshapedMultiplied, stftObj->getFFTBinSize(), fSupport.size());
 
        // Calculate distrIdcsTFHI by repeating the tSize vector, multiplying by fftSize and adding to fSupIdcsTFHI
        std::vector<Eigen::MatrixXi> distrIdcsTFHI =
                calculateDistrIdcs(tSize, nFreqCorrs, nInharmCoeffs, fSupport.size(), fSupIdcsTFHI);
 
        // Extract subMatrices from distr based on indices
        std::vector<Eigen::MatrixXf> submatrices_cell =
                extractSubmatrices(/*TODO::Make this back to dist*/ distr, distrIdcsTFHI);
        // Compute totalErgTFI by summing matrices across the 3rd dimension
        std::vector<Eigen::MatrixXf> totalErgTFI =
                computeTotalErgTFI(submatrices_cell, tSize, nFreqCorrs, nInharmCoeffs);
 
        // Calculte scoreTI matrix
        Eigen::MatrixXf scoreTI(tSize, nInharmCoeffs);
        // Iterate over each 2D matrix in the vector
        for (int i = 0; i < totalErgTFI.size(); ++i)
        {
            scoreTI.col(i) = totalErgTFI[i].rowwise().maxCoeff();
        }
 
        // If the maximum score is within 1 % of the first score, just choose the first
        // score.That is, just choose the inharmonicity coefficient at index 1.
        auto [maxScoreTI, inharmCoeffIdcsT] = maxRowValues(scoreTI);
        // threshold it
        float threshold = 0.01f;
        for (int i = 0; i < maxScoreTI.size(); ++i)
        {
            if ((maxScoreTI(i) - scoreTI(i, 0)) / scoreTI(i, 0) <= threshold)
            {
                inharmCoeffIdcsT(i) = 0;
            }
        }
        // Basic time and freq vectors
        Eigen::VectorXi basicTimeIndices = Eigen::VectorXi::LinSpaced(tSize, 1, tSize);
        Eigen::VectorXi basicFreqCorrIndices = Eigen::VectorXi::LinSpaced(nFreqCorrs, 1, nFreqCorrs);
 
        Eigen::MatrixXi adjustedTimeIndices =
                (basicTimeIndices.array() + tSize * nFreqCorrs * (inharmCoeffIdcsT.array())).replicate(1, nFreqCorrs);
        Eigen::MatrixXi freqCorrRepeatIndices =
                ((basicFreqCorrIndices.array() - 1) * tSize).replicate(1, tSize).transpose();
        Eigen::MatrixXi combinedIndices = adjustedTimeIndices.array() + freqCorrRepeatIndices.array();
        Eigen::MatrixXi colIdcs =
                Eigen::Map<Eigen::MatrixXi>(combinedIndices.data(), combinedIndices.size(), 1).array() - 1;
 
        // Processing to populate total energy matrix using the flattened indices
        Eigen::MatrixXf totalErgTF(tSize, nFreqCorrs);
        for (int i = 0; i < colIdcs.size(); i++)
        {
            // Extract correct subMatrix
            int counter = std::floor(colIdcs(i) / (tSize * nFreqCorrs));
            auto& subMatrix = totalErgTFI[counter];
            // Extract element
            int linearIndex = colIdcs(i) % (tSize * nFreqCorrs);
            // Calculate row and column indices from linear index
            int row = linearIndex % subMatrix.rows();
            int col = linearIndex / subMatrix.rows();
 
            // Store the fetched value into fSupIdcsHTF
            totalErgTF(i) = subMatrix(row, col);
        }
 
        // Adjusted indices calculation for harmonic and frequency correction.
        Eigen::VectorXi harmonicAdjustedIndices =
                (tSize * nFreqCorrs * harmRep.nHarms * (inharmCoeffIdcsT.array())).array() + basicTimeIndices.array();
        std::vector<Eigen::MatrixXi> freqCorrectionMatrices(harmRep.nHarms,
                                                            harmonicAdjustedIndices.replicate(1, nFreqCorrs));
 
        std::vector<Eigen::MatrixXi> replicated_freq_indices(harmRep.nHarms);
        Eigen::MatrixXi freqCorrectionRepeatedIndices = basicFreqCorrIndices.transpose().replicate(tSize, 1);
 
        // Replicate freq_indices row-wise tSize times for each submatrix
        for (int i = 0; i < harmRep.nHarms; i++)
        {
            Eigen::MatrixXi submatrix = Eigen::MatrixXi::Zero(tSize, nFreqCorrs);
            for (int j = 0; j < tSize; j++)
            {
                submatrix.row(j) = basicFreqCorrIndices;
            }
            replicated_freq_indices[i] = submatrix;
        }
        // Generate indeces for harmonic repetitions
        Eigen::VectorXi harm_indices = Eigen::VectorXi::LinSpaced(harmRep.nHarms, 1, harmRep.nHarms);
        std::vector<Eigen::MatrixXi> replicated_harm_indices(harmRep.nHarms);
        // Replicate each element of harm_indices accordingly
        for (int i = 0; i < harmRep.nHarms; i++)
        {
            Eigen::MatrixXi submatrix = Eigen::MatrixXi::Zero(tSize, nFreqCorrs);
            for (int j = 0; j < tSize; j++)
            {
                submatrix.row(j).fill((harm_indices(i) - 1) *
                                      nFreqCorrs); // Fill each row with the current element of harm_indices
            }
            replicated_harm_indices[i] = submatrix;
        }
 
        // Combine frequency and harmonic indeces
        std::vector<Eigen::MatrixXi> combinedHarmonicFreqIndices(harmRep.nHarms);
        for (size_t i = 0; i < combinedHarmonicFreqIndices.size(); ++i)
        {
            // Add corresponding matrices element-wise
            Eigen::MatrixXi sum =
                    (replicated_harm_indices[i].array() + replicated_freq_indices[i].array() - 1) * tSize +
                    freqCorrectionMatrices[i].array();
            combinedHarmonicFreqIndices[i] = sum.array() - 1 /*for c++ indexing!*/;
        }
 
        // Reshape the matrices into a single column vector
        Eigen::VectorXi colIdcsNew(harmRep.nHarms * tSize * nFreqCorrs);
        int index = 0;
        for (size_t i = 0; i < replicated_harm_indices.size(); ++i)
        {
            auto subMatrix = combinedHarmonicFreqIndices[i];
            for (int i = 0; i < subMatrix.cols(); i++)
            {
                // Copy the current column into the vector
                colIdcsNew.segment(index, subMatrix.rows()) = subMatrix.col(i);
                index += subMatrix.rows(); // Increment index
            }
        }
 
        // Calculate fSupIdcsHTF
        std::vector<Eigen::MatrixXi> fSupIdcsHTF(harmRep.nHarms, Eigen::MatrixXi::Zero(tSize, nFreqCorrs));
        int setCounter = 0;
        for (int i = 0; i < colIdcsNew.size(); i++)
        {
            // Extract correct subMatrix
            int counter = std::floor(colIdcsNew(i) / (tSize * nFreqCorrs));
            auto& subMatrix = fSupIdcsTFHI[counter];
            // Extract element
            int linearIndex = colIdcsNew(i) % (tSize * nFreqCorrs);
            // Calculate row and column indices from linear index
            int row = linearIndex % subMatrix.rows();
            int col = linearIndex / subMatrix.rows();
 
            // Calculate correct subMatrix of fSupIdcsHTF to store in
            if (i % (tSize * nFreqCorrs) == 0 && i > 0)
                ++setCounter;
            // Store the fetched value into fSupIdcsHTF
            fSupIdcsHTF[setCounter](row, col) = subMatrix(row, col);
        }
        // Permute fSupIdcsHTF
        std::vector<Eigen::MatrixXi> permfSupIdcsHTF(nFreqCorrs, Eigen::MatrixXi::Zero(harmRep.nHarms, tSize));
        for (int c = 0; c < harmRep.nHarms; c++) // for each submatrix
        {
            for (int i = 0; i < tSize; i++) // for each row
            {
                for (int j = 0; j < nFreqCorrs; j++) // for each column
                {
                    // Retrieve element of original
                    auto elem = fSupIdcsHTF[c](i, j);
                    // Permute it
                    permfSupIdcsHTF[j](c, i) = elem;
                }
            }
        }
 
        // Find the indices of the maximum element in each column
        Eigen::VectorXi freqCorrIdcsT(totalErgTF.rows());
        for (int i = 0; i < totalErgTF.rows(); ++i)
        {
            float maxVal = totalErgTF.row(i).maxCoeff(); // Find the maximum value in the current row
            for (int j = 0; j < totalErgTF.cols(); ++j)
            {
                if (totalErgTF(i, j) == maxVal)
                {
                    freqCorrIdcsT(i) = j; // Store the index of the maximum element in maxIndices
                    break;                // Exit the loop once the maximum is found
                }
            }
        }
 
        // Calculate colIdcs
        Eigen::MatrixXi harmonicIndicesRepeated = harm_indices.replicate(1, tSize);
        Eigen::VectorXi freqCorrectionAdjustment = tSize * (freqCorrIdcsT.transpose().array()) - 1;
        Eigen::VectorXi combinedIndexCalculation = (basicTimeIndices + freqCorrectionAdjustment) * harmRep.nHarms;
        Eigen::MatrixXi combinedIndicesReplicated = combinedIndexCalculation.transpose().replicate(harmRep.nHarms, 1);
        Eigen::MatrixXi finalIndicesMatrix = harmonicIndicesRepeated + combinedIndicesReplicated;
        // Reshape the elements into a column vector
        Eigen::VectorXi columnVector =
                Eigen::Map<Eigen::VectorXi>(finalIndicesMatrix.data(), finalIndicesMatrix.size());
        Eigen::VectorXi colIdcsNewNew = columnVector.array() - 1; // for C++
 
        // Calculate partialFreqs
        Eigen::MatrixXf partialFreqs(harmRep.nHarms, tSize);
        partialFreqs.setZero();
        Eigen::VectorXi idxVec(colIdcsNewNew.size());
        for (int i = 0; i < colIdcsNewNew.size(); i++)
        {
            // Extract correct subMatrix
            int counter = std::floor(colIdcsNewNew(i) / (tSize * harmRep.nHarms));
            auto& subMatrix = permfSupIdcsHTF[counter];
            // Extract element
            int linearIndex = colIdcsNewNew(i) % (tSize * harmRep.nHarms);
            // Calculate row and column indices from linear index
            int row = linearIndex % subMatrix.rows();
            int col = linearIndex / subMatrix.rows();
 
            // Store the fetched value into fSupIdcsHTF
            idxVec(i) = subMatrix(row, col);
        }
 
        // Ensure the indices are within the valid range for fSupport
        jassert(idxVec.maxCoeff() < fSupport.size() && idxVec.minCoeff() >= 0);
 
        // Assign elements from fSupport to partialFreqs
        Eigen::VectorXf fSupportElems = idxVec.unaryExpr([&](int idx) { return fSupport(idx); });
        partialFreqs = Eigen::Map<Eigen::MatrixXf>(fSupportElems.data(), harmRep.nHarms, tSize);
 
        // Calculate partialAmps
        Eigen::VectorXi scaled_sequence = Eigen::VectorXi::LinSpaced(tSize, 0, tSize - 1)
                                                  .transpose()
                                                  .replicate(harmRep.nHarms, 1)
                                                  .reshaped(harmRep.nHarms * tSize, 1)
                                                  .array() *
                                          fSupport.size();
        Eigen::VectorXi summed_values = idxVec + scaled_sequence;
 
        // Access elements of distr using summed_values and store in partialAmps
        Eigen::MatrixXf partialAmps =
                summed_values.unaryExpr([&](int idx) { return distr(idx); }).reshaped(harmRep.nHarms, tSize);
 
        value.conservativeResize(partialFreqs.rows() + partialAmps.rows(), partialFreqs.cols());
        value << partialFreqs, partialAmps;
 
        return value;
    }
 
    void HarmonicRepresentation::replaceNanWithMedian(Eigen::VectorXf& estPitches)
    {
        // Create a vector to store non-NaN values
        std::vector<float> nonNanValues;
 
        // Find non-NaN values using Eigen's array operations
        nonNanValues.reserve(estPitches.size());
        for (int i = 0; i < estPitches.size(); ++i)
        {
            if (!std::isnan(estPitches(i)))
            {
                nonNanValues.push_back(estPitches(i));
            }
        }
 
        // Calculate median of non-NaN values
        float medianValue = 0.0f;
        if (!nonNanValues.empty())
        {
            std::sort(nonNanValues.begin(), nonNanValues.end());
            size_t middle = nonNanValues.size() / 2;
            if (nonNanValues.size() % 2 == 0)
            {
                medianValue = (nonNanValues[middle - 1] + nonNanValues[middle]) / 2.0f;
            }
            else
            {
                medianValue = nonNanValues[middle];
            }
        }
 
        // Replace NaN values in estPitches with median
        estPitches = estPitches.unaryExpr([medianValue](float v) { return std::isnan(v) ? medianValue : v; });
    }
 
    Eigen::VectorXf HarmonicRepresentation::fevalbp(const Eigen::MatrixXf& bp, const Eigen::VectorXf& x_v)
    {
        // Check for NaN values in bp and x_v
        jassert(!bp.array().isNaN().any());
        jassert(!x_v.array().isNaN().any());
 
        Eigen::VectorXf y_v(x_v.size());
 
        // Find indices where x_v is less than the first time in bp
        Eigen::VectorXi pos1 = findIndices(x_v, -std::numeric_limits<float>::infinity(), bp(0, 0));
        // Check if pos1 is not empty
        if (pos1.size() > 0)
        {
            // Set the values of y_v at indices in pos1 to the corresponding value from bp
            for (int i = 0; i < pos1.size(); ++i)
            {
                y_v(pos1(i)) = bp(0, 1);
            }
        }
 
        // Find indices where x_v is greater than the last time in bp
        Eigen::VectorXi pos2 = findIndices(x_v, bp(bp.rows() - 1, 0), std::numeric_limits<float>::infinity());
        // Check if pos1 is not empty
        if (pos2.size() > 0)
        {
            // Set the values of y_v at indices in pos1 to the corresponding value from bp
            for (int i = 0; i < pos2.size(); ++i)
            {
                y_v(pos2(i)) = bp(bp.rows() - 1, 1);
            }
        }
 
        // Find indices where x_v is within the range of times in bp
        Eigen::ArrayXi pos = findIndices(x_v, bp(0, 0), bp(bp.rows() - 1, 0));
 
        if (pos.size() > 1)
        {
            // If there are more than one index with the above property(which is probably
            // almost always the case, no ? ) estimate the pitch at the times in x_v with the
            // above property by interpolating on the(time, pitch) pairs in bp
            y_v(pos, 0) = LinearInterpolator::interpolate2DEigen(bp.col(0), bp.col(1), x_v(pos));
        }
        else
        {
            // Vector of length equal to the number of times at which we would like to have values
            y_v.setZero();
 
            for (int n = 0; n < x_v.size(); n++)
            {
                float x = x_v(n);
 
                // Calculate the differences between x and the elements of the first column of bp
                Eigen::VectorXf diff = bp.col(0).array() - x;
 
                // Find the index of the minimum absolute difference
                int minimum_pos;
                float minimum_value = diff.array().abs().minCoeff(&minimum_pos);
 
                // calculate the number of rows in bp
                // (this shouldn't need to be calculated every time)
                int L = bp.rows();
 
                // (Shouldn't L be only of 1 dimension?)
                if ((bp(minimum_pos, 0) == x) || (L == 1) || ((bp(minimum_pos, 0) < x) && (minimum_pos == L)) ||
                    ((bp(minimum_pos, 0) > x) && (minimum_pos == 0)))
                {
                    // If the closest time in bp is exactly equal to x OR
                    // if the number of rows in bp is equal to 1 OR
                    // if the closest time is less than x and its index is the last one
                    // in bp OR
                    // if the closest time is greater than x and its index is the first one in bp
                    // SET y_v AT THIS INDEX n TO THE PITCH AT THE TIME CLOSEST TO x IN bp
                    y_v(n) = bp(minimum_pos, 1);
                }
                else if (bp(minimum_pos, 0) < x)
                {
                    // Otherwise, if the time closest to x is less than x(but is not the last
                    // row - value in bp) linearly interpolate between the frequency at the index
                    // closest to x and the one afterward to obtain the frequency to put in y_v(n)
                    y_v(n) = (bp(minimum_pos + 1, 1) - bp(minimum_pos, 1)) /
                                     (bp(minimum_pos + 1, 0) - bp(minimum_pos, 0)) * (x - bp(minimum_pos, 0)) +
                             bp(minimum_pos, 1);
                }
                else if (bp(minimum_pos, 0) > x)
                {
                    // If the time closest to x is greater than x (but is not the first row-value in bp)
                    // linearly interpolate between the index one before the index closest to x and the one
                    // closest to x to obtain the frequency to put in y_v(n)
                    y_v(n) = (bp(minimum_pos, 1) - bp(minimum_pos - 1, 1)) /
                                     (bp(minimum_pos, 0) - bp(minimum_pos - 1, 0)) * (x - bp(minimum_pos - 1, 0)) +
                             bp(minimum_pos - 1, 1);
                }
                else
                {
                    // If none of the above cases match, throw an error
                    throw std::runtime_error("Not a valid case");
                }
            }
        }
 
        return y_v;
    }
 
    Eigen::VectorXi HarmonicRepresentation::findIndices(const Eigen::VectorXf& x_v, float threshMin, float threshMax)
    {
        // Adjust threshMin and threshMax for 'less than' and 'greater than' cases
        // If only one threshold is provided, the function works for 'less than' or 'greater than' depending on the
        // threshold
        Eigen::Array<bool, Eigen::Dynamic, 1> mask = (x_v.array() >= threshMin) && (x_v.array() <= threshMax);
 
        // Count true elements in the mask to allocate space for indices
        Eigen::VectorXi indices;
        Eigen::Index numIndices = mask.cast<int>().sum();
        indices.resize(numIndices);
 
        int index = 0;
        for (int i = 0; i < mask.size(); ++i)
        {
            if (mask(i))
            {
                indices(index++) = i;
            }
        }
 
        return indices;
    }
 
    std::vector<Eigen::MatrixXf> HarmonicRepresentation::calculateReshaped2D(const Eigen::MatrixXf& multiplied,
                                                                             int tSize, int nFreqCorrs,
                                                                             int nInharmCoeffs)
    {
        std::vector<Eigen::MatrixXf> reshapedMultiplied;
        reshapedMultiplied.reserve(harmRep.nHarms * nInharmCoeffs);
 
        int counter = 0;
        int expectedElements = tSize * nFreqCorrs * harmRep.nHarms * nInharmCoeffs;
 
        jassert(multiplied.size() == expectedElements); // Validity check for sufficient data in 'multiplied'
 
        for (int i = 0; i < harmRep.nHarms; ++i)
        {
            for (int j = 0; j < nInharmCoeffs; ++j)
            {
                // Extract elements from multiplied
                Eigen::MatrixXf c_mat =
                        Eigen::Map<const Eigen::MatrixXf>(multiplied.data() + counter, tSize, nFreqCorrs);
                // Store matrix in the vector
                reshapedMultiplied.push_back(c_mat);
                // Update counter
                counter += tSize * nFreqCorrs;
            }
        }
 
        return reshapedMultiplied;
    }
 
    std::vector<Eigen::MatrixXi> HarmonicRepresentation::applyTransformation(
            const std::vector<Eigen::MatrixXf>& reshapedMultiplied2_d, float binSize, int fSize)
    {
        std::vector<Eigen::MatrixXi> fSupIdcsTFHI;
 
        for (const auto& mat : reshapedMultiplied2_d)
        {
            // Apply the transformation to the current matrix
            Eigen::MatrixXf scaledArray = (1.0f / binSize) * mat.array();
            // Round each element of the scaled array
            Eigen::MatrixXf roundedArray = scaledArray.array().round();
            // Add 1 to each element of the rounded array
            Eigen::MatrixXi transformedMat = (roundedArray.array() /* + 1*/).cast<int>();
            // Clip values exceeding fSize
            transformedMat = transformedMat.array().min(fSize);
            // Store the transformed matrix
            fSupIdcsTFHI.push_back(transformedMat);
        }
 
        return fSupIdcsTFHI;
    }
 
    std::vector<Eigen::MatrixXi> HarmonicRepresentation::calculateDistrIdcs(
            int tSize, int nFreqCorrs, int nInharmCoeffs, int fSize, std::vector<Eigen::MatrixXi>& fSupIdcsTFHI2_d)
    {
        std::vector<Eigen::MatrixXi> s_cell;
        s_cell.reserve(harmRep.nHarms * nInharmCoeffs); // Pre-allocate memory
 
        // Part a) Distribute elements from (0:tSize-1) into smaller matrices
        for (int i = 0; i < harmRep.nHarms; ++i)
        {
            for (int j = 0; j < nInharmCoeffs; ++j)
            {
                // Create a matrix with dimensions (tSize, nFreqCorrs)
                Eigen::MatrixXi s_mat = Eigen::MatrixXi::Zero(tSize, nFreqCorrs);
                // Fill it with elements from (0:tSize-1)
                for (int k = 0; k < tSize; ++k)
                {
                    for (int l = 0; l < nFreqCorrs; ++l)
                    {
                        s_mat(k, l) = k * fSize;
                    }
                }
                // Store the matrix in the cell array
                s_cell.push_back(s_mat);
            }
        }
 
        // Part b) Add all matrices
        std::vector<Eigen::MatrixXi> distrIdcsTFHI;
        distrIdcsTFHI.reserve(fSupIdcsTFHI2_d.size()); // Pre-allocate memory
 
        for (size_t i = 0; i < fSupIdcsTFHI2_d.size(); ++i)
        {
            // Add matrices from cells element-wise
            Eigen::MatrixXi combinedMatrix = fSupIdcsTFHI2_d[i] + s_cell[i];
            distrIdcsTFHI.push_back(combinedMatrix);
        }
 
        return distrIdcsTFHI;
    }
 
    std::vector<Eigen::MatrixXf> HarmonicRepresentation::extractSubmatrices(
            const Eigen::MatrixXf& distr, const std::vector<Eigen::MatrixXi>& distrIdcsTFHI2_d)
    {
        std::vector<Eigen::MatrixXf> submatrices_cell;
 
        for (const auto& indices : distrIdcsTFHI2_d)
        {
            // Initialize a submatrix with the same dimensions as indices
            Eigen::MatrixXf submatrix(indices.rows(), indices.cols());
 
            // Iterate over each element in indices and access the corresponding element in distr
            for (int i = 0; i < indices.rows(); ++i)
            {
                for (int j = 0; j < indices.cols(); ++j)
                {
                    // Calculate the linear index from the row and column indices of indices
                    int linearIndex = indices(i, j);
 
                    // Access the corresponding element in distr and assign it to the submatrix
                    submatrix(i, j) = distr(linearIndex);
                }
            }
 
            // Store the submatrix in the result vector
            submatrices_cell.push_back(submatrix);
        }
 
        return submatrices_cell;
    }
 
    std::vector<Eigen::MatrixXf> HarmonicRepresentation::computeTotalErgTFI(
            const std::vector<Eigen::MatrixXf>& submatrices_cell, int tSize, int nFreqCorrs, int nInharmCoeffs)
    {
        std::vector<Eigen::MatrixXf> totalErgTFI(nInharmCoeffs);
 
        // Iterate over each pair of matrices in the cell
        for (int i = 0; i < nInharmCoeffs; ++i)
        {
            // Initialize a matrix to store the sum of the current pair
            Eigen::MatrixXf currentSum = Eigen::MatrixXf::Zero(tSize, nFreqCorrs);
            // Iterate over each matrix in the current pair
            for (int j = 0; j < harmRep.nHarms; ++j)
            {
                // Add the current matrix to the currentSum
                currentSum += submatrices_cell[i * harmRep.nHarms + j];
            }
            // Assign the currentSum matrix to the corresponding index in totalErgTFI2_d
            totalErgTFI[i] = currentSum;
        }
 
        return totalErgTFI;
    }
 
    std::pair<Eigen::VectorXf, Eigen::VectorXi> HarmonicRepresentation::maxRowValues(const Eigen::MatrixXf& matrix)
    {
        // Vector to store max values
        Eigen::VectorXf maxValues(matrix.rows());
        // Vector to store corresponding column indices
        Eigen::VectorXi maxIndices(matrix.rows());
 
        for (int row = 0; row < matrix.rows(); ++row)
        {
            float maxVal;
            int maxIdx;
            maxVal = matrix.row(row).maxCoeff(&maxIdx); // Get max value and index in current row
            maxValues(row) = maxVal;
            maxIndices(row) = maxIdx;
        }
 
        return {maxValues, maxIndices};
    }
 
} // namespace krotos