annotate data/model/FFTModel.cpp @ 1281:6463046bd21d 3.0-integration

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