ERB_FFTSpectrogram.cpp

Go to the documentation of this file.

//=====================================================================================
//=====================================================================================
 
//=====================================================================================
krotos::ERB_FFTSpectrogram::ERB_FFTSpectrogram()
{
    // create the vector of N ERB-scale frequencies spaced between m_lo and m_hi
    // TODO: Tsovolas!
    m_cfArray = erbSpace(m_lo, m_hi, m_numChans);
}
 
//=====================================================================================
void krotos::ERB_FFTSpectrogram::setSampleRate(float fs)
{
    m_sampleRate = fs;
 
    // calculate hopSize in samples from sampling rate
    m_hopSizeSamples = static_cast<int>(round(m_sampleRate * m_hopSizeSec));
 
    // configure earFilter
    m_earFilter.configure(fs);
 
    // limit to 1/2 erb below Nyquist
    m_hi = std::min(m_hi, m_sampleRate / 2.0f - cferb(m_sampleRate / 2.0f) / 2.0f);
 
    // winSize param
    m_wSize = static_cast<int>(std::pow(2.0f, getNextPowerOfTwo(m_ERD * m_sampleRate * 2.0f)));
 
    // create gammaTone window
    m_gtWin = gtWindow(m_wSize, (float)m_wSize / (m_ERD * m_sampleRate), 4);
}
 
std::vector<std::vector<float>> krotos::ERB_FFTSpectrogram::filterSpectrum(std::vector<float>& inputSignal)
{
    // apply the mid outer ear filter to the signal
    std::vector<float> filteredSignal = applyEarFilter(inputSignal);
    int length = static_cast<int>(filteredSignal.size());
 
    // pad signal with zeros to align analysis point with window power centroid
    // square gammatone window
    std::vector<float> gtWinSquared;
    gtWinSquared.reserve(m_gtWin.size());
 
    std::transform(m_gtWin.begin(), m_gtWin.end(), std::back_inserter(gtWinSquared), [](float n) { return n * n; });
    // calculate centroid and spread - take ceil to provide zeropad offset
    int zpOffest = static_cast<int>(std::ceil(centroid(gtWinSquared) - 1.0f));
    // zeroPad the signal
    filteredSignal.insert(filteredSignal.begin(), zpOffest, 0.0f);
    // calculate last hopIndex where the window will fit within the end of the
    // signal
    int lastIndex = (length - (m_wSize - zpOffest)) / m_hopSizeSamples * m_hopSizeSamples + 1;
    std::vector<float> startSamples;
    for (int i = 0; i <= lastIndex; i += m_hopSizeSamples)
    {
        startSamples.push_back(static_cast<float>(i));
    }
 
    // build the filterBank
    buildFilterBank(static_cast<int>(startSamples.size()));
 
    // matrix of windowed slices of signal
    std::vector<std::vector<float>> fr(m_wSize, std::vector<float>(startSamples.size()));
 
    // obtain windowed segments of the signal
    for (int i = 0; i < startSamples.size(); ++i)
    {
        for (int j = 0; j < m_wSize; ++j)
        {
            fr[j][i] = filteredSignal[static_cast<int>(startSamples[i]) + j] * m_gtWin[j];
        }
    }
 
    // power spectrum
    std::vector<std::vector<float>> powerSpec = calculatePowerSpectrum(fr);
 
    std::vector<std::vector<float>> tmp(powerSpec[0].size(), std::vector<float>(m_wfunct[0].size(), 0.0f));
 
    for (size_t i = 0; i < m_wfunct[0].size(); ++i)
    {
        for (size_t j = 0; j < powerSpec[0].size(); ++j)
        {
            for (size_t k = 0; k < m_wfunct.size(); ++k)
            {
                tmp[j][i] += static_cast<float>(m_wfunct[k][i]) * powerSpec[k][j];
            }
            tmp[j][i] = std::pow(tmp[j][i], 0.25f);
        }
    }
 
    // write to file to import to MATLAB and compare against the timbretoolbox
    // implementation std::string filename = "matrix_data.txt";
    // writeMatrixToFile(tmp, filename);
 
    return tmp;
}
 
