Add new pitch extractor PitchHPS (Harmonic Product Spectrum)

src/algorithms/tonal/pitchhps.cpp

@@ -0,0 +1,131 @@
+/*
+ * Copyright (C) 2006-2021  Music Technology Group - Universitat Pompeu Fabra
+ *
+ * This file is part of Essentia
+ *
+ * Essentia is free software: you can redistribute it and/or modify it under
+ * the terms of the GNU Affero General Public License as published by the Free
+ * Software Foundation (FSF), either version 3 of the License, or (at your
+ * option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the Affero GNU General Public License
+ * version 3 along with this program.  If not, see http://www.gnu.org/licenses/
+ */
+
+#include "pitchhps.h"
+#include "essentiamath.h"
+#include <complex>
+
+using namespace std;
+using namespace essentia;
+using namespace standard;
+
+const char* PitchHPS::name = "PitchHPS";
+const char* PitchHPS::category = "Pitch";
+const char* PitchHPS::description = DOC("This algorithm estimates the fundamental frequency given the spectrum of a monophonic music signal. It is an implementation of the Harmonic Product Spectrum algorithm [1], computed in the frequency domain. It is recommended to window the input spectrum with a Hann window. The raw spectrum can be computed with the Spectrum algorithm.\n"
+"\n"
+"An exception is thrown if an empty spectrum is provided.\n"
+"\n"
+"Note that a null “pitch” is never output by the algorithm.\n"
+"\n"
+"References:\n"
+"  [1] Noll, A. M. (1970). Pitch Determination of Human Speech by the\n"
+"  Harmonic Product Spectrum, the Harmonic Sum Spectrum, and a Maximum\n"
+"  Likelihood Estimate. Symposium on Computer Processing in Communication,\n"
+"  Ed., 19, 779–797.");
+
+void PitchHPS::configure() {
+  // compute buffer sizes
+  _frameSize = parameter("frameSize").toInt();
+  _sampleRate = parameter("sampleRate").toReal();
+  _numHarmonics = parameter("numHarmonics").toInt();
+  _magnitudeThreshold = parameter("magnitudeThreshold").toReal();
+  _minFrequency = parameter("minFrequency").toReal();
+  _maxFrequency = parameter("maxFrequency").toReal();
+
+
+  _tauMax = min(int(ceil(_sampleRate / _minFrequency)), _frameSize/2);
+  _tauMin = min(int(floor(_sampleRate / _maxFrequency)), _frameSize/2);
+
+  if (_tauMax <= _tauMin) {
+    throw EssentiaException("PitchHPS: maxFrequency is lower than minFrequency, or they are too close, or they are out of the interval of detectable frequencies with respect to the specified frameSize. Minimum detectable frequency is ", _sampleRate / (_frameSize/2), " Hz");
+  }
+
+  // configure peak detection algorithm
+  _peakDetect->configure("range", _frameSize/2+1,
+                        "minPosition", _tauMin,
+                        "maxPosition", _tauMax,
+                        "orderBy", "amplitude");
+}
+
+void PitchHPS::compute() {
+  const vector<Real>& spectrum = _spectrum.get();
+  if (spectrum.empty()) {
+    throw EssentiaException("PitchHPS: Cannot compute pitch detection on empty spectrum.");
+  }
+  Real& pitch = _pitch.get();
+//  Real& pitchConfidence = _pitchConfidence.get();
+
+  if ((int)spectrum.size() != _frameSize/2+1) {//_sqrMag.size()/2+1) {
+    Algorithm::configure( "frameSize", int(2*(spectrum.size()-1)) );
+  }
+
+  vector<Real> filteredSpectrum(spectrum);
+
+  if (_tauMin < _tauMax) {
+    double frequencyResolution = _sampleRate / spectrum.size();
+
+    size_t minBin = static_cast<size_t>(std::floor(_minFrequency / frequencyResolution));
+    size_t maxBin = static_cast<size_t>(std::floor(_maxFrequency * _numHarmonics / frequencyResolution));
+
+    std::fill(filteredSpectrum.begin(), filteredSpectrum.begin() + minBin, 0);
+    std::fill(filteredSpectrum.begin() + maxBin + 1, filteredSpectrum.end(), 0);
+
+//    Real minAmplitude = filteredSpectrum[argmax(filteredSpectrum)] * _magnitudeThreshold;
+//    int maxFreqPos = min(int(_numHarmonics * _tauMax), int(filteredSpectrum.size()));
+//
+//    for (int i = _tauMin; i < maxFreqPos; i++) {
+//      if (filteredSpectrum[i] < minAmplitude) {
+//        filteredSpectrum[i] = 0.0;
+//      }
+//    }
+  }
+
+
+  vector<Real> hps(filteredSpectrum);
+
+  for (int h=2; h < _numHarmonics; h++) {
+
+      vector<Real> downsampled(filteredSpectrum.size()/h, 1.0);
+      for (int bin=0; bin < downsampled.size(); bin++) {
+          downsampled[bin] = filteredSpectrum[bin*h];
+      }
+      for (int bin=0; bin < downsampled.size(); bin++) {
+          hps[bin] *= downsampled[bin];
+      }
+  }
+
+  for (int bin=hps.size()/_numHarmonics; bin < hps.size(); bin++) {
+      hps[bin] = 0;
+  }
+  vector<Real> _positions;
+  vector<Real> _amplitudes;
+
+  _peakDetect->input("array").set(hps);
+  _peakDetect->output("positions").set(_positions);
+  _peakDetect->output("amplitudes").set(_amplitudes);
+  _peakDetect->compute();
+
+  if (_positions.size() == 0) {
+    pitch = 0.0;
+//    pitchConfidence = 0.0;
+  } else {
+    pitch = _positions[0] * _sampleRate / _frameSize;
+//    pitchConfidence = 1.0;
+  }
+}

