Adding polyphase FIR filters

bench/FIRFilterBench.cpp

@@ -75,6 +75,50 @@ static auto makeChowFIRChoices()
 
 auto chowFIRChoices = makeChowFIRChoices();
 
+static auto makeChowPolyphaseDecimFIR (int order)
+{
+    const auto coefsData = bench_utils::makeRandomVector<float> (order);
+
+    chowdsp::FIRPolyphaseDecimator<float> filter {};
+    filter.prepare (2, 1, blockSize, coefsData);
+
+    return filter;
+}
+
+static auto makeChowPolyphaseDecimFIRChoices()
+{
+    std::unordered_map<int, chowdsp::FIRPolyphaseDecimator<float>> choices;
+
+    for (int order = startOrder; order <= endOrder; order *= orderMult)
+        choices.emplace (std::make_pair (order, makeChowPolyphaseDecimFIR (order)));
+
+    return choices;
+}
+
+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);
@@ -114,10 +158,29 @@ static void ChowFIR (benchmark::State& state)
     auto& fir = chowFIRChoices[(int) state.range (0)];
     for (auto _ : state)
     {
-        auto&& audioBlock = juce::dsp::AudioBlock<float> { audioBuffer };
-        fir.processBlock (audioBlock);
+        fir.processBlock (audioBuffer);
     }
 }
 BENCHMARK (ChowFIR)->MinTime (1)->RangeMultiplier (orderMult)->Range (startOrder, endOrder);
 
+static void ChowPolyphaseDecimFIR (benchmark::State& state)
+{
+    auto& fir = chowPolyphaseDecimFIRChoices[(int) state.range (0)];
+    for (auto _ : state)
+    {
+        fir.processBlock (audioBuffer, audioBuffer);
+    }
+}
+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_FIRFilter.h

@@ -56,21 +56,38 @@
 
     /** Process block of samples */
     void processBlock (FloatType* block, const int numSamples, const int channel = 0) noexcept