std::vector<float> krotos::ERB_FFTSpectrogram::applyEarFilter(std::vector<float>& inputSignal)
{
    jassert(m_sampleRate != -1.0f);
 
    // TODO: rewrite to not copy memory around !!
    std::vector<float> earFilteredSignal = m_earFilter.processSignal(inputSignal);
 
    return earFilteredSignal;
}
 
void krotos::ERB_FFTSpectrogram::buildFilterBank(int startSamplesSize)
{
    // array of kernel bandwidth coeffs
    std::vector<float> b = cferb(m_cfArray);
 
    // ERB to gammatone b parameter
    float scalar = 0.982f; // keep here for now
    transform(b.begin(), b.end(), b.begin(), [scalar](float value) { return value / scalar; });
 
    std::vector<float> bb(b.size());
    transform(b.begin(), b.end(), bb.begin(), [=](float x) {
        float diffSquared = x * x - m_b0 * m_b0;
        return std::sqrt(diffSquared);
    });
 
    // matrix of kernels(array of gammatone power tfs sampled at fft spectrum
    // frequencies).
    std::vector<std::vector<float>> fSupport(m_wSize / 2, std::vector<float>(m_numChans));
    for (int i = 0; i < m_wSize / 2; i++)
    {
        for (int j = 0; j < m_numChans; j++)
        {
            fSupport[i][j] = static_cast<float>(i + 1) * m_sampleRate / m_wSize;
        }
    }
 
    std::vector<std::vector<float>> cf(m_wSize / 2, std::vector<float>(m_cfArray.size()));
    for (int i = 0; i < m_wSize / 2; i++)
    {
        for (int j = 0; j < m_cfArray.size(); j++)
        {
            cf[i][j] = m_cfArray[j];
        }
    }
 
    std::vector<std::vector<float>> repeated_bb(m_wSize / 2, std::vector<float>(bb.size()));
    for (int i = 0; i < m_wSize / 2; ++i)
    {
        repeated_bb[i] = bb;
    }
 
    // power transfer functions- need to cast to double to avoid underflow -- DAMN
    // IT!!!
    size_t rows = fSupport.size();
    size_t cols = cf[0].size();
    m_wfunct.resize(rows, std::vector<double>(cols));
 
    for (size_t i = 0; i < rows; ++i)
    {
        for (size_t j = 0; j < cols; ++j)
        {
            double term = std::pow((1.0 / (startSamplesSize * (fSupport[i][j] - cf[i][j]) + repeated_bb[i][j])), 4.0);
            m_wfunct[i][j] = std::abs(std::pow(term, 2.0));
        }
    }
 
    // adjust so weight == ERB
    std::vector<double> adjustweight(m_cfArray.size());
    for (size_t j = 0; j < m_wfunct[0].size(); ++j)
    {
        double sum_wfunct = 0.0;
        for (size_t i = 0; i < m_wfunct.size(); ++i)
        {
            sum_wfunct += m_wfunct[i][j];
        }
 
        adjustweight[j] = cferb(static_cast<double>(m_cfArray[j])) / sum_wfunct;
    }
 
    // TODO::Ask if this would work:  Optimize it by using SIMD (Single
    // Instruction, Multiple Data) instructions to perform parallel computations on
    // multiple elements simultaneously. perhaps using: Intel's SIMD-oriented Fast
    // Mersenne Twister (SIMD-optimized PRNG) or Intel Intrinsics.??????????
    size_t adjustweight_size = adjustweight.size();
    for (size_t i = 0; i < m_wfunct.size(); ++i)
    {
        for (size_t j = 0; j < m_wfunct[0].size(); ++j)
        {
            m_wfunct[i][j] *= adjustweight[j % adjustweight_size];
        }
    }
 
    // find the maximum value in wfunct
    double maxVal = 0.0;
    for (const auto& row : m_wfunct)
    {
        auto maxElement = std::max_element(row.begin(), row.end());
        if (maxElement != row.end() && *maxElement > maxVal)
        {
            maxVal = *maxElement;
        }
    }
 
    // divide each element of wfunct by the maximum value
    for (auto& row : m_wfunct)
    {
        std::transform(row.begin(), row.end(), row.begin(), [maxVal](double val) { return val / maxVal; });
    }
}
 
