Chris@0: 
Chris@0: #include "TuningDifference.h"
Chris@0: 
Chris@1: #include <iostream>
Chris@1: 
Chris@4: #include <cmath>
Chris@4: #include <cstdio>
Chris@4: 
Chris@13: #include <algorithm>
Chris@1: 
Chris@13: using namespace std;
Chris@13: 
Chris@13: static double pitchToFrequency(int pitch,
Chris@13: 			       double centsOffset = 0.,
Chris@13: 			       double concertA = 440.)
Chris@13: {
Chris@13:     double p = double(pitch) + (centsOffset / 100.);
Chris@13:     return concertA * pow(2.0, (p - 69.0) / 12.0); 
Chris@13: }
Chris@13: 
Chris@13: static double frequencyForCentsAbove440(double cents)
Chris@13: {
Chris@13:     return pitchToFrequency(69, cents, 440.);
Chris@13: }
Chris@5: 
Chris@0: TuningDifference::TuningDifference(float inputSampleRate) :
Chris@13:     Plugin(inputSampleRate),
Chris@13:     m_bpo(60),
Chris@15:     m_refChroma(new Chromagram(paramsForTuningFrequency(440.)))
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: TuningDifference::~TuningDifference()
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: TuningDifference::getIdentifier() const
Chris@0: {
Chris@1:     return "tuning-difference";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: TuningDifference::getName() const
Chris@0: {
Chris@1:     return "Tuning Difference";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: TuningDifference::getDescription() const
Chris@0: {
Chris@0:     // Return something helpful here!
Chris@0:     return "";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: TuningDifference::getMaker() const
Chris@0: {
Chris@0:     // Your name here
Chris@0:     return "";
Chris@0: }
Chris@0: 
Chris@0: int
Chris@0: TuningDifference::getPluginVersion() const
Chris@0: {
Chris@0:     // Increment this each time you release a version that behaves
Chris@0:     // differently from the previous one
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: TuningDifference::getCopyright() const
Chris@0: {
Chris@0:     // This function is not ideally named.  It does not necessarily
Chris@0:     // need to say who made the plugin -- getMaker does that -- but it
Chris@0:     // should indicate the terms under which it is distributed.  For
Chris@0:     // example, "Copyright (year). All Rights Reserved", or "GPL"
Chris@0:     return "";
Chris@0: }
Chris@0: 
Chris@0: TuningDifference::InputDomain
Chris@0: TuningDifference::getInputDomain() const
Chris@0: {
Chris@13:     return TimeDomain;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: TuningDifference::getPreferredBlockSize() const
Chris@0: {
Chris@13:     return 0;
Chris@0: }
Chris@0: 
Chris@0: size_t 
Chris@0: TuningDifference::getPreferredStepSize() const
Chris@0: {
Chris@1:     return 0;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: TuningDifference::getMinChannelCount() const
Chris@0: {
Chris@1:     return 2;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: TuningDifference::getMaxChannelCount() const
Chris@0: {
Chris@1:     return 2;
Chris@0: }
Chris@0: 
Chris@0: TuningDifference::ParameterList
Chris@0: TuningDifference::getParameterDescriptors() const
Chris@0: {
Chris@0:     ParameterList list;
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: float
Chris@1: TuningDifference::getParameter(string) const
Chris@0: {
Chris@0:     return 0;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@1: TuningDifference::setParameter(string, float) 
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: TuningDifference::ProgramList
Chris@0: TuningDifference::getPrograms() const
Chris@0: {
Chris@0:     ProgramList list;
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: TuningDifference::getCurrentProgram() const
Chris@0: {
Chris@0:     return ""; // no programs
Chris@0: }
Chris@0: 
Chris@0: void
Chris@1: TuningDifference::selectProgram(string)
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: TuningDifference::OutputList
Chris@0: TuningDifference::getOutputDescriptors() const
Chris@0: {
Chris@0:     OutputList list;
Chris@0: 
Chris@1:     OutputDescriptor d;
Chris@1:     d.identifier = "cents";
Chris@1:     d.name = "Tuning Difference";
Chris@1:     d.description = "Difference in averaged frequency profile between channels 1 and 2, in cents. A positive value means channel 2 is higher.";
Chris@1:     d.unit = "cents";
Chris@1:     d.hasFixedBinCount = true;
Chris@1:     d.binCount = 1;
Chris@1:     d.hasKnownExtents = false;
Chris@1:     d.isQuantized = false;
Chris@1:     d.sampleType = OutputDescriptor::VariableSampleRate;
Chris@1:     d.hasDuration = false;
Chris@13:     m_outputs[d.identifier] = list.size();
Chris@1:     list.push_back(d);
Chris@0: 
Chris@1:     d.identifier = "tuningfreq";
Chris@1:     d.name = "Relative Tuning Frequency";
Chris@1:     d.description = "Tuning frequency of channel 2, if channel 1 is assumed to contain the same music as it at a tuning frequency of A=440Hz.";
Chris@4:     d.unit = "hz";
Chris@1:     d.hasFixedBinCount = true;
Chris@1:     d.binCount = 1;
Chris@1:     d.hasKnownExtents = false;
Chris@1:     d.isQuantized = false;
Chris@1:     d.sampleType = OutputDescriptor::VariableSampleRate;
Chris@1:     d.hasDuration = false;
Chris@13:     m_outputs[d.identifier] = list.size();
Chris@1:     list.push_back(d);
Chris@1: 
Chris@13:     d.identifier = "reffeature";
Chris@13:     d.name = "Reference Feature";
Chris@13:     d.description = "Chroma feature from reference audio.";
Chris@0:     d.unit = "";
Chris@0:     d.hasFixedBinCount = true;
Chris@13:     d.binCount = m_bpo;
Chris@4:     d.hasKnownExtents = false;
Chris@4:     d.isQuantized = false;
Chris@4:     d.sampleType = OutputDescriptor::FixedSampleRate;
Chris@5:     d.sampleRate = 1;
Chris@4:     d.hasDuration = false;
Chris@13:     m_outputs[d.identifier] = list.size();
Chris@13:     list.push_back(d);
Chris@13: 
Chris@13:     d.identifier = "otherfeature";
Chris@13:     d.name = "Other Feature";
Chris@13:     d.description = "Chroma feature from other audio, before rotation.";
Chris@13:     d.unit = "";
Chris@13:     d.hasFixedBinCount = true;
Chris@13:     d.binCount = m_bpo;
Chris@13:     d.hasKnownExtents = false;
Chris@13:     d.isQuantized = false;
Chris@13:     d.sampleType = OutputDescriptor::FixedSampleRate;
Chris@13:     d.sampleRate = 1;
Chris@13:     d.hasDuration = false;
Chris@13:     m_outputs[d.identifier] = list.size();
Chris@13:     list.push_back(d);
Chris@13: 
Chris@13:     d.identifier = "rotfeature";
Chris@13:     d.name = "Other Feature at Rotated Frequency";
Chris@13:     d.description = "Chroma feature from reference audio calculated with the tuning frequency obtained from rotation matching.";
Chris@13:     d.unit = "";
Chris@13:     d.hasFixedBinCount = true;
Chris@13:     d.binCount = m_bpo;
Chris@13:     d.hasKnownExtents = false;
Chris@13:     d.isQuantized = false;
Chris@13:     d.sampleType = OutputDescriptor::FixedSampleRate;
Chris@13:     d.sampleRate = 1;
Chris@13:     d.hasDuration = false;
Chris@13:     m_outputs[d.identifier] = list.size();
Chris@4:     list.push_back(d);
Chris@4: 
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: bool
Chris@0: TuningDifference::initialise(size_t channels, size_t stepSize, size_t blockSize)
Chris@0: {
Chris@0:     if (channels < getMinChannelCount() ||
Chris@0: 	channels > getMaxChannelCount()) return false;
Chris@0: 
Chris@13:     if (stepSize != blockSize) return false;
Chris@0: 
Chris@1:     m_blockSize = blockSize;
Chris@1: 
Chris@1:     reset();
Chris@1:     
Chris@0:     return true;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: TuningDifference::reset()
Chris@0: {
Chris@13:     if (m_frameCount > 0) {
Chris@15: 	m_refChroma.reset(new Chromagram(paramsForTuningFrequency(440.)));
Chris@13: 	m_frameCount = 0;
Chris@13:     }
Chris@15:     m_refTotals = TFeature(m_bpo, 0.0);
Chris@13:     m_other.clear();
Chris@13: }
Chris@13: 
Chris@13: template<typename T>
Chris@13: void addTo(vector<T> &a, const vector<T> &b)
Chris@13: {
Chris@13:     transform(a.begin(), a.end(), b.begin(), a.begin(), plus<T>());
Chris@13: }
Chris@13: 
Chris@13: template<typename T>
Chris@13: T distance(const vector<T> &a, const vector<T> &b)
Chris@13: {
Chris@13:     return inner_product(a.begin(), a.end(), b.begin(), T(),
Chris@13: 			 plus<T>(), [](T x, T y) { return fabs(x - y); });
Chris@13: }
Chris@13: 
Chris@13: TuningDifference::TFeature
Chris@15: TuningDifference::computeFeatureFromTotals(const TFeature &totals) const
Chris@13: {
Chris@16:     if (m_frameCount == 0) return totals;
Chris@16:     
Chris@13:     TFeature feature(m_bpo);
Chris@13:     double sum = 0.0;
Chris@16: 
Chris@15:     for (int i = 0; i < m_bpo; ++i) {
Chris@13: 	double value = totals[i] / m_frameCount;
Chris@15: 	feature[i] += value;
Chris@13: 	sum += value;
Chris@13:     }
Chris@13: 
Chris@13:     for (int i = 0; i < m_bpo; ++i) {
Chris@13: 	feature[i] /= sum;
Chris@13:     }
Chris@13: 
Chris@13:     cerr << "computeFeatureFromTotals: feature values:" << endl;
Chris@13:     for (auto v: feature) cerr << v << " ";
Chris@13:     cerr << endl;
Chris@13:     
Chris@13:     return feature;
Chris@13: }
Chris@13: 
Chris@15: Chromagram::Parameters
Chris@13: TuningDifference::paramsForTuningFrequency(double hz) const
Chris@13: {
Chris@15:     Chromagram::Parameters params(m_inputSampleRate);
Chris@15:     params.lowestOctave = 0;
Chris@15:     params.octaves = 6;
Chris@15:     params.bpo = m_bpo;
Chris@15:     params.tuningFrequency = hz;
Chris@15:     return params;
Chris@13: }
Chris@13: 
Chris@13: TuningDifference::TFeature
Chris@13: TuningDifference::computeFeatureFromSignal(const Signal &signal, double hz) const
Chris@13: {
Chris@15:     Chromagram chromagram(paramsForTuningFrequency(hz));
Chris@13: 
Chris@15:     TFeature totals(m_bpo, 0.0);
Chris@16: 
Chris@16:     cerr << "computeFeatureFromSignal: hz = " << hz << ", frame count = " << m_frameCount << endl;
Chris@13:     
Chris@13:     for (int i = 0; i < m_frameCount; ++i) {
Chris@13: 	Signal::const_iterator first = signal.begin() + i * m_blockSize;
Chris@13: 	Signal::const_iterator last = first + m_blockSize;
Chris@13: 	if (last > signal.end()) last = signal.end();
Chris@15: 	CQBase::RealSequence input(first, last);
Chris@13: 	input.resize(m_blockSize);
Chris@15: 	CQBase::RealBlock block = chromagram.process(input);
Chris@13: 	for (const auto &v: block) addTo(totals, v);
Chris@13:     }
Chris@13: 
Chris@13:     return computeFeatureFromTotals(totals);
Chris@0: }
Chris@0: 
Chris@0: TuningDifference::FeatureSet
Chris@13: TuningDifference::process(const float *const *inputBuffers, Vamp::RealTime)
Chris@0: {
Chris@15:     CQBase::RealBlock block;
Chris@15:     CQBase::RealSequence input;
Chris@13: 
Chris@15:     input = CQBase::RealSequence
Chris@13: 	(inputBuffers[0], inputBuffers[0] + m_blockSize);
Chris@15:     block = m_refChroma->process(input);
Chris@13:     for (const auto &v: block) addTo(m_refTotals, v);
Chris@13: 
Chris@13:     m_other.insert(m_other.end(),
Chris@13: 		   inputBuffers[1], inputBuffers[1] + m_blockSize);
Chris@1:     
Chris@1:     ++m_frameCount;
Chris@0:     return FeatureSet();
Chris@0: }
Chris@0: 
Chris@13: double
Chris@13: TuningDifference::featureDistance(const TFeature &other, int rotation) const
Chris@13: {
Chris@13:     if (rotation == 0) {
Chris@13: 	return distance(m_refFeature, other);
Chris@13:     } else {
Chris@15: 	// A positive rotation pushes the tuning frequency up for this
Chris@15: 	// chroma, negative one pulls it down. If a positive rotation
Chris@15: 	// makes this chroma match an un-rotated reference, then this
Chris@15: 	// chroma must have initially been lower than the reference.
Chris@13: 	TFeature r(other);
Chris@15: 	if (rotation < 0) {
Chris@15: 	    rotate(r.begin(), r.begin() - rotation, r.end());
Chris@13: 	} else {
Chris@15: 	    rotate(r.begin(), r.end() - rotation, r.end());
Chris@13: 	}
Chris@13: 	return distance(m_refFeature, r);
Chris@13:     }
Chris@13: }
Chris@13: 
Chris@13: int
Chris@13: TuningDifference::findBestRotation(const TFeature &other) const
Chris@13: {
Chris@13:     map<double, int> dists;
Chris@13: 
Chris@13:     int maxSemis = 6;
Chris@13:     int maxRotation = (m_bpo * maxSemis) / 12;
Chris@13: 
Chris@13:     for (int r = -maxRotation; r <= maxRotation; ++r) {
Chris@13: 	double dist = featureDistance(other, r);
Chris@13: 	dists[dist] = r;
Chris@13: 	cerr << "rotation " << r << ": score " << dist << endl;
Chris@13:     }
Chris@13: 
Chris@13:     int best = dists.begin()->second;
Chris@13: 
Chris@13:     cerr << "best is " << best << endl;
Chris@13:     return best;
Chris@13: }
Chris@13: 
Chris@16: double
Chris@16: TuningDifference::findFineFrequency(int coarseCents, double coarseScore)
Chris@16: {
Chris@16:     int coarseResolution = 1200 / m_bpo;
Chris@16:     int searchDistance = coarseResolution/2 - 1;
Chris@16: 
Chris@16:     double bestScore = coarseScore;
Chris@16:     double bestHz = frequencyForCentsAbove440(coarseCents);
Chris@16: 
Chris@16:     cerr << "corresponding coarse Hz " << bestHz << " scores " << coarseScore << endl;
Chris@16:     cerr << "searchDistance = " << searchDistance << endl;
Chris@16:     
Chris@16:     for (int sign = -1; sign <= 1; sign += 2) {
Chris@16: 	for (int offset = 1; offset <= searchDistance; ++offset) {
Chris@16: 
Chris@16: 	    int fineCents = coarseCents + sign * offset;
Chris@16: 
Chris@16: 	    cerr << "trying with fineCents = " << fineCents << "..." << endl;
Chris@16: 	    
Chris@16: 	    double fineHz = frequencyForCentsAbove440(fineCents);
Chris@16: 	    TFeature fineFeature = computeFeatureFromSignal(m_other, fineHz);
Chris@16: 	    double fineScore = featureDistance(fineFeature);
Chris@16: 
Chris@16: 	    cerr << "fine offset = " << offset << ", cents = " << fineCents
Chris@16: 		 << ", Hz = " << fineHz << ", score " << fineScore
Chris@16: 		 << " (best score so far " << bestScore << ")" << endl;
Chris@16: 	    
Chris@16: 	    if (fineScore < bestScore) {
Chris@16: 		cerr << "is good!" << endl;
Chris@16: 		bestScore = fineScore;
Chris@17: 		bestHz = fineHz;
Chris@16: 	    } else {
Chris@16: 		break;
Chris@16: 	    }
Chris@16: 	}
Chris@16:     }
Chris@16: 
Chris@16:     return bestHz;
Chris@16: }
Chris@16: 
Chris@0: TuningDifference::FeatureSet
Chris@0: TuningDifference::getRemainingFeatures()
Chris@0: {
Chris@13:     FeatureSet fs;
Chris@13:     if (m_frameCount == 0) return fs;
Chris@13: 
Chris@13:     m_refFeature = computeFeatureFromTotals(m_refTotals);
Chris@13:     TFeature otherFeature = computeFeatureFromSignal(m_other, 440.);
Chris@1: 
Chris@1:     Feature f;
Chris@1: 
Chris@4:     f.values.clear();
Chris@13:     for (auto v: m_refFeature) f.values.push_back(v);
Chris@13:     fs[m_outputs["reffeature"]].push_back(f);
Chris@4: 
Chris@4:     f.values.clear();
Chris@13:     for (auto v: otherFeature) f.values.push_back(v);
Chris@13:     fs[m_outputs["otherfeature"]].push_back(f); 
Chris@13:    
Chris@13:     int rotation = findBestRotation(otherFeature);
Chris@13: 
Chris@16:     int coarseCents = -(rotation * 1200) / m_bpo;
Chris@13: 
Chris@13:     cerr << "rotation " << rotation << " -> cents " << coarseCents << endl;
Chris@13: 
Chris@13:     double coarseHz = frequencyForCentsAbove440(coarseCents);
Chris@13: 
Chris@13:     TFeature coarseFeature = computeFeatureFromSignal(m_other, coarseHz);
Chris@13:     double coarseScore = featureDistance(coarseFeature);
Chris@13: 
Chris@13:     cerr << "corresponding Hz " << coarseHz << " scores " << coarseScore << endl;
Chris@4: 
Chris@4:     f.values.clear();
Chris@13:     for (auto v: coarseFeature) f.values.push_back(v);
Chris@13:     fs[m_outputs["rotfeature"]].push_back(f);
Chris@16: 
Chris@16:     double fineHz = findFineFrequency(coarseCents, coarseScore);
Chris@16: 
Chris@16:     cerr << "overall best Hz = " << fineHz << endl;
Chris@4:     
Chris@1:     return fs;
Chris@0: }
Chris@0: