Mercurial > hg > svcore
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 |