Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add new pitch extractor PitchHPS (Harmonic Product Spectrum) #1429

Open
wants to merge 14 commits into
base: master
Choose a base branch
from
131 changes: 131 additions & 0 deletions src/algorithms/tonal/pitchhps.cpp
Original file line number Diff line number Diff line change
@@ -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;
}
}
107 changes: 107 additions & 0 deletions src/algorithms/tonal/pitchhps.h
Original file line number Diff line number Diff line change
@@ -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
114 changes: 114 additions & 0 deletions test/src/unittests/tonal/test_pitchhps.py
Original file line number Diff line number Diff line change
@@ -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)