+    {
+        processBlock (block, block, numSamples, channel);
+    }
+
+    /** Process block of samples out-of-place */
+    void processBlock (const FloatType* blockIn, FloatType* blockOut, const int numSamples, const int channel = 0) noexcept
     {
         auto* z = state[channel].data();
         const auto* h = coefficients.data();
         ScopedValue zPtrLocal { zPtr[channel] };
 
         for (int n = 0; n < numSamples; ++n)
-            block[n] = processSampleInternal (block[n], z, h, zPtrLocal.get(), order, paddedOrder);
+            blockOut[n] = processSampleInternal (blockIn[n], z, h, zPtrLocal.get(), order, paddedOrder);
     }
 
     /** Process block of samples */
     void processBlock (const BufferView<FloatType>& block) noexcept
     {
-        const auto numSamples = block.getNumSamples();
-        for (const auto [channel, data] : buffer_iters::channels (block))
-            processBlock (data.data(), numSamples, channel);
+        processBlock (block, block);
+    }
+
+    /** Process block of samples out-of-place */
+    void processBlock (const BufferView<const FloatType>& blockIn, const BufferView<FloatType>& blockOut) noexcept
+    {
+        jassert (blockIn.getNumChannels() == blockOut.getNumChannels());
+        jassert (blockIn.getNumSamples() == blockOut.getNumSamples());
+
+        const auto numChannels = blockIn.getNumChannels();
+        const auto numSamples = blockIn.getNumSamples();
+
+        for (int ch = 0; ch < numChannels; ++ch)
+            processBlock (blockIn.getReadPointer (ch), blockOut.getWritePointer (ch), numSamples, ch);
     }
 
     /**

modules/dsp/chowdsp_filters/Other/chowdsp_FIRPolyphaseDecimator.h

@@ -0,0 +1,118 @@
+#pragma once
+
+namespace chowdsp
+{
+/**
+ * A polyphase FIR decimation filter.
+ *
+ * Reference: http://www.ws.binghamton.edu/fowler/fowler%20personal%20page/EE521_files/IV-05%20Polyphase%20FIlters%20Revised.pdf
+ */
+template <typename T>
+class FIRPolyphaseDecimator
+{
+public:
+    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 maxBlockSizeIn, const nonstd::span<const T> coeffs)
+    {
+        jassert (maxBlockSizeIn % decimationFactor == 0);
+        const auto numCoeffs = coeffs.size();
+        const auto coeffsPerFilter = Math::ceiling_divide (numCoeffs, (size_t) decimationFactor);
+
+        std::vector<float> oneFilterCoeffs (coeffsPerFilter);
+
+#if JUCE_TEENSY
+        bufferPtrs.resize ((size_t) decimationFactor);
+#endif
+
+        buffers.clear();
+        buffers.reserve ((size_t) decimationFactor);
+        overlapState = std::vector<T> ((size_t) numChannels * (size_t) decimationFactor, T {});
+        filters.clear();
+        filters.reserve ((size_t) decimationFactor);
+        for (int i = 0; i < decimationFactor; ++i)
+        {
+            buffers.emplace_back (numChannels, maxBlockSizeIn / decimationFactor);
+
+            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) decimationFactor;
+                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 / decimationFactor.
+     */
+    void processBlock (const T* inBlock, T* outBlock, const int numSamplesIn, const int channel = 0) noexcept
+    {
+        const auto decimationFactor = (int) filters.size();
+        jassert (numSamplesIn % decimationFactor == 0);
+        const auto numSamplesOut = numSamplesIn / decimationFactor;
+
+        // set up sub-buffer pointers
+#if ! JUCE_TEENSY
+        auto* bufferPtrs = static_cast<T**> (alloca (sizeof (T*) * (size_t) decimationFactor));
+#endif
+        for (size_t filterIndex = 0; filterIndex < (size_t) decimationFactor; ++filterIndex)
+            bufferPtrs[filterIndex] = buffers[filterIndex].getWritePointer (channel);
+        auto* channelOverlapState = overlapState.data() + (size_t) channel * (size_t) decimationFactor;
+
+        // fill sub-buffers
+        for (int n = 0; n < numSamplesOut; ++n)
+            bufferPtrs[0][n] = inBlock[n * decimationFactor];
+
+        for (size_t filterIndex = 1; filterIndex < (size_t) decimationFactor; ++filterIndex)
+        {
+            bufferPtrs[(size_t) decimationFactor - filterIndex][0] = channelOverlapState[filterIndex];
+            for (int n = 0; n < numSamplesOut - 1; ++n)
+                bufferPtrs[(size_t) decimationFactor - filterIndex][n + 1] = inBlock[n * decimationFactor + (int) filterIndex];
+            channelOverlapState[filterIndex] = inBlock[(numSamplesOut - 1) * decimationFactor + (int) filterIndex];
+        }
+
+        // process sub-buffers
+        for (size_t filterIndex = 0; filterIndex < (size_t) decimationFactor; ++filterIndex)
+            filters[filterIndex].processBlock (bufferPtrs[filterIndex], numSamplesOut, channel);
+
+        // sum the sub-buffers
+        juce::FloatVectorOperations::copy (outBlock, bufferPtrs[0], numSamplesOut);
+        for (size_t filterIndex = 1; filterIndex < (size_t) decimationFactor; ++filterIndex)
+            juce::FloatVectorOperations::add (outBlock, bufferPtrs[filterIndex], numSamplesOut);
+    }
+
+    /**
+     * 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());
+        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 {};
+    std::vector<T> overlapState {};
+
+#if JUCE_TEENSY
+    std::vector<T*> bufferPtrs {};
+#endif
+};
+} // namespace chowdsp

modules/dsp/chowdsp_filters/Other/chowdsp_FIRPolyphaseInterpolator.h

@@ -0,0 +1,100 @@
+#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);
+
+#if JUCE_TEENSY
+        bufferPtrs.resize ((size_t) interpolationFactor);
+#endif
+
+        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
+#if ! JUCE_TEENSY
+        auto* bufferPtrs = static_cast<T**> (alloca (sizeof (T*) * (size_t) interpolationFactor));
+#endif
+        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 {};
+
+#if JUCE_TEENSY
+    std::vector<T*> bufferPtrs {};
+#endif
+};
+} // namespace chowdsp

modules/dsp/chowdsp_filters/Utils/chowdsp_CoefficientCalculators.h

@@ -480,7 +480,7 @@ namespace CoefficientCalculators
             a_copy[2] = ipow<2> (expBw);
 
             const auto f0 = juce::MathConstants<NumericType>::twoPi / w0;
-            const auto w_N_Sq = ipow<2> ((T) 0.49 * f0);
+            const auto w_N_Sq = ipow<2> ((T) (NumericType) 0.49 * f0);
             const auto C = ipow<2> (Beta * f0); // = (2*Beta)^2 w_N^2
             const auto Gn2 = g * (ipow<2> (1 - Alpha * w_N_Sq) + C) / (ipow<2> (Alpha - w_N_Sq) + C);
 

modules/dsp/chowdsp_filters/chowdsp_filters.h

@@ -46,6 +46,8 @@ BEGIN_JUCE_MODULE_DECLARATION
 #include "Other/chowdsp_FractionalOrderFilter.h"
 #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_dsp_juce_test/FIRFilterTest.cpp

@@ -3,7 +3,7 @@
 #include <juce_dsp/juce_dsp.h>
 #include <chowdsp_filters/chowdsp_filters.h>
 
-TEST_CASE ("FIR Filter Test", "[dsp][filters]")
+TEST_CASE ("FIR Filter Test", "[dsp][filters][fir]")
 {
     SECTION ("FIR Filter Test")
     {

tests/dsp_tests/chowdsp_filters_test/ButterworthFilterTest.cpp

@@ -27,7 +27,7 @@ void testFilter (Filter& filt, std::vector<float> freqs, std::vector<float> mags
         testFrequency (freqs[i], mags[i], errs[i], "Incorrect gain at " + messages[i] + " frequency.");
 }
 
-TEMPLATE_TEST_CASE ("Butterworth Filter Test", "[dsp][filters]", float)
+TEMPLATE_TEST_CASE ("Butterworth Filter Test", "[dsp][filters][anti-aliasing]", float)
 {
     using T = TestType;
     using FilterType = chowdsp::ButterworthFilterType;

tests/dsp_tests/chowdsp_filters_test/CMakeLists.txt

@@ -19,4 +19,6 @@ target_sources(chowdsp_filters_test
         ButterQsTest.cpp
         LinearTransformsTest.cpp
         CrossoverFilterTest.cpp
+        FIRPolyphaseDecimatorTest.cpp
+        FIRPolyphaseInterpolatorTest.cpp
 )

tests/dsp_tests/chowdsp_filters_test/ChebyshevIIFilterTest.cpp

@@ -30,7 +30,7 @@ void testFilter (Filter& filt, std::vector<float> freqs, std::vector<float> mags
         testFrequency (freqs[i], mags[i], errs[i], "Incorrect gain at " + messages[i] + " frequency.");
 }
 
-TEMPLATE_TEST_CASE ("Chebyshev II Filter Test", "[dsp][filters]", float)
+TEMPLATE_TEST_CASE ("Chebyshev II Filter Test", "[dsp][filters][anti-aliasing]", float)
 {
     using T = TestType;
     using FilterType = chowdsp::ChebyshevFilterType;