Line data    Source code

       1             : /**
       2             :  * \file SD1OverdriveFilter.cpp
       3             :  */
       4             : 
       5             : #include "SD1OverdriveFilter.h"
       6             : #include <ATK/Utility/fmath.h>
       7             : #include <ATK/Utility/ScalarNewtonRaphson.h>
       8             : 
       9             : #include <boost/math/special_functions/sign.hpp>
      10             : 
      11             : #include <stdexcept>
      12             : 
      13             : namespace ATK
      14             : {
      15             :   template<typename DataType_>
      16             :   class SD1OverdriveFilter<DataType_>::SD1OverdriveFunction
      17             :   {
      18             :   public:
      19             :     using DataType = DataType_;
      20             :   protected:
      21             :     const DataType R;
      22             :     const DataType R1;
      23             :     const DataType C;
      24             :     const DataType Q;
      25             :     DataType drive = 0.5;
      26             :     const DataType is;
      27             :     const DataType vt;
      28             : 
      29             :     DataType ieq{0};
      30             :     DataType i{0};
      31             : 
      32             :     DataType expdiode_y1_p{1};
      33             :     DataType expdiode_y1_m{1};
      34             : 
      35             :   public:
      36           1 :     SD1OverdriveFunction(DataType dt, DataType R, DataType C, DataType R1, DataType Q, DataType is, DataType vt)
      37           1 :       :R(R), R1(R1), C(2 * C / dt), Q(Q), is(is), vt(vt)
      38             :     {
      39           1 :     }
      40             :     
      41           2 :     void set_drive(DataType drive)
      42             :     {
      43           2 :       this->drive = (R1 + drive * Q);
      44           2 :     }
      45             : 
      46        7751 :     std::pair<DataType, DataType> operator()(const DataType* ATK_RESTRICT input, DataType* ATK_RESTRICT output, DataType y1)
      47             :     {
      48        7751 :       auto x1 = input[0];
      49        7751 :       y1 -= x1;
      50        7751 :       expdiode_y1_p = fmath::exp(y1 / vt);
      51        7751 :       expdiode_y1_m = 1 / expdiode_y1_p;
      52             : 
      53        7751 :       DataType diode1 = is * (expdiode_y1_p - 2 * expdiode_y1_m + 1);
      54        7751 :       DataType diode1_derivative = is * (expdiode_y1_p + 2 * expdiode_y1_m) / vt;
      55             : 
      56        7751 :       i = (C * x1 - ieq) / (1 + R * C);
      57             : 
      58        7751 :       return std::make_pair(y1 / drive + diode1 - i, 1 / drive + diode1_derivative);
      59             :     }
      60             :     
      61        4000 :     void update_state(const DataType* ATK_RESTRICT input, DataType* ATK_RESTRICT output)
      62             :     {
      63        4000 :       auto x1 = input[0];
      64        4000 :       ieq = 2 * C * (x1 - i * R) - ieq;
      65        4000 :     }
      66             : 
      67        4000 :     DataType estimate(const DataType* ATK_RESTRICT input, DataType* ATK_RESTRICT output)
      68             :     {
      69        4000 :       auto x0 = input[-1];
      70        4000 :       auto x1 = input[0];
      71        4000 :       auto y0 = output[-1];
      72        4000 :       return affine_estimate(x0, x1, y0);
      73             :     }
      74             : 
      75        4000 :     DataType affine_estimate(DataType x0, DataType x1, DataType y0)
      76             :     {
      77        4000 :       y0 -= x0;
      78        4000 :       auto sinh = is * (expdiode_y1_p - 2 * expdiode_y1_m + 1);
      79        4000 :       auto cosh = is * (expdiode_y1_p + 2 * expdiode_y1_m);
      80        4000 :       auto i = (C * x1 - ieq) / (1 + R * C);
      81             : 
      82        4000 :       return (i - (sinh - y0 / vt * cosh)) / (cosh / vt + (1 / drive)) + x1;
      83             :     }
      84             :   };
      85             :   
      86             :   template <typename DataType>
      87           4 :   SD1OverdriveFilter<DataType>::SD1OverdriveFilter()
      88           4 :     :TypedBaseFilter<DataType>(1, 1)
      89             :   {
      90           4 :     input_delay = 1;
      91           4 :     output_delay = 1;
      92           4 :   }
      93             : 
      94             :   template <typename DataType>
      95           4 :   SD1OverdriveFilter<DataType>::~SD1OverdriveFilter()
      96             :   {
      97           4 :   }
      98             : 
      99             :   template <typename DataType>
     100           1 :   void SD1OverdriveFilter<DataType>::setup()
     101             :   {
     102           1 :     Parent::setup();
     103           1 :     optimizer = std::make_unique<ScalarNewtonRaphson<SD1OverdriveFunction, num_iterations, true>>(SD1OverdriveFunction(static_cast<DataType>(1. / input_sampling_rate),
     104             :       static_cast<DataType>(4.7e3), static_cast<DataType>(0.047e-6), static_cast<DataType>(33e3),
     105             :       static_cast<DataType>(1e6), static_cast<DataType>(1e-12), static_cast<DataType>(26e-3)));
     106             : 
     107           1 :     optimizer->get_function().set_drive(drive);
     108           1 :   }
     109             : 
     110             :   template <typename DataType_>
     111           4 :   void SD1OverdriveFilter<DataType_>::set_drive(DataType_ drive)
     112             :   {
     113           4 :     if(drive < 0 || drive > 1)
     114             :     {
     115           2 :       throw std::out_of_range("Drive must be a value between 0 and 1");
     116             :     }
     117           2 :     this->drive = drive;
     118           2 :     if(optimizer)
     119             :     {
     120           1 :         optimizer->get_function().set_drive(drive);
     121             :     }
     122           2 :   }
     123             : 
     124             :   template <typename DataType_>
     125           1 :   DataType_ SD1OverdriveFilter<DataType_>::get_drive() const
     126             :   {
     127           1 :     return drive;
     128             :   }
     129             : 
     130             :   template <typename DataType>
     131           1 :   void SD1OverdriveFilter<DataType>::process_impl(gsl::index size) const
     132             :   {
     133           1 :     const DataType* ATK_RESTRICT input = converted_inputs[0];
     134           1 :     DataType* ATK_RESTRICT output = outputs[0];
     135        4001 :     for(gsl::index i = 0; i < size; ++i)
     136             :     {
     137        4000 :       optimizer->optimize(input + i, output + i);
     138        4000 :       optimizer->get_function().update_state(input + i, output + i);
     139             :     }
     140           1 :   }
     141             : 
     142             : #if ATK_ENABLE_INSTANTIATION
     143             :   template class SD1OverdriveFilter<float>;
     144             : #endif
     145             :   template class SD1OverdriveFilter<double>;
     146             : }