Add variant of BBD for multi-channel

Author: jpcima

sources/bbd/bbd_line.cpp

@@ -4,119 +4,143 @@
 
 static constexpr unsigned interp_size = 128;
 
-void BBD_Line::setup(double fs, unsigned ns, const BBD_Filter_Spec &fsin, const BBD_Filter_Spec &fsout)
+template <unsigned Channels>
+void BBD_Line<Channels>::setup(double fs, unsigned ns, const BBD_Filter_Spec &fsin, const BBD_Filter_Spec &fsout)
 {
-    mem_.reserve(8192);
+    mem_.reserve(8192 * Channels);
 
     const BBD_Filter_Coef &fin = BBD::compute_filter_cached(fs, interp_size, fsin);
     const BBD_Filter_Coef &fout = BBD::compute_filter_cached(fs, interp_size, fsout);
     fin_ = &fin;
     fout_ = &fout;
 
-    unsigned Min = fin.M;
-    unsigned Mout = fout.M;
-    Xin_.reset(new cdouble[Min]);
-    Xout_.reset(new cdouble[Mout]);
-    Xout_mem_.reset(new cdouble[Mout]);
-    Gin_.reset(new cdouble[Min]);
-    Gout_.reset(new cdouble[Mout]);
+    assert(fin.M <= Mmax);
+    assert(fout.M <= Mmax);
 
     set_delay_size(ns);
     clear();
 }
 
-void BBD_Line::set_delay_size(unsigned ns)
+template <unsigned Channels>
+void BBD_Line<Channels>::set_delay_size(unsigned ns)
 {
     mem_.clear();
-    mem_.resize(ns);
+    mem_.resize(ns * Channels);
     imem_ = 0;
     ns_ = ns;
 }
 
-void BBD_Line::clear()
+template <unsigned Channels>
+void BBD_Line<Channels>::clear()
 {
     std::fill(mem_.begin(), mem_.end(), 0);
     imem_ = 0;
     pclk_ = 0;
     ptick_ = 0;
-    ybbd_old_ = 0;
-    unsigned Min = fin_->M;
-    unsigned Mout = fout_->M;
-    std::fill(&Xin_[0], &Xin_[Min], 0);
-    std::fill(&Xout_[0], &Xout_[Mout], 0);
-    std::fill(&Xout_mem_[0], &Xout_mem_[Mout], 0);
+    ybbd1_.fill(0);
+    for (unsigned c = 0; c < Channels; ++c) {
+        Xin_[c].fill(0);
+        Xout1_[c].fill(0);
+    }
 }
 
