Polyphase interpolator

bench/FIRFilterBench.cpp

@@ -97,6 +97,28 @@ static auto makeChowPolyphaseDecimFIRChoices()
 
 auto chowPolyphaseDecimFIRChoices = makeChowPolyphaseDecimFIRChoices();
 
+static auto makeChowPolyphaseInterpFIR (int order)
+{
+    const auto coefsData = bench_utils::makeRandomVector<float> (order);
+
+    chowdsp::FIRPolyphaseInterpolator<float> filter {};
+    filter.prepare (2, 1, blockSize / 2, coefsData);
+
+    return filter;
+}
+
+static auto makeChowPolyphaseInterpFIRChoices()
+{
+    std::unordered_map<int, chowdsp::FIRPolyphaseInterpolator<float>> choices;
+
+    for (int order = startOrder; order <= endOrder; order *= orderMult)
+        choices.emplace (std::make_pair (order, makeChowPolyphaseInterpFIR (order)));
+
+    return choices;
+}
+
+auto chowPolyphaseInterpFIRChoices = makeChowPolyphaseInterpFIRChoices();
+
 static auto makeAudioBuffer()
 {
     auto bufferData = bench_utils::makeRandomVector<float> (blockSize);
@@ -141,14 +163,24 @@ static void ChowFIR (benchmark::State& state)
 }
 BENCHMARK (ChowFIR)->MinTime (1)->RangeMultiplier (orderMult)->Range (startOrder, endOrder);
 
-static void ChowPolyphaseFIR (benchmark::State& state)
+static void ChowPolyphaseDecimFIR (benchmark::State& state)
 {
     auto& fir = chowPolyphaseDecimFIRChoices[(int) state.range (0)];
     for (auto _ : state)
     {
         fir.processBlock (audioBuffer, audioBuffer);
     }
 }
-BENCHMARK (ChowPolyphaseFIR)->MinTime (1)->RangeMultiplier (orderMult)->Range (startOrder, endOrder);
+BENCHMARK (ChowPolyphaseDecimFIR)->MinTime (1)->RangeMultiplier (orderMult)->Range (startOrder, endOrder);
+
+static void ChowPolyphaseInterpFIR (benchmark::State& state)
+{
+    auto& fir = chowPolyphaseInterpFIRChoices[(int) state.range (0)];
+    for (auto _ : state)
+    {
+        fir.processBlock (chowdsp::BufferView { audioBuffer, 0, blockSize / 2 }, audioBuffer);
+    }
+}
+BENCHMARK (ChowPolyphaseInterpFIR)->MinTime (1)->RangeMultiplier (orderMult)->Range (startOrder, endOrder);
 
 BENCHMARK_MAIN();

modules/dsp/chowdsp_filters/Other/chowdsp_FIRPolyphaseDecimator.h

@@ -14,9 +14,9 @@ class FIRPolyphaseDecimator
     FIRPolyphaseDecimator () = default;
 
     /** Prepares the filter to process a stream of data, with the given configuration and filter coefficients. */
-    void prepare (int decimationFactor, int numChannels, int maxBlockSize, const nonstd::span<const T> coeffs)
+    void prepare (int decimationFactor, int numChannels, int maxBlockSizeIn, const nonstd::span<const T> coeffs)
     {
-        jassert (maxBlockSize % decimationFactor == 0);
+        jassert (maxBlockSizeIn % decimationFactor == 0);
         const auto numCoeffs = coeffs.size();
         const auto coeffsPerFilter = Math::ceiling_divide (numCoeffs, (size_t) decimationFactor);
 
@@ -29,7 +29,7 @@ class FIRPolyphaseDecimator
         filters.reserve ((size_t) decimationFactor);
         for (int i = 0; i < decimationFactor; ++i)
         {
-            buffers.emplace_back (numChannels, maxBlockSize / decimationFactor);
+            buffers.emplace_back (numChannels, maxBlockSizeIn / decimationFactor);
 
             auto& filter = filters.emplace_back (coeffsPerFilter);
             filter.prepare (numChannels);
@@ -44,7 +44,12 @@ class FIRPolyphaseDecimator
         }
     }
 
-    /** Processes a block of data. */
+    /**
+     * Processes a block of data.
+     *
+     * inBlock should have a size of numSamplesIn, and outBlock should have a size of
+     * numSamplesIn / decimationFactor.
+     */
     void processBlock (const T* inBlock, T* outBlock, const int numSamplesIn, const int channel = 0) noexcept
     {
         const auto decimationFactor = (int) filters.size();
@@ -79,7 +84,11 @@ class FIRPolyphaseDecimator
             juce::FloatVectorOperations::add (outBlock, bufferPtrs[filterIndex], numSamplesOut);
     }
 
