c@5: //
c@17: //  NoveltyCurveProcessor.cpp
c@5: //  Tempogram
c@5: //
c@5: //  Created by Carl Bussey on 10/07/2014.
c@5: //  Copyright (c) 2014 Carl Bussey. All rights reserved.
c@5: //
c@5: 
c@11: //Spectrogram dimensions should be flipped?
c@11: 
c@14: #include "NoveltyCurveProcessor.h"
c@5: using namespace std;
c@5: 
c@14: NoveltyCurveProcessor::NoveltyCurveProcessor(const float &samplingFrequency, const size_t &fftLength, const size_t &numberOfBlocks, const size_t &compressionConstant) :
c@5:     m_samplingFrequency(samplingFrequency),
c@7:     m_fftLength(fftLength),
c@7:     m_blockSize(fftLength/2 + 1),
c@5:     m_numberOfBlocks(numberOfBlocks),
c@5:     m_compressionConstant(compressionConstant),
c@5:     m_numberOfBands(5),
c@13:     m_pBandBoundaries(0),
c@13:     m_pBandSum(0)
c@5: {
c@5:     initialise();
c@5: }
c@5: 
c@14: NoveltyCurveProcessor::~NoveltyCurveProcessor(){
c@5:     cleanup();
c@5: }
c@5: 
c@9: //allocate all space and set variable
c@5: void
c@14: NoveltyCurveProcessor::initialise(){
c@5:     
c@13:     // for bandwise processing, the band is split into 5 bands. m_pBandBoundaries contains the upper and lower bin boundaries for each band.
c@13:     m_pBandBoundaries = new int[m_numberOfBands+1];
c@13:     m_pBandBoundaries[0] = 0;
c@13:     for (unsigned int band = 1; band < m_numberOfBands; band++){
c@13:         float lowFreq = 500*pow(2.5, (int)band-1);
c@13:         m_pBandBoundaries[band] = m_fftLength*lowFreq/m_samplingFrequency;
c@5:     }
c@13:     m_pBandBoundaries[m_numberOfBands] = m_blockSize;
c@5:     
c@13:     m_pBandSum = new float [m_numberOfBands];
c@5: }
c@5: 
c@9: //delete space allocated in initialise()
c@5: void
c@14: NoveltyCurveProcessor::cleanup(){
c@13:     delete []m_pBandBoundaries;
c@13:     m_pBandBoundaries = 0;
c@13:     delete []m_pBandSum;
c@13:     m_pBandSum = 0;
c@5: }
c@5: 
c@9: //calculate max of spectrogram
c@14: float NoveltyCurveProcessor::calculateMax(const vector< vector<float> > &spectrogram) const
c@13: {
c@5:     float max = 0;
c@5:     
c@13:     for (unsigned int j = 0; j < m_numberOfBlocks; j++){
c@13:         for (unsigned int i = 0; i < m_blockSize; i++){
c@7:             max = max > fabs(spectrogram[i][j]) ? max : fabs(spectrogram[i][j]);
c@5:         }
c@5:     }
c@5:     
c@5:     return max;
c@5: }
c@5: 
c@9: //subtract local average of novelty curve
c@9: //uses m_hannWindow as filter
c@14: void NoveltyCurveProcessor::subtractLocalAverage(vector<float> &noveltyCurve, const size_t &smoothLength) const
c@13: {
c@5:     vector<float> localAverage(m_numberOfBlocks);
c@5:     
c@13:     float * m_hannWindow = new float[smoothLength];
c@13:     WindowFunction::hanning(m_hannWindow, smoothLength, true);
c@9:     
c@13:     FIRFilter filter(m_numberOfBlocks, smoothLength);
c@15:     filter.process(&noveltyCurve[0], m_hannWindow, &localAverage[0], FIRFilter::middle);
c@5:     
c@7:     assert(noveltyCurve.size() == m_numberOfBlocks);
c@13:     for (unsigned int i = 0; i < m_numberOfBlocks; i++){
c@5:         noveltyCurve[i] -= localAverage[i];
c@5:         noveltyCurve[i] = noveltyCurve[i] >= 0 ? noveltyCurve[i] : 0;
c@5:     }
c@9:     
c@11:     delete []m_hannWindow;
c@13:     m_hannWindow = 0;
c@5: }
c@5: 
c@9: //smoothed differentiator filter. Flips upper half of hanning window about y-axis to create coefficients.
c@14: void NoveltyCurveProcessor::smoothedDifferentiator(vector< vector<float> > &spectrogram, const size_t &smoothLength) const
c@13: {
c@7:     
c@7:     float * diffHannWindow = new float [smoothLength];
c@7:     WindowFunction::hanning(diffHannWindow, smoothLength, true);
c@7:     
c@7:     if(smoothLength%2) diffHannWindow[(smoothLength+1)/2 - 1] = 0;
c@20:     for(int i = (smoothLength+1)/2; i < (int)smoothLength; i++){
c@7:         diffHannWindow[i] = -diffHannWindow[i];
c@7:     }
c@7:     
c@11:     FIRFilter smoothFilter(m_numberOfBlocks, smoothLength);
c@7:     
c@20:     for (int i = 0; i < (int)m_blockSize; i++){
c@15:         smoothFilter.process(&spectrogram[i][0], diffHannWindow, &spectrogram[i][0], FIRFilter::middle);
c@7:     }
c@7: }
c@7: 
c@9: //half rectification (set negative to zero)
c@14: void NoveltyCurveProcessor::halfWaveRectify(vector< vector<float> > &spectrogram) const
c@13: {
c@20:     for (int block = 0; block < (int)m_numberOfBlocks; block++){
c@20:         for (int k = 0; k < (int)m_blockSize; k++){
c@7:             if (spectrogram[k][block] < 0.0) spectrogram[k][block] = 0.0;
c@7:         }
c@7:     }
c@7: }
c@7: 
c@9: //process method
c@5: vector<float>
c@14: NoveltyCurveProcessor::spectrogramToNoveltyCurve(Spectrogram spectrogram) const
c@13: {
c@13:     std::vector<float> noveltyCurve(m_numberOfBlocks);
c@7:     
c@14:     //cout << spectrogram[0].size() << " : " << m_numberOfBlocks << endl;
c@7:     assert(spectrogram.size() == m_blockSize);
c@7:     assert(spectrogram[0].size() == m_numberOfBlocks);
c@5:     
c@9:     //normalise and log spectrogram
c@5:     float normaliseScale = calculateMax(spectrogram);
c@20:     for (int block = 0; block < (int)m_numberOfBlocks; block++){
c@20:         for (int k = 0; k < (int)m_blockSize; k++){
c@7:             if(normaliseScale != 0.0) spectrogram[k][block] /= normaliseScale; //normalise
c@7:             spectrogram[k][block] = log(1+m_compressionConstant*spectrogram[k][block]);
c@7:         }
c@7:     }
c@7: 
c@9:     //smooted differentiator
c@7:     smoothedDifferentiator(spectrogram, 5); //make smoothLength a parameter!
c@9:     //halfwave rectification
c@7:     halfWaveRectify(spectrogram);
c@7:     
c@9:     //bandwise processing
c@20:     for (int block = 0; block < (int)m_numberOfBlocks; block++){
c@20:         for (int band = 0; band < (int)m_numberOfBands; band++){
c@13:             int k = m_pBandBoundaries[band];
c@13:             int bandEnd = m_pBandBoundaries[band+1];
c@13:             m_pBandSum[band] = 0;
c@5:             
c@7:             while(k < bandEnd){
c@13:                 m_pBandSum[band] += spectrogram[k][block];
c@7:                 k++;
c@5:             }
c@5:         }
c@5:         float total = 0;
c@20:         for(int band = 0; band < (int)m_numberOfBands; band++){
c@13:             total += m_pBandSum[band];
c@5:         }
c@13:         noveltyCurve[block] = total/m_numberOfBands;
c@5:     }
c@5:     
c@9:     //subtract local averages
c@13:     subtractLocalAverage(noveltyCurve, 65);
c@5:     
c@13:     return noveltyCurve;
c@7: }