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 | 
