Move the filters to core

6 files changed, 902 insertions, 0 deletions

diff --git a/core/filters/biquad.cpp b/core/filters/biquad.cpp
new file mode 100644
index 00000000..fefdc8e1
--- /dev/null
+++ b/core/filters/biquad.cpp
@@ -0,0 +1,163 @@
+
+#include "config.h"
+
+#include "biquad.h"
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+
+#include "opthelpers.h"
+
+
+template<typename Real>
+void BiquadFilterR<Real>::setParams(BiquadType type, Real f0norm, Real gain, Real rcpQ)
+{
+    // Limit gain to -100dB
+    assert(gain > 0.00001f);
+
+    const Real w0{al::MathDefs<Real>::Tau() * f0norm};
+    const Real sin_w0{std::sin(w0)};
+    const Real cos_w0{std::cos(w0)};
+    const Real alpha{sin_w0/2.0f * rcpQ};
+
+    Real sqrtgain_alpha_2;
+    Real a[3]{ 1.0f, 0.0f, 0.0f };
+    Real b[3]{ 1.0f, 0.0f, 0.0f };
+
+    /* Calculate filter coefficients depending on filter type */
+    switch(type)
+    {
+        case BiquadType::HighShelf:
+            sqrtgain_alpha_2 = 2.0f * std::sqrt(gain) * alpha;
+            b[0] =       gain*((gain+1.0f) + (gain-1.0f)*cos_w0 + sqrtgain_alpha_2);
+            b[1] = -2.0f*gain*((gain-1.0f) + (gain+1.0f)*cos_w0                   );
+            b[2] =       gain*((gain+1.0f) + (gain-1.0f)*cos_w0 - sqrtgain_alpha_2);
+            a[0] =             (gain+1.0f) - (gain-1.0f)*cos_w0 + sqrtgain_alpha_2;
+            a[1] =  2.0f*     ((gain-1.0f) - (gain+1.0f)*cos_w0                   );
+            a[2] =             (gain+1.0f) - (gain-1.0f)*cos_w0 - sqrtgain_alpha_2;
+            break;
+        case BiquadType::LowShelf:
+            sqrtgain_alpha_2 = 2.0f * std::sqrt(gain) * alpha;
+            b[0] =       gain*((gain+1.0f) - (gain-1.0f)*cos_w0 + sqrtgain_alpha_2);
+            b[1] =  2.0f*gain*((gain-1.0f) - (gain+1.0f)*cos_w0                   );
+            b[2] =       gain*((gain+1.0f) - (gain-1.0f)*cos_w0 - sqrtgain_alpha_2);
+            a[0] =             (gain+1.0f) + (gain-1.0f)*cos_w0 + sqrtgain_alpha_2;
+            a[1] = -2.0f*     ((gain-1.0f) + (gain+1.0f)*cos_w0                   );
+            a[2] =             (gain+1.0f) + (gain-1.0f)*cos_w0 - sqrtgain_alpha_2;
+            break;
+        case BiquadType::Peaking:
+            b[0] =  1.0f + alpha * gain;
+            b[1] = -2.0f * cos_w0;
+            b[2] =  1.0f - alpha * gain;
+            a[0] =  1.0f + alpha / gain;
+            a[1] = -2.0f * cos_w0;
+            a[2] =  1.0f - alpha / gain;
+            break;
+
+        case BiquadType::LowPass:
+            b[0] = (1.0f - cos_w0) / 2.0f;
+            b[1] =  1.0f - cos_w0;
+            b[2] = (1.0f - cos_w0) / 2.0f;
+            a[0] =  1.0f + alpha;
+            a[1] = -2.0f * cos_w0;
+            a[2] =  1.0f - alpha;
+            break;
+        case BiquadType::HighPass:
+            b[0] =  (1.0f + cos_w0) / 2.0f;
+            b[1] = -(1.0f + cos_w0);
+            b[2] =  (1.0f + cos_w0) / 2.0f;
+            a[0] =   1.0f + alpha;
+            a[1] =  -2.0f * cos_w0;
+            a[2] =   1.0f - alpha;
+            break;
+        case BiquadType::BandPass:
+            b[0] =  alpha;
+            b[1] =  0.0f;
+            b[2] = -alpha;
+            a[0] =  1.0f + alpha;
+            a[1] = -2.0f * cos_w0;
+            a[2] =  1.0f - alpha;
+            break;
+    }
+
+    mA1 = a[1] / a[0];
+    mA2 = a[2] / a[0];
+    mB0 = b[0] / a[0];
+    mB1 = b[1] / a[0];
+    mB2 = b[2] / a[0];
+}
+
+template<typename Real>
+void BiquadFilterR<Real>::process(const al::span<const Real> src, Real *dst)
+{
+    const Real b0{mB0};
+    const Real b1{mB1};
+    const Real b2{mB2};
+    const Real a1{mA1};
+    const Real a2{mA2};
+    Real z1{mZ1};
+    Real z2{mZ2};
+
+    /* Processing loop is Transposed Direct Form II. This requires less storage
+     * compared to Direct Form I (only two delay components, instead of a four-
+     * sample history; the last two inputs and outputs), and works better for
+     * floating-point which favors summing similarly-sized values while being
+     * less bothered by overflow.
+     *
+     * See: http://www.earlevel.com/main/2003/02/28/biquads/
+     */
+    auto proc_sample = [b0,b1,b2,a1,a2,&z1,&z2](Real input) noexcept -> Real
+    {
+        const Real output{input*b0 + z1};
+        z1 = input*b1 - output*a1 + z2;
+        z2 = input*b2 - output*a2;
+        return output;
+    };
+    std::transform(src.cbegin(), src.cend(), dst, proc_sample);
+
+    mZ1 = z1;
+    mZ2 = z2;
+}
+
+template<typename Real>
+void BiquadFilterR<Real>::dualProcess(BiquadFilterR &other, const al::span<const Real> src,
+    Real *dst)
+{
+    const Real b00{mB0};
+    const Real b01{mB1};
+    const Real b02{mB2};
+    const Real a01{mA1};
+    const Real a02{mA2};
+    const Real b10{other.mB0};
+    const Real b11{other.mB1};
+    const Real b12{other.mB2};
+    const Real a11{other.mA1};
+    const Real a12{other.mA2};
+    Real z01{mZ1};
+    Real z02{mZ2};
+    Real z11{other.mZ1};
+    Real z12{other.mZ2};
+
+    auto proc_sample = [b00,b01,b02,a01,a02,b10,b11,b12,a11,a12,&z01,&z02,&z11,&z12](Real input) noexcept -> Real
+    {
+        const Real tmpout{input*b00 + z01};
+        z01 = input*b01 - tmpout*a01 + z02;
+        z02 = input*b02 - tmpout*a02;
+        input = tmpout;
+
+        const Real output{input*b10 + z11};
+        z11 = input*b11 - output*a11 + z12;
+        z12 = input*b12 - output*a12;
+        return output;
+    };
+    std::transform(src.cbegin(), src.cend(), dst, proc_sample);
+
+    mZ1 = z01;
+    mZ2 = z02;
+    other.mZ1 = z11;
+    other.mZ2 = z12;
+}
+
+template class BiquadFilterR<float>;
+template class BiquadFilterR<double>;
diff --git a/core/filters/biquad.h b/core/filters/biquad.h
new file mode 100644
index 00000000..b2e2cfdb
--- /dev/null
+++ b/core/filters/biquad.h
@@ -0,0 +1,144 @@
+#ifndef CORE_FILTERS_BIQUAD_H
+#define CORE_FILTERS_BIQUAD_H
+
+#include <algorithm>
+#include <cmath>
+#include <cstddef>
+#include <utility>
+
+#include "alspan.h"
+#include "math_defs.h"
+
+
+/* Filters implementation is based on the "Cookbook formulae for audio
+ * EQ biquad filter coefficients" by Robert Bristow-Johnson
+ * http://www.musicdsp.org/files/Audio-EQ-Cookbook.txt
+ */
+/* Implementation note: For the shelf and peaking filters, the specified gain
+ * is for the centerpoint of the transition band. This better fits EFX filter
+ * behavior, which expects the shelf's reference frequency to reach the given
+ * gain. To set the gain for the shelf or peak itself, use the square root of
+ * the desired linear gain (or halve the dB gain).
+ */
+
+enum class BiquadType {
+    /** EFX-style low-pass filter, specifying a gain and reference frequency. */
+    HighShelf,
+    /** EFX-style high-pass filter, specifying a gain and reference frequency. */
+    LowShelf,
+    /** Peaking filter, specifying a gain and reference frequency. */
+    Peaking,
+
+    /** Low-pass cut-off filter, specifying a cut-off frequency. */
+    LowPass,
+    /** High-pass cut-off filter, specifying a cut-off frequency. */
+    HighPass,
+    /** Band-pass filter, specifying a center frequency. */
+    BandPass,
+};
+
+template<typename Real>
+class BiquadFilterR {
+    /* Last two delayed components for direct form II. */
+    Real mZ1{0.0f}, mZ2{0.0f};
+    /* Transfer function coefficients "b" (numerator) */
+    Real mB0{1.0f}, mB1{0.0f}, mB2{0.0f};
+    /* Transfer function coefficients "a" (denominator; a0 is pre-applied). */
+    Real mA1{0.0f}, mA2{0.0f};
+
+    void setParams(BiquadType type, Real f0norm, Real gain, Real rcpQ);
+
+    /**
+     * Calculates the rcpQ (i.e. 1/Q) coefficient for shelving filters, using
+     * the reference gain and shelf slope parameter.
+     * \param gain 0 < gain
+     * \param slope 0 < slope <= 1
+     */
+    static Real rcpQFromSlope(Real gain, Real slope)
+    { return std::sqrt((gain + 1.0f/gain)*(1.0f/slope - 1.0f) + 2.0f); }
+
+    /**
+     * Calculates the rcpQ (i.e. 1/Q) coefficient for filters, using the
+     * normalized reference frequency and bandwidth.
+     * \param f0norm 0 < f0norm < 0.5.
+     * \param bandwidth 0 < bandwidth
+     */
+    static Real rcpQFromBandwidth(Real f0norm, Real bandwidth)
+    {
+        const Real w0{al::MathDefs<Real>::Tau() * f0norm};
+        return 2.0f*std::sinh(std::log(Real{2.0f})/2.0f*bandwidth*w0/std::sin(w0));
+    }
+
+public:
+    void clear() noexcept { mZ1 = mZ2 = 0.0f; }
+
+    /**
+     * Sets the filter state for the specified filter type and its parameters.
+     *
+     * \param type The type of filter to apply.
+     * \param f0norm The normalized reference frequency (ref / sample_rate).
+     * This is the center point for the Shelf, Peaking, and BandPass filter
+     * types, or the cutoff frequency for the LowPass and HighPass filter
+     * types.
+     * \param gain The gain for the reference frequency response. Only used by
+     * the Shelf and Peaking filter types.
+     * \param slope Slope steepness of the transition band.
+     */
+    void setParamsFromSlope(BiquadType type, Real f0norm, Real gain, Real slope)
+    {
+        gain = std::max<Real>(gain, 0.001f); /* Limit -60dB */
+        setParams(type, f0norm, gain, rcpQFromSlope(gain, slope));
+    }
+
+    /**
+     * Sets the filter state for the specified filter type and its parameters.
+     *
+     * \param type The type of filter to apply.
+     * \param f0norm The normalized reference frequency (ref / sample_rate).
+     * This is the center point for the Shelf, Peaking, and BandPass filter
+     * types, or the cutoff frequency for the LowPass and HighPass filter
+     * types.
+     * \param gain The gain for the reference frequency response. Only used by
+     * the Shelf and Peaking filter types.
+     * \param bandwidth Normalized bandwidth of the transition band.
+     */
+    void setParamsFromBandwidth(BiquadType type, Real f0norm, Real gain, Real bandwidth)
+    { setParams(type, f0norm, gain, rcpQFromBandwidth(f0norm, bandwidth)); }
+
+    void copyParamsFrom(const BiquadFilterR &other)
+    {
+        mB0 = other.mB0;
+        mB1 = other.mB1;
+        mB2 = other.mB2;
+        mA1 = other.mA1;
+        mA2 = other.mA2;
+    }
+
+    void process(const al::span<const Real> src, Real *dst);
+    /** Processes this filter and the other at the same time. */
+    void dualProcess(BiquadFilterR &other, const al::span<const Real> src, Real *dst);
+
+    /* Rather hacky. It's just here to support "manual" processing. */
+    std::pair<Real,Real> getComponents() const noexcept { return {mZ1, mZ2}; }
+    void setComponents(Real z1, Real z2) noexcept { mZ1 = z1; mZ2 = z2; }
+    Real processOne(const Real in, Real &z1, Real &z2) const noexcept
+    {
+        const Real out{in*mB0 + z1};
+        z1 = in*mB1 - out*mA1 + z2;
+        z2 = in*mB2 - out*mA2;
+        return out;
+    }
+};
+
+template<typename Real>
+struct DualBiquadR {
+    BiquadFilterR<Real> &f0, &f1;
+
+    void process(const al::span<const Real> src, Real *dst)
+    { f0.dualProcess(f1, src, dst); }
+};
+
+using BiquadFilter = BiquadFilterR<float>;
+using DualBiquad = DualBiquadR<float>;
+
+#endif /* CORE_FILTERS_BIQUAD_H */
diff --git a/core/filters/nfc.cpp b/core/filters/nfc.cpp
new file mode 100644
index 00000000..9a28517c
--- /dev/null
+++ b/core/filters/nfc.cpp
@@ -0,0 +1,383 @@
+
+#include "config.h"
+
+#include "nfc.h"
+
+#include <algorithm>
+
+#include "opthelpers.h"
+
+
+/* Near-field control filters are the basis for handling the near-field effect.
+ * The near-field effect is a bass-boost present in the directional components
+ * of a recorded signal, created as a result of the wavefront curvature (itself
+ * a function of sound distance). Proper reproduction dictates this be
+ * compensated for using a bass-cut given the playback speaker distance, to
+ * avoid excessive bass in the playback.
+ *
+ * For real-time rendered audio, emulating the near-field effect based on the
+ * sound source's distance, and subsequently compensating for it at output
+ * based on the speaker distances, can create a more realistic perception of
+ * sound distance beyond a simple 1/r attenuation.
+ *
+ * These filters do just that. Each one applies a low-shelf filter, created as
+ * the combination of a bass-boost for a given sound source distance (near-
+ * field emulation) along with a bass-cut for a given control/speaker distance
+ * (near-field compensation).
+ *
+ * Note that it is necessary to apply a cut along with the boost, since the
+ * boost alone is unstable in higher-order ambisonics as it causes an infinite
+ * DC gain (even first-order ambisonics requires there to be no DC offset for
+ * the boost to work). Consequently, ambisonics requires a control parameter to
+ * be used to avoid an unstable boost-only filter. NFC-HOA defines this control
+ * as a reference delay, calculated with:
+ *
+ * reference_delay = control_distance / speed_of_sound
+ *
+ * This means w0 (for input) or w1 (for output) should be set to:
+ *
+ * wN = 1 / (reference_delay * sample_rate)
+ *
+ * when dealing with NFC-HOA content. For FOA input content, which does not
+ * specify a reference_delay variable, w0 should be set to 0 to apply only
+ * near-field compensation for output. It's important that w1 be a finite,
+ * positive, non-0 value or else the bass-boost will become unstable again.
+ * Also, w0 should not be too large compared to w1, to avoid excessively loud
+ * low frequencies.
+ */
+
+namespace {
+
+constexpr float B[5][4] = {
+    {    0.0f                             },
+    {    1.0f                             },
+    {    3.0f,     3.0f                   },
+    { 3.6778f,  6.4595f, 2.3222f          },
+    { 4.2076f, 11.4877f, 5.7924f, 9.1401f }
+};
+
+NfcFilter1 NfcFilterCreate1(const float w0, const float w1) noexcept
+{
+    NfcFilter1 nfc{};
+    float b_00, g_0;
+    float r;
+
+    nfc.base_gain = 1.0f;
+    nfc.gain = 1.0f;
+
+    /* Calculate bass-boost coefficients. */
+    r = 0.5f * w0;
+    b_00 = B[1][0] * r;
+    g_0 = 1.0f + b_00;
+
+    nfc.gain *= g_0;
+    nfc.b1 = 2.0f * b_00 / g_0;
+
+    /* Calculate bass-cut coefficients. */
+    r = 0.5f * w1;
+    b_00 = B[1][0] * r;
+    g_0 = 1.0f + b_00;
+
+    nfc.base_gain /= g_0;
+    nfc.gain /= g_0;
+    nfc.a1 = 2.0f * b_00 / g_0;
+
+    return nfc;
+}
+
+void NfcFilterAdjust1(NfcFilter1 *nfc, const float w0) noexcept
+{
+    const float r{0.5f * w0};
+    const float b_00{B[1][0] * r};
+    const float g_0{1.0f + b_00};
+
+    nfc->gain = nfc->base_gain * g_0;
+    nfc->b1 = 2.0f * b_00 / g_0;
+}
+
+
+NfcFilter2 NfcFilterCreate2(const float w0, const float w1) noexcept
+{
+    NfcFilter2 nfc{};
+    float b_10, b_11, g_1;
+    float r;
+
+    nfc.base_gain = 1.0f;
+    nfc.gain = 1.0f;
+
+    /* Calculate bass-boost coefficients. */
+    r = 0.5f * w0;
+    b_10 = B[2][0] * r;
+    b_11 = B[2][1] * r * r;
+    g_1 = 1.0f + b_10 + b_11;
+
+    nfc.gain *= g_1;
+    nfc.b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc.b2 = 4.0f * b_11 / g_1;
+
+    /* Calculate bass-cut coefficients. */
+    r = 0.5f * w1;
+    b_10 = B[2][0] * r;
+    b_11 = B[2][1] * r * r;
+    g_1 = 1.0f + b_10 + b_11;
+
+    nfc.base_gain /= g_1;
+    nfc.gain /= g_1;
+    nfc.a1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc.a2 = 4.0f * b_11 / g_1;
+
+    return nfc;
+}
+
+void NfcFilterAdjust2(NfcFilter2 *nfc, const float w0) noexcept
+{
+    const float r{0.5f * w0};
+    const float b_10{B[2][0] * r};
+    const float b_11{B[2][1] * r * r};
+    const float g_1{1.0f + b_10 + b_11};
+
+    nfc->gain = nfc->base_gain * g_1;
+    nfc->b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc->b2 = 4.0f * b_11 / g_1;
+}
+
+
+NfcFilter3 NfcFilterCreate3(const float w0, const float w1) noexcept
+{
+    NfcFilter3 nfc{};
+    float b_10, b_11, g_1;
+    float b_00, g_0;
+    float r;
+
+    nfc.base_gain = 1.0f;
+    nfc.gain = 1.0f;
+
+    /* Calculate bass-boost coefficients. */
+    r = 0.5f * w0;
+    b_10 = B[3][0] * r;
+    b_11 = B[3][1] * r * r;
+    b_00 = B[3][2] * r;
+    g_1 = 1.0f + b_10 + b_11;
+    g_0 = 1.0f + b_00;
+
+    nfc.gain *= g_1 * g_0;
+    nfc.b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc.b2 = 4.0f * b_11 / g_1;
+    nfc.b3 = 2.0f * b_00 / g_0;
+
+    /* Calculate bass-cut coefficients. */
+    r = 0.5f * w1;
+    b_10 = B[3][0] * r;
+    b_11 = B[3][1] * r * r;
+    b_00 = B[3][2] * r;
+    g_1 = 1.0f + b_10 + b_11;
+    g_0 = 1.0f + b_00;
+
+    nfc.base_gain /= g_1 * g_0;
+    nfc.gain /= g_1 * g_0;
+    nfc.a1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc.a2 = 4.0f * b_11 / g_1;
+    nfc.a3 = 2.0f * b_00 / g_0;
+
+    return nfc;
+}
+
+void NfcFilterAdjust3(NfcFilter3 *nfc, const float w0) noexcept
+{
+    const float r{0.5f * w0};
+    const float b_10{B[3][0] * r};
+    const float b_11{B[3][1] * r * r};
+    const float b_00{B[3][2] * r};
+    const float g_1{1.0f + b_10 + b_11};
+    const float g_0{1.0f + b_00};
+
+    nfc->gain = nfc->base_gain * g_1 * g_0;
+    nfc->b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc->b2 = 4.0f * b_11 / g_1;
+    nfc->b3 = 2.0f * b_00 / g_0;
+}
+
+
+NfcFilter4 NfcFilterCreate4(const float w0, const float w1) noexcept
+{
+    NfcFilter4 nfc{};
+    float b_10, b_11, g_1;
+    float b_00, b_01, g_0;
+    float r;
+
+    nfc.base_gain = 1.0f;
+    nfc.gain = 1.0f;
+
+    /* Calculate bass-boost coefficients. */
+    r = 0.5f * w0;
+    b_10 = B[4][0] * r;
+    b_11 = B[4][1] * r * r;
+    b_00 = B[4][2] * r;
+    b_01 = B[4][3] * r * r;
+    g_1 = 1.0f + b_10 + b_11;
+    g_0 = 1.0f + b_00 + b_01;
+
+    nfc.gain *= g_1 * g_0;
+    nfc.b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc.b2 = 4.0f * b_11 / g_1;
+    nfc.b3 = (2.0f*b_00 + 4.0f*b_01) / g_0;
+    nfc.b4 = 4.0f * b_01 / g_0;
+
+    /* Calculate bass-cut coefficients. */
+    r = 0.5f * w1;
+    b_10 = B[4][0] * r;
+    b_11 = B[4][1] * r * r;
+    b_00 = B[4][2] * r;
+    b_01 = B[4][3] * r * r;
+    g_1 = 1.0f + b_10 + b_11;
+    g_0 = 1.0f + b_00 + b_01;
+
+    nfc.base_gain /= g_1 * g_0;
+    nfc.gain /= g_1 * g_0;
+    nfc.a1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc.a2 = 4.0f * b_11 / g_1;
+    nfc.a3 = (2.0f*b_00 + 4.0f*b_01) / g_0;
+    nfc.a4 = 4.0f * b_01 / g_0;
+
+    return nfc;
+}
+
+void NfcFilterAdjust4(NfcFilter4 *nfc, const float w0) noexcept
+{
+    const float r{0.5f * w0};
+    const float b_10{B[4][0] * r};
+    const float b_11{B[4][1] * r * r};
+    const float b_00{B[4][2] * r};
+    const float b_01{B[4][3] * r * r};
+    const float g_1{1.0f + b_10 + b_11};
+    const float g_0{1.0f + b_00 + b_01};
+
+    nfc->gain = nfc->base_gain * g_1 * g_0;
+    nfc->b1 = (2.0f*b_10 + 4.0f*b_11) / g_1;
+    nfc->b2 = 4.0f * b_11 / g_1;
+    nfc->b3 = (2.0f*b_00 + 4.0f*b_01) / g_0;
+    nfc->b4 = 4.0f * b_01 / g_0;
+}
+
+} // namespace
+
+void NfcFilter::init(const float w1) noexcept
+{
+    first = NfcFilterCreate1(0.0f, w1);
+    second = NfcFilterCreate2(0.0f, w1);
+    third = NfcFilterCreate3(0.0f, w1);
+    fourth = NfcFilterCreate4(0.0f, w1);
+}
+
+void NfcFilter::adjust(const float w0) noexcept
+{
+    NfcFilterAdjust1(&first, w0);
+    NfcFilterAdjust2(&second, w0);
+    NfcFilterAdjust3(&third, w0);
+    NfcFilterAdjust4(&fourth, w0);
+}
+
+
+void NfcFilter::process1(const al::span<const float> src, float *RESTRICT dst)
+{
+    const float gain{first.gain};
+    const float b1{first.b1};
+    const float a1{first.a1};
+    float z1{first.z[0]};
+    auto proc_sample = [gain,b1,a1,&z1](const float in) noexcept -> float
+    {
+        const float y{in*gain - a1*z1};
+        const float out{y + b1*z1};
+        z1 += y;
+        return out;
+    };
+    std::transform(src.cbegin(), src.cend(), dst, proc_sample);
+    first.z[0] = z1;
+}
+
+void NfcFilter::process2(const al::span<const float> src, float *RESTRICT dst)
+{
+    const float gain{second.gain};
+    const float b1{second.b1};
+    const float b2{second.b2};
+    const float a1{second.a1};
+    const float a2{second.a2};
+    float z1{second.z[0]};
+    float z2{second.z[1]};
+    auto proc_sample = [gain,b1,b2,a1,a2,&z1,&z2](const float in) noexcept -> float
+    {
+        const float y{in*gain - a1*z1 - a2*z2};
+        const float out{y + b1*z1 + b2*z2};
+        z2 += z1;
+        z1 += y;
+        return out;
+    };
+    std::transform(src.cbegin(), src.cend(), dst, proc_sample);
+    second.z[0] = z1;
+    second.z[1] = z2;
+}
+
+void NfcFilter::process3(const al::span<const float> src, float *RESTRICT dst)
+{
+    const float gain{third.gain};
+    const float b1{third.b1};
+    const float b2{third.b2};
+    const float b3{third.b3};
+    const float a1{third.a1};
+    const float a2{third.a2};
+    const float a3{third.a3};
+    float z1{third.z[0]};
+    float z2{third.z[1]};
+    float z3{third.z[2]};
+    auto proc_sample = [gain,b1,b2,b3,a1,a2,a3,&z1,&z2,&z3](const float in) noexcept -> float
+    {
+        float y{in*gain - a1*z1 - a2*z2};
+        float out{y + b1*z1 + b2*z2};
+        z2 += z1;
+        z1 += y;
+
+        y = out - a3*z3;
+        out = y + b3*z3;
+        z3 += y;
+        return out;
+    };
+    std::transform(src.cbegin(), src.cend(), dst, proc_sample);
+    third.z[0] = z1;
+    third.z[1] = z2;
+    third.z[2] = z3;
+}
+
+void NfcFilter::process4(const al::span<const float> src, float *RESTRICT dst)
+{
+    const float gain{fourth.gain};
+    const float b1{fourth.b1};
+    const float b2{fourth.b2};
+    const float b3{fourth.b3};
+    const float b4{fourth.b4};
+    const float a1{fourth.a1};
+    const float a2{fourth.a2};
+    const float a3{fourth.a3};
+    const float a4{fourth.a4};
+    float z1{fourth.z[0]};
+    float z2{fourth.z[1]};
+    float z3{fourth.z[2]};
+    float z4{fourth.z[3]};
+    auto proc_sample = [gain,b1,b2,b3,b4,a1,a2,a3,a4,&z1,&z2,&z3,&z4](const float in) noexcept -> float
+    {
+        float y{in*gain - a1*z1 - a2*z2};
+        float out{y + b1*z1 + b2*z2};
+        z2 += z1;
+        z1 += y;
+
+        y = out - a3*z3 - a4*z4;
+        out = y + b3*z3 + b4*z4;
+        z4 += z3;
+        z3 += y;
+        return out;
+    };
+    std::transform(src.cbegin(), src.cend(), dst, proc_sample);
+    fourth.z[0] = z1;
+    fourth.z[1] = z2;
+    fourth.z[2] = z3;
+    fourth.z[3] = z4;
+}
diff --git a/core/filters/nfc.h b/core/filters/nfc.h
new file mode 100644
index 00000000..33f67a5f
--- /dev/null
+++ b/core/filters/nfc.h
@@ -0,0 +1,63 @@
+#ifndef CORE_FILTERS_NFC_H
+#define CORE_FILTERS_NFC_H
+
+#include <cstddef>
+
+#include "alspan.h"
+
+
+struct NfcFilter1 {
+    float base_gain, gain;
+    float b1, a1;
+    float z[1];
+};
+struct NfcFilter2 {
+    float base_gain, gain;
+    float b1, b2, a1, a2;
+    float z[2];
+};
+struct NfcFilter3 {
+    float base_gain, gain;
+    float b1, b2, b3, a1, a2, a3;
+    float z[3];
+};
+struct NfcFilter4 {
+    float base_gain, gain;
+    float b1, b2, b3, b4, a1, a2, a3, a4;
+    float z[4];
+};
+
+class NfcFilter {
+    NfcFilter1 first;
+    NfcFilter2 second;
+    NfcFilter3 third;
+    NfcFilter4 fourth;
+
+public:
+    /* NOTE:
+     * w0 = speed_of_sound / (source_distance * sample_rate);
+     * w1 = speed_of_sound / (control_distance * sample_rate);
+     *
+     * Generally speaking, the control distance should be approximately the
+     * average speaker distance, or based on the reference delay if outputing
+     * NFC-HOA. It must not be negative, 0, or infinite. The source distance
+     * should not be too small relative to the control distance.
+     */
+
+    void init(const float w1) noexcept;
+    void adjust(const float w0) noexcept;
+
+    /* Near-field control filter for first-order ambisonic channels (1-3). */
+    void process1(const al::span<const float> src, float *RESTRICT dst);
+
+    /* Near-field control filter for second-order ambisonic channels (4-8). */
+    void process2(const al::span<const float> src, float *RESTRICT dst);
+
+    /* Near-field control filter for third-order ambisonic channels (9-15). */
+    void process3(const al::span<const float> src, float *RESTRICT dst);
+
+    /* Near-field control filter for fourth-order ambisonic channels (16-24). */
+    void process4(const al::span<const float> src, float *RESTRICT dst);
+};
+
+#endif /* CORE_FILTERS_NFC_H */
diff --git a/core/filters/splitter.cpp b/core/filters/splitter.cpp
new file mode 100644
index 00000000..5cc670b7
--- /dev/null
+++ b/core/filters/splitter.cpp
@@ -0,0 +1,113 @@
+
+#include "config.h"
+
+#include "splitter.h"
+
+#include <algorithm>
+#include <cmath>
+#include <limits>
+
+#include "math_defs.h"
+#include "opthelpers.h"
+
+
+template<typename Real>
+void BandSplitterR<Real>::init(Real f0norm)
+{
+    const Real w{f0norm * al::MathDefs<Real>::Tau()};
+    const Real cw{std::cos(w)};
+    if(cw > std::numeric_limits<float>::epsilon())
+        mCoeff = (std::sin(w) - 1.0f) / cw;
+    else
+        mCoeff = cw * -0.5f;
+
+    mLpZ1 = 0.0f;
+    mLpZ2 = 0.0f;
+    mApZ1 = 0.0f;
+}
+
+template<typename Real>
+void BandSplitterR<Real>::process(const al::span<const Real> input, Real *hpout, Real *lpout)
+{
+    const Real ap_coeff{mCoeff};
+    const Real lp_coeff{mCoeff*0.5f + 0.5f};
+    Real lp_z1{mLpZ1};
+    Real lp_z2{mLpZ2};
+    Real ap_z1{mApZ1};
+    auto proc_sample = [ap_coeff,lp_coeff,&lp_z1,&lp_z2,&ap_z1,&lpout](const Real in) noexcept -> Real
+    {
+        /* Low-pass sample processing. */
+        Real d{(in - lp_z1) * lp_coeff};
+        Real lp_y{lp_z1 + d};
+        lp_z1 = lp_y + d;
+
+        d = (lp_y - lp_z2) * lp_coeff;
+        lp_y = lp_z2 + d;
+        lp_z2 = lp_y + d;
+
+        *(lpout++) = lp_y;
+
+        /* All-pass sample processing. */
+        Real ap_y{in*ap_coeff + ap_z1};
+        ap_z1 = in - ap_y*ap_coeff;
+
+        /* High-pass generated from removing low-passed output. */
+        return ap_y - lp_y;
+    };
+    std::transform(input.cbegin(), input.cend(), hpout, proc_sample);
+    mLpZ1 = lp_z1;
+    mLpZ2 = lp_z2;
+    mApZ1 = ap_z1;
+}
+
+template<typename Real>
+void BandSplitterR<Real>::processHfScale(const al::span<Real> samples, const Real hfscale)
+{
+    const Real ap_coeff{mCoeff};
+    const Real lp_coeff{mCoeff*0.5f + 0.5f};
+    Real lp_z1{mLpZ1};
+    Real lp_z2{mLpZ2};
+    Real ap_z1{mApZ1};
+    auto proc_sample = [hfscale,ap_coeff,lp_coeff,&lp_z1,&lp_z2,&ap_z1](const Real in) noexcept -> Real
+    {
+        /* Low-pass sample processing. */
+        Real d{(in - lp_z1) * lp_coeff};
+        Real lp_y{lp_z1 + d};
+        lp_z1 = lp_y + d;
+
+        d = (lp_y - lp_z2) * lp_coeff;
+        lp_y = lp_z2 + d;
+        lp_z2 = lp_y + d;
+
+        /* All-pass sample processing. */
+        Real ap_y{in*ap_coeff + ap_z1};
+        ap_z1 = in - ap_y*ap_coeff;
+
+        /* High-pass generated by removing the low-passed signal, which is then
+         * scaled and added back to the low-passed signal.
+         */
+        return (ap_y-lp_y)*hfscale + lp_y;
+    };
+    std::transform(samples.begin(), samples.end(), samples.begin(), proc_sample);
+    mLpZ1 = lp_z1;
+    mLpZ2 = lp_z2;
+    mApZ1 = ap_z1;
+}
+
+template<typename Real>
+void BandSplitterR<Real>::applyAllpass(const al::span<Real> samples) const
+{
+    const Real coeff{mCoeff};
+    Real z1{0.0f};
+    auto proc_sample = [coeff,&z1](const Real in) noexcept -> Real
+    {
+        const Real out{in*coeff + z1};
+        z1 = in - out*coeff;
+        return out;
+    };
+    std::transform(samples.begin(), samples.end(), samples.begin(), proc_sample);
+}
+
+
+template class BandSplitterR<float>;
+template class BandSplitterR<double>;
diff --git a/core/filters/splitter.h b/core/filters/splitter.h
new file mode 100644
index 00000000..ba548c10
--- /dev/null
+++ b/core/filters/splitter.h
@@ -0,0 +1,36 @@
+#ifndef CORE_FILTERS_SPLITTER_H
+#define CORE_FILTERS_SPLITTER_H
+
+#include <cstddef>
+
+#include "alspan.h"
+
+
+/* Band splitter. Splits a signal into two phase-matching frequency bands. */
+template<typename Real>
+class BandSplitterR {
+    Real mCoeff{0.0f};
+    Real mLpZ1{0.0f};
+    Real mLpZ2{0.0f};
+    Real mApZ1{0.0f};
+
+public:
+    BandSplitterR() = default;
+    BandSplitterR(const BandSplitterR&) = default;
+    BandSplitterR(Real f0norm) { init(f0norm); }
+
+    void init(Real f0norm);
+    void clear() noexcept { mLpZ1 = mLpZ2 = mApZ1 = 0.0f; }
+    void process(const al::span<const Real> input, Real *hpout, Real *lpout);
+
+    void processHfScale(const al::span<Real> samples, const Real hfscale);
+
+    /* The all-pass portion of the band splitter. Applies the same phase shift
+     * without splitting the signal. Note that each use of this method is
+     * indepedent, it does not track history between calls.
+     */
+    void applyAllpass(const al::span<Real> samples) const;
+};
+using BandSplitter = BandSplitterR<float>;
+
+#endif /* CORE_FILTERS_SPLITTER_H */