c@41: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@41: 
c@41: /*
c@67:  * SimilarityPlugin.cpp
c@41:  *
c@118:  * Copyright 2009 Centre for Digital Music, Queen Mary, University of London.
c@135: 
c@135:     This program is free software; you can redistribute it and/or
c@135:     modify it under the terms of the GNU General Public License as
c@135:     published by the Free Software Foundation; either version 2 of the
c@135:     License, or (at your option) any later version.  See the file
c@135:     COPYING included with this distribution for more information.
c@41:  */
c@41: 
c@41: #include <iostream>
c@44: #include <cstdio>
c@190: #include <algorithm>
c@41: 
c@41: #include "SimilarityPlugin.h"
c@42: #include "base/Pitch.h"
c@41: #include "dsp/mfcc/MFCC.h"
c@42: #include "dsp/chromagram/Chromagram.h"
c@41: #include "dsp/rateconversion/Decimator.h"
c@47: #include "dsp/rhythm/BeatSpectrum.h"
c@47: #include "maths/KLDivergence.h"
c@47: #include "maths/CosineDistance.h"
c@49: #include "maths/MathUtilities.h"
c@41: 
c@41: using std::string;
c@41: using std::vector;
c@41: using std::cerr;
c@41: using std::endl;
c@41: using std::ostringstream;
c@41: 
c@47: const float
c@47: SimilarityPlugin::m_noRhythm = 0.009;
c@47: 
c@47: const float
c@47: SimilarityPlugin::m_allRhythm = 0.991;
c@47: 
c@41: SimilarityPlugin::SimilarityPlugin(float inputSampleRate) :
c@41:     Plugin(inputSampleRate),
c@42:     m_type(TypeMFCC),
c@41:     m_mfcc(0),
c@47:     m_rhythmfcc(0),
c@42:     m_chromagram(0),
c@41:     m_decimator(0),
c@42:     m_featureColumnSize(20),
c@48:     m_rhythmWeighting(0.5f),
c@47:     m_rhythmClipDuration(4.f), // seconds
c@47:     m_rhythmClipOrigin(40.f), // seconds
c@47:     m_rhythmClipFrameSize(0),
c@47:     m_rhythmClipFrames(0),
c@47:     m_rhythmColumnSize(20),
c@41:     m_blockSize(0),
c@47:     m_channels(0),
c@47:     m_processRate(0),
c@47:     m_frameNo(0),
c@47:     m_done(false)
c@41: {
c@47:     int rate = lrintf(m_inputSampleRate);
c@47:     int internalRate = 22050;
c@47:     int decimationFactor = rate / internalRate;
c@47:     if (decimationFactor < 1) decimationFactor = 1;
c@47: 
c@47:     // must be a power of two
c@47:     while (decimationFactor & (decimationFactor - 1)) ++decimationFactor;
c@47: 
c@47:     m_processRate = rate / decimationFactor; // may be 22050, 24000 etc
c@41: }
c@41: 
c@41: SimilarityPlugin::~SimilarityPlugin()
c@41: {
c@41:     delete m_mfcc;
c@47:     delete m_rhythmfcc;
c@42:     delete m_chromagram;
c@41:     delete m_decimator;
c@41: }
c@41: 
c@41: string
c@41: SimilarityPlugin::getIdentifier() const
c@41: {
c@41:     return "qm-similarity";
c@41: }
c@41: 
c@41: string
c@41: SimilarityPlugin::getName() const
c@41: {
c@41:     return "Similarity";
c@41: }
c@41: 
c@41: string
c@41: SimilarityPlugin::getDescription() const
c@41: {
c@42:     return "Return a distance matrix for similarity between the input audio channels";
c@41: }
c@41: 
c@41: string
c@41: SimilarityPlugin::getMaker() const
c@41: {
c@50:     return "Queen Mary, University of London";
c@41: }
c@41: 
c@41: int
c@41: SimilarityPlugin::getPluginVersion() const
c@41: {
c@41:     return 1;
c@41: }
c@41: 
c@41: string
c@41: SimilarityPlugin::getCopyright() const
c@41: {
c@118:     return "Plugin by Mark Levy, Kurt Jacobson and Chris Cannam.  Copyright (c) 2009 QMUL - All Rights Reserved";
c@41: }
c@41: 
c@41: size_t
c@41: SimilarityPlugin::getMinChannelCount() const
c@41: {
c@43:     return 1;
c@41: }
c@41: 
c@41: size_t
c@41: SimilarityPlugin::getMaxChannelCount() const
c@41: {
c@41:     return 1024;
c@41: }
c@41: 
c@41: int
c@41: SimilarityPlugin::getDecimationFactor() const
c@41: {
c@41:     int rate = lrintf(m_inputSampleRate);
c@47:     return rate / m_processRate;
c@41: }
c@41: 
c@41: size_t
c@41: SimilarityPlugin::getPreferredStepSize() const
c@41: {
c@42:     if (m_blockSize == 0) calculateBlockSize();
c@47: 
c@47:     // there is also an assumption to this effect in process()
c@47:     // (referring to m_fftSize/2 instead of a literal post-decimation
c@47:     // step size):
c@45:     return m_blockSize/2;
c@41: }
c@41: 
c@41: size_t
c@41: SimilarityPlugin::getPreferredBlockSize() const
c@41: {
c@42:     if (m_blockSize == 0) calculateBlockSize();
c@42:     return m_blockSize;
c@42: }
c@42: 
c@42: void
c@42: SimilarityPlugin::calculateBlockSize() const
c@42: {
c@42:     if (m_blockSize != 0) return;
c@42:     int decimationFactor = getDecimationFactor();
c@66:     m_blockSize = 2048 * decimationFactor;
c@41: }
c@41: 
c@41: SimilarityPlugin::ParameterList SimilarityPlugin::getParameterDescriptors() const
c@41: {
c@41:     ParameterList list;
c@42: 
c@42:     ParameterDescriptor desc;
c@42:     desc.identifier = "featureType";
c@42:     desc.name = "Feature Type";
c@48:     desc.description = "Audio feature used for similarity measure.  Timbral: use the first 20 MFCCs (19 plus C0).  Chromatic: use 12 bin-per-octave chroma.  Rhythmic: compare beat spectra of short regions.";
c@42:     desc.unit = "";
c@42:     desc.minValue = 0;
c@48:     desc.maxValue = 4;
c@48:     desc.defaultValue = 1;
c@42:     desc.isQuantized = true;
c@42:     desc.quantizeStep = 1;
c@48:     desc.valueNames.push_back("Timbre");
c@48:     desc.valueNames.push_back("Timbre and Rhythm");
c@48:     desc.valueNames.push_back("Chroma");
c@48:     desc.valueNames.push_back("Chroma and Rhythm");
c@48:     desc.valueNames.push_back("Rhythm only");
c@42:     list.push_back(desc);	
c@48: /*
c@47:     desc.identifier = "rhythmWeighting";
c@47:     desc.name = "Influence of Rhythm";
c@47:     desc.description = "Proportion of similarity measure made up from rhythmic similarity component, from 0 (entirely timbral or chromatic) to 100 (entirely rhythmic).";
c@47:     desc.unit = "%";
c@47:     desc.minValue = 0;
c@47:     desc.maxValue = 100;
c@47:     desc.defaultValue = 0;
c@48:     desc.isQuantized = false;
c@47:     desc.valueNames.clear();
c@47:     list.push_back(desc);	
c@48: */
c@41:     return list;
c@41: }
c@41: 
c@41: float
c@41: SimilarityPlugin::getParameter(std::string param) const
c@41: {
c@42:     if (param == "featureType") {
c@48: 
c@48:         if (m_rhythmWeighting > m_allRhythm) {
c@48:             return 4;
c@48:         }
c@48: 
c@48:         switch (m_type) {
c@48: 
c@48:         case TypeMFCC:
c@48:             if (m_rhythmWeighting < m_noRhythm) return 0;
c@48:             else return 1;
c@48:             break;
c@48: 
c@48:         case TypeChroma:
c@48:             if (m_rhythmWeighting < m_noRhythm) return 2;
c@48:             else return 3;
c@48:             break;
c@48:         }            
c@48: 
c@48:         return 1;
c@48: 
c@48: //    } else if (param == "rhythmWeighting") {
c@48: //        return nearbyint(m_rhythmWeighting * 100.0);
c@42:     }
c@42: 
c@41:     std::cerr << "WARNING: SimilarityPlugin::getParameter: unknown parameter \""
c@41:               << param << "\"" << std::endl;
c@41:     return 0.0;
c@41: }
c@41: 
c@41: void
c@41: SimilarityPlugin::setParameter(std::string param, float value)
c@41: {
c@42:     if (param == "featureType") {
c@48: 
c@42:         int v = int(value + 0.1);
c@48: 
c@48:         Type newType = m_type;
c@48: 
c@48:         switch (v) {
c@48:         case 0: newType = TypeMFCC; m_rhythmWeighting = 0.0f; break;
c@48:         case 1: newType = TypeMFCC; m_rhythmWeighting = 0.5f; break;
c@48:         case 2: newType = TypeChroma; m_rhythmWeighting = 0.0f; break;
c@48:         case 3: newType = TypeChroma; m_rhythmWeighting = 0.5f; break;
c@48:         case 4: newType = TypeMFCC; m_rhythmWeighting = 1.f; break;
c@48:         }
c@48: 
c@48:         if (newType != m_type) m_blockSize = 0;
c@48: 
c@48:         m_type = newType;
c@42:         return;
c@48: 
c@48: //    } else if (param == "rhythmWeighting") {
c@48: //        m_rhythmWeighting = value / 100;
c@48: //        return;
c@42:     }
c@42: 
c@41:     std::cerr << "WARNING: SimilarityPlugin::setParameter: unknown parameter \""
c@41:               << param << "\"" << std::endl;
c@41: }
c@41: 
c@41: SimilarityPlugin::OutputList
c@41: SimilarityPlugin::getOutputDescriptors() const
c@41: {
c@41:     OutputList list;
c@41: 	
c@41:     OutputDescriptor similarity;
c@43:     similarity.identifier = "distancematrix";
c@43:     similarity.name = "Distance Matrix";
c@43:     similarity.description = "Distance matrix for similarity metric.  Smaller = more similar.  Should be symmetrical.";
c@41:     similarity.unit = "";
c@41:     similarity.hasFixedBinCount = true;
c@41:     similarity.binCount = m_channels;
c@41:     similarity.hasKnownExtents = false;
c@41:     similarity.isQuantized = false;
c@41:     similarity.sampleType = OutputDescriptor::FixedSampleRate;
c@41:     similarity.sampleRate = 1;
c@41: 	
c@43:     m_distanceMatrixOutput = list.size();
c@41:     list.push_back(similarity);
c@41: 	
c@43:     OutputDescriptor simvec;
c@43:     simvec.identifier = "distancevector";
c@43:     simvec.name = "Distance from First Channel";
c@43:     simvec.description = "Distance vector for similarity of each channel to the first channel.  Smaller = more similar.";
c@43:     simvec.unit = "";
c@43:     simvec.hasFixedBinCount = true;
c@43:     simvec.binCount = m_channels;
c@43:     simvec.hasKnownExtents = false;
c@43:     simvec.isQuantized = false;
c@43:     simvec.sampleType = OutputDescriptor::FixedSampleRate;
c@43:     simvec.sampleRate = 1;
c@43: 	
c@43:     m_distanceVectorOutput = list.size();
c@43:     list.push_back(simvec);
c@43: 	
c@44:     OutputDescriptor sortvec;
c@44:     sortvec.identifier = "sorteddistancevector";
c@44:     sortvec.name = "Ordered Distances from First Channel";
c@44:     sortvec.description = "Vector of the order of other channels in similarity to the first, followed by distance vector for similarity of each to the first.  Smaller = more similar.";
c@44:     sortvec.unit = "";
c@44:     sortvec.hasFixedBinCount = true;
c@44:     sortvec.binCount = m_channels;
c@44:     sortvec.hasKnownExtents = false;
c@44:     sortvec.isQuantized = false;
c@44:     sortvec.sampleType = OutputDescriptor::FixedSampleRate;
c@44:     sortvec.sampleRate = 1;
c@44: 	
c@44:     m_sortedVectorOutput = list.size();
c@44:     list.push_back(sortvec);
c@44: 	
c@41:     OutputDescriptor means;
c@41:     means.identifier = "means";
c@42:     means.name = "Feature Means";
c@43:     means.description = "Means of the feature bins.  Feature time (sec) corresponds to input channel.  Number of bins depends on selected feature type.";
c@41:     means.unit = "";
c@41:     means.hasFixedBinCount = true;
c@43:     means.binCount = m_featureColumnSize;
c@41:     means.hasKnownExtents = false;
c@41:     means.isQuantized = false;
c@43:     means.sampleType = OutputDescriptor::FixedSampleRate;
c@43:     means.sampleRate = 1;
c@41: 	
c@43:     m_meansOutput = list.size();
c@41:     list.push_back(means);
c@41:     
c@41:     OutputDescriptor variances;
c@41:     variances.identifier = "variances";
c@42:     variances.name = "Feature Variances";
c@43:     variances.description = "Variances of the feature bins.  Feature time (sec) corresponds to input channel.  Number of bins depends on selected feature type.";
c@41:     variances.unit = "";
c@41:     variances.hasFixedBinCount = true;
c@43:     variances.binCount = m_featureColumnSize;
c@41:     variances.hasKnownExtents = false;
c@41:     variances.isQuantized = false;
c@43:     variances.sampleType = OutputDescriptor::FixedSampleRate;
c@43:     variances.sampleRate = 1;
c@41: 	
c@43:     m_variancesOutput = list.size();
c@41:     list.push_back(variances);
c@41:     
c@47:     OutputDescriptor beatspectrum;
c@47:     beatspectrum.identifier = "beatspectrum";
c@47:     beatspectrum.name = "Beat Spectra";
c@47:     beatspectrum.description = "Rhythmic self-similarity vectors (beat spectra) for the input channels.  Feature time (sec) corresponds to input channel.  Not returned if rhythm weighting is zero.";
c@47:     beatspectrum.unit = "";
c@47:     if (m_rhythmClipFrames > 0) {
c@47:         beatspectrum.hasFixedBinCount = true;
c@47:         beatspectrum.binCount = m_rhythmClipFrames / 2;
c@47:     } else {
c@47:         beatspectrum.hasFixedBinCount = false;
c@47:     }
c@47:     beatspectrum.hasKnownExtents = false;
c@47:     beatspectrum.isQuantized = false;
c@47:     beatspectrum.sampleType = OutputDescriptor::FixedSampleRate;
c@47:     beatspectrum.sampleRate = 1;
c@47: 	
c@47:     m_beatSpectraOutput = list.size();
c@47:     list.push_back(beatspectrum);
c@47:     
c@41:     return list;
c@41: }
c@41: 
c@68: bool
c@68: SimilarityPlugin::initialise(size_t channels, size_t stepSize, size_t blockSize)
c@68: {
c@68:     if (channels < getMinChannelCount()) return false;
c@68: 
c@68:     // Using more than getMaxChannelCount is not actually a problem
c@68:     // for us.  Using "incorrect" step and block sizes would be fine
c@68:     // for timbral or chroma similarity, but will break rhythmic
c@68:     // similarity, so we'd better enforce these.
c@68: 
c@68:     if (stepSize != getPreferredStepSize()) {
c@68:         std::cerr << "SimilarityPlugin::initialise: supplied step size "
c@68:                   << stepSize << " differs from required step size "
c@68:                   << getPreferredStepSize() << std::endl;
c@68:         return false;
c@68:     }
c@68: 
c@68:     if (blockSize != getPreferredBlockSize()) {
c@68:         std::cerr << "SimilarityPlugin::initialise: supplied block size "
c@68:                   << blockSize << " differs from required block size "
c@68:                   << getPreferredBlockSize() << std::endl;
c@68:         return false;
c@68:     }        
c@68:     
c@68:     m_blockSize = blockSize;
c@68:     m_channels = channels;
c@68: 
c@68:     m_lastNonEmptyFrame = std::vector<int>(m_channels);
c@68:     for (int i = 0; i < m_channels; ++i) m_lastNonEmptyFrame[i] = -1;
c@68: 
c@68:     m_emptyFrameCount = std::vector<int>(m_channels);
c@68:     for (int i = 0; i < m_channels; ++i) m_emptyFrameCount[i] = 0;
c@68: 
c@68:     m_frameNo = 0;
c@68: 
c@68:     int decimationFactor = getDecimationFactor();
c@68:     if (decimationFactor > 1) {
c@68:         m_decimator = new Decimator(m_blockSize, decimationFactor);
c@68:     }
c@68: 
c@68:     if (m_type == TypeMFCC) {
c@68: 
c@68:         m_featureColumnSize = 20;
c@68: 
c@68:         MFCCConfig config(m_processRate);
c@68:         config.fftsize = 2048;
c@68:         config.nceps = m_featureColumnSize - 1;
c@68:         config.want_c0 = true;
c@68:         config.logpower = 1;
c@68:         m_mfcc = new MFCC(config);
c@68:         m_fftSize = m_mfcc->getfftlength();
c@68:         m_rhythmClipFrameSize = m_fftSize / 4;
c@68: 
c@68: //        std::cerr << "MFCC FS = " << config.FS << ", FFT size = " << m_fftSize<< std::endl;
c@68: 
c@68:     } else if (m_type == TypeChroma) {
c@68: 
c@68:         m_featureColumnSize = 12;
c@68: 
c@68:         // For simplicity, aim to have the chroma fft size equal to
c@68:         // 2048, the same as the mfcc fft size (so the input block
c@68:         // size does not depend on the feature type and we can use the
c@68:         // same processing parameters for rhythm etc).  This is also
c@68:         // why getPreferredBlockSize can confidently return 2048 * the
c@68:         // decimation factor.
c@68:         
c@68:         // The fft size for a chromagram is the filterbank Q value
c@68:         // times the sample rate, divided by the minimum frequency,
c@68:         // rounded up to the nearest power of two.
c@68: 
c@68:         double q = 1.0 / (pow(2.0, (1.0 / 12.0)) - 1.0);
c@68:         double fmin = (q * m_processRate) / 2048.0;
c@68: 
c@68:         // Round fmin up to the nearest MIDI pitch multiple of 12.
c@68:         // So long as fmin is greater than 12 to start with, this
c@68:         // should not change the resulting fft size.
c@68: 
c@68:         int pmin = Pitch::getPitchForFrequency(float(fmin));
c@68:         pmin = ((pmin / 12) + 1) * 12;
c@68:         fmin = Pitch::getFrequencyForPitch(pmin);
c@68: 
c@68:         float fmax = Pitch::getFrequencyForPitch(pmin + 36);
c@68: 
c@68:         ChromaConfig config;
c@68:         config.FS = m_processRate;
c@68:         config.min = fmin;
c@68:         config.max = fmax;
c@68:         config.BPO = 12;
c@68:         config.CQThresh = 0.0054;
c@68:         // We don't normalise the chromagram's columns individually;
c@68:         // we normalise the mean at the end instead
c@68:         config.normalise = MathUtilities::NormaliseNone;
c@68:         m_chromagram = new Chromagram(config);
c@68:         m_fftSize = m_chromagram->getFrameSize();
c@68:         
c@68:         if (m_fftSize != 2048) {
c@68:             std::cerr << "WARNING: SimilarityPlugin::initialise: Internal processing FFT size " << m_fftSize << " != expected size 2048 in chroma mode" << std::endl;
c@68:         }
c@68: 
c@68: //        std::cerr << "fftsize = " << m_fftSize << std::endl;
c@68: 
c@68:         m_rhythmClipFrameSize = m_fftSize / 4;
c@68: 
c@68: //        std::cerr << "m_rhythmClipFrameSize = " << m_rhythmClipFrameSize << std::endl;
c@68: //        std::cerr << "min = "<< config.min << ", max = " << config.max << std::endl;
c@68: 
c@68:     } else {
c@68: 
c@68:         std::cerr << "SimilarityPlugin::initialise: internal error: unknown type " << m_type << std::endl;
c@68:         return false;
c@68:     }
c@68:     
c@68:     if (needRhythm()) {
c@68:         m_rhythmClipFrames =
c@68:             int(ceil((m_rhythmClipDuration * m_processRate) 
c@68:                      / m_rhythmClipFrameSize));
c@68: //        std::cerr << "SimilarityPlugin::initialise: rhythm clip requires "
c@68: //                  << m_rhythmClipFrames << " frames of size "
c@68: //                  << m_rhythmClipFrameSize << " at process rate "
c@68: //                  << m_processRate << " ( = "
c@68: //                  << (float(m_rhythmClipFrames * m_rhythmClipFrameSize) / m_processRate) << " sec )"
c@68: //                  << std::endl;
c@68: 
c@68:         MFCCConfig config(m_processRate);
c@68:         config.fftsize = m_rhythmClipFrameSize;
c@68:         config.nceps = m_rhythmColumnSize - 1;
c@68:         config.want_c0 = true;
c@68:         config.logpower = 1;
c@68:         config.window = RectangularWindow; // because no overlap
c@68:         m_rhythmfcc = new MFCC(config);
c@68:     }
c@68: 
c@68:     for (int i = 0; i < m_channels; ++i) {
c@68: 
c@68:         m_values.push_back(FeatureMatrix());
c@68: 
c@68:         if (needRhythm()) {
c@68:             m_rhythmValues.push_back(FeatureColumnQueue());
c@68:         }
c@68:     }
c@68: 
c@68:     m_done = false;
c@68: 
c@68:     return true;
c@68: }
c@68: 
c@68: void
c@68: SimilarityPlugin::reset()
c@68: {
c@178:     for (int i = 0; i < int(m_values.size()); ++i) {
c@68:         m_values[i].clear();
c@68:     }
c@68: 
c@178:     for (int i = 0; i < int(m_rhythmValues.size()); ++i) {
c@68:         m_rhythmValues[i].clear();
c@68:     }
c@68: 
c@178:     for (int i = 0; i < int(m_lastNonEmptyFrame.size()); ++i) {
c@68:         m_lastNonEmptyFrame[i] = -1;
c@68:     }
c@68: 
c@178:     for (int i = 0; i < int(m_emptyFrameCount.size()); ++i) {
c@68:         m_emptyFrameCount[i] = 0;
c@68:     }
c@68: 
c@68:     m_done = false;
c@68: }
c@68: 
c@41: SimilarityPlugin::FeatureSet
c@41: SimilarityPlugin::process(const float *const *inputBuffers, Vamp::RealTime /* timestamp */)
c@41: {
c@47:     if (m_done) {
c@47:         return FeatureSet();
c@47:     }
c@47: 
c@41:     double *dblbuf = new double[m_blockSize];
c@41:     double *decbuf = dblbuf;
c@42:     if (m_decimator) decbuf = new double[m_fftSize];
c@42: 
c@47:     double *raw = new double[std::max(m_featureColumnSize,
c@47:                                       m_rhythmColumnSize)];
c@41: 
c@43:     float threshold = 1e-10;
c@43: 
c@47:     bool someRhythmFrameNeeded = false;
c@47: 
c@178:     for (int c = 0; c < m_channels; ++c) {
c@41: 
c@43:         bool empty = true;
c@43: 
c@41:         for (int i = 0; i < m_blockSize; ++i) {
c@43:             float val = inputBuffers[c][i];
c@43:             if (fabs(val) > threshold) empty = false;
c@43:             dblbuf[i] = val;
c@41:         }
c@41: 
c@47:         if (empty) {
c@47:             if (needRhythm() && ((m_frameNo % 2) == 0)) {
c@47:                 for (int i = 0; i < m_fftSize / m_rhythmClipFrameSize; ++i) {
c@178:                     if (int(m_rhythmValues[c].size()) < m_rhythmClipFrames) {
c@47:                         FeatureColumn mf(m_rhythmColumnSize);
c@47:                         for (int i = 0; i < m_rhythmColumnSize; ++i) {
c@47:                             mf[i] = 0.0;
c@47:                         }
c@47:                         m_rhythmValues[c].push_back(mf);
c@47:                     }
c@47:                 }
c@47:             }
c@60:             m_emptyFrameCount[c]++;
c@47:             continue;
c@47:         }
c@47: 
c@44:         m_lastNonEmptyFrame[c] = m_frameNo;
c@43: 
c@41:         if (m_decimator) {
c@41:             m_decimator->process(dblbuf, decbuf);
c@41:         }
c@42: 
c@47:         if (needTimbre()) {
c@47: 
c@66:             FeatureColumn mf(m_featureColumnSize);
c@66: 
c@47:             if (m_type == TypeMFCC) {
c@47:                 m_mfcc->process(decbuf, raw);
c@66:                 for (int i = 0; i < m_featureColumnSize; ++i) {
c@66:                     mf[i] = raw[i];
c@66:                 }
c@47:             } else if (m_type == TypeChroma) {
c@66:                 double *chroma = m_chromagram->process(decbuf);
c@66:                 for (int i = 0; i < m_featureColumnSize; ++i) {
c@66:                     mf[i] = chroma[i];
c@66:                 }
c@47:             }
c@41:         
c@47:             m_values[c].push_back(mf);
c@44:         }
c@41: 
c@47: //        std::cerr << "needRhythm = " << needRhythm() << ", frame = " << m_frameNo << std::endl;
c@47: 
c@47:         if (needRhythm() && ((m_frameNo % 2) == 0)) {
c@47: 
c@47:             // The incoming frames are overlapping; we only use every
c@47:             // other one, because we don't want the overlap (it would
c@47:             // screw up the rhythm)
c@47: 
c@47:             int frameOffset = 0;
c@47: 
c@47:             while (frameOffset + m_rhythmClipFrameSize <= m_fftSize) {
c@47: 
c@47:                 bool needRhythmFrame = true;
c@47: 
c@178:                 if (int(m_rhythmValues[c].size()) >= m_rhythmClipFrames) {
c@47: 
c@47:                     needRhythmFrame = false;
c@47: 
c@47:                     // assumes hopsize = framesize/2
c@47:                     float current = m_frameNo * (m_fftSize/2) + frameOffset;
c@47:                     current = current / m_processRate;
c@47:                     if (current - m_rhythmClipDuration < m_rhythmClipOrigin) {
c@47:                         needRhythmFrame = true;
c@47:                         m_rhythmValues[c].pop_front();
c@47:                     }
c@47: 
c@53: //                    if (needRhythmFrame) {
c@53: //                        std::cerr << "at current = " <<current << " (frame = " << m_frameNo << "), have " << m_rhythmValues[c].size() << ", need rhythm = " << needRhythmFrame << std::endl;
c@53: //                    }
c@47: 
c@47:                 }
c@47:                 
c@47:                 if (needRhythmFrame) {
c@47: 
c@47:                     someRhythmFrameNeeded = true;
c@47: 
c@47:                     m_rhythmfcc->process(decbuf + frameOffset, raw);
c@47:                     
c@47:                     FeatureColumn mf(m_rhythmColumnSize);
c@47:                     for (int i = 0; i < m_rhythmColumnSize; ++i) {
c@47:                         mf[i] = raw[i];
c@47:                     }
c@47: 
c@47:                     m_rhythmValues[c].push_back(mf);
c@47:                 }
c@47: 
c@47:                 frameOffset += m_rhythmClipFrameSize;
c@47:             }
c@47:         }
c@47:     }
c@47: 
c@47:     if (!needTimbre() && !someRhythmFrameNeeded && ((m_frameNo % 2) == 0)) {
c@53: //        std::cerr << "done!" << std::endl;
c@47:         m_done = true;
c@41:     }
c@41: 
c@41:     if (m_decimator) delete[] decbuf;
c@41:     delete[] dblbuf;
c@47:     delete[] raw;
c@41: 	
c@44:     ++m_frameNo;
c@44: 
c@41:     return FeatureSet();
c@41: }
c@41: 
c@47: SimilarityPlugin::FeatureMatrix
c@47: SimilarityPlugin::calculateTimbral(FeatureSet &returnFeatures)
c@41: {
c@47:     FeatureMatrix m(m_channels); // means
c@47:     FeatureMatrix v(m_channels); // variances
c@41:     
c@41:     for (int i = 0; i < m_channels; ++i) {
c@41: 
c@42:         FeatureColumn mean(m_featureColumnSize), variance(m_featureColumnSize);
c@41: 
c@42:         for (int j = 0; j < m_featureColumnSize; ++j) {
c@41: 
c@43:             mean[j] = 0.0;
c@43:             variance[j] = 0.0;
c@41:             int count;
c@41: 
c@44:             // We want to take values up to, but not including, the
c@44:             // last non-empty frame (which may be partial)
c@43: 
c@60:             int sz = m_lastNonEmptyFrame[i] - m_emptyFrameCount[i];
c@44:             if (sz < 0) sz = 0;
c@178:             if (sz >= int(m_values[i].size())) {
c@178:                 sz = int(m_values[i].size())-1;
c@178:             }
c@43: 
c@41:             count = 0;
c@43:             for (int k = 0; k < sz; ++k) {
c@42:                 double val = m_values[i][k][j];
c@130:                 if (ISNAN(val) || ISINF(val)) continue;
c@41:                 mean[j] += val;
c@41:                 ++count;
c@41:             }
c@41:             if (count > 0) mean[j] /= count;
c@41: 
c@41:             count = 0;
c@43:             for (int k = 0; k < sz; ++k) {
c@42:                 double val = ((m_values[i][k][j] - mean[j]) *
c@42:                               (m_values[i][k][j] - mean[j]));
c@130:                 if (ISNAN(val) || ISINF(val)) continue;
c@41:                 variance[j] += val;
c@41:                 ++count;
c@41:             }
c@41:             if (count > 0) variance[j] /= count;
c@41:         }
c@41: 
c@41:         m[i] = mean;
c@41:         v[i] = variance;
c@41:     }
c@41: 
c@47:     FeatureMatrix distances(m_channels);
c@42: 
c@48:     if (m_type == TypeMFCC) {
c@48: 
c@48:         // "Despite the fact that MFCCs extracted from music are
c@48:         // clearly not Gaussian, [14] showed, somewhat surprisingly,
c@48:         // that a similarity function comparing single Gaussians
c@48:         // modelling MFCCs for each track can perform as well as
c@48:         // mixture models.  A great advantage of using single
c@48:         // Gaussians is that a simple closed form exists for the KL
c@48:         // divergence." -- Mark Levy, "Lightweight measures for
c@48:         // timbral similarity of musical audio"
c@48:         // (http://www.elec.qmul.ac.uk/easaier/papers/mlevytimbralsimilarity.pdf)
c@48: 
c@48:         KLDivergence kld;
c@48: 
c@48:         for (int i = 0; i < m_channels; ++i) {
c@48:             for (int j = 0; j < m_channels; ++j) {
c@48:                 double d = kld.distanceGaussian(m[i], v[i], m[j], v[j]);
c@48:                 distances[i].push_back(d);
c@48:             }
c@48:         }
c@48: 
c@48:     } else {
c@48: 
c@49:         // We use the KL divergence for distributions of discrete
c@49:         // variables, as chroma are histograms already.  Or at least,
c@49:         // they will be when we've normalised them like this:
c@49:         for (int i = 0; i < m_channels; ++i) {
c@49:             MathUtilities::normalise(m[i], MathUtilities::NormaliseUnitSum);
c@49:         }
c@48: 
c@48:         KLDivergence kld;
c@48: 
c@48:         for (int i = 0; i < m_channels; ++i) {
c@48:             for (int j = 0; j < m_channels; ++j) {
c@48:                 double d = kld.distanceDistribution(m[i], m[j], true);
c@48:                 distances[i].push_back(d);
c@48:             }
c@41:         }
c@41:     }
c@47:     
c@44:     Feature feature;
c@44:     feature.hasTimestamp = true;
c@44: 
c@44:     char labelBuffer[100];
c@43: 
c@41:     for (int i = 0; i < m_channels; ++i) {
c@41: 
c@41:         feature.timestamp = Vamp::RealTime(i, 0);
c@41: 
c@44:         sprintf(labelBuffer, "Means for channel %d", i+1);
c@44:         feature.label = labelBuffer;
c@44: 
c@41:         feature.values.clear();
c@42:         for (int k = 0; k < m_featureColumnSize; ++k) {
c@41:             feature.values.push_back(m[i][k]);
c@41:         }
c@41: 
c@43:         returnFeatures[m_meansOutput].push_back(feature);
c@41: 
c@44:         sprintf(labelBuffer, "Variances for channel %d", i+1);
c@44:         feature.label = labelBuffer;
c@44: 
c@41:         feature.values.clear();
c@42:         for (int k = 0; k < m_featureColumnSize; ++k) {
c@41:             feature.values.push_back(v[i][k]);
c@41:         }
c@41: 
c@43:         returnFeatures[m_variancesOutput].push_back(feature);
c@47:     }
c@47: 
c@47:     return distances;
c@47: }
c@47: 
c@47: SimilarityPlugin::FeatureMatrix
c@47: SimilarityPlugin::calculateRhythmic(FeatureSet &returnFeatures)
c@47: {
c@47:     if (!needRhythm()) return FeatureMatrix();
c@47: 
c@60: //        std::cerr << "SimilarityPlugin::initialise: rhythm clip for channel 0 contains "
c@60: //                  << m_rhythmValues[0].size() << " frames of size "
c@60: //                  << m_rhythmClipFrameSize << " at process rate "
c@60: //                  << m_processRate << " ( = "
c@60: //                  << (float(m_rhythmValues[0].size() * m_rhythmClipFrameSize) / m_processRate) << " sec )"
c@60: //                  << std::endl;
c@60: 
c@47:     BeatSpectrum bscalc;
c@47:     CosineDistance cd;
c@47: 
c@47:     // Our rhythm feature matrix is a deque of vectors for practical
c@47:     // reasons, but BeatSpectrum::process wants a vector of vectors
c@47:     // (which is what FeatureMatrix happens to be).
c@47: 
c@47:     FeatureMatrixSet bsinput(m_channels);
c@47:     for (int i = 0; i < m_channels; ++i) {
c@178:         for (int j = 0; j < int(m_rhythmValues[i].size()); ++j) {
c@47:             bsinput[i].push_back(m_rhythmValues[i][j]);
c@47:         }
c@47:     }
c@47: 
c@47:     FeatureMatrix bs(m_channels);
c@47:     for (int i = 0; i < m_channels; ++i) {
c@47:         bs[i] = bscalc.process(bsinput[i]);
c@47:     }
c@47: 
c@47:     FeatureMatrix distances(m_channels);
c@47:     for (int i = 0; i < m_channels; ++i) {
c@47:         for (int j = 0; j < m_channels; ++j) {
c@47:             double d = cd.distance(bs[i], bs[j]);
c@47:             distances[i].push_back(d);
c@47:         }
c@47:     }
c@47: 
c@47:     Feature feature;
c@47:     feature.hasTimestamp = true;
c@47: 
c@47:     char labelBuffer[100];
c@47: 
c@47:     for (int i = 0; i < m_channels; ++i) {
c@47: 
c@47:         feature.timestamp = Vamp::RealTime(i, 0);
c@47: 
c@47:         sprintf(labelBuffer, "Beat spectrum for channel %d", i+1);
c@47:         feature.label = labelBuffer;
c@47: 
c@47:         feature.values.clear();
c@178:         for (int j = 0; j < int(bs[i].size()); ++j) {
c@47:             feature.values.push_back(bs[i][j]);
c@47:         }
c@47: 
c@47:         returnFeatures[m_beatSpectraOutput].push_back(feature);
c@47:     }
c@47: 
c@47:     return distances;
c@47: }
c@47: 
c@47: double
c@47: SimilarityPlugin::getDistance(const FeatureMatrix &timbral,
c@47:                               const FeatureMatrix &rhythmic,
c@47:                               int i, int j)
c@47: {
c@47:     double distance = 1.0;
c@47:     if (needTimbre()) distance *= timbral[i][j];
c@47:     if (needRhythm()) distance *= rhythmic[i][j];
c@47:     return distance;
c@47: }
c@47: 
c@47: SimilarityPlugin::FeatureSet
c@47: SimilarityPlugin::getRemainingFeatures()
c@47: {
c@47:     FeatureSet returnFeatures;
c@47: 
c@47:     // We want to return a matrix of the distances between channels,
c@47:     // but Vamp doesn't have a matrix return type so we will actually
c@47:     // return a series of vectors
c@47: 
c@47:     FeatureMatrix timbralDistances, rhythmicDistances;
c@47: 
c@47:     if (needTimbre()) {
c@47:         timbralDistances = calculateTimbral(returnFeatures);
c@47:     }
c@47: 
c@47:     if (needRhythm()) {
c@47:         rhythmicDistances = calculateRhythmic(returnFeatures);
c@47:     }
c@47:     
c@47:     // We give all features a timestamp, otherwise hosts will tend to
c@47:     // stamp them at the end of the file, which is annoying
c@47: 
c@47:     Feature feature;
c@47:     feature.hasTimestamp = true;
c@47: 
c@47:     Feature distanceVectorFeature;
c@47:     distanceVectorFeature.label = "Distance from first channel";
c@47:     distanceVectorFeature.hasTimestamp = true;
c@47:     distanceVectorFeature.timestamp = Vamp::RealTime::zeroTime;
c@47: 
c@47:     std::map<double, int> sorted;
c@47: 
c@47:     char labelBuffer[100];
c@47: 
c@47:     for (int i = 0; i < m_channels; ++i) {
c@47: 
c@47:         feature.timestamp = Vamp::RealTime(i, 0);
c@41: 
c@41:         feature.values.clear();
c@41:         for (int j = 0; j < m_channels; ++j) {
c@47:             double dist = getDistance(timbralDistances, rhythmicDistances, i, j);
c@47:             feature.values.push_back(dist);
c@41:         }
c@43: 
c@44:         sprintf(labelBuffer, "Distances from channel %d", i+1);
c@44:         feature.label = labelBuffer;
c@41: 		
c@43:         returnFeatures[m_distanceMatrixOutput].push_back(feature);
c@43: 
c@47:         double fromFirst = 
c@47:             getDistance(timbralDistances, rhythmicDistances, 0, i);
c@44: 
c@47:         distanceVectorFeature.values.push_back(fromFirst);
c@47:         sorted[fromFirst] = i;
c@41:     }
c@41: 
c@43:     returnFeatures[m_distanceVectorOutput].push_back(distanceVectorFeature);
c@43: 
c@44:     feature.label = "Order of channels by similarity to first channel";
c@44:     feature.values.clear();
c@44:     feature.timestamp = Vamp::RealTime(0, 0);
c@44: 
c@44:     for (std::map<double, int>::iterator i = sorted.begin();
c@44:          i != sorted.end(); ++i) {
c@45:         feature.values.push_back(i->second + 1);
c@44:     }
c@44: 
c@44:     returnFeatures[m_sortedVectorOutput].push_back(feature);
c@44: 
c@44:     feature.label = "Ordered distances of channels from first channel";
c@44:     feature.values.clear();
c@44:     feature.timestamp = Vamp::RealTime(1, 0);
c@44: 
c@44:     for (std::map<double, int>::iterator i = sorted.begin();
c@44:          i != sorted.end(); ++i) {
c@44:         feature.values.push_back(i->first);
c@44:     }
c@44: 
c@44:     returnFeatures[m_sortedVectorOutput].push_back(feature);
c@44: 
c@41:     return returnFeatures;
c@41: }