annotate data/model/FFTModel.cpp @ 1483:7459f4c4d7c3

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