FrequencyDomainAnalysisFramework.cpp

Go to the documentation of this file.

//=====================================================================================
//=====================================================================================
 
#include "FrequencyDomainAnalysisFramework.h"
#include "Utilities.h"
 
//=====================================================================================
krotos::FrequencyDomainAnalysisFramework::FrequencyDomainAnalysisFramework(int analysisFrameSize, int sampleRate)
    : m_numMelCoefficients(40), m_numMelCepstralCoefficients(13), m_spectralFrameSize(analysisFrameSize)
{
    setSamplingFrequency(sampleRate);
 
    initialiseSpectralAnalysisAlgorithms();
}
 
//=====================================================================================
void krotos::FrequencyDomainAnalysisFramework::setSamplingFrequency(int fs)
{
    m_samplingFrequency = fs;
 
    // TODO:RD-52 that has to do with making sure that Mel Spectrum works @Chris
    // Scott
    m_melFrequencySpectrum.setSamplingFrequency(m_samplingFrequency);
}
 
//=====================================================================================
void krotos::FrequencyDomainAnalysisFramework::setNumberOfMelFrequencyCoefficients(int numcoeffs)
{
    jassert(m_spectralFrameSize > -1); // audio buffer size has not been set to a valid value
 
    m_numMelCoefficients = numcoeffs;
 
    m_melFrequencySpectrum.setParameters(m_spectralFrameSize, m_numMelCoefficients);
}
 
//=====================================================================================
void krotos::FrequencyDomainAnalysisFramework::setNumberOfMelFrequencyCepstralCoefficients(int numcoeffs)
{
    jassert(m_spectralFrameSize > -1); // audio buffer size has not been set to a valid value
 
    m_numMelCepstralCoefficients = numcoeffs;
 
    m_dct.setInputOutputSize(m_numMelCoefficients, m_numMelCepstralCoefficients);
}
 
//=====================================================================================
void krotos::FrequencyDomainAnalysisFramework::setFrequencyVector(const std::vector<float>& freqVector)
{
 
    if (!m_freqVector.empty())
    {
        m_freqVector.clear();
    }
 
    m_freqVector = freqVector;
    ;
}
 
//=====================================================================================
int krotos::FrequencyDomainAnalysisFramework::getSamplingFrequency()
{
    jassert(m_samplingFrequency > -1); // Sampling frequency has not been set to a valid value
    return m_samplingFrequency;
}
 
//=====================================================================================
int krotos::FrequencyDomainAnalysisFramework::getNumberOfMelFrequencyCepstralCoefficients()
{
    return m_numMelCepstralCoefficients;
}
 
//=====================================================================================
void krotos::FrequencyDomainAnalysisFramework::initialiseSpectralAnalysisAlgorithms()
{
    // Ensure your audio analysis buffer size is a positive power of two
    assert(isPositivePowerOfTwo(m_spectralFrameSize));
 
    // set prev magnitude spectrum for spectral flux
    m_prevMagSpecSpectralFlux.resize(m_spectralFrameSize);
 
    if (m_prevMagSpecSpectralFlux.size() > 0)
    {
        FloatVectorOperations::fill(&m_prevMagSpecSpectralFlux[0], 0.f, m_prevMagSpecSpectralFlux.size());
    }
 
    // set buffers for calculating complex spectral difference
    m_prevMagSpecCSD.resize(m_spectralFrameSize);
    m_prevPhase1CSD.resize(m_spectralFrameSize);
    m_prevPhase2CSD.resize(m_spectralFrameSize);
 
    if (m_prevMagSpecCSD.size() > 0)
    {
        FloatVectorOperations::fill(&m_prevMagSpecCSD[0], 0.f, m_prevMagSpecCSD.size());
    }
    if (m_prevPhase1CSD.size() > 0)
    {
        FloatVectorOperations::fill(&m_prevPhase1CSD[0], 0.f, m_prevPhase1CSD.size());
    }
    if (m_prevPhase2CSD.size() > 0)
    {
        FloatVectorOperations::fill(&m_prevPhase2CSD[0], 0.f, m_prevPhase2CSD.size());
    }
 
    setNumberOfMelFrequencyCoefficients(m_numMelCoefficients);
 
    setNumberOfMelFrequencyCepstralCoefficients(m_numMelCepstralCoefficients);
}
 