src/algorithms/tonal/pitchhps.h

@@ -0,0 +1,107 @@
+/*
+ * Copyright (C) 2006-2021  Music Technology Group - Universitat Pompeu Fabra
+ *
+ * This file is part of Essentia
+ *
+ * Essentia is free software: you can redistribute it and/or modify it under
+ * the terms of the GNU Affero General Public License as published by the Free
+ * Software Foundation (FSF), either version 3 of the License, or (at your
+ * option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful, but WITHOUT
+ * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+ * FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+ * details.
+ *
+ * You should have received a copy of the Affero GNU General Public License
+ * version 3 along with this program.  If not, see http://www.gnu.org/licenses/
+ */
+
+#ifndef ESSENTIA_PITCHHPS_H
+#define ESSENTIA_PITCHHPS_H
+
+#include "algorithmfactory.h"
+
+namespace essentia {
+namespace standard {
+
+class PitchHPS : public Algorithm {
+
+ private:
+  Input<std::vector<Real> > _spectrum;
+  Output<Real> _pitch;
+//  Output<Real> _pitchConfidence;
+
+  Algorithm* _peakDetect;
+
+  std::vector<Real> _positions;   /** peak positions */
+  std::vector<Real> _amplitudes;  /** peak amplitudes */
+  Real _sampleRate;               /** sampling rate of the audio signal */
+  int _frameSize;
+  int _tauMin;
+  int _tauMax;
+  int _numHarmonics;
+  int _minFrequency;
+  int _maxFrequency;
+  Real _magnitudeThreshold;
+
+ public:
+  PitchHPS() {
+    declareInput(_spectrum, "spectrum", "the input spectrum (preferably created with a hann window)");
+    declareOutput(_pitch, "pitch", "detected pitch [Hz]");
+//    declareOutput(_pitchConfidence, "pitchConfidence", "confidence with which the pitch was detected [0,1]");
+
+    _peakDetect = AlgorithmFactory::create("PeakDetection");
+  }
+
+  ~PitchHPS() {
+    delete _peakDetect;
+  };
+
+  void declareParameters() {
+    declareParameter("frameSize", "number of samples in the input spectrum", "[2,inf)", 2048);
+    declareParameter("sampleRate", "sampling rate of the input spectrum [Hz]", "(0,inf)", 44100.);
+    declareParameter("minFrequency", "the minimum allowed frequency [Hz]", "(0,inf)", 20.0);
+    declareParameter("maxFrequency", "the maximum allowed frequency [Hz]", "(0,inf)", 22050.0);
+    declareParameter("numHarmonics", "number of harmonics to consider", "[1,inf)", 5);
+    declareParameter("magnitudeThreshold", "threshold ratio for the amplitude filtering [0,1]", "[0,1]", 0.2);
+  }
+
+  void configure();
+  void compute();
+
+  static const char* name;
+  static const char* category;
+  static const char* description;
+
+}; // class PitchHPS
+
+} // namespace standard
+} // namespace essentia
+
+
+#include "streamingalgorithmwrapper.h"
+
+namespace essentia {
+namespace streaming {
+
+class PitchHPS : public StreamingAlgorithmWrapper {
+
+ protected:
+  Sink<std::vector<Real> > _spectrum;
+  Source<Real> _pitch;
+//  Source<Real> _pitchConfidence;
+
+ public:
+  PitchHPS() {
+    declareAlgorithm("PitchHPS");
+    declareInput(_spectrum, TOKEN, "spectrum");
+    declareOutput(_pitch, TOKEN, "pitch");
+//    declareOutput(_pitchConfidence, TOKEN, "pitchConfidence");
+  }
+};
+
+} // namespace streaming
+} // namespace essentia
+
+#endif // ESSENTIA_PITCHHPS_H