-void BBD_Line::process(unsigned n, const float *input, float *output, const float *clock)
+template <unsigned Channels>
+void BBD_Line<Channels>::process(unsigned n, const float *const inputs[Channels], float *const outputs[Channels], const float *clock)
 {
     unsigned ns = ns_;
     float *mem = mem_.data();
     unsigned imem = imem_;
     double pclk = pclk_;
     unsigned ptick = ptick_;
-    double ybbd_old = ybbd_old_;
+    std::array<double, Channels> ybbd1 = ybbd1_;
 
-    const BBD_Filter_Coef &fin = *fin_, &fout = *fout_;
-    unsigned Min = fin.M, Mout = fout.M;
-    cdouble *Xin = Xin_.get(), *Xout = Xout_.get();
-    cdouble *Xout_mem = Xout_mem_.get();
-    cdouble *Gin = Gin_.get(), *Gout = Gout_.get();
-    const cdouble *Pin = fin.P.get(), *Pout = fout.P.get();
+    const BBD_Filter_Coef &fin = *fin_;
+    const BBD_Filter_Coef &fout = *fout_;
+    unsigned Min = fin.M;
+    unsigned Mout = fout.M;
+    const cdouble *const Pin = fin.P.get();
+    const cdouble *const Pout = fout.P.get();
 
     for (unsigned i = 0; i < n; ++i) {
-        double fclk = clock[i];
-
-        for (unsigned m = 0; m < Mout; ++m)
-            Xout[m] = 0;
+        std::array<cdouble, Mmax> Xout_[Channels] = {};
 
+        double fclk = clock[i];
         if (fclk > 0) {
             double pclk_old = pclk;
             pclk += fclk;
+
             unsigned tick_count = (unsigned)pclk;
             pclk -= tick_count;
+
             for (unsigned tick = 0; tick < tick_count; ++tick) {
                 double d = (1 - pclk_old + tick) * (1 / fclk);
                 d -= (unsigned)d;
+
                 if ((ptick & 1) == 0) {
+                    cdouble Gin[Mmax];
                     fin.interpolate_G(d, Gin);
-                    cdouble s = 0;
-                    for (unsigned m = 0; m < Min; ++m)
-                        s += Gin[m] * Xin[m];
-                    mem[imem] = s.real();
+
+                    for (unsigned c = 0; c < Channels; ++c) {
+                        const std::array<cdouble, Mmax> &Xin = Xin_[c];
+                        cdouble s = 0;
+                        for (unsigned m = 0; m < Min; ++m)
+                            s += Gin[m] * Xin[m];
+                        mem[imem * Channels + c] = s.real();
+                    }
+
                     imem = ((imem + 1) < ns) ? (imem + 1) : 0;
                 }
                 else {
+                    cdouble Gout[Mmax];
                     fout.interpolate_G(d, Gout);
-                    double ybbd = mem[imem];
-                    double delta = ybbd - ybbd_old;
-                    ybbd_old = ybbd;
-                    for (unsigned m = 0; m < Mout; ++m)
-                        Xout[m] += Gout[m] * delta;
+
+                    for (unsigned c = 0; c < Channels; ++c) {
+                        std::array<cdouble, Mmax> &Xout = Xout_[c];
+                        double ybbd = mem[imem * Channels + c];
+                        double delta = ybbd - ybbd1[c];
+                        ybbd1[c] = ybbd;
+                        for (unsigned m = 0; m < Mout; ++m)
+                            Xout[m] += Gout[m] * delta;
+                    }
                 }
+
                 ++ptick;
             }
         }
 
-        for (unsigned m = 0; m < Min; ++m)
-            Xin[m] = Pin[m] * Xin[m] + cdouble(input[i]);
+        for (unsigned c = 0; c < Channels; ++c) {
+            const float *input = inputs[c];
+            float *output = outputs[c];
 
-        cdouble y = fout.H * ybbd_old;
-        for (unsigned m = 0; m < Mout; ++m) {
-            cdouble xout = Pout[m] * Xout_mem[m] + Xout[m];
-            Xout_mem[m] = xout;
-            y += xout;
-        }
+            std::array<cdouble, Mmax> &Xin = Xin_[c];
+            std::array<cdouble, Mmax> &Xout1 = Xout1_[c];
+            const std::array<cdouble, Mmax> &Xout = Xout_[c];
 
-        output[i] = y.real();
+            for (unsigned m = 0; m < Min; ++m)
+                Xin[m] = Pin[m] * Xin[m] + cdouble(input[i]);
+
+            cdouble y = fout.H * ybbd1[c];
+            for (unsigned m = 0; m < Mout; ++m) {
+                cdouble xout = Pout[m] * Xout1[m] + Xout[m];
+                Xout1[m] = xout;
+                y += xout;
+            }
+
+            output[i] = y.real();
+        }
     }
 
     imem_ = imem;
     pclk_ = pclk;
     ptick_ = ptick;
-    ybbd_old_ = ybbd_old;
+    ybbd1_ = ybbd1;
 }
 
-void BBD_Line::process(unsigned n, float *inout, const float *clock)
+template <unsigned Channels>
+void BBD_Line<Channels>::process(unsigned n, float *const inouts[Channels], const float *clock)
 {
-    process(n, inout, inout, clock);
+    process(n, inouts, inouts, clock);
 }