-    /** Processes a block of data. */
+    /**
+     * Processes a block of data.
+     *
+     * bufferOut should have a size of bufferIn.getNumSamples() / decimationFactor.
+     */
     void processBlock (const BufferView<const T>& bufferIn, const BufferView<T>& bufferOut) noexcept
     {
         jassert (bufferIn.getNumChannels() == bufferOut.getNumChannels());

modules/dsp/chowdsp_filters/Other/chowdsp_FIRPolyphaseInterpolator.h

@@ -0,0 +1,90 @@
+#pragma once
+
+namespace chowdsp
+{
+/**
+ * A polyphase FIR interpolation filter.
+ *
+ * Reference: http://www.ws.binghamton.edu/fowler/fowler%20personal%20page/EE521_files/IV-05%20Polyphase%20FIlters%20Revised.pdf
+ */
+template <typename T>
+class FIRPolyphaseInterpolator
+{
+public:
+    FIRPolyphaseInterpolator() = default;
+
+    /** Prepares the filter to process a stream of data, with the given configuration and filter coefficients. */
+    void prepare (int interpolationFactor, int numChannels, int maxBlockSizeIn, const nonstd::span<const T> coeffs)
+    {
+        const auto numCoeffs = coeffs.size();
+        const auto coeffsPerFilter = Math::ceiling_divide (numCoeffs, (size_t) interpolationFactor);
+
+        std::vector<float> oneFilterCoeffs (coeffsPerFilter);
+
+        buffers.clear();
+        buffers.reserve ((size_t) interpolationFactor);
+        filters.clear();
+        filters.reserve ((size_t) interpolationFactor);
+        for (int i = 0; i < interpolationFactor; ++i)
+        {
+            buffers.emplace_back (numChannels, maxBlockSizeIn);
+
+            auto& filter = filters.emplace_back (coeffsPerFilter);
+            filter.prepare (numChannels);
+
+            std::fill (oneFilterCoeffs.begin(), oneFilterCoeffs.end(), T{});
+            for (size_t j = 0; j < coeffsPerFilter; ++j)
+            {
+                const auto index = (size_t) i + j * (size_t) interpolationFactor;
+                oneFilterCoeffs[j] = index >= coeffs.size() ? T{} : coeffs[index];
+            }
+            filter.setCoefficients (oneFilterCoeffs.data());
+        }
+    }
+
+    /**
+     * Processes a block of data.
+     *
+     * inBlock should have a size of numSamplesIn, and outBlock should have a size of
+     * numSamplesIn * interpolationFactor.
+     */
+    void processBlock (const T* inBlock, T* outBlock, const int numSamplesIn, const int channel = 0) noexcept
+    {
+        const auto interpolationFactor = (int) filters.size();
+
+        // set up sub-buffer pointers
+        auto* bufferPtrs = static_cast<T**> (alloca (sizeof (T*) * (size_t) interpolationFactor));
+        for (size_t filterIndex = 0; filterIndex < (size_t) interpolationFactor; ++filterIndex)
+            bufferPtrs[filterIndex] = buffers[filterIndex].getWritePointer (channel);
+
+        // process sub-buffers
+        for (size_t filterIndex = 0; filterIndex < (size_t) interpolationFactor; ++filterIndex)
+            filters[filterIndex].processBlock (inBlock, bufferPtrs[filterIndex], numSamplesIn, channel);
+
+        // fill output buffer
+        for (size_t filterIndex = 0; filterIndex < (size_t) interpolationFactor; ++filterIndex)
+            for (int n = 0; n < numSamplesIn; ++n)
+                outBlock[n * interpolationFactor + (int) filterIndex] = bufferPtrs[filterIndex][n];
+    }
+
+    /**
+     * Processes a block of data.
+     *
+     * bufferOut should have a size of bufferIn.getNumSamples() * interpolationFactor.
+     */
+    void processBlock (const BufferView<const T>& bufferIn, const BufferView<T>& bufferOut) noexcept
+    {
+        jassert (bufferIn.getNumChannels() == bufferOut.getNumChannels());
+        const auto numSamples = bufferIn.getNumSamples();
+        jassert (numSamples == bufferOut.getNumSamples() * (int) filters.size());
+
+        for (auto [ch, dataIn] : buffer_iters::channels (bufferIn))
+            processBlock (dataIn.data(), bufferOut.getWritePointer (ch), numSamples, ch);
+    }
+
+private:
+    std::vector<FIRFilter<T>> filters{};
+
+    std::vector<Buffer<T>> buffers{};
+};
+}

modules/dsp/chowdsp_filters/chowdsp_filters.h

@@ -47,6 +47,7 @@ BEGIN_JUCE_MODULE_DECLARATION
 #include "Other/chowdsp_HilbertFilter.h"
 #include "Other/chowdsp_FIRFilter.h"
 #include "Other/chowdsp_FIRPolyphaseDecimator.h"
+#include "Other/chowdsp_FIRPolyphaseInterpolator.h"
 #include "Other/chowdsp_WernerFilter.h"
 #include "Other/chowdsp_ARPFilter.h"
 #include "Other/chowdsp_LinkwitzRileyFilter.h"

tests/dsp_tests/chowdsp_filters_test/CMakeLists.txt

@@ -20,4 +20,5 @@ target_sources(chowdsp_filters_test
         LinearTransformsTest.cpp
         CrossoverFilterTest.cpp
         FIRPolyphaseDecimatorTest.cpp
+        FIRPolyphaseInterpolatorTest.cpp
 )

tests/dsp_tests/chowdsp_filters_test/FIRPolyphaseInterpolatorTest.cpp

@@ -0,0 +1,62 @@
+#include <CatchUtils.h>
+#include <chowdsp_filters/chowdsp_filters.h>
+
+static void interpolationFilterCompare (int filterOrder, int interpolationFactor, int numChannels)
+{
+    const auto numSamples = 12;
+    const auto halfSamples = numSamples / 2;
+
+    chowdsp::Buffer<float> bufferIn{ numChannels, numSamples };
+    for (auto [ch, data] : chowdsp::buffer_iters::channels (bufferIn))
+        for (auto [n, x] : chowdsp::enumerate (data))
+            x = static_cast<float> (n + (size_t) ch);
+
+    std::vector coeffs ((size_t) filterOrder, 0.0f);
+    for (auto [k, h] : chowdsp::enumerate (coeffs))
+        h = static_cast<float> (k);
+
+    chowdsp::FIRFilter<float> referenceFilter{ filterOrder };
+    referenceFilter.prepare (numChannels);
+    referenceFilter.setCoefficients (coeffs.data());
+
+    chowdsp::Buffer<float> referenceBufferIn{ numChannels, numSamples * interpolationFactor };
+    chowdsp::Buffer<float> referenceBufferOut{ numChannels, numSamples * interpolationFactor };
+    referenceBufferIn.clear();
+    for (auto [ch, data] : chowdsp::buffer_iters::channels (std::as_const (bufferIn)))
+        for (auto [n, x] : chowdsp::enumerate (data))
+            referenceBufferIn.getWritePointer (ch)[(int) n * interpolationFactor] = x;
+    referenceFilter.processBlock (chowdsp::BufferView{ referenceBufferIn, 0, halfSamples * interpolationFactor },
+                                  chowdsp::BufferView{ referenceBufferOut, 0, halfSamples * interpolationFactor });
+    referenceFilter.processBlock (chowdsp::BufferView{ referenceBufferIn, halfSamples * interpolationFactor, halfSamples * interpolationFactor },
+                                  chowdsp::BufferView{ referenceBufferOut, halfSamples * interpolationFactor, halfSamples * interpolationFactor });
+
+    chowdsp::FIRPolyphaseInterpolator<float> interpolatorFilter;
+    interpolatorFilter.prepare (interpolationFactor, numChannels, numSamples, coeffs);
+    chowdsp::Buffer<float> testBufferOut{ numChannels, numSamples * interpolationFactor };
+    interpolatorFilter.processBlock (chowdsp::BufferView{ bufferIn, 0, halfSamples },
+                                  chowdsp::BufferView{ testBufferOut, 0, halfSamples * interpolationFactor });
+    interpolatorFilter.processBlock (chowdsp::BufferView{ bufferIn, halfSamples, halfSamples },
+                                  chowdsp::BufferView{ testBufferOut, halfSamples * interpolationFactor, halfSamples * interpolationFactor });
+
+    for (int ch = 0; ch < numChannels; ++ch)
+    {
+        for (int n = 0; n < numSamples * interpolationFactor; ++n)
+        {
+            const auto ref = referenceBufferOut.getReadPointer (ch)[n];
+            const auto test = testBufferOut.getReadPointer (ch)[n];
+            REQUIRE (test == Catch::Approx { ref }.margin (1.0e-6));
+        }
+    }
+}
+
+TEST_CASE ("FIR Polyphase Interpolator Test", "[dsp][filters][fir][anti-aliasing]")
+{
+    interpolationFilterCompare (10, 2, 1);
+    interpolationFilterCompare (9, 2, 1);
+
+    interpolationFilterCompare (16, 3, 4);
+    interpolationFilterCompare (19, 3, 4);
+
+    interpolationFilterCompare (32, 4, 2);
+    interpolationFilterCompare (33, 4, 2);
+}