annotate data/model/FFTModel.cpp @ 1765:2e2497cba59e by-id

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