Skip to content

Commit

Permalink
Fix cubic-curve edges
Browse files Browse the repository at this point in the history
  • Loading branch information
@geraintluff
geraintluff committed Feb 24, 2022
1 parent 4f497b8 commit cb94481
Showing 1 changed file with 142 additions and 85 deletions.
227 changes: 142 additions & 85 deletions envelopes.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,12 @@
#include <random>
#include <vector>
#include <iterator>
#include <algorithm> // std::stable_sort

namespace signalsmith {
namespace envelopes {
/** @defgroup Envelopes Envelopes and LFOs
@brief Envelopes, mostly fast and approximate
@brief LFOs, envelopes and filters for manipulating them
@{
@file
Expand Down Expand Up @@ -176,6 +177,7 @@ namespace envelopes {
};

/** Rectangular moving average filter (FIR).
\diagram{box-filter-example.svg}
A filter of length 1 has order 0 (i.e. does nothing). */
template<typename Sample=double>
class BoxFilter {
Expand Down Expand Up @@ -516,103 +518,158 @@ namespace envelopes {
return value = std::max<Sample>(v, value + (v - peak)*stepMultiplier);
}
};

/* Smooth interpolation (optionally monotonic) between points, using cubic segments.
To produce a sharp corner, use a repeated point. */

/// A real-valued cubic curve
template<typename Sample=double>
class CubicSegments {
// Gradient - must have y0 != y1
class CubicSegment {
Sample xStart, a0, a1, a2, a3;

// Only use with y0 != y1
static inline double gradient(Sample x0, Sample x1, Sample y0, Sample y1) {
return (y1 - y0)/(x1 - x0);
}
public:
class Segment {
Sample xOffset, a, b, c, d;
public:
Segment(Sample xOffset, Sample a, Sample b, Sample c, Sample d) : xOffset(xOffset), a(a), b(b), c(c), d(d) {}

Sample operator ()(Sample x) const {
x -= xOffset;
return a + x*(b + x*(c + d*x));
}
/// Differentiate
Segment dx() const {
return Segment(xOffset, b, 2*c, 3*d, 0);
// Ensure a gradient produces monotonic segments
static inline void ensureMonotonic(Sample &curveGrad, Sample gradA, Sample gradB) {
if ((gradA <= 0 && gradB >= 0) || (gradA >= 0 && gradB <= 0)) {
curveGrad = 0; // point is a local minimum/maximum
} else {
if (std::abs(curveGrad) > std::abs(gradA*3)) {
curveGrad = gradA*3;
}
if (std::abs(curveGrad) > std::abs(gradB*3)) {
curveGrad = gradB*3;
}
}

static Segment hermite(Sample x0, Sample x1, Sample y0, Sample y1, Sample g0, Sample g1) {
Sample xScale = 1/(x1 - x0);
return Segment{
x0, y0, g0,
(3*(y1 - y0)*xScale - 2*g0 - g1)*xScale,
(2*(y0 - y1)*xScale + g0 + g1)*(xScale*xScale)
};
}
// When we have duplicate x-values (either side) make up a gradient
static inline void chooseGradient(Sample &curveGrad, Sample grad1, Sample curveGradOther, Sample y0, Sample y1, bool monotonic) {
curveGrad = 2*grad1 - curveGradOther;
if (y0 != y1 && (y1 > y0) != (grad1 >= 0)) { // not duplicate y, but a local min/max
curveGrad = 0;
} else if (monotonic) {
if (grad1 >= 0) {
curveGrad = std::max<Sample>(0, curveGrad);
} else {
curveGrad = std::min<Sample>(0, curveGrad);
}
}
};
}
public:
CubicSegment() : CubicSegment(0, 0, 0, 0, 0) {}
CubicSegment(Sample xStart, Sample a0, Sample a1, Sample a2, Sample a3) : xStart(xStart), a0(a0), a1(a1), a2(a2), a3(a3) {}

/** Cubic segment valid between `x1` and `x2`, which is smooth when applied to an adjacent set of points.
If `x0 == x1` or `x2 == x3` it will choose a gradient which minimises the second differential.
Sample operator ()(Sample x) const {
x -= xStart;
return a0 + x*(a1 + x*(a2 + x*a3));
}
/// Differentiate
CubicSegment dx() const {
return {xStart, a1, 2*a2, 3*a3, 0};
}
/// The reference x-value, used as the centre of the cubic expansion
Sample start() {
return xStart;
}

/// Cubic segment based on start/end values and gradients
static CubicSegment hermite(Sample x0, Sample x1, Sample y0, Sample y1, Sample g0, Sample g1) {
Sample xScale = 1/(x1 - x0);
return {
x0, y0, g0,
(3*(y1 - y0)*xScale - 2*g0 - g1)*xScale,
(2*(y0 - y1)*xScale + g0 + g1)*(xScale*xScale)
};
}

/** Cubic segment (valid between `x1` and `x2`), which is smooth when applied to an adjacent set of points.
If `x0 == x1` or `x2 == x3` it will choose a gradient which continues in a quadratic curve, or 0 if the point is a local minimum/maximum.
*/
static Segment fromPoints(Sample x0, Sample x1, Sample x2, Sample x3, Sample y0, Sample y1, Sample y2, Sample y3, bool monotonic=false) {
if (x1 == x2) {
return Segment(0, y1, 0, 0, 0);
} else if (x0 == x1) {
Sample grad1 = gradient(x1, x2, y1, y2);
if (x2 == x3) {
// Duplicate points on both sides - return linear segment
return Segment(x1, y1, grad1, 0, 0);
}
Sample grad2 = gradient(x2, x3, y2, y3);
Sample curveGrad2 = (grad1 + grad2)*Sample(0.5);
if (monotonic) {
if ((grad1 <= 0 && grad2 >= 0) || (grad1 >= 0 && grad2 <= 0)) {
curveGrad2 = 0; // Point 2 is a local minimum/maximum
} else if (std::abs(curveGrad2) > std::abs(grad1)*3) {
curveGrad2 = grad1*Sample(1.0/3);
}
}
Sample curveGrad1 = 2*grad1 - curveGrad2;
return Segment::hermite(x1, x2, y1, y2, curveGrad1, curveGrad2);
} else if (x2 == x3) {
static CubicSegment smooth(Sample x0, Sample x1, Sample x2, Sample x3, Sample y0, Sample y1, Sample y2, Sample y3, bool monotonic=false) {
if (x1 == x2) return {0, y1, 0, 0, 0}; // zero-width segment, just return constant

Sample grad1 = gradient(x1, x2, y1, y2);
Sample curveGrad1 = grad1;
if (x0 != x1) { // we have a defined x0-x1 gradient
Sample grad0 = gradient(x0, x1, y0, y1);
Sample grad1 = gradient(x1, x2, y1, y2);
Sample curveGrad1 = (grad0 + grad1)*Sample(0.5);
if (monotonic) {
if ((grad0 <= 0 && grad1 >= 0) || (grad0 >= 0 && grad1 <= 0)) {
curveGrad1 = 0; // Point 1 is a local minimum/maximum
} else if (std::abs(curveGrad1) > std::abs(grad0)*3) {
curveGrad1 = grad1*Sample(1.0/3);
}
}
Sample curveGrad2 = 2*grad1 - curveGrad1;
return Segment::hermite(x1, x2, y1, y2, curveGrad1, curveGrad2);
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
}
Sample grad0 = gradient(x0, x1, y0, y1);
Sample grad1 = gradient(x1, x2, y1, y2);
Sample grad2 = gradient(x2, x3, y2, y3);
Sample curveGrad1 = (grad0 + grad1)*Sample(0.5);
Sample curveGrad2 = (grad1 + grad2)*Sample(0.5);
if (monotonic) {
if ((grad0 <= 0 && grad1 >= 0) || (grad0 >= 0 && grad1 <= 0)) {
curveGrad1 = 0; // Point 1 is a local minimum/maximum
}
if (std::abs(curveGrad1) > std::abs(grad0)*3) {
curveGrad1 = grad0*Sample(1.0/3);
}
if (std::abs(curveGrad1) > std::abs(grad1)*3) {
curveGrad1 = grad1*Sample(1.0/3);
}
if ((grad1 <= 0 && grad2 >= 0) || (grad1 >= 0 && grad2 <= 0)) {
curveGrad2 = 0; // Point 2 is a local minimum/maximum
}
if (std::abs(curveGrad2) > std::abs(grad1)*3) {
curveGrad2 = grad1*Sample(1.0/3);
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);
}
if (std::abs(curveGrad2) > std::abs(grad2)*3) {
curveGrad2 = grad2*Sample(1.0/3);
} else {
chooseGradient(curveGrad2, grad1, curveGrad1, y2, y3, monotonic);
}
return hermite(x1, x2, y1, y2, curveGrad1, curveGrad2);
}
};

/** Smooth interpolation (optionally monotonic) between points, using cubic segments.
\diagram{cubic-segments-example.svg,Example curve including a repeated point at (2.5\, 3.5) and an instantaneous jump at x=4. The curve is flat beyond the first/last points.}
To produce a sharp corner, use a repeated point. The gradient is flat at the edges, unless you use repeated points at the start/end.*/
template<typename Sample=double>
class CubicSegmentCurve {
struct Point {
Sample x, y;
bool operator <(const Point &other) const {
return x < other.x;
}
};
std::vector<Point> points;
Point first{0, 0}, last{0, 0};

using Segment = CubicSegment<Sample>;
std::vector<Segment> segments{1};
public:
/// Clear existing points and segments
void clear() {
points.resize(0);
segments.resize(0);
first = last = {0, 0};
}

/// Add a new point, but does not recalculate the segments. `corner` just writes the point twice, for convenience.
CubicSegmentCurve & add(Sample x, Sample y, bool corner=false) {
points.push_back({x, y});
if (corner) points.push_back({x, y});
return *this;
}

/// Recalculates the segments.
void update(bool monotonic=false) {
if (points.empty()) add(0, 0);
std::stable_sort(points.begin(), points.end()); // Ensure ascending order
segments.resize(0);
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));
}
}
return Segment::hermite(x1, x2, y1, y2, curveGrad1, curveGrad2);
}

/// Reads a value out from the curve.
Sample operator ()(Sample x) {
if (x <= first.x) return first.y;
if (x >= last.x) return last.y;
size_t index = 1;
while (index < segments.size() && segments[index].start() <= x) {
++index;
}
return segments[index - 1](x);
}
};

Expand Down

0 comments on commit cb94481

Please sign in to comment.