std::vector<float> krotos::ERB_FFTSpectrogram::gtWindow(int n, float b, int order)
{
    std::vector<float> y(n);
 
    for (int i = 0; i < n; ++i)
    {
        float t = static_cast<float>(i) / n;
        y[i] = std::pow(b, (float)order) * std::pow(t, (float)order - 1.0f) * std::exp(-2.0f * float(M_PI) * b * t);
    }
 
    std::reverse(y.begin(), y.end());
 
    // normalize the window
    float maxVal = *std::max_element(y.begin(), y.end());
    for (auto& val : y)
    {
        val /= maxVal;
    }
 
    return y;
}
 
std::vector<float> krotos::ERB_FFTSpectrogram::cferb(std::vector<float> cf)
{
    std::vector<float> bw;
 
    for (float cfValue : cf)
    {
        cfValue = (24.7f * (1.0f + 4.37f * cfValue / 1000.0f));
        bw.push_back(cfValue);
    }
 
    return bw;
}
 
float krotos::ERB_FFTSpectrogram::cferb(const float cf)
{
    float bw = 24.7f * (1.0f + 4.37f * cf / 1000.0f);
    return bw;
}
 
double krotos::ERB_FFTSpectrogram::cferb(const double cf)
{
    double bw = 24.7 * (1.0 + 4.37 * cf / 1000.0);
    return bw;
}
 
std::vector<float> krotos::ERB_FFTSpectrogram::erbSpace(float low, float hi, int N)
{
    // change the following parameters if you wish to use a different
    // ERB scale. Must change in makeERBCoeffs too.
    float EarQ = 9.26449f; // Glasberg and Moore Parameters
    float minBW = 24.7f;
 
    std::vector<float> y(N);
 
    float a = EarQ * minBW;
    for (int i = 0; i < N; ++i)
    {
        float cf = -a + std::exp(i * (-std::log(hi + a) + std::log(low + a)) / (N - 1)) * (hi + a);
        y[N - i - 1] = cf;
    }
 
    return y;
}
 
float krotos::ERB_FFTSpectrogram::centroid(const std::vector<float>& x)
{
    std::size_t n = x.size();
    float sum = 0.0f;
    float totalWeight = 0.0f;
 
    // compute the weighted sum and total weight
    for (std::size_t i = 0; i < n; i++)
    {
        sum += static_cast<float>(i + 1) * x[i];
        totalWeight += x[i];
    }
 
    // normalize and return the centroid
    return (sum + std::numeric_limits<float>::epsilon()) / (totalWeight + std::numeric_limits<float>::epsilon());
}
 
std::vector<std::vector<float>> krotos::ERB_FFTSpectrogram::calculatePowerSpectrum(
        const std::vector<std::vector<float>>& inFrames)
{
    size_t numRows = inFrames.size();
    size_t numCols = inFrames[0].size();
    size_t fftSize = numRows;
 
    std::vector<std::vector<float>> powerSpectrum(numRows, std::vector<float>(numCols));
 
    // create the dsp::FFT object
    int fftOrder = static_cast<int>(log2(fftSize));
    dsp::FFT fft(fftOrder);
    m_fftSize = static_cast<int>(fftSize);
 
    // create input and output buffers for the FFT
    std::vector<float> inputBuffer(2 * fftSize, 0.0f);
 
    for (size_t col = 0; col < numCols; ++col)
    {
        // copy column data to the input buffer
        for (size_t row = 0; row < numRows; ++row)
        {
            inputBuffer[row] = inFrames[row][col];
        }
 
        // perform frequency-only forward transform
        fft.performFrequencyOnlyForwardTransform(inputBuffer.data());
 
        // get the magnitude squared of the FFT coefficients
        for (size_t row = 0; row < numRows; ++row)
        {
            float magnitudeSquared = std::norm(inputBuffer[row]);
            powerSpectrum[row][col] = magnitudeSquared;
        }
 
        // zero inputBuffer to prepare for next frame
        std::fill(inputBuffer.begin(), inputBuffer.end(), 0.0f);
    }
 
    return powerSpectrum;
}
 
