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@13:     m_refCQ(new CQSpectrogram(paramsForTuningFrequency(440.),
Chris@13: 			      CQSpectrogram::InterpolateHold))
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@13: 	m_refCQ.reset(new CQSpectrogram(paramsForTuningFrequency(440.),
Chris@13: 					CQSpectrogram::InterpolateHold));
Chris@13: 	m_frameCount = 0;
Chris@13:     }
Chris@13:     m_refTotals = Chroma(m_refCQ->getTotalBins(), 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@13: TuningDifference::computeFeatureFromTotals(const Chroma &totals) const
Chris@13: {
Chris@13:     TFeature feature(m_bpo);
Chris@13:     double sum = 0.0;
Chris@13:     
Chris@13:     for (int i = 0; i < (int)totals.size(); ++i) {
Chris@13: 	double value = totals[i] / m_frameCount;
Chris@13: 	feature[i % m_bpo] += 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@13: CQParameters
Chris@13: TuningDifference::paramsForTuningFrequency(double hz) const
Chris@13: {
Chris@13:     return CQParameters(m_inputSampleRate,
Chris@13: 			pitchToFrequency(36, hz),
Chris@13: 			pitchToFrequency(96, hz),
Chris@13: 			m_bpo);
Chris@13: }
Chris@13: 
Chris@13: TuningDifference::TFeature
Chris@13: TuningDifference::computeFeatureFromSignal(const Signal &signal, double hz) const
Chris@13: {
Chris@13:     CQSpectrogram cq(paramsForTuningFrequency(hz),
Chris@13: 		     CQSpectrogram::InterpolateHold);
Chris@13: 
Chris@13:     Chroma totals(m_refCQ->getTotalBins(), 0.0);
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@13: 	CQSpectrogram::RealSequence input(first, last);
Chris@13: 	input.resize(m_blockSize);
Chris@13: 	CQSpectrogram::RealBlock block = cq.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@13:     CQSpectrogram::RealBlock block;
Chris@13:     CQSpectrogram::RealSequence input;
Chris@13: 
Chris@13:     input = CQSpectrogram::RealSequence
Chris@13: 	(inputBuffers[0], inputBuffers[0] + m_blockSize);
Chris@13:     block = m_refCQ->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@13: 	TFeature r(other);
Chris@13: 	if (rotation > 0) {
Chris@13: 	    rotate(r.begin(), r.begin() + rotation, r.end());
Chris@13: 	} else {
Chris@13: 	    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@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@13:     int coarseCents = -(rotation * 100) / (m_bpo / 12);
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@4:     
Chris@1:     return fs;
Chris@0: }
Chris@0: