annotate data/model/FFTModel.cpp @ 1858:14747f24ad04

Use cancellable serialiser; add some more profiling points
author Chris Cannam
date Thu, 14 May 2020 16:36:48 +0100
parents 915d316a5609
children
rev   line source
Chris@152 1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
Chris@152 2
Chris@152 3 /*
Chris@152 4 Sonic Visualiser
Chris@152 5 An audio file viewer and annotation editor.
Chris@152 6 Centre for Digital Music, Queen Mary, University of London.
Chris@152 7 This file copyright 2006 Chris Cannam.
Chris@152 8
Chris@152 9 This program is free software; you can redistribute it and/or
Chris@152 10 modify it under the terms of the GNU General Public License as
Chris@152 11 published by the Free Software Foundation; either version 2 of the
Chris@152 12 License, or (at your option) any later version. See the file
Chris@152 13 COPYING included with this distribution for more information.
Chris@152 14 */
Chris@152 15
Chris@152 16 #include "FFTModel.h"
Chris@152 17 #include "DenseTimeValueModel.h"
Chris@152 18
Chris@183 19 #include "base/Profiler.h"
Chris@275 20 #include "base/Pitch.h"
Chris@1256 21 #include "base/HitCount.h"
Chris@1428 22 #include "base/Debug.h"
Chris@1573 23 #include "base/MovingMedian.h"
Chris@183 24
Chris@402 25 #include <algorithm>
Chris@402 26
Chris@152 27 #include <cassert>
Chris@1090 28 #include <deque>
Chris@152 29
Chris@1090 30 using namespace std;
Chris@1090 31
Chris@1256 32 static HitCount inSmallCache("FFTModel: Small FFT cache");
Chris@1256 33 static HitCount inSourceCache("FFTModel: Source data cache");
Chris@1256 34
Chris@1744 35 FFTModel::FFTModel(ModelId modelId,
Chris@152 36 int channel,
Chris@152 37 WindowType windowType,
Chris@929 38 int windowSize,
Chris@929 39 int windowIncrement,
Chris@1090 40 int fftSize) :
Chris@1744 41 m_model(modelId),
Chris@1780 42 m_sampleRate(0),
Chris@1090 43 m_channel(channel),
Chris@1090 44 m_windowType(windowType),
Chris@1090 45 m_windowSize(windowSize),
Chris@1090 46 m_windowIncrement(windowIncrement),
Chris@1090 47 m_fftSize(fftSize),
Chris@1091 48 m_windower(windowType, windowSize),
Chris@1093 49 m_fft(fftSize),
Chris@1780 50 m_maximumFrequency(0.0),
Chris@1371 51 m_cacheWriteIndex(0),
Chris@1093 52 m_cacheSize(3)
Chris@152 53 {
Chris@1837 54 clearCaches();
Chris@1371 55
Chris@1091 56 if (m_windowSize > m_fftSize) {
Chris@1428 57 SVCERR << "ERROR: FFTModel::FFTModel: window size (" << m_windowSize
Chris@1680 58 << ") may not exceed FFT size (" << m_fftSize << ")" << endl;
Chris@1680 59 throw invalid_argument("FFTModel window size may not exceed FFT size");
Chris@1091 60 }
Chris@1133 61
Chris@1270 62 m_fft.initFloat();
Chris@1270 63
Chris@1744 64 auto model = ModelById::getAs<DenseTimeValueModel>(m_model);
Chris@1744 65 if (model) {
Chris@1780 66 m_sampleRate = model->getSampleRate();
Chris@1780 67
Chris@1752 68 connect(model.get(), SIGNAL(modelChanged(ModelId)),
Chris@1752 69 this, SIGNAL(modelChanged(ModelId)));
Chris@1752 70 connect(model.get(), SIGNAL(modelChangedWithin(ModelId, sv_frame_t, sv_frame_t)),
Chris@1752 71 this, SIGNAL(modelChangedWithin(ModelId, sv_frame_t, sv_frame_t)));
Chris@1784 72 } else {
Chris@1784 73 m_error = QString("Model #%1 is not available").arg(m_model.untyped);
Chris@1744 74 }
Chris@152 75 }
Chris@152 76
Chris@152 77 FFTModel::~FFTModel()
Chris@152 78 {
Chris@152 79 }
Chris@152 80
Chris@1837 81 void
Chris@1837 82 FFTModel::clearCaches()
Chris@1837 83 {
Chris@1837 84 m_cached.clear();
Chris@1837 85 while (m_cached.size() < m_cacheSize) {
Chris@1837 86 m_cached.push_back({ -1, complexvec_t(m_fftSize / 2 + 1) });
Chris@1837 87 }
Chris@1837 88 m_cacheWriteIndex = 0;
Chris@1837 89 m_savedData.range = { 0, 0 };
Chris@1837 90 }
Chris@1837 91
Chris@1744 92 bool
Chris@1744 93 FFTModel::isOK() const
Chris@360 94 {
Chris@1744 95 auto model = ModelById::getAs<DenseTimeValueModel>(m_model);
Chris@1784 96 if (!model) {
Chris@1784 97 m_error = QString("Model #%1 is not available").arg(m_model.untyped);
Chris@1784 98 return false;
Chris@1784 99 }
Chris@1784 100 if (!model->isOK()) {
Chris@1784 101 m_error = QString("Model #%1 is not OK").arg(m_model.untyped);
Chris@1784 102 return false;
Chris@1784 103 }
Chris@1784 104 return true;
Chris@1744 105 }
Chris@1744 106
Chris@1744 107 int
Chris@1744 108 FFTModel::getCompletion() const
Chris@1744 109 {
Chris@1744 110 int c = 100;
Chris@1744 111 auto model = ModelById::getAs<DenseTimeValueModel>(m_model);
Chris@1744 112 if (model) {
Chris@1744 113 if (model->isReady(&c)) return 100;
Chris@360 114 }
Chris@1744 115 return c;
Chris@1744 116 }
Chris@1744 117
Chris@1780 118 void
Chris@1780 119 FFTModel::setMaximumFrequency(double freq)
Chris@1780 120 {
Chris@1780 121 m_maximumFrequency = freq;
Chris@1837 122 clearCaches();
Chris@360 123 }
Chris@360 124
Chris@1091 125 int
Chris@1091 126 FFTModel::getWidth() const
Chris@1091 127 {
Chris@1744 128 auto model = ModelById::getAs<DenseTimeValueModel>(m_model);
Chris@1744 129 if (!model) return 0;
Chris@1744 130 return int((model->getEndFrame() - model->getStartFrame())
Chris@1091 131 / m_windowIncrement) + 1;
Chris@1091 132 }
Chris@1091 133
Chris@1091 134 int
Chris@1091 135 FFTModel::getHeight() const
Chris@1091 136 {
Chris@1780 137 int height = m_fftSize / 2 + 1;
Chris@1780 138 if (m_maximumFrequency != 0.0) {
Chris@1780 139 int maxBin = int(ceil(m_maximumFrequency * m_fftSize) / m_sampleRate);
Chris@1780 140 if (maxBin >= 0 && maxBin < height) {
Chris@1780 141 return maxBin + 1;
Chris@1780 142 }
Chris@1780 143 }
Chris@1780 144 return height;
Chris@1091 145 }
Chris@1091 146
Chris@152 147 QString
Chris@929 148 FFTModel::getBinName(int n) const
Chris@152 149 {
Chris@1790 150 return tr("%1 Hz").arg(getBinValue(n));
Chris@1790 151 }
Chris@1790 152
Chris@1790 153 float
Chris@1790 154 FFTModel::getBinValue(int n) const
Chris@1790 155 {
Chris@1790 156 return float((m_sampleRate * n) / m_fftSize);
Chris@152 157 }
Chris@152 158
Chris@1091 159 FFTModel::Column
Chris@1091 160 FFTModel::getColumn(int x) const
Chris@1091 161 {
Chris@1091 162 auto cplx = getFFTColumn(x);
Chris@1091 163 Column col;
Chris@1154 164 col.reserve(cplx.size());
Chris@1091 165 for (auto c: cplx) col.push_back(abs(c));
Chris@1319 166 return col;
Chris@1091 167 }
Chris@1091 168
Chris@1200 169 FFTModel::Column
Chris@1200 170 FFTModel::getPhases(int x) const
Chris@1200 171 {
Chris@1200 172 auto cplx = getFFTColumn(x);
Chris@1200 173 Column col;
Chris@1200 174 col.reserve(cplx.size());
Chris@1201 175 for (auto c: cplx) {
Chris@1201 176 col.push_back(arg(c));
Chris@1201 177 }
Chris@1319 178 return col;
Chris@1200 179 }
Chris@1200 180
Chris@1091 181 float
Chris@1091 182 FFTModel::getMagnitudeAt(int x, int y) const
Chris@1091 183 {
Chris@1569 184 if (x < 0 || x >= getWidth() || y < 0 || y >= getHeight()) {
Chris@1569 185 return 0.f;
Chris@1569 186 }
Chris@1093 187 auto col = getFFTColumn(x);
Chris@1093 188 return abs(col[y]);
Chris@1091 189 }
Chris@1091 190
Chris@1091 191 float
Chris@1091 192 FFTModel::getMaximumMagnitudeAt(int x) const
Chris@1091 193 {
Chris@1091 194 Column col(getColumn(x));
Chris@1092 195 float max = 0.f;
Chris@1154 196 int n = int(col.size());
Chris@1154 197 for (int i = 0; i < n; ++i) {
Chris@1092 198 if (col[i] > max) max = col[i];
Chris@1092 199 }
Chris@1092 200 return max;
Chris@1091 201 }
Chris@1091 202
Chris@1091 203 float
Chris@1091 204 FFTModel::getPhaseAt(int x, int y) const
Chris@1091 205 {
Chris@1093 206 if (x < 0 || x >= getWidth() || y < 0 || y >= getHeight()) return 0.f;
Chris@1091 207 return arg(getFFTColumn(x)[y]);
Chris@1091 208 }
Chris@1091 209
Chris@1091 210 void
Chris@1091 211 FFTModel::getValuesAt(int x, int y, float &re, float &im) const
Chris@1091 212 {
Chris@1780 213 if (x < 0 || x >= getWidth() || y < 0 || y >= getHeight()) {
Chris@1780 214 re = 0.f;
Chris@1780 215 im = 0.f;
Chris@1780 216 return;
Chris@1780 217 }
Chris@1091 218 auto col = getFFTColumn(x);
Chris@1091 219 re = col[y].real();
Chris@1091 220 im = col[y].imag();
Chris@1091 221 }
Chris@1091 222
Chris@1091 223 bool
Chris@1091 224 FFTModel::getMagnitudesAt(int x, float *values, int minbin, int count) const
Chris@1091 225 {
Chris@1839 226 if (count == 0) {
Chris@1839 227 count = getHeight() - minbin;
Chris@1839 228 }
Chris@1091 229 auto col = getFFTColumn(x);
Chris@1091 230 for (int i = 0; i < count; ++i) {
Chris@1091 231 values[i] = abs(col[minbin + i]);
Chris@1091 232 }
Chris@1091 233 return true;
Chris@1091 234 }
Chris@1091 235
Chris@1091 236 bool
Chris@1091 237 FFTModel::getPhasesAt(int x, float *values, int minbin, int count) const
Chris@1091 238 {
Chris@1091 239 if (count == 0) count = getHeight();
Chris@1091 240 auto col = getFFTColumn(x);
Chris@1091 241 for (int i = 0; i < count; ++i) {
Chris@1091 242 values[i] = arg(col[minbin + i]);
Chris@1091 243 }
Chris@1091 244 return true;
Chris@1091 245 }
Chris@1091 246
Chris@1091 247 bool
Chris@1091 248 FFTModel::getValuesAt(int x, float *reals, float *imags, int minbin, int count) const
Chris@1091 249 {
Chris@1091 250 if (count == 0) count = getHeight();
Chris@1091 251 auto col = getFFTColumn(x);
Chris@1091 252 for (int i = 0; i < count; ++i) {
Chris@1091 253 reals[i] = col[minbin + i].real();
Chris@1091 254 }
Chris@1091 255 for (int i = 0; i < count; ++i) {
Chris@1091 256 imags[i] = col[minbin + i].imag();
Chris@1091 257 }
Chris@1091 258 return true;
Chris@1091 259 }
Chris@1091 260
Chris@1837 261 floatvec_t
Chris@1091 262 FFTModel::getSourceSamples(int column) const
Chris@1091 263 {
Chris@1094 264 // m_fftSize may be greater than m_windowSize, but not the reverse
Chris@1094 265
Chris@1094 266 // cerr << "getSourceSamples(" << column << ")" << endl;
Chris@1094 267
Chris@1091 268 auto range = getSourceSampleRange(column);
Chris@1094 269 auto data = getSourceData(range);
Chris@1094 270
Chris@1091 271 int off = (m_fftSize - m_windowSize) / 2;
Chris@1094 272
Chris@1094 273 if (off == 0) {
Chris@1094 274 return data;
Chris@1094 275 } else {
Chris@1094 276 vector<float> pad(off, 0.f);
Chris@1837 277 floatvec_t padded;
Chris@1094 278 padded.reserve(m_fftSize);
Chris@1094 279 padded.insert(padded.end(), pad.begin(), pad.end());
Chris@1094 280 padded.insert(padded.end(), data.begin(), data.end());
Chris@1094 281 padded.insert(padded.end(), pad.begin(), pad.end());
Chris@1094 282 return padded;
Chris@1094 283 }
Chris@1094 284 }
Chris@1094 285
Chris@1837 286 floatvec_t
Chris@1094 287 FFTModel::getSourceData(pair<sv_frame_t, sv_frame_t> range) const
Chris@1094 288 {
Chris@1094 289 // cerr << "getSourceData(" << range.first << "," << range.second
Chris@1094 290 // << "): saved range is (" << m_savedData.range.first
Chris@1094 291 // << "," << m_savedData.range.second << ")" << endl;
Chris@1094 292
Chris@1100 293 if (m_savedData.range == range) {
Chris@1256 294 inSourceCache.hit();
Chris@1100 295 return m_savedData.data;
Chris@1100 296 }
Chris@1094 297
Chris@1270 298 Profiler profiler("FFTModel::getSourceData (cache miss)");
Chris@1270 299
Chris@1094 300 if (range.first < m_savedData.range.second &&
Chris@1094 301 range.first >= m_savedData.range.first &&
Chris@1094 302 range.second > m_savedData.range.second) {
Chris@1094 303
Chris@1256 304 inSourceCache.partial();
Chris@1256 305
Chris@1100 306 sv_frame_t discard = range.first - m_savedData.range.first;
Chris@1100 307
Chris@1837 308 floatvec_t data;
Chris@1457 309 data.reserve(range.second - range.first);
Chris@1094 310
Chris@1457 311 data.insert(data.end(),
Chris@1457 312 m_savedData.data.begin() + discard,
Chris@1457 313 m_savedData.data.end());
Chris@1100 314
Chris@1837 315 floatvec_t rest = getSourceDataUncached
Chris@1457 316 ({ m_savedData.range.second, range.second });
Chris@1457 317
Chris@1457 318 data.insert(data.end(), rest.begin(), rest.end());
Chris@1094 319
Chris@1457 320 m_savedData = { range, data };
Chris@1457 321 return data;
Chris@1095 322
Chris@1095 323 } else {
Chris@1095 324
Chris@1256 325 inSourceCache.miss();
Chris@1256 326
Chris@1095 327 auto data = getSourceDataUncached(range);
Chris@1095 328 m_savedData = { range, data };
Chris@1095 329 return data;
Chris@1094 330 }
Chris@1095 331 }
Chris@1094 332
Chris@1837 333 floatvec_t
Chris@1095 334 FFTModel::getSourceDataUncached(pair<sv_frame_t, sv_frame_t> range) const
Chris@1095 335 {
Chris@1457 336 Profiler profiler("FFTModel::getSourceDataUncached");
Chris@1688 337
Chris@1744 338 auto model = ModelById::getAs<DenseTimeValueModel>(m_model);
Chris@1744 339 if (!model) return {};
Chris@1457 340
Chris@1091 341 decltype(range.first) pfx = 0;
Chris@1091 342 if (range.first < 0) {
Chris@1091 343 pfx = -range.first;
Chris@1091 344 range = { 0, range.second };
Chris@1091 345 }
Chris@1096 346
Chris@1744 347 auto data = model->getData(m_channel,
Chris@1744 348 range.first,
Chris@1744 349 range.second - range.first);
Chris@1773 350 /*
Chris@1281 351 if (data.empty()) {
Chris@1281 352 SVDEBUG << "NOTE: empty source data for range (" << range.first << ","
Chris@1281 353 << range.second << ") (model end frame "
Chris@1744 354 << model->getEndFrame() << ")" << endl;
Chris@1281 355 }
Chris@1773 356 */
Chris@1281 357
Chris@1096 358 // don't return a partial frame
Chris@1096 359 data.resize(range.second - range.first, 0.f);
Chris@1096 360
Chris@1096 361 if (pfx > 0) {
Chris@1096 362 vector<float> pad(pfx, 0.f);
Chris@1096 363 data.insert(data.begin(), pad.begin(), pad.end());
Chris@1096 364 }
Chris@1096 365
Chris@1091 366 if (m_channel == -1) {
Chris@1744 367 int channels = model->getChannelCount();
Chris@1429 368 if (channels > 1) {
Chris@1096 369 int n = int(data.size());
Chris@1096 370 float factor = 1.f / float(channels);
Chris@1100 371 // use mean instead of sum for fft model input
Chris@1429 372 for (int i = 0; i < n; ++i) {
Chris@1429 373 data[i] *= factor;
Chris@1429 374 }
Chris@1429 375 }
Chris@1091 376 }
Chris@1094 377
Chris@1094 378 return data;
Chris@1091 379 }
Chris@1091 380
Chris@1837 381 const complexvec_t &
Chris@1093 382 FFTModel::getFFTColumn(int n) const
Chris@1091 383 {
Chris@1258 384 // The small cache (i.e. the m_cached deque) is for cases where
Chris@1258 385 // values are looked up individually, and for e.g. peak-frequency
Chris@1258 386 // spectrograms where values from two consecutive columns are
Chris@1257 387 // needed at once. This cache gets essentially no hits when
Chris@1258 388 // scrolling through a magnitude spectrogram, but 95%+ hits with a
Chris@1569 389 // peak-frequency spectrogram or spectrum.
Chris@1257 390 for (const auto &incache : m_cached) {
Chris@1093 391 if (incache.n == n) {
Chris@1256 392 inSmallCache.hit();
Chris@1837 393 return incache.col;
Chris@1093 394 }
Chris@1093 395 }
Chris@1256 396 inSmallCache.miss();
Chris@1258 397
Chris@1258 398 Profiler profiler("FFTModel::getFFTColumn (cache miss)");
Chris@1093 399
Chris@1093 400 auto samples = getSourceSamples(n);
Chris@1567 401 m_windower.cut(samples.data() + (m_fftSize - m_windowSize) / 2);
Chris@1270 402 breakfastquay::v_fftshift(samples.data(), m_fftSize);
Chris@1270 403
Chris@1837 404 complexvec_t &col = m_cached[m_cacheWriteIndex].col;
Chris@1837 405
Chris@1837 406 // expand to large enough for fft destination, if truncated previously
Chris@1837 407 col.resize(m_fftSize / 2 + 1);
Chris@1270 408
Chris@1270 409 m_fft.forwardInterleaved(samples.data(),
Chris@1270 410 reinterpret_cast<float *>(col.data()));
Chris@1093 411
Chris@1837 412 // keep only the number of elements we need - so that we can
Chris@1837 413 // return a const ref without having to resize on a cache hit
Chris@1837 414 col.resize(getHeight());
Chris@1837 415
Chris@1371 416 m_cached[m_cacheWriteIndex].n = n;
Chris@1371 417
Chris@1371 418 m_cacheWriteIndex = (m_cacheWriteIndex + 1) % m_cacheSize;
Chris@1093 419
Chris@1837 420 return col;
Chris@1091 421 }
Chris@1091 422
Chris@275 423 bool
Chris@1045 424 FFTModel::estimateStableFrequency(int x, int y, double &frequency)
Chris@275 425 {
Chris@275 426 if (!isOK()) return false;
Chris@275 427
Chris@1090 428 frequency = double(y * getSampleRate()) / m_fftSize;
Chris@275 429
Chris@275 430 if (x+1 >= getWidth()) return false;
Chris@275 431
Chris@275 432 // At frequency f, a phase shift of 2pi (one cycle) happens in 1/f sec.
Chris@275 433 // At hopsize h and sample rate sr, one hop happens in h/sr sec.
Chris@275 434 // At window size w, for bin b, f is b*sr/w.
Chris@275 435 // thus 2pi phase shift happens in w/(b*sr) sec.
Chris@275 436 // We need to know what phase shift we expect from h/sr sec.
Chris@275 437 // -> 2pi * ((h/sr) / (w/(b*sr)))
Chris@275 438 // = 2pi * ((h * b * sr) / (w * sr))
Chris@275 439 // = 2pi * (h * b) / w.
Chris@275 440
Chris@1038 441 double oldPhase = getPhaseAt(x, y);
Chris@1038 442 double newPhase = getPhaseAt(x+1, y);
Chris@275 443
Chris@929 444 int incr = getResolution();
Chris@275 445
Chris@1090 446 double expectedPhase = oldPhase + (2.0 * M_PI * y * incr) / m_fftSize;
Chris@275 447
Chris@1038 448 double phaseError = princarg(newPhase - expectedPhase);
Chris@275 449
Chris@275 450 // The new frequency estimate based on the phase error resulting
Chris@275 451 // from assuming the "native" frequency of this bin
Chris@275 452
Chris@275 453 frequency =
Chris@1090 454 (getSampleRate() * (expectedPhase + phaseError - oldPhase)) /
Chris@1045 455 (2.0 * M_PI * incr);
Chris@275 456
Chris@275 457 return true;
Chris@275 458 }
Chris@275 459
Chris@275 460 FFTModel::PeakLocationSet
Chris@1191 461 FFTModel::getPeaks(PeakPickType type, int x, int ymin, int ymax) const
Chris@275 462 {
Chris@551 463 Profiler profiler("FFTModel::getPeaks");
Chris@1575 464
Chris@275 465 FFTModel::PeakLocationSet peaks;
Chris@275 466 if (!isOK()) return peaks;
Chris@275 467
Chris@275 468 if (ymax == 0 || ymax > getHeight() - 1) {
Chris@275 469 ymax = getHeight() - 1;
Chris@275 470 }
Chris@275 471
Chris@275 472 if (type == AllPeaks) {
Chris@551 473 int minbin = ymin;
Chris@551 474 if (minbin > 0) minbin = minbin - 1;
Chris@551 475 int maxbin = ymax;
Chris@551 476 if (maxbin < getHeight() - 1) maxbin = maxbin + 1;
Chris@551 477 const int n = maxbin - minbin + 1;
Chris@1218 478 float *values = new float[n];
Chris@551 479 getMagnitudesAt(x, values, minbin, maxbin - minbin + 1);
Chris@929 480 for (int bin = ymin; bin <= ymax; ++bin) {
Chris@551 481 if (bin == minbin || bin == maxbin) continue;
Chris@551 482 if (values[bin - minbin] > values[bin - minbin - 1] &&
Chris@551 483 values[bin - minbin] > values[bin - minbin + 1]) {
Chris@275 484 peaks.insert(bin);
Chris@275 485 }
Chris@275 486 }
Chris@1218 487 delete[] values;
Chris@275 488 return peaks;
Chris@275 489 }
Chris@275 490
Chris@551 491 Column values = getColumn(x);
Chris@1154 492 int nv = int(values.size());
Chris@275 493
Chris@500 494 float mean = 0.f;
Chris@1154 495 for (int i = 0; i < nv; ++i) mean += values[i];
Chris@1154 496 if (nv > 0) mean = mean / float(values.size());
Chris@1038 497
Chris@275 498 // For peak picking we use a moving median window, picking the
Chris@275 499 // highest value within each continuous region of values that
Chris@275 500 // exceed the median. For pitch adaptivity, we adjust the window
Chris@275 501 // size to a roughly constant pitch range (about four tones).
Chris@275 502
Chris@1040 503 sv_samplerate_t sampleRate = getSampleRate();
Chris@275 504
Chris@1090 505 vector<int> inrange;
Chris@1576 506 double dist = 0.5;
Chris@500 507
Chris@929 508 int medianWinSize = getPeakPickWindowSize(type, sampleRate, ymin, dist);
Chris@929 509 int halfWin = medianWinSize/2;
Chris@275 510
Chris@1573 511 MovingMedian<float> window(medianWinSize);
Chris@1573 512
Chris@929 513 int binmin;
Chris@275 514 if (ymin > halfWin) binmin = ymin - halfWin;
Chris@275 515 else binmin = 0;
Chris@275 516
Chris@929 517 int binmax;
Chris@1154 518 if (ymax + halfWin < nv) binmax = ymax + halfWin;
Chris@1154 519 else binmax = nv - 1;
Chris@275 520
Chris@929 521 int prevcentre = 0;
Chris@500 522
Chris@929 523 for (int bin = binmin; bin <= binmax; ++bin) {
Chris@275 524
Chris@275 525 float value = values[bin];
Chris@275 526
Chris@280 527 // so-called median will actually be the dist*100'th percentile
Chris@280 528 medianWinSize = getPeakPickWindowSize(type, sampleRate, bin, dist);
Chris@275 529 halfWin = medianWinSize/2;
Chris@275 530
Chris@1573 531 int actualSize = std::min(medianWinSize, bin - binmin + 1);
Chris@1573 532 window.resize(actualSize);
Chris@1573 533 window.setPercentile(dist * 100.0);
Chris@1573 534 window.push(value);
Chris@275 535
Chris@275 536 if (type == MajorPitchAdaptivePeaks) {
Chris@1154 537 if (ymax + halfWin < nv) binmax = ymax + halfWin;
Chris@1154 538 else binmax = nv - 1;
Chris@275 539 }
Chris@275 540
Chris@1573 541 float median = window.get();
Chris@275 542
Chris@929 543 int centrebin = 0;
Chris@500 544 if (bin > actualSize/2) centrebin = bin - actualSize/2;
Chris@500 545
Chris@500 546 while (centrebin > prevcentre || bin == binmin) {
Chris@275 547
Chris@500 548 if (centrebin > prevcentre) ++prevcentre;
Chris@500 549
Chris@500 550 float centre = values[prevcentre];
Chris@500 551
Chris@500 552 if (centre > median) {
Chris@500 553 inrange.push_back(centrebin);
Chris@500 554 }
Chris@500 555
Chris@1154 556 if (centre <= median || centrebin+1 == nv) {
Chris@500 557 if (!inrange.empty()) {
Chris@929 558 int peakbin = 0;
Chris@500 559 float peakval = 0.f;
Chris@929 560 for (int i = 0; i < (int)inrange.size(); ++i) {
Chris@500 561 if (i == 0 || values[inrange[i]] > peakval) {
Chris@500 562 peakval = values[inrange[i]];
Chris@500 563 peakbin = inrange[i];
Chris@500 564 }
Chris@500 565 }
Chris@500 566 inrange.clear();
Chris@500 567 if (peakbin >= ymin && peakbin <= ymax) {
Chris@500 568 peaks.insert(peakbin);
Chris@275 569 }
Chris@275 570 }
Chris@275 571 }
Chris@500 572
Chris@500 573 if (bin == binmin) break;
Chris@275 574 }
Chris@275 575 }
Chris@275 576
Chris@275 577 return peaks;
Chris@275 578 }
Chris@275 579
Chris@929 580 int
Chris@1040 581 FFTModel::getPeakPickWindowSize(PeakPickType type, sv_samplerate_t sampleRate,
Chris@1576 582 int bin, double &dist) const
Chris@275 583 {
Chris@1576 584 dist = 0.5; // dist is percentile / 100.0
Chris@275 585 if (type == MajorPeaks) return 10;
Chris@275 586 if (bin == 0) return 3;
Chris@280 587
Chris@1091 588 double binfreq = (sampleRate * bin) / m_fftSize;
Chris@1038 589 double hifreq = Pitch::getFrequencyForPitch(73, 0, binfreq);
Chris@280 590
Chris@1091 591 int hibin = int(lrint((hifreq * m_fftSize) / sampleRate));
Chris@275 592 int medianWinSize = hibin - bin;
Chris@1576 593
Chris@1575 594 if (medianWinSize < 3) {
Chris@1575 595 medianWinSize = 3;
Chris@1575 596 }
Chris@1576 597
Chris@1576 598 // We want to avoid the median window size changing too often, as
Chris@1576 599 // it requires a reallocation. So snap to a nearby round number.
Chris@1576 600
Chris@1575 601 if (medianWinSize > 20) {
Chris@1575 602 medianWinSize = (1 + medianWinSize / 10) * 10;
Chris@1575 603 }
Chris@1576 604 if (medianWinSize > 200) {
Chris@1576 605 medianWinSize = (1 + medianWinSize / 100) * 100;
Chris@1576 606 }
Chris@1576 607 if (medianWinSize > 2000) {
Chris@1576 608 medianWinSize = (1 + medianWinSize / 1000) * 1000;
Chris@1576 609 }
Chris@1576 610 if (medianWinSize > 20000) {
Chris@1576 611 medianWinSize = 20000;
Chris@1575 612 }
Chris@280 613
Chris@1576 614 if (medianWinSize < 100) {
Chris@1576 615 dist = 1.0 - (4.0 / medianWinSize);
Chris@1576 616 } else {
Chris@1576 617 dist = 1.0 - (8.0 / medianWinSize);
Chris@1576 618 }
Chris@1576 619 if (dist < 0.5) dist = 0.5;
Chris@1575 620
Chris@275 621 return medianWinSize;
Chris@275 622 }
Chris@275 623
Chris@275 624 FFTModel::PeakSet
Chris@929 625 FFTModel::getPeakFrequencies(PeakPickType type, int x,
Chris@1191 626 int ymin, int ymax) const
Chris@275 627 {
Chris@551 628 Profiler profiler("FFTModel::getPeakFrequencies");
Chris@551 629
Chris@275 630 PeakSet peaks;
Chris@275 631 if (!isOK()) return peaks;
Chris@275 632 PeakLocationSet locations = getPeaks(type, x, ymin, ymax);
Chris@275 633
Chris@1040 634 sv_samplerate_t sampleRate = getSampleRate();
Chris@929 635 int incr = getResolution();
Chris@275 636
Chris@275 637 // This duplicates some of the work of estimateStableFrequency to
Chris@275 638 // allow us to retrieve the phases in two separate vertical
Chris@275 639 // columns, instead of jumping back and forth between columns x and
Chris@275 640 // x+1, which may be significantly slower if re-seeking is needed
Chris@275 641
Chris@1090 642 vector<float> phases;
Chris@275 643 for (PeakLocationSet::iterator i = locations.begin();
Chris@275 644 i != locations.end(); ++i) {
Chris@275 645 phases.push_back(getPhaseAt(x, *i));
Chris@275 646 }
Chris@275 647
Chris@929 648 int phaseIndex = 0;
Chris@275 649 for (PeakLocationSet::iterator i = locations.begin();
Chris@275 650 i != locations.end(); ++i) {
Chris@1038 651 double oldPhase = phases[phaseIndex];
Chris@1038 652 double newPhase = getPhaseAt(x+1, *i);
Chris@1090 653 double expectedPhase = oldPhase + (2.0 * M_PI * *i * incr) / m_fftSize;
Chris@1038 654 double phaseError = princarg(newPhase - expectedPhase);
Chris@1038 655 double frequency =
Chris@275 656 (sampleRate * (expectedPhase + phaseError - oldPhase))
Chris@275 657 / (2 * M_PI * incr);
Chris@1045 658 peaks[*i] = frequency;
Chris@275 659 ++phaseIndex;
Chris@275 660 }
Chris@275 661
Chris@275 662 return peaks;
Chris@275 663 }
Chris@275 664