Line data    Source code

       1             : /**
       2             :  * \file DiodeClipperFilter.cpp
       3             :  */
       4             : 
       5             : #include "DiodeClipperFilter.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 DiodeClipperFilter<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           1 :     SimpleOverdriveFunction(DataType dt, DataType R, DataType C, DataType is, DataType vt)
      32           1 :     :A(dt / (2 * C * R)), B(dt / (2 * C)), is(is), vt(vt)
      33             :     {
      34           1 :     }
      35             :     
      36        6152 :     std::pair<DataType, DataType> operator()(const DataType* ATK_RESTRICT input, DataType* ATK_RESTRICT output, DataType y1)
      37             :     {
      38        6152 :       auto x1 = input[0];
      39        6152 :       auto y0 = output[-1];
      40        6152 :       DataType expdiode_y1_p = fmath::exp(y1 / vt);
      41        6152 :       DataType expdiode_y1_m = 1 / expdiode_y1_p;
      42             :       DataType expdiode_y0_p;
      43             :       DataType expdiode_y0_m;
      44             :           
      45        6152 :             if(y0 == oldy0)
      46             :             {
      47        2161 :               expdiode_y0_p = oldexpy0;
      48        2161 :               expdiode_y0_m = oldinvexpy0;
      49             :             }
      50        3991 :             else if(y0 == oldy1)
      51             :             {
      52         116 :               expdiode_y0_p = oldexpy1;
      53         116 :               expdiode_y0_m = oldinvexpy1;
      54             :             }
      55             :             else
      56             :             {
      57        3875 :               expdiode_y0_p = fmath::exp(y0 / vt);
      58        3875 :               expdiode_y0_m = 1 / expdiode_y0_p;
      59             :             }
      60             :           
      61        6152 :       oldy0 = y0;
      62        6152 :       oldexpy0 = expdiode_y0_p;
      63        6152 :       oldinvexpy0 = expdiode_y0_m;
      64             : 
      65        6152 :       oldy1 = y1;
      66        6152 :       oldexpy1 = expdiode_y1_p;
      67        6152 :       oldinvexpy1 = expdiode_y1_m;
      68             :       
      69        6152 :       std::pair<DataType, DataType> diode = std::make_pair(is * (expdiode_y1_p - expdiode_y1_m), is * (expdiode_y1_p + expdiode_y1_m) / vt);
      70        6152 :       auto old_diode = is * (expdiode_y0_p - expdiode_y0_m);
      71        6152 :       return std::make_pair(y0 - y1 + 2 * A * x1 - (A * (y1 + y0) + B * (diode.first + old_diode)), -1 - A - B * diode.second);
      72             :     }
      73             : 
      74        4000 :     DataType estimate(const DataType* ATK_RESTRICT input, DataType* ATK_RESTRICT output)
      75             :     {
      76        4000 :       auto x0 = input[-1];
      77        4000 :       auto x1 = input[0];
      78        4000 :       auto y0 = output[-1];
      79        4000 :       return affine_estimate(x0, x1, y0);
      80             :     }
      81             : 
      82        4000 :     DataType affine_estimate(DataType x0, DataType x1, DataType y0)
      83             :     {
      84        4000 :       auto sinh = is * (oldexpy1 - oldinvexpy1);
      85        4000 :       auto cosh = is * (oldexpy1 + oldinvexpy1);
      86        4000 :       return (y0 + A * (2 * x1 - y0) - B * (2 * sinh - y0 / vt * cosh) ) / (B * cosh / vt + 1 + A);
      87             :     }
      88             :   };
      89             :   
      90             :   template <typename DataType>
      91           1 :   DiodeClipperFilter<DataType>::DiodeClipperFilter()
      92           1 :   :TypedBaseFilter<DataType>(1, 1)
      93             :   {
      94           1 :     input_delay = 1;
      95           1 :     output_delay = 1;
      96           1 :   }
      97             : 
      98             :   template <typename DataType>
      99           1 :   DiodeClipperFilter<DataType>::~DiodeClipperFilter()
     100             :   {
     101           1 :   }
     102             : 
     103             :   template <typename DataType>
     104           1 :   void DiodeClipperFilter<DataType>::setup()
     105             :   {
     106           1 :     Parent::setup();
     107           1 :     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           1 :   }
     110             : 
     111             :   template <typename DataType>
     112           1 :   void DiodeClipperFilter<DataType>::process_impl(gsl::index size) const
     113             :   {
     114           1 :     const DataType* ATK_RESTRICT input = converted_inputs[0];
     115           1 :     DataType* ATK_RESTRICT output = outputs[0];
     116        4001 :     for(gsl::index i = 0; i < size; ++i)
     117             :     {
     118        4000 :       optimizer->optimize(input + i, output + i);
     119             :     }
     120           1 :   }
     121             : 
     122             : #if ATK_ENABLE_INSTANTIATION
     123             :   template class DiodeClipperFilter<float>;
     124             : #endif
     125             :   template class DiodeClipperFilter<double>;
     126             : }