annotate ConstrainedHarmonicPeak.cpp @ 9:5fb59edfab99

no pitch output in case of a) edge case or b) interpolated is outside freq range
author matthiasm
date Thu, 10 Apr 2014 18:37:55 +0100
parents e9b629578488
children f82a28c2209f
rev   line source
Chris@0 1
Chris@0 2 #include "ConstrainedHarmonicPeak.h"
Chris@0 3
Chris@0 4 #include <cmath>
Chris@0 5 #include <cstdio>
Chris@0 6
Chris@0 7 using std::cerr;
Chris@0 8 using std::endl;
Chris@0 9 using std::vector;
Chris@0 10
Chris@0 11 ConstrainedHarmonicPeak::ConstrainedHarmonicPeak(float inputSampleRate) :
Chris@0 12 Plugin(inputSampleRate),
Chris@4 13 m_fftSize(0),
Chris@0 14 m_minFreq(0),
Chris@6 15 m_maxFreq(22050),
Chris@1 16 m_harmonics(5)
Chris@0 17 {
Chris@0 18 }
Chris@0 19
Chris@0 20 ConstrainedHarmonicPeak::~ConstrainedHarmonicPeak()
Chris@0 21 {
Chris@0 22 }
Chris@0 23
Chris@0 24 string
Chris@0 25 ConstrainedHarmonicPeak::getIdentifier() const
Chris@0 26 {
Chris@0 27 return "constrainedharmonicpeak";
Chris@0 28 }
Chris@0 29
Chris@0 30 string
Chris@0 31 ConstrainedHarmonicPeak::getName() const
Chris@0 32 {
Chris@0 33 return "Frequency-Constrained Harmonic Peak";
Chris@0 34 }
Chris@0 35
Chris@0 36 string
Chris@0 37 ConstrainedHarmonicPeak::getDescription() const
Chris@0 38 {
Chris@6 39 return "Return the interpolated peak frequency of a harmonic product spectrum within a given frequency range";
Chris@0 40 }
Chris@0 41
Chris@0 42 string
Chris@0 43 ConstrainedHarmonicPeak::getMaker() const
Chris@0 44 {
Chris@0 45 return "Queen Mary, University of London";
Chris@0 46 }
Chris@0 47
Chris@0 48 int
Chris@0 49 ConstrainedHarmonicPeak::getPluginVersion() const
Chris@0 50 {
Chris@0 51 return 1;
Chris@0 52 }
Chris@0 53
Chris@0 54 string
Chris@0 55 ConstrainedHarmonicPeak::getCopyright() const
Chris@0 56 {
Chris@0 57 return "GPL";
Chris@0 58 }
Chris@0 59
Chris@0 60 ConstrainedHarmonicPeak::InputDomain
Chris@0 61 ConstrainedHarmonicPeak::getInputDomain() const
Chris@0 62 {
Chris@0 63 return FrequencyDomain;
Chris@0 64 }
Chris@0 65
Chris@0 66 size_t
Chris@0 67 ConstrainedHarmonicPeak::getPreferredBlockSize() const
Chris@0 68 {
Chris@0 69 return 2048;
Chris@0 70 }
Chris@0 71
Chris@0 72 size_t
Chris@0 73 ConstrainedHarmonicPeak::getPreferredStepSize() const
Chris@0 74 {
Chris@0 75 return 512;
Chris@0 76 }
Chris@0 77
Chris@0 78 size_t
Chris@0 79 ConstrainedHarmonicPeak::getMinChannelCount() const
Chris@0 80 {
Chris@0 81 return 1;
Chris@0 82 }
Chris@0 83
Chris@0 84 size_t
Chris@0 85 ConstrainedHarmonicPeak::getMaxChannelCount() const
Chris@0 86 {
Chris@0 87 return 1;
Chris@0 88 }
Chris@0 89
Chris@0 90 ConstrainedHarmonicPeak::ParameterList
Chris@0 91 ConstrainedHarmonicPeak::getParameterDescriptors() const
Chris@0 92 {
Chris@0 93 ParameterList list;
Chris@0 94
Chris@0 95 ParameterDescriptor d;
Chris@0 96 d.identifier = "minfreq";
Chris@0 97 d.name = "Minimum frequency";
Chris@6 98 d.description = "Minimum frequency for peak finding. Will be rounded down to the nearest spectral bin.";
Chris@0 99 d.unit = "Hz";
Chris@0 100 d.minValue = 0;
Chris@0 101 d.maxValue = m_inputSampleRate/2;
Chris@0 102 d.defaultValue = 0;
Chris@0 103 d.isQuantized = false;
Chris@0 104 list.push_back(d);
Chris@0 105
Chris@0 106 d.identifier = "maxfreq";
Chris@0 107 d.name = "Maximum frequency";
Chris@6 108 d.description = "Maximum frequency for peak finding. Will be rounded up to the nearest spectral bin.";
Chris@0 109 d.unit = "Hz";
Chris@0 110 d.minValue = 0;
Chris@0 111 d.maxValue = m_inputSampleRate/2;
Chris@6 112 d.defaultValue = 22050;
Chris@0 113 d.isQuantized = false;
Chris@0 114 list.push_back(d);
Chris@0 115
Chris@1 116 d.identifier = "harmonics";
Chris@1 117 d.name = "Harmonics";
Chris@1 118 d.description = "Maximum number of harmonics to consider";
Chris@1 119 d.unit = "";
Chris@1 120 d.minValue = 1;
Chris@1 121 d.maxValue = 20;
Chris@1 122 d.defaultValue = 5;
Chris@1 123 d.isQuantized = true;
Chris@1 124 d.quantizeStep = 1;
Chris@1 125 list.push_back(d);
Chris@1 126
Chris@0 127 return list;
Chris@0 128 }
Chris@0 129
Chris@0 130 float
Chris@0 131 ConstrainedHarmonicPeak::getParameter(string identifier) const
Chris@0 132 {
Chris@0 133 if (identifier == "minfreq") {
Chris@0 134 return m_minFreq;
Chris@0 135 } else if (identifier == "maxfreq") {
Chris@0 136 return m_maxFreq;
Chris@1 137 } else if (identifier == "harmonics") {
Chris@1 138 return m_harmonics;
Chris@0 139 }
Chris@0 140 return 0;
Chris@0 141 }
Chris@0 142
Chris@0 143 void
Chris@0 144 ConstrainedHarmonicPeak::setParameter(string identifier, float value)
Chris@0 145 {
Chris@0 146 if (identifier == "minfreq") {
Chris@0 147 m_minFreq = value;
Chris@0 148 } else if (identifier == "maxfreq") {
Chris@0 149 m_maxFreq = value;
Chris@1 150 } else if (identifier == "harmonics") {
Chris@1 151 m_harmonics = int(round(value));
Chris@0 152 }
Chris@0 153 }
Chris@0 154
Chris@0 155 ConstrainedHarmonicPeak::ProgramList
Chris@0 156 ConstrainedHarmonicPeak::getPrograms() const
Chris@0 157 {
Chris@0 158 ProgramList list;
Chris@0 159 return list;
Chris@0 160 }
Chris@0 161
Chris@0 162 string
Chris@0 163 ConstrainedHarmonicPeak::getCurrentProgram() const
Chris@0 164 {
Chris@0 165 return ""; // no programs
Chris@0 166 }
Chris@0 167
Chris@0 168 void
Chris@0 169 ConstrainedHarmonicPeak::selectProgram(string name)
Chris@0 170 {
Chris@0 171 }
Chris@0 172
Chris@0 173 ConstrainedHarmonicPeak::OutputList
Chris@0 174 ConstrainedHarmonicPeak::getOutputDescriptors() const
Chris@0 175 {
Chris@0 176 OutputList list;
Chris@0 177
Chris@0 178 OutputDescriptor d;
Chris@0 179 d.identifier = "peak";
Chris@0 180 d.name = "Peak frequency";
Chris@6 181 d.description = "Interpolated frequency of the harmonic spectral peak within the given frequency range";
Chris@0 182 d.unit = "Hz";
Chris@0 183 d.sampleType = OutputDescriptor::OneSamplePerStep;
Chris@0 184 d.hasDuration = false;
Chris@0 185 list.push_back(d);
Chris@0 186
Chris@0 187 return list;
Chris@0 188 }
Chris@0 189
Chris@0 190 bool
Chris@0 191 ConstrainedHarmonicPeak::initialise(size_t channels, size_t stepSize, size_t blockSize)
Chris@0 192 {
Chris@0 193 if (channels < getMinChannelCount() ||
Chris@0 194 channels > getMaxChannelCount()) {
Chris@0 195 cerr << "ConstrainedHarmonicPeak::initialise: ERROR: channels " << channels
Chris@0 196 << " out of acceptable range " << getMinChannelCount()
Chris@0 197 << " -> " << getMaxChannelCount() << endl;
Chris@0 198 return false;
Chris@0 199 }
Chris@0 200
Chris@4 201 m_fftSize = blockSize;
Chris@0 202
Chris@0 203 return true;
Chris@0 204 }
Chris@0 205
Chris@0 206 void
Chris@0 207 ConstrainedHarmonicPeak::reset()
Chris@0 208 {
Chris@0 209 }
Chris@0 210
Chris@1 211 double
Chris@1 212 ConstrainedHarmonicPeak::findInterpolatedPeak(const double *in,
Chris@1 213 int peakbin,
Chris@1 214 int bins)
Chris@1 215 {
Chris@1 216 // duplicate with SimpleCepstrum plugin
Chris@1 217 // after jos,
Chris@1 218 // https://ccrma.stanford.edu/~jos/sasp/Quadratic_Interpolation_Spectral_Peaks.html
Chris@1 219
Chris@1 220 if (peakbin < 1 || peakbin > bins - 2) {
Chris@1 221 return peakbin;
Chris@1 222 }
Chris@1 223
Chris@1 224 double alpha = in[peakbin-1];
Chris@1 225 double beta = in[peakbin];
Chris@1 226 double gamma = in[peakbin+1];
Chris@1 227
Chris@1 228 double denom = (alpha - 2*beta + gamma);
Chris@1 229
Chris@1 230 if (denom == 0) {
Chris@1 231 // flat
Chris@1 232 return peakbin;
Chris@1 233 }
Chris@1 234
Chris@1 235 double p = ((alpha - gamma) / denom) / 2.0;
Chris@1 236
Chris@1 237 return double(peakbin) + p;
Chris@1 238 }
Chris@1 239
Chris@0 240 ConstrainedHarmonicPeak::FeatureSet
Chris@0 241 ConstrainedHarmonicPeak::process(const float *const *inputBuffers, Vamp::RealTime timestamp)
Chris@0 242 {
Chris@0 243 FeatureSet fs;
Chris@0 244
Chris@6 245 // This could be better. The procedure here is
Chris@6 246 //
Chris@6 247 // 1 Produce a harmonic product spectrum within a limited
Chris@6 248 // frequency range by effectively summing the dB values of the
Chris@6 249 // bins at each multiple of the bin numbers (up to a given
Chris@6 250 // number of harmonics) in the range under consideration
Chris@6 251 // 2 Find the peak bin
Chris@6 252 // 3 Calculate the peak location by quadratic interpolation
Chris@6 253 // from the peak bin and its two neighbouring bins
Chris@6 254 //
Chris@6 255 // Problems with this:
Chris@6 256 //
Chris@6 257 // 1 Harmonics might not be located at integer multiples of the
Chris@6 258 // original bin frequency
Chris@6 259 // 2 Quadratic interpolation works "correctly" for dB-valued
Chris@6 260 // magnitude spectra but might not produce the right results in
Chris@6 261 // the dB-summed hps, especially in light of the first problem
Chris@6 262 // 3 Interpolation might not make sense at all if there are
Chris@6 263 // multiple nearby frequencies interfering across the three
Chris@6 264 // bins used for interpolation (we may be unable to identify
Chris@6 265 // the right frequency at all, but it's possible interpolation
Chris@6 266 // will make our guess worse rather than better)
Chris@6 267 //
Chris@6 268 // Possible improvements:
Chris@6 269 //
Chris@6 270 // 1 Find the higher harmonics by looking for the peak bin within
Chris@6 271 // a range around the nominal peak location
Chris@6 272 // 2 Once a peak has been identified as the peak of the HPS, use
Chris@6 273 // the original spectrum (not the HPS) to obtain the values for
Chris@6 274 // interpolation? (would help with problem 2 but might make
Chris@6 275 // problem 3 worse)
Chris@6 276
Chris@4 277 int hs = m_fftSize/2;
Chris@1 278
Chris@1 279 double *mags = new double[hs+1];
Chris@1 280 for (int i = 0; i <= hs; ++i) {
Chris@4 281 mags[i] = sqrt(inputBuffers[0][i*2] * inputBuffers[0][i*2] +
Chris@4 282 inputBuffers[0][i*2+1] * inputBuffers[0][i*2+1]);
Chris@1 283 }
Chris@1 284
Chris@1 285 // bin freq is bin * samplerate / fftsize
Chris@1 286
Chris@4 287 int minbin = int(floor((m_minFreq * m_fftSize) / m_inputSampleRate));
Chris@4 288 int maxbin = int(ceil((m_maxFreq * m_fftSize) / m_inputSampleRate));
Chris@1 289 if (minbin > hs) minbin = hs;
Chris@1 290 if (maxbin > hs) maxbin = hs;
Chris@1 291 if (maxbin <= minbin) return fs;
Chris@1 292
Chris@1 293 double *hps = new double[maxbin - minbin + 1];
Chris@1 294
Chris@1 295 // HPS in dB after MzHarmonicSpectrum
Chris@1 296
Chris@1 297 for (int bin = minbin; bin <= maxbin; ++bin) {
Chris@1 298
Chris@1 299 int i = bin - minbin;
Chris@1 300 hps[i] = 1.0;
Chris@1 301
Chris@1 302 int contributing = 0;
Chris@1 303
Chris@1 304 for (int j = 1; j <= m_harmonics; ++j) {
Chris@1 305 if (j * bin > hs) break;
Chris@1 306 hps[i] *= mags[j * bin];
Chris@1 307 ++contributing;
Chris@1 308 }
Chris@1 309
Chris@1 310 if (hps[i] <= 0.0) {
Chris@1 311 hps[i] = -120.0;
Chris@1 312 } else {
Chris@1 313 hps[i] = 20.0 / contributing * log10(hps[i]);
Chris@1 314 }
Chris@1 315 }
Chris@1 316
Chris@1 317 double maxdb = -120.0;
Chris@1 318 int maxidx = 0;
Chris@1 319 for (int i = 0; i <= maxbin - minbin; ++i) {
matthiasm@9 320 if (hps[i] > maxdb) {
matthiasm@9 321 maxdb = hps[i];
matthiasm@9 322 maxidx = i;
matthiasm@9 323 }
matthiasm@9 324 }
matthiasm@9 325
matthiasm@9 326 if (maxidx == 0 || maxidx == maxbin - minbin) { // edge cases are useless
matthiasm@9 327 return fs;
Chris@1 328 }
Chris@1 329
Chris@1 330 double interpolated = findInterpolatedPeak(hps, maxidx, maxbin - minbin + 1);
matthiasm@9 331
Chris@1 332 interpolated = interpolated + minbin;
Chris@1 333
Chris@4 334 double freq = interpolated * m_inputSampleRate / m_fftSize;
matthiasm@9 335
matthiasm@9 336 if (freq < m_minFreq || freq > m_maxFreq) {
matthiasm@9 337 return fs;
matthiasm@9 338 }
Chris@1 339
Chris@1 340 Feature f;
Chris@1 341 f.values.push_back(freq);
Chris@1 342 fs[0].push_back(f);
Chris@1 343
Chris@0 344 return fs;
Chris@0 345 }
Chris@0 346
Chris@0 347 ConstrainedHarmonicPeak::FeatureSet
Chris@0 348 ConstrainedHarmonicPeak::getRemainingFeatures()
Chris@0 349 {
Chris@0 350 FeatureSet fs;
Chris@0 351 return fs;
Chris@0 352 }
Chris@0 353