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