Chris@31: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
Chris@31: 
Chris@31: /*
Chris@31:   Silvet
Chris@31: 
Chris@31:   A Vamp plugin for note transcription.
Chris@31:   Centre for Digital Music, Queen Mary University of London.
Chris@31:     
Chris@31:   This program is free software; you can redistribute it and/or
Chris@31:   modify it under the terms of the GNU General Public License as
Chris@31:   published by the Free Software Foundation; either version 2 of the
Chris@31:   License, or (at your option) any later version.  See the file
Chris@31:   COPYING included with this distribution for more information.
Chris@31: */
Chris@31: 
Chris@31: #include "Silvet.h"
Chris@34: #include "EM.h"
Chris@31: 
Chris@32: #include "maths/MedianFilter.h"
Chris@55: #include "maths/MathUtilities.h"
Chris@31: #include "dsp/rateconversion/Resampler.h"
Chris@31: 
Chris@32: #include "constant-q-cpp/cpp-qm-dsp/CQInterpolated.h"
Chris@31: 
Chris@31: #include <vector>
Chris@31: 
Chris@32: #include <cstdio>
Chris@32: 
Chris@31: using std::vector;
Chris@48: using std::cout;
Chris@31: using std::cerr;
Chris@31: using std::endl;
Chris@40: using Vamp::RealTime;
Chris@31: 
Chris@31: static int processingSampleRate = 44100;
Chris@31: static int processingBPO = 60;
Chris@32: static int processingHeight = 545;
Chris@38: static int processingNotes = 88;
Chris@31: 
Chris@31: Silvet::Silvet(float inputSampleRate) :
Chris@31:     Plugin(inputSampleRate),
Chris@31:     m_resampler(0),
Chris@31:     m_cq(0)
Chris@31: {
Chris@31: }
Chris@31: 
Chris@31: Silvet::~Silvet()
Chris@31: {
Chris@31:     delete m_resampler;
Chris@31:     delete m_cq;
Chris@41:     for (int i = 0; i < (int)m_postFilter.size(); ++i) {
Chris@41:         delete m_postFilter[i];
Chris@41:     }
Chris@31: }
Chris@31: 
Chris@31: string
Chris@31: Silvet::getIdentifier() const
Chris@31: {
Chris@31:     return "silvet";
Chris@31: }
Chris@31: 
Chris@31: string
Chris@31: Silvet::getName() const
Chris@31: {
Chris@31:     return "Silvet Note Transcription";
Chris@31: }
Chris@31: 
Chris@31: string
Chris@31: Silvet::getDescription() const
Chris@31: {
Chris@31:     // Return something helpful here!
Chris@31:     return "";
Chris@31: }
Chris@31: 
Chris@31: string
Chris@31: Silvet::getMaker() const
Chris@31: {
Chris@31:     // Your name here
Chris@31:     return "";
Chris@31: }
Chris@31: 
Chris@31: int
Chris@31: Silvet::getPluginVersion() const
Chris@31: {
Chris@31:     return 1;
Chris@31: }
Chris@31: 
Chris@31: string
Chris@31: Silvet::getCopyright() const
Chris@31: {
Chris@31:     // This function is not ideally named.  It does not necessarily
Chris@31:     // need to say who made the plugin -- getMaker does that -- but it
Chris@31:     // should indicate the terms under which it is distributed.  For
Chris@31:     // example, "Copyright (year). All Rights Reserved", or "GPL"
Chris@31:     return "";
Chris@31: }
Chris@31: 
Chris@31: Silvet::InputDomain
Chris@31: Silvet::getInputDomain() const
Chris@31: {
Chris@31:     return TimeDomain;
Chris@31: }
Chris@31: 
Chris@31: size_t
Chris@31: Silvet::getPreferredBlockSize() const
Chris@31: {
Chris@31:     return 0;
Chris@31: }
Chris@31: 
Chris@31: size_t 
Chris@31: Silvet::getPreferredStepSize() const
Chris@31: {
Chris@31:     return 0;
Chris@31: }
Chris@31: 
Chris@31: size_t
Chris@31: Silvet::getMinChannelCount() const
Chris@31: {
Chris@31:     return 1;
Chris@31: }
Chris@31: 
Chris@31: size_t
Chris@31: Silvet::getMaxChannelCount() const
Chris@31: {
Chris@31:     return 1;
Chris@31: }
Chris@31: 
Chris@31: Silvet::ParameterList
Chris@31: Silvet::getParameterDescriptors() const
Chris@31: {
Chris@31:     ParameterList list;
Chris@31:     return list;
Chris@31: }
Chris@31: 
Chris@31: float
Chris@31: Silvet::getParameter(string identifier) const
Chris@31: {
Chris@31:     return 0;
Chris@31: }
Chris@31: 
Chris@31: void
Chris@31: Silvet::setParameter(string identifier, float value) 
Chris@31: {
Chris@31: }
Chris@31: 
Chris@31: Silvet::ProgramList
Chris@31: Silvet::getPrograms() const
Chris@31: {
Chris@31:     ProgramList list;
Chris@31:     return list;
Chris@31: }
Chris@31: 
Chris@31: string
Chris@31: Silvet::getCurrentProgram() const
Chris@31: {
Chris@31:     return ""; 
Chris@31: }
Chris@31: 
Chris@31: void
Chris@31: Silvet::selectProgram(string name)
Chris@31: {
Chris@31: }
Chris@31: 
Chris@31: Silvet::OutputList
Chris@31: Silvet::getOutputDescriptors() const
Chris@31: {
Chris@31:     OutputList list;
Chris@31: 
Chris@31:     OutputDescriptor d;
Chris@51:     d.identifier = "notes";
Chris@51:     d.name = "Note transcription";
Chris@51:     d.description = "Overall note transcription across all instruments";
Chris@41:     d.unit = "Hz";
Chris@31:     d.hasFixedBinCount = true;
Chris@31:     d.binCount = 2;
Chris@41:     d.binNames.push_back("Frequency");
Chris@31:     d.binNames.push_back("Velocity");
Chris@31:     d.hasKnownExtents = false;
Chris@31:     d.isQuantized = false;
Chris@31:     d.sampleType = OutputDescriptor::VariableSampleRate;
Chris@51:     d.sampleRate = m_inputSampleRate / (m_cq ? m_cq->getColumnHop() : 62);
Chris@31:     d.hasDuration = true;
Chris@32:     m_notesOutputNo = list.size();
Chris@32:     list.push_back(d);
Chris@32: 
Chris@51:     d.identifier = "cq";
Chris@51:     d.name = "Raw constant-Q";
Chris@51:     d.description = "Unfiltered constant-Q time-frequency distribution";
Chris@51:     d.unit = "";
Chris@51:     d.hasFixedBinCount = true;
Chris@51:     d.binCount = processingHeight + 55;
Chris@51:     d.binNames.clear();
Chris@51:     if (m_cq) {
Chris@51:         char name[20];
Chris@51:         for (int i = 0; i < processingHeight + 55; ++i) {
Chris@51:             float freq = m_cq->getBinFrequency(i);
Chris@51:             sprintf(name, "%.1f Hz", freq);
Chris@51:             d.binNames.push_back(name);
Chris@51:         }
Chris@51:     }
Chris@51:     d.hasKnownExtents = false;
Chris@51:     d.isQuantized = false;
Chris@51:     d.sampleType = OutputDescriptor::FixedSampleRate;
Chris@51:     d.sampleRate = m_inputSampleRate / (m_cq ? m_cq->getColumnHop() : 62);
Chris@51:     d.hasDuration = false;
Chris@51:     m_cqOutputNo = list.size();
Chris@51:     list.push_back(d);
Chris@51: 
Chris@32:     d.identifier = "inputgrid";
Chris@51:     d.name = "Filtered constant-Q";
Chris@51:     d.description = "Filtered constant-Q time-frequency distribution used as input to the PLCA step";
Chris@32:     d.unit = "";
Chris@32:     d.hasFixedBinCount = true;
Chris@32:     d.binCount = processingHeight;
Chris@32:     d.binNames.clear();
Chris@32:     if (m_cq) {
Chris@32:         char name[20];
Chris@32:         for (int i = 0; i < processingHeight; ++i) {
Chris@32:             float freq = m_cq->getBinFrequency(i + 55);
Chris@32:             sprintf(name, "%.1f Hz", freq);
Chris@32:             d.binNames.push_back(name);
Chris@32:         }
Chris@32:     }
Chris@32:     d.hasKnownExtents = false;
Chris@32:     d.isQuantized = false;
Chris@32:     d.sampleType = OutputDescriptor::FixedSampleRate;
Chris@32:     d.sampleRate = 25;
Chris@32:     d.hasDuration = false;
Chris@51:     m_fcqOutputNo = list.size();
Chris@31:     list.push_back(d);
Chris@31: 
Chris@51:     d.identifier = "pitches";
Chris@51:     d.name = "Pitch activation";
Chris@51:     d.description = "Estimated pitch activation matrix";
Chris@38:     d.unit = "";
Chris@38:     d.hasFixedBinCount = true;
Chris@55:     d.binCount = processingNotes;
Chris@38:     d.binNames.clear();
Chris@55:     for (int i = 0; i < processingNotes; ++i) {
Chris@38:         d.binNames.push_back(noteName(i));
Chris@38:     }
Chris@38:     d.hasKnownExtents = false;
Chris@38:     d.isQuantized = false;
Chris@38:     d.sampleType = OutputDescriptor::FixedSampleRate;
Chris@38:     d.sampleRate = 25;
Chris@38:     d.hasDuration = false;
Chris@38:     m_pitchOutputNo = list.size();
Chris@38:     list.push_back(d);
Chris@38: 
Chris@31:     return list;
Chris@31: }
Chris@31: 
Chris@38: std::string
Chris@38: Silvet::noteName(int i) const
Chris@38: {
Chris@38:     static const char *names[] = {
Chris@38:         "A", "A#", "B", "C", "C#", "D", "D#", "E", "F", "F#", "G", "G#"
Chris@38:     };
Chris@38: 
Chris@38:     const char *n = names[i % 12];
Chris@38: 
Chris@38:     int oct = (i + 9) / 12; 
Chris@38:     
Chris@38:     char buf[20];
Chris@38:     sprintf(buf, "%s%d", n, oct);
Chris@38: 
Chris@38:     return buf;
Chris@38: }
Chris@38: 
Chris@41: float
Chris@41: Silvet::noteFrequency(int note) const
Chris@41: {
Chris@41:     return float(27.5 * pow(2.0, note / 12.0));
Chris@41: }
Chris@41: 
Chris@31: bool
Chris@31: Silvet::initialise(size_t channels, size_t stepSize, size_t blockSize)
Chris@31: {
Chris@31:     if (channels < getMinChannelCount() ||
Chris@31: 	channels > getMaxChannelCount()) return false;
Chris@31: 
Chris@31:     if (stepSize != blockSize) {
Chris@31: 	cerr << "Silvet::initialise: Step size must be the same as block size ("
Chris@31: 	     << stepSize << " != " << blockSize << ")" << endl;
Chris@31: 	return false;
Chris@31:     }
Chris@31: 
Chris@31:     m_blockSize = blockSize;
Chris@31: 
Chris@31:     reset();
Chris@31: 
Chris@31:     return true;
Chris@31: }
Chris@31: 
Chris@31: void
Chris@31: Silvet::reset()
Chris@31: {
Chris@31:     delete m_resampler;
Chris@31:     delete m_cq;
Chris@31: 
Chris@31:     if (m_inputSampleRate != processingSampleRate) {
Chris@31: 	m_resampler = new Resampler(m_inputSampleRate, processingSampleRate);
Chris@31:     } else {
Chris@31: 	m_resampler = 0;
Chris@31:     }
Chris@31: 
Chris@32:     m_cq = new CQInterpolated
Chris@32: 	(processingSampleRate, 27.5, processingSampleRate / 3, processingBPO,
Chris@32:          CQInterpolated::Linear);
Chris@31: 
Chris@41:     for (int i = 0; i < (int)m_postFilter.size(); ++i) {
Chris@41:         delete m_postFilter[i];
Chris@41:     }
Chris@41:     m_postFilter.clear();
Chris@41:     for (int i = 0; i < processingNotes; ++i) {
Chris@41:         m_postFilter.push_back(new MedianFilter<double>(3));
Chris@41:     }
Chris@41:     m_pianoRoll.clear();
Chris@32:     m_columnCount = 0;
Chris@32:     m_reducedColumnCount = 0;
Chris@40:     m_startTime = RealTime::zeroTime;
Chris@31: }
Chris@31: 
Chris@31: Silvet::FeatureSet
Chris@31: Silvet::process(const float *const *inputBuffers, Vamp::RealTime timestamp)
Chris@31: {
Chris@40:     if (m_columnCount == 0) {
Chris@40:         m_startTime = timestamp;
Chris@40:     }
Chris@40:     
Chris@31:     vector<double> data;
Chris@40:     for (int i = 0; i < m_blockSize; ++i) {
Chris@40:         data.push_back(inputBuffers[0][i]);
Chris@40:     }
Chris@31: 
Chris@31:     if (m_resampler) {
Chris@31: 	data = m_resampler->process(data.data(), data.size());
Chris@31:     }
Chris@31: 
Chris@32:     Grid cqout = m_cq->process(data);
Chris@51:     FeatureSet fs = transcribe(cqout);
Chris@51: 
Chris@51:     for (int i = 0; i < (int)cqout.size(); ++i) {
Chris@51:         Feature f;
Chris@51:         for (int j = 0; j < (int)cqout[i].size(); ++j) {
Chris@51:             f.values.push_back(float(cqout[i][j]));
Chris@51:         }
Chris@51:         fs[m_cqOutputNo].push_back(f);
Chris@51:     }
Chris@51: 
Chris@51:     return fs;
Chris@34: }
Chris@34: 
Chris@34: Silvet::FeatureSet
Chris@34: Silvet::getRemainingFeatures()
Chris@34: {
Chris@34:     Grid cqout = m_cq->getRemainingBlocks();
Chris@51:     FeatureSet fs = transcribe(cqout);
Chris@51: 
Chris@51:     for (int i = 0; i < (int)cqout.size(); ++i) {
Chris@51:         Feature f;
Chris@51:         for (int j = 0; j < (int)cqout[i].size(); ++j) {
Chris@51:             f.values.push_back(float(cqout[i][j]));
Chris@51:         }
Chris@51:         fs[m_cqOutputNo].push_back(f);
Chris@51:     }
Chris@51: 
Chris@51:     return fs;
Chris@34: }
Chris@34: 
Chris@34: Silvet::FeatureSet
Chris@34: Silvet::transcribe(const Grid &cqout)
Chris@34: {
Chris@32:     Grid filtered = preProcess(cqout);
Chris@31: 
Chris@32:     FeatureSet fs;
Chris@32: 
Chris@32:     for (int i = 0; i < (int)filtered.size(); ++i) {
Chris@32:         Feature f;
Chris@32:         for (int j = 0; j < processingHeight; ++j) {
Chris@32:             f.values.push_back(float(filtered[i][j]));
Chris@32:         }
Chris@51:         fs[m_fcqOutputNo].push_back(f);
Chris@32:     }
Chris@32: 
Chris@34:     int width = filtered.size();
Chris@34: 
Chris@34:     int iterations = 12;
Chris@34: 
Chris@34:     for (int i = 0; i < width; ++i) {
Chris@37: 
Chris@37:         double sum = 0.0;
Chris@37:         for (int j = 0; j < processingHeight; ++j) {
Chris@37:             sum += filtered[i][j];
Chris@37:         }
Chris@37: 
Chris@37:         if (sum < 1e-5) continue;
Chris@37: 
Chris@34:         EM em;
Chris@34:         for (int j = 0; j < iterations; ++j) {
Chris@34:             em.iterate(filtered[i]);
Chris@34:         }
Chris@37: 
Chris@38:         vector<double> pitches = em.getPitchDistribution();
Chris@41:         
Chris@55:         for (int j = 0; j < processingNotes; ++j) {
Chris@41:             pitches[j] *= sum;
Chris@41:         }
Chris@41: 
Chris@38:         Feature f;
Chris@55:         for (int j = 0; j < processingNotes; ++j) {
Chris@41:             f.values.push_back(float(pitches[j]));
Chris@38:         }
Chris@38:         fs[m_pitchOutputNo].push_back(f);
Chris@38: 
Chris@41:         FeatureList noteFeatures = postProcess(pitches);
Chris@45: 
Chris@41:         for (FeatureList::const_iterator fi = noteFeatures.begin();
Chris@41:              fi != noteFeatures.end(); ++fi) {
Chris@41:             fs[m_notesOutputNo].push_back(*fi);
Chris@40:         }
Chris@34:     }
Chris@34: 
Chris@32:     return fs;
Chris@31: }
Chris@31: 
Chris@32: Silvet::Grid
Chris@32: Silvet::preProcess(const Grid &in)
Chris@32: {
Chris@32:     int width = in.size();
Chris@32: 
Chris@32:     // reduce to 100 columns per second, or one column every 441 samples
Chris@32: 
Chris@32:     int spacing = processingSampleRate / 100;
Chris@32: 
Chris@32:     Grid out;
Chris@32: 
Chris@33:     //!!! nb we count the CQ latency in terms of processing hops, but
Chris@33:     //!!! actually it isn't guaranteed to be an exact number (in fact
Chris@33:     //!!! it probably isn't) so this is imprecise -- fix
Chris@33:     int latentColumns = m_cq->getLatency() / m_cq->getColumnHop();
Chris@33: 
Chris@32:     for (int i = 0; i < width; ++i) {
Chris@32: 
Chris@33:         if (m_columnCount < latentColumns) {
Chris@33:             ++m_columnCount;
Chris@33:             continue;
Chris@33:         }
Chris@33: 
Chris@32:         int prevSampleNo = (m_columnCount - 1) * m_cq->getColumnHop();
Chris@32:         int sampleNo = m_columnCount * m_cq->getColumnHop();
Chris@32: 
Chris@32:         bool select = (sampleNo / spacing != prevSampleNo / spacing);
Chris@32: 
Chris@32:         if (select) {
Chris@32:             vector<double> inCol = in[i];
Chris@32:             vector<double> outCol(processingHeight);
Chris@32: 
Chris@32:             // we reverse the column as we go (the CQ output is
Chris@32:             // "upside-down", with high frequencies at the start of
Chris@32:             // each column, and we want it the other way around) and
Chris@32:             // then ignore the first 55 (lowest-frequency) bins,
Chris@32:             // giving us 545 bins instead of 600
Chris@32: 
Chris@32:             for (int j = 0; j < processingHeight; ++j) {
Chris@46:                 int ix = inCol.size() - j - 55;
Chris@46:                 outCol[j] = inCol[ix];
Chris@46:             }
Chris@32: 
Chris@46:             vector<double> noiseLevel1 = 
Chris@46:                 MedianFilter<double>::filter(40, outCol);
Chris@46:             for (int j = 0; j < processingHeight; ++j) {
Chris@46:                 noiseLevel1[j] = std::min(outCol[j], noiseLevel1[j]);
Chris@46:             }
Chris@32: 
Chris@46:             vector<double> noiseLevel2 = 
Chris@46:                 MedianFilter<double>::filter(40, noiseLevel1);
Chris@46:             for (int j = 0; j < processingHeight; ++j) {
Chris@46:                 outCol[j] = std::max(outCol[j] - noiseLevel2[j], 0.0);
Chris@32:             }
Chris@32: 
Chris@32:             // then we only use every fourth filtered column, for 25
Chris@32:             // columns per second in the eventual grid
Chris@46:             //!!! why, if we're filtering the time columns, don't we just
Chris@46:             // reduce to this frame rate before filtering at all?
Chris@32: 
Chris@32:             if (m_reducedColumnCount % 4 == 0) {
Chris@32:                 out.push_back(outCol);
Chris@32:             }
Chris@32: 
Chris@32:             ++m_reducedColumnCount;
Chris@32:         }
Chris@32: 
Chris@32:         ++m_columnCount;
Chris@32:     }
Chris@32: 
Chris@32:     return out;
Chris@32: }
Chris@32:     
Chris@41: Vamp::Plugin::FeatureList
Chris@41: Silvet::postProcess(const vector<double> &pitches)        
Chris@41: {        
Chris@41:     vector<double> filtered;
Chris@41: 
Chris@41:     for (int j = 0; j < processingNotes; ++j) {
Chris@55:         m_postFilter[j]->push(pitches[j]);
Chris@41:         filtered.push_back(m_postFilter[j]->get());
Chris@41:     }
Chris@41: 
Chris@41:     // Threshold for level and reduce number of candidate pitches
Chris@41: 
Chris@41:     int polyphony = 5;
Chris@41:     double threshold = 4.8;
Chris@41: 
Chris@41:     typedef std::multimap<double, int> ValueIndexMap;
Chris@41: 
Chris@41:     ValueIndexMap strengths;
Chris@41:     for (int j = 0; j < processingNotes; ++j) {
Chris@41:         strengths.insert(ValueIndexMap::value_type(filtered[j], j));
Chris@41:     }
Chris@41: 
Chris@55:     map<int, double> active;
Chris@41:     ValueIndexMap::const_iterator si = strengths.end();
Chris@45:     while (int(active.size()) < polyphony) {
Chris@41:         --si;
Chris@41:         if (si->first < threshold) break;
Chris@41:         cerr << si->second << " : " << si->first << endl;
Chris@55:         active[si->second] = si->first;
Chris@45:         if (si == strengths.begin()) break;
Chris@41:     }
Chris@41: 
Chris@41:     // Minimum duration pruning, and conversion to notes. We can only
Chris@41:     // report notes that have just ended (i.e. that are absent in the
Chris@41:     // latest active set but present in the last set in the piano
Chris@41:     // roll) -- any notes that ended earlier will have been reported
Chris@41:     // already, and if they haven't ended, we don't know their
Chris@41:     // duration.
Chris@41: 
Chris@41:     int width = m_pianoRoll.size();
Chris@41: 
Chris@41:     int durationThreshold = 2; // columns
Chris@41: 
Chris@41:     FeatureList noteFeatures;
Chris@41: 
Chris@41:     if (width < durationThreshold + 1) {
Chris@41:         m_pianoRoll.push_back(active);
Chris@41:         return noteFeatures;
Chris@41:     }
Chris@41: 
Chris@41:     // we have 25 columns per second
Chris@41:     double columnDuration = 1.0 / 25.0;
Chris@41:     
Chris@55:     for (map<int, double>::const_iterator ni = m_pianoRoll[width-1].begin();
Chris@41:          ni != m_pianoRoll[width-1].end(); ++ni) {
Chris@41: 
Chris@55:         int note = ni->first;
Chris@41:         
Chris@41:         if (active.find(note) != active.end()) {
Chris@41:             // the note is still playing
Chris@41:             continue;
Chris@41:         }
Chris@41: 
Chris@41:         // the note was playing but just ended
Chris@41:         int end = width;
Chris@41:         int start = end-1;
Chris@41: 
Chris@57:         double maxStrength = 0.0;
Chris@55: 
Chris@41:         while (m_pianoRoll[start].find(note) != m_pianoRoll[start].end()) {
Chris@57:             double strength = m_pianoRoll[start][note];
Chris@57:             if (strength > maxStrength) {
Chris@57:                 maxStrength = strength;
Chris@57:             }
Chris@41:             --start;
Chris@41:         }
Chris@41:         ++start;
Chris@41: 
Chris@41:         int duration = width - start;
Chris@41:         cerr << "duration " << duration << " for just-ended note " << note << endl;
Chris@41:         if (duration < durationThreshold) {
Chris@41:             // spurious
Chris@41:             continue;
Chris@41:         }
Chris@41: 
Chris@57:         int velocity = maxStrength * 2;
Chris@55:         if (velocity > 127) velocity = 127;
Chris@55: 
Chris@41:         Feature nf;
Chris@41:         nf.hasTimestamp = true;
Chris@41:         nf.timestamp = RealTime::fromSeconds(columnDuration * start);
Chris@41:         nf.hasDuration = true;
Chris@41:         nf.duration = RealTime::fromSeconds(columnDuration * duration);
Chris@41:         nf.values.push_back(noteFrequency(note));
Chris@55:         nf.values.push_back(velocity);
Chris@41:         nf.label = noteName(note);
Chris@41:         noteFeatures.push_back(nf);
Chris@41:     }
Chris@41: 
Chris@41:     m_pianoRoll.push_back(active);
Chris@41: 
Chris@41:     cerr << "returning " << noteFeatures.size() << " complete notes" << endl;
Chris@41: 
Chris@41:     return noteFeatures;
Chris@41: }
Chris@41: