LCOV - code coverage report
Current view: top level - Preamplifier - FollowerTransistorClassAFilter.cpp (source / functions) Hit Total Coverage
Test: coverage.info.cleaned Lines: 128 131 97.7 %
Date: 2021-02-18 20:07:22 Functions: 13 15 86.7 %

          Line data    Source code
       1             : /**
       2             : * \file FollowerTransistorClassAFilter.cpp
       3             : */
       4             : 
       5             : #include <ATK/Preamplifier/FollowerTransistorClassAFilter.h>
       6             : #include <ATK/Preamplifier/TransistorFunction.h>
       7             : #include <ATK/Utility/fmath.h>
       8             : #include <ATK/Utility/SimplifiedVectorizedNewtonRaphson.h>
       9             : #include <ATK/Utility/VectorizedNewtonRaphson.h>
      10             : 
      11             : namespace
      12             : {
      13             :   static constexpr size_t vector_size = 5;
      14             : }
      15             : 
      16             : namespace ATK
      17             : {
      18             :   template <typename DataType_>
      19             :   class FollowerTransistorClassAFilter<DataType_>::TransistorClassAFunction
      20             :   {
      21             :     const DataType_ Rp;
      22             :     const DataType_ Rg1;
      23             :     const DataType_ Rg2;
      24             :     const DataType_ Ro;
      25             :     const DataType_ Rk1;
      26             :     const DataType_ Rk2;
      27             :     const DataType_ Vbias;
      28             :     const DataType_ Cg;
      29             :     const DataType_ Co;
      30             :     const DataType_ Ck;
      31             : 
      32             :     DataType_ ickeq;
      33             :     DataType_ icgeq;
      34             :     DataType_ icoeq;
      35             : 
      36             :     TransistorFunction<DataType_>& transistor_function_1;
      37             :     TransistorFunction<DataType_>& transistor_function_2;
      38             :   public:
      39             :     using DataType = DataType_;
      40             :     using Vector = Eigen::Matrix<DataType, vector_size, 1>;
      41             :     using Matrix = Eigen::Matrix<DataType, vector_size, vector_size>;
      42             : 
      43             :     std::pair<DataType, DataType> exp_y0;
      44             :   
      45             :     template<typename T>
      46           1 :     TransistorClassAFunction(DataType dt, DataType Rp, DataType Rg1, DataType Rg2, DataType Ro, DataType Rk1, DataType Rk2, DataType Vbias, DataType Cg, DataType Co, DataType Ck, TransistorFunction<DataType_>& transistor_function_1, TransistorFunction<DataType_>& transistor_function_2, const T& default_output)
      47           1 :       :Rp(1 / Rp), Rg1(1 / Rg1), Rg2(1 / Rg2), Ro(1 / Ro), Rk1(1 / Rk1), Rk2(1 / Rk2), Vbias(Vbias), Cg(2 / dt * Cg), Co(2 / dt * Co), Ck(2 / dt * Ck), ickeq(2 / dt * Ck * default_output[1]), icgeq(2 / dt * -Cg * default_output[3]), icoeq(2 / dt * Co * (default_output[0] - default_output[4])), transistor_function_1(transistor_function_1), transistor_function_2(transistor_function_2)
      48             :     {
      49           1 :     }
      50             : 
      51        4800 :     Vector estimate(gsl::index i, const DataType* const * ATK_RESTRICT input, DataType* const * ATK_RESTRICT output)
      52             :     {
      53        4800 :       return id_estimate(i, input, output);
      54             :     }
      55             : 
      56        4800 :     Vector id_estimate(gsl::index i, const DataType* const * ATK_RESTRICT input, DataType* const * ATK_RESTRICT output)
      57             :     {
      58        4800 :       Vector y0 = Vector::Zero();
      59       28800 :       for (int j = 0; j < vector_size; ++j)
      60             :       {
      61       24000 :         y0.data()[j] = output[j][i - 1];
      62             :       }
      63             : 
      64        4800 :       return y0;
      65             :     }
      66             : 
      67        4800 :     void update_state(gsl::index i, const DataType* const * ATK_RESTRICT input, DataType* const * ATK_RESTRICT output)
      68             :     {
      69        4800 :       ickeq = 2 * Ck * output[1][i] - ickeq;
      70        4800 :       icgeq = 2 * Cg * (input[0][i] - output[3][i]) - icgeq;
      71        4800 :       icoeq = 2 * Co * (output[0][i] - output[4][i]) - icoeq;
      72        4800 :     }
      73             : 
      74       15766 :     Vector operator()(gsl::index i, const DataType* const * ATK_RESTRICT input, DataType* const * ATK_RESTRICT output, const Vector& y1)
      75             :     {
      76       15766 :       std::pair<DataType, DataType> exp_y1 = std::make_pair(fmath::exp((y1(3) - y1(1)) / transistor_function_1.Vt), fmath::exp((y1(3) - y1(2)) / transistor_function_1.Vt));
      77             : 
      78       15766 :       auto Ib1 = transistor_function_1.Lb(exp_y1);
      79       15766 :       auto Ic1 = transistor_function_1.Lc(exp_y1);
      80             : 
      81       15766 :       auto Ib1_Vbe = transistor_function_1.Lb_Vbe(exp_y1);
      82       15766 :       auto Ib1_Vbc = transistor_function_1.Lb_Vbc(exp_y1);
      83             : 
      84       15766 :       auto Ic1_Vbe = transistor_function_1.Lc_Vbe(exp_y1);
      85       15766 :       auto Ic1_Vbc = transistor_function_1.Lc_Vbc(exp_y1);
      86             : 
      87       15766 :       std::pair<DataType, DataType> exp_y2 = std::make_pair(fmath::exp((y1(2) - y1(4)) / transistor_function_2.Vt), fmath::exp((y1(2) - Vbias) / transistor_function_2.Vt));
      88             : 
      89       15766 :       auto Ib2 = transistor_function_2.Lb(exp_y2);
      90       15766 :       auto Ic2 = transistor_function_2.Lc(exp_y2);
      91             : 
      92       15766 :       auto Ib2_Vbe = transistor_function_2.Lb_Vbe(exp_y2);
      93       15766 :       auto Ib2_Vbc = transistor_function_2.Lb_Vbc(exp_y2);
      94             : 
      95       15766 :       auto Ic2_Vbe = transistor_function_2.Lc_Vbe(exp_y2);
      96       15766 :       auto Ic2_Vbc = transistor_function_2.Lc_Vbc(exp_y2);
      97             : 
      98       15766 :       auto f1 = Ib1 + icgeq + y1(3) * Rg2 + (y1(3) - Vbias) * Rg1 + (y1(3) - input[0][i]) * Cg;
      99       15766 :       auto f2 = Ib1 + Ic1 + ickeq - y1(1) * (Rk1 + Ck);
     100       15766 :       auto f3 = Ic1 + Ib2 + (y1(2) - Vbias) * Rp;
     101       15766 :       auto f4 = Ib2 + Ic2 - icoeq - y1(4) * (Co + Rk2) + y1(0) * Co;
     102       15766 :       auto f5 = y1(0) * Ro - icoeq + (y1(0) - y1(4)) * Co;
     103             : 
     104       15766 :       auto a10 = -Ib1_Vbe - Ib1_Vbc;
     105       15766 :       auto a20 = Ib1_Vbc;
     106       15766 :       auto a30 = Ib1_Vbc + Rg2 + Rg1 + Cg;
     107       15766 :       auto a11 = -Ib1_Vbe - Ib1_Vbc - Ic1_Vbe - Ic1_Vbc - (Rk1 + Ck);
     108       15766 :       auto a21 = Ib1_Vbc + Ic1_Vbc;
     109       15766 :       auto a31 = Ib1_Vbe + Ic1_Vbe;
     110       15766 :       auto a12 = -Ic1_Vbe;
     111       15766 :       auto a22 = -Ic1_Vbc + Ib2_Vbe + Ib2_Vbc + Rp;
     112       15766 :       auto a32 = Ic1_Vbe;
     113       15766 :       auto a42 = -Ib2_Vbe;
     114       15766 :       auto a03 = Co;
     115       15766 :       auto a23 = Ib2_Vbe + Ib2_Vbc + Ic2_Vbe + Ic2_Vbc;
     116       15766 :       auto a43 = -Co - Rk2 - Ib2_Vbe - Ic2_Vbe;
     117       15766 :       auto a04 = Ro + Co;
     118       15766 :       auto a44 = -Co;
     119             : 
     120       15766 :       auto invdet = 1 / (-a03*a10*a21*a32*a44 + a03*a10*a22*a31*a44 + a03*a11*a20*a32*a44 - a03*a11*a22*a30*a44 - a03*a12*a20*a31*a44 + a03*a12*a21*a30*a44 + a04*a10*a21*a32*a43 - a04*a10*a22*a31*a43 + a04*a10*a23*a31*a42 - a04*a11*a20*a32*a43 + a04*a11*a22*a30*a43 - a04*a11*a23*a30*a42 + a04*a12*a20*a31*a43 - a04*a12*a21*a30*a43);
     121             : 
     122       15766 :       Vector F;
     123       15766 :       F << f1, f2, f3, f4, f5;
     124             : 
     125       15766 :       Matrix M;
     126       31532 :       M << a23*a44*(-a11*a32 + a12*a31), a23*a44*(a10*a32 - a12*a30), a23*a44*(-a10*a31 + a11*a30), a44*(-a10*a21*a32 + a10*a22*a31 + a11*a20*a32 - a11*a22*a30 - a12*a20*a31 + a12*a21*a30), a10*a21*a32*a43 - a10*a22*a31*a43 + a10*a23*a31*a42 - a11*a20*a32*a43 + a11*a22*a30*a43 - a11*a23*a30*a42 + a12*a20*a31*a43 - a12*a21*a30*a43,
     127       15766 :         -a03*a21*a32*a44 + a03*a22*a31*a44 + a04*a21*a32*a43 - a04*a22*a31*a43 + a04*a23*a31*a42, a03*a20*a32*a44 - a03*a22*a30*a44 - a04*a20*a32*a43 + a04*a22*a30*a43 - a04*a23*a30*a42, -a03*a20*a31*a44 + a03*a21*a30*a44 + a04*a20*a31*a43 - a04*a21*a30*a43, a04*a42*(-a20*a31 + a21*a30), a03*a42*(a20*a31 - a21*a30),
     128       15766 :         a03*a11*a32*a44 - a03*a12*a31*a44 - a04*a11*a32*a43 + a04*a12*a31*a43, -a03*a10*a32*a44 + a03*a12*a30*a44 + a04*a10*a32*a43 - a04*a12*a30*a43, a03*a10*a31*a44 - a03*a11*a30*a44 - a04*a10*a31*a43 + a04*a11*a30*a43, a04*a42*(a10*a31 - a11*a30), a03*a42*(-a10*a31 + a11*a30),
     129       15766 :         -a03*a11*a22*a44 + a03*a12*a21*a44 + a04*a11*a22*a43 - a04*a11*a23*a42 - a04*a12*a21*a43, a03*a10*a22*a44 - a03*a12*a20*a44 - a04*a10*a22*a43 + a04*a10*a23*a42 + a04*a12*a20*a43, -a03*a10*a21*a44 + a03*a11*a20*a44 + a04*a10*a21*a43 - a04*a11*a20*a43, a04*a42*(-a10*a21 + a11*a20), a03*a42*(a10*a21 - a11*a20),
     130       31532 :         a04*a23*(a11*a32 - a12*a31), a04*a23*(-a10*a32 + a12*a30), a04*a23*(a10*a31 - a11*a30), a04*(a10*a21*a32 - a10*a22*a31 - a11*a20*a32 + a11*a22*a30 + a12*a20*a31 - a12*a21*a30), a03*(-a10*a21*a32 + a10*a22*a31 + a11*a20*a32 - a11*a22*a30 - a12*a20*a31 + a12*a21*a30);
     131             : 
     132       15766 :       return M * F * invdet;
     133             :     }
     134             :   };
     135             : 
     136             :   template <typename DataType_>
     137             :   class TransistorClassAInitialFunction
     138             :   {
     139             :     const DataType_ Rp;
     140             :     const DataType_ Rg1;
     141             :     const DataType_ Rg2;
     142             :     const DataType_ Ro;
     143             :     const DataType_ Rk1;
     144             :     const DataType_ Rk2;
     145             :     const DataType_ Vbias;
     146             : 
     147             :     TransistorFunction<DataType_>& transistor_function_1;
     148             :     TransistorFunction<DataType_>& transistor_function_2;
     149             :   public:
     150             :     using DataType = DataType_;
     151             :     using Vector = Eigen::Matrix<DataType, 4, 1>;
     152             :     using Matrix = Eigen::Matrix<DataType, 4, 4>;
     153             : 
     154           1 :     TransistorClassAInitialFunction(DataType Rp, DataType Rg1, DataType Rg2, DataType Ro, DataType Rk1, DataType Rk2, DataType Vbias, TransistorFunction<DataType_>& transistor_function_1, TransistorFunction<DataType_>& transistor_function_2)
     155           1 :       :Rp(1 / Rp), Rg1(1 / Rg1), Rg2(1 / Rg2), Ro(1 / Ro), Rk1(1 / Rk1), Rk2(1 / Rk2), Vbias(Vbias), transistor_function_1(transistor_function_1), transistor_function_2(transistor_function_2)
     156             :     {
     157           1 :     }
     158             : 
     159          21 :     Vector operator()(const Vector& y1)
     160             :     {
     161          21 :       std::pair<DataType, DataType> exp_y1 = std::make_pair(fmath::exp((y1(2) - y1(0)) / transistor_function_1.Vt), fmath::exp((y1(2) - y1(1)) / transistor_function_1.Vt));
     162             : 
     163          21 :       auto Ib1 = transistor_function_1.Lb(exp_y1);
     164          21 :       auto Ic1 = transistor_function_1.Lc(exp_y1);
     165             : 
     166          21 :       auto Ib1_Vbe = transistor_function_1.Lb_Vbe(exp_y1);
     167          21 :       auto Ib1_Vbc = transistor_function_1.Lb_Vbc(exp_y1);
     168             : 
     169          21 :       auto Ic1_Vbe = transistor_function_1.Lc_Vbe(exp_y1);
     170          21 :       auto Ic1_Vbc = transistor_function_1.Lc_Vbc(exp_y1);
     171             : 
     172          21 :       std::pair<DataType, DataType> exp_y2 = std::make_pair(fmath::exp((y1(1) - y1(3)) / transistor_function_2.Vt), fmath::exp((y1(1) - Vbias) / transistor_function_2.Vt));
     173             : 
     174          21 :       auto Ib2 = transistor_function_2.Lb(exp_y2);
     175          21 :       auto Ic2 = transistor_function_2.Lc(exp_y2);
     176             : 
     177          21 :       auto Ib2_Vbe = transistor_function_2.Lb_Vbe(exp_y2);
     178          21 :       auto Ib2_Vbc = transistor_function_2.Lb_Vbc(exp_y2);
     179             : 
     180          21 :       auto Ic2_Vbe = transistor_function_2.Lc_Vbe(exp_y2);
     181          21 :       auto Ic2_Vbc = transistor_function_2.Lc_Vbc(exp_y2);
     182             : 
     183          21 :       auto f1 = Ib1 + Ic1 - y1(0) * Rk1;
     184          21 :       auto f2 = Ic1 + Ib2 + (y1(1) - Vbias) * Rp;
     185          21 :       auto f3 = Ib1 + y1(2) * Rg2 + (y1(2) - Vbias) * Rg1;
     186          21 :       auto f4 = Ib2 + Ic2 - y1(3) * Rk2;
     187             : 
     188          21 :       Vector F(f1,
     189             :         f2,
     190             :         f3,
     191             :         f4);
     192             : 
     193          21 :       Matrix M;
     194          42 :       M << -Ib1_Vbe - Ic1_Vbe - Rk1, -Ib1_Vbc - Ic1_Vbc, Ib1_Vbe + Ib1_Vbc + Ic1_Vbe + Ic1_Vbc, 0,
     195          21 :         -Ic1_Vbe, -Ic1_Vbc + Ib2_Vbe + Ib2_Vbc + Rp, Ic1_Vbe + Ic1_Vbc, -Ib2_Vbe,
     196          21 :         -Ib1_Vbe, -Ib1_Vbc, Ib1_Vbc + Ib1_Vbe + Rg2 + Rg1, 0,
     197          42 :         0, Ib2_Vbe + Ib2_Vbc + Ic2_Vbe + Ic2_Vbc, 0, -Rk2 - Ib2_Vbe - Ic2_Vbe;
     198             : 
     199          42 :       return M.inverse() * F;
     200             :     }
     201             :   };
     202             : 
     203             :   template <typename DataType>
     204           1 :   FollowerTransistorClassAFilter<DataType>::FollowerTransistorClassAFilter(DataType Rp, DataType Rg1, DataType Rg2, DataType Ro, DataType Rk1, DataType Rk2, DataType Vbias, DataType Cg, DataType Co, DataType Ck, TransistorFunction<DataType>&& transistor_function_1, TransistorFunction<DataType>&& transistor_function_2)
     205           1 :     :Parent(1, vector_size), Rp(Rp), Rg1(Rg1), Rg2(Rg2), Ro(Ro), Rk1(Rk1), Rk2(Rk2), Vbias(Vbias), Cg(Cg), Co(Co), Ck(Ck), transistor_function_1(std::move(transistor_function_1)), transistor_function_2(std::move(transistor_function_2))
     206             :   {
     207           1 :     input_delay = output_delay = 1;
     208           1 :   }
     209             : 
     210             :   template <typename DataType>
     211           0 :   FollowerTransistorClassAFilter<DataType>::FollowerTransistorClassAFilter(FollowerTransistorClassAFilter&& other)
     212           0 :     :Parent(std::move(other)), Rp(other.Rp), Rg1(other.Rg1), Rg2(other.Rg2), Ro(other.Ro), Rk1(other.Rk1), Rk2(other.Rk2), Vbias(other.Vbias), Cg(other.Cg), Co(other.Co), Ck(other.Ck), transistor_function_1(std::move(other.transistor_function_1)), transistor_function_2(std::move(other.transistor_function_2))
     213             :   {
     214           0 :   }
     215             : 
     216             :   template <typename DataType>
     217           1 :   FollowerTransistorClassAFilter<DataType>::~FollowerTransistorClassAFilter()
     218             :   {
     219           1 :   }
     220             : 
     221             :   template<typename DataType_>
     222           1 :   void FollowerTransistorClassAFilter<DataType_>::setup()
     223             :   {
     224           1 :     Parent::setup();
     225           1 :     optimizer = std::make_unique<VectorizedNewtonRaphson<TransistorClassAFunction, vector_size, nb_max_iter, true>>(TransistorClassAFunction(static_cast<DataType>(1. / input_sampling_rate),
     226           1 :       Rp, Rg1, Rg2, Ro, Rk1, Rk2, //R
     227           1 :       Vbias, // Vbias
     228           1 :       Cg, Co, Ck, // C
     229           1 :       transistor_function_1, // transistor
     230           1 :       transistor_function_2, // transistor
     231           1 :       default_output));
     232           1 :   }
     233             : 
     234             :   template<typename DataType_>
     235           1 :   void FollowerTransistorClassAFilter<DataType_>::full_setup()
     236             :   {
     237             :     // setup default_output
     238           1 :     SimplifiedVectorizedNewtonRaphson<TransistorClassAInitialFunction<DataType_>, 4, 50> custom(TransistorClassAInitialFunction<DataType_>(
     239           1 :       Rp, Rg1, Rg2, Ro, Rk1, Rk2, //R
     240           1 :       Vbias, // Vbias
     241           1 :       transistor_function_1, // transistor
     242           1 :       transistor_function_2 // transistor
     243           1 :       ), typename TransistorClassAInitialFunction<DataType_>::Vector(Vbias / 20, Vbias * 6 / 10, Vbias / 10, Vbias / 2));
     244             : 
     245           1 :     auto stable = custom.optimize();
     246             : 
     247           1 :     default_output[0] = 0;
     248           1 :     default_output[1] = stable(0);
     249           1 :     default_output[2] = stable(1);
     250           1 :     default_output[3] = stable(2);
     251           1 :     default_output[4] = stable(3);
     252             : 
     253           1 :     Parent::full_setup();
     254           1 :   }
     255             : 
     256             :   template<typename DataType_>
     257           1 :   void FollowerTransistorClassAFilter<DataType_>::process_impl(gsl::index size) const
     258             :   {
     259           1 :     assert(input_sampling_rate == output_sampling_rate);
     260             : 
     261        4801 :     for (gsl::index i = 0; i < size; ++i)
     262             :     {
     263        4800 :       optimizer->optimize(i, converted_inputs.data(), outputs.data());
     264        4800 :       optimizer->get_function().update_state(i, converted_inputs.data(), outputs.data());
     265             :     }
     266           1 :   }
     267             : 
     268             :   template<typename DataType_>
     269           1 :   FollowerTransistorClassAFilter<DataType_> FollowerTransistorClassAFilter<DataType_>::build_standard_filter(DataType_ Rp, DataType_ Rg1, DataType_ Rg2, DataType_ Ro,
     270             :     DataType_ Rk1, DataType_ Rk2, DataType_ Vbias, DataType_ Cg, DataType_ Co, DataType_ Ck, TransistorFunction<DataType_> function1, TransistorFunction<DataType_> function2)
     271             :   {
     272             :     return FollowerTransistorClassAFilter<DataType_>(Rp, Rg1, Rg2, Ro, Rk1, Rk2, Vbias, Cg, Co, Ck,
     273           1 :       std::move(function1), std::move(function2));
     274             :   }
     275             : 
     276             : #if ATK_ENABLE_INSTANTIATION
     277             :   template class FollowerTransistorClassAFilter<float>;
     278             : #endif
     279             :   template class FollowerTransistorClassAFilter<double>;
     280             : }

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