test/src/unittests/tonal/test_pitchhps.py

@@ -0,0 +1,114 @@
+#!/usr/bin/env python
+
+# Copyright (C) 2006-2021  Music Technology Group - Universitat Pompeu Fabra
+#
+# This file is part of Essentia
+#
+# Essentia is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the Free
+# Software Foundation (FSF), either version 3 of the License, or (at your
+# option) any later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the Affero GNU General Public License
+# version 3 along with this program. If not, see http://www.gnu.org/licenses/
+
+
+
+from essentia_test import *
+from numpy import sin, pi, mean, random
+
+class TestPitchHPS(TestCase):
+
+    def testEmpty(self):
+        self.assertComputeFails(PitchHPS(), [])
+
+    def testZero(self):
+        pitch = PitchHPS()(zeros(1024))
+        self.assertEqual(pitch, 0)
+
+
+    def testSine(self):
+        sr = 44100
+        size = sr*1;
+        freq = 440
+        signal = [sin(2.0*pi*freq*i/sr) for i in range(size)]
+        self.runTest(signal, sr, freq)
+
+    def testBandLimitedSquare(self):
+        sr = 44100
+        size = sr*1;
+        freq = 660
+        w = 2.0*pi*freq
+        nharms = 10
+        signal = zeros(size)
+        for i in range(size):
+            for harm in range(nharms):
+                signal[i] += .5/(2.*harm+1)*sin((2*harm+1)*i*w/sr)
+
+        self.runTest(signal, sr, freq)
+
+    def testBandLimitedSaw(self):
+        sr = 44100
+        size = sr*1;
+        freq = 660
+        w = 2.0*pi*freq
+        nharms = 10
+        signal = zeros(size)
+        for i in range(1,size):
+            for harm in range(1,nharms+1):
+                signal[i] += 1./harm*sin(harm*i*w/sr)
+        self.runTest(signal, sr, freq, 1.1)
+
+    def testBandLimitedSawMasked(self):
+        sr = 44100
+        size = sr*1;
+        freq = 440
+        w = 2.0*pi*freq
+        subw = 2.0*pi*(freq-100)
+        nharms = 10
+        signal = zeros(size)
+        for i in range(1,size):
+            # masking noise:
+            whitenoise = 2*(random.rand(1)-0.5)
+            signal[i] += 2*whitenoise
+            for harm in range(1,nharms):
+                signal[i] += 1./harm*sin(i*harm*w/sr)
+        signal = 5*LowPass()(signal)
+        for i in range(1,size):
+            for harm in range(1,nharms+1):
+                signal[i] += .1/harm*sin(i*harm*w/sr)
+            signal[i] += 0.5*sin(i*subw/sr)
+        max_signal = max(signal) + 1
+        signal = signal/max_signal
+        self.runTest(signal, sr, freq, 1.5)
+
+
+    def runTest(self, signal, sr, freq, pitch_precision = 1):
+        frameSize = 1024
+        hopsize = frameSize
+
+        frames = FrameGenerator(signal, frameSize=frameSize, hopSize=hopsize)
+        win = Windowing(type='hann')
+        pitchDetect = PitchHPS(frameSize=frameSize, sampleRate = sr)
+        pitch = []
+        for frame in frames:
+            spec = Spectrum()(win(frame))
+            f = pitchDetect(spec)
+            pitch += [f]
+
+        self.assertAlmostEqual(mean(f), freq, pitch_precision)
+
+    def testInvalidParam(self):
+        self.assertConfigureFails(PitchHPS(), {'frameSize' : 1})
+        self.assertConfigureFails(PitchHPS(), {'sampleRate' : 0})
+
+
+suite = allTests(TestPitchHPS)
+
+if __name__ == '__main__':
+    TextTestRunner(verbosity=2).run(suite)