c@41: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@41: 
c@41: /*
c@41:  * SegmenterPlugin.cpp
c@41:  *
c@41:  * Copyright 2008 Centre for Digital Music, Queen Mary, University of London.
c@41:  * All rights reserved.
c@41:  */
c@41: 
c@41: #include <iostream>
c@41: #include <sstream>
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@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@41: SimilarityPlugin::SimilarityPlugin(float inputSampleRate) :
c@41:     Plugin(inputSampleRate),
c@42:     m_type(TypeMFCC),
c@41:     m_mfcc(0),
c@42:     m_chromagram(0),
c@41:     m_decimator(0),
c@42:     m_featureColumnSize(20),
c@41:     m_blockSize(0),
c@41:     m_channels(0)
c@41: {
c@41: 	
c@41: }
c@41: 
c@41: SimilarityPlugin::~SimilarityPlugin()
c@41: {
c@41:     delete m_mfcc;
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@41:     return "Chris Cannam, 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@41:     return "Copyright (c) 2008 - All Rights Reserved";
c@41: }
c@41: 
c@41: size_t
c@41: SimilarityPlugin::getMinChannelCount() const
c@41: {
c@41:     return 2;
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: bool
c@41: SimilarityPlugin::initialise(size_t channels, size_t stepSize, size_t blockSize)
c@41: {
c@41:     if (channels < getMinChannelCount() ||
c@41: 	channels > getMaxChannelCount()) return false;
c@41: 
c@41:     if (stepSize != getPreferredStepSize()) {
c@41:         std::cerr << "SimilarityPlugin::initialise: supplied step size "
c@41:                   << stepSize << " differs from required step size "
c@41:                   << getPreferredStepSize() << std::endl;
c@41:         return false;
c@41:     }
c@41: 
c@41:     if (blockSize != getPreferredBlockSize()) {
c@41:         std::cerr << "SimilarityPlugin::initialise: supplied block size "
c@41:                   << blockSize << " differs from required block size "
c@41:                   << getPreferredBlockSize() << std::endl;
c@41:         return false;
c@41:     }        
c@41:     
c@41:     m_blockSize = blockSize;
c@41:     m_channels = channels;
c@41: 
c@41:     int decimationFactor = getDecimationFactor();
c@41:     if (decimationFactor > 1) {
c@42:         m_decimator = new Decimator(m_blockSize, decimationFactor);
c@41:     }
c@41: 
c@42:     if (m_type == TypeMFCC) {
c@42: 
c@42:         m_featureColumnSize = 20;
c@42: 
c@42:         MFCCConfig config;
c@42:         config.FS = lrintf(m_inputSampleRate) / decimationFactor;
c@42:         config.fftsize = 2048;
c@42:         config.nceps = m_featureColumnSize - 1;
c@42:         config.want_c0 = true;
c@42:         m_mfcc = new MFCC(config);
c@42:         m_fftSize = m_mfcc->getfftlength();
c@42: 
c@42:     } else if (m_type == TypeChroma) {
c@42: 
c@42:         m_featureColumnSize = 12;
c@42: 
c@42:         ChromaConfig config;
c@42:         config.FS = lrintf(m_inputSampleRate) / decimationFactor;
c@42:         config.min = Pitch::getFrequencyForPitch(24, 0, 440);
c@42:         config.max = Pitch::getFrequencyForPitch(96, 0, 440);
c@42:         config.BPO = 12;
c@42:         config.CQThresh = 0.0054;
c@42:         config.isNormalised = true;
c@42:         m_chromagram = new Chromagram(config);
c@42:         m_fftSize = m_chromagram->getFrameSize();
c@42: 
c@42:         std::cerr << "min = "<< config.min << ", max = " << config.max << std::endl;
c@42: 
c@42:     } else {
c@42: 
c@42:         std::cerr << "SimilarityPlugin::initialise: internal error: unknown type " << m_type << std::endl;
c@42:         return false;
c@42:     }
c@41:     
c@41:     for (int i = 0; i < m_channels; ++i) {
c@42:         m_values.push_back(FeatureMatrix());
c@41:     }
c@41: 
c@41:     return true;
c@41: }
c@41: 
c@41: void
c@41: SimilarityPlugin::reset()
c@41: {
c@41:     //!!!
c@41: }
c@41: 
c@41: int
c@41: SimilarityPlugin::getDecimationFactor() const
c@41: {
c@41:     int rate = lrintf(m_inputSampleRate);
c@41:     int internalRate = 22050;
c@41:     int decimationFactor = rate / internalRate;
c@41:     if (decimationFactor < 1) decimationFactor = 1;
c@41: 
c@41:     // must be a power of two
c@41:     while (decimationFactor & (decimationFactor - 1)) ++decimationFactor;
c@41: 
c@41:     return decimationFactor;
c@41: }
c@41: 
c@41: size_t
c@41: SimilarityPlugin::getPreferredStepSize() const
c@41: {
c@42:     if (m_blockSize == 0) calculateBlockSize();
c@42:     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@42:     if (m_type == TypeChroma) {
c@42:         ChromaConfig config;
c@42:         config.FS = lrintf(m_inputSampleRate) / decimationFactor;
c@42:         config.min = Pitch::getFrequencyForPitch(24, 0, 440);
c@42:         config.max = Pitch::getFrequencyForPitch(96, 0, 440);
c@42:         config.BPO = 12;
c@42:         config.CQThresh = 0.0054;
c@42:         config.isNormalised = false;
c@42:         Chromagram *c = new Chromagram(config);
c@42:         size_t sz = c->getFrameSize();
c@42:         delete c;
c@42:         m_blockSize = sz * decimationFactor;
c@42:     } else {
c@42:         m_blockSize = 2048 * decimationFactor;
c@42:     }
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@42:     desc.description = "";//!!!
c@42:     desc.unit = "";
c@42:     desc.minValue = 0;
c@42:     desc.maxValue = 1;
c@42:     desc.defaultValue = 0;
c@42:     desc.isQuantized = true;
c@42:     desc.quantizeStep = 1;
c@42:     desc.valueNames.push_back("Timbral (MFCC)");
c@42:     desc.valueNames.push_back("Chromatic (Chroma)");
c@42:     list.push_back(desc);	
c@42: 	
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@42:         if (m_type == TypeMFCC) return 0;
c@42:         else if (m_type == TypeChroma) return 1;
c@42:         else return 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@42:         int v = int(value + 0.1);
c@42:         Type prevType = m_type;
c@42:         if (v == 0) m_type = TypeMFCC;
c@42:         else if (v == 1) m_type = TypeChroma;
c@42:         if (m_type != prevType) m_blockSize = 0;
c@42:         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@41:     similarity.identifier = "distance";
c@41:     similarity.name = "Distance";
c@42:     similarity.description = "Distance Metric for Similarity (smaller = more similar)";
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@41:     list.push_back(similarity);
c@41: 	
c@41:     OutputDescriptor means;
c@41:     means.identifier = "means";
c@42:     means.name = "Feature Means";
c@41:     means.description = "";
c@41:     means.unit = "";
c@41:     means.hasFixedBinCount = true;
c@41:     means.binCount = m_channels;
c@41:     means.hasKnownExtents = false;
c@41:     means.isQuantized = false;
c@41:     means.sampleType = OutputDescriptor::VariableSampleRate;
c@41:     means.sampleRate = m_inputSampleRate / getPreferredStepSize();
c@41: 	
c@41:     list.push_back(means);
c@41:     
c@41:     OutputDescriptor variances;
c@41:     variances.identifier = "variances";
c@42:     variances.name = "Feature Variances";
c@41:     variances.description = "";
c@41:     variances.unit = "";
c@41:     variances.hasFixedBinCount = true;
c@41:     variances.binCount = m_channels;
c@41:     variances.hasKnownExtents = false;
c@41:     variances.isQuantized = false;
c@41:     variances.sampleType = OutputDescriptor::VariableSampleRate;
c@41:     variances.sampleRate = m_inputSampleRate / getPreferredStepSize();
c@41: 	
c@41:     list.push_back(variances);
c@41:     
c@41:     return list;
c@41: }
c@41: 
c@41: SimilarityPlugin::FeatureSet
c@41: SimilarityPlugin::process(const float *const *inputBuffers, Vamp::RealTime /* timestamp */)
c@41: {
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@42:     double *raw = 0;
c@42:     bool ownRaw = false;
c@42: 
c@42:     if (m_type == TypeMFCC) {
c@42:         raw = new double[m_featureColumnSize];
c@42:         ownRaw = true;
c@42:     }
c@41: 
c@41:     for (size_t c = 0; c < m_channels; ++c) {
c@41: 
c@41:         for (int i = 0; i < m_blockSize; ++i) {
c@41:             dblbuf[i] = inputBuffers[c][i];
c@41:         }
c@41: 
c@41:         if (m_decimator) {
c@41:             m_decimator->process(dblbuf, decbuf);
c@41:         }
c@42: 
c@42:         if (m_type == TypeMFCC) {
c@42:             m_mfcc->process(m_fftSize, decbuf, raw);
c@42:         } else if (m_type == TypeChroma) {
c@42:             raw = m_chromagram->process(decbuf);
c@42:         }                
c@41:         
c@42:         FeatureColumn mf(m_featureColumnSize);
c@42:         for (int i = 0; i < m_featureColumnSize; ++i) mf[i] = raw[i];
c@41: 
c@42:         m_values[c].push_back(mf);
c@41:     }
c@41: 
c@41:     if (m_decimator) delete[] decbuf;
c@41:     delete[] dblbuf;
c@42: 
c@42:     if (ownRaw) delete[] raw;
c@41: 	
c@41:     return FeatureSet();
c@41: }
c@41: 
c@41: SimilarityPlugin::FeatureSet
c@41: SimilarityPlugin::getRemainingFeatures()
c@41: {
c@42:     std::vector<FeatureColumn> m(m_channels);
c@42:     std::vector<FeatureColumn> v(m_channels);
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@41:             mean[j] = variance[j] = 0.0;
c@41:             int count;
c@41: 
c@42: //            std::cout << i << "," << j << ":" << std::endl;
c@42: 
c@41:             count = 0;
c@42:             for (int k = 0; k < m_values[i].size(); ++k) {
c@42:                 double val = m_values[i][k][j];
c@42: //                std::cout << val << " ";
c@41:                 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@42: 
c@42: //            std::cout << std::endl;
c@41: 
c@41:             count = 0;
c@42:             for (int k = 0; k < m_values[i].size(); ++k) {
c@42:                 double val = ((m_values[i][k][j] - mean[j]) *
c@42:                               (m_values[i][k][j] - mean[j]));
c@41:                 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@42:     // we want to return a matrix of the distances between channels,
c@41:     // but Vamp doesn't have a matrix return type so we actually
c@41:     // return a series of vectors
c@41: 
c@41:     std::vector<std::vector<double> > distances;
c@41: 
c@42:     // "Despite the fact that MFCCs extracted from music are clearly
c@42:     // not Gaussian, [14] showed, somewhat surprisingly, that a
c@42:     // similarity function comparing single Gaussians modelling MFCCs
c@42:     // for each track can perform as well as mixture models.  A great
c@42:     // advantage of using single Gaussians is that a simple closed
c@42:     // form exists for the KL divergence." -- Mark Levy, "Lightweight
c@42:     // measures for timbral similarity of musical audio"
c@42:     // (http://www.elec.qmul.ac.uk/easaier/papers/mlevytimbralsimilarity.pdf)
c@42:     //
c@42:     // This code calculates a symmetrised distance metric based on the
c@42:     // KL divergence of Gaussian models of the MFCC values.
c@42: 
c@41:     for (int i = 0; i < m_channels; ++i) {
c@41:         distances.push_back(std::vector<double>());
c@41:         for (int j = 0; j < m_channels; ++j) {
c@42:             double d = -2.0 * m_featureColumnSize;
c@42:             for (int k = 0; k < m_featureColumnSize; ++k) {
c@42:                 // m[i][k] is the mean of feature bin k for channel i
c@42:                 // v[i][k] is the variance of feature bin k for channel i
c@41:                 d += v[i][k] / v[j][k] + v[j][k] / v[i][k];
c@41:                 d += (m[i][k] - m[j][k])
c@41:                     * (1.0 / v[i][k] + 1.0 / v[j][k])
c@41:                     * (m[i][k] - m[j][k]);
c@41:             }
c@41:             d /= 2.0;
c@41:             distances[i].push_back(d);
c@41:         }
c@41:     }
c@41: 
c@41:     FeatureSet returnFeatures;
c@41: 
c@41:     for (int i = 0; i < m_channels; ++i) {
c@41: 
c@41:         Feature feature;
c@41:         feature.hasTimestamp = true; // otherwise hosts will tend to stamp them at the end of the file, which is annoying
c@41:         feature.timestamp = Vamp::RealTime(i, 0);
c@41: 
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@41:         returnFeatures[1].push_back(feature);
c@41: 
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@41:         returnFeatures[2].push_back(feature);
c@41: 
c@41:         feature.values.clear();
c@41:         for (int j = 0; j < m_channels; ++j) {
c@41:             feature.values.push_back(distances[i][j]);
c@41:         }
c@41:         ostringstream oss;
c@41:         oss << "Distance from " << (i + 1);
c@41:         feature.label = oss.str();
c@41: 		
c@41:         returnFeatures[0].push_back(feature);
c@41:     }
c@41: 
c@41:     return returnFeatures;
c@41: }