+
+template class BBD_Line<1>;
+template class BBD_Line<2>;

sources/bbd/bbd_line.h

@@ -2,17 +2,18 @@
 #include "bbd_filter.h"
 #include <algorithm>
 #include <vector>
-#include <memory>
+#include <array>
 #include <complex>
 typedef std::complex<double> cdouble;
 
+template <unsigned Channels = 1>
 class BBD_Line {
 public:
     void setup(double fs, unsigned ns, const BBD_Filter_Spec &fsin, const BBD_Filter_Spec &fsout);
     void set_delay_size(unsigned ns);
     void clear();
-    void process(unsigned n, const float *input, float *output, const float *clock);
-    void process(unsigned n, float *inout, const float *clock);
+    void process(unsigned n, const float *const inputs[Channels], float *const outputs[Channels], const float *clock);
+    void process(unsigned n, float *const inouts[Channels], const float *clock);
 
     const BBD_Filter_Coef &filter_in() const noexcept { return *fin_; }
     const BBD_Filter_Coef &filter_out() const noexcept { return *fout_; }
@@ -28,12 +29,11 @@ class BBD_Line {
     unsigned imem_; // delay memory index
     double pclk_; // clock phase
     unsigned ptick_; // clock tick counter
-    double ybbd_old_;
     const BBD_Filter_Coef *fin_;
     const BBD_Filter_Coef *fout_;
-    std::unique_ptr<cdouble[]> Xin_;
-    std::unique_ptr<cdouble[]> Xout_;
-    std::unique_ptr<cdouble[]> Xout_mem_; // sample memory of output filter
-    std::unique_ptr<cdouble[]> Gin_;
-    std::unique_ptr<cdouble[]> Gout_;
+    std::array<double, Channels> ybbd1_;
+
+    enum { Mmax = 8 }; // the maximum filter order
+    std::array<cdouble, Mmax> Xin_[Channels];
+    std::array<cdouble, Mmax> Xout1_[Channels]; // sample memory of output filter
 };

sources/dsp/Delay3Phase.cpp

@@ -66,20 +66,20 @@ void Delay3Phase::AnalogDelay::process(const float *input, const float *const mo
     ///
     float avgDelay = 5e-3f;
     float varDelay = 1e-3f;
-    float delayClockLo = BBD_Line::hz_rate_for_delay(avgDelay - varDelay, NumAnalogBBDStages) * sampleTime;
-    float delayClockHi = BBD_Line::hz_rate_for_delay(avgDelay + varDelay, NumAnalogBBDStages) * sampleTime;
+    float delayClockLo = BBD_Line<1>::hz_rate_for_delay(avgDelay - varDelay, NumAnalogBBDStages) * sampleTime;
+    float delayClockHi = BBD_Line<1>::hz_rate_for_delay(avgDelay + varDelay, NumAnalogBBDStages) * sampleTime;
 
     ///
     for (unsigned l = 0; l < 3; ++l) {
-        BBD_Line &line = fDelayLine[l];
+        BBD_Line<1> &line = fDelayLine[l];
         float *lineOutput = lineOutputs[l];
         const float *mod = mods[l];
 
         float clock[BufferLimit];
         for (unsigned i = 0; i < count; ++i)
             clock[i] = delayClockLo + (delayClockHi - delayClockLo) * 0.5f * (1.0f + mod[i]);
 
-        line.process(count, input, lineOutput, clock);
+        line.process(count, &input, &lineOutput, clock);
     }
 
     float *outputL = outputs[0];

sources/dsp/Delay3Phase.h

@@ -23,7 +23,7 @@ class Delay3Phase {
 
     private:
         float fSampleTime = 0;
-        BBD_Line fDelayLine[3];
+        BBD_Line<1> fDelayLine[3];
     };
 
 private: