annotate ConstrainedHarmonicPeak.cpp @ 10:f82a28c2209f

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