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