comparison data/model/FFTModel.cpp @ 1090:420fc961c0c4 simple-fft-model

Gut the old code, but don't replace it yet (so nothing will link yet)
author Chris Cannam
date Fri, 12 Jun 2015 14:51:46 +0100
parents 0fd3661bcfff
children bdebff3265ae
comparison
equal deleted inserted replaced
1089:655cd4e68e9a 1090:420fc961c0c4
13 COPYING included with this distribution for more information. 13 COPYING included with this distribution for more information.
14 */ 14 */
15 15
16 #include "FFTModel.h" 16 #include "FFTModel.h"
17 #include "DenseTimeValueModel.h" 17 #include "DenseTimeValueModel.h"
18 #include "AggregateWaveModel.h"
19 18
20 #include "base/Profiler.h" 19 #include "base/Profiler.h"
21 #include "base/Pitch.h" 20 #include "base/Pitch.h"
22 21
23 #include <algorithm> 22 #include <algorithm>
24 23
25 #include <cassert> 24 #include <cassert>
25 #include <deque>
26 26
27 #ifndef __GNUC__ 27 #ifndef __GNUC__
28 #include <alloca.h> 28 #include <alloca.h>
29 #endif 29 #endif
30
31 using namespace std;
30 32
31 FFTModel::FFTModel(const DenseTimeValueModel *model, 33 FFTModel::FFTModel(const DenseTimeValueModel *model,
32 int channel, 34 int channel,
33 WindowType windowType, 35 WindowType windowType,
34 int windowSize, 36 int windowSize,
35 int windowIncrement, 37 int windowIncrement,
36 int fftSize, 38 int fftSize) :
37 bool polar, 39 m_model(model),
38 StorageAdviser::Criteria criteria, 40 m_channel(channel),
39 sv_frame_t fillFromFrame) : 41 m_windowType(windowType),
40 //!!! ZoomConstraint! 42 m_windowSize(windowSize),
41 m_server(0), 43 m_windowIncrement(windowIncrement),
42 m_xshift(0), 44 m_fftSize(fftSize),
43 m_yshift(0) 45 m_windower(windowType, windowSize)
44 { 46 {
45 setSourceModel(const_cast<DenseTimeValueModel *>(model)); //!!! hmm.
46
47 m_server = getServer(model,
48 channel,
49 windowType,
50 windowSize,
51 windowIncrement,
52 fftSize,
53 polar,
54 criteria,
55 fillFromFrame);
56
57 if (!m_server) return; // caller should check isOK()
58
59 int xratio = windowIncrement / m_server->getWindowIncrement();
60 int yratio = m_server->getFFTSize() / fftSize;
61
62 while (xratio > 1) {
63 if (xratio & 0x1) {
64 cerr << "ERROR: FFTModel: Window increment ratio "
65 << windowIncrement << " / "
66 << m_server->getWindowIncrement()
67 << " must be a power of two" << endl;
68 assert(!(xratio & 0x1));
69 }
70 ++m_xshift;
71 xratio >>= 1;
72 }
73
74 while (yratio > 1) {
75 if (yratio & 0x1) {
76 cerr << "ERROR: FFTModel: FFT size ratio "
77 << m_server->getFFTSize() << " / " << fftSize
78 << " must be a power of two" << endl;
79 assert(!(yratio & 0x1));
80 }
81 ++m_yshift;
82 yratio >>= 1;
83 }
84 } 47 }
85 48
86 FFTModel::~FFTModel() 49 FFTModel::~FFTModel()
87 { 50 {
88 if (m_server) FFTDataServer::releaseInstance(m_server);
89 } 51 }
90 52
91 void 53 void
92 FFTModel::sourceModelAboutToBeDeleted() 54 FFTModel::sourceModelAboutToBeDeleted()
93 { 55 {
94 if (m_sourceModel) { 56 if (m_model) {
95 cerr << "FFTModel[" << this << "]::sourceModelAboutToBeDeleted(" << m_sourceModel << ")" << endl; 57 cerr << "FFTModel[" << this << "]::sourceModelAboutToBeDeleted(" << m_model << ")" << endl;
96 if (m_server) { 58 m_model = 0;
97 FFTDataServer::releaseInstance(m_server); 59 }
98 m_server = 0;
99 }
100 FFTDataServer::modelAboutToBeDeleted(m_sourceModel);
101 }
102 }
103
104 FFTDataServer *
105 FFTModel::getServer(const DenseTimeValueModel *model,
106 int channel,
107 WindowType windowType,
108 int windowSize,
109 int windowIncrement,
110 int fftSize,
111 bool polar,
112 StorageAdviser::Criteria criteria,
113 sv_frame_t fillFromFrame)
114 {
115 // Obviously, an FFT model of channel C (where C != -1) of an
116 // aggregate model is the same as the FFT model of the appropriate
117 // channel of whichever model that aggregate channel is drawn
118 // from. We should use that model here, in case we already have
119 // the data for it or will be wanting the same data again later.
120
121 // If the channel is -1 (i.e. mixture of all channels), then we
122 // can't do this shortcut unless the aggregate model only has one
123 // channel or contains exactly all of the channels of a single
124 // other model. That isn't very likely -- if it were the case,
125 // why would we be using an aggregate model?
126
127 if (channel >= 0) {
128
129 const AggregateWaveModel *aggregate =
130 dynamic_cast<const AggregateWaveModel *>(model);
131
132 if (aggregate && channel < aggregate->getComponentCount()) {
133
134 AggregateWaveModel::ModelChannelSpec spec =
135 aggregate->getComponent(channel);
136
137 return getServer(spec.model,
138 spec.channel,
139 windowType,
140 windowSize,
141 windowIncrement,
142 fftSize,
143 polar,
144 criteria,
145 fillFromFrame);
146 }
147 }
148
149 // The normal case
150
151 return FFTDataServer::getFuzzyInstance(model,
152 channel,
153 windowType,
154 windowSize,
155 windowIncrement,
156 fftSize,
157 polar,
158 criteria,
159 fillFromFrame);
160 }
161
162 sv_samplerate_t
163 FFTModel::getSampleRate() const
164 {
165 return isOK() ? m_server->getModel()->getSampleRate() : 0;
166 }
167
168 FFTModel::Column
169 FFTModel::getColumn(int x) const
170 {
171 Profiler profiler("FFTModel::getColumn", false);
172
173 Column result;
174
175 result.clear();
176 int h = getHeight();
177 result.reserve(h);
178
179 #ifdef __GNUC__
180 float magnitudes[h];
181 #else
182 float *magnitudes = (float *)alloca(h * sizeof(float));
183 #endif
184
185 if (m_server->getMagnitudesAt(x << m_xshift, magnitudes)) {
186
187 for (int y = 0; y < h; ++y) {
188 result.push_back(magnitudes[y]);
189 }
190
191 } else {
192 for (int i = 0; i < h; ++i) result.push_back(0.f);
193 }
194
195 return result;
196 } 60 }
197 61
198 QString 62 QString
199 FFTModel::getBinName(int n) const 63 FFTModel::getBinName(int n) const
200 { 64 {
207 bool 71 bool
208 FFTModel::estimateStableFrequency(int x, int y, double &frequency) 72 FFTModel::estimateStableFrequency(int x, int y, double &frequency)
209 { 73 {
210 if (!isOK()) return false; 74 if (!isOK()) return false;
211 75
212 sv_samplerate_t sampleRate = m_server->getModel()->getSampleRate(); 76 frequency = double(y * getSampleRate()) / m_fftSize;
213
214 int fftSize = m_server->getFFTSize() >> m_yshift;
215 frequency = double(y * sampleRate) / fftSize;
216 77
217 if (x+1 >= getWidth()) return false; 78 if (x+1 >= getWidth()) return false;
218 79
219 // At frequency f, a phase shift of 2pi (one cycle) happens in 1/f sec. 80 // At frequency f, a phase shift of 2pi (one cycle) happens in 1/f sec.
220 // At hopsize h and sample rate sr, one hop happens in h/sr sec. 81 // At hopsize h and sample rate sr, one hop happens in h/sr sec.
228 double oldPhase = getPhaseAt(x, y); 89 double oldPhase = getPhaseAt(x, y);
229 double newPhase = getPhaseAt(x+1, y); 90 double newPhase = getPhaseAt(x+1, y);
230 91
231 int incr = getResolution(); 92 int incr = getResolution();
232 93
233 double expectedPhase = oldPhase + (2.0 * M_PI * y * incr) / fftSize; 94 double expectedPhase = oldPhase + (2.0 * M_PI * y * incr) / m_fftSize;
234 95
235 double phaseError = princarg(newPhase - expectedPhase); 96 double phaseError = princarg(newPhase - expectedPhase);
236
237 // bool stable = (fabsf(phaseError) < (1.1f * (m_windowIncrement * M_PI) / m_fftSize));
238 97
239 // The new frequency estimate based on the phase error resulting 98 // The new frequency estimate based on the phase error resulting
240 // from assuming the "native" frequency of this bin 99 // from assuming the "native" frequency of this bin
241 100
242 frequency = 101 frequency =
243 (sampleRate * (expectedPhase + phaseError - oldPhase)) / 102 (getSampleRate() * (expectedPhase + phaseError - oldPhase)) /
244 (2.0 * M_PI * incr); 103 (2.0 * M_PI * incr);
245 104
246 return true; 105 return true;
247 } 106 }
248 107
291 // exceed the median. For pitch adaptivity, we adjust the window 150 // exceed the median. For pitch adaptivity, we adjust the window
292 // size to a roughly constant pitch range (about four tones). 151 // size to a roughly constant pitch range (about four tones).
293 152
294 sv_samplerate_t sampleRate = getSampleRate(); 153 sv_samplerate_t sampleRate = getSampleRate();
295 154
296 std::deque<float> window; 155 deque<float> window;
297 std::vector<int> inrange; 156 vector<int> inrange;
298 float dist = 0.5; 157 float dist = 0.5;
299 158
300 int medianWinSize = getPeakPickWindowSize(type, sampleRate, ymin, dist); 159 int medianWinSize = getPeakPickWindowSize(type, sampleRate, ymin, dist);
301 int halfWin = medianWinSize/2; 160 int halfWin = medianWinSize/2;
302 161
329 if (type == MajorPitchAdaptivePeaks) { 188 if (type == MajorPitchAdaptivePeaks) {
330 if (ymax + halfWin < values.size()) binmax = ymax + halfWin; 189 if (ymax + halfWin < values.size()) binmax = ymax + halfWin;
331 else binmax = values.size()-1; 190 else binmax = values.size()-1;
332 } 191 }
333 192
334 std::deque<float> sorted(window); 193 deque<float> sorted(window);
335 std::sort(sorted.begin(), sorted.end()); 194 sort(sorted.begin(), sorted.end());
336 float median = sorted[int(float(sorted.size()) * dist)]; 195 float median = sorted[int(float(sorted.size()) * dist)];
337 196
338 int centrebin = 0; 197 int centrebin = 0;
339 if (bin > actualSize/2) centrebin = bin - actualSize/2; 198 if (bin > actualSize/2) centrebin = bin - actualSize/2;
340 199
378 { 237 {
379 percentile = 0.5; 238 percentile = 0.5;
380 if (type == MajorPeaks) return 10; 239 if (type == MajorPeaks) return 10;
381 if (bin == 0) return 3; 240 if (bin == 0) return 3;
382 241
383 int fftSize = m_server->getFFTSize() >> m_yshift; 242 double binfreq = (getSampleRate() * bin) / m_fftSize;
384 double binfreq = (sampleRate * bin) / fftSize;
385 double hifreq = Pitch::getFrequencyForPitch(73, 0, binfreq); 243 double hifreq = Pitch::getFrequencyForPitch(73, 0, binfreq);
386 244
387 int hibin = int(lrint((hifreq * fftSize) / sampleRate)); 245 int hibin = int(lrint((hifreq * m_fftSize) / getSampleRate()));
388 int medianWinSize = hibin - bin; 246 int medianWinSize = hibin - bin;
389 if (medianWinSize < 3) medianWinSize = 3; 247 if (medianWinSize < 3) medianWinSize = 3;
390 248
391 percentile = 0.5f + float(binfreq / sampleRate); 249 percentile = 0.5f + float(binfreq / getSampleRate());
392 250
393 return medianWinSize; 251 return medianWinSize;
394 } 252 }
395 253
396 FFTModel::PeakSet 254 FFTModel::PeakSet
402 PeakSet peaks; 260 PeakSet peaks;
403 if (!isOK()) return peaks; 261 if (!isOK()) return peaks;
404 PeakLocationSet locations = getPeaks(type, x, ymin, ymax); 262 PeakLocationSet locations = getPeaks(type, x, ymin, ymax);
405 263
406 sv_samplerate_t sampleRate = getSampleRate(); 264 sv_samplerate_t sampleRate = getSampleRate();
407 int fftSize = m_server->getFFTSize() >> m_yshift;
408 int incr = getResolution(); 265 int incr = getResolution();
409 266
410 // This duplicates some of the work of estimateStableFrequency to 267 // This duplicates some of the work of estimateStableFrequency to
411 // allow us to retrieve the phases in two separate vertical 268 // allow us to retrieve the phases in two separate vertical
412 // columns, instead of jumping back and forth between columns x and 269 // columns, instead of jumping back and forth between columns x and
413 // x+1, which may be significantly slower if re-seeking is needed 270 // x+1, which may be significantly slower if re-seeking is needed
414 271
415 std::vector<float> phases; 272 vector<float> phases;
416 for (PeakLocationSet::iterator i = locations.begin(); 273 for (PeakLocationSet::iterator i = locations.begin();
417 i != locations.end(); ++i) { 274 i != locations.end(); ++i) {
418 phases.push_back(getPhaseAt(x, *i)); 275 phases.push_back(getPhaseAt(x, *i));
419 } 276 }
420 277
421 int phaseIndex = 0; 278 int phaseIndex = 0;
422 for (PeakLocationSet::iterator i = locations.begin(); 279 for (PeakLocationSet::iterator i = locations.begin();
423 i != locations.end(); ++i) { 280 i != locations.end(); ++i) {
424 double oldPhase = phases[phaseIndex]; 281 double oldPhase = phases[phaseIndex];
425 double newPhase = getPhaseAt(x+1, *i); 282 double newPhase = getPhaseAt(x+1, *i);
426 double expectedPhase = oldPhase + (2.0 * M_PI * *i * incr) / fftSize; 283 double expectedPhase = oldPhase + (2.0 * M_PI * *i * incr) / m_fftSize;
427 double phaseError = princarg(newPhase - expectedPhase); 284 double phaseError = princarg(newPhase - expectedPhase);
428 double frequency = 285 double frequency =
429 (sampleRate * (expectedPhase + phaseError - oldPhase)) 286 (sampleRate * (expectedPhase + phaseError - oldPhase))
430 / (2 * M_PI * incr); 287 / (2 * M_PI * incr);
431 // bool stable = (fabsf(phaseError) < (1.1f * (incr * M_PI) / fftSize));
432 // if (stable)
433 peaks[*i] = frequency; 288 peaks[*i] = frequency;
434 ++phaseIndex; 289 ++phaseIndex;
435 } 290 }
436 291
437 return peaks; 292 return peaks;
438 } 293 }
439 294
440 FFTModel::FFTModel(const FFTModel &model) :
441 DenseThreeDimensionalModel(),
442 m_server(model.m_server),
443 m_xshift(model.m_xshift),
444 m_yshift(model.m_yshift)
445 {
446 FFTDataServer::claimInstance(m_server);
447 }
448