annotate data/model/FFTModel.cpp @ 1717:417528c41e66

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