Fix API. Normalize gain for ChannelSynthesizer out.

Author: vitalsong

benchs/subband.cpp

@@ -18,7 +18,7 @@ static void BM_Channelizer(benchmark::State& state) {
 
 static void BM_ChannelSynthesizer(benchmark::State& state) {
     auto chan = dsplib::ChannelSynthesizer(N, D, M);
-    dsplib::arr_cmplx x = dsplib::randn(N);
+    dsplib::arr_cmplx x = dsplib::complex(dsplib::randn(N));
     for (auto _ : state) {
         auto y = chan.process(x);
         benchmark::DoNotOptimize(y);

include/dsplib/subband.h

@@ -9,6 +9,8 @@ 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
@@ -26,26 +28,24 @@ class 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);
+    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]
-     * @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);
+    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
+     * @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);
 
@@ -60,6 +60,12 @@ class Channelizer
         return this->process(x);
     }
 
+    /**
+     * @brief Processing frame size
+     * @return int frame len
+     */
+    [[nodiscard]] int frame_len() const noexcept;
+
 private:
     std::shared_ptr<ChannelizerImpl> d_;
 };
@@ -77,32 +83,33 @@ class ChannelSynthesizer
 public:
     ChannelSynthesizer(const ChannelSynthesizer&) = delete;
 
+    //TODO: remove num_taps if filter is calculated
+    //TODO: params order (bands, decim, ntaps/coeffs)
+
     /**
      * @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);
+    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]
-     * @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);
+    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
+     * @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);
 
@@ -117,6 +124,12 @@ class ChannelSynthesizer
         return this->process(x);
     }
 
+    /**
+     * @brief Processing frame size
+     * @return int frame len
+     */
+    [[nodiscard]] int frame_len() const noexcept;
+
 private:
     std::shared_ptr<ChannelSynthesizerImpl> d_;
 };

lib/subband.cpp

@@ -1,5 +1,5 @@
 #include "dsplib/subband.h"
-#include "dsplib/throw.h"
+#include "dsplib/assert.h"
 #include "dsplib/resample.h"
 #include "dsplib/fft.h"
 #include "dsplib/ifft.h"
@@ -124,7 +124,6 @@ class DFTFilterBank
 public:
     virtual ~DFTFilterBank() = default;
 
-protected:
     explicit DFTFilterBank(int num_bands, int decim_factor, int num_taps, bool synthesis = false)
       : nbands_{num_bands}
       , ntaps_{num_taps}
@@ -139,19 +138,21 @@ class DFTFilterBank
     const int decim_;
     const int d_;
 
+    //TODO: size optimization
     CircBuffer buf_;
     CircBuffer gsi_;
 };
 
 arr_real _design_filter(int num_bands, int num_taps) {
     //TODO: cache last coeffs
-    return design_multirate_fir(1, num_bands, std::ceil(num_taps / 2.0));
+    assert(num_taps % 2 == 0);
+    return design_multirate_fir(1, num_bands, num_taps / 2, 80);
 }
 
 }   // namespace
 
 //--------------------------------------------------------------------------------------------------------
-class ChannelizerImpl : protected DFTFilterBank
+class ChannelizerImpl : public DFTFilterBank
 {
 public:
     explicit ChannelizerImpl(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor, int num_taps)
@@ -173,8 +174,6 @@ class ChannelizerImpl : protected DFTFilterBank
                 convert[k + i * d_] = gsi[k];
             }
         }
-
-        //TODO: move or use `convert` array
         buf_.push(convert, true);
 
         // calculate outputs of polyphase filters
@@ -198,7 +197,7 @@ class ChannelizerImpl : protected DFTFilterBank
 };
 
 //--------------------------------------------------------------------------------------------------------------
-class ChannelSynthesizerImpl : private DFTFilterBank
+class ChannelSynthesizerImpl : public DFTFilterBank
 {
 public:
     explicit ChannelSynthesizerImpl(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor,
@@ -212,7 +211,7 @@ class ChannelSynthesizerImpl : private DFTFilterBank
     arr_real process(const dsplib::arr_cmplx& x) {
         DSPLIB_ASSERT(x.size() == nbands_, "input vector size error");
 
-        const auto xx = ifft_(x * nbands_);
+        const auto xx = ifft_(x) * nbands_;
         buf_.push(xx, true);
 
         // calculate outputs of polyphase filters
@@ -237,7 +236,10 @@ class ChannelSynthesizerImpl : private DFTFilterBank
                 out[d_ - k - 1] += gsi[k + i * d_];
             }
         }
-        return out;
+
+        //normalize ouput
+        //TODO: more precision, gain error ~ 1 dB
+        return out * (nbands_ / decim_);
     }
 
 private:
