diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..3f27406 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,15 @@ +cmake_minimum_required(VERSION 3.24) + +add_library(signalsmith-stretch INTERFACE) +set_target_properties(signalsmith-stretch PROPERTIES INTERFACE_INCLUDE_DIRECTORIES ${CMAKE_CURRENT_SOURCE_DIR}/include) + +include(FetchContent) +FetchContent_Declare( + signalsmith-linear + GIT_REPOSITORY https://github.com/Signalsmith-Audio/linear.git + GIT_TAG 355e014e330fde0c441004ba33e3964f04d4d28c + GIT_SHALLOW ON +) +FetchContent_MakeAvailable(signalsmith-linear) + +target_link_libraries(signalsmith-stretch INTERFACE signalsmith-linear) diff --git a/cmd/CMakeLists.txt b/cmd/CMakeLists.txt index d924577..21c6f05 100644 --- a/cmd/CMakeLists.txt +++ b/cmd/CMakeLists.txt @@ -1,8 +1,12 @@ cmake_minimum_required(VERSION 3.24) project(example-cmd VERSION 1.0.0) set(CMAKE_CXX_STANDARD 11) + add_executable(stretch ${CMAKE_CURRENT_SOURCE_DIR}/main.cpp) +add_subdirectory(".." "build-stretch") # Needs a (relative) build directory since it's not a subdir +target_link_libraries(stretch signalsmith-stretch) + include(FetchContent) set(FETCHCONTENT_QUIET False) FetchContent_Declare( @@ -11,3 +15,4 @@ FetchContent_Declare( SOURCE_DIR "${CMAKE_CURRENT_SOURCE_DIR}/inputs" ) FetchContent_MakeAvailable(inputs) + diff --git a/cmd/Makefile b/cmd/Makefile index 3321386..ad2649b 100644 --- a/cmd/Makefile +++ b/cmd/Makefile @@ -1,15 +1,11 @@ all: out/stretch -dev: clean out/stretch - mkdir -p out/examples - inputs/run-all.sh out/examples/u2- out/stretch --semitones=2 - cd ../linear/tests; make - -out/stretch: ../signalsmith-stretch.h main.cpp util/*.h util/*.hxx ../dsp/*.h +out/stretch: main.cpp ../signalsmith-stretch.h util/*.h util/*.hxx ../dsp/*.h mkdir -p out g++ -std=c++11 -O3 -g \ -Wall -Wextra -Wfatal-errors -Wpedantic -pedantic-errors \ -framework Accelerate -DSIGNALSMITH_USE_ACCELERATE \ + -I ../include \ main.cpp -o out/stretch # Uses input files from: https://signalsmith-audio.co.uk/code/stretch/inputs.zip @@ -37,8 +33,4 @@ clean: cmake: # CMAKE_BUILD_TYPE is needed for single-config generators (e.g. Makefiles) cmake -B build -DCMAKE_BUILD_TYPE=Release - cmake --build build --config Release - - # Copy to out/ so that `make clean cmake examples` works - mkdir -p out - cp build/Release/* out/ \ No newline at end of file + cmake --build build --config Release \ No newline at end of file diff --git a/cmd/main.cpp b/cmd/main.cpp index 6674209..596d37d 100644 --- a/cmd/main.cpp +++ b/cmd/main.cpp @@ -1,12 +1,13 @@ -#include "util/stopwatch.h" -#include "util/memory-tracker.hxx" - +// helper for debugging #include #define LOG_EXPR(expr) std::cout << #expr << " = " << (expr) << "\n"; -#include "../signalsmith-stretch.h" -#include "util/simple-args.h" -#include "util/wav.h" +#include "signalsmith-stretch/signalsmith-stretch.h" + +#include "./util/stopwatch.h" +#include "./util/memory-tracker.hxx" +#include "./util/simple-args.h" +#include "./util/wav.h" int main(int argc, char* argv[]) { signalsmith::stretch::SignalsmithStretch stretch; // optional cheaper RNG for performance comparison diff --git a/dsp/LICENSE.txt b/dsp/LICENSE.txt deleted file mode 100644 index 4e91ac9..0000000 --- a/dsp/LICENSE.txt +++ /dev/null @@ -1,21 +0,0 @@ -MIT License - -Copyright (c) 2021 Geraint Luff / Signalsmith Audio Ltd. - -Permission is hereby granted, free of charge, to any person obtaining a copy -of this software and associated documentation files (the "Software"), to deal -in the Software without restriction, including without limitation the rights -to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -copies of the Software, and to permit persons to whom the Software is -furnished to do so, subject to the following conditions: - -The above copyright notice and this permission notice shall be included in all -copies or substantial portions of the Software. - -THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE -SOFTWARE. \ No newline at end of file diff --git a/dsp/README.md b/dsp/README.md deleted file mode 100644 index 6363d74..0000000 --- a/dsp/README.md +++ /dev/null @@ -1,40 +0,0 @@ -# Signalsmith Audio's DSP Library - -A C++11 header-only library, providing classes/templates for (mostly audio) signal-processing tasks. - -More detail is in the [main project page](https://signalsmith-audio.co.uk/code/dsp/), and the [Doxygen docs](https://signalsmith-audio.co.uk/code/dsp/html/modules.html). - -## Basic use - -``` -git clone https://signalsmith-audio.co.uk/code/dsp.git -``` - -Just include the header file(s) you need, and start using classes: - -```cpp -#include "dsp/delay.h" - -using Delay = signalsmith::delay::Delay; -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, 6, 1) -``` - -### Development / contributing - -Tests (and source-scripts for the above docs) are available in a separate repo: - -``` -git clone https://signalsmith-audio.co.uk/code/dsp-doc.git -``` - -The goal (where possible) is to measure/test the actual audio characteristics of the tools (e.g. frequency responses and aliasing levels). - -### License - -This code is [MIT licensed](LICENSE.txt). If you'd prefer something else, get in touch. diff --git a/dsp/common.h b/dsp/common.h deleted file mode 100644 index 4526309..0000000 --- a/dsp/common.h +++ /dev/null @@ -1,47 +0,0 @@ -#ifndef SIGNALSMITH_DSP_COMMON_H -#define SIGNALSMITH_DSP_COMMON_H - -#if defined(__FAST_MATH__) && (__apple_build_version__ >= 16000000) && (__apple_build_version__ <= 16000099) -# error Apple Clang 16.0.0 generates incorrect SIMD for ARM. If you HAVE to use this version of Clang, turn off -ffast-math. -#endif - -#ifndef M_PI -#define M_PI 3.14159265358979323846264338327950288 -#endif - -namespace signalsmith { - /** @defgroup Common Common - @brief Definitions and helper classes used by the rest of the library - - @{ - @file - */ - -#define SIGNALSMITH_DSP_VERSION_MAJOR 1 -#define SIGNALSMITH_DSP_VERSION_MINOR 6 -#define SIGNALSMITH_DSP_VERSION_PATCH 1 -#define SIGNALSMITH_DSP_VERSION_STRING "1.6.1" - - /** Version compatability check. - \code{.cpp} - static_assert(signalsmith::version(1, 4, 1), "version check"); - \endcode - ... or use the equivalent `SIGNALSMITH_DSP_VERSION_CHECK`. - Major versions are not compatible with each other. Minor and patch versions are backwards-compatible. - */ - constexpr bool versionCheck(int major, int minor, int patch=0) { - return major == SIGNALSMITH_DSP_VERSION_MAJOR - && (SIGNALSMITH_DSP_VERSION_MINOR > minor - || (SIGNALSMITH_DSP_VERSION_MINOR == minor && SIGNALSMITH_DSP_VERSION_PATCH >= patch)); - } - -/// Check the library version is compatible (semver). -#define SIGNALSMITH_DSP_VERSION_CHECK(major, minor, patch) \ - static_assert(::signalsmith::versionCheck(major, minor, patch), "signalsmith library version is " SIGNALSMITH_DSP_VERSION_STRING); - -/** @} */ -} // 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 == 6 && SIGNALSMITH_DSP_VERSION_PATCH == 1, "multiple versions of the Signalsmith DSP library"); -#endif // include guard diff --git a/dsp/curves.h b/dsp/curves.h deleted file mode 100644 index 4a88e24..0000000 --- a/dsp/curves.h +++ /dev/null @@ -1,371 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_CURVES_H -#define SIGNALSMITH_DSP_CURVES_H - -#include -#include // std::stable_sort - -namespace signalsmith { -namespace curves { - /** @defgroup Curves Curves - @brief User-defined mapping functions - - @{ - @file - */ - - /// Linear map for real values. - template - class Linear { - Sample a1, a0; - public: - Linear() : Linear(0, 1) {} - Linear(Sample a0, Sample a1) : a1(a1), a0(a0) {} - /// Construct by from/to value pairs - Linear(Sample x0, Sample x1, Sample y0, Sample y1) : a1((x0 == x1) ? 0 : (y1 - y0)/(x1 - x0)), a0(y0 - x0*a1) {} - - Sample operator ()(Sample x) const { - return a0 + x*a1; - } - - Sample dx() const { - return a1; - } - - /// Returns the inverse map (with some numerical error) - Linear inverse() const { - Sample invA1 = 1/a1; - return Linear(-a0*invA1, invA1); - } - }; - - /// A real-valued cubic curve. It has a "start" point where accuracy is highest. - template - class Cubic { - Sample xStart, a0, a1, a2, a3; - - // Only use with y0 != y1 - static inline Sample gradient(Sample x0, Sample x1, Sample y0, Sample y1) { - return (y1 - y0)/(x1 - x0); - } - // 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; - } - } - } - // 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(0, curveGrad); - } else { - curveGrad = std::min(0, curveGrad); - } - } - } - public: - Cubic() : Cubic(0, 0, 0, 0, 0) {} - Cubic(Sample xStart, Sample a0, Sample a1, Sample a2, Sample a3) : xStart(xStart), a0(a0), a1(a1), a2(a2), a3(a3) {} - - Sample operator ()(Sample x) const { - x -= xStart; - return a0 + x*(a1 + x*(a2 + x*a3)); - } - /// The reference x-value, used as the centre of the cubic expansion - Sample start() const { - return xStart; - } - /// Differentiate - 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) { - 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 Cubic 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; - 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 && (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); - } 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); - } - }; - - /** Smooth interpolation (optionally monotonic) between points, using cubic segments. - \diagram{cubic-segments-example.svg,Example curve including a repeated point and an instantaneous jump. 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 - 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; - } - }; - std::vector points; - Point first{0, 0}, last{0, 0}; - - std::vector> _segments{1}; - // Not public because it's only valid inside the bounds - const Cubic & 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); - _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, 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(0, std::min(p1.curveGrad, p1.lineGrad*monotonicFactor)); - p2.curveGrad = std::max(0, std::min(p2.curveGrad, p1.lineGrad*monotonicFactor)); - } else { - p1.curveGrad = std::min(0, std::max(p1.curveGrad, p1.lineGrad*monotonicFactor)); - p2.curveGrad = std::min(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(); - 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 + (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 = 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; - std::vector & segments() { - return _segments; - } - const std::vector & segments() const { - return _segments; - } - }; - - /** A warped-range map, based on 1/x - \diagram{curves-reciprocal-example.svg}*/ - template - class Reciprocal { - 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 barkScale() { - return {1, 10, 24, 60, 1170, 13500}; - } - /// Returns a map from 0-1 to the given (non-negative) Hz range. - static Reciprocal 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) {} - Reciprocal(Sample x0, Sample x1, Sample x2, Sample y0, Sample y1, Sample y2) { - Sample kx = (x1 - x0)/(x2 - x1); - Sample ky = (y1 - y0)/(y2 - y1); - a = (kx*x2)*y0 - (ky*x0)*y2; - b = ky*y2 - kx*y0; - c = kx*x2 - ky*x0; - d = ky - kx; - } - - Sample operator ()(double x) const { - return (a + b*x)/(c + d*x); - } - Reciprocal inverse() const { - return Reciprocal(-a, c, b, -d); - } - 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 diff --git a/dsp/delay.h b/dsp/delay.h deleted file mode 100644 index 04c5631..0000000 --- a/dsp/delay.h +++ /dev/null @@ -1,717 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_DELAY_H -#define SIGNALSMITH_DSP_DELAY_H - -#include -#include -#include // for std::ceil() -#include - -#include -#include "./fft.h" -#include "./windows.h" - -namespace signalsmith { -namespace delay { - /** @defgroup Delay Delay utilities - @brief Standalone templated classes for delays - - You can set up a `Buffer` or `MultiBuffer`, and get interpolated samples using a `Reader` (separately on each channel in the multi-channel case) - or you can use `Delay`/`MultiDelay` which include their own buffers. - - Interpolation quality is chosen using a template class, from @ref Interpolators. - - @{ - @file - */ - - /** @brief Single-channel delay buffer - - Access is used with `buffer[]`, relative to the internal read/write position ("head"). This head is moved using `++buffer` (or `buffer += n`), such that `buffer[1] == (buffer + 1)[0]` in a similar way iterators/pointers. - - Operations like `buffer - 10` or `buffer++` return a View, which holds a fixed position in the buffer (based on the read/write position at the time). - - The capacity includes both positive and negative indices. For example, a capacity of 100 would support using any of the ranges: - - * `buffer[-99]` to buffer[0]` - * `buffer[-50]` to buffer[49]` - * `buffer[0]` to buffer[99]` - - Although buffers are usually used with historical samples accessed using negative indices e.g. `buffer[-10]`, you could equally use it flipped around (moving the head backwards through the buffer using `--buffer`). - */ - template - class Buffer { - unsigned bufferIndex; - unsigned bufferMask; - std::vector buffer; - public: - Buffer(int minCapacity=0) { - resize(minCapacity); - } - // We shouldn't accidentally copy a delay buffer - Buffer(const Buffer &other) = delete; - Buffer & operator =(const Buffer &other) = delete; - // But moving one is fine - Buffer(Buffer &&other) = default; - Buffer & operator =(Buffer &&other) = default; - - void resize(int minCapacity, Sample value=Sample()) { - int bufferLength = 1; - while (bufferLength < minCapacity) bufferLength *= 2; - buffer.assign(bufferLength, value); - bufferMask = unsigned(bufferLength - 1); - bufferIndex = 0; - } - void reset(Sample value=Sample()) { - buffer.assign(buffer.size(), value); - } - - /// Holds a view for a particular position in the buffer - template - class View { - using CBuffer = typename std::conditional::type; - using CSample = typename std::conditional::type; - CBuffer *buffer = nullptr; - unsigned bufferIndex = 0; - public: - View(CBuffer &buffer, int offset=0) : buffer(&buffer), bufferIndex(buffer.bufferIndex + (unsigned)offset) {} - View(const View &other, int offset=0) : buffer(other.buffer), bufferIndex(other.bufferIndex + (unsigned)offset) {} - View & operator =(const View &other) { - buffer = other.buffer; - bufferIndex = other.bufferIndex; - return *this; - } - - CSample & operator[](int offset) { - return buffer->buffer[(bufferIndex + (unsigned)offset)&buffer->bufferMask]; - } - const Sample & operator[](int offset) const { - return buffer->buffer[(bufferIndex + (unsigned)offset)&buffer->bufferMask]; - } - - /// Write data into the buffer - template - void write(Data &&data, int length) { - for (int i = 0; i < length; ++i) { - (*this)[i] = data[i]; - } - } - /// Read data out from the buffer - template - void read(int length, Data &&data) const { - for (int i = 0; i < length; ++i) { - data[i] = (*this)[i]; - } - } - - View operator +(int offset) const { - return View(*this, offset); - } - View operator -(int offset) const { - return View(*this, -offset); - } - }; - using MutableView = View; - using ConstView = View; - - MutableView view(int offset=0) { - return MutableView(*this, offset); - } - ConstView view(int offset=0) const { - return ConstView(*this, offset); - } - ConstView constView(int offset=0) const { - return ConstView(*this, offset); - } - - Sample & operator[](int offset) { - return buffer[(bufferIndex + (unsigned)offset)&bufferMask]; - } - const Sample & operator[](int offset) const { - return buffer[(bufferIndex + (unsigned)offset)&bufferMask]; - } - - /// Write data into the buffer - template - void write(Data &&data, int length) { - for (int i = 0; i < length; ++i) { - (*this)[i] = data[i]; - } - } - /// Read data out from the buffer - template - void read(int length, Data &&data) const { - for (int i = 0; i < length; ++i) { - data[i] = (*this)[i]; - } - } - - Buffer & operator ++() { - ++bufferIndex; - return *this; - } - Buffer & operator +=(int i) { - bufferIndex += (unsigned)i; - return *this; - } - Buffer & operator --() { - --bufferIndex; - return *this; - } - Buffer & operator -=(int i) { - bufferIndex -= (unsigned)i; - return *this; - } - - MutableView operator ++(int) { - MutableView view(*this); - ++bufferIndex; - return view; - } - MutableView operator +(int i) { - return MutableView(*this, i); - } - ConstView operator +(int i) const { - return ConstView(*this, i); - } - MutableView operator --(int) { - MutableView view(*this); - --bufferIndex; - return view; - } - MutableView operator -(int i) { - return MutableView(*this, -i); - } - ConstView operator -(int i) const { - return ConstView(*this, -i); - } - }; - - /** @brief Multi-channel delay buffer - - This behaves similarly to the single-channel `Buffer`, with the following differences: - - * `buffer[c]` returns a view for a single channel, which behaves like the single-channel `Buffer::View`. - * The constructor and `.resize()` take an additional first `channel` argument. - */ - template - class MultiBuffer { - int channels, stride; - Buffer buffer; - public: - using ConstChannel = typename Buffer::ConstView; - using MutableChannel = typename Buffer::MutableView; - - MultiBuffer(int channels=0, int capacity=0) : channels(channels), stride(capacity), buffer(channels*capacity) {} - - void resize(int nChannels, int capacity, Sample value=Sample()) { - channels = nChannels; - stride = capacity; - buffer.resize(channels*capacity, value); - } - void reset(Sample value=Sample()) { - buffer.reset(value); - } - - /// A reference-like multi-channel result for a particular sample index - template - class Stride { - using CChannel = typename std::conditional::type; - using CSample = typename std::conditional::type; - CChannel view; - int channels, stride; - public: - Stride(CChannel view, int channels, int stride) : view(view), channels(channels), stride(stride) {} - Stride(const Stride &other) : view(other.view), channels(other.channels), stride(other.stride) {} - - CSample & operator[](int channel) { - return view[channel*stride]; - } - const Sample & operator[](int channel) const { - return view[channel*stride]; - } - - /// Reads from the buffer into a multi-channel result - template - void get(Data &&result) const { - for (int c = 0; c < channels; ++c) { - result[c] = view[c*stride]; - } - } - /// Writes from multi-channel data into the buffer - template - void set(Data &&data) { - for (int c = 0; c < channels; ++c) { - view[c*stride] = data[c]; - } - } - template - Stride & operator =(const Data &data) { - set(data); - return *this; - } - Stride & operator =(const Stride &data) { - set(data); - return *this; - } - }; - - Stride at(int offset) { - return {buffer.view(offset), channels, stride}; - } - Stride at(int offset) const { - return {buffer.view(offset), channels, stride}; - } - - /// Holds a particular position in the buffer - template - class View { - using CChannel = typename std::conditional::type; - CChannel view; - int channels, stride; - public: - View(CChannel view, int channels, int stride) : view(view), channels(channels), stride(stride) {} - - CChannel operator[](int channel) { - return view + channel*stride; - } - ConstChannel operator[](int channel) const { - return view + channel*stride; - } - - Stride at(int offset) { - return {view + offset, channels, stride}; - } - Stride at(int offset) const { - return {view + offset, channels, stride}; - } - }; - using MutableView = View; - using ConstView = View; - - MutableView view(int offset=0) { - return MutableView(buffer.view(offset), channels, stride); - } - ConstView view(int offset=0) const { - return ConstView(buffer.view(offset), channels, stride); - } - ConstView constView(int offset=0) const { - return ConstView(buffer.view(offset), channels, stride); - } - - MutableChannel operator[](int channel) { - return buffer + channel*stride; - } - ConstChannel operator[](int channel) const { - return buffer + channel*stride; - } - - MultiBuffer & operator ++() { - ++buffer; - return *this; - } - MultiBuffer & operator +=(int i) { - buffer += i; - return *this; - } - MutableView operator ++(int) { - return MutableView(buffer++, channels, stride); - } - MutableView operator +(int i) { - return MutableView(buffer + i, channels, stride); - } - ConstView operator +(int i) const { - return ConstView(buffer + i, channels, stride); - } - MultiBuffer & operator --() { - --buffer; - return *this; - } - MultiBuffer & operator -=(int i) { - buffer -= i; - return *this; - } - MutableView operator --(int) { - return MutableView(buffer--, channels, stride); - } - MutableView operator -(int i) { - return MutableView(buffer - i, channels, stride); - } - ConstView operator -(int i) const { - return ConstView(buffer - i, channels, stride); - } - }; - - /** \defgroup Interpolators Interpolators - \ingroup Delay - @{ */ - /// Nearest-neighbour interpolator - /// \diagram{delay-random-access-nearest.svg,aliasing and maximum amplitude/delay errors for different input frequencies} - template - struct InterpolatorNearest { - static constexpr int inputLength = 1; - static constexpr Sample latency = -0.5; // Because we're truncating, which rounds down too often - - template - static Sample fractional(const Data &data, Sample) { - return data[0]; - } - }; - /// Linear interpolator - /// \diagram{delay-random-access-linear.svg,aliasing and maximum amplitude/delay errors for different input frequencies} - template - struct InterpolatorLinear { - static constexpr int inputLength = 2; - static constexpr int latency = 0; - - template - static Sample fractional(const Data &data, Sample fractional) { - Sample a = data[0], b = data[1]; - return a + fractional*(b - a); - } - }; - /// Spline cubic interpolator - /// \diagram{delay-random-access-cubic.svg,aliasing and maximum amplitude/delay errors for different input frequencies} - template - struct InterpolatorCubic { - static constexpr int inputLength = 4; - static constexpr int latency = 1; - - template - static Sample fractional(const Data &data, Sample fractional) { - // Cubic interpolation - Sample a = data[0], b = data[1], c = data[2], d = data[3]; - Sample cbDiff = c - b; - Sample k1 = (c - a)*0.5; - Sample k3 = k1 + (d - b)*0.5 - cbDiff*2; - Sample k2 = cbDiff - k3 - k1; - return b + fractional*(k1 + fractional*(k2 + fractional*k3)); // 16 ops total, not including the indexing - } - }; - - // Efficient Algorithms and Structures for Fractional Delay Filtering Based on Lagrange Interpolation - // Franck 2009 https://www.aes.org/e-lib/browse.cfm?elib=14647 - namespace _franck_impl { - template - struct ProductRange { - using Array = std::array; - static constexpr int mid = (low + high)/2; - using Left = ProductRange; - using Right = ProductRange; - - Left left; - Right right; - - const Sample total; - ProductRange(Sample x) : left(x), right(x), total(left.total*right.total) {} - - template - Sample calculateResult(Sample extraFactor, const Data &data, const Array &invFactors) { - return left.calculateResult(extraFactor*right.total, data, invFactors) - + right.calculateResult(extraFactor*left.total, data, invFactors); - } - }; - template - struct ProductRange { - using Array = std::array; - - const Sample total; - ProductRange(Sample x) : total(x - index) {} - - template - Sample calculateResult(Sample extraFactor, const Data &data, const Array &invFactors) { - return extraFactor*data[index]*invFactors[index]; - } - }; - } - /** Fixed-order Lagrange interpolation. - \diagram{interpolator-LagrangeN.svg,aliasing and amplitude/delay errors for different sizes} - */ - template - struct InterpolatorLagrangeN { - static constexpr int inputLength = n + 1; - static constexpr int latency = (n - 1)/2; - - using Array = std::array; - Array invDivisors; - - InterpolatorLagrangeN() { - for (int j = 0; j <= n; ++j) { - double divisor = 1; - for (int k = 0; k < j; ++k) divisor *= (j - k); - for (int k = j + 1; k <= n; ++k) divisor *= (j - k); - invDivisors[j] = 1/divisor; - } - } - - template - Sample fractional(const Data &data, Sample fractional) const { - constexpr int mid = n/2; - using Left = _franck_impl::ProductRange; - using Right = _franck_impl::ProductRange; - - Sample x = fractional + latency; - - Left left(x); - Right right(x); - - return left.calculateResult(right.total, data, invDivisors) + right.calculateResult(left.total, data, invDivisors); - } - }; - template - using InterpolatorLagrange3 = InterpolatorLagrangeN; - template - using InterpolatorLagrange7 = InterpolatorLagrangeN; - template - using InterpolatorLagrange19 = InterpolatorLagrangeN; - - /** Fixed-size Kaiser-windowed sinc interpolation. - \diagram{interpolator-KaiserSincN.svg,aliasing and amplitude/delay errors for different sizes} - If `minimumPhase` is enabled, a minimum-phase version of the kernel is used: - \diagram{interpolator-KaiserSincN-min.svg,aliasing and amplitude/delay errors for minimum-phase mode} - */ - template - struct InterpolatorKaiserSincN { - static constexpr int inputLength = n; - static constexpr Sample latency = minimumPhase ? 0 : (n*Sample(0.5) - 1); - - int subSampleSteps; - std::vector coefficients; - - InterpolatorKaiserSincN() : InterpolatorKaiserSincN(0.5 - 0.45/std::sqrt(n)) {} - InterpolatorKaiserSincN(double passFreq) : InterpolatorKaiserSincN(passFreq, 1 - passFreq) {} - InterpolatorKaiserSincN(double passFreq, double stopFreq) { - subSampleSteps = 2*n; // Heuristic again. Really it depends on the bandwidth as well. - double kaiserBandwidth = (stopFreq - passFreq)*(n + 1.0/subSampleSteps); - kaiserBandwidth += 1.25/kaiserBandwidth; // We want to place the first zero, but (because using this to window a sinc essentially integrates it in the freq-domain), our ripples (and therefore zeroes) are out of phase. This is a heuristic fix. - double sincScale = M_PI*(passFreq + stopFreq); - - double centreIndex = n*subSampleSteps*0.5, scaleFactor = 1.0/subSampleSteps; - std::vector windowedSinc(subSampleSteps*n + 1); - - ::signalsmith::windows::Kaiser::withBandwidth(kaiserBandwidth, false).fill(windowedSinc, windowedSinc.size()); - - for (size_t i = 0; i < windowedSinc.size(); ++i) { - double x = (i - centreIndex)*scaleFactor; - int intX = std::round(x); - if (intX != 0 && std::abs(x - intX) < 1e-6) { - // Exact 0s - windowedSinc[i] = 0; - } else if (std::abs(x) > 1e-6) { - double p = x*sincScale; - windowedSinc[i] *= std::sin(p)/p; - } - } - - if (minimumPhase) { - signalsmith::fft::FFT fft(windowedSinc.size()*2, 1); - windowedSinc.resize(fft.size(), 0); - std::vector> spectrum(fft.size()); - std::vector> cepstrum(fft.size()); - fft.fft(windowedSinc, spectrum); - for (size_t i = 0; i < fft.size(); ++i) { - spectrum[i] = std::log(std::abs(spectrum[i]) + 1e-30); - } - fft.fft(spectrum, cepstrum); - for (size_t i = 1; i < fft.size()/2; ++i) { - cepstrum[i] *= 0; - } - for (size_t i = fft.size()/2 + 1; i < fft.size(); ++i) { - cepstrum[i] *= 2; - } - Sample scaling = Sample(1)/fft.size(); - fft.ifft(cepstrum, spectrum); - - for (size_t i = 0; i < fft.size(); ++i) { - Sample phase = spectrum[i].imag()*scaling; - Sample mag = std::exp(spectrum[i].real()*scaling); - spectrum[i] = {mag*std::cos(phase), mag*std::sin(phase)}; - } - fft.ifft(spectrum, cepstrum); - windowedSinc.resize(subSampleSteps*n + 1); - windowedSinc.shrink_to_fit(); - for (size_t i = 0; i < windowedSinc.size(); ++i) { - windowedSinc[i] = cepstrum[i].real()*scaling; - } - } - - // Re-order into FIR fractional-delay blocks - coefficients.resize(n*(subSampleSteps + 1)); - for (int k = 0; k <= subSampleSteps; ++k) { - for (int i = 0; i < n; ++i) { - coefficients[k*n + i] = windowedSinc[(subSampleSteps - k) + i*subSampleSteps]; - } - } - } - - template - Sample fractional(const Data &data, Sample fractional) const { - Sample subSampleDelay = fractional*subSampleSteps; - int lowIndex = subSampleDelay; - if (lowIndex >= subSampleSteps) lowIndex = subSampleSteps - 1; - Sample subSampleFractional = subSampleDelay - lowIndex; - int highIndex = lowIndex + 1; - - Sample sumLow = 0, sumHigh = 0; - const Sample *coeffLow = coefficients.data() + lowIndex*n; - const Sample *coeffHigh = coefficients.data() + highIndex*n; - for (int i = 0; i < n; ++i) { - sumLow += data[i]*coeffLow[i]; - sumHigh += data[i]*coeffHigh[i]; - } - return sumLow + (sumHigh - sumLow)*subSampleFractional; - } - }; - - template - using InterpolatorKaiserSinc20 = InterpolatorKaiserSincN; - template - using InterpolatorKaiserSinc8 = InterpolatorKaiserSincN; - template - using InterpolatorKaiserSinc4 = InterpolatorKaiserSincN; - - template - using InterpolatorKaiserSinc20Min = InterpolatorKaiserSincN; - template - using InterpolatorKaiserSinc8Min = InterpolatorKaiserSincN; - template - using InterpolatorKaiserSinc4Min = InterpolatorKaiserSincN; - /// @} - - /** @brief A delay-line reader which uses an external buffer - - This is useful if you have multiple delay-lines reading from the same buffer. - */ - template class Interpolator=InterpolatorLinear> - class Reader : public Interpolator /* so we can get the empty-base-class optimisation */ { - using Super = Interpolator; - public: - Reader () {} - /// Pass in a configured interpolator - Reader (const Interpolator &interpolator) : Super(interpolator) {} - - template - Sample read(const Buffer &buffer, Sample delaySamples) const { - int startIndex = delaySamples; - Sample remainder = delaySamples - startIndex; - - // Delay buffers use negative indices, but interpolators use positive ones - using View = decltype(buffer - startIndex); - struct Flipped { - View view; - Sample operator [](int i) const { - return view[-i]; - } - }; - return Super::fractional(Flipped{buffer - startIndex}, remainder); - } - }; - - /** @brief A single-channel delay-line containing its own buffer.*/ - template class Interpolator=InterpolatorLinear> - class Delay : private Reader { - using Super = Reader; - Buffer buffer; - public: - static constexpr Sample latency = Super::latency; - - Delay(int capacity=0) : buffer(1 + capacity + Super::inputLength) {} - /// Pass in a configured interpolator - Delay(const Interpolator &interp, int capacity=0) : Super(interp), buffer(1 + capacity + Super::inputLength) {} - - void reset(Sample value=Sample()) { - buffer.reset(value); - } - void resize(int minCapacity, Sample value=Sample()) { - buffer.resize(minCapacity + Super::inputLength, value); - } - - /** Read a sample from `delaySamples` >= 0 in the past. - The interpolator may add its own latency on top of this (see `Delay::latency`). The default interpolation (linear) has 0 latency. - */ - Sample read(Sample delaySamples) const { - return Super::read(buffer, delaySamples); - } - /// Writes a sample. Returns the same object, so that you can say `delay.write(v).read(delay)`. - Delay & write(Sample value) { - ++buffer; - buffer[0] = value; - return *this; - } - }; - - /** @brief A multi-channel delay-line with its own buffer. */ - template class Interpolator=InterpolatorLinear> - class MultiDelay : private Reader { - using Super = Reader; - int channels; - MultiBuffer multiBuffer; - public: - static constexpr Sample latency = Super::latency; - - MultiDelay(int channels=0, int capacity=0) : channels(channels), multiBuffer(channels, 1 + capacity + Super::inputLength) {} - - void reset(Sample value=Sample()) { - multiBuffer.reset(value); - } - void resize(int nChannels, int capacity, Sample value=Sample()) { - channels = nChannels; - multiBuffer.resize(channels, capacity + Super::inputLength, value); - } - - /// A single-channel delay-line view, similar to a `const Delay` - struct ChannelView { - static constexpr Sample latency = Super::latency; - - const Super &reader; - typename MultiBuffer::ConstChannel channel; - - Sample read(Sample delaySamples) const { - return reader.read(channel, delaySamples); - } - }; - ChannelView operator [](int channel) const { - return ChannelView{*this, multiBuffer[channel]}; - } - - /// A multi-channel result, lazily calculating samples - struct DelayView { - Super &reader; - typename MultiBuffer::ConstView view; - Sample delaySamples; - - // Calculate samples on-the-fly - Sample operator [](int c) const { - return reader.read(view[c], delaySamples); - } - }; - DelayView read(Sample delaySamples) { - return DelayView{*this, multiBuffer.constView(), delaySamples}; - } - /// Reads into the provided output structure - template - void read(Sample delaySamples, Output &output) { - for (int c = 0; c < channels; ++c) { - output[c] = Super::read(multiBuffer[c], delaySamples); - } - } - /// Reads separate delays for each channel - template - void readMulti(const Delays &delays, Output &output) { - for (int c = 0; c < channels; ++c) { - output[c] = Super::read(multiBuffer[c], delays[c]); - } - } - template - MultiDelay & write(const Data &data) { - ++multiBuffer; - for (int c = 0; c < channels; ++c) { - multiBuffer[c][0] = data[c]; - } - return *this; - } - }; - -/** @} */ -}} // signalsmith::delay:: -#endif // include guard diff --git a/dsp/envelopes.h b/dsp/envelopes.h deleted file mode 100644 index 37c2b4d..0000000 --- a/dsp/envelopes.h +++ /dev/null @@ -1,523 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_ENVELOPES_H -#define SIGNALSMITH_DSP_ENVELOPES_H - -#include -#include -#include -#include - -namespace signalsmith { -namespace envelopes { - /** @defgroup Envelopes Envelopes and LFOs - @brief LFOs, envelopes and filters for manipulating them - - @{ - @file - */ - - /** An LFO based on cubic segments. - You can randomise the rate and/or the depth. Randomising the depth past `0.5` means it no longer neatly alternates sides: - \diagram{cubic-lfo-example.svg,Some example LFO curves.} - Without randomisation, it is approximately sine-like: - \diagram{cubic-lfo-spectrum-pure.svg} - */ - class CubicLfo { - float ratio = 0; - float ratioStep = 0; - - float valueFrom = 0, valueTo = 1, valueRange = 1; - float targetLow = 0, targetHigh = 1; - float targetRate = 0; - float rateRandom = 0.5, depthRandom = 0; - bool freshReset = true; - - std::default_random_engine randomEngine; - std::uniform_real_distribution randomUnit; - float random() { - return randomUnit(randomEngine); - } - float randomRate() { - return targetRate*exp(rateRandom*(random() - 0.5)); - } - float randomTarget(float previous) { - float randomOffset = depthRandom*random()*(targetLow - targetHigh); - if (previous < (targetLow + targetHigh)*0.5f) { - return targetHigh + randomOffset; - } else { - return targetLow - randomOffset; - } - } - public: - CubicLfo() : randomEngine(std::random_device()()), randomUnit(0, 1) { - reset(); - } - CubicLfo(long seed) : randomUnit(0, 1) { - randomEngine.seed(seed); - reset(); - } - - /// Resets the LFO state, starting with random phase. - void reset() { - ratio = random(); - ratioStep = randomRate(); - if (random() < 0.5) { - valueFrom = targetLow; - valueTo = targetHigh; - } else { - valueFrom = targetHigh; - valueTo = targetLow; - } - valueRange = valueTo - valueFrom; - freshReset = true; - } - /** Smoothly updates the LFO parameters. - - If called directly after `.reset()`, oscillation will immediately start within the specified range. Otherwise, it will remain smooth and fit within the new range after at most one cycle: - \diagram{cubic-lfo-changes.svg} - - The LFO will complete a full oscillation in (approximately) `1/rate` samples. `rateVariation` can be any number, but 0-1 is a good range. - - `depthVariation` must be in the range [0, 1], where ≤ 0.5 produces random amplitude but still alternates up/down. - \diagram{cubic-lfo-spectrum.svg,Spectra for the two types of randomisation - note the jump as depth variation goes past 50%} - */ - void set(float low, float high, float rate, float rateVariation=0, float depthVariation=0) { - rate *= 2; // We want to go up and down during this period - targetRate = rate; - targetLow = std::min(low, high); - targetHigh = std::max(low, high); - rateRandom = rateVariation; - depthRandom = std::min(1, std::max(0, depthVariation)); - - // If we haven't called .next() yet, don't bother being smooth. - if (freshReset) return reset(); - - // Only update the current rate if it's outside our new random-variation range - float maxRandomRatio = exp((float)0.5*rateRandom); - if (ratioStep > rate*maxRandomRatio || ratioStep < rate/maxRandomRatio) { - ratioStep = randomRate(); - } - } - - /// Returns the next output sample - float next() { - freshReset = false; - float result = ratio*ratio*(3 - 2*ratio)*valueRange + valueFrom; - - ratio += ratioStep; - while (ratio >= 1) { - ratio -= 1; - ratioStep = randomRate(); - valueFrom = valueTo; - valueTo = randomTarget(valueFrom); - valueRange = valueTo - valueFrom; - } - return result; - } - }; - - /** Variable-width rectangular sum */ - template - class BoxSum { - int bufferLength, index; - std::vector buffer; - Sample sum = 0, wrapJump = 0; - public: - BoxSum(int maxLength) { - resize(maxLength); - } - - /// Sets the maximum size (and reset contents) - void resize(int maxLength) { - bufferLength = maxLength + 1; - buffer.resize(bufferLength); - if (maxLength != 0) buffer.shrink_to_fit(); - reset(); - } - - /// Resets (with an optional "fill" value) - void reset(Sample value=Sample()) { - index = 0; - sum = 0; - for (size_t i = 0; i < buffer.size(); ++i) { - buffer[i] = sum; - sum += value; - } - wrapJump = sum; - sum = 0; - } - - Sample read(int width) { - int readIndex = index - width; - double result = sum; - if (readIndex < 0) { - result += wrapJump; - readIndex += bufferLength; - } - return result - buffer[readIndex]; - } - - void write(Sample value) { - ++index; - if (index == bufferLength) { - index = 0; - wrapJump = sum; - sum = 0; - } - sum += value; - buffer[index] = sum; - } - - Sample readWrite(Sample value, int width) { - write(value); - return read(width); - } - }; - - /** Rectangular moving average filter (FIR). - \diagram{box-filter-example.svg} - A filter of length 1 has order 0 (i.e. does nothing). */ - template - class BoxFilter { - BoxSum boxSum; - int _length, _maxLength; - Sample multiplier; - public: - BoxFilter(int maxLength) : boxSum(maxLength) { - resize(maxLength); - } - /// Sets the maximum size (and current size, and resets) - void resize(int maxLength) { - _maxLength = maxLength; - boxSum.resize(maxLength); - set(maxLength); - } - /// Sets the current size (expanding/allocating only if needed) - void set(int length) { - _length = length; - multiplier = Sample(1)/length; - if (length > _maxLength) resize(length); - } - - /// Resets (with an optional "fill" value) - void reset(Sample fill=Sample()) { - boxSum.reset(fill); - } - - Sample operator()(Sample v) { - return boxSum.readWrite(v, _length)*multiplier; - } - }; - - /** FIR filter made from a stack of `BoxFilter`s. - This filter has a non-negative impulse (monotonic step response), making it useful for smoothing positive-only values. It provides an optimal set of box-lengths, chosen to minimise peaks in the stop-band: - \diagram{box-stack-long.svg,Impulse responses for various stack sizes at length N=1000} - Since the underlying box-averages must have integer width, the peaks are slightly higher for shorter lengths with higher numbers of layers: - \diagram{box-stack-short-freq.svg,Frequency responses for various stack sizes at length N=30} - */ - template - class BoxStackFilter { - struct Layer { - double ratio = 0, lengthError = 0; - int length = 0; - BoxFilter filter{0}; - Layer() {} - }; - int _size; - std::vector layers; - - template - void setupLayers(const Iterable &ratios) { - layers.resize(0); - double sum = 0; - for (auto ratio : ratios) { - Layer layer; - layer.ratio = ratio; - layers.push_back(layer); - sum += ratio; - } - double factor = 1/sum; - for (auto &l : layers) { - l.ratio *= factor; - } - } - public: - BoxStackFilter(int maxSize, int layers=4) { - resize(maxSize, layers); - } - - /// Returns an optimal set of length ratios (heuristic for larger depths) - static std::vector optimalRatios(int layerCount) { - // Coefficients up to 6, found through numerical search - static double hardcoded[] = {1, 0.58224186169, 0.41775813831, 0.404078562416, 0.334851475794, 0.261069961789, 0.307944914938, 0.27369945234, 0.22913263601, 0.189222996712, 0.248329349789, 0.229253789144, 0.201191468123, 0.173033035122, 0.148192357821, 0.205275202874, 0.198413552119, 0.178256637764, 0.157821404506, 0.138663023387, 0.121570179349 /*, 0.178479592135, 0.171760666359, 0.158434068954, 0.143107825806, 0.125907148711, 0.11853946895, 0.103771229086, 0.155427880834, 0.153063152848, 0.142803459422, 0.131358358458, 0.104157805178, 0.119338029601, 0.0901675284678, 0.103683785192, 0.143949349747, 0.139813248378, 0.132051305252, 0.122216776152, 0.112888320989, 0.102534988632, 0.0928386714364, 0.0719750997699, 0.0817322396428, 0.130587011572, 0.127244563184, 0.121228748787, 0.113509941974, 0.105000272288, 0.0961938290157, 0.0880639725438, 0.0738389766046, 0.0746781936619, 0.0696544903682 */}; - if (layerCount <= 0) { - return {}; - } else if (layerCount <= 6) { - double *start = &hardcoded[layerCount*(layerCount - 1)/2]; - return std::vector(start, start + layerCount); - } - std::vector result(layerCount); - - double invN = 1.0/layerCount, sqrtN = std::sqrt(layerCount); - double p = 1 - invN; - double k = 1 + 4.5/sqrtN + 0.08*sqrtN; - - double sum = 0; - for (int i = 0; i < layerCount; ++i) { - double x = i*invN; - double power = -x*(1 - p*std::exp(-x*k)); - double length = std::pow(2, power); - result[i] = length; - sum += length; - } - double factor = 1/sum; - for (auto &r : result) r *= factor; - return result; - } - /** Approximate (optimal) bandwidth for a given number of layers - \diagram{box-stack-bandwidth.svg,Approximate main lobe width (bandwidth)} - */ - static constexpr double layersToBandwidth(int layers) { - return 1.58*(layers + 0.1); - } - /** Approximate (optimal) peak in the stop-band - \diagram{box-stack-peak.svg,Heuristic stop-band peak} - */ - static constexpr double layersToPeakDb(int layers) { - return 5 - layers*18; - } - - /// Sets size using an optimal (heuristic at larger sizes) set of length ratios - void resize(int maxSize, int layerCount) { - resize(maxSize, optimalRatios(layerCount)); - } - /// Sets the maximum (and current) impulse response length and explicit length ratios - template - auto resize(int maxSize, List ratios) -> decltype(void(std::begin(ratios)), void(std::end(ratios))) { - setupLayers(ratios); - for (auto &layer : layers) layer.filter.resize(0); // .set() will expand it later - _size = -1; - set(maxSize); - reset(); - } - void resize(int maxSize, std::initializer_list ratios) { - resize &>(maxSize, ratios); - } - - /// Sets the impulse response length (does not reset if `size` ≤ `maxSize`) - void set(int size) { - if (layers.size() == 0) return; // meaningless - - if (_size == size) return; - _size = size; - int order = size - 1; - int totalOrder = 0; - - for (auto &layer : layers) { - double layerOrderFractional = layer.ratio*order; - int layerOrder = int(layerOrderFractional); - layer.length = layerOrder + 1; - layer.lengthError = layerOrder - layerOrderFractional; - totalOrder += layerOrder; - } - // Round some of them up, so the total is correct - this is O(N²), but `layers.size()` is small - while (totalOrder < order) { - int minIndex = 0; - double minError = layers[0].lengthError; - for (size_t i = 1; i < layers.size(); ++i) { - if (layers[i].lengthError < minError) { - minError = layers[i].lengthError; - minIndex = i; - } - } - layers[minIndex].length++; - layers[minIndex].lengthError += 1; - totalOrder++; - } - for (auto &layer : layers) layer.filter.set(layer.length); - } - - /// Resets the filter - void reset(Sample fill=Sample()) { - for (auto &layer : layers) layer.filter.reset(fill); - } - - Sample operator()(Sample v) { - for (auto &layer : layers) { - v = layer.filter(v); - } - return v; - } - }; - - /** Peak-hold filter. - \diagram{peak-hold.svg} - - The size is variable, and can be changed instantly with `.set()`, or by using `.push()`/`.pop()` in an unbalanced way. - - This has complexity O(1) every sample when the length remains constant (balanced `.push()`/`.pop()`, or using `filter(v)`), and amortised O(1) complexity otherwise. To avoid allocations while running, it pre-allocates a vector (not a `std::deque`) which determines the maximum length. - */ - template - class PeakHold { - static constexpr Sample lowest = std::numeric_limits::lowest(); - int bufferMask; - std::vector buffer; - int backIndex = 0, middleStart = 0, workingIndex = 0, middleEnd = 0, frontIndex = 0; - Sample frontMax = lowest, workingMax = lowest, middleMax = lowest; - - public: - PeakHold(int maxLength) { - resize(maxLength); - } - int size() { - return frontIndex - backIndex; - } - void resize(int maxLength) { - int bufferLength = 1; - while (bufferLength < maxLength) bufferLength *= 2; - buffer.resize(bufferLength); - bufferMask = bufferLength - 1; - - frontIndex = backIndex + maxLength; - reset(); - } - void reset(Sample fill=lowest) { - int prevSize = size(); - buffer.assign(buffer.size(), fill); - frontMax = workingMax = middleMax = lowest; - middleEnd = workingIndex = frontIndex = 0; - middleStart = middleEnd - (prevSize/2); - backIndex = frontIndex - prevSize; - } - /** Sets the size immediately. - Must be `0 <= newSize <= maxLength` (see constructor and `.resize()`). - - Shrinking doesn't destroy information, and if you expand again (with `preserveCurrentPeak=false`), you will get the same output as before shrinking. Expanding when `preserveCurrentPeak` is enabled is destructive, re-writing its history such that the current output value is unchanged.*/ - void set(int newSize, bool preserveCurrentPeak=false) { - while (size() < newSize) { - Sample &backPrev = buffer[backIndex&bufferMask]; - --backIndex; - Sample &back = buffer[backIndex&bufferMask]; - back = preserveCurrentPeak ? backPrev : std::max(back, backPrev); - } - while (size() > newSize) { - pop(); - } - } - - void push(Sample v) { - buffer[frontIndex&bufferMask] = v; - ++frontIndex; - frontMax = std::max(frontMax, v); - } - void pop() { - if (backIndex == middleStart) { - // Move along the maximums - workingMax = lowest; - middleMax = frontMax; - frontMax = lowest; - - int prevFrontLength = frontIndex - middleEnd; - int prevMiddleLength = middleEnd - middleStart; - if (prevFrontLength <= prevMiddleLength + 1) { - // Swap over simply - middleStart = middleEnd; - middleEnd = frontIndex; - workingIndex = middleEnd; - } else { - // The front is longer than the middle - only happens if unbalanced - // We don't move *all* of the front over, keeping half the surplus in the front - int middleLength = (frontIndex - middleStart)/2; - middleStart = middleEnd; - middleEnd += middleLength; - - // Working index is close enough that it will be finished by the time the back is empty - int backLength = middleStart - backIndex; - int workingLength = std::min(backLength, middleEnd - middleStart); - workingIndex = middleStart + workingLength; - - // Since the front was not completely consumed, we re-calculate the front's maximum - for (int i = middleEnd; i != frontIndex; ++i) { - frontMax = std::max(frontMax, buffer[i&bufferMask]); - } - // The index might not start at the end of the working block - compute the last bit immediately - for (int i = middleEnd - 1; i != workingIndex - 1; --i) { - buffer[i&bufferMask] = workingMax = std::max(workingMax, buffer[i&bufferMask]); - } - } - - // Is the new back (previous middle) empty? Only happens if unbalanced - if (backIndex == middleStart) { - // swap over again (front's empty, no change) - workingMax = lowest; - middleMax = frontMax; - frontMax = lowest; - middleStart = workingIndex = middleEnd; - - if (backIndex == middleStart) { - --backIndex; // Only happens if you pop from an empty list - fail nicely - } - } - - buffer[frontIndex&bufferMask] = lowest; // In case of length 0, when everything points at this value - } - - ++backIndex; - if (workingIndex != middleStart) { - --workingIndex; - buffer[workingIndex&bufferMask] = workingMax = std::max(workingMax, buffer[workingIndex&bufferMask]); - } - } - Sample read() { - Sample backMax = buffer[backIndex&bufferMask]; - return std::max(backMax, std::max(middleMax, frontMax)); - } - - // For simple use as a constant-length filter - Sample operator ()(Sample v) { - push(v); - pop(); - return read(); - } - }; - - /** Peak-decay filter with a linear shape and fixed-time return to constant value. - \diagram{peak-decay-linear.svg} - This is equivalent to a `BoxFilter` which resets itself whenever the output would be less than the input. - */ - template - class PeakDecayLinear { - static constexpr Sample lowest = std::numeric_limits::lowest(); - PeakHold peakHold; - Sample value = lowest; - Sample stepMultiplier = 1; - public: - PeakDecayLinear(int maxLength) : peakHold(maxLength) { - set(maxLength); - } - void resize(int maxLength) { - peakHold.resize(maxLength); - reset(); - } - void set(double length) { - peakHold.set(std::ceil(length)); - // Overshoot slightly but don't exceed 1 - stepMultiplier = Sample(1.0001)/std::max(1.0001, length); - } - void reset(Sample start=lowest) { - peakHold.reset(start); - set(peakHold.size()); - value = start; - } - - Sample operator ()(Sample v) { - Sample peak = peakHold.read(); - peakHold(v); - return value = std::max(v, value + (v - peak)*stepMultiplier); - } - }; - -/** @} */ -}} // signalsmith::envelopes:: -#endif // include guard diff --git a/dsp/fft.h b/dsp/fft.h deleted file mode 100644 index 5897586..0000000 --- a/dsp/fft.h +++ /dev/null @@ -1,523 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_FFT_V5 -#define SIGNALSMITH_FFT_V5 - -#include "./perf.h" - -#include -#include -#include - -namespace signalsmith { namespace fft { - /** @defgroup FFT FFT (complex and real) - @brief Fourier transforms (complex and real) - - @{ - @file - */ - - namespace _fft_impl { - - template - SIGNALSMITH_INLINE V complexReal(const std::complex &c) { - return ((V*)(&c))[0]; - } - template - SIGNALSMITH_INLINE V complexImag(const std::complex &c) { - return ((V*)(&c))[1]; - } - - // Complex multiplication has edge-cases around Inf/NaN - handling those properly makes std::complex non-inlineable, so we use our own - template - SIGNALSMITH_INLINE std::complex complexMul(const std::complex &a, const std::complex &b) { - V aReal = complexReal(a), aImag = complexImag(a); - V bReal = complexReal(b), bImag = complexImag(b); - return conjugateSecond ? std::complex{ - bReal*aReal + bImag*aImag, - bReal*aImag - bImag*aReal - } : std::complex{ - aReal*bReal - aImag*bImag, - aReal*bImag + aImag*bReal - }; - } - - template - SIGNALSMITH_INLINE std::complex complexAddI(const std::complex &a, const std::complex &b) { - V aReal = complexReal(a), aImag = complexImag(a); - V bReal = complexReal(b), bImag = complexImag(b); - return flipped ? std::complex{ - aReal + bImag, - aImag - bReal - } : std::complex{ - aReal - bImag, - aImag + bReal - }; - } - - // Use SFINAE to get an iterator from std::begin(), if supported - otherwise assume the value itself is an iterator - template - struct GetIterator { - static T get(const T &t) { - return t; - } - }; - template - struct GetIterator()))> { - static auto get(const T &t) -> decltype(std::begin(t)) { - return std::begin(t); - } - }; - } - - /** Floating-point FFT implementation. - It is fast for 2^a * 3^b. - Here are the peak and RMS errors for `float`/`double` computation: - \diagram{fft-errors.svg Simulated errors for pure-tone harmonic inputs\, compared to a theoretical upper bound from "Roundoff error analysis of the fast Fourier transform" (G. Ramos, 1971)} - */ - template - class FFT { - using complex = std::complex; - size_t _size; - std::vector workingVector; - - enum class StepType { - generic, step2, step3, step4 - }; - struct Step { - StepType type; - size_t factor; - size_t startIndex; - size_t innerRepeats; - size_t outerRepeats; - size_t twiddleIndex; - }; - std::vector factors; - std::vector plan; - std::vector twiddleVector; - - struct PermutationPair {size_t from, to;}; - std::vector permutation; - - void addPlanSteps(size_t factorIndex, size_t start, size_t length, size_t repeats) { - if (factorIndex >= factors.size()) return; - - size_t factor = factors[factorIndex]; - if (factorIndex + 1 < factors.size()) { - if (factors[factorIndex] == 2 && factors[factorIndex + 1] == 2) { - ++factorIndex; - factor = 4; - } - } - - size_t subLength = length/factor; - Step mainStep{StepType::generic, factor, start, subLength, repeats, twiddleVector.size()}; - - if (factor == 2) mainStep.type = StepType::step2; - if (factor == 3) mainStep.type = StepType::step3; - if (factor == 4) mainStep.type = StepType::step4; - - // Twiddles - bool foundStep = false; - for (const Step &existingStep : plan) { - if (existingStep.factor == mainStep.factor && existingStep.innerRepeats == mainStep.innerRepeats) { - foundStep = true; - mainStep.twiddleIndex = existingStep.twiddleIndex; - break; - } - } - if (!foundStep) { - for (size_t i = 0; i < subLength; ++i) { - for (size_t f = 0; f < factor; ++f) { - double phase = 2*M_PI*i*f/length; - complex twiddle = {V(std::cos(phase)), V(-std::sin(phase))}; - twiddleVector.push_back(twiddle); - } - } - } - - if (repeats == 1 && sizeof(complex)*subLength > 65536) { - for (size_t i = 0; i < factor; ++i) { - addPlanSteps(factorIndex + 1, start + i*subLength, subLength, 1); - } - } else { - addPlanSteps(factorIndex + 1, start, subLength, repeats*factor); - } - plan.push_back(mainStep); - } - void setPlan() { - factors.resize(0); - size_t size = _size, factor = 2; - while (size > 1) { - if (size%factor == 0) { - factors.push_back(factor); - size /= factor; - } else if (factor > sqrt(size)) { - factor = size; - } else { - ++factor; - } - } - - plan.resize(0); - twiddleVector.resize(0); - addPlanSteps(0, 0, _size, 1); - twiddleVector.shrink_to_fit(); - - permutation.resize(0); - permutation.reserve(_size); - permutation.push_back(PermutationPair{0, 0}); - size_t indexLow = 0, indexHigh = factors.size(); - size_t inputStepLow = _size, outputStepLow = 1; - size_t inputStepHigh = 1, outputStepHigh = _size; - while (outputStepLow*inputStepHigh < _size) { - size_t f, inputStep, outputStep; - if (outputStepLow <= inputStepHigh) { - f = factors[indexLow++]; - inputStep = (inputStepLow /= f); - outputStep = outputStepLow; - outputStepLow *= f; - } else { - f = factors[--indexHigh]; - inputStep = inputStepHigh; - inputStepHigh *= f; - outputStep = (outputStepHigh /= f); - } - size_t oldSize = permutation.size(); - for (size_t i = 1; i < f; ++i) { - for (size_t j = 0; j < oldSize; ++j) { - PermutationPair pair = permutation[j]; - pair.from += i*inputStep; - pair.to += i*outputStep; - permutation.push_back(pair); - } - } - } - } - - template - void fftStepGeneric(RandomAccessIterator &&origData, const Step &step) { - complex *working = workingVector.data(); - const size_t stride = step.innerRepeats; - - for (size_t outerRepeat = 0; outerRepeat < step.outerRepeats; ++outerRepeat) { - RandomAccessIterator data = origData; - - const complex *twiddles = twiddleVector.data() + step.twiddleIndex; - const size_t factor = step.factor; - for (size_t repeat = 0; repeat < step.innerRepeats; ++repeat) { - for (size_t i = 0; i < step.factor; ++i) { - working[i] = _fft_impl::complexMul(data[i*stride], twiddles[i]); - } - for (size_t f = 0; f < factor; ++f) { - complex sum = working[0]; - for (size_t i = 1; i < factor; ++i) { - double phase = 2*M_PI*f*i/factor; - complex twiddle = {V(std::cos(phase)), V(-std::sin(phase))}; - sum += _fft_impl::complexMul(working[i], twiddle); - } - data[f*stride] = sum; - } - ++data; - twiddles += factor; - } - origData += step.factor*step.innerRepeats; - } - } - - template - SIGNALSMITH_INLINE void fftStep2(RandomAccessIterator &&origData, const Step &step) { - const size_t stride = step.innerRepeats; - const complex *origTwiddles = twiddleVector.data() + step.twiddleIndex; - for (size_t outerRepeat = 0; outerRepeat < step.outerRepeats; ++outerRepeat) { - const complex* twiddles = origTwiddles; - for (RandomAccessIterator data = origData; data < origData + stride; ++data) { - complex A = data[0]; - complex B = _fft_impl::complexMul(data[stride], twiddles[1]); - - data[0] = A + B; - data[stride] = A - B; - twiddles += 2; - } - origData += 2*stride; - } - } - - template - SIGNALSMITH_INLINE void fftStep3(RandomAccessIterator &&origData, const Step &step) { - constexpr complex factor3 = {-0.5, inverse ? 0.8660254037844386 : -0.8660254037844386}; - const size_t stride = step.innerRepeats; - const complex *origTwiddles = twiddleVector.data() + step.twiddleIndex; - - for (size_t outerRepeat = 0; outerRepeat < step.outerRepeats; ++outerRepeat) { - const complex* twiddles = origTwiddles; - for (RandomAccessIterator data = origData; data < origData + stride; ++data) { - complex A = data[0]; - complex B = _fft_impl::complexMul(data[stride], twiddles[1]); - complex C = _fft_impl::complexMul(data[stride*2], twiddles[2]); - - complex realSum = A + (B + C)*factor3.real(); - complex imagSum = (B - C)*factor3.imag(); - - data[0] = A + B + C; - data[stride] = _fft_impl::complexAddI(realSum, imagSum); - data[stride*2] = _fft_impl::complexAddI(realSum, imagSum); - - twiddles += 3; - } - origData += 3*stride; - } - } - - template - SIGNALSMITH_INLINE void fftStep4(RandomAccessIterator &&origData, const Step &step) { - const size_t stride = step.innerRepeats; - const complex *origTwiddles = twiddleVector.data() + step.twiddleIndex; - - for (size_t outerRepeat = 0; outerRepeat < step.outerRepeats; ++outerRepeat) { - const complex* twiddles = origTwiddles; - for (RandomAccessIterator data = origData; data < origData + stride; ++data) { - complex A = data[0]; - complex C = _fft_impl::complexMul(data[stride], twiddles[2]); - complex B = _fft_impl::complexMul(data[stride*2], twiddles[1]); - complex D = _fft_impl::complexMul(data[stride*3], twiddles[3]); - - complex sumAC = A + C, sumBD = B + D; - complex diffAC = A - C, diffBD = B - D; - - data[0] = sumAC + sumBD; - data[stride] = _fft_impl::complexAddI(diffAC, diffBD); - data[stride*2] = sumAC - sumBD; - data[stride*3] = _fft_impl::complexAddI(diffAC, diffBD); - - twiddles += 4; - } - origData += 4*stride; - } - } - - template - void permute(InputIterator input, OutputIterator data) { - for (auto pair : permutation) { - data[pair.from] = input[pair.to]; - } - } - - template - void run(InputIterator &&input, OutputIterator &&data) { - permute(input, data); - - for (const Step &step : plan) { - switch (step.type) { - case StepType::generic: - fftStepGeneric(data + step.startIndex, step); - break; - case StepType::step2: - fftStep2(data + step.startIndex, step); - break; - case StepType::step3: - fftStep3(data + step.startIndex, step); - break; - case StepType::step4: - fftStep4(data + step.startIndex, step); - break; - } - } - } - - static bool validSize(size_t size) { - constexpr static bool filter[32] = { - 1, 1, 1, 1, 1, 0, 1, 0, 1, 1, // 0-9 - 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, // 10-19 - 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, // 20-29 - 0, 0 - }; - return filter[size]; - } - public: - static size_t fastSizeAbove(size_t size) { - size_t power2 = 1; - while (size >= 32) { - size = (size - 1)/2 + 1; - power2 *= 2; - } - while (size < 32 && !validSize(size)) { - ++size; - } - return power2*size; - } - static size_t fastSizeBelow(size_t size) { - size_t power2 = 1; - while (size >= 32) { - size /= 2; - power2 *= 2; - } - while (size > 1 && !validSize(size)) { - --size; - } - return power2*size; - } - - FFT(size_t size, int fastDirection=0) : _size(0) { - if (fastDirection > 0) size = fastSizeAbove(size); - if (fastDirection < 0) size = fastSizeBelow(size); - this->setSize(size); - } - - size_t setSize(size_t size) { - if (size != _size) { - _size = size; - workingVector.resize(size); - setPlan(); - } - return _size; - } - size_t setFastSizeAbove(size_t size) { - return setSize(fastSizeAbove(size)); - } - size_t setFastSizeBelow(size_t size) { - return setSize(fastSizeBelow(size)); - } - const size_t & size() const { - return _size; - } - - template - void fft(InputIterator &&input, OutputIterator &&output) { - auto inputIter = _fft_impl::GetIterator::get(input); - auto outputIter = _fft_impl::GetIterator::get(output); - return run(inputIter, outputIter); - } - - template - void ifft(InputIterator &&input, OutputIterator &&output) { - auto inputIter = _fft_impl::GetIterator::get(input); - auto outputIter = _fft_impl::GetIterator::get(output); - return run(inputIter, outputIter); - } - }; - - struct FFTOptions { - static constexpr int halfFreqShift = 1; - }; - - template - class RealFFT { - static constexpr bool modified = (optionFlags&FFTOptions::halfFreqShift); - - using complex = std::complex; - std::vector complexBuffer1, complexBuffer2; - std::vector twiddlesMinusI; - std::vector modifiedRotations; - FFT complexFft; - public: - static size_t fastSizeAbove(size_t size) { - return FFT::fastSizeAbove((size + 1)/2)*2; - } - static size_t fastSizeBelow(size_t size) { - return FFT::fastSizeBelow(size/2)*2; - } - - RealFFT(size_t size=0, int fastDirection=0) : complexFft(0) { - if (fastDirection > 0) size = fastSizeAbove(size); - if (fastDirection < 0) size = fastSizeBelow(size); - this->setSize(std::max(size, 2)); - } - - size_t setSize(size_t size) { - complexBuffer1.resize(size/2); - complexBuffer2.resize(size/2); - - size_t hhSize = size/4 + 1; - twiddlesMinusI.resize(hhSize); - for (size_t i = 0; i < hhSize; ++i) { - V rotPhase = -2*M_PI*(modified ? i + 0.5 : i)/size; - twiddlesMinusI[i] = {std::sin(rotPhase), -std::cos(rotPhase)}; - } - if (modified) { - modifiedRotations.resize(size/2); - for (size_t i = 0; i < size/2; ++i) { - V rotPhase = -2*M_PI*i/size; - modifiedRotations[i] = {std::cos(rotPhase), std::sin(rotPhase)}; - } - } - - return complexFft.setSize(size/2); - } - size_t setFastSizeAbove(size_t size) { - return setSize(fastSizeAbove(size)); - } - size_t setFastSizeBelow(size_t size) { - return setSize(fastSizeBelow(size)); - } - size_t size() const { - return complexFft.size()*2; - } - - template - void fft(InputIterator &&input, OutputIterator &&output) { - size_t hSize = complexFft.size(); - for (size_t i = 0; i < hSize; ++i) { - if (modified) { - complexBuffer1[i] = _fft_impl::complexMul({input[2*i], input[2*i + 1]}, modifiedRotations[i]); - } else { - complexBuffer1[i] = {input[2*i], input[2*i + 1]}; - } - } - - complexFft.fft(complexBuffer1.data(), complexBuffer2.data()); - - if (!modified) output[0] = { - complexBuffer2[0].real() + complexBuffer2[0].imag(), - complexBuffer2[0].real() - complexBuffer2[0].imag() - }; - for (size_t i = modified ? 0 : 1; i <= hSize/2; ++i) { - size_t conjI = modified ? (hSize - 1 - i) : (hSize - i); - - complex odd = (complexBuffer2[i] + conj(complexBuffer2[conjI]))*(V)0.5; - complex evenI = (complexBuffer2[i] - conj(complexBuffer2[conjI]))*(V)0.5; - complex evenRotMinusI = _fft_impl::complexMul(evenI, twiddlesMinusI[i]); - - output[i] = odd + evenRotMinusI; - output[conjI] = conj(odd - evenRotMinusI); - } - } - - template - void ifft(InputIterator &&input, OutputIterator &&output) { - size_t hSize = complexFft.size(); - if (!modified) complexBuffer1[0] = { - input[0].real() + input[0].imag(), - input[0].real() - input[0].imag() - }; - for (size_t i = modified ? 0 : 1; i <= hSize/2; ++i) { - size_t conjI = modified ? (hSize - 1 - i) : (hSize - i); - complex v = input[i], v2 = input[conjI]; - - complex odd = v + conj(v2); - complex evenRotMinusI = v - conj(v2); - complex evenI = _fft_impl::complexMul(evenRotMinusI, twiddlesMinusI[i]); - - complexBuffer1[i] = odd + evenI; - complexBuffer1[conjI] = conj(odd - evenI); - } - - complexFft.ifft(complexBuffer1.data(), complexBuffer2.data()); - - for (size_t i = 0; i < hSize; ++i) { - complex v = complexBuffer2[i]; - if (modified) v = _fft_impl::complexMul(v, modifiedRotations[i]); - output[2*i] = v.real(); - output[2*i + 1] = v.imag(); - } - } - }; - - template - struct ModifiedRealFFT : public RealFFT { - using RealFFT::RealFFT; - }; - -/// @} -}} // namespace -#endif // include guard diff --git a/dsp/filters.h b/dsp/filters.h deleted file mode 100644 index 6e6a3c6..0000000 --- a/dsp/filters.h +++ /dev/null @@ -1,436 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_FILTERS_H -#define SIGNALSMITH_DSP_FILTERS_H - -#include "./perf.h" - -#include -#include - -namespace signalsmith { -namespace filters { - /** @defgroup Filters Basic filters - @brief Classes for some common filter types - - @{ - @file - */ - - /** Filter design methods. - These differ mostly in how they handle frequency-warping near Nyquist: - \diagram{filters-lowpass.svg} - \diagram{filters-highpass.svg} - \diagram{filters-peak.svg} - \diagram{filters-bandpass.svg} - \diagram{filters-notch.svg} - \diagram{filters-high-shelf.svg} - \diagram{filters-low-shelf.svg} - \diagram{filters-allpass.svg} - */ - enum class BiquadDesign { - bilinear, ///< Bilinear transform, adjusting for centre frequency but not bandwidth - cookbook, ///< RBJ's "Audio EQ Cookbook". Based on `bilinear`, adjusting bandwidth (for peak/notch/bandpass) to preserve the ratio between upper/lower boundaries. This performs oddly near Nyquist. - oneSided, ///< Based on `bilinear`, adjusting bandwidth to preserve the lower boundary (leaving the upper one loose). - vicanek ///< From Martin Vicanek's [Matched Second Order Digital Filters](https://vicanek.de/articles/BiquadFits.pdf). Falls back to `oneSided` for shelf and allpass filters. This takes the poles from the impulse-invariant approach, and then picks the zeros to create a better match. This means that Nyquist is not 0dB for peak/notch (or -Inf for lowpass), but it is a decent match to the analogue prototype. - }; - - /** A standard biquad. - - This is not guaranteed to be stable if modulated at audio rate. - - The default highpass/lowpass bandwidth (`defaultBandwidth`) produces a Butterworth filter when bandwidth-compensation is disabled. - - Bandwidth compensation defaults to `BiquadDesign::oneSided` (or `BiquadDesign::cookbook` if `cookbookBandwidth` is enabled) for all filter types aside from highpass/lowpass (which use `BiquadDesign::bilinear`).*/ - template - class BiquadStatic { - static constexpr BiquadDesign bwDesign = cookbookBandwidth ? BiquadDesign::cookbook : BiquadDesign::oneSided; - Sample a1 = 0, a2 = 0, b0 = 1, b1 = 0, b2 = 0; - Sample x1 = 0, x2 = 0, y1 = 0, y2 = 0; - - enum class Type {highpass, lowpass, highShelf, lowShelf, bandpass, notch, peak, allpass}; - - struct FreqSpec { - double scaledFreq; - double w0, sinW0, cosW0; - double inv2Q; - - FreqSpec(double freq, BiquadDesign design) { - scaledFreq = std::max(1e-6, std::min(0.4999, freq)); - if (design == BiquadDesign::cookbook) { - scaledFreq = std::min(0.45, scaledFreq); - } - w0 = 2*M_PI*scaledFreq; - cosW0 = std::cos(w0); - sinW0 = std::sin(w0); - } - - void oneSidedCompQ() { - // Ratio between our (digital) lower boundary f1 and centre f0 - double f1Factor = std::sqrt(inv2Q*inv2Q + 1) - inv2Q; - // Bilinear means discrete-time freq f = continuous-time freq tan(pi*xf/pi) - double ctF1 = std::tan(M_PI*scaledFreq*f1Factor), invCtF0 = (1 + cosW0)/sinW0; - double ctF1Factor = ctF1*invCtF0; - inv2Q = 0.5/ctF1Factor - 0.5*ctF1Factor; - } - }; - SIGNALSMITH_INLINE static FreqSpec octaveSpec(double scaledFreq, double octaves, BiquadDesign design) { - FreqSpec spec(scaledFreq, design); - - if (design == BiquadDesign::cookbook) { - // Approximately preserves bandwidth between halfway points - octaves *= spec.w0/spec.sinW0; - } - spec.inv2Q = std::sinh(std::log(2)*0.5*octaves); // 1/(2Q) - if (design == BiquadDesign::oneSided) spec.oneSidedCompQ(); - return spec; - } - SIGNALSMITH_INLINE static FreqSpec qSpec(double scaledFreq, double q, BiquadDesign design) { - FreqSpec spec(scaledFreq, design); - - spec.inv2Q = 0.5/q; - if (design == BiquadDesign::oneSided) spec.oneSidedCompQ(); - return spec; - } - - SIGNALSMITH_INLINE double dbToSqrtGain(double db) { - return std::pow(10, db*0.025); - } - - SIGNALSMITH_INLINE BiquadStatic & configure(Type type, FreqSpec calc, double sqrtGain, BiquadDesign design) { - double w0 = calc.w0; - - if (design == BiquadDesign::vicanek) { - if (type == Type::notch) { // Heuristic for notches near Nyquist - calc.inv2Q *= (1 - calc.scaledFreq*0.5); - } - double Q = (type == Type::peak ? 0.5*sqrtGain : 0.5)/calc.inv2Q; - double q = (type == Type::peak ? 1/sqrtGain : 1)*calc.inv2Q; - double expmqw = std::exp(-q*w0); - double da1, da2; - if (q <= 1) { - a1 = da1 = -2*expmqw*std::cos(std::sqrt(1 - q*q)*w0); - } else { - a1 = da1 = -2*expmqw*std::cosh(std::sqrt(q*q - 1)*w0); - } - a2 = da2 = expmqw*expmqw; - double sinpd2 = std::sin(w0/2); - double p0 = 1 - sinpd2*sinpd2, p1 = sinpd2*sinpd2, p2 = 4*p0*p1; - double A0 = 1 + da1 + da2, A1 = 1 - da1 + da2, A2 = -4*da2; - A0 *= A0; - A1 *= A1; - if (type == Type::lowpass) { - double R1 = (A0*p0 + A1*p1 + A2*p2)*Q*Q; - double B0 = A0, B1 = (R1 - B0*p0)/p1; - b0 = 0.5*(std::sqrt(B0) + std::sqrt(std::max(0.0, B1))); - b1 = std::sqrt(B0) - b0; - b2 = 0; - return *this; - } else if (type == Type::highpass) { - b2 = b0 = std::sqrt(A0*p0 + A1*p1 + A2*p2)*Q/(4*p1); - b1 = -2*b0; - return *this; - } else if (type == Type::bandpass) { - double R1 = A0*p0 + A1*p1 + A2*p2; - double R2 = -A0 + A1 + 4*(p0 - p1)*A2; - double B2 = (R1 - R2*p1)/(4*p1*p1); - double B1 = R2 + 4*(p1 - p0)*B2; - b1 = -0.5*std::sqrt(std::max(0.0, B1)); - b0 = 0.5*(std::sqrt(std::max(0.0, B2 + 0.25*B1)) - b1); - b2 = -b0 - b1; - return *this; - } else if (type == Type::notch) { - // The Vicanek paper doesn't cover notches (band-stop), but we know where the zeros should be: - b0 = 1; - double db1 = -2*std::cos(w0); // might be higher precision - b1 = db1; - b2 = 1; - // Scale so that B0 == A0 to get 0dB at f=0 - double scale = std::sqrt(A0)/(b0 + db1 + b2); - b0 *= scale; - b1 *= scale; - b2 *= scale; - return *this; - } else if (type == Type::peak) { - double G2 = (sqrtGain*sqrtGain)*(sqrtGain*sqrtGain); - double R1 = (A0*p0 + A1*p1 + A2*p2)*G2; - double R2 = (-A0 + A1 + 4*(p0 - p1)*A2)*G2; - double B0 = A0; - double B2 = (R1 - R2*p1 - B0)/(4*p1*p1); - double B1 = R2 + B0 + 4*(p1 - p0)*B2; - double W = 0.5*(std::sqrt(B0) + std::sqrt(std::max(0.0, B1))); - b0 = 0.5*(W + std::sqrt(std::max(0.0, W*W + B2))); - b1 = 0.5*(std::sqrt(B0) - std::sqrt(std::max(0.0, B1))); - b2 = -B2/(4*b0); - return *this; - } - // All others fall back to `oneSided` - design = BiquadDesign::oneSided; - calc.oneSidedCompQ(); - } - - double alpha = calc.sinW0*calc.inv2Q; - double A = sqrtGain, sqrtA2alpha = 2*std::sqrt(A)*alpha; - - double a0; - if (type == Type::highpass) { - b1 = -1 - calc.cosW0; - b0 = b2 = (1 + calc.cosW0)*0.5; - a0 = 1 + alpha; - a1 = -2*calc.cosW0; - a2 = 1 - alpha; - } else if (type == Type::lowpass) { - b1 = 1 - calc.cosW0; - b0 = b2 = b1*0.5; - a0 = 1 + alpha; - a1 = -2*calc.cosW0; - a2 = 1 - alpha; - } else if (type == Type::highShelf) { - b0 = A*((A+1)+(A-1)*calc.cosW0+sqrtA2alpha); - b2 = A*((A+1)+(A-1)*calc.cosW0-sqrtA2alpha); - b1 = -2*A*((A-1)+(A+1)*calc.cosW0); - a0 = (A+1)-(A-1)*calc.cosW0+sqrtA2alpha; - a2 = (A+1)-(A-1)*calc.cosW0-sqrtA2alpha; - a1 = 2*((A-1)-(A+1)*calc.cosW0); - } else if (type == Type::lowShelf) { - b0 = A*((A+1)-(A-1)*calc.cosW0+sqrtA2alpha); - b2 = A*((A+1)-(A-1)*calc.cosW0-sqrtA2alpha); - b1 = 2*A*((A-1)-(A+1)*calc.cosW0); - a0 = (A+1)+(A-1)*calc.cosW0+sqrtA2alpha; - a2 = (A+1)+(A-1)*calc.cosW0-sqrtA2alpha; - a1 = -2*((A-1)+(A+1)*calc.cosW0); - } else if (type == Type::bandpass) { - b0 = alpha; - b1 = 0; - b2 = -alpha; - a0 = 1 + alpha; - a1 = -2*calc.cosW0; - a2 = 1 - alpha; - } else if (type == Type::notch) { - b0 = 1; - b1 = -2*calc.cosW0; - b2 = 1; - a0 = 1 + alpha; - a1 = b1; - a2 = 1 - alpha; - } else if (type == Type::peak) { - b0 = 1 + alpha*A; - b1 = -2*calc.cosW0; - b2 = 1 - alpha*A; - a0 = 1 + alpha/A; - a1 = b1; - a2 = 1 - alpha/A; - } else if (type == Type::allpass) { - a0 = b2 = 1 + alpha; - a1 = b1 = -2*calc.cosW0; - a2 = b0 = 1 - alpha; - } else { - // reset to neutral - a1 = a2 = b1 = b2 = 0; - a0 = b0 = 1; - } - double invA0 = 1/a0; - b0 *= invA0; - b1 *= invA0; - b2 *= invA0; - a1 *= invA0; - a2 *= invA0; - return *this; - } - public: - static constexpr double defaultQ = 0.7071067811865476; // sqrt(0.5) - static constexpr double defaultBandwidth = 1.8999686269529916; // equivalent to above Q - - Sample operator ()(Sample x0) { - Sample y0 = x0*b0 + x1*b1 + x2*b2 - y1*a1 - y2*a2; - y2 = y1; - y1 = y0; - x2 = x1; - x1 = x0; - return y0; - } - - void reset() { - x1 = x2 = y1 = y2 = 0; - } - - std::complex response(Sample scaledFreq) const { - Sample w = scaledFreq*Sample(2*M_PI); - std::complex invZ = {std::cos(w), -std::sin(w)}, invZ2 = invZ*invZ; - return (b0 + invZ*b1 + invZ2*b2)/(Sample(1) + invZ*a1 + invZ2*a2); - } - Sample responseDb(Sample scaledFreq) const { - Sample w = scaledFreq*Sample(2*M_PI); - std::complex invZ = {std::cos(w), -std::sin(w)}, invZ2 = invZ*invZ; - Sample energy = std::norm(b0 + invZ*b1 + invZ2*b2)/std::norm(Sample(1) + invZ*a1 + invZ2*a2); - return 10*std::log10(energy); - } - - /// @name Lowpass - /// @{ - BiquadStatic & lowpass(double scaledFreq, double octaves=defaultBandwidth, BiquadDesign design=BiquadDesign::bilinear) { - return configure(Type::lowpass, octaveSpec(scaledFreq, octaves, design), 0, design); - } - BiquadStatic & lowpassQ(double scaledFreq, double q, BiquadDesign design=BiquadDesign::bilinear) { - return configure(Type::lowpass, qSpec(scaledFreq, q, design), 0, design); - } - /// @deprecated use `BiquadDesign` instead - void lowpass(double scaledFreq, double octaves, bool correctBandwidth) { - lowpass(scaledFreq, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @deprecated By the time you care about `design`, you should care about the bandwidth - BiquadStatic & lowpass(double scaledFreq, BiquadDesign design) { - return lowpass(scaledFreq, defaultBandwidth, design); - } - /// @} - - /// @name Highpass - /// @{ - BiquadStatic & highpass(double scaledFreq, double octaves=defaultBandwidth, BiquadDesign design=BiquadDesign::bilinear) { - return configure(Type::highpass, octaveSpec(scaledFreq, octaves, design), 0, design); - } - BiquadStatic & highpassQ(double scaledFreq, double q, BiquadDesign design=BiquadDesign::bilinear) { - return configure(Type::highpass, qSpec(scaledFreq, q, design), 0, design); - } - /// @deprecated use `BiquadDesign` instead - void highpass(double scaledFreq, double octaves, bool correctBandwidth) { - highpass(scaledFreq, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @deprecated By the time you care about `design`, you should care about the bandwidth - BiquadStatic & highpass(double scaledFreq, BiquadDesign design) { - return highpass(scaledFreq, defaultBandwidth, design); - } - /// @} - - /// @name Bandpass - /// @{ - BiquadStatic & bandpass(double scaledFreq, double octaves=defaultBandwidth, BiquadDesign design=bwDesign) { - return configure(Type::bandpass, octaveSpec(scaledFreq, octaves, design), 0, design); - } - BiquadStatic & bandpassQ(double scaledFreq, double q, BiquadDesign design=bwDesign) { - return configure(Type::bandpass, qSpec(scaledFreq, q, design), 0, design); - } - /// @deprecated use `BiquadDesign` instead - void bandpass(double scaledFreq, double octaves, bool correctBandwidth) { - bandpass(scaledFreq, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @deprecated By the time you care about `design`, you should care about the bandwidth - BiquadStatic & bandpass(double scaledFreq, BiquadDesign design) { - return bandpass(scaledFreq, defaultBandwidth, design); - } - /// @} - - /// @name Notch - /// @{ - BiquadStatic & notch(double scaledFreq, double octaves=defaultBandwidth, BiquadDesign design=bwDesign) { - return configure(Type::notch, octaveSpec(scaledFreq, octaves, design), 0, design); - } - BiquadStatic & notchQ(double scaledFreq, double q, BiquadDesign design=bwDesign) { - return configure(Type::notch, qSpec(scaledFreq, q, design), 0, design); - } - /// @deprecated use `BiquadDesign` instead - void notch(double scaledFreq, double octaves, bool correctBandwidth) { - notch(scaledFreq, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @deprecated By the time you care about `design`, you should care about the bandwidth - BiquadStatic & notch(double scaledFreq, BiquadDesign design) { - return notch(scaledFreq, defaultBandwidth, design); - } - /// @deprecated alias for `.notch()` - void bandStop(double scaledFreq, double octaves=1, bool correctBandwidth=true) { - notch(scaledFreq, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @} - - /// @name Peak - /// @{ - BiquadStatic & peak(double scaledFreq, double gain, double octaves=1, BiquadDesign design=bwDesign) { - return configure(Type::peak, octaveSpec(scaledFreq, octaves, design), std::sqrt(gain), design); - } - BiquadStatic & peakDb(double scaledFreq, double db, double octaves=1, BiquadDesign design=bwDesign) { - return configure(Type::peak, octaveSpec(scaledFreq, octaves, design), dbToSqrtGain(db), design); - } - BiquadStatic & peakQ(double scaledFreq, double gain, double q, BiquadDesign design=bwDesign) { - return configure(Type::peak, qSpec(scaledFreq, q, design), std::sqrt(gain), design); - } - BiquadStatic & peakDbQ(double scaledFreq, double db, double q, BiquadDesign design=bwDesign) { - return configure(Type::peak, qSpec(scaledFreq, q, design), dbToSqrtGain(db), design); - } - /// @deprecated By the time you care about `design`, you should care about the bandwidth - BiquadStatic & peak(double scaledFreq, double gain, BiquadDesign design) { - return peak(scaledFreq, gain, 1, design); - } - /// @} - - /// @name High shelf - /// @{ - BiquadStatic & highShelf(double scaledFreq, double gain, double octaves=defaultBandwidth, BiquadDesign design=bwDesign) { - return configure(Type::highShelf, octaveSpec(scaledFreq, octaves, design), std::sqrt(gain), design); - } - BiquadStatic & highShelfDb(double scaledFreq, double db, double octaves=defaultBandwidth, BiquadDesign design=bwDesign) { - return configure(Type::highShelf, octaveSpec(scaledFreq, octaves, design), dbToSqrtGain(db), design); - } - BiquadStatic & highShelfQ(double scaledFreq, double gain, double q, BiquadDesign design=bwDesign) { - return configure(Type::highShelf, qSpec(scaledFreq, q, design), std::sqrt(gain), design); - } - BiquadStatic & highShelfDbQ(double scaledFreq, double db, double q, BiquadDesign design=bwDesign) { - return configure(Type::highShelf, qSpec(scaledFreq, q, design), dbToSqrtGain(db), design); - } - /// @deprecated use `BiquadDesign` instead - BiquadStatic & highShelf(double scaledFreq, double gain, double octaves, bool correctBandwidth) { - return highShelf(scaledFreq, gain, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @deprecated use `BiquadDesign` instead - BiquadStatic & highShelfDb(double scaledFreq, double db, double octaves, bool correctBandwidth) { - return highShelfDb(scaledFreq, db, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @} - - /// @name Low shelf - /// @{ - BiquadStatic & lowShelf(double scaledFreq, double gain, double octaves=2, BiquadDesign design=bwDesign) { - return configure(Type::lowShelf, octaveSpec(scaledFreq, octaves, design), std::sqrt(gain), design); - } - BiquadStatic & lowShelfDb(double scaledFreq, double db, double octaves=2, BiquadDesign design=bwDesign) { - return configure(Type::lowShelf, octaveSpec(scaledFreq, octaves, design), dbToSqrtGain(db), design); - } - BiquadStatic & lowShelfQ(double scaledFreq, double gain, double q, BiquadDesign design=bwDesign) { - return configure(Type::lowShelf, qSpec(scaledFreq, q, design), std::sqrt(gain), design); - } - BiquadStatic & lowShelfDbQ(double scaledFreq, double db, double q, BiquadDesign design=bwDesign) { - return configure(Type::lowShelf, qSpec(scaledFreq, q, design), dbToSqrtGain(db), design); - } - /// @deprecated use `BiquadDesign` instead - BiquadStatic & lowShelf(double scaledFreq, double gain, double octaves, bool correctBandwidth) { - return lowShelf(scaledFreq, gain, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @deprecated use `BiquadDesign` instead - BiquadStatic & lowShelfDb(double scaledFreq, double db, double octaves, bool correctBandwidth) { - return lowShelfDb(scaledFreq, db, octaves, correctBandwidth ? bwDesign : BiquadDesign::bilinear); - } - /// @} - - /// @name Allpass - /// @{ - BiquadStatic & allpass(double scaledFreq, double octaves=1, BiquadDesign design=bwDesign) { - return configure(Type::allpass, octaveSpec(scaledFreq, octaves, design), 0, design); - } - BiquadStatic & allpassQ(double scaledFreq, double q, BiquadDesign design=bwDesign) { - return configure(Type::allpass, qSpec(scaledFreq, q, design), 0, design); - } - /// @} - - BiquadStatic & addGain(double factor) { - b0 *= factor; - b1 *= factor; - b2 *= factor; - return *this; - } - BiquadStatic & addGainDb(double db) { - return addGain(std::pow(10, db*0.05)); - } - }; - - /** @} */ -}} // signalsmith::filters:: -#endif // include guard diff --git a/dsp/mix.h b/dsp/mix.h deleted file mode 100644 index c4b25b1..0000000 --- a/dsp/mix.h +++ /dev/null @@ -1,218 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_MULTI_CHANNEL_H -#define SIGNALSMITH_DSP_MULTI_CHANNEL_H - -#include - -namespace signalsmith { -namespace mix { - /** @defgroup Mix Multichannel mixing - @brief Utilities for stereo/multichannel mixing operations - - @{ - @file - */ - - /** @defgroup Matrices Orthogonal matrices - @brief Some common matrices used for audio - @ingroup Mix - @{ */ - - /// @brief Hadamard: high mixing levels, N log(N) operations - template - class Hadamard { - public: - static_assert(size >= 0, "Size must be positive (or -1 for dynamic)"); - /// Applies the matrix, scaled so it's orthogonal - template - static void inPlace(Data &&data) { - unscaledInPlace(data); - - Sample factor = scalingFactor(); - for (int c = 0; c < size; ++c) { - data[c] *= factor; - } - } - - /// Scaling factor applied to make it orthogonal - static Sample scalingFactor() { - /// TODO: test for C++20, or whatever makes this constexpr. Maybe a `#define` in `common.h`? - return std::sqrt(Sample(1)/(size ? size : 1)); - } - - /// Skips the scaling, so it's a matrix full of `1`s - template - static void unscaledInPlace(Data &&data) { - if (size <= 1) return; - constexpr int hSize = size/2; - - Hadamard::template unscaledInPlace(data); - Hadamard::template unscaledInPlace(data); - - for (int i = 0; i < hSize; ++i) { - Sample a = data[i + startIndex], b = data[i + startIndex + hSize]; - data[i + startIndex] = (a + b); - data[i + startIndex + hSize] = (a - b); - } - } - }; - /// @brief Hadamard with dynamic size - template - class Hadamard { - int size; - public: - Hadamard(int size) : size(size) {} - - /// Applies the matrix, scaled so it's orthogonal - template - void inPlace(Data &&data) const { - unscaledInPlace(data); - - Sample factor = scalingFactor(); - for (int c = 0; c < size; ++c) { - data[c] *= factor; - } - } - - /// Scaling factor applied to make it orthogonal - Sample scalingFactor() const { - return std::sqrt(Sample(1)/(size ? size : 1)); - } - - /// Skips the scaling, so it's a matrix full of `1`s - template - void unscaledInPlace(Data &&data) const { - int hSize = size/2; - while (hSize > 0) { - for (int startIndex = 0; startIndex < size; startIndex += hSize*2) { - for (int i = startIndex; i < startIndex + hSize; ++i) { - Sample a = data[i], b = data[i + hSize]; - data[i] = (a + b); - data[i + hSize] = (a - b); - } - } - hSize /= 2; - } - } - }; - /// @brief Householder: moderate mixing, 2N operations - template - class Householder { - public: - static_assert(size >= 0, "Size must be positive (or -1 for dynamic)"); - template - static void inPlace(Data &&data) { - if (size < 1) return; - /// TODO: test for C++20, which makes `std::complex::operator/` constexpr - const Sample factor = Sample(-2)/Sample(size ? size : 1); - - Sample sum = data[0]; - for (int i = 1; i < size; ++i) { - sum += data[i]; - } - sum *= factor; - for (int i = 0; i < size; ++i) { - data[i] += sum; - } - } - /// @deprecated The matrix is already orthogonal, but this is here for compatibility with Hadamard - constexpr static Sample scalingFactor() { - return 1; - } - }; - /// @brief Householder with dynamic size - template - class Householder { - int size; - public: - Householder(int size) : size(size) {} - - template - void inPlace(Data &&data) const { - if (size < 1) return; - const Sample factor = Sample(-2)/Sample(size ? size : 1); - - Sample sum = data[0]; - for (int i = 1; i < size; ++i) { - sum += data[i]; - } - sum *= factor; - for (int i = 0; i < size; ++i) { - data[i] += sum; - } - } - /// @deprecated The matrix is already orthogonal, but this is here for compatibility with Hadamard - constexpr static Sample scalingFactor() { - return 1; - } - }; - /// @} - - /** @brief Upmix/downmix a stereo signal to an (even) multi-channel signal - - When spreading out, it rotates the input by various amounts (e.g. a four-channel signal would produce `(left, right, mid side)`), such that energy is preserved for each pair. - - When mixing together, it uses the opposite rotations, such that upmix → downmix produces the same stereo signal (when scaled by `.scalingFactor1()`. - */ - template - class StereoMultiMixer { - static_assert((channels/2)*2 == channels, "StereoMultiMixer must have an even number of channels"); - static_assert(channels > 0, "StereoMultiMixer must have a positive number of channels"); - static constexpr int hChannels = channels/2; - std::array coeffs; - public: - StereoMultiMixer() { - coeffs[0] = 1; - coeffs[1] = 0; - for (int i = 1; i < hChannels; ++i) { - double phase = M_PI*i/channels; - coeffs[2*i] = std::cos(phase); - coeffs[2*i + 1] = std::sin(phase); - } - } - - template - void stereoToMulti(In &input, Out &output) const { - output[0] = input[0]; - output[1] = input[1]; - for (int i = 2; i < channels; i += 2) { - output[i] = input[0]*coeffs[i] + input[1]*coeffs[i + 1]; - output[i + 1] = input[1]*coeffs[i] - input[0]*coeffs[i + 1]; - } - } - template - void multiToStereo(In &input, Out &output) const { - output[0] = input[0]; - output[1] = input[1]; - for (int i = 2; i < channels; i += 2) { - output[0] += input[i]*coeffs[i] - input[i + 1]*coeffs[i + 1]; - output[1] += input[i + 1]*coeffs[i] + input[i]*coeffs[i + 1]; - } - } - /// Scaling factor for the downmix, if channels are phase-aligned - static constexpr Sample scalingFactor1() { - return 2/Sample(channels); - } - /// Scaling factor for the downmix, if channels are independent - static Sample scalingFactor2() { - return std::sqrt(scalingFactor1()); - } - }; - - /// A cheap (polynomial) almost-energy-preserving crossfade - /// Maximum energy error: 1.06%, average 0.64%, curves overshoot by 0.3% - /// See: http://signalsmith-audio.co.uk/writing/2021/cheap-energy-crossfade/ - template - void cheapEnergyCrossfade(Sample x, Result &toCoeff, Result &fromCoeff) { - Sample x2 = 1 - x; - // Other powers p can be approximated by: k = -6.0026608 + p*(6.8773512 - 1.5838104*p) - Sample A = x*x2, B = A*(1 + (Sample)1.4186*A); - Sample C = (B + x), D = (B + x2); - toCoeff = C*C; - fromCoeff = D*D; - } - -/** @} */ -}} // signalsmith::delay:: -#endif // include guard diff --git a/dsp/perf.h b/dsp/perf.h deleted file mode 100644 index 779c7b7..0000000 --- a/dsp/perf.h +++ /dev/null @@ -1,84 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_PERF_H -#define SIGNALSMITH_DSP_PERF_H - -#include - -#if defined(__SSE__) || defined(_M_X64) -# include -#else -# include // for uintptr_t -#endif - -namespace signalsmith { -namespace perf { - /** @defgroup Performance Performance helpers - @brief Nothing serious, just some `#defines` and helpers - - @{ - @file - */ - - /// *Really* insist that a function/method is inlined (mostly for performance in DEBUG builds) - #ifndef SIGNALSMITH_INLINE - #ifdef __GNUC__ - #define SIGNALSMITH_INLINE __attribute__((always_inline)) inline - #elif defined(__MSVC__) - #define SIGNALSMITH_INLINE __forceinline inline - #else - #define SIGNALSMITH_INLINE inline - #endif - #endif - - /** @brief Complex-multiplication (with optional conjugate second-arg), without handling NaN/Infinity - The `std::complex` multiplication has edge-cases around NaNs which slow things down and prevent auto-vectorisation. Flags like `-ffast-math` sort this out anyway, but this helps with Debug builds. - */ - template - SIGNALSMITH_INLINE static std::complex mul(const std::complex &a, const std::complex &b) { - return conjugateSecond ? std::complex{ - b.real()*a.real() + b.imag()*a.imag(), - b.real()*a.imag() - b.imag()*a.real() - } : std::complex{ - a.real()*b.real() - a.imag()*b.imag(), - a.real()*b.imag() + a.imag()*b.real() - }; - } - -#if defined(__SSE__) || defined(_M_X64) - class StopDenormals { - unsigned int controlStatusRegister; - public: - StopDenormals() : controlStatusRegister(_mm_getcsr()) { - _mm_setcsr(controlStatusRegister|0x8040); // Flush-to-Zero and Denormals-Are-Zero - } - ~StopDenormals() { - _mm_setcsr(controlStatusRegister); - } - }; -#elif (defined (__ARM_NEON) || defined (__ARM_NEON__)) - class StopDenormals { - uintptr_t status; - public: - StopDenormals() { - uintptr_t asmStatus; - asm volatile("mrs %0, fpcr" : "=r"(asmStatus)); - status = asmStatus = asmStatus|0x01000000U; // Flush to Zero - asm volatile("msr fpcr, %0" : : "ri"(asmStatus)); - } - ~StopDenormals() { - uintptr_t asmStatus = status; - asm volatile("msr fpcr, %0" : : "ri"(asmStatus)); - } - }; -#else -# if __cplusplus >= 202302L -# warning "The `StopDenormals` class doesn't do anything for this architecture" -# endif - class StopDenormals {}; // FIXME: add for other architectures -#endif - -/** @} */ -}} // signalsmith::perf:: - -#endif // include guard diff --git a/dsp/rates.h b/dsp/rates.h deleted file mode 100644 index b5da8c7..0000000 --- a/dsp/rates.h +++ /dev/null @@ -1,184 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_RATES_H -#define SIGNALSMITH_DSP_RATES_H - -#include "./windows.h" -#include "./delay.h" - -namespace signalsmith { -namespace rates { - /** @defgroup Rates Multi-rate processing - @brief Classes for oversampling/upsampling/downsampling etc. - - @{ - @file - */ - - /// @brief Fills a container with a Kaiser-windowed sinc for an FIR lowpass. - /// \diagram{rates-kaiser-sinc.svg,33-point results for various pass/stop frequencies} - template - void fillKaiserSinc(Data &&data, int length, double passFreq, double stopFreq) { - if (length <= 0) return; - double kaiserBandwidth = (stopFreq - passFreq)*length; - kaiserBandwidth += 1.25/kaiserBandwidth; // heuristic for transition band, see `InterpolatorKaiserSincN` - auto kaiser = signalsmith::windows::Kaiser::withBandwidth(kaiserBandwidth); - kaiser.fill(data, length); - - double centreIndex = (length - 1)*0.5; - double sincScale = M_PI*(passFreq + stopFreq); - double ampScale = (passFreq + stopFreq); - for (int i = 0; i < length; ++i) { - double x = (i - centreIndex), px = x*sincScale; - double sinc = (std::abs(px) > 1e-6) ? std::sin(px)*ampScale/px : ampScale; - data[i] *= sinc; - } - }; - /// @brief If only the centre frequency is specified, a heuristic is used to balance ripples and transition width - /// \diagram{rates-kaiser-sinc-heuristic.svg,The transition width is set to: 0.9/sqrt(length)} - template - void fillKaiserSinc(Data &&data, int length, double centreFreq) { - double halfWidth = 0.45/std::sqrt(length); - if (halfWidth > centreFreq) halfWidth = (halfWidth + centreFreq)*0.5; - fillKaiserSinc(data, length, centreFreq - halfWidth, centreFreq + halfWidth); - } - - /** 2x FIR oversampling for block-based processing. - - \diagram{rates-oversampler2xfir-responses-45.svg,Upsample response for various lengths} - - The oversampled signal is stored inside this object, with channels accessed via `oversampler[c]`. For example, you might do: - \code{.cpp} - // Upsample from multi-channel input (inputBuffers[c][i] is a sample) - oversampler.up(inputBuffers, bufferLength) - - // Modify the contents at the higher rate - for (int c = 0; c < 2; ++c) { - float *channel = oversampler[c]; - for (int i = 0; i < bufferLength*2; ++i) { - channel[i] = std::abs(channel[i]); - } - } - - // Downsample into the multi-channel output - oversampler.down(outputBuffers, bufferLength); - \endcode - - The performance depends not just on the length, but also on where you end the passband, allowing a wider/narrower transition band. Frequencies above this limit (relative to the lower sample-rate) may alias when upsampling and downsampling. - - \diagram{rates-oversampler2xfir-lengths.svg,Resample error rates for different passband thresholds} - - Since both upsample and downsample are stateful, channels are meaningful. If your input channel-count doesn't match your output, you can size it to the larger of the two, and use `.upChannel()` and `.downChannel()` to only process the channels which exist.*/ - template - struct Oversampler2xFIR { - Oversampler2xFIR() : Oversampler2xFIR(0, 0) {} - Oversampler2xFIR(int channels, int maxBlock, int halfLatency=16, double passFreq=0.43) { - resize(channels, maxBlock, halfLatency, passFreq); - } - - void resize(int nChannels, int maxBlockLength) { - resize(nChannels, maxBlockLength, oneWayLatency); - } - void resize(int nChannels, int maxBlockLength, int halfLatency, double passFreq=0.43) { - oneWayLatency = halfLatency; - kernelLength = oneWayLatency*2; - channels = nChannels; - halfSampleKernel.resize(kernelLength); - fillKaiserSinc(halfSampleKernel, kernelLength, passFreq, 1 - passFreq); - inputStride = kernelLength + maxBlockLength; - inputBuffer.resize(channels*inputStride); - stride = (maxBlockLength + kernelLength)*2; - buffer.resize(stride*channels); - } - - void reset() { - inputBuffer.assign(inputBuffer.size(), 0); - buffer.assign(buffer.size(), 0); - } - - /// @brief Round-trip latency (or equivalently: upsample latency at the higher rate). - /// This will be twice the value passed into the constructor or `.resize()`. - int latency() const { - return kernelLength; - } - - /// Upsamples from a multi-channel input into the internal buffer - template - void up(Data &&data, int lowSamples) { - for (int c = 0; c < channels; ++c) { - upChannel(c, data[c], lowSamples); - } - } - - /// Upsamples a single-channel input into the internal buffer - template - void upChannel(int c, Data &&data, int lowSamples) { - Sample *inputChannel = inputBuffer.data() + c*inputStride; - for (int i = 0; i < lowSamples; ++i) { - inputChannel[kernelLength + i] = data[i]; - } - Sample *output = (*this)[c]; - for (int i = 0; i < lowSamples; ++i) { - output[2*i] = inputChannel[i + oneWayLatency]; - Sample *offsetInput = inputChannel + (i + 1); - Sample sum = 0; - for (int o = 0; o < kernelLength; ++o) { - sum += offsetInput[o]*halfSampleKernel[o]; - } - output[2*i + 1] = sum; - } - // Copy the end of the buffer back to the beginning - for (int i = 0; i < kernelLength; ++i) { - inputChannel[i] = inputChannel[lowSamples + i]; - } - } - - /// Downsamples from the internal buffer to a multi-channel output - template - void down(Data &&data, int lowSamples) { - for (int c = 0; c < channels; ++c) { - downChannel(c, data[c], lowSamples); - } - } - - /// Downsamples a single channel from the internal buffer to a single-channel output - template - void downChannel(int c, Data &&data, int lowSamples) { - Sample *input = buffer.data() + c*stride; // no offset for latency - for (int i = 0; i < lowSamples; ++i) { - Sample v1 = input[2*i + kernelLength]; - Sample sum = 0; - for (int o = 0; o < kernelLength; ++o) { - Sample v2 = input[2*(i + o) + 1]; - sum += v2*halfSampleKernel[o]; - } - Sample v2 = sum; - Sample v = (v1 + v2)*Sample(0.5); - data[i] = v; - } - // Copy the end of the buffer back to the beginning - for (int i = 0; i < kernelLength*2; ++i) { - input[i] = input[lowSamples*2 + i]; - } - } - - /// Gets the samples for a single (higher-rate) channel. The valid length depends how many input samples were passed into `.up()`/`.upChannel()`. - Sample * operator[](int c) { - return buffer.data() + kernelLength*2 + stride*c; - } - const Sample * operator[](int c) const { - return buffer.data() + kernelLength*2 + stride*c; - } - - private: - int oneWayLatency, kernelLength; - int channels; - int stride, inputStride; - std::vector inputBuffer; - std::vector halfSampleKernel; - std::vector buffer; - }; - -/** @} */ -}} // namespace -#endif // include guard diff --git a/dsp/spectral.h b/dsp/spectral.h deleted file mode 100644 index effc75e..0000000 --- a/dsp/spectral.h +++ /dev/null @@ -1,496 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_SPECTRAL_H -#define SIGNALSMITH_DSP_SPECTRAL_H - -#include "./perf.h" -#include "./fft.h" -#include "./windows.h" -#include "./delay.h" - -#include - -namespace signalsmith { -namespace spectral { - /** @defgroup Spectral Spectral Processing - @brief Tools for frequency-domain manipulation of audio signals - - @{ - @file - */ - - /** @brief An FFT with built-in windowing and round-trip scaling - - This uses a Modified Real FFT, which applies half-bin shift before the transform. The result therefore has `N/2` bins, centred at the frequencies: `(i + 0.5)/N`. - - This avoids the awkward (real-valued) bands for DC-offset and Nyquist. - */ - template - class WindowedFFT { - using MRFFT = signalsmith::fft::ModifiedRealFFT; - using Complex = std::complex; - MRFFT mrfft{2}; - - std::vector fftWindow; - std::vector timeBuffer; - int offsetSamples = 0; - public: - /// Returns a fast FFT size <= `size` - static int fastSizeAbove(int size, int divisor=1) { - return MRFFT::fastSizeAbove(size/divisor)*divisor; - } - /// Returns a fast FFT size >= `size` - static int fastSizeBelow(int size, int divisor=1) { - return MRFFT::fastSizeBelow(1 + (size - 1)/divisor)*divisor; - } - - WindowedFFT() {} - WindowedFFT(int size, int rotateSamples=0) { - setSize(size, rotateSamples); - } - template - 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 & setSizeWindow(int size, int rotateSamples=0) { - mrfft.setSize(size); - fftWindow.assign(size, 1); - timeBuffer.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 - 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) { - Sample r = (i + windowOffset)*invSize; - fftWindow[i] = fn(r); - } - } - /// Sets the size (using the default Blackman-Harris window) - 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.01168*std::cos(phase*3); - }, Sample(0.5), rotateSamples); - } - - const std::vector & window() const { - return this->fftWindow; - } - int size() const { - return mrfft.size(); - } - - /// Performs an FFT, with windowing and rotation (if enabled) - template - void fft(Input &&input, Output &&output) { - int fftSize = size(); - const Sample norm = (withScaling ? 1/(Sample)fftSize : 1); - for (int i = 0; i < offsetSamples; ++i) { - // Inverted polarity since we're using the MRFFT - timeBuffer[i + fftSize - offsetSamples] = -input[i]*norm*(withWindow ? fftWindow[i] : Sample(1)); - } - for (int i = offsetSamples; i < fftSize; ++i) { - timeBuffer[i - offsetSamples] = input[i]*norm*(withWindow ? fftWindow[i] : Sample(1)); - } - mrfft.fft(timeBuffer, output); - } - /// Performs an FFT (no windowing or rotation) - template - void fftRaw(Input &&input, Output &&output) { - mrfft.fft(input, output); - } - - /// Inverse FFT, with windowing, 1/N scaling and rotation (if enabled) - template - void ifft(Input &&input, Output &&output) { - mrfft.ifft(input, timeBuffer); - int fftSize = mrfft.size(); - const Sample norm = (withScaling ? 1/(Sample)fftSize : 1); - - for (int i = 0; i < offsetSamples; ++i) { - // Inverted polarity since we're using the MRFFT - output[i] = -timeBuffer[i + fftSize - offsetSamples]*norm*(withWindow ? fftWindow[i] : Sample(1)); - } - for (int i = offsetSamples; i < fftSize; ++i) { - output[i] = timeBuffer[i - offsetSamples]*norm*(withWindow ? fftWindow[i] : Sample(1)); - } - } - /// Performs an IFFT (no windowing, scaling or rotation) - template - void ifftRaw(Input &&input, Output &&output) { - mrfft.ifft(input, output); - } - }; - - /** STFT synthesis, built on a `MultiBuffer`. - - Any window length and block interval is supported, but the FFT size may be rounded up to a faster size (by zero-padding). It uses a heuristically-optimal Kaiser window modified for perfect-reconstruction. - - \diagram{stft-aliasing-simulated.svg,Simulated bad-case aliasing (random phase-shift for each band) for overlapping ratios} - - There is a "latest valid index", and you can read the output up to one `historyLength` behind this (see `.resize()`). You can read up to one window-length _ahead_ to get partially-summed future output. - - \diagram{stft-buffer-validity.svg} - - You move the valid index along using `.ensureValid()`, passing in a functor which provides spectra (using `.analyse()` and/or direct modification through `.spectrum[c]`): - - \code - void processSample(...) { - stft.ensureValid([&](int) { - // Here, we introduce (1 - windowSize) of latency - stft.analyse(inputBuffer.view(1 - windowSize)) - }); - // read as a MultiBuffer - auto result = stft.at(0); - ++stft; // also moves the latest valid index - } - - void processBlock(...) { - // assuming `historyLength` == blockSize - stft.ensureValid(blockSize, [&](int blockStartIndex) { - int inputStart = blockStartIndex + (1 - windowSize); - stft.analyse(inputBuffer.view(inputStart)); - }); - auto earliestValid = stft.at(0); - auto latestValid = stft.at(blockSize); - stft += blockSize; - } - \endcode - - The index passed to this functor will be greater than the previous valid index, and `<=` the index you pass in. Therefore, if you call `.ensureValid()` every sample, it can only ever be `0`. - */ - template - class STFT : public signalsmith::delay::MultiBuffer { - - using Super = signalsmith::delay::MultiBuffer; - using Complex = std::complex; - - int channels = 0, _windowSize = 0, _fftSize = 0, _interval = 1; - int validUntilIndex = 0; - - class MultiSpectrum { - int channels, stride; - std::vector buffer; - public: - MultiSpectrum() : MultiSpectrum(0, 0) {} - MultiSpectrum(int channels, int bands) : channels(channels), stride(bands), buffer(channels*bands, 0) {} - - void resize(int nChannels, int nBands) { - channels = nChannels; - stride = nBands; - buffer.assign(channels*stride, 0); - } - - void reset() { - buffer.assign(buffer.size(), 0); - } - - void swap(MultiSpectrum &other) { - using std::swap; - swap(buffer, other.buffer); - } - - Complex * operator [](int channel) { - return buffer.data() + channel*stride; - } - const Complex * operator [](int channel) const { - return buffer.data() + channel*stride; - } - }; - std::vector timeBuffer; - - bool rotate = false; - void resizeInternal(int newChannels, int windowSize, int newInterval, int historyLength, int zeroPadding) { - Super::resize(newChannels, - windowSize /* for output summing */ - + newInterval /* so we can read `windowSize` ahead (we'll be at most `interval-1` from the most recent block */ - + historyLength); - - int fftSize = fft.fastSizeAbove(windowSize + zeroPadding); - - this->channels = newChannels; - _windowSize = windowSize; - this->_fftSize = fftSize; - this->_interval = newInterval; - validUntilIndex = -1; - - setWindow(windowShape, rotate); - - 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}*/ - // TODO: these should both be set before resize() - void setWindow(Window shape, bool rotateToZero=false) { - windowShape = shape; - rotate = rotateToZero; - - 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; - WindowedFFT fft; - - STFT() {} - /// Parameters passed straight to `.resize()` - STFT(int channels, int windowSize, int interval, int historyLength=0, int zeroPadding=0) { - resize(channels, windowSize, interval, historyLength, zeroPadding); - } - - /// Sets the channel-count, FFT size and interval. - void resize(int nChannels, int windowSize, int interval, int historyLength=0, int zeroPadding=0) { - resizeInternal(nChannels, windowSize, interval, historyLength, zeroPadding); - } - - int windowSize() const { - return _windowSize; - } - int fftSize() const { - return _fftSize; - } - int interval() const { - return _interval; - } - /// Returns the (analysis and synthesis) window - decltype(fft.window()) window() const { - return fft.window(); - } - /// Calculates the effective window for the partially-summed future output (relative to the most recent block) - std::vector partialSumWindow(bool includeLatestBlock=true) const { - const auto &w = window(); - std::vector result(_windowSize, 0); - 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; - } - } - return result; - } - - /// Resets everything - since we clear the output sum, it will take `windowSize` samples to get proper output. - void reset() { - Super::reset(); - spectrum.reset(); - validUntilIndex = -1; - } - - /** Generates valid output up to the specified index (or 0), using the callback as many times as needed. - - The callback should be a functor accepting a single integer argument, which is the index for which a spectrum is required. - - The block created from these spectra will start at this index in the output, plus `.latency()`. - */ - template - void ensureValid(int i, AnalysisFn fn) { - while (validUntilIndex < i) { - int blockIndex = validUntilIndex + 1; - fn(blockIndex); - - auto output = this->view(blockIndex); - for (int c = 0; c < channels; ++c) { - auto channel = output[c]; - - // Clear out the future sum, a window-length and an interval ahead - for (int wi = _windowSize; wi < _windowSize + _interval; ++wi) { - channel[wi] = 0; - } - - // Add in the IFFT'd result - fft.ifft(spectrum[c], timeBuffer); - for (int wi = 0; wi < _windowSize; ++wi) { - channel[wi] += timeBuffer[wi]; - } - } - validUntilIndex += _interval; - } - } - /// The same as above, assuming index 0 - template - 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 - - Results can be read/edited using `.spectrum`. */ - template - void analyse(Data &&data) { - for (int c = 0; c < channels; ++c) { - fft.fft(data[c], spectrum[c]); - } - } - template - void analyse(int c, Data &&data) { - fft.fft(data, spectrum[c]); - } - /// Analyse without windowing or zero-rotation - template - void analyseRaw(Data &&data) { - for (int c = 0; c < channels; ++c) { - fft.fftRaw(data[c], spectrum[c]); - } - } - template - void analyseRaw(int c, Data &&data) { - fft.fftRaw(data, spectrum[c]); - } - - int bands() const { - return _fftSize/2; - } - - /** Internal latency (between the block-index requested in `.ensureValid()` and its position in the output) - - Currently unused, but it's in here to allow for a future implementation which spreads the FFT calculations out across each interval.*/ - int latency() { - return 0; - } - - // @name Shift the underlying buffer (moving the "valid" index accordingly) - // @{ - STFT & operator ++() { - Super::operator ++(); - validUntilIndex--; - return *this; - } - STFT & operator +=(int i) { - Super::operator +=(i); - validUntilIndex -= i; - return *this; - } - STFT & operator --() { - Super::operator --(); - validUntilIndex++; - return *this; - } - STFT & operator -=(int i) { - Super::operator -=(i); - validUntilIndex += i; - return *this; - } - // @} - - typename Super::MutableView operator ++(int postIncrement) { - auto result = Super::operator ++(postIncrement); - validUntilIndex--; - return result; - } - typename Super::MutableView operator --(int postIncrement) { - auto result = Super::operator --(postIncrement); - validUntilIndex++; - return result; - } - }; - - /** STFT processing, with input/output. - Before calling `.ensureValid(index)`, you should make sure the input is filled up to `index`. - */ - template - class ProcessSTFT : public STFT { - using Super = STFT; - public: - signalsmith::delay::MultiBuffer input; - - ProcessSTFT(int inChannels, int outChannels, int windowSize, int interval, int historyLength=0) { - resize(inChannels, outChannels, windowSize, interval, historyLength); - } - - /** Alter the spectrum, using input up to this point, for the output block starting from this point. - Sub-classes should replace this with whatever processing is desired. */ - virtual void processSpectrum(int /*blockIndex*/) {} - - /// Sets the input/output channels, FFT size and interval. - void resize(int inChannels, int outChannels, int windowSize, int interval, int historyLength=0) { - Super::resize(outChannels, windowSize, interval, historyLength); - input.resize(inChannels, windowSize + interval + historyLength); - } - void reset(Sample value=Sample()) { - Super::reset(value); - input.reset(value); - } - - /// Internal latency, including buffering samples for analysis. - int latency() { - return Super::latency() + (this->windowSize() - 1); - } - - void ensureValid(int i=0) { - Super::ensureValid(i, [&](int blockIndex) { - this->analyse(input.view(blockIndex - this->windowSize() + 1)); - this->processSpectrum(blockIndex); - }); - } - - // @name Shift the output, input, and valid index. - // @{ - ProcessSTFT & operator ++() { - Super::operator ++(); - ++input; - return *this; - } - ProcessSTFT & operator +=(int i) { - Super::operator +=(i); - input += i; - return *this; - } - ProcessSTFT & operator --() { - Super::operator --(); - --input; - return *this; - } - ProcessSTFT & operator -=(int i) { - Super::operator -=(i); - input -= i; - return *this; - } - // @} - }; - -/** @} */ -}} // signalsmith::spectral:: -#endif // include guard diff --git a/dsp/windows.h b/dsp/windows.h deleted file mode 100644 index 59f2f33..0000000 --- a/dsp/windows.h +++ /dev/null @@ -1,219 +0,0 @@ -#include "./common.h" - -#ifndef SIGNALSMITH_DSP_WINDOWS_H -#define SIGNALSMITH_DSP_WINDOWS_H - -#include -#include - -namespace signalsmith { -namespace windows { - /** @defgroup Windows Window functions - @brief Windows for spectral analysis - - These are generally double-precision, because they are mostly calculated during setup/reconfiguring, not real-time code. - - @{ - @file - */ - - /** @brief The Kaiser window (almost) maximises the energy in the main-lobe compared to the side-lobes. - - Kaiser windows can be constructing using the shape-parameter (beta) or using the static `with???()` methods.*/ - class Kaiser { - // I_0(x)=\sum_{k=0}^{N}\frac{x^{2k}}{(k!)^2\cdot4^k} - inline static double bessel0(double x) { - const double significanceLimit = 1e-4; - double result = 0; - double term = 1; - double m = 0; - while (term > significanceLimit) { - result += term; - ++m; - term *= (x*x)/(4*m*m); - } - - return result; - } - double beta; - double invB0; - - static double heuristicBandwidth(double bandwidth) { - // Good peaks - //return bandwidth + 8/((bandwidth + 3)*(bandwidth + 3)); - // Good average - //return bandwidth + 14/((bandwidth + 2.5)*(bandwidth + 2.5)); - // Compromise - return bandwidth + 8/((bandwidth + 3)*(bandwidth + 3)) + 0.25*std::max(3 - bandwidth, 0.0); - } - public: - /// Set up a Kaiser window with a given shape. `beta` is `pi*alpha` (since there is ambiguity about shape parameters) - Kaiser(double beta) : beta(beta), invB0(1/bessel0(beta)) {} - - /// @name Bandwidth methods - /// @{ - static Kaiser withBandwidth(double bandwidth, bool heuristicOptimal=false) { - return Kaiser(bandwidthToBeta(bandwidth, heuristicOptimal)); - } - - /** Returns the Kaiser shape where the main lobe has the specified bandwidth (as a factor of 1/window-length). - \diagram{kaiser-windows.svg,You can see that the main lobe matches the specified bandwidth.} - If `heuristicOptimal` is enabled, the main lobe width is _slightly_ wider, improving both the peak and total energy - see `bandwidthToEnergyDb()` and `bandwidthToPeakDb()`. - \diagram{kaiser-windows-heuristic.svg, The main lobe extends to ±bandwidth/2.} */ - static double bandwidthToBeta(double bandwidth, bool heuristicOptimal=false) { - if (heuristicOptimal) { // Heuristic based on numerical search - bandwidth = heuristicBandwidth(bandwidth); - } - bandwidth = std::max(bandwidth, 2.0); - double alpha = std::sqrt(bandwidth*bandwidth*0.25 - 1); - return alpha*M_PI; - } - - static double betaToBandwidth(double beta) { - double alpha = beta*(1.0/M_PI); - return 2*std::sqrt(alpha*alpha + 1); - } - /// @} - - /// @name Performance methods - /// @{ - /** @brief Total energy ratio (in dB) between side-lobes and the main lobe. - \diagram{windows-kaiser-sidelobe-energy.svg,Measured main/side lobe energy ratio. You can see that the heuristic improves performance for all bandwidth values.} - This function uses an approximation which is accurate to ±0.5dB for 2 ⩽ bandwidth ≤ 10, or 1 ⩽ bandwidth ≤ 10 when `heuristicOptimal`is enabled. - */ - static double bandwidthToEnergyDb(double bandwidth, bool heuristicOptimal=false) { - // Horrible heuristic fits - if (heuristicOptimal) { - if (bandwidth < 3) bandwidth += (3 - bandwidth)*0.5; - return 12.9 + -3/(bandwidth + 0.4) - 13.4*bandwidth + (bandwidth < 3)*-9.6*(bandwidth - 3); - } - return 10.5 + 15/(bandwidth + 0.4) - 13.25*bandwidth + (bandwidth < 2)*13*(bandwidth - 2); - } - static double energyDbToBandwidth(double energyDb, bool heuristicOptimal=false) { - double bw = 1; - while (bw < 20 && bandwidthToEnergyDb(bw, heuristicOptimal) > energyDb) { - bw *= 2; - } - double step = bw/2; - while (step > 0.0001) { - if (bandwidthToEnergyDb(bw, heuristicOptimal) > energyDb) { - bw += step; - } else { - bw -= step; - } - step *= 0.5; - } - return bw; - } - /** @brief Peak ratio (in dB) between side-lobes and the main lobe. - \diagram{windows-kaiser-sidelobe-peaks.svg,Measured main/side lobe peak ratio. You can see that the heuristic improves performance, except in the bandwidth range 1-2 where peak ratio was sacrificed to improve total energy ratio.} - This function uses an approximation which is accurate to ±0.5dB for 2 ⩽ bandwidth ≤ 9, or 0.5 ⩽ bandwidth ≤ 9 when `heuristicOptimal`is enabled. - */ - static double bandwidthToPeakDb(double bandwidth, bool heuristicOptimal=false) { - // Horrible heuristic fits - if (heuristicOptimal) { - return 14.2 - 20/(bandwidth + 1) - 13*bandwidth + (bandwidth < 3)*-6*(bandwidth - 3) + (bandwidth < 2.25)*5.8*(bandwidth - 2.25); - } - return 10 + 8/(bandwidth + 2) - 12.75*bandwidth + (bandwidth < 2)*4*(bandwidth - 2); - } - static double peakDbToBandwidth(double peakDb, bool heuristicOptimal=false) { - double bw = 1; - while (bw < 20 && bandwidthToPeakDb(bw, heuristicOptimal) > peakDb) { - bw *= 2; - } - double step = bw/2; - while (step > 0.0001) { - if (bandwidthToPeakDb(bw, heuristicOptimal) > peakDb) { - bw += step; - } else { - bw -= step; - } - step *= 0.5; - } - return bw; - } - /** @} */ - - /** Equivalent noise bandwidth (ENBW), a measure of frequency resolution. - \diagram{windows-kaiser-enbw.svg,Measured ENBW\, with and without the heuristic bandwidth adjustment.} - This approximation is accurate to ±0.05 up to a bandwidth of 22. - */ - static double bandwidthToEnbw(double bandwidth, bool heuristicOptimal=false) { - if (heuristicOptimal) bandwidth = heuristicBandwidth(bandwidth); - double b2 = std::max(bandwidth - 2, 0); - return 1 + b2*(0.2 + b2*(-0.005 + b2*(-0.000005 + b2*0.0000022))); - } - - /// Return the window's value for position in the range [0, 1] - double operator ()(double unit) { - double r = 2*unit - 1; - double arg = std::sqrt(1 - r*r); - return bessel0(beta*arg)*invB0; - } - - /// Fills an arbitrary container with a Kaiser window - template - void fill(Data &&data, int size) const { - double invSize = 1.0/size; - for (int i = 0; i < size; ++i) { - double r = (2*i + 1)*invSize - 1; - double arg = std::sqrt(1 - r*r); - data[i] = bessel0(beta*arg)*invB0; - } - } - }; - - /** @brief The Approximate Confined Gaussian window is (almost) optimal - - ACG windows can be constructing using the shape-parameter (sigma) or using the static `with???()` methods.*/ - class ApproximateConfinedGaussian { - double gaussianFactor; - - double gaussian(double x) const { - return std::exp(-x*x*gaussianFactor); - } - public: - /// Heuristic map from bandwidth to the appropriately-optimal sigma - static double bandwidthToSigma(double bandwidth) { - return 0.3/std::sqrt(bandwidth); - } - static ApproximateConfinedGaussian withBandwidth(double bandwidth) { - return ApproximateConfinedGaussian(bandwidthToSigma(bandwidth)); - } - - ApproximateConfinedGaussian(double sigma) : gaussianFactor(0.0625/(sigma*sigma)) {} - - /// Fills an arbitrary container - template - void fill(Data &&data, int size) const { - double invSize = 1.0/size; - double offsetScale = gaussian(1)/(gaussian(3) + gaussian(-1)); - double norm = 1/(gaussian(0) - 2*offsetScale*(gaussian(2))); - for (int i = 0; i < size; ++i) { - double r = (2*i + 1)*invSize - 1; - data[i] = norm*(gaussian(r) - offsetScale*(gaussian(r - 2) + gaussian(r + 2))); - } - } - }; - - /** Forces STFT perfect-reconstruction (WOLA) on an existing window, for a given STFT interval. - For example, here are perfect-reconstruction versions of the approximately-optimal @ref Kaiser windows: - \diagram{kaiser-windows-heuristic-pr.svg,Note the lower overall energy\, and the pointy top for 2x bandwidth. Spectral performance is about the same\, though.} - */ - template - void forcePerfectReconstruction(Data &&data, int windowLength, int interval) { - for (int i = 0; i < interval; ++i) { - double sum2 = 0; - for (int index = i; index < windowLength; index += interval) { - sum2 += data[index]*data[index]; - } - double factor = 1/std::sqrt(sum2); - for (int index = i; index < windowLength; index += interval) { - data[index] *= factor; - } - } - } - -/** @} */ -}} // signalsmith::windows -#endif // include guard diff --git a/include/signalsmith-stretch/signalsmith-stretch.h b/include/signalsmith-stretch/signalsmith-stretch.h new file mode 100644 index 0000000..d897b5e --- /dev/null +++ b/include/signalsmith-stretch/signalsmith-stretch.h @@ -0,0 +1 @@ +#include "../../signalsmith-stretch.h" \ No newline at end of file diff --git a/signalsmith-stretch.h b/signalsmith-stretch.h index 8ce65b5..72a2318 100644 --- a/signalsmith-stretch.h +++ b/signalsmith-stretch.h @@ -3,10 +3,10 @@ //#include "dsp/spectral.h" //#include "dsp/delay.h" -#include "dsp/perf.h" -SIGNALSMITH_DSP_VERSION_CHECK(1, 6, 0); // Check version is compatible +//#include "dsp/perf.h" +//SIGNALSMITH_DSP_VERSION_CHECK(1, 6, 0); // Check version is compatible -#include "linear/stft.h" +#include "signalsmith-linear/stft.h" // https://github.com/Signalsmith-Audio/linear #include #include #include @@ -14,6 +14,19 @@ SIGNALSMITH_DSP_VERSION_CHECK(1, 6, 0); // Check version is compatible namespace signalsmith { namespace stretch { +namespace _impl { + template + static std::complex mul(const std::complex &a, const std::complex &b) { + return conjugateSecond ? std::complex{ + b.real()*a.real() + b.imag()*a.imag(), + b.real()*a.imag() - b.imag()*a.real() + } : std::complex{ + a.real()*b.real() - a.imag()*b.imag(), + a.real()*b.imag() + a.imag()*b.real() + }; + } +} + template struct SignalsmithStretch { static constexpr size_t version[3] = {1, 1, 1}; @@ -385,9 +398,9 @@ private: for (int b = 0; b < bands; ++b) { auto &bin = bins[b]; - bin.output = signalsmith::perf::mul(bin.output, rot); - bin.prevInput = signalsmith::perf::mul(bin.prevInput, rot); - rot = signalsmith::perf::mul(rot, rotStep); + bin.output = _impl::mul(bin.output, rot); + bin.prevInput = _impl::mul(bin.prevInput, rot); + rot = _impl::mul(rot, rotStep); } } } @@ -426,8 +439,8 @@ private: auto &outputBin = bins[b]; Complex prevInput = getFractional<&Band::prevInput>(c, lowIndex, fracIndex); - Complex freqTwist = signalsmith::perf::mul(prediction.input, prevInput); - Complex phase = signalsmith::perf::mul(outputBin.output, freqTwist); + Complex freqTwist = _impl::mul(prediction.input, prevInput); + Complex phase = _impl::mul(outputBin.output, freqTwist); outputBin.output = phase/(std::max(prevEnergy, prediction.energy) + noiseFloor); } } @@ -457,17 +470,17 @@ private: if (b > 0) { Sample binTimeFactor = randomTimeFactor ? timeFactorDist(randomEngine) : timeFactor; Complex downInput = getFractional<&Band::input>(maxChannel, mapPoint.inputBin - binTimeFactor); - Complex shortVerticalTwist = signalsmith::perf::mul(prediction.input, downInput); + Complex shortVerticalTwist = _impl::mul(prediction.input, downInput); auto &downBin = bins[b - 1]; - phase += signalsmith::perf::mul(downBin.output, shortVerticalTwist); + phase += _impl::mul(downBin.output, shortVerticalTwist); if (b >= longVerticalStep) { Complex longDownInput = getFractional<&Band::input>(maxChannel, mapPoint.inputBin - longVerticalStep*binTimeFactor); - Complex longVerticalTwist = signalsmith::perf::mul(prediction.input, longDownInput); + Complex longVerticalTwist = _impl::mul(prediction.input, longDownInput); auto &longDownBin = bins[b - longVerticalStep]; - phase += signalsmith::perf::mul(longDownBin.output, longVerticalTwist); + phase += _impl::mul(longDownBin.output, longVerticalTwist); } } // Downwards vertical steps @@ -477,20 +490,20 @@ private: Sample binTimeFactor = randomTimeFactor ? timeFactorDist(randomEngine) : timeFactor; Complex downInput = getFractional<&Band::input>(maxChannel, upMapPoint.inputBin - binTimeFactor); - Complex shortVerticalTwist = signalsmith::perf::mul(upPrediction.input, downInput); + Complex shortVerticalTwist = _impl::mul(upPrediction.input, downInput); auto &upBin = bins[b + 1]; - phase += signalsmith::perf::mul(upBin.output, shortVerticalTwist); + phase += _impl::mul(upBin.output, shortVerticalTwist); if (b < bands - longVerticalStep) { auto &longUpPrediction = predictions[b + longVerticalStep]; auto &longUpMapPoint = outputMap[b + longVerticalStep]; Complex longDownInput = getFractional<&Band::input>(maxChannel, longUpMapPoint.inputBin - longVerticalStep*binTimeFactor); - Complex longVerticalTwist = signalsmith::perf::mul(longUpPrediction.input, longDownInput); + Complex longVerticalTwist = _impl::mul(longUpPrediction.input, longDownInput); auto &longUpBin = bins[b + longVerticalStep]; - phase += signalsmith::perf::mul(longUpBin.output, longVerticalTwist); + phase += _impl::mul(longUpBin.output, longVerticalTwist); } } @@ -502,8 +515,8 @@ private: auto &channelBin = bandsForChannel(c)[b]; auto &channelPrediction = predictionsForChannel(c)[b]; - Complex channelTwist = signalsmith::perf::mul(channelPrediction.input, prediction.input); - Complex channelPhase = signalsmith::perf::mul(outputBin.output, channelTwist); + Complex channelTwist = _impl::mul(channelPrediction.input, prediction.input); + Complex channelPhase = _impl::mul(outputBin.output, channelTwist); channelBin.output = channelPrediction.makeOutput(channelPhase); } }