LCOV - code coverage report
Current view: top level - Preamplifier - Triode2Filter.cpp (source / functions) Hit Total Coverage
Test: coverage.info.cleaned Lines: 92 95 96.8 %
Date: 2021-02-18 20:07:22 Functions: 13 90 14.4 %

          Line data    Source code
       1             : /**
       2             :  * \file Triode2Filter.cpp
       3             :  */
       4             : 
       5             : #include <ATK/Preamplifier/DempwolfTriodeFunction.h>
       6             : #include <ATK/Preamplifier/EnhancedKorenTriodeFunction.h>
       7             : #include <ATK/Preamplifier/KorenTriodeFunction.h>
       8             : #include <ATK/Preamplifier/LeachTriodeFunction.h>
       9             : #include <ATK/Preamplifier/MunroPiazzaTriodeFunction.h>
      10             : #include <ATK/Preamplifier/ModifiedMunroPiazzaTriodeFunction.h>
      11             : #include <ATK/Preamplifier/Triode2Filter.h>
      12             : #include <ATK/Utility/SimplifiedVectorizedNewtonRaphson.h>
      13             : #include <ATK/Utility/VectorizedNewtonRaphson.h>
      14             : 
      15             : #include <cassert>
      16             : 
      17             : namespace ATK
      18             : {
      19             :   template <typename DataType_, typename TriodeFunction>
      20             :   class Triode2Filter<DataType_, TriodeFunction>::CommonCathodeTriode2Function
      21             :   {
      22             :     const DataType_ Rp;
      23             :     const DataType_ Rg;
      24             :     const DataType_ Ro;
      25             :     const DataType_ Rk;
      26             :     const DataType_ Vbias;
      27             :     const DataType_ Co;
      28             :     const DataType_ Ck;
      29             :     const DataType_ Cpg;
      30             : 
      31             :     DataType_ ickeq;
      32             :     DataType_ icoeq;
      33             :     DataType_ icpgeq;
      34             : 
      35             :     TriodeFunction& tube_function;
      36             : 
      37             :   public:
      38             :     using DataType = DataType_;
      39             :     using Vector = Eigen::Matrix<DataType, 4, 1>;
      40             :     using Matrix = Eigen::Matrix<DataType, 4, 4>;
      41             : 
      42             :     template<typename T>
      43           3 :     CommonCathodeTriode2Function(DataType dt, DataType Rp, DataType Rg, DataType Ro, DataType Rk, DataType Vbias, DataType Co, DataType Ck, TriodeFunction& tube_function, const T& default_output)
      44           3 :       :Rp(1/Rp), Rg(1/Rg), Ro(1/Ro), Rk(1/Rk), Vbias(Vbias), Co(2 / dt * Co), Ck(2 / dt * Ck), Cpg(2 / dt * tube_function.Cpg), ickeq(2 / dt * Ck * default_output[1]), icoeq(-2 / dt * Co * default_output[2]), icpgeq(Cpg * (default_output[3] - default_output[4]) ), tube_function(tube_function)
      45             :     {
      46           3 :     }
      47             : 
      48       39600 :     Vector estimate(gsl::index i, const DataType* const * ATK_RESTRICT input, DataType* const * ATK_RESTRICT output)
      49             :     {
      50       39600 :       return id_estimate(i, input, output);
      51             :     }
      52             :     
      53       39600 :     Vector id_estimate(gsl::index i, const DataType* const * ATK_RESTRICT input, DataType* const * ATK_RESTRICT output)
      54             :     {
      55       39600 :       Vector y0 = Vector::Zero();
      56      198000 :       for (int j = 0; j < 4; ++j)
      57             :       {
      58      158400 :         y0.data()[j] = output[j][i - 1];
      59             :       }
      60             :       
      61       39600 :       return y0;
      62             :     }
      63             :     
      64       39600 :     void update_state(gsl::index i, const DataType* const * ATK_RESTRICT input, DataType* const * ATK_RESTRICT output)
      65             :     {
      66       39600 :       ickeq = 2 * Ck * output[1][i] - ickeq;
      67       39600 :       icoeq = -2 * Co * output[2][i] - icoeq;
      68       39600 :       icpgeq = 2 * Cpg * (output[3][i] - output[4][i]) - icpgeq;
      69       39600 :     }
      70             : 
      71      102497 :     Vector operator()(gsl::index i, const DataType* const * ATK_RESTRICT input, DataType* const * ATK_RESTRICT output, const Vector& y1)
      72             :     {
      73      102497 :       auto Ib = tube_function.Lb(y1(3) - y1(0), y1(2) - y1(0));
      74      102497 :       auto Ic = tube_function.Lc(y1(3) - y1(0), y1(2) - y1(0));
      75             : 
      76      102497 :       auto Ib_Vbe = tube_function.Lb_Vbe(y1(3) - y1(0), y1(2) - y1(0));
      77      102497 :       auto Ib_Vce = tube_function.Lb_Vce(y1(3) - y1(0), y1(2) - y1(0));
      78             : 
      79      102497 :       auto Ic_Vbe = tube_function.Lc_Vbe(y1(3) - y1(0), y1(2) - y1(0));
      80      102497 :       auto Ic_Vce = tube_function.Lc_Vce(y1(3) - y1(0), y1(2) - y1(0));
      81             : 
      82      102497 :       auto f1 = Ib + Ic + ickeq - y1(0) * (Rk + Ck);
      83      102497 :       auto f2 = icoeq + (y1(1) + y1(2)) * Ro + y1(1) * Co;
      84             : 
      85      102497 :       auto g1 = (y1(2) - Vbias) * Rp + (Ic + (y1(1) + y1(2)) * Ro) + (y1(2) - y1(3)) * Cpg - icpgeq;
      86      102497 :       auto g2 = (y1(3) - input[0][i]) * Rg + Ib + icpgeq - (y1(2) - y1(3)) * Cpg;
      87             :       
      88      102497 :       Vector F(f1,
      89             :            f2,
      90             :            g1,
      91             :            g2);
      92             : 
      93      102497 :       Matrix M;
      94      204994 :       M << -(Ib_Vbe + Ic_Vbe + Ib_Vce + Ic_Vce) - (Rk + Ck), 0, (Ib_Vce + Ic_Vce), (Ib_Vbe + Ic_Vbe),
      95      102497 :             0, Ro + Co, Ro, 0,
      96      102497 :       -(Ic_Vbe + Ic_Vce), Ro, Rp + Ro + Ic_Vce + Cpg, Ic_Vbe - Cpg,
      97      204994 :             -(Ib_Vbe + Ib_Vce), 0, Ib_Vce - Cpg, Ib_Vbe + Rg + Cpg;
      98             : 
      99      204994 :       return M.inverse() * F;
     100             :     }
     101             : 
     102             :   };
     103             :   
     104             :   template <typename DataType_, typename TriodeFunction>
     105             :   class CommonCathodeTriode2InitialFunction
     106             :   {
     107             :     const DataType_ Rp;
     108             :     const DataType_ Rg;
     109             :     const DataType_ Ro;
     110             :     const DataType_ Rk;
     111             :     const DataType_ Vbias;
     112             : 
     113             :     TriodeFunction& tube_function;
     114             : 
     115             :   public:
     116             :     using DataType = DataType_;
     117             :     using Vector = Eigen::Matrix<DataType, 3, 1>;
     118             :     using Matrix = Eigen::Matrix<DataType, 3, 3>;
     119             : 
     120           3 :     CommonCathodeTriode2InitialFunction(DataType Rp, DataType Rg, DataType Ro, DataType Rk, DataType Vbias, TriodeFunction& tube_function)
     121           3 :       :Rp(Rp), Rg(Rg), Ro(Ro), Rk(Rk), Vbias(Vbias), tube_function(tube_function)
     122             :     {
     123           3 :     }
     124             : 
     125          42 :     Vector operator()(const Vector& y1)
     126             :     {
     127          42 :       auto Ib = tube_function.Lb(y1(1) - y1(2), y1(0) - y1(2));
     128          42 :       auto Ic = tube_function.Lc(y1(1) - y1(2), y1(0) - y1(2));
     129             : 
     130          42 :       auto Ib_Vbe = tube_function.Lb_Vbe(y1(1) - y1(2), y1(0) - y1(2));
     131          42 :       auto Ib_Vce = tube_function.Lb_Vce(y1(1) - y1(2), y1(0) - y1(2));
     132             : 
     133          42 :       auto Ic_Vbe = tube_function.Lc_Vbe(y1(1) - y1(2), y1(0) - y1(2));
     134          42 :       auto Ic_Vce = tube_function.Lc_Vce(y1(1) - y1(2), y1(0) - y1(2));
     135             : 
     136          42 :       Vector F(y1(0) - Vbias + Ic * Rp,
     137          42 :         Ib * Rg + y1(1),
     138          42 :         y1(2) - (Ib + Ic) * Rk);
     139             : 
     140          42 :       Matrix M;
     141          84 :       M << 1 + Rp * Ic_Vce, Rp * Ic_Vbe, -Rp * (Ic_Vbe + Ic_Vce),
     142          42 :         Ib_Vce * Rg, 1 + Rg * Ib_Vbe, -Rg * (Ib_Vbe + Ib_Vce),
     143          84 :         -(Ic_Vce + Ib_Vce) * Rk, -(Ic_Vbe + Ib_Vbe) * Rk, 1 + (Ic_Vbe + Ic_Vce + Ib_Vbe + Ib_Vce) * Rk;
     144             : 
     145          84 :       return M.inverse() * F;
     146             :     }
     147             :   };
     148             : 
     149             :   template <typename DataType, typename TriodeFunction>
     150           2 :   Triode2Filter<DataType, TriodeFunction>::Triode2Filter(DataType Rp, DataType Rg, DataType Ro, DataType Rk, DataType Vbias, DataType Co, DataType Ck, TriodeFunction&& tube_function)
     151           2 :   :Parent(1, 5), Rp(Rp), Rg(Rg), Ro(Ro), Rk(Rk), Vbias(Vbias), Co(Co), Ck(Ck), tube_function(std::move(tube_function))
     152             :   {
     153           2 :     input_delay = output_delay = 1;
     154           2 :   }
     155             : 
     156             :   template <typename DataType, typename TriodeFunction>
     157           0 :   Triode2Filter<DataType, TriodeFunction>::Triode2Filter(Triode2Filter&& other)
     158           0 :   :Parent(std::move(other)), Rp(other.Rp), Rg(other.Rg), Ro(other.Ro), Rk(other.Rk), Vbias(other.Vbias), Co(other.Co), Ck(other.Ck), tube_function(std::move(other.tube_function))
     159             :   {
     160           0 :   }
     161             : 
     162             :   template<typename DataType,  typename TriodeFunction>
     163           2 :   Triode2Filter<DataType, TriodeFunction>::~Triode2Filter()
     164             :   {
     165           2 :   }
     166             : 
     167             :   template<typename DataType,  typename TriodeFunction>
     168           3 :   void Triode2Filter<DataType, TriodeFunction>::setup()
     169             :   {
     170           3 :     Parent::setup();
     171           3 :     optimizer = std::make_unique<VectorizedNewtonRaphson<CommonCathodeTriode2Function, 4, iterations, true>>(CommonCathodeTriode2Function(static_cast<DataType>(1. / input_sampling_rate),
     172           3 :       Rp, Rg, Ro, Rk, //R
     173           3 :       Vbias, // Vbias
     174           3 :       Co, Ck, // C
     175           3 :       tube_function, // tube
     176           3 :       default_output));
     177           3 :   }
     178             : 
     179             :   template<typename DataType, typename TriodeFunction>
     180           3 :   void Triode2Filter<DataType, TriodeFunction>::full_setup()
     181             :   {
     182           3 :     Eigen::Matrix<DataType, 3, 1> y0;
     183           3 :     y0 << Vbias, 0, 0;
     184             :     
     185             :     // setup default_output
     186           3 :     SimplifiedVectorizedNewtonRaphson<CommonCathodeTriode2InitialFunction<DataType, TriodeFunction>, 3, 20> custom(CommonCathodeTriode2InitialFunction<DataType, TriodeFunction>(
     187           3 :       Rp, Rg, Ro, Rk, //R
     188           3 :       Vbias, // Vbias
     189           3 :       tube_function // tube
     190           3 :       ), std::move(y0));
     191             : 
     192           3 :     auto stable = custom.optimize();
     193             : 
     194           3 :     default_output[0] = 0;
     195           3 :     default_output[1] = stable(2);
     196           3 :     default_output[2] = -stable(0);
     197           3 :     default_output[3] = stable(0);
     198           3 :     default_output[4] = stable(1);
     199             : 
     200           3 :     Parent::full_setup();
     201           3 :   }
     202             : 
     203             :   template<typename DataType, typename TriodeFunction>
     204           2 :   void Triode2Filter<DataType, TriodeFunction>::process_impl(gsl::index size) const
     205             :   {
     206           2 :     assert(input_sampling_rate == output_sampling_rate);
     207             : 
     208       39602 :     for(gsl::index i = 0; i < size; ++i)
     209             :     {
     210       39600 :       optimizer->optimize(i, converted_inputs.data(), outputs.data() + 1);
     211       39600 :       outputs[0][i] = outputs[2][i] + outputs[3][i];
     212       39600 :       optimizer->get_function().update_state(i, converted_inputs.data(), outputs.data());
     213             :     }
     214           2 :   }
     215             : 
     216             :   template<typename DataType, typename TriodeFunction>
     217           2 :   Triode2Filter<DataType, TriodeFunction> Triode2Filter<DataType, TriodeFunction>::build_standard_filter(DataType Rp, DataType Rg, DataType Ro, DataType Rk, DataType Vbias, DataType Co, DataType Ck, TriodeFunction function)
     218             :   {
     219             :     return Triode2Filter<DataType, TriodeFunction>(Rp, Rg, Ro, Rk, //R
     220             :                                                   Vbias, // Vbias
     221             :                                                   Co, Ck, // C
     222           2 :                                                   std::move(function) // tube
     223           2 :                                                   );
     224             :   }
     225             : 
     226             : #if ATK_ENABLE_INSTANTIATION
     227             :   template class Triode2Filter<float, LeachTriodeFunction<float> >;
     228             :   template class Triode2Filter<float, KorenTriodeFunction<float> >;
     229             :   template class Triode2Filter<float, EnhancedKorenTriodeFunction<float> >;
     230             :   template class Triode2Filter<float, MunroPiazzaTriodeFunction<float> >;
     231             :   template class Triode2Filter<float, ModifiedMunroPiazzaTriodeFunction<float> >;
     232             :   template class Triode2Filter<float, DempwolfTriodeFunction<float> >;
     233             : #endif
     234             :   template class Triode2Filter<double, LeachTriodeFunction<double> >;
     235             :   template class Triode2Filter<double, KorenTriodeFunction<double> >;
     236             :   template class Triode2Filter<double, EnhancedKorenTriodeFunction<double> >;
     237             :   template class Triode2Filter<double, MunroPiazzaTriodeFunction<double> >;
     238             :   template class Triode2Filter<double, ModifiedMunroPiazzaTriodeFunction<double> >;
     239             :   template class Triode2Filter<double, DempwolfTriodeFunction<double> >;
     240             : }

Generated by: LCOV version TK-3.3.0-4-gdba42eea