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