Add .seek() method to setup/move input position
Also update DSP library to v1.6.0 (for `STFT::nextInvalid()`)
This commit is contained in:
parent
c3153785b0
commit
ebaf93d494
@ -22,7 +22,7 @@ Delay delayLine(1024);
|
||||
You can add a compile-time version-check to make sure you have a compatible version of the library:
|
||||
```cpp
|
||||
#include "dsp/envelopes.h"
|
||||
SIGNALSMITH_DSP_VERSION_CHECK(1, 4, 4)
|
||||
SIGNALSMITH_DSP_VERSION_CHECK(1, 6, 0)
|
||||
```
|
||||
|
||||
### Development / contributing
|
||||
|
||||
@ -14,9 +14,9 @@ namespace signalsmith {
|
||||
*/
|
||||
|
||||
#define SIGNALSMITH_DSP_VERSION_MAJOR 1
|
||||
#define SIGNALSMITH_DSP_VERSION_MINOR 4
|
||||
#define SIGNALSMITH_DSP_VERSION_PATCH 4
|
||||
#define SIGNALSMITH_DSP_VERSION_STRING "1.4.4"
|
||||
#define SIGNALSMITH_DSP_VERSION_MINOR 6
|
||||
#define SIGNALSMITH_DSP_VERSION_PATCH 0
|
||||
#define SIGNALSMITH_DSP_VERSION_STRING "1.6.0"
|
||||
|
||||
/** Version compatability check.
|
||||
\code{.cpp}
|
||||
@ -39,5 +39,5 @@ namespace signalsmith {
|
||||
} // signalsmith::
|
||||
#else
|
||||
// If we've already included it, check it's the same version
|
||||
static_assert(SIGNALSMITH_DSP_VERSION_MAJOR == 1 && SIGNALSMITH_DSP_VERSION_MINOR == 4 && SIGNALSMITH_DSP_VERSION_PATCH == 4, "multiple versions of the Signalsmith DSP library");
|
||||
static_assert(SIGNALSMITH_DSP_VERSION_MAJOR == 1 && SIGNALSMITH_DSP_VERSION_MINOR == 6 && SIGNALSMITH_DSP_VERSION_PATCH == 0, "multiple versions of the Signalsmith DSP library");
|
||||
#endif // include guard
|
||||
|
||||
180
dsp/curves.h
180
dsp/curves.h
@ -91,6 +91,10 @@ namespace curves {
|
||||
Cubic dx() const {
|
||||
return {xStart, a1, 2*a2, 3*a3, 0};
|
||||
}
|
||||
Sample dx(Sample x) const {
|
||||
x -= xStart;
|
||||
return a1 + x*(2*a2 + x*(3*a3));
|
||||
}
|
||||
|
||||
/// Cubic segment based on start/end values and gradients
|
||||
static Cubic hermite(Sample x0, Sample x1, Sample y0, Sample y1, Sample g0, Sample g1) {
|
||||
@ -110,25 +114,28 @@ namespace curves {
|
||||
|
||||
Sample grad1 = gradient(x1, x2, y1, y2);
|
||||
Sample curveGrad1 = grad1;
|
||||
bool chooseGrad1 = false;
|
||||
if (x0 != x1) { // we have a defined x0-x1 gradient
|
||||
Sample grad0 = gradient(x0, x1, y0, y1);
|
||||
curveGrad1 = (grad0 + grad1)*Sample(0.5);
|
||||
if (monotonic) ensureMonotonic(curveGrad1, grad0, grad1);
|
||||
} else if (y0 != y1) {
|
||||
if ((y1 > y0) != (grad1 >= 0)) curveGrad1 = 0; // set to 0 if it's a min/max
|
||||
} else if (y0 != y1 && (y1 > y0) != (grad1 >= 0)) {
|
||||
curveGrad1 = 0; // set to 0 if it's a min/max
|
||||
} else {
|
||||
curveGrad1 = 0;
|
||||
chooseGrad1 = true;
|
||||
}
|
||||
Sample curveGrad2;
|
||||
if (x2 != x3) { // we have a defined x1-x2 gradient
|
||||
Sample grad2 = gradient(x2, x3, y2, y3);
|
||||
curveGrad2 = (grad1 + grad2)*Sample(0.5);
|
||||
if (monotonic) ensureMonotonic(curveGrad2, grad1, grad2);
|
||||
|
||||
if (x0 == x1) { // If the other gradient isn't defined, make one up
|
||||
chooseGradient(curveGrad1, grad1, curveGrad2, y0, y1, monotonic);
|
||||
}
|
||||
} else {
|
||||
chooseGradient(curveGrad2, grad1, curveGrad1, y2, y3, monotonic);
|
||||
}
|
||||
if (chooseGrad1) {
|
||||
chooseGradient(curveGrad1, grad1, curveGrad2, y0, y1, monotonic);
|
||||
}
|
||||
return hermite(x1, x2, y1, y2, curveGrad1, curveGrad2);
|
||||
}
|
||||
};
|
||||
@ -140,6 +147,12 @@ namespace curves {
|
||||
class CubicSegmentCurve {
|
||||
struct Point {
|
||||
Sample x, y;
|
||||
Sample lineGrad = 0, curveGrad = 0;
|
||||
bool hasCurveGrad = false;
|
||||
|
||||
Point() : Point(0, 0) {}
|
||||
Point(Sample x, Sample y) : x(x), y(y) {}
|
||||
|
||||
bool operator <(const Point &other) const {
|
||||
return x < other.x;
|
||||
}
|
||||
@ -148,7 +161,25 @@ namespace curves {
|
||||
Point first{0, 0}, last{0, 0};
|
||||
|
||||
std::vector<Cubic<Sample>> _segments{1};
|
||||
// Not public because it's only valid inside the bounds
|
||||
const Cubic<Sample> & findSegment(Sample x) const {
|
||||
// Binary search
|
||||
size_t low = 0, high = _segments.size();
|
||||
while (true) {
|
||||
size_t mid = (low + high)/2;
|
||||
if (low == mid) break;
|
||||
if (_segments[mid].start() <= x) {
|
||||
low = mid;
|
||||
} else {
|
||||
high = mid;
|
||||
}
|
||||
}
|
||||
return _segments[low];
|
||||
}
|
||||
public:
|
||||
Sample lowGrad = 0;
|
||||
Sample highGrad = 0;
|
||||
|
||||
/// Clear existing points and segments
|
||||
void clear() {
|
||||
points.resize(0);
|
||||
@ -164,52 +195,118 @@ namespace curves {
|
||||
}
|
||||
|
||||
/// Recalculates the segments.
|
||||
void update(bool monotonic=false) {
|
||||
void update(bool monotonic=false, bool extendGrad=true, Sample monotonicFactor=3) {
|
||||
if (points.empty()) add(0, 0);
|
||||
std::stable_sort(points.begin(), points.end()); // Ensure ascending order
|
||||
_segments.resize(0);
|
||||
|
||||
// Calculate the point-to-point gradients
|
||||
for (size_t i = 1; i < points.size(); ++i) {
|
||||
auto &prev = points[i - 1];
|
||||
auto &next = points[i];
|
||||
if (prev.x != next.x) {
|
||||
prev.lineGrad = (next.y - prev.y)/(next.x - prev.x);
|
||||
} else {
|
||||
prev.lineGrad = 0;
|
||||
}
|
||||
}
|
||||
|
||||
for (auto &p : points) p.hasCurveGrad = false;
|
||||
points[0].curveGrad = lowGrad;
|
||||
points[0].hasCurveGrad = true;
|
||||
points.back().curveGrad = highGrad;
|
||||
points.back().hasCurveGrad = true;
|
||||
|
||||
// Calculate curve gradient where we know it
|
||||
for (size_t i = 1; i + 1 < points.size(); ++i) {
|
||||
auto &p0 = points[i - 1];
|
||||
auto &p1 = points[i];
|
||||
auto &p2 = points[i + 1];
|
||||
if (p0.x != p1.x && p1.x != p2.x) {
|
||||
p1.curveGrad = (p0.lineGrad + p1.lineGrad)*Sample(0.5);
|
||||
p1.hasCurveGrad = true;
|
||||
}
|
||||
}
|
||||
|
||||
for (size_t i = 1; i < points.size(); ++i) {
|
||||
Point &p1 = points[i - 1];
|
||||
Point &p2 = points[i];
|
||||
if (p1.x == p2.x) continue;
|
||||
if (p1.hasCurveGrad) {
|
||||
if (!p2.hasCurveGrad) {
|
||||
p2.curveGrad = 2*p1.lineGrad - p1.curveGrad;
|
||||
}
|
||||
} else if (p2.hasCurveGrad) {
|
||||
p1.curveGrad = 2*p1.lineGrad - p2.curveGrad;
|
||||
} else {
|
||||
p1.curveGrad = p2.curveGrad = p1.lineGrad;
|
||||
}
|
||||
}
|
||||
|
||||
if (monotonic) {
|
||||
for (size_t i = 1; i < points.size(); ++i) {
|
||||
Point &p1 = points[i - 1];
|
||||
Point &p2 = points[i];
|
||||
if (p1.x != p2.x) {
|
||||
if (p1.lineGrad >= 0) {
|
||||
p1.curveGrad = std::max<Sample>(0, std::min(p1.curveGrad, p1.lineGrad*monotonicFactor));
|
||||
p2.curveGrad = std::max<Sample>(0, std::min(p2.curveGrad, p1.lineGrad*monotonicFactor));
|
||||
} else {
|
||||
p1.curveGrad = std::min<Sample>(0, std::max(p1.curveGrad, p1.lineGrad*monotonicFactor));
|
||||
p2.curveGrad = std::min<Sample>(0, std::max(p2.curveGrad, p1.lineGrad*monotonicFactor));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (size_t i = 1; i < points.size(); ++i) {
|
||||
Point &p1 = points[i - 1];
|
||||
Point &p2 = points[i];
|
||||
if (p1.x != p2.x) {
|
||||
_segments.push_back(Segment::hermite(p1.x, p2.x, p1.y, p2.y, p1.curveGrad, p2.curveGrad));
|
||||
}
|
||||
}
|
||||
|
||||
first = points[0];
|
||||
last = points.back();
|
||||
for (size_t i = 1; i < points.size(); ++i) {
|
||||
Point p1 = points[i - 1];
|
||||
Point p2 = points[i];
|
||||
if (p1.x != p2.x) {
|
||||
Point p0 = (i > 1) ? points[i - 2] : Point{p1.x, p2.y};
|
||||
Point p3 = (i + 1 < points.size()) ? points[i + 1] : Point{p2.x, p1.y};
|
||||
_segments.push_back(Segment::smooth(p0.x, p1.x, p2.x, p3.x, p0.y, p1.y, p2.y, p3.y, monotonic));
|
||||
if (extendGrad && _segments.size()) {
|
||||
if (points[0].x != points[1].x || points[0].y == points[1].y) {
|
||||
lowGrad = _segments[0].dx(first.x);
|
||||
}
|
||||
auto &last = points.back(), &last2 = points[points.size() - 1];
|
||||
if (last.x != last2.x || last.y == last2.y) {
|
||||
highGrad = _segments.back().dx(last.x);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/// Reads a value out from the curve.
|
||||
Sample operator()(Sample x) const {
|
||||
if (x <= first.x) return first.y;
|
||||
if (x >= last.x) return last.y;
|
||||
|
||||
// Binary search
|
||||
size_t low = 0, high = _segments.size();
|
||||
while (true) {
|
||||
size_t mid = (low + high)/2;
|
||||
if (low == mid) break;
|
||||
if (_segments[mid].start() <= x) {
|
||||
low = mid;
|
||||
} else {
|
||||
high = mid;
|
||||
}
|
||||
}
|
||||
return _segments[low](x);
|
||||
if (x <= first.x) return first.y + (x - first.x)*lowGrad;
|
||||
if (x >= last.x) return last.y + (x - last.x)*highGrad;
|
||||
return findSegment(x)(x);
|
||||
}
|
||||
|
||||
|
||||
CubicSegmentCurve dx() const {
|
||||
CubicSegmentCurve result{*this};
|
||||
result.first.y = result.last.y = 0;
|
||||
result.first.y = lowGrad;
|
||||
result.last.y = highGrad;
|
||||
result.lowGrad = result.highGrad = 0;
|
||||
for (auto &s : result._segments) {
|
||||
s = s.dx();
|
||||
}
|
||||
return result;
|
||||
}
|
||||
Sample dx(Sample x) const {
|
||||
if (x < first.x) return lowGrad;
|
||||
if (x >= last.x) return highGrad;
|
||||
return findSegment(x).dx(x);
|
||||
}
|
||||
|
||||
using Segment = Cubic<Sample>;
|
||||
std::vector<Segment> & segments() {
|
||||
return _segments;
|
||||
}
|
||||
const std::vector<Segment> & segments() const {
|
||||
return _segments;
|
||||
}
|
||||
@ -222,6 +319,21 @@ namespace curves {
|
||||
Sample a, b, c, d; // (a + bx)/(c + dx)
|
||||
Reciprocal(Sample a, Sample b, Sample c, Sample d) : a(a), b(b), c(c), d(d) {}
|
||||
public:
|
||||
/** Decent approximation to the Bark scale
|
||||
|
||||
The Bark index goes from 1-24, but this map is valid from approximately 0.25 - 27.5.
|
||||
You can get the bandwidth by `barkScale.dx(barkIndex)`.
|
||||
\diagram{curves-reciprocal-approx-bark.svg}*/
|
||||
static Reciprocal<Sample> barkScale() {
|
||||
return {1, 10, 24, 60, 1170, 13500};
|
||||
}
|
||||
/// Returns a map from 0-1 to the given (non-negative) Hz range.
|
||||
static Reciprocal<Sample> barkRange(Sample lowHz, Sample highHz) {
|
||||
Reciprocal bark = barkScale();
|
||||
Sample lowBark = bark.inverse(lowHz), highBark = bark.inverse(highHz);
|
||||
return Reciprocal(lowBark, (lowBark + highBark)/2, highBark).then(bark);
|
||||
}
|
||||
|
||||
Reciprocal() : Reciprocal(0, 0.5, 1) {}
|
||||
/// If no x-range given, default to the unit range
|
||||
Reciprocal(Sample y0, Sample y1, Sample y2) : Reciprocal(0, 0.5, 1, y0, y1, y2) {}
|
||||
@ -243,13 +355,17 @@ namespace curves {
|
||||
Sample inverse(Sample y) const {
|
||||
return (c*y - a)/(b - d*y);
|
||||
}
|
||||
Sample dx(Sample x) const {
|
||||
Sample l = (c + d*x);
|
||||
return (b*c - a*d)/(l*l);
|
||||
}
|
||||
|
||||
/// Combine two `Reciprocal`s together in sequence
|
||||
Reciprocal then(const Reciprocal &other) const {
|
||||
return Reciprocal(other.a*c + other.b*a, other.a*d + other.b*b, other.c*c + other.d*a, other.c*d + other.d*b);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
/** @} */
|
||||
}} // namespace
|
||||
#endif // include guard
|
||||
|
||||
@ -518,77 +518,6 @@ namespace envelopes {
|
||||
}
|
||||
};
|
||||
|
||||
/** Rolling Statistics
|
||||
\diagram{envelopes-order-statistics}
|
||||
|
||||
This keeps a history of recent samples
|
||||
*/
|
||||
template<typename Sample>
|
||||
class RollingStats {
|
||||
size_t index = 0;
|
||||
std::vector<Sample> buffer, sorted;
|
||||
|
||||
struct SkipListEntry {
|
||||
Sample value;
|
||||
int distance;
|
||||
int nextEntry;
|
||||
int finerEntry;
|
||||
|
||||
int findBefore(const Sample &v, int thisIndex) const {
|
||||
if (nextEntry == -1) return thisIndex;
|
||||
SkipListEntry &next = skipList[nextEntry];
|
||||
if (
|
||||
}
|
||||
};
|
||||
std::vector<SkipListEntry> skipList;
|
||||
const int findBefore(const Sample &v) const {
|
||||
if (skipList.empty() || skipList[0].value >= v) return -1;
|
||||
|
||||
return skipList[0].findBefore(v, 0);
|
||||
}
|
||||
void addEntry(const Sample &v) {
|
||||
if (skipList.empty()) {
|
||||
skipList.push_back({v, 1, -1, -1});
|
||||
return;
|
||||
}
|
||||
SkipListEntry *entry = &skipList[0];
|
||||
while (entry) {
|
||||
if (
|
||||
}
|
||||
}
|
||||
|
||||
public:
|
||||
RollingStats(int length, Sample initial=Sample()) {
|
||||
buffer.assign(length, initial);
|
||||
sorted.assign(length, initial);
|
||||
|
||||
for (int i = 0; i < length; ++i) {
|
||||
addEntry(initial);
|
||||
}
|
||||
}
|
||||
|
||||
void reset(Sample initial=Sample()) {
|
||||
buffer.assign(buffer.size(), initial);
|
||||
sorted.assign(buffer.size(), initial);
|
||||
index = 0;
|
||||
}
|
||||
|
||||
void operator()(const Sample &v) {
|
||||
buffer[index] = v;
|
||||
sorted = buffer;
|
||||
std::sort(sorted.begin(), sorted.end());
|
||||
if (++index >= buffer.size()) index = 0;
|
||||
}
|
||||
|
||||
const Sample & operator[](int index) const {
|
||||
return sorted[index];
|
||||
}
|
||||
|
||||
Sample percentile(double percentile) {
|
||||
return sorted[(sorted.size() - 1e-10)*percentile];
|
||||
}
|
||||
};
|
||||
|
||||
/** @} */
|
||||
}} // signalsmith::envelopes::
|
||||
#endif // include guard
|
||||
|
||||
124
dsp/spectral.h
124
dsp/spectral.h
@ -33,7 +33,7 @@ namespace spectral {
|
||||
|
||||
std::vector<Sample> fftWindow;
|
||||
std::vector<Sample> timeBuffer;
|
||||
std::vector<Complex> freqBuffer;
|
||||
int offsetSamples = 0;
|
||||
public:
|
||||
/// Returns a fast FFT size <= `size`
|
||||
static int fastSizeAbove(int size, int divisor=1) {
|
||||
@ -45,26 +45,27 @@ namespace spectral {
|
||||
}
|
||||
|
||||
WindowedFFT() {}
|
||||
WindowedFFT(int size) {
|
||||
setSize(size);
|
||||
WindowedFFT(int size, int rotateSamples=0) {
|
||||
setSize(size, rotateSamples);
|
||||
}
|
||||
template<class WindowFn>
|
||||
WindowedFFT(int size, WindowFn fn, Sample windowOffset=0.5) {
|
||||
setSize(size, fn, windowOffset);
|
||||
WindowedFFT(int size, WindowFn fn, Sample windowOffset=0.5, int rotateSamples=0) {
|
||||
setSize(size, fn, windowOffset, rotateSamples);
|
||||
}
|
||||
|
||||
/// Sets the size, returning the window for modification (initially all 1s)
|
||||
std::vector<Sample> & setSizeWindow(int size) {
|
||||
std::vector<Sample> & setSizeWindow(int size, int rotateSamples=0) {
|
||||
mrfft.setSize(size);
|
||||
fftWindow.resize(size, 1);
|
||||
fftWindow.assign(size, 1);
|
||||
timeBuffer.resize(size);
|
||||
freqBuffer.resize(size);
|
||||
offsetSamples = rotateSamples;
|
||||
if (offsetSamples < 0) offsetSamples += size; // TODO: for a negative rotation, the other half of the result is inverted
|
||||
return fftWindow;
|
||||
}
|
||||
/// Sets the FFT size, with a user-defined functor for the window
|
||||
template<class WindowFn>
|
||||
void setSize(int size, WindowFn fn, Sample windowOffset=0.5) {
|
||||
setSizeWindow(size);
|
||||
void setSize(int size, WindowFn fn, Sample windowOffset=0.5, int rotateSamples=0) {
|
||||
setSizeWindow(size, rotateSamples);
|
||||
|
||||
Sample invSize = 1/(Sample)size;
|
||||
for (int i = 0; i < size; ++i) {
|
||||
@ -73,12 +74,12 @@ namespace spectral {
|
||||
}
|
||||
}
|
||||
/// Sets the size (using the default Blackman-Harris window)
|
||||
void setSize(int size) {
|
||||
void setSize(int size, int rotateSamples=0) {
|
||||
setSize(size, [](double x) {
|
||||
double phase = 2*M_PI*x;
|
||||
// Blackman-Harris
|
||||
return 0.35875 + 0.48829*std::cos(phase) + 0.14128*std::cos(phase*2) + 0.1168*std::cos(phase*3);
|
||||
});
|
||||
return 0.35875 - 0.48829*std::cos(phase) + 0.14128*std::cos(phase*2) - 0.01168*std::cos(phase*3);
|
||||
}, Sample(0.5), rotateSamples);
|
||||
}
|
||||
|
||||
const std::vector<Sample> & window() const {
|
||||
@ -91,17 +92,17 @@ namespace spectral {
|
||||
/// Performs an FFT (with windowing)
|
||||
template<class Input, class Output>
|
||||
void fft(Input &&input, Output &&output) {
|
||||
struct WindowedInput {
|
||||
const Input &input;
|
||||
std::vector<Sample> &window;
|
||||
SIGNALSMITH_INLINE Sample operator [](int i) {
|
||||
return input[i]*window[i];
|
||||
}
|
||||
};
|
||||
|
||||
mrfft.fft(WindowedInput{input, fftWindow}, output);
|
||||
int fftSize = size();
|
||||
for (int i = 0; i < offsetSamples; ++i) {
|
||||
// Inverted polarity since we're using the MRFFT
|
||||
timeBuffer[i + fftSize - offsetSamples] = -input[i]*fftWindow[i];
|
||||
}
|
||||
for (int i = offsetSamples; i < fftSize; ++i) {
|
||||
timeBuffer[i - offsetSamples] = input[i]*fftWindow[i];
|
||||
}
|
||||
mrfft.fft(timeBuffer, output);
|
||||
}
|
||||
/// Performs an FFT (no windowing)
|
||||
/// Performs an FFT (no windowing or rotation)
|
||||
template<class Input, class Output>
|
||||
void fftRaw(Input &&input, Output &&output) {
|
||||
mrfft.fft(input, output);
|
||||
@ -110,14 +111,19 @@ namespace spectral {
|
||||
/// Inverse FFT, with windowing and 1/N scaling
|
||||
template<class Input, class Output>
|
||||
void ifft(Input &&input, Output &&output) {
|
||||
mrfft.ifft(input, output);
|
||||
int size = mrfft.size();
|
||||
Sample norm = 1/(Sample)size;
|
||||
for (int i = 0; i < size; ++i) {
|
||||
output[i] *= norm*fftWindow[i];
|
||||
mrfft.ifft(input, timeBuffer);
|
||||
int fftSize = mrfft.size();
|
||||
Sample norm = 1/(Sample)fftSize;
|
||||
|
||||
for (int i = 0; i < offsetSamples; ++i) {
|
||||
// Inverted polarity since we're using the MRFFT
|
||||
output[i] = -timeBuffer[i + fftSize - offsetSamples]*norm*fftWindow[i];
|
||||
}
|
||||
for (int i = offsetSamples; i < fftSize; ++i) {
|
||||
output[i] = timeBuffer[i - offsetSamples]*norm*fftWindow[i];
|
||||
}
|
||||
}
|
||||
/// Performs an IFFT (no windowing)
|
||||
/// Performs an IFFT (no windowing or rotation)
|
||||
template<class Input, class Output>
|
||||
void ifftRaw(Input &&input, Output &&output) {
|
||||
mrfft.ifft(input, output);
|
||||
@ -213,37 +219,48 @@ namespace spectral {
|
||||
this->_fftSize = fftSize;
|
||||
this->_interval = newInterval;
|
||||
validUntilIndex = -1;
|
||||
|
||||
auto &window = fft.setSizeWindow(fftSize);
|
||||
if (windowShape == Window::kaiser) {
|
||||
using Kaiser = ::signalsmith::windows::Kaiser;
|
||||
/// Roughly optimal Kaiser for STFT analysis (forced to perfect reconstruction)
|
||||
auto kaiser = Kaiser::withBandwidth(windowSize/double(_interval), true);
|
||||
kaiser.fill(window, windowSize);
|
||||
} else {
|
||||
using Confined = ::signalsmith::windows::ApproximateConfinedGaussian;
|
||||
auto confined = Confined::withBandwidth(windowSize/double(_interval));
|
||||
confined.fill(window, windowSize);
|
||||
}
|
||||
::signalsmith::windows::forcePerfectReconstruction(window, windowSize, _interval);
|
||||
|
||||
// TODO: fill extra bits of an input buffer with NaN/Infinity, to break this, and then fix by adding zero-padding to WindowedFFT (as opposed to zero-valued window sections)
|
||||
for (int i = windowSize; i < fftSize; ++i) {
|
||||
window[i] = 0;
|
||||
}
|
||||
setWindow(windowShape);
|
||||
|
||||
spectrum.resize(channels, fftSize/2);
|
||||
timeBuffer.resize(fftSize);
|
||||
}
|
||||
public:
|
||||
enum class Window {kaiser, acg};
|
||||
/// \deprecated use `.setWindow()` which actually updates the window when you change it
|
||||
Window windowShape = Window::kaiser;
|
||||
// for convenience
|
||||
static constexpr Window kaiser = Window::kaiser;
|
||||
static constexpr Window acg = Window::acg;
|
||||
|
||||
/** Swaps between the default (Kaiser) shape and Approximate Confined Gaussian (ACG).
|
||||
\diagram{stft-windows.svg,Default (Kaiser) windows and partial cumulative sum}
|
||||
The ACG has better rolloff since its edges go to 0:
|
||||
\diagram{stft-windows-acg.svg,ACG windows and partial cumulative sum}
|
||||
However, it generally has worse performance in terms of total sidelobe energy, affecting worst-case aliasing levels for (most) higher overlap ratios:
|
||||
\diagram{stft-aliasing-simulated-acg.svg,Simulated bad-case aliasing for ACG windows - compare with above}*/
|
||||
enum class Window {kaiser, acg};
|
||||
Window windowShape = Window::kaiser;
|
||||
// TODO: these should both be set before resize()
|
||||
void setWindow(Window shape, bool rotateToZero=false) {
|
||||
windowShape = shape;
|
||||
|
||||
auto &window = fft.setSizeWindow(_fftSize, rotateToZero ? _windowSize/2 : 0);
|
||||
if (windowShape == Window::kaiser) {
|
||||
using Kaiser = ::signalsmith::windows::Kaiser;
|
||||
/// Roughly optimal Kaiser for STFT analysis (forced to perfect reconstruction)
|
||||
auto kaiser = Kaiser::withBandwidth(_windowSize/double(_interval), true);
|
||||
kaiser.fill(window, _windowSize);
|
||||
} else {
|
||||
using Confined = ::signalsmith::windows::ApproximateConfinedGaussian;
|
||||
auto confined = Confined::withBandwidth(_windowSize/double(_interval));
|
||||
confined.fill(window, _windowSize);
|
||||
}
|
||||
::signalsmith::windows::forcePerfectReconstruction(window, _windowSize, _interval);
|
||||
|
||||
// TODO: fill extra bits of an input buffer with NaN/Infinity, to break this, and then fix by adding zero-padding to WindowedFFT (as opposed to zero-valued window sections)
|
||||
for (int i = _windowSize; i < _fftSize; ++i) {
|
||||
window[i] = 0;
|
||||
}
|
||||
}
|
||||
|
||||
using Spectrum = MultiSpectrum;
|
||||
Spectrum spectrum;
|
||||
@ -274,10 +291,11 @@ namespace spectral {
|
||||
return fft.window();
|
||||
}
|
||||
/// Calculates the effective window for the partially-summed future output (relative to the most recent block)
|
||||
std::vector<Sample> partialSumWindow() const {
|
||||
std::vector<Sample> partialSumWindow(bool includeLatestBlock=true) const {
|
||||
const auto &w = window();
|
||||
std::vector<Sample> result(_windowSize, 0);
|
||||
for (int offset = 0; offset < _windowSize; offset += _interval) {
|
||||
int firstOffset = (includeLatestBlock ? 0 : _interval);
|
||||
for (int offset = firstOffset; offset < _windowSize; offset += _interval) {
|
||||
for (int i = 0; i < _windowSize - offset; ++i) {
|
||||
Sample value = w[i + offset];
|
||||
result[i] += value*value;
|
||||
@ -328,6 +346,10 @@ namespace spectral {
|
||||
void ensureValid(AnalysisFn fn) {
|
||||
return ensureValid(0, fn);
|
||||
}
|
||||
/// Returns the next invalid index (a.k.a. the index of the next block)
|
||||
int nextInvalid() const {
|
||||
return validUntilIndex + 1;
|
||||
}
|
||||
|
||||
/** Analyse a multi-channel input, for any type where `data[channel][index]` returns samples
|
||||
|
||||
@ -342,7 +364,7 @@ namespace spectral {
|
||||
void analyse(int c, Data &&data) {
|
||||
fft.fft(data, spectrum[c]);
|
||||
}
|
||||
/// Analyse without windowing
|
||||
/// Analyse without windowing or zero-rotation
|
||||
template<class Data>
|
||||
void analyseRaw(Data &&data) {
|
||||
for (int c = 0; c < channels; ++c) {
|
||||
|
||||
@ -4,7 +4,7 @@
|
||||
#include "dsp/spectral.h"
|
||||
#include "dsp/delay.h"
|
||||
#include "dsp/perf.h"
|
||||
SIGNALSMITH_DSP_VERSION_CHECK(1, 3, 3); // Check version is compatible
|
||||
SIGNALSMITH_DSP_VERSION_CHECK(1, 6, 0); // Check version is compatible
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
#include <functional>
|
||||
@ -37,6 +37,7 @@ struct SignalsmithStretch {
|
||||
prevInputOffset = -1;
|
||||
channelBands.assign(channelBands.size(), Band());
|
||||
silenceCounter = 2*stft.windowSize();
|
||||
didSeek = false;
|
||||
}
|
||||
|
||||
// Configures using a default preset
|
||||
@ -86,7 +87,24 @@ struct SignalsmithStretch {
|
||||
void setFreqMap(std::function<Sample(Sample)> inputToOutput) {
|
||||
customFreqMap = inputToOutput;
|
||||
}
|
||||
|
||||
|
||||
// Provide previous input ("pre-roll"), without affecting the speed calculation. You should ideally feed it one block-length + one interval
|
||||
template<class Inputs>
|
||||
void seek(Inputs &&inputs, int inputSamples, double playbackRate) {
|
||||
inputBuffer.reset();
|
||||
for (int c = 0; c < channels; ++c) {
|
||||
auto &&inputChannel = inputs[c];
|
||||
auto &&bufferChannel = inputBuffer[c];
|
||||
int startIndex = std::max<int>(0, inputSamples - stft.windowSize() - stft.interval());
|
||||
for (int i = startIndex; i < inputSamples; ++i) {
|
||||
bufferChannel[i] = inputChannel[i];
|
||||
}
|
||||
}
|
||||
inputBuffer += inputSamples;
|
||||
didSeek = true;
|
||||
seekTimeFactor = (playbackRate*stft.interval() > 1) ? 1/playbackRate : stft.interval();
|
||||
}
|
||||
|
||||
template<class Inputs, class Outputs>
|
||||
void process(Inputs &&inputs, int inputSamples, Outputs &&outputs, int outputSamples) {
|
||||
Sample totalEnergy = 0;
|
||||
@ -128,7 +146,7 @@ struct SignalsmithStretch {
|
||||
for (int c = 0; c < channels; ++c) {
|
||||
auto &&inputChannel = inputs[c];
|
||||
auto &&bufferChannel = inputBuffer[c];
|
||||
int startIndex = std::max<int>(0, inputSamples - stft.windowSize());
|
||||
int startIndex = std::max<int>(0, inputSamples - stft.windowSize() - stft.interval());
|
||||
for (int i = startIndex; i < inputSamples; ++i) {
|
||||
bufferChannel[i] = inputChannel[i];
|
||||
}
|
||||
@ -150,7 +168,7 @@ struct SignalsmithStretch {
|
||||
int inputInterval = inputOffset - prevInputOffset;
|
||||
prevInputOffset = inputOffset;
|
||||
|
||||
bool newSpectrum = (inputInterval > 0);
|
||||
bool newSpectrum = didSeek || (inputInterval > 0);
|
||||
if (newSpectrum) {
|
||||
for (int c = 0; c < channels; ++c) {
|
||||
// Copy from the history buffer, if needed
|
||||
@ -174,7 +192,7 @@ struct SignalsmithStretch {
|
||||
}
|
||||
}
|
||||
|
||||
if (inputInterval != stft.interval()) { // make sure the previous input is the correct distance in the past
|
||||
if (didSeek || inputInterval != stft.interval()) { // make sure the previous input is the correct distance in the past
|
||||
int prevIntervalOffset = inputOffset - stft.interval();
|
||||
for (int c = 0; c < channels; ++c) {
|
||||
// Copy from the history buffer, if needed
|
||||
@ -199,8 +217,9 @@ struct SignalsmithStretch {
|
||||
}
|
||||
}
|
||||
|
||||
Sample timeFactor = stft.interval()/std::max<Sample>(1, inputInterval);
|
||||
Sample timeFactor = didSeek ? seekTimeFactor : stft.interval()/std::max<Sample>(1, inputInterval);
|
||||
processSpectrum(newSpectrum, timeFactor);
|
||||
didSeek = false;
|
||||
|
||||
for (int c = 0; c < channels; ++c) {
|
||||
auto channelBands = bandsForChannel(c);
|
||||
@ -247,6 +266,8 @@ private:
|
||||
int channels = 0, bands = 0;
|
||||
int prevInputOffset = -1;
|
||||
std::vector<Sample> timeBuffer;
|
||||
bool didSeek = false;
|
||||
Sample seekTimeFactor = 1;
|
||||
|
||||
std::vector<Complex> rotCentreSpectrum, rotPrevInterval;
|
||||
Sample bandToFreq(Sample b) const {
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user