libs/rubberband/src/StretchCalculator.cpp

   1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
   2
   3 /*
   4     Rubber Band
   5     An audio time-stretching and pitch-shifting library.
   6     Copyright 2007-2008 Chris Cannam.
   7
   8     This program is free software; you can redistribute it and/or
   9     modify it under the terms of the GNU General Public License as
  10     published by the Free Software Foundation; either version 2 of the
  11     License, or (at your option) any later version.  See the file
  12     COPYING included with this distribution for more information.
  13 */
  14
  15 #include "StretchCalculator.h"
  16
  17 #include <algorithm>
  18 #include <math.h>
  19 #include <algorithm>
  20 #include <iostream>
  21 #include <deque>
  22 #include <set>
  23 #include <cassert>
  24 #include <algorithm>
  25
  26 #include "sysutils.h"
  27
  28 namespace RubberBand
  29 {
  30
  31 StretchCalculator::StretchCalculator(size_t sampleRate,
  32                                      size_t inputIncrement,
  33                                      bool useHardPeaks) :
  34     m_sampleRate(sampleRate),
  35     m_increment(inputIncrement),
  36     m_prevDf(0),
  37     m_divergence(0),
  38     m_recovery(0),
  39     m_prevRatio(1.0),
  40     m_transientAmnesty(0),
  41     m_useHardPeaks(useHardPeaks)
  42 {
  43 //    std::cerr << "StretchCalculator::StretchCalculator: useHardPeaks = " << useHardPeaks << std::endl;
  44 }
  45
  46 StretchCalculator::~StretchCalculator()
  47 {
  48 }
  49
  50 std::vector<int>
  51 StretchCalculator::calculate(double ratio, size_t inputDuration,
  52                              const std::vector<float> &phaseResetDf,
  53                              const std::vector<float> &stretchDf)
  54 {
  55     assert(phaseResetDf.size() == stretchDf.size());
  56
  57     m_lastPeaks = findPeaks(phaseResetDf);
  58     std::vector<Peak> &peaks = m_lastPeaks;
  59     size_t totalCount = phaseResetDf.size();
  60
  61     std::vector<int> increments;
  62
  63     size_t outputDuration = lrint(inputDuration * ratio);
  64
  65     if (m_debugLevel > 0) {
  66         std::cerr << "StretchCalculator::calculate(): inputDuration " << inputDuration << ", ratio " << ratio << ", outputDuration " << outputDuration;
  67     }
  68
  69     outputDuration = lrint((phaseResetDf.size() * m_increment) * ratio);
  70
  71     if (m_debugLevel > 0) {
  72         std::cerr << " (rounded up to " << outputDuration << ")";
  73         std::cerr << ", df size " << phaseResetDf.size() << std::endl;
  74     }
  75
  76     std::vector<size_t> fixedAudioChunks;
  77     for (size_t i = 0; i < peaks.size(); ++i) {
  78         fixedAudioChunks.push_back
  79             (lrint((double(peaks[i].chunk) * outputDuration) / totalCount));
  80     }
  81
  82     if (m_debugLevel > 1) {
  83         std::cerr << "have " << peaks.size() << " fixed positions" << std::endl;
  84     }
  85
  86     size_t totalInput = 0, totalOutput = 0;
  87
  88     // For each region between two consecutive time sync points, we
  89     // want to take the number of output chunks to be allocated and
  90     // the detection function values within the range, and produce a
  91     // series of increments that sum to the number of output chunks,
  92     // such that each increment is displaced from the input increment
  93     // by an amount inversely proportional to the magnitude of the
  94     // stretch detection function at that input step.
  95
  96     size_t regionTotalChunks = 0;
  97
  98     for (size_t i = 0; i <= peaks.size(); ++i) {
  99
 100         size_t regionStart, regionStartChunk, regionEnd, regionEndChunk;
 101         bool phaseReset = false;
 102
 103         if (i == 0) {
 104             regionStartChunk = 0;
 105             regionStart = 0;
 106         } else {
 107             regionStartChunk = peaks[i-1].chunk;
 108             regionStart = fixedAudioChunks[i-1];
 109             phaseReset = peaks[i-1].hard;
 110         }
 111
 112         if (i == peaks.size()) {
 113             regionEndChunk = totalCount;
 114             regionEnd = outputDuration;
 115         } else {
 116             regionEndChunk = peaks[i].chunk;
 117             regionEnd = fixedAudioChunks[i];
 118         }
 119
 120         size_t regionDuration = regionEnd - regionStart;
 121         regionTotalChunks += regionDuration;
 122
 123         std::vector<float> dfRegion;
 124
 125         for (size_t j = regionStartChunk; j != regionEndChunk; ++j) {
 126             dfRegion.push_back(stretchDf[j]);
 127         }
 128
 129         if (m_debugLevel > 1) {
 130             std::cerr << "distributeRegion from " << regionStartChunk << " to " << regionEndChunk << " (chunks " << regionStart << " to " << regionEnd << ")" << std::endl;
 131         }
 132
 133         dfRegion = smoothDF(dfRegion);
 134
 135         std::vector<int> regionIncrements = distributeRegion
 136             (dfRegion, regionDuration, ratio, phaseReset);
 137
 138         size_t totalForRegion = 0;
 139
 140         for (size_t j = 0; j < regionIncrements.size(); ++j) {
 141
 142             int incr = regionIncrements[j];
 143
 144             if (j == 0 && phaseReset) increments.push_back(-incr);
 145             else increments.push_back(incr);
 146
 147             if (incr > 0) totalForRegion += incr;
 148             else totalForRegion += -incr;
 149
 150             totalInput += m_increment;
 151         }
 152
 153         if (totalForRegion != regionDuration) {
 154             std::cerr << "*** WARNING: distributeRegion returned wrong duration " << totalForRegion << ", expected " << regionDuration << std::endl;
 155         }
 156
 157         totalOutput += totalForRegion;
 158     }
 159
 160     if (m_debugLevel > 0) {
 161         std::cerr << "total input increment = " << totalInput << " (= " << totalInput / m_increment << " chunks), output = " << totalOutput << ", ratio = " << double(totalOutput)/double(totalInput) << ", ideal output " << size_t(ceil(totalInput * ratio)) << std::endl;
 162         std::cerr << "(region total = " << regionTotalChunks << ")" << std::endl;
 163     }
 164
 165     return increments;
 166 }
 167
 168 int
 169 StretchCalculator::calculateSingle(double ratio,
 170                                    float df,
 171                                    size_t increment)
 172 {
 173     if (increment == 0) increment = m_increment;
 174
 175     bool isTransient = false;
 176
 177     // We want to ensure, as close as possible, that the phase reset
 178     // points appear at _exactly_ the right audio frame numbers.
 179
 180     // In principle, the threshold depends on chunk size: larger chunk
 181     // sizes need higher thresholds.  Since chunk size depends on
 182     // ratio, I suppose we could in theory calculate the threshold
 183     // from the ratio directly.  For the moment we're happy if it
 184     // works well in common situations.
 185
 186     float transientThreshold = 0.35f;
 187     if (ratio > 1) transientThreshold = 0.25f;
 188
 189     if (m_useHardPeaks && df > m_prevDf * 1.1f && df > transientThreshold) {
 190         isTransient = true;
 191     }
 192
 193     if (m_debugLevel > 2) {
 194         std::cerr << "df = " << df << ", prevDf = " << m_prevDf
 195                   << ", thresh = " << transientThreshold << std::endl;
 196     }
 197
 198     m_prevDf = df;
 199
 200     bool ratioChanged = (ratio != m_prevRatio);
 201     m_prevRatio = ratio;
 202
 203     if (isTransient && m_transientAmnesty == 0) {
 204         if (m_debugLevel > 1) {
 205             std::cerr << "StretchCalculator::calculateSingle: transient"
 206                       << std::endl;
 207         }
 208         m_divergence += increment - (increment * ratio);
 209
 210         // as in offline mode, 0.05 sec approx min between transients
 211         m_transientAmnesty =
 212             lrint(ceil(double(m_sampleRate) / (20 * double(increment))));
 213
 214         m_recovery = m_divergence / ((m_sampleRate / 10.0) / increment);
 215         return -int(increment);
 216     }
 217
 218     if (ratioChanged) {
 219         m_recovery = m_divergence / ((m_sampleRate / 10.0) / increment);
 220     }
 221
 222     if (m_transientAmnesty > 0) --m_transientAmnesty;
 223
 224     int incr = lrint(increment * ratio - m_recovery);
 225     if (m_debugLevel > 2 || (m_debugLevel > 1 && m_divergence != 0)) {
 226         std::cerr << "divergence = " << m_divergence << ", recovery = " << m_recovery << ", incr = " << incr << ", ";
 227     }
 228     if (incr < lrint((increment * ratio) / 2)) {
 229         incr = lrint((increment * ratio) / 2);
 230     } else if (incr > lrint(increment * ratio * 2)) {
 231         incr = lrint(increment * ratio * 2);
 232     }
 233
 234     double divdiff = (increment * ratio) - incr;
 235
 236     if (m_debugLevel > 2 || (m_debugLevel > 1 && m_divergence != 0)) {
 237         std::cerr << "divdiff = " << divdiff << std::endl;
 238     }
 239
 240     double prevDivergence = m_divergence;
 241     m_divergence -= divdiff;
 242     if ((prevDivergence < 0 && m_divergence > 0) ||
 243         (prevDivergence > 0 && m_divergence < 0)) {
 244         m_recovery = m_divergence / ((m_sampleRate / 10.0) / increment);
 245     }
 246
 247     return incr;
 248 }
 249
 250 void
 251 StretchCalculator::reset()
 252 {
 253     m_prevDf = 0;
 254     m_divergence = 0;
 255 }
 256
 257 std::vector<StretchCalculator::Peak>
 258 StretchCalculator::findPeaks(const std::vector<float> &rawDf)
 259 {
 260     std::vector<float> df = smoothDF(rawDf);
 261
 262     // We distinguish between "soft" and "hard" peaks.  A soft peak is
 263     // simply the result of peak-picking on the smoothed onset
 264     // detection function, and it represents any (strong-ish) onset.
 265     // We aim to ensure always that soft peaks are placed at the
 266     // correct position in time.  A hard peak is where there is a very
 267     // rapid rise in detection function, and it presumably represents
 268     // a more broadband, noisy transient.  For these we perform a
 269     // phase reset (if in the appropriate mode), and we locate the
 270     // reset at the first point where we notice enough of a rapid
 271     // rise, rather than necessarily at the peak itself, in order to
 272     // preserve the shape of the transient.
 273
 274     std::set<size_t> hardPeakCandidates;
 275     std::set<size_t> softPeakCandidates;
 276
 277     if (m_useHardPeaks) {
 278
 279         // 0.05 sec approx min between hard peaks
 280         size_t hardPeakAmnesty = lrint(ceil(double(m_sampleRate) /
 281                                             (20 * double(m_increment))));
 282         size_t prevHardPeak = 0;
 283
 284         if (m_debugLevel > 1) {
 285             std::cerr << "hardPeakAmnesty = " << hardPeakAmnesty << std::endl;
 286         }
 287
 288         for (size_t i = 1; i + 1 < df.size(); ++i) {
 289
 290             if (df[i] < 0.1) continue;
 291             if (df[i] <= df[i-1] * 1.1) continue;
 292             if (df[i] < 0.22) continue;
 293
 294             if (!hardPeakCandidates.empty() &&
 295                 i < prevHardPeak + hardPeakAmnesty) {
 296                 continue;
 297             }
 298
 299             bool hard = (df[i] > 0.4);
 300
 301             if (hard && (m_debugLevel > 1)) {
 302                 std::cerr << "hard peak at " << i << ": " << df[i]
 303                           << " > absolute " << 0.4
 304                           << std::endl;
 305             }
 306
 307             if (!hard) {
 308                 hard = (df[i] > df[i-1] * 1.4);
 309
 310                 if (hard && (m_debugLevel > 1)) {
 311                     std::cerr << "hard peak at " << i << ": " << df[i]
 312                               << " > prev " << df[i-1] << " * 1.4"
 313                               << std::endl;
 314                 }
 315             }
 316
 317             if (!hard && i > 1) {
 318                 hard = (df[i]   > df[i-1] * 1.2 &&
 319                         df[i-1] > df[i-2] * 1.2);
 320
 321                 if (hard && (m_debugLevel > 1)) {
 322                     std::cerr << "hard peak at " << i << ": " << df[i]
 323                               << " > prev " << df[i-1] << " * 1.2 and "
 324                               << df[i-1] << " > prev " << df[i-2] << " * 1.2"
 325                               << std::endl;
 326                 }
 327             }
 328
 329             if (!hard && i > 2) {
 330                 // have already established that df[i] > df[i-1] * 1.1
 331                 hard = (df[i] > 0.3 &&
 332                         df[i-1] > df[i-2] * 1.1 &&
 333                         df[i-2] > df[i-3] * 1.1);
 334
 335                 if (hard && (m_debugLevel > 1)) {
 336                     std::cerr << "hard peak at " << i << ": " << df[i]
 337                               << " > prev " << df[i-1] << " * 1.1 and "
 338                               << df[i-1] << " > prev " << df[i-2] << " * 1.1 and "
 339                               << df[i-2] << " > prev " << df[i-3] << " * 1.1"
 340                               << std::endl;
 341                 }
 342             }
 343
 344             if (!hard) continue;
 345
 346 //            (df[i+1] > df[i] && df[i+1] > df[i-1] * 1.8) ||
 347 //                df[i] > 0.4) {
 348
 349             size_t peakLocation = i;
 350
 351             if (i + 1 < rawDf.size() &&
 352                 rawDf[i + 1] > rawDf[i] * 1.4) {
 353
 354                 ++peakLocation;
 355
 356                 if (m_debugLevel > 1) {
 357                     std::cerr << "pushing hard peak forward to " << peakLocation << ": " << df[peakLocation] << " > " << df[peakLocation-1] << " * " << 1.4 << std::endl;
 358                 }
 359             }
 360
 361             hardPeakCandidates.insert(peakLocation);
 362             prevHardPeak = peakLocation;
 363         }
 364     }
 365
 366     size_t medianmaxsize = lrint(ceil(double(m_sampleRate) /
 367                                  double(m_increment))); // 1 sec ish
 368
 369     if (m_debugLevel > 1) {
 370         std::cerr << "mediansize = " << medianmaxsize << std::endl;
 371     }
 372     if (medianmaxsize < 7) {
 373         medianmaxsize = 7;
 374         if (m_debugLevel > 1) {
 375             std::cerr << "adjusted mediansize = " << medianmaxsize << std::endl;
 376         }
 377     }
 378
 379     int minspacing = lrint(ceil(double(m_sampleRate) /
 380                                 (20 * double(m_increment)))); // 0.05 sec ish
 381
 382     std::deque<float> medianwin;
 383     std::vector<float> sorted;
 384     int softPeakAmnesty = 0;
 385
 386     for (size_t i = 0; i < medianmaxsize/2; ++i) {
 387         medianwin.push_back(0);
 388     }
 389     for (size_t i = 0; i < medianmaxsize/2 && i < df.size(); ++i) {
 390         medianwin.push_back(df[i]);
 391     }
 392
 393     size_t lastSoftPeak = 0;
 394
 395     for (size_t i = 0; i < df.size(); ++i) {
 396
 397         size_t mediansize = medianmaxsize;
 398
 399         if (medianwin.size() < mediansize) {
 400             mediansize = medianwin.size();
 401         }
 402
 403         size_t middle = medianmaxsize / 2;
 404         if (middle >= mediansize) middle = mediansize-1;
 405
 406         size_t nextDf = i + mediansize - middle;
 407
 408         if (mediansize < 2) {
 409             if (mediansize > medianmaxsize) { // absurd, but never mind that
 410                 medianwin.pop_front();
 411             }
 412             if (nextDf < df.size()) {
 413                 medianwin.push_back(df[nextDf]);
 414             } else {
 415                 medianwin.push_back(0);
 416             }
 417             continue;
 418         }
 419
 420         if (m_debugLevel > 2) {
 421 //            std::cerr << "have " << mediansize << " in median buffer" << std::endl;
 422         }
 423
 424         sorted.clear();
 425         for (size_t j = 0; j < mediansize; ++j) {
 426             sorted.push_back(medianwin[j]);
 427         }
 428         std::sort(sorted.begin(), sorted.end());
 429
 430         size_t n = 90; // percentile above which we pick peaks
 431         size_t index = (sorted.size() * n) / 100;
 432         if (index >= sorted.size()) index = sorted.size()-1;
 433         if (index == sorted.size()-1 && index > 0) --index;
 434         float thresh = sorted[index];
 435
 436 //        if (m_debugLevel > 2) {
 437 //            std::cerr << "medianwin[" << middle << "] = " << medianwin[middle] << ", thresh = " << thresh << std::endl;
 438 //            if (medianwin[middle] == 0.f) {
 439 //                std::cerr << "contents: ";
 440 //                for (size_t j = 0; j < medianwin.size(); ++j) {
 441 //                    std::cerr << medianwin[j] << " ";
 442 //                }
 443 //                std::cerr << std::endl;
 444 //            }
 445 //        }
 446
 447         if (medianwin[middle] > thresh &&
 448             medianwin[middle] > medianwin[middle-1] &&
 449             medianwin[middle] > medianwin[middle+1] &&
 450             softPeakAmnesty == 0) {
 451
 452             size_t maxindex = middle;
 453             float maxval = medianwin[middle];
 454
 455             for (size_t j = middle+1; j < mediansize; ++j) {
 456                 if (medianwin[j] > maxval) {
 457                     maxval = medianwin[j];
 458                     maxindex = j;
 459                 } else if (medianwin[j] < medianwin[middle]) {
 460                     break;
 461                 }
 462             }
 463
 464             size_t peak = i + maxindex - middle;
 465
 466 //            std::cerr << "i = " << i << ", maxindex = " << maxindex << ", middle = " << middle << ", so peak at " << peak << std::endl;
 467
 468             if (softPeakCandidates.empty() || lastSoftPeak != peak) {
 469
 470                 if (m_debugLevel > 1) {
 471                     std::cerr << "soft peak at " << peak << " ("
 472                               << peak * m_increment << "): "
 473                               << medianwin[middle] << " > "
 474                               << thresh << " and "
 475                               << medianwin[middle]
 476                               << " > " << medianwin[middle-1] << " and "
 477                               << medianwin[middle]
 478                               << " > " << medianwin[middle+1]
 479                               << std::endl;
 480                 }
 481
 482                 if (peak >= df.size()) {
 483                     if (m_debugLevel > 2) {
 484                         std::cerr << "peak is beyond end"  << std::endl;
 485                     }
 486                 } else {
 487                     softPeakCandidates.insert(peak);
 488                     lastSoftPeak = peak;
 489                 }
 490             }
 491
 492             softPeakAmnesty = minspacing + maxindex - middle;
 493             if (m_debugLevel > 2) {
 494                 std::cerr << "amnesty = " << softPeakAmnesty << std::endl;
 495             }
 496
 497         } else if (softPeakAmnesty > 0) --softPeakAmnesty;
 498
 499         if (mediansize >= medianmaxsize) {
 500             medianwin.pop_front();
 501         }
 502         if (nextDf < df.size()) {
 503             medianwin.push_back(df[nextDf]);
 504         } else {
 505             medianwin.push_back(0);
 506         }
 507     }
 508
 509     std::vector<Peak> peaks;
 510
 511     while (!hardPeakCandidates.empty() || !softPeakCandidates.empty()) {
 512
 513         bool haveHardPeak = !hardPeakCandidates.empty();
 514         bool haveSoftPeak = !softPeakCandidates.empty();
 515
 516         size_t hardPeak = (haveHardPeak ? *hardPeakCandidates.begin() : 0);
 517         size_t softPeak = (haveSoftPeak ? *softPeakCandidates.begin() : 0);
 518
 519         Peak peak;
 520         peak.hard = false;
 521         peak.chunk = softPeak;
 522
 523         bool ignore = false;
 524
 525         if (haveHardPeak &&
 526             (!haveSoftPeak || hardPeak <= softPeak)) {
 527
 528             if (m_debugLevel > 2) {
 529                 std::cerr << "Hard peak: " << hardPeak << std::endl;
 530             }
 531
 532             peak.hard = true;
 533             peak.chunk = hardPeak;
 534             hardPeakCandidates.erase(hardPeakCandidates.begin());
 535
 536         } else {
 537             if (m_debugLevel > 2) {
 538                 std::cerr << "Soft peak: " << softPeak << std::endl;
 539             }
 540             if (!peaks.empty() &&
 541                 peaks[peaks.size()-1].hard &&
 542                 peaks[peaks.size()-1].chunk + 3 >= softPeak) {
 543                 if (m_debugLevel > 2) {
 544                     std::cerr << "(ignoring, as we just had a hard peak)"
 545                               << std::endl;
 546                 }
 547                 ignore = true;
 548             }
 549         }
 550
 551         if (haveSoftPeak && peak.chunk == softPeak) {
 552             softPeakCandidates.erase(softPeakCandidates.begin());
 553         }
 554
 555         if (!ignore) {
 556             peaks.push_back(peak);
 557         }
 558     }
 559
 560     return peaks;
 561 }
 562
 563 std::vector<float>
 564 StretchCalculator::smoothDF(const std::vector<float> &df)
 565 {
 566     std::vector<float> smoothedDF;
 567
 568     for (size_t i = 0; i < df.size(); ++i) {
 569         // three-value moving mean window for simple smoothing
 570         float total = 0.f, count = 0;
 571         if (i > 0) { total += df[i-1]; ++count; }
 572         total += df[i]; ++count;
 573         if (i+1 < df.size()) { total += df[i+1]; ++count; }
 574         float mean = total / count;
 575         smoothedDF.push_back(mean);
 576     }
 577
 578     return smoothedDF;
 579 }
 580
 581 std::vector<int>
 582 StretchCalculator::distributeRegion(const std::vector<float> &dfIn,
 583                                     size_t duration, float ratio, bool phaseReset)
 584 {
 585     std::vector<float> df(dfIn);
 586     std::vector<int> increments;
 587
 588     // The peak for the stretch detection function may appear after
 589     // the peak that we're using to calculate the start of the region.
 590     // We don't want that.  If we find a peak in the first half of
 591     // the region, we should set all the values up to that point to
 592     // the same value as the peak.
 593
 594     // (This might not be subtle enough, especially if the region is
 595     // long -- we want a bound that corresponds to acoustic perception
 596     // of the audible bounce.)
 597
 598     for (size_t i = 1; i < df.size()/2; ++i) {
 599         if (df[i] < df[i-1]) {
 600             if (m_debugLevel > 1) {
 601                 std::cerr << "stretch peak offset: " << i-1 << " (peak " << df[i-1] << ")" << std::endl;
 602             }
 603             for (size_t j = 0; j < i-1; ++j) {
 604                 df[j] = df[i-1];
 605             }
 606             break;
 607         }
 608     }
 609
 610     float maxDf = 0;
 611
 612     for (size_t i = 0; i < df.size(); ++i) {
 613         if (i == 0 || df[i] > maxDf) maxDf = df[i];
 614     }
 615
 616     // We want to try to ensure the last 100ms or so (if possible) are
 617     // tending back towards the maximum df, so that the stretchiness
 618     // reduces at the end of the stretched region.
 619
 620     int reducedRegion = lrint((0.1 * m_sampleRate) / m_increment);
 621     if (reducedRegion > int(df.size()/5)) reducedRegion = df.size()/5;
 622
 623     for (int i = 0; i < reducedRegion; ++i) {
 624         size_t index = df.size() - reducedRegion + i;
 625         df[index] = df[index] + ((maxDf - df[index]) * i) / reducedRegion;
 626     }
 627
 628     long toAllot = long(duration) - long(m_increment * df.size());
 629
 630     if (m_debugLevel > 1) {
 631         std::cerr << "region of " << df.size() << " chunks, output duration " << duration << ", toAllot " << toAllot << std::endl;
 632     }
 633
 634     size_t totalIncrement = 0;
 635
 636     // We place limits on the amount of displacement per chunk.  if
 637     // ratio < 0, no increment should be larger than increment*ratio
 638     // or smaller than increment*ratio/2; if ratio > 0, none should be
 639     // smaller than increment*ratio or larger than increment*ratio*2.
 640     // We need to enforce this in the assignment of displacements to
 641     // allotments, not by trying to respond if something turns out
 642     // wrong.
 643
 644     // Note that the ratio is only provided to this function for the
 645     // purposes of establishing this bound to the displacement.
 646
 647     // so if
 648     // maxDisplacement / totalDisplacement > increment * ratio*2 - increment
 649     // (for ratio > 1)
 650     // or
 651     // maxDisplacement / totalDisplacement < increment * ratio/2
 652     // (for ratio < 1)
 653
 654     // then we need to adjust and accommodate
 655
 656     bool acceptableSquashRange = false;
 657
 658     double totalDisplacement = 0;
 659     double maxDisplacement = 0; // min displacement will be 0 by definition
 660
 661     maxDf = 0;
 662     float adj = 0;
 663
 664     while (!acceptableSquashRange) {
 665
 666         acceptableSquashRange = true;
 667         calculateDisplacements(df, maxDf, totalDisplacement, maxDisplacement,
 668                                adj);
 669
 670         if (m_debugLevel > 1) {
 671             std::cerr << "totalDisplacement " << totalDisplacement << ", max " << maxDisplacement << " (maxDf " << maxDf << ", df count " << df.size() << ")" << std::endl;
 672         }
 673
 674         if (totalDisplacement == 0) {
 675 // Not usually a problem, in fact
 676 //            std::cerr << "WARNING: totalDisplacement == 0 (duration " << duration << ", " << df.size() << " values in df)" << std::endl;
 677             if (!df.empty() && adj == 0) {
 678                 acceptableSquashRange = false;
 679                 adj = 1;
 680             }
 681             continue;
 682         }
 683
 684         int extremeIncrement = m_increment + lrint((toAllot * maxDisplacement) / totalDisplacement);
 685         if (ratio < 1.0) {
 686             if (extremeIncrement > lrint(ceil(m_increment * ratio))) {
 687                 std::cerr << "ERROR: extreme increment " << extremeIncrement << " > " << m_increment * ratio << " (this should not happen)" << std::endl;
 688             } else if (extremeIncrement < (m_increment * ratio) / 2) {
 689                 if (m_debugLevel > 0) {
 690                     std::cerr << "WARNING: extreme increment " << extremeIncrement << " < " << (m_increment * ratio) / 2 << std::endl;
 691                 }
 692                 acceptableSquashRange = false;
 693             }
 694         } else {
 695             if (extremeIncrement > m_increment * ratio * 2) {
 696                 if (m_debugLevel > 0) {
 697                     std::cerr << "WARNING: extreme increment " << extremeIncrement << " > " << m_increment * ratio * 2 << std::endl;
 698                 }
 699                 acceptableSquashRange = false;
 700             } else if (extremeIncrement < lrint(floor(m_increment * ratio))) {
 701                 std::cerr << "ERROR: extreme increment " << extremeIncrement << " < " << m_increment * ratio << " (I thought this couldn't happen?)" << std::endl;
 702             }
 703         }
 704
 705         if (!acceptableSquashRange) {
 706             // Need to make maxDisplacement smaller as a proportion of
 707             // the total displacement, yet ensure that the
 708             // displacements still sum to the total.
 709             adj += maxDf/10;
 710         }
 711     }
 712
 713     for (size_t i = 0; i < df.size(); ++i) {
 714
 715         double displacement = maxDf - df[i];
 716         if (displacement < 0) displacement -= adj;
 717         else displacement += adj;
 718
 719         if (i == 0 && phaseReset) {
 720             if (df.size() == 1) {
 721                 increments.push_back(duration);
 722                 totalIncrement += duration;
 723             } else {
 724                 increments.push_back(m_increment);
 725                 totalIncrement += m_increment;
 726             }
 727             totalDisplacement -= displacement;
 728             continue;
 729         }
 730
 731         double theoreticalAllotment = 0;
 732
 733         if (totalDisplacement != 0) {
 734             theoreticalAllotment = (toAllot * displacement) / totalDisplacement;
 735         }
 736         int allotment = lrint(theoreticalAllotment);
 737         if (i + 1 == df.size()) allotment = toAllot;
 738
 739         int increment = m_increment + allotment;
 740
 741         if (increment <= 0) {
 742             // this is a serious problem, the allocation is quite
 743             // wrong if it allows increment to diverge so far from the
 744             // input increment
 745             std::cerr << "*** WARNING: increment " << increment << " <= 0, rounding to zero" << std::endl;
 746             increment = 0;
 747             allotment = increment - m_increment;
 748         }
 749
 750         increments.push_back(increment);
 751         totalIncrement += increment;
 752
 753         toAllot -= allotment;
 754         totalDisplacement -= displacement;
 755
 756         if (m_debugLevel > 2) {
 757             std::cerr << "df " << df[i] << ", smoothed " << df[i] << ", disp " << displacement << ", allot " << theoreticalAllotment << ", incr " << increment << ", remain " << toAllot << std::endl;
 758         }
 759     }
 760
 761     if (m_debugLevel > 2) {
 762         std::cerr << "total increment: " << totalIncrement << ", left over: " << toAllot << " to allot, displacement " << totalDisplacement << std::endl;
 763     }
 764
 765     if (totalIncrement != duration) {
 766         std::cerr << "*** WARNING: calculated output duration " << totalIncrement << " != expected " << duration << std::endl;
 767     }
 768
 769     return increments;
 770 }
 771
 772 void
 773 StretchCalculator::calculateDisplacements(const std::vector<float> &df,
 774                                           float &maxDf,
 775                                           double &totalDisplacement,
 776                                           double &maxDisplacement,
 777                                           float adj) const
 778 {
 779     totalDisplacement = maxDisplacement = 0;
 780
 781     maxDf = 0;
 782
 783     for (size_t i = 0; i < df.size(); ++i) {
 784         if (i == 0 || df[i] > maxDf) maxDf = df[i];
 785     }
 786
 787     for (size_t i = 0; i < df.size(); ++i) {
 788         double displacement = maxDf - df[i];
 789         if (displacement < 0) displacement -= adj;
 790         else displacement += adj;
 791         totalDisplacement += displacement;
 792         if (i == 0 || displacement > maxDisplacement) {
 793             maxDisplacement = displacement;
 794         }
 795     }
 796 }
 797
 798 }
 799