annotate data/model/FFTModel.cpp @ 1111:457a1a619c5f simple-fft-model

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