1#ifndef SIGNALSMITH_STRETCH_H
2#define SIGNALSMITH_STRETCH_H
3
4#include "dsp/spectral.h"
5#include "dsp/delay.h"
6#include "dsp/perf.h"
7SIGNALSMITH_DSP_VERSION_CHECK(1, 6, 0); // Check version is compatible
8#include <vector>
9#include <algorithm>
10#include <functional>
11#include <random>
12
13namespace signalsmith { namespace stretch {
14
15template<typename Sample=float>
16struct SignalsmithStretch {
17
18 SignalsmithStretch() : randomEngine(std::random_device{}()) {}
19 SignalsmithStretch(long seed) : randomEngine(seed) {}
20
21 int blockSamples() const {
22 return stft.windowSize();
23 }
24 int intervalSamples() const {
25 return stft.interval();
26 }
27 int inputLatency() const {
28 return stft.windowSize()/2;
29 }
30 int outputLatency() const {
31 return stft.windowSize() - inputLatency();
32 }
33
34 void reset() {
35 stft.reset();
36 inputBuffer.reset();
37 prevInputOffset = -1;
38 channelBands.assign(channelBands.size(), Band());
39 silenceCounter = 2*stft.windowSize();
40 didSeek = false;
41 flushed = true;
42 }
43
44 // Configures using a default preset
45 void presetDefault(int nChannels, Sample sampleRate) {
46 configure(nChannels, blockSamples: sampleRate*0.12, intervalSamples: sampleRate*0.03);
47 }
48 void presetCheaper(int nChannels, Sample sampleRate) {
49 configure(nChannels, blockSamples: sampleRate*0.1, intervalSamples: sampleRate*0.04);
50 }
51
52 // Manual setup
53 void configure(int nChannels, int blockSamples, int intervalSamples) {
54 channels = nChannels;
55 stft.resize(channels, blockSamples, intervalSamples);
56 bands = stft.bands();
57 inputBuffer.resize(channels, blockSamples + intervalSamples + 1);
58 timeBuffer.assign(stft.fftSize(), 0);
59 channelBands.assign(bands*channels, Band());
60
61 // Various phase rotations
62 rotCentreSpectrum.resize(bands);
63 rotPrevInterval.assign(bands, 0);
64 timeShiftPhases(shiftSamples: blockSamples*Sample(-0.5), output&: rotCentreSpectrum);
65 timeShiftPhases(shiftSamples: -intervalSamples, output&: rotPrevInterval);
66 peaks.reserve(bands);
67 energy.resize(bands);
68 smoothedEnergy.resize(bands);
69 outputMap.resize(bands);
70 channelPredictions.resize(channels*bands);
71 }
72
73 /// Frequency multiplier, and optional tonality limit (as multiple of sample-rate)
74 void setTransposeFactor(Sample multiplier, Sample tonalityLimit=0) {
75 freqMultiplier = multiplier;
76 if (tonalityLimit > 0) {
77 freqTonalityLimit = tonalityLimit/std::sqrt(multiplier); // compromise between input and output limits
78 } else {
79 freqTonalityLimit = 1;
80 }
81 customFreqMap = nullptr;
82 }
83 void setTransposeSemitones(Sample semitones, Sample tonalityLimit=0) {
84 setTransposeFactor(multiplier: std::pow(2, semitones/12), tonalityLimit);
85 customFreqMap = nullptr;
86 }
87 // Sets a custom frequency map - should be monotonically increasing
88 void setFreqMap(std::function<Sample(Sample)> inputToOutput) {
89 customFreqMap = inputToOutput;
90 }
91
92 // Provide previous input ("pre-roll"), without affecting the speed calculation. You should ideally feed it one block-length + one interval
93 template<class Inputs>
94 void seek(Inputs &&inputs, int inputSamples, double playbackRate) {
95 inputBuffer.reset();
96 Sample totalEnergy = 0;
97 for (int c = 0; c < channels; ++c) {
98 auto &&inputChannel = inputs[c];
99 auto &&bufferChannel = inputBuffer[c];
100 int startIndex = std::max<int>(0, inputSamples - stft.windowSize() - stft.interval());
101 for (int i = startIndex; i < inputSamples; ++i) {
102 Sample s = inputChannel[i];
103 totalEnergy += s*s;
104 bufferChannel[i] = s;
105 }
106 }
107 if (totalEnergy >= noiseFloor) {
108 silenceCounter = 0;
109 silenceFirst = true;
110 }
111 inputBuffer += inputSamples;
112 didSeek = true;
113 seekTimeFactor = (playbackRate*stft.interval() > 1) ? 1/playbackRate : stft.interval();
114 }
115
116 template<class Inputs, class Outputs>
117 void process(Inputs &&inputs, int inputSamples, Outputs &&outputs, int outputSamples) {
118 Sample totalEnergy = 0;
119 for (int c = 0; c < channels; ++c) {
120 auto &&inputChannel = inputs[c];
121 for (int i = 0; i < inputSamples; ++i) {
122 Sample s = inputChannel[i];
123 totalEnergy += s*s;
124 }
125 }
126 if (totalEnergy < noiseFloor) {
127 if (silenceCounter >= 2*stft.windowSize()) {
128 if (silenceFirst) {
129 silenceFirst = false;
130 for (auto &b : channelBands) {
131 b.input = b.prevInput = b.output = b.prevOutput = 0;
132 b.inputEnergy = 0;
133 }
134 }
135
136 if (inputSamples > 0) {
137 // copy from the input, wrapping around if needed
138 for (int outputIndex = 0; outputIndex < outputSamples; ++outputIndex) {
139 int inputIndex = outputIndex%inputSamples;
140 for (int c = 0; c < channels; ++c) {
141 outputs[c][outputIndex] = inputs[c][inputIndex];
142 }
143 }
144 } else {
145 for (int c = 0; c < channels; ++c) {
146 auto &&outputChannel = outputs[c];
147 for (int outputIndex = 0; outputIndex < outputSamples; ++outputIndex) {
148 outputChannel[outputIndex] = 0;
149 }
150 }
151 }
152
153 // Store input in history buffer
154 for (int c = 0; c < channels; ++c) {
155 auto &&inputChannel = inputs[c];
156 auto &&bufferChannel = inputBuffer[c];
157 int startIndex = std::max<int>(0, inputSamples - stft.windowSize() - stft.interval());
158 for (int i = startIndex; i < inputSamples; ++i) {
159 bufferChannel[i] = inputChannel[i];
160 }
161 }
162 inputBuffer += inputSamples;
163 return;
164 } else {
165 silenceCounter += inputSamples;
166 }
167 } else {
168 silenceCounter = 0;
169 silenceFirst = true;
170 }
171
172 for (int outputIndex = 0; outputIndex < outputSamples; ++outputIndex) {
173 stft.ensureValid(outputIndex, [&](int outputOffset) {
174 // Time to process a spectrum! Where should it come from in the input?
175 int inputOffset = std::round(outputOffset*Sample(inputSamples)/outputSamples) - stft.windowSize();
176 int inputInterval = inputOffset - prevInputOffset;
177 prevInputOffset = inputOffset;
178
179 bool newSpectrum = didSeek || (inputInterval > 0);
180 if (newSpectrum) {
181 for (int c = 0; c < channels; ++c) {
182 // Copy from the history buffer, if needed
183 auto &&bufferChannel = inputBuffer[c];
184 for (int i = 0; i < -inputOffset; ++i) {
185 timeBuffer[i] = bufferChannel[i + inputOffset];
186 }
187 // Copy the rest from the input
188 auto &&inputChannel = inputs[c];
189 for (int i = std::max<int>(a: 0, b: -inputOffset); i < stft.windowSize(); ++i) {
190 timeBuffer[i] = inputChannel[i + inputOffset];
191 }
192 stft.analyse(c, timeBuffer);
193 }
194 flushed = false; // TODO: first block after a flush should be gain-compensated
195
196 for (int c = 0; c < channels; ++c) {
197 auto channelBands = bandsForChannel(channel: c);
198 auto &&spectrumBands = stft.spectrum[c];
199 for (int b = 0; b < bands; ++b) {
200 channelBands[b].input = signalsmith::perf::mul(spectrumBands[b], rotCentreSpectrum[b]);
201 }
202 }
203
204 if (didSeek || inputInterval != stft.interval()) { // make sure the previous input is the correct distance in the past
205 int prevIntervalOffset = inputOffset - stft.interval();
206 for (int c = 0; c < channels; ++c) {
207 // Copy from the history buffer, if needed
208 auto &&bufferChannel = inputBuffer[c];
209 for (int i = 0; i < std::min(-prevIntervalOffset, stft.windowSize()); ++i) {
210 timeBuffer[i] = bufferChannel[i + prevIntervalOffset];
211 }
212 // Copy the rest from the input
213 auto &&inputChannel = inputs[c];
214 for (int i = std::max<int>(a: 0, b: -prevIntervalOffset); i < stft.windowSize(); ++i) {
215 timeBuffer[i] = inputChannel[i + prevIntervalOffset];
216 }
217 stft.analyse(c, timeBuffer);
218 }
219 for (int c = 0; c < channels; ++c) {
220 auto channelBands = bandsForChannel(channel: c);
221 auto &&spectrumBands = stft.spectrum[c];
222 for (int b = 0; b < bands; ++b) {
223 channelBands[b].prevInput = signalsmith::perf::mul(spectrumBands[b], rotCentreSpectrum[b]);
224 }
225 }
226 }
227 }
228
229 Sample timeFactor = didSeek ? seekTimeFactor : stft.interval()/std::max<Sample>(1, inputInterval);
230 processSpectrum(newSpectrum, timeFactor);
231 didSeek = false;
232
233 for (int c = 0; c < channels; ++c) {
234 auto channelBands = bandsForChannel(channel: c);
235 auto &&spectrumBands = stft.spectrum[c];
236 for (int b = 0; b < bands; ++b) {
237 spectrumBands[b] = signalsmith::perf::mul<true>(channelBands[b].output, rotCentreSpectrum[b]);
238 }
239 }
240 });
241
242 for (int c = 0; c < channels; ++c) {
243 auto &&outputChannel = outputs[c];
244 auto &&stftChannel = stft[c];
245 outputChannel[outputIndex] = stftChannel[outputIndex];
246 }
247 }
248
249 // Store input in history buffer
250 for (int c = 0; c < channels; ++c) {
251 auto &&inputChannel = inputs[c];
252 auto &&bufferChannel = inputBuffer[c];
253 int startIndex = std::max<int>(0, inputSamples - stft.windowSize());
254 for (int i = startIndex; i < inputSamples; ++i) {
255 bufferChannel[i] = inputChannel[i];
256 }
257 }
258 inputBuffer += inputSamples;
259 stft += outputSamples;
260 prevInputOffset -= inputSamples;
261 }
262
263 // Read the remaining output, providing no further input. `outputSamples` should ideally be at least `.outputLatency()`
264 template<class Outputs>
265 void flush(Outputs &&outputs, int outputSamples) {
266 int plainOutput = std::min<int>(outputSamples, stft.windowSize());
267 int foldedBackOutput = std::min<int>(outputSamples, stft.windowSize() - plainOutput);
268 for (int c = 0; c < channels; ++c) {
269 auto &&outputChannel = outputs[c];
270 auto &&stftChannel = stft[c];
271 for (int i = 0; i < plainOutput; ++i) {
272 // TODO: plain output should be gain-
273 outputChannel[i] = stftChannel[i];
274 }
275 for (int i = 0; i < foldedBackOutput; ++i) {
276 outputChannel[outputSamples - 1 - i] -= stftChannel[plainOutput + i];
277 }
278 for (int i = 0; i < plainOutput + foldedBackOutput; ++i) {
279 stftChannel[i] = 0;
280 }
281 }
282 // Skip the output we just used/cleared
283 stft += plainOutput + foldedBackOutput;
284 // Reset the phase-vocoder stuff, so the next block gets a fresh start
285 for (int c = 0; c < channels; ++c) {
286 auto channelBands = bandsForChannel(channel: c);
287 for (int b = 0; b < bands; ++b) {
288 channelBands[b].prevInput = channelBands[b].prevOutput = 0;
289 }
290 }
291 flushed = true;
292 }
293private:
294 using Complex = std::complex<Sample>;
295 static constexpr Sample noiseFloor{1e-15};
296 static constexpr Sample maxCleanStretch{2}; // time-stretch ratio before we start randomising phases
297 int silenceCounter = 0;
298 bool silenceFirst = true;
299
300 Sample freqMultiplier = 1, freqTonalityLimit = 0.5;
301 std::function<Sample(Sample)> customFreqMap = nullptr;
302
303 signalsmith::spectral::STFT<Sample> stft{0, 1, 1};
304 signalsmith::delay::MultiBuffer<Sample> inputBuffer;
305 int channels = 0, bands = 0;
306 int prevInputOffset = -1;
307 std::vector<Sample> timeBuffer;
308 bool didSeek = false, flushed = true;
309 Sample seekTimeFactor = 1;
310
311 std::vector<Complex> rotCentreSpectrum, rotPrevInterval;
312 Sample bandToFreq(Sample b) const {
313 return (b + Sample(0.5))/stft.fftSize();
314 }
315 Sample freqToBand(Sample f) const {
316 return f*stft.fftSize() - Sample(0.5);
317 }
318 void timeShiftPhases(Sample shiftSamples, std::vector<Complex> &output) const {
319 for (int b = 0; b < bands; ++b) {
320 Sample phase = bandToFreq(b)*shiftSamples*Sample(-2*M_PI);
321 output[b] = {std::cos(phase), std::sin(phase)};
322 }
323 }
324
325 struct Band {
326 Complex input, prevInput{0};
327 Complex output, prevOutput{0};
328 Sample inputEnergy;
329 };
330 std::vector<Band> channelBands;
331 Band * bandsForChannel(int channel) {
332 return channelBands.data() + channel*bands;
333 }
334 template<Complex Band::*member>
335 Complex getBand(int channel, int index) {
336 if (index < 0 || index >= bands) return 0;
337 return channelBands[index + channel*bands].*member;
338 }
339 template<Complex Band::*member>
340 Complex getFractional(int channel, int lowIndex, Sample fractional) {
341 Complex low = getBand<member>(channel, lowIndex);
342 Complex high = getBand<member>(channel, lowIndex + 1);
343 return low + (high - low)*fractional;
344 }
345 template<Complex Band::*member>
346 Complex getFractional(int channel, Sample inputIndex) {
347 int lowIndex = std::floor(inputIndex);
348 Sample fracIndex = inputIndex - lowIndex;
349 return getFractional<member>(channel, lowIndex, fracIndex);
350 }
351 template<Sample Band::*member>
352 Sample getBand(int channel, int index) {
353 if (index < 0 || index >= bands) return 0;
354 return channelBands[index + channel*bands].*member;
355 }
356 template<Sample Band::*member>
357 Sample getFractional(int channel, int lowIndex, Sample fractional) {
358 Sample low = getBand<member>(channel, lowIndex);
359 Sample high = getBand<member>(channel, lowIndex + 1);
360 return low + (high - low)*fractional;
361 }
362 template<Sample Band::*member>
363 Sample getFractional(int channel, Sample inputIndex) {
364 int lowIndex = std::floor(inputIndex);
365 Sample fracIndex = inputIndex - lowIndex;
366 return getFractional<member>(channel, lowIndex, fracIndex);
367 }
368
369 struct Peak {
370 Sample input, output;
371 };
372 std::vector<Peak> peaks;
373 std::vector<Sample> energy, smoothedEnergy;
374 struct PitchMapPoint {
375 Sample inputBin, freqGrad;
376 };
377 std::vector<PitchMapPoint> outputMap;
378
379 struct Prediction {
380 Sample energy = 0;
381 Complex input;
382 Complex shortVerticalTwist, longVerticalTwist;
383
384 Complex makeOutput(Complex phase) {
385 Sample phaseNorm = std::norm(phase);
386 if (phaseNorm <= noiseFloor) {
387 phase = input; // prediction is too weak, fall back to the input
388 phaseNorm = std::norm(input) + noiseFloor;
389 }
390 return phase*std::sqrt(energy/phaseNorm);
391 }
392 };
393 std::vector<Prediction> channelPredictions;
394 Prediction * predictionsForChannel(int c) {
395 return channelPredictions.data() + c*bands;
396 }
397
398 std::default_random_engine randomEngine;
399
400 void processSpectrum(bool newSpectrum, Sample timeFactor) {
401 timeFactor = std::max<Sample>(timeFactor, 1/maxCleanStretch);
402 bool randomTimeFactor = (timeFactor > maxCleanStretch);
403 std::uniform_real_distribution<Sample> timeFactorDist(maxCleanStretch*2*randomTimeFactor - timeFactor, timeFactor);
404
405 if (newSpectrum) {
406 for (int c = 0; c < channels; ++c) {
407 auto bins = bandsForChannel(channel: c);
408 for (int b = 0; b < bands; ++b) {
409 auto &bin = bins[b];
410 bin.prevOutput = signalsmith::perf::mul(bin.prevOutput, rotPrevInterval[b]);
411 bin.prevInput = signalsmith::perf::mul(bin.prevInput, rotPrevInterval[b]);
412 }
413 }
414 }
415
416 Sample smoothingBins = Sample(stft.fftSize())/stft.interval();
417 int longVerticalStep = std::round(smoothingBins);
418 if (customFreqMap || freqMultiplier != 1) {
419 findPeaks(smoothingBins);
420 updateOutputMap();
421 } else { // we're not pitch-shifting, so no need to find peaks etc.
422 for (int c = 0; c < channels; ++c) {
423 Band *bins = bandsForChannel(channel: c);
424 for (int b = 0; b < bands; ++b) {
425 bins[b].inputEnergy = std::norm(bins[b].input);
426 }
427 }
428 for (int b = 0; b < bands; ++b) {
429 outputMap[b] = {Sample(b), 1};
430 }
431 }
432
433 // Preliminary output prediction from phase-vocoder
434 for (int c = 0; c < channels; ++c) {
435 Band *bins = bandsForChannel(channel: c);
436 auto *predictions = predictionsForChannel(c);
437 for (int b = 0; b < bands; ++b) {
438 auto mapPoint = outputMap[b];
439 int lowIndex = std::floor(mapPoint.inputBin);
440 Sample fracIndex = mapPoint.inputBin - lowIndex;
441
442 Prediction &prediction = predictions[b];
443 Sample prevEnergy = prediction.energy;
444 prediction.energy = getFractional<&Band::inputEnergy>(c, lowIndex, fracIndex);
445 prediction.energy *= std::max<Sample>(0, mapPoint.freqGrad); // scale the energy according to local stretch factor
446 prediction.input = getFractional<&Band::input>(c, lowIndex, fracIndex);
447
448 auto &outputBin = bins[b];
449 Complex prevInput = getFractional<&Band::prevInput>(c, lowIndex, fracIndex);
450 Complex freqTwist = signalsmith::perf::mul<true>(prediction.input, prevInput);
451 Complex phase = signalsmith::perf::mul(outputBin.prevOutput, freqTwist);
452 outputBin.output = phase/(std::max(prevEnergy, prediction.energy) + noiseFloor);
453
454 if (b > 0) {
455 Sample binTimeFactor = randomTimeFactor ? timeFactorDist(randomEngine) : timeFactor;
456 Complex downInput = getFractional<&Band::input>(c, mapPoint.inputBin - binTimeFactor);
457 prediction.shortVerticalTwist = signalsmith::perf::mul<true>(prediction.input, downInput);
458 if (b >= longVerticalStep) {
459 Complex longDownInput = getFractional<&Band::input>(c, mapPoint.inputBin - longVerticalStep*binTimeFactor);
460 prediction.longVerticalTwist = signalsmith::perf::mul<true>(prediction.input, longDownInput);
461 } else {
462 prediction.longVerticalTwist = 0;
463 }
464 } else {
465 prediction.shortVerticalTwist = prediction.longVerticalTwist = 0;
466 }
467 }
468 }
469
470 // Re-predict using phase differences between frequencies
471 for (int b = 0; b < bands; ++b) {
472 // Find maximum-energy channel and calculate that
473 int maxChannel = 0;
474 Sample maxEnergy = predictionsForChannel(c: 0)[b].energy;
475 for (int c = 1; c < channels; ++c) {
476 Sample e = predictionsForChannel(c)[b].energy;
477 if (e > maxEnergy) {
478 maxChannel = c;
479 maxEnergy = e;
480 }
481 }
482
483 auto *predictions = predictionsForChannel(c: maxChannel);
484 auto &prediction = predictions[b];
485 auto *bins = bandsForChannel(channel: maxChannel);
486 auto &outputBin = bins[b];
487
488 Complex phase = 0;
489
490 // Upwards vertical steps
491 if (b > 0) {
492 auto &downBin = bins[b - 1];
493 phase += signalsmith::perf::mul(downBin.output, prediction.shortVerticalTwist);
494
495 if (b >= longVerticalStep) {
496 auto &longDownBin = bins[b - longVerticalStep];
497 phase += signalsmith::perf::mul(longDownBin.output, prediction.longVerticalTwist);
498 }
499 }
500 // Downwards vertical steps
501 if (b < bands - 1) {
502 auto &upPrediction = predictions[b + 1];
503 auto &upBin = bins[b + 1];
504 phase += signalsmith::perf::mul<true>(upBin.output, upPrediction.shortVerticalTwist);
505
506 if (b < bands - longVerticalStep) {
507 auto &longUpPrediction = predictions[b + longVerticalStep];
508 auto &longUpBin = bins[b + longVerticalStep];
509 phase += signalsmith::perf::mul<true>(longUpBin.output, longUpPrediction.longVerticalTwist);
510 }
511 }
512
513 outputBin.output = prediction.makeOutput(phase);
514
515 // All other bins are locked in phase
516 for (int c = 0; c < channels; ++c) {
517 if (c != maxChannel) {
518 auto &channelBin = bandsForChannel(channel: c)[b];
519 auto &channelPrediction = predictionsForChannel(c)[b];
520
521 Complex channelTwist = signalsmith::perf::mul<true>(channelPrediction.input, prediction.input);
522 Complex channelPhase = signalsmith::perf::mul(outputBin.output, channelTwist);
523 channelBin.output = channelPrediction.makeOutput(channelPhase);
524 }
525 }
526 }
527
528 if (newSpectrum) {
529 for (auto &bin : channelBands) {
530 bin.prevOutput = bin.output;
531 bin.prevInput = bin.input;
532 }
533 } else {
534 for (auto &bin : channelBands) bin.prevOutput = bin.output;
535 }
536 }
537
538 // Produces smoothed energy across all channels
539 void smoothEnergy(Sample smoothingBins) {
540 Sample smoothingSlew = 1/(1 + smoothingBins*Sample(0.5));
541 for (auto &e : energy) e = 0;
542 for (int c = 0; c < channels; ++c) {
543 Band *bins = bandsForChannel(channel: c);
544 for (int b = 0; b < bands; ++b) {
545 Sample e = std::norm(bins[b].input);
546 bins[b].inputEnergy = e; // Used for interpolating prediction energy
547 energy[b] += e;
548 }
549 }
550 for (int b = 0; b < bands; ++b) {
551 smoothedEnergy[b] = energy[b];
552 }
553 Sample e = 0;
554 for (int repeat = 0; repeat < 2; ++repeat) {
555 for (int b = bands - 1; b >= 0; --b) {
556 e += (smoothedEnergy[b] - e)*smoothingSlew;
557 smoothedEnergy[b] = e;
558 }
559 for (int b = 0; b < bands; ++b) {
560 e += (smoothedEnergy[b] - e)*smoothingSlew;
561 smoothedEnergy[b] = e;
562 }
563 }
564 }
565
566 Sample mapFreq(Sample freq) const {
567 if (customFreqMap) return customFreqMap(freq);
568 if (freq > freqTonalityLimit) {
569 Sample diff = freq - freqTonalityLimit;
570 return freqTonalityLimit*freqMultiplier + diff;
571 }
572 return freq*freqMultiplier;
573 }
574
575 // Identifies spectral peaks using energy across all channels
576 void findPeaks(Sample smoothingBins) {
577 smoothEnergy(smoothingBins);
578
579 peaks.resize(0);
580
581 int start = 0;
582 while (start < bands) {
583 if (energy[start] > smoothedEnergy[start]) {
584 int end = start;
585 Sample bandSum = 0, energySum = 0;
586 while (end < bands && energy[end] > smoothedEnergy[end]) {
587 bandSum += end*energy[end];
588 energySum += energy[end];
589 ++end;
590 }
591 Sample avgBand = bandSum/energySum;
592 Sample avgFreq = bandToFreq(b: avgBand);
593 peaks.emplace_back(Peak{avgBand, freqToBand(f: mapFreq(freq: avgFreq))});
594
595 start = end;
596 }
597 ++start;
598 }
599 }
600
601 void updateOutputMap() {
602 if (peaks.empty()) {
603 for (int b = 0; b < bands; ++b) {
604 outputMap[b] = {Sample(b), 1};
605 }
606 return;
607 }
608 Sample bottomOffset = peaks[0].input - peaks[0].output;
609 for (int b = 0; b < std::min<int>(bands, std::ceil(peaks[0].output)); ++b) {
610 outputMap[b] = {b + bottomOffset, 1};
611 }
612 // Interpolate between points
613 for (size_t p = 1; p < peaks.size(); ++p) {
614 const Peak &prev = peaks[p - 1], &next = peaks[p];
615 Sample rangeScale = 1/(next.output - prev.output);
616 Sample outOffset = prev.input - prev.output;
617 Sample outScale = next.input - next.output - prev.input + prev.output;
618 Sample gradScale = outScale*rangeScale;
619 int startBin = std::max<int>(0, std::ceil(prev.output));
620 int endBin = std::min<int>(bands, std::ceil(next.output));
621 for (int b = startBin; b < endBin; ++b) {
622 Sample r = (b - prev.output)*rangeScale;
623 Sample h = r*r*(3 - 2*r);
624 Sample outB = b + outOffset + h*outScale;
625
626 Sample gradH = 6*r*(1 - r);
627 Sample gradB = 1 + gradH*gradScale;
628
629 outputMap[b] = {outB, gradB};
630 }
631 }
632 Sample topOffset = peaks.back().input - peaks.back().output;
633 for (int b = std::max<int>(0, peaks.back().output); b < bands; ++b) {
634 outputMap[b] = {b + topOffset, 1};
635 }
636 }
637};
638
639}} // namespace
640#endif // include guard
641

source code of qtmultimedia/src/3rdparty/signalsmith-stretch/signalsmith-stretch.h