void krotos::ERB_FFTSpectrogram::drawSpectrogram(juce::Image& image, const std::vector<std::vector<float>>& stftMatrix)
{
    auto numOfFrames = stftMatrix.size();
 
    auto imageWidth = image.getWidth();
    auto imageHeight = image.getHeight();
 
    jassert(imageWidth > 0 && imageHeight > 0);
 
    float factor = (static_cast<float>(numOfFrames)) / (static_cast<float>(imageWidth));
 
    for (int x = 0; x < imageWidth; x++)
    {
        int frameIndex = static_cast<int>(factor * static_cast<float>(x));
 
        jassert(frameIndex <= numOfFrames);
 
        for (auto y = 1; y < imageHeight; ++y)
        {
            auto maxLevel = juce::FloatVectorOperations::findMinAndMax(stftMatrix.at(frameIndex).data(), m_fftSize);
 
            // proportion of the image that this y pixel corresponds to ?
            auto skewedProportionY = 1.0f - std::exp(std::log((float)y / (float)imageHeight) * m_scalingConstant);
            // proportion to fft Index
            auto fftDataIndex = (size_t)juce::jlimit(
                    /*lowerLimit*/ 0,
                    /*upperLimit*/ m_fftSize,
                    /*value to constrain*/ (int)(skewedProportionY * m_fftSize));
 
            // maps from Magnitude of spectrum to 0...1
            auto level = juce::jmap(/*sourceValue*/ stftMatrix.at(frameIndex).at(fftDataIndex),
                                    /*sourceRangeMin*/ 0.0f,
                                    /*sourceRangeMax*/ juce::jmax(maxLevel.getEnd(), 1e-5f),
                                    /*targetRangeMin*/ 0.0f,
                                    /*targetRangeMax*/ 1.0f);
 
            image.setPixelAt(x, y, juce::Colour::fromHSV(level, 1.0f, level, 1.0f));
        }
    }
}
 
float krotos::ERB_FFTSpectrogram::getNextPowerOfTwo(float value)
{
    if (value <= 0.0f)
        return 0.0f;
 
    float exponent = std::ceilf(std::log2f(value));
 
    return exponent;
}
 
std::vector<float> krotos::ERB_FFTSpectrogram::erbs2Hz(const std::vector<float>& erbsVector)
{
    std::vector<float> hzVector(erbsVector.size());
 
    std::transform(erbsVector.begin(), erbsVector.end(), hzVector.begin(),
                   [](float hz) { return std::pow(2.0f, (hz / 6.44f) + 7.84f) - 229.0f; });
 
    return hzVector;
}
 
Eigen::VectorXf krotos::ERB_FFTSpectrogram::erbs2Hz(const Eigen::VectorXf& erbsVector)
{
    return Eigen::pow(2.0f, erbsVector.array() / 6.44f + 7.84f) - 229.0f;
}
 
std::vector<float> krotos::ERB_FFTSpectrogram::hz2Erbs(const std::vector<float>& hzVector)
{
    std::vector<float> erbsVector(hzVector.size());
 
    std::transform(hzVector.begin(), hzVector.end(), erbsVector.begin(),
                   [](float hz) { return 6.44f * (std::log2(229.0f + hz) - 7.84f); });
 
    return erbsVector;
}
 
float krotos::ERB_FFTSpectrogram::hz2Erbs(const float hz) { return 6.44f * (std::log2(229.0f + hz) - 7.84f); }
 
Eigen::VectorXf krotos::ERB_FFTSpectrogram::hz2Erbs(const Eigen::VectorXf& hzVector)
{
    return 6.44f * (Eigen::log2(229.0f + hzVector.array()) - 7.84f);
}