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 Chris@0: #include 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@1: return 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@0: ConstrainedHarmonicPeak::selectProgram(string name) 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@0: ConstrainedHarmonicPeak::initialise(size_t channels, size_t stepSize, 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@4: m_fftSize = 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@0: ConstrainedHarmonicPeak::process(const float *const *inputBuffers, Vamp::RealTime timestamp) 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@4: int minbin = int(floor((m_minFreq * m_fftSize) / m_inputSampleRate)); Chris@4: int maxbin = int(ceil((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@1: f.values.push_back(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: