Chris@7: /*
Chris@7:     CHP
Chris@7:     Copyright (c) 2014 Queen Mary, University of London
Chris@7: 
Chris@7:     Permission is hereby granted, free of charge, to any person
Chris@7:     obtaining a copy of this software and associated documentation
Chris@7:     files (the "Software"), to deal in the Software without
Chris@7:     restriction, including without limitation the rights to use, copy,
Chris@7:     modify, merge, publish, distribute, sublicense, and/or sell copies
Chris@7:     of the Software, and to permit persons to whom the Software is
Chris@7:     furnished to do so, subject to the following conditions:
Chris@7: 
Chris@7:     The above copyright notice and this permission notice shall be
Chris@7:     included in all copies or substantial portions of the Software.
Chris@7: 
Chris@7:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
Chris@7:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
Chris@7:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
Chris@7:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
Chris@7:     CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
Chris@7:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
Chris@7:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
Chris@7: */
Chris@0: 
Chris@0: #include "ConstrainedHarmonicPeak.h"
Chris@0: 
Chris@0: #include <cmath>
Chris@0: #include <cstdio>
Chris@14: #include <climits>
Chris@0: 
Chris@0: using std::cerr;
Chris@0: using std::endl;
Chris@0: using std::vector;
Chris@0: 
Chris@0: ConstrainedHarmonicPeak::ConstrainedHarmonicPeak(float inputSampleRate) :
Chris@0:     Plugin(inputSampleRate),
Chris@4:     m_fftSize(0),
Chris@0:     m_minFreq(0),
Chris@6:     m_maxFreq(22050),
Chris@1:     m_harmonics(5)
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: ConstrainedHarmonicPeak::~ConstrainedHarmonicPeak()
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: ConstrainedHarmonicPeak::getIdentifier() const
Chris@0: {
Chris@0:     return "constrainedharmonicpeak";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: ConstrainedHarmonicPeak::getName() const
Chris@0: {
Chris@0:     return "Frequency-Constrained Harmonic Peak";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: ConstrainedHarmonicPeak::getDescription() const
Chris@0: {
Chris@6:     return "Return the interpolated peak frequency of a harmonic product spectrum within a given frequency range";
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: ConstrainedHarmonicPeak::getMaker() const
Chris@0: {
Chris@0:     return "Queen Mary, University of London";
Chris@0: }
Chris@0: 
Chris@0: int
Chris@0: ConstrainedHarmonicPeak::getPluginVersion() const
Chris@0: {
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: ConstrainedHarmonicPeak::getCopyright() const
Chris@0: {
Chris@0:     return "GPL";
Chris@0: }
Chris@0: 
Chris@0: ConstrainedHarmonicPeak::InputDomain
Chris@0: ConstrainedHarmonicPeak::getInputDomain() const
Chris@0: {
Chris@0:     return FrequencyDomain;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: ConstrainedHarmonicPeak::getPreferredBlockSize() const
Chris@0: {
Chris@0:     return 2048;
Chris@0: }
Chris@0: 
Chris@0: size_t 
Chris@0: ConstrainedHarmonicPeak::getPreferredStepSize() const
Chris@0: {
Chris@0:     return 512;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: ConstrainedHarmonicPeak::getMinChannelCount() const
Chris@0: {
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: size_t
Chris@0: ConstrainedHarmonicPeak::getMaxChannelCount() const
Chris@0: {
Chris@0:     return 1;
Chris@0: }
Chris@0: 
Chris@0: ConstrainedHarmonicPeak::ParameterList
Chris@0: ConstrainedHarmonicPeak::getParameterDescriptors() const
Chris@0: {
Chris@0:     ParameterList list;
Chris@0: 
Chris@0:     ParameterDescriptor d;
Chris@0:     d.identifier = "minfreq";
Chris@0:     d.name = "Minimum frequency";
Chris@6:     d.description = "Minimum frequency for peak finding. Will be rounded down to the nearest spectral bin.";
Chris@0:     d.unit = "Hz";
Chris@0:     d.minValue = 0;
Chris@0:     d.maxValue = m_inputSampleRate/2;
Chris@0:     d.defaultValue = 0;
Chris@0:     d.isQuantized = false;
Chris@0:     list.push_back(d);
Chris@0: 
Chris@0:     d.identifier = "maxfreq";
Chris@0:     d.name = "Maximum frequency";
Chris@6:     d.description = "Maximum frequency for peak finding. Will be rounded up to the nearest spectral bin.";
Chris@0:     d.unit = "Hz";
Chris@0:     d.minValue = 0;
Chris@0:     d.maxValue = m_inputSampleRate/2;
Chris@6:     d.defaultValue = 22050;
Chris@0:     d.isQuantized = false;
Chris@0:     list.push_back(d);
Chris@0: 
Chris@1:     d.identifier = "harmonics";
Chris@1:     d.name = "Harmonics";
Chris@1:     d.description = "Maximum number of harmonics to consider";
Chris@1:     d.unit = "";
Chris@1:     d.minValue = 1;
Chris@1:     d.maxValue = 20;
Chris@1:     d.defaultValue = 5;
Chris@1:     d.isQuantized = true;
Chris@1:     d.quantizeStep = 1;
Chris@1:     list.push_back(d);
Chris@1: 
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: float
Chris@0: ConstrainedHarmonicPeak::getParameter(string identifier) const
Chris@0: {
Chris@0:     if (identifier == "minfreq") {
Chris@0: 	return m_minFreq;
Chris@0:     } else if (identifier == "maxfreq") {
Chris@0: 	return m_maxFreq;
Chris@1:     } else if (identifier == "harmonics") {
Chris@14: 	return float(m_harmonics);
Chris@0:     }
Chris@0:     return 0;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: ConstrainedHarmonicPeak::setParameter(string identifier, float value) 
Chris@0: {
Chris@0:     if (identifier == "minfreq") {
Chris@0: 	m_minFreq = value;
Chris@0:     } else if (identifier == "maxfreq") {
Chris@0: 	m_maxFreq = value;
Chris@1:     } else if (identifier == "harmonics") {
Chris@1: 	m_harmonics = int(round(value));
Chris@0:     }
Chris@0: }
Chris@0: 
Chris@0: ConstrainedHarmonicPeak::ProgramList
Chris@0: ConstrainedHarmonicPeak::getPrograms() const
Chris@0: {
Chris@0:     ProgramList list;
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: string
Chris@0: ConstrainedHarmonicPeak::getCurrentProgram() const
Chris@0: {
Chris@0:     return ""; // no programs
Chris@0: }
Chris@0: 
Chris@0: void
Chris@14: ConstrainedHarmonicPeak::selectProgram(string)
Chris@0: {
Chris@0: }
Chris@0: 
Chris@0: ConstrainedHarmonicPeak::OutputList
Chris@0: ConstrainedHarmonicPeak::getOutputDescriptors() const
Chris@0: {
Chris@0:     OutputList list;
Chris@0: 
Chris@0:     OutputDescriptor d;
Chris@0:     d.identifier = "peak";
Chris@0:     d.name = "Peak frequency";
Chris@6:     d.description = "Interpolated frequency of the harmonic spectral peak within the given frequency range";
Chris@0:     d.unit = "Hz";
Chris@0:     d.sampleType = OutputDescriptor::OneSamplePerStep;
Chris@0:     d.hasDuration = false;
Chris@0:     list.push_back(d);
Chris@0: 
Chris@0:     return list;
Chris@0: }
Chris@0: 
Chris@0: bool
Chris@14: ConstrainedHarmonicPeak::initialise(size_t channels, size_t, size_t blockSize)
Chris@0: {
Chris@0:     if (channels < getMinChannelCount() ||
Chris@0: 	channels > getMaxChannelCount()) {
Chris@0: 	cerr << "ConstrainedHarmonicPeak::initialise: ERROR: channels " << channels
Chris@0: 	     << " out of acceptable range " << getMinChannelCount()
Chris@0: 	     << " -> " << getMaxChannelCount() << endl;
Chris@0: 	return false;
Chris@0:     }
Chris@0: 
Chris@14:     if (blockSize > INT_MAX) {
Chris@14:         return false; // arf
Chris@14:     }
Chris@14:     
Chris@14:     m_fftSize = int(blockSize);
Chris@0: 
Chris@0:     return true;
Chris@0: }
Chris@0: 
Chris@0: void
Chris@0: ConstrainedHarmonicPeak::reset()
Chris@0: {
Chris@0: }
Chris@0: 
Chris@1: double
Chris@1: ConstrainedHarmonicPeak::findInterpolatedPeak(const double *in, 
Chris@1: 					      int peakbin,
Chris@1: 					      int bins)
Chris@1: {
Chris@1:     // duplicate with SimpleCepstrum plugin
Chris@1:     // after jos, 
Chris@1:     // https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
Chris@1: 
Chris@1:     if (peakbin < 1 || peakbin > bins - 2) {
Chris@1:         return peakbin;
Chris@1:     }
Chris@1: 
Chris@1:     double alpha = in[peakbin-1];
Chris@1:     double beta  = in[peakbin];
Chris@1:     double gamma = in[peakbin+1];
Chris@1: 
Chris@1:     double denom = (alpha - 2*beta + gamma);
Chris@1: 
Chris@1:     if (denom == 0) {
Chris@1:         // flat
Chris@1:         return peakbin;
Chris@1:     }
Chris@1: 
Chris@1:     double p = ((alpha - gamma) / denom) / 2.0;
Chris@1: 
Chris@1:     return double(peakbin) + p;
Chris@1: }
Chris@1: 
Chris@0: ConstrainedHarmonicPeak::FeatureSet
Chris@14: ConstrainedHarmonicPeak::process(const float *const *inputBuffers, Vamp::RealTime)
Chris@0: {
Chris@0:     FeatureSet fs;
Chris@0: 
Chris@6:     // This could be better. The procedure here is
Chris@6:     // 
Chris@6:     // 1 Produce a harmonic product spectrum within a limited
Chris@6:     //   frequency range by effectively summing the dB values of the
Chris@6:     //   bins at each multiple of the bin numbers (up to a given
Chris@6:     //   number of harmonics) in the range under consideration
Chris@6:     // 2 Find the peak bin
Chris@6:     // 3 Calculate the peak location by quadratic interpolation
Chris@6:     //   from the peak bin and its two neighbouring bins
Chris@6:     //
Chris@6:     // Problems with this: 
Chris@6:     //
Chris@6:     // 1 Harmonics might not be located at integer multiples of the
Chris@6:     //   original bin frequency
Chris@6:     // 2 Quadratic interpolation works "correctly" for dB-valued
Chris@6:     //   magnitude spectra but might not produce the right results in
Chris@6:     //   the dB-summed hps, especially in light of the first problem
Chris@6:     // 3 Interpolation might not make sense at all if there are
Chris@6:     //   multiple nearby frequencies interfering across the three
Chris@6:     //   bins used for interpolation (we may be unable to identify
Chris@6:     //   the right frequency at all, but it's possible interpolation
Chris@6:     //   will make our guess worse rather than better)
Chris@6:     //
Chris@6:     // Possible improvements:
Chris@6:     // 
Chris@6:     // 1 Find the higher harmonics by looking for the peak bin within
Chris@6:     //   a range around the nominal peak location
Chris@6:     // 2 Once a peak has been identified as the peak of the HPS, use
Chris@6:     //   the original spectrum (not the HPS) to obtain the values for
Chris@6:     //   interpolation? (would help with problem 2 but might make
Chris@6:     //   problem 3 worse)
Chris@6: 
Chris@4:     int hs = m_fftSize/2;
Chris@1: 
Chris@1:     double *mags = new double[hs+1];
Chris@1:     for (int i = 0; i <= hs; ++i) {
Chris@4: 	mags[i] = sqrt(inputBuffers[0][i*2] * inputBuffers[0][i*2] +
Chris@4: 		       inputBuffers[0][i*2+1] * inputBuffers[0][i*2+1]);
Chris@1:     }
Chris@1: 
Chris@1:     // bin freq is bin * samplerate / fftsize
Chris@1: 
Chris@14:     int minbin = int(floor((double(m_minFreq) * m_fftSize) / m_inputSampleRate));
Chris@14:     int maxbin = int(ceil((double(m_maxFreq) * m_fftSize) / m_inputSampleRate));
Chris@1:     if (minbin > hs) minbin = hs;
Chris@1:     if (maxbin > hs) maxbin = hs;
Chris@1:     if (maxbin <= minbin) return fs;
Chris@1: 
Chris@1:     double *hps = new double[maxbin - minbin + 1];
Chris@1: 
Chris@1:     // HPS in dB after MzHarmonicSpectrum
Chris@1: 
Chris@1:     for (int bin = minbin; bin <= maxbin; ++bin) {
Chris@1: 
Chris@1: 	int i = bin - minbin;
Chris@1: 	hps[i] = 1.0;
Chris@1: 
Chris@1: 	int contributing = 0;
Chris@1: 
Chris@1: 	for (int j = 1; j <= m_harmonics; ++j) {
Chris@1: 	    if (j * bin > hs) break;
Chris@1: 	    hps[i] *= mags[j * bin];
Chris@1: 	    ++contributing;
Chris@1: 	}
Chris@1: 
Chris@1: 	if (hps[i] <= 0.0) {
Chris@1: 	    hps[i] = -120.0;
Chris@1: 	} else {
Chris@1: 	    hps[i] = 20.0 / contributing * log10(hps[i]);
Chris@1: 	}
Chris@1:     }
Chris@1: 
Chris@1:     double maxdb = -120.0;
Chris@1:     int maxidx = 0;
Chris@1:     for (int i = 0; i <= maxbin - minbin; ++i) {
matthiasm@9:         if (hps[i] > maxdb) {
matthiasm@9:             maxdb = hps[i];
matthiasm@9:             maxidx = i;
matthiasm@9:         }
matthiasm@9:     }
matthiasm@9: 
matthiasm@9:     if (maxidx == 0 || maxidx == maxbin - minbin) { // edge cases are useless
matthiasm@9:         return fs;
Chris@1:     }
Chris@1: 
Chris@1:     double interpolated = findInterpolatedPeak(hps, maxidx, maxbin - minbin + 1);
matthiasm@9:     
Chris@1:     interpolated = interpolated + minbin;
Chris@1: 
Chris@4:     double freq = interpolated * m_inputSampleRate / m_fftSize;
matthiasm@9:     
matthiasm@9:     if (freq < m_minFreq || freq > m_maxFreq) {
matthiasm@9:         return fs;
matthiasm@9:     }
Chris@1: 
Chris@1:     Feature f;
Chris@14:     f.values.push_back(float(freq));
Chris@1:     fs[0].push_back(f);
Chris@1: 
Chris@0:     return fs;
Chris@0: }
Chris@0: 
Chris@0: ConstrainedHarmonicPeak::FeatureSet
Chris@0: ConstrainedHarmonicPeak::getRemainingFeatures()
Chris@0: {
Chris@0:     FeatureSet fs;
Chris@0:     return fs;
Chris@0: }
Chris@0: