Line data    Source code

       1             : /**
       2             :  * \file SimpleOverdriveFilter.cpp
       3             :  */
       4             : 
       5             : #include "SimpleOverdriveFilter.h"
       6             : #include <ATK/Utility/fmath.h>
       7             : #include <ATK/Utility/ScalarNewtonRaphson.h>
       8             : 
       9             : #include <boost/math/special_functions/sign.hpp>
      10             : 
      11             : namespace ATK
      12             : {
      13             :   template<typename DataType_>
      14             :   class SimpleOverdriveFilter<DataType_>::SimpleOverdriveFunction
      15             :   {
      16             :   public:
      17             :     using DataType = DataType_;
      18             :   protected:
      19             :     DataType A;
      20             :     DataType B;
      21             :     DataType is;
      22             :     DataType vt;
      23             : 
      24             :     DataType oldy0{0};
      25             :     DataType oldexpy0{1};
      26             :     DataType oldinvexpy0{1};
      27             :     DataType oldy1{0};
      28             :     DataType oldexpy1{1};
      29             :     DataType oldinvexpy1{1};
      30             :   public:
      31           3 :     SimpleOverdriveFunction(DataType dt, DataType R, DataType C, DataType is, DataType vt)
      32           3 :     :A(dt / (2 * C) + R), B(dt / (2 * C) - R), is(is), vt(vt)
      33             :     {
      34           3 :     }
      35             :     
      36        8082 :     std::pair<DataType, DataType> operator()(const DataType* ATK_RESTRICT input, DataType* ATK_RESTRICT output, DataType y1)
      37             :     {
      38        8082 :       auto x0 = input[-1];
      39        8082 :       auto x1 = input[0];
      40        8082 :       auto y0 = output[-1];
      41        8082 :       DataType expdiode_y1_p = fmath::exp(y1 / vt);
      42        8082 :       DataType expdiode_y1_m = 1 / expdiode_y1_p;
      43             :       DataType expdiode_y0_p;
      44             :       DataType expdiode_y0_m;
      45             :           
      46        8082 :             if(y0 == oldy0)
      47             :             {
      48        4091 :               expdiode_y0_p = oldexpy0;
      49        4091 :               expdiode_y0_m = oldinvexpy0;
      50             :             }
      51        3991 :             else if(y0 == oldy1)
      52             :             {
      53          34 :               expdiode_y0_p = oldexpy1;
      54          34 :               expdiode_y0_m = oldinvexpy1;
      55             :             }
      56             :             else
      57             :             {
      58        3957 :               expdiode_y0_p = fmath::exp(y0 / vt);
      59        3957 :               expdiode_y0_m = 1 / expdiode_y0_p;
      60             :             }
      61             :           
      62        8082 :       oldy0 = y0;
      63        8082 :       oldexpy0 = expdiode_y0_p;
      64        8082 :       oldinvexpy0 = expdiode_y0_m;
      65             : 
      66        8082 :       oldy1 = y1;
      67        8082 :       oldexpy1 = expdiode_y1_p;
      68        8082 :       oldinvexpy1 = expdiode_y1_m;
      69             : 
      70        8082 :       std::pair<DataType, DataType> diode = std::make_pair(is * (expdiode_y1_p - expdiode_y1_m), is * (expdiode_y1_p + expdiode_y1_m) / vt);
      71        8082 :       return std::make_pair(A * diode.first + y1 + x0 - x1 + B * is * (expdiode_y0_p - expdiode_y0_m) - y0, A * diode.second + 1);
      72             :     }
      73             : 
      74        5000 :     DataType estimate(const DataType* ATK_RESTRICT input, DataType* ATK_RESTRICT output)
      75             :     {
      76        5000 :       auto x0 = input[-1];
      77        5000 :       auto x1 = input[0];
      78        5000 :       auto y0 = output[-1];
      79        5000 :       return affine_estimate(x0, x1, y0);
      80             :     }
      81             : 
      82        5000 :     DataType affine_estimate(DataType x0, DataType x1, DataType y0)
      83             :     {
      84        5000 :       auto sinh = (oldexpy1 - oldinvexpy1);
      85        5000 :       auto cosh = (oldexpy1 + oldinvexpy1);
      86        5000 :       return (x1 - x0 + y0 - is * sinh * B - (sinh - y0 / vt * cosh) * is * A) / (is * cosh * A / vt + 1);
      87             :     }
      88             :   };
      89             :   
      90             :   template <typename DataType>
      91           2 :   SimpleOverdriveFilter<DataType>::SimpleOverdriveFilter()
      92           2 :   :TypedBaseFilter<DataType>(1, 1)
      93             :   {
      94           2 :     input_delay = 1;
      95           2 :     output_delay = 1;
      96           2 :   }
      97             : 
      98             :   template <typename DataType>
      99           2 :   SimpleOverdriveFilter<DataType>::~SimpleOverdriveFilter()
     100             :   {
     101           2 :   }
     102             : 
     103             :   template <typename DataType>
     104           3 :   void SimpleOverdriveFilter<DataType>::setup()
     105             :   {
     106           3 :     Parent::setup();
     107           3 :     optimizer = std::make_unique<ScalarNewtonRaphson<SimpleOverdriveFunction>>(SimpleOverdriveFunction(static_cast<DataType>(1. / input_sampling_rate),
     108             :       10000, static_cast<DataType>(22e-9), static_cast<DataType>(1e-12), static_cast<DataType>(26e-3)));
     109           3 :   }
     110             : 
     111             :   template <typename DataType>
     112           2 :   void SimpleOverdriveFilter<DataType>::process_impl(gsl::index size) const
     113             :   {
     114           2 :     const DataType* ATK_RESTRICT input = converted_inputs[0];
     115           2 :     DataType* ATK_RESTRICT output = outputs[0];
     116        5002 :     for(gsl::index i = 0; i < size; ++i)
     117             :     {
     118        5000 :       optimizer->optimize(input + i, output + i);
     119             :     }
     120           2 :   }
     121             : 
     122             : #if ATK_ENABLE_INSTANTIATION
     123             :   template class SimpleOverdriveFilter<float>;
     124             : #endif
     125             :   template class SimpleOverdriveFilter<double>;
     126             :   }