//=====================================================================================
void krotos::FrequencyDomainAnalysisFramework::setSpectralFrame(std::vector<float> spectralFrame)
{
    m_spectralFrame = spectralFrame;
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getFreqOfLoudestPartial()
{
    jassert(m_samplingFrequency > -1); // Sampling frequency has not been set to a valid value
 
    float currentMax{0.0f};
    int currentMaxIndex{0};
    for (int i = 0; i < m_spectralFrame.size(); i++)
    {
        if (std::fabsf(m_spectralFrame[i]) > currentMax)
        {
            currentMax = m_spectralFrame[i];
            currentMaxIndex = i;
        }
    }
 
    return (float)currentMaxIndex * (float)m_samplingFrequency / (float)m_spectralFrameSize;
}
 
//=====================================================================================
std::vector<float> krotos::FrequencyDomainAnalysisFramework::getMelSpectrum()
{
    return m_melFrequencySpectrum.calculateMelSpectrum(m_spectralFrame);
}
 
//=====================================================================================
std::vector<float> krotos::FrequencyDomainAnalysisFramework::getMFCC()
{
    // To-do: check if mel spectrum has already been computed
    auto melSpec = m_melFrequencySpectrum.calculateMelSpectrum(m_spectralFrame);
 
    // log transform
    for (auto& value : melSpec)
        value = std::logf(std::max(value, 1e-5f));
 
    return m_dct.calculateDCT(melSpec);
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getSpectralCentroid()
{
    jassert(m_samplingFrequency > -1); // Sampling frequency has not been set to a valid value
 
    float weightedSum = 0;
    float sum = 0;
 
    for (int i = 0; i < m_spectralFrame.size(); i++) // Careful :: never use m_specFrameSize they can be different!!!
    {
        float centreFrequency = m_freqVector[i];
 
        weightedSum += centreFrequency * m_spectralFrame[i];
 
        sum += m_spectralFrame[i];
    }
 
    if (sum == 0.0)
    {
        return 0.0;
    }
    else
    {
        return weightedSum / sum;
    }
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getSpectralFlatness()
{
    float logSum = 0;
    float sum = 0;
    const float eps = 2.2204e-16f; // MATLAB epsilon
 
    for (int i = 0; i < m_spectralFrame.size(); i++)
    {
        // we add 1 here to stop zeros making the logs -inf
        //@Chris Scott is gonna rectify the calculation here to reflect Python
        float sample = m_spectralFrame[i] + eps;
 
        logSum += logf(sample);
        sum += sample;
    }
 
    logSum = logSum / (float)m_spectralFrame.size();
    sum = sum / (float)m_spectralFrame.size();
 
    logSum = expf(logSum);
 
    if (sum == 0)
    {
        return 0.0;
    }
    else
    {
        return logSum / (sum + eps);
    }
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getSpectralSpread()
{
    jassert(m_samplingFrequency > -1); // Sampling frequency has not been set to a valid value
 
    float resolution = (float)m_samplingFrequency / (float)m_spectralFrameSize;
    float centroid = getSpectralCentroid(); // Calculate the spectral centroid
 
    float sum = 0;
    float magnitudeSum = 0;
 
    for (int i = 0; i < m_spectralFrame.size(); i++)
    {
        float centreFrequency = resolution * (float)i;
        float deviation = centreFrequency - centroid;      // Calculate the deviation from the centroid
        sum += deviation * deviation * m_spectralFrame[i]; // Accumulate the squared deviation multiplied by
                                                           // the magnitude
        magnitudeSum += m_spectralFrame[i];                // Accumulate the magnitudes
    }
 
    if (magnitudeSum == 0.0)
    {
        return 0.0;
    }
    else
    {
        float spread = std::sqrt(sum / magnitudeSum); // Calculate the square root of the
                                                      // average squared deviation
        return spread;
    }
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getSpectralSlope()
{
    jassert(m_samplingFrequency > -1); // Sampling frequency has not been set to a valid value
 
    float resolution = (float)m_samplingFrequency / (float)m_spectralFrameSize;
 
    float sumXY = 0;
    float sumX = 0;
    float sumY = 0;
    float sumX2 = 0;
 
    for (int i = 0; i < m_spectralFrame.size(); i++)
    {
        float centreFrequency = resolution * (float)i;
        float magnitude = m_spectralFrame[i];
 
        sumXY += centreFrequency * magnitude;
        sumX += centreFrequency;
        sumY += magnitude;
        sumX2 += centreFrequency * centreFrequency;
    }
 
    float n = (float)m_spectralFrame.size();
    float slope = (n * sumXY - sumX * sumY) / ((n * sumX2 - sumX * sumX) * sumY);
 
    return slope;
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getSpectralRolloff()
{
    // sum the spectrum
    float spectralSum = 0;
 
    for (int i = 0; i < m_spectralFrameSize; i++)
    {
        spectralSum += m_spectralFrame[i];
    }
 
    // calculate the rolloff threshold;
    float rolloffThreshold = spectralSum * 0.85f;
 
    // define the rolloff index
    float rolloffIndex = 0;
 
    // reset the spectral sum
    spectralSum = 0;
 
    for (int i = 0; i < m_spectralFrameSize; i++)
    {
        // sum the spectrum again
        spectralSum += m_spectralFrame[i];
 
        // if we exceed the rolloff threshold
        if (spectralSum > rolloffThreshold)
        {
            // store the index and break
            rolloffIndex = (float)i;
            break;
        }
    }
 
    // calclulate spectral rolloff normalised by the length of the magnitude
    // spectrum
    float spectralRolloff = rolloffIndex / (float)m_spectralFrameSize;
 
    return spectralRolloff;
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getSpectralFlux()
{
    float sum = 0;
 
    for (int i = 0; i < m_spectralFrameSize; i++)
    {
        float difference = fabs(m_spectralFrame[i] - m_prevMagSpecSpectralFlux[i]);
 
        sum += difference;
 
        m_prevMagSpecSpectralFlux[i] = m_spectralFrame[i];
    }
 
    return sum;
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getFrameEnergy()
{
    float energy = 0.0f;
 
    for (int i = 0; i < m_spectralFrameSize; ++i)
    {
        float amplitude = m_spectralFrame[i];
        energy += amplitude * amplitude;
    }
 
    return energy;
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getHighFreqencyContentOnsetDectionFunction()
{
    float sum = 0;
 
    for (int i = 0; i < m_spectralFrameSize; i++)
    {
        sum += (m_spectralFrame[i] * (float)i);
    }
 
    return sum;
}
 
//=====================================================================================
// float
// krotos::FrequencyDomainAnalysisFramework::getComplexSpectralDifference()
//{
//    checkFFTHasBeenCalculated();
//
//    float phaseDeviation;
//    float complexSpectralDifference;
//
//    complexSpectralDifference = 0;
//
//    for (int i = 0;i < m_spectralFrameSize;i++)
//    {
//        float phaseValue = atan2(m_fft.imag[i],m_fft.real[i]);
//
//        phaseDeviation = phaseValue - (2*m_prevPhase1CSD[i]) +
//        m_prevPhase2CSD[i];
//
//        //phaseDeviation = phaseWrap(phaseDeviation);
//
//        float magnitudeDifference = m_fft.magnitudeSpectrum[i] -
//        m_prevMagSpecCSD[i];
//
//        float phaseDifference =
//        -m_fft.magnitudeSpectrum[i]*sin(phaseDeviation);
//
//        complexSpectralDifference +=
//        sqrt((magnitudeDifference*magnitudeDifference) +
//        (phaseDifference*phaseDifference));
//
//        m_prevPhase2CSD[i] = m_prevPhase1CSD[i];
//        m_prevPhase1CSD[i] = phaseValue;
//        m_prevMagSpecCSD[i] = m_fft.magnitudeSpectrum[i];
//    }
//
//    return complexSpectralDifference;
//}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::calculateSpectralMoment(int n)
{
    float denominatorSum = 0;
    float numeratorSum = 0;
 
    for (int i = 0; i < m_spectralFrameSize; i++)
    {
        float k = (float)i;
        numeratorSum += powf(k, static_cast<float>(n)) * m_spectralFrame[i];
 
        denominatorSum += m_spectralFrame[i];
    }
 
    if (denominatorSum == 0)
    {
        return 0.0;
    }
    else
    {
        return numeratorSum / denominatorSum;
    }
}
 
//=====================================================================================
float krotos::FrequencyDomainAnalysisFramework::getSpectralKurtosis()
{
    float u1 = calculateSpectralMoment(1);
    float u2 = calculateSpectralMoment(2);
    float u3 = calculateSpectralMoment(3);
    float u4 = calculateSpectralMoment(4);
 
    float spread = sqrt(u2 - u1 * u1);
 
    float x = -3 * powf(u1, 4) + 6 * u1 * u2 + 4 * u1 * u3 + u4;
    float y = powf(spread, 4);
 
    float kurtosis;
 
    if (y == 0)
    {
        kurtosis = 0.0;
    }
    else
    {
        kurtosis = (x / y) - 3;
    }
 
    return kurtosis;
}