Channelizer/ChannelSynthesizer (#64)

* Add polyphase FFT filter banks
* Add `zadoff_chu` function
* Update tests

Author: vitalsong

CMakeLists.txt

@@ -45,6 +45,7 @@ set(DSPLIB_SOURCES
     lib/fft.cpp
     lib/ifft.cpp
     lib/czt.cpp
+    lib/subband.cpp
 )
 
 if (DSPLIB_EXCLUDE_FFT)

benchs/CMakeLists.txt

@@ -43,6 +43,7 @@ set(SOURCES
     resample.cpp
     slice.cpp
     fir.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::complex(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

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

include/dsplib/subband.h

@@ -0,0 +1,136 @@
+#pragma once
+
+#include <dsplib/array.h>
+
+#include <cassert>
+#include <memory>
+
+namespace dsplib {
+
+class ChannelizerImpl;
+
+//TODO: use decim factor in matlab style (M/D)
+
+/**
+ * @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]
+     */
+    explicit Channelizer(const arr_real& filter, int num_bands, int decim_factor);
+
+    /**
+     * @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]
+     */
+    explicit Channelizer(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor);
+
+    /**
+     * @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 (expected to be even)
+     */
+    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);
+    }
+
+    /**
+     * @brief Processing frame size
+     * @return int frame len
+     */
+    [[nodiscard]] int frame_len() const noexcept;
+
+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.
+ * @todo remove num_taps if filter is calculated
+ * @todo params order (bands, decim, ntaps/coeffs)
+ */
+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]
+     */
+    explicit ChannelSynthesizer(const arr_real& filter, int num_bands, int decim_factor);
+
+    /**
+     * @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]
+     */
+    explicit ChannelSynthesizer(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor);
+
+    /**
+     * @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 (expected to be even)
+     */
+    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);
+    }
+
+    /**
+     * @brief Processing frame size
+     * @return int frame len
+     */
+    [[nodiscard]] int frame_len() const noexcept;
+
+private:
+    std::shared_ptr<ChannelSynthesizerImpl> d_;
+};
+
+}   // namespace dsplib

include/dsplib/utils.h

@@ -240,4 +240,12 @@ arr_cmplx circshift(const arr_cmplx& x, int k);
 arr_real fftshift(const arr_real& x);
 arr_cmplx fftshift(const arr_cmplx& x);
 
+/**
+ * @brief generate Zadoff-Chu sequence
+ * @param r root, from 1 to n-1
+ * @param n length
+ * @return Zadoff-Chu sequence
+ */
+arr_cmplx zadoff_chu(int r, int n);
+
 }   // namespace dsplib

lib/corr.cpp

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

lib/czt.cpp

@@ -15,7 +15,7 @@ class CztPlanImpl
       : _n{n}
       , _m{m} {
         assert(abs(abs(w) - 1.0) < 2 * eps());
-        auto t = abs2(arange(1 - n, max(m, n))) / 2;
+        const auto t = abs2(arange(1 - n, max(m, n))) / 2;
         arr_cmplx chirp(t.size());
         const auto w_a = angle(w);
         for (int i = 0; i < t.size(); ++i) {