@@ -247,40 +249,53 @@ class ChannelSynthesizerImpl : private DFTFilterBank
 };
 
 //--------------------------------------------------------------------------------------------------------------
-Channelizer::Channelizer(const arr_real& filter, int num_bands, int decim_factor, int num_taps) {
-    const auto fptr = std::make_shared<dsplib::arr_real>(filter);
-    d_ = std::make_unique<ChannelizerImpl>(fptr, num_bands, decim_factor, num_taps);
+Channelizer::Channelizer(const arr_real& filter, int num_bands, int decim_factor)
+  : Channelizer(std::make_shared<dsplib::arr_real>(filter), num_bands, decim_factor) {
 }
 
-Channelizer::Channelizer(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor, int num_taps) {
+Channelizer::Channelizer(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor) {
+    const int num_taps = filter->size() / num_bands;
+    DSPLIB_ASSERT(num_taps % 2 == 0, "`num_taps` expected to be even");
+    DSPLIB_ASSERT(num_bands % decim_factor == 0, "only integer ratio M/D supported");
+    DSPLIB_ASSERT(filter->size() == num_bands * num_taps, "filter size must be equal `num_bands * num_taps`");
     d_ = std::make_unique<ChannelizerImpl>(std::move(filter), num_bands, decim_factor, num_taps);
 }
 
 Channelizer::Channelizer(int num_bands, int decim_factor, int num_taps)
-  : Channelizer(_design_filter(num_bands, num_taps), num_bands, decim_factor, num_taps) {
+  : Channelizer(_design_filter(num_bands, num_taps), num_bands, decim_factor) {
 }
 
 arr_cmplx Channelizer::process(const arr_real& x) {
     return d_->process(x);
 }
 
+int Channelizer::frame_len() const noexcept {
+    return (d_->nbands_ / d_->decim_);
+}
+
 //--------------------------------------------------------------------------------------------------------------
-ChannelSynthesizer::ChannelSynthesizer(const arr_real& filter, int num_bands, int decim_factor, int num_taps) {
-    const auto fptr = std::make_shared<dsplib::arr_real>(filter);
-    d_ = std::make_unique<ChannelSynthesizerImpl>(fptr, num_bands, decim_factor, num_taps);
+ChannelSynthesizer::ChannelSynthesizer(const arr_real& filter, int num_bands, int decim_factor)
+  : ChannelSynthesizer(std::make_shared<dsplib::arr_real>(filter), num_bands, decim_factor) {
 }
 
-ChannelSynthesizer::ChannelSynthesizer(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor,
-                                       int num_taps) {
+ChannelSynthesizer::ChannelSynthesizer(std::shared_ptr<const arr_real> filter, int num_bands, int decim_factor) {
+    const int num_taps = filter->size() / num_bands;
+    DSPLIB_ASSERT(num_taps % 2 == 0, "`num_taps` expected to be even");
+    DSPLIB_ASSERT(num_bands % decim_factor == 0, "only integer ratio M/D supported");
+    DSPLIB_ASSERT(filter->size() == num_bands * num_taps, "filter size must be equal `num_bands * num_taps`");
     d_ = std::make_unique<ChannelSynthesizerImpl>(std::move(filter), num_bands, decim_factor, num_taps);
 }
 
 ChannelSynthesizer::ChannelSynthesizer(int num_bands, int decim_factor, int num_taps)
-  : ChannelSynthesizer(_design_filter(num_bands, num_taps), num_bands, decim_factor, num_taps) {
+  : ChannelSynthesizer(_design_filter(num_bands, num_taps), num_bands, decim_factor) {
 }
 
 arr_real ChannelSynthesizer::process(const arr_cmplx& x) {
     return d_->process(x);
 }
 
+int ChannelSynthesizer::frame_len() const noexcept {
+    return (d_->nbands_ / d_->decim_);
+}
+
 }   // namespace dsplib