Polyphase FFT filter bank

benchs/CMakeLists.txt

@@ -42,6 +42,7 @@ set(SOURCES
     adaptive.cpp
     resample.cpp
     slice.cpp
+    subband.cpp
 )
 
 add_executable(${PROJECT_NAME} ${SOURCES})

benchs/subband.cpp

@@ -0,0 +1,29 @@
+#include <benchmark/benchmark.h>
+#include <dsplib.h>
+
+#include "dsplib/subband.h"
+
+constexpr int M = 4;
+constexpr int N = 1024;
+constexpr int D = 2;
+
+static void BM_Channelizer(benchmark::State& state) {
+    auto chan = dsplib::Channelizer(N, D, M);
+    auto x = dsplib::randn(N / D);
+    for (auto _ : state) {
+        auto y = chan.process(x);
+        benchmark::DoNotOptimize(y);
+    }
+}
+
+static void BM_ChannelSynthesizer(benchmark::State& state) {
+    auto chan = dsplib::ChannelSynthesizer(N, D, M);
+    dsplib::arr_cmplx x = dsplib::randn(N);
+    for (auto _ : state) {
+        auto y = chan.process(x);
+        benchmark::DoNotOptimize(y);
+    }
+}
+
+BENCHMARK(BM_Channelizer)->Unit(benchmark::kMicrosecond)->MinTime(5);
+BENCHMARK(BM_ChannelSynthesizer)->Unit(benchmark::kMicrosecond)->MinTime(5);

include/dsplib.h

@@ -24,6 +24,7 @@
 #include <dsplib/resample.h>
 #include <dsplib/spectrum.h>
 #include <dsplib/stft.h>
+#include <dsplib/subband.h>
 
 #include <dsplib/audio/noise-gate.h>
 #include <dsplib/audio/compressor.h>

include/dsplib/subband.h

@@ -0,0 +1,124 @@
+#pragma once
+
+#include <dsplib/array.h>
+
+#include <cassert>
+#include <memory>
+
+namespace dsplib {
+
+class ChannelizerImpl;
+
+/**
+ * @brief Polyphase FFT analysis filter bank
+ * @details Separates a broadband input signal into multiple narrow subbands
+ * @see Matlab dsp.Channelizer
+ * @todo Add `Range` parameter [OneSided, Twosided, Centered]
+ */
+class Channelizer
+{
+public:
+    Channelizer(const Channelizer&) = delete;
+
+    /**
+     * @brief Construct Channelizer
+     * 
+     * @param filter Multirate FIR coeffs [num_bands * ntaps]
+     * @param num_bands Number of frequency bands
+     * @param decim_factor Decimation factor [1 : M-1]
+     * @param num_taps Number of filter coefficients per frequency band
+     */
+    explicit Channelizer(const arr_real& filter, int num_bands, int decim_factor, int num_taps);
+
+    /**
+     * @brief Construct Channelizer
+     * @details Use this function to save memory if multiple Channelizer/ChannelSynthesizer objects can exist at the same time
+     * @param filter Pointer to multirate FIR coeffs
+     * @param num_bands Number of frequency bands
+     * @param decim_factor Decimation factor [1 : M-1]
+     * @param num_taps Number of filter coefficients per frequency band
+     */
+    explicit Channelizer(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor, int num_taps);
+
+    /**
+     * @brief Construct Channelizer
+     * @details The filter will be calculated using the `design_multirate_fir(1, num_bands, ceil(num_taps / 2.0))`
+     * @param num_bands Number of frequency bands
+     * @param decim_factor Decimation factor [1 : M-1]
+     * @param num_taps Number of filter coefficients per frequency band
+     */
+    explicit Channelizer(int num_bands, int decim_factor, int num_taps);
+
+    /**
+     * @brief Filter bank processing (analysis)
+     * @param x Input broadband signal [num_bands / decim_factor]
+     * @return arr_cmplx Subband signal [num_bands]
+     */
+    [[nodiscard]] arr_cmplx process(const arr_real& x);
+
+    arr_cmplx operator()(const arr_real& x) {
+        return this->process(x);
+    }
+
+private:
+    std::shared_ptr<ChannelizerImpl> d_;
+};
+
+class ChannelSynthesizerImpl;
+
+/**
+ * @brief Polyphase FFT synthesis filter bank
+ * @details Ideally ChannelSynthesizer(Channelizer(x)) == x.
+ * @see Matlab dsp.ChannelSynthesizer
+ * @warning This implementation differs from matlab in the decimation parameter. For decim_factor=1, the results should be identical.
+ */
+class ChannelSynthesizer
+{
+public:
+    ChannelSynthesizer(const ChannelSynthesizer&) = delete;
+
+    /**
+     * @brief Construct ChannelSynthesizer
+     * 
+     * @param filter Multirate FIR coeffs [num_bands * ntaps]
+     * @param num_bands Number of frequency bands
+     * @param decim_factor Decimation factor [1 : M-1]
+     * @param num_taps Number of filter coefficients per frequency band
+     */
+    explicit ChannelSynthesizer(const arr_real& filter, int num_bands, int decim_factor, int num_taps);
+
+    /**
+     * @brief Construct ChannelSynthesizer
+     * @details Use this function to save memory if multiple Channelizer/ChannelSynthesizer objects can exist at the same time
+     * @param filter Pointer to multirate FIR coeffs
+     * @param num_bands Number of frequency bands
+     * @param decim_factor Decimation factor [1 : M-1]
+     * @param num_taps Number of filter coefficients per frequency band
+     */
+    explicit ChannelSynthesizer(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor, int num_taps);
+
+    /**
+     * @brief Construct ChannelSynthesizer
+     * @details The filter will be calculated using the `design_multirate_fir(1, num_bands, ceil(num_taps / 2.0))`
+     * @param num_bands Number of frequency bands
+     * @param decim_factor Decimation factor [1 : M-1]
+     * @param num_taps Number of filter coefficients per frequency band
+     */
+    explicit ChannelSynthesizer(int num_bands, int decim_factor, int num_taps);
+
+    /**
+     * @brief Filter bank processing (synthesis)
+     * @param x Input subband signal [num_bands]
+     * @return arr_cmplx Restored broadband signal [num_bands / decim_factor]
+     */
+    [[nodiscard]] arr_real process(const arr_cmplx& x);
+
+    arr_real operator()(const arr_cmplx& x) {
+        return this->process(x);
+    }
+
+private:
+    std::shared_ptr<ChannelSynthesizerImpl> d_;
+};
+
+}   // namespace dsplib

lib/corr.cpp

@@ -82,6 +82,7 @@ real_t corr(const arr_real& x, const arr_real& y, Correlation type) {
     case Correlation::Kendall:
         return _kendall_corr(x, y);
     }
+    return 0;
 }
 
 }   // namespace dsplib