annotate data/model/FFTModel.cpp @ 1258:200c60de27ca 3.0-integration

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