Line data Source code
1 : /**
2 : * \file BesselFilter.hxx
3 : */
4 :
5 : #include "BesselFilter.h"
6 : #include <ATK/EQ/helpers.h>
7 :
8 : #include <boost/math/constants/constants.hpp>
9 :
10 : namespace BesselUtilities
11 : {
12 : template<typename DataType>
13 90 : void create_bessel_analog_coefficients(int order, EQUtilities::ZPK<DataType>& zpk)
14 : {
15 90 : zpk.k = 1;
16 90 : auto& z = zpk.z;
17 90 : auto& p = zpk.p;
18 90 : z.clear(); // no zeros for this filter type
19 90 : p.clear();
20 90 : switch (order)
21 : {
22 26 : case 0:
23 26 : break;
24 30 : case 1:
25 30 : p.push_back(1);
26 30 : break;
27 1 : case 2:
28 1 : p.push_back(std::complex<DataType>(-.8660254037844386467637229, .4999999999999999999999996));
29 1 : p.push_back(std::complex<DataType>(-.8660254037844386467637229, -.4999999999999999999999996));
30 1 : break;
31 22 : case 3:
32 22 : p.push_back(std::complex<DataType>(-.9416000265332067855971980, 0));
33 22 : p.push_back(std::complex<DataType>(-.7456403858480766441810907, .7113666249728352680992154));
34 22 : p.push_back(std::complex<DataType>(-.7456403858480766441810907, -.7113666249728352680992154));
35 22 : break;
36 1 : case 4:
37 1 : p.push_back(std::complex<DataType>(-.6572111716718829545787781, .8301614350048733772399715));
38 1 : p.push_back(std::complex<DataType>(-.6572111716718829545787781, -.8301614350048733772399715));
39 1 : p.push_back(std::complex<DataType>(-.9047587967882449459642637, .2709187330038746636700923));
40 1 : p.push_back(std::complex<DataType>(-.9047587967882449459642637, -.2709187330038746636700923));
41 1 : break;
42 1 : case 5:
43 1 : p.push_back(std::complex<DataType>(-.9264420773877602247196260, 0));
44 1 : p.push_back(std::complex<DataType>(-.8515536193688395541722677, .4427174639443327209850002));
45 1 : p.push_back(std::complex<DataType>(-.8515536193688395541722677, -.4427174639443327209850002));
46 1 : p.push_back(std::complex<DataType>(-.5905759446119191779319432, .9072067564574549539291747));
47 1 : p.push_back(std::complex<DataType>(-.5905759446119191779319432, -.9072067564574549539291747));
48 1 : break;
49 1 : case 6:
50 1 : p.push_back(std::complex<DataType>(-.9093906830472271808050953, .1856964396793046769246397));
51 1 : p.push_back(std::complex<DataType>(-.9093906830472271808050953, -.1856964396793046769246397));
52 1 : p.push_back(std::complex<DataType>(-.7996541858328288520243325, .5621717346937317988594118));
53 1 : p.push_back(std::complex<DataType>(-.7996541858328288520243325, -.5621717346937317988594118));
54 1 : p.push_back(std::complex<DataType>(-.5385526816693109683073792, .9616876881954277199245657));
55 1 : p.push_back(std::complex<DataType>(-.5385526816693109683073792, -.9616876881954277199245657));
56 1 : break;
57 1 : case 7:
58 1 : p.push_back(std::complex<DataType>(-.9194871556490290014311619, 0));
59 1 : p.push_back(std::complex<DataType>(-.8800029341523374639772340, .3216652762307739398381830));
60 1 : p.push_back(std::complex<DataType>(-.8800029341523374639772340, -.3216652762307739398381830));
61 1 : p.push_back(std::complex<DataType>(-.7527355434093214462291616, .6504696305522550699212995));
62 1 : p.push_back(std::complex<DataType>(-.7527355434093214462291616, -.6504696305522550699212995));
63 1 : p.push_back(std::complex<DataType>(-.4966917256672316755024763, 1.002508508454420401230220));
64 1 : p.push_back(std::complex<DataType>(-.4966917256672316755024763, -1.002508508454420401230220));
65 1 : break;
66 1 : case 8:
67 1 : p.push_back(std::complex<DataType>(-.9096831546652910216327629, .1412437976671422927888150));
68 1 : p.push_back(std::complex<DataType>(-.9096831546652910216327629, -.1412437976671422927888150));
69 1 : p.push_back(std::complex<DataType>(-.8473250802359334320103023, .4259017538272934994996429));
70 1 : p.push_back(std::complex<DataType>(-.8473250802359334320103023, -.4259017538272934994996429));
71 1 : p.push_back(std::complex<DataType>(-.7111381808485399250796172, .7186517314108401705762571));
72 1 : p.push_back(std::complex<DataType>(-.7111381808485399250796172, -.7186517314108401705762571));
73 1 : p.push_back(std::complex<DataType>(-.4621740412532122027072175, 1.034388681126901058116589));
74 1 : p.push_back(std::complex<DataType>(-.4621740412532122027072175, -1.034388681126901058116589));
75 1 : break;
76 1 : case 9:
77 1 : p.push_back(std::complex<DataType>(-.9154957797499037686769223, 0));
78 1 : p.push_back(std::complex<DataType>(-.8911217017079759323183848, .2526580934582164192308115));
79 1 : p.push_back(std::complex<DataType>(-.8911217017079759323183848, -.2526580934582164192308115));
80 1 : p.push_back(std::complex<DataType>(-.8148021112269012975514135, .5085815689631499483745341));
81 1 : p.push_back(std::complex<DataType>(-.8148021112269012975514135, -.5085815689631499483745341));
82 1 : p.push_back(std::complex<DataType>(-.6743622686854761980403401, .7730546212691183706919682));
83 1 : p.push_back(std::complex<DataType>(-.6743622686854761980403401, -.7730546212691183706919682));
84 1 : p.push_back(std::complex<DataType>(-.4331415561553618854685942, 1.060073670135929666774323));
85 1 : p.push_back(std::complex<DataType>(-.4331415561553618854685942, -1.060073670135929666774323));
86 1 : break;
87 1 : case 10:
88 1 : p.push_back(std::complex<DataType>(-.9091347320900502436826431, .1139583137335511169927714));
89 1 : p.push_back(std::complex<DataType>(-.9091347320900502436826431, -.1139583137335511169927714));
90 1 : p.push_back(std::complex<DataType>(-.8688459641284764527921864, .3430008233766309973110589));
91 1 : p.push_back(std::complex<DataType>(-.8688459641284764527921864, -.3430008233766309973110589));
92 1 : p.push_back(std::complex<DataType>(-.7837694413101441082655890, .5759147538499947070009852));
93 1 : p.push_back(std::complex<DataType>(-.7837694413101441082655890, -.5759147538499947070009852));
94 1 : p.push_back(std::complex<DataType>(-.6417513866988316136190854, .8175836167191017226233947));
95 1 : p.push_back(std::complex<DataType>(-.6417513866988316136190854, -.8175836167191017226233947));
96 1 : p.push_back(std::complex<DataType>(-.4083220732868861566219785, 1.081274842819124562037210));
97 1 : p.push_back(std::complex<DataType>(-.4083220732868861566219785, -1.081274842819124562037210));
98 1 : break;
99 1 : case 11:
100 1 : p.push_back(std::complex<DataType>(-.9129067244518981934637318, 0));
101 1 : p.push_back(std::complex<DataType>(-.8963656705721166099815744, .2080480375071031919692341));
102 1 : p.push_back(std::complex<DataType>(-.8963656705721166099815744, -.2080480375071031919692341));
103 1 : p.push_back(std::complex<DataType>(-.8453044014712962954184557, .4178696917801248292797448));
104 1 : p.push_back(std::complex<DataType>(-.8453044014712962954184557, -.4178696917801248292797448));
105 1 : p.push_back(std::complex<DataType>(-.7546938934722303128102142, .6319150050721846494520941));
106 1 : p.push_back(std::complex<DataType>(-.7546938934722303128102142, -.6319150050721846494520941));
107 1 : p.push_back(std::complex<DataType>(-.6126871554915194054182909, .8547813893314764631518509));
108 1 : p.push_back(std::complex<DataType>(-.6126871554915194054182909, -.8547813893314764631518509));
109 1 : p.push_back(std::complex<DataType>(-.3868149510055090879155425, 1.099117466763120928733632));
110 1 : p.push_back(std::complex<DataType>(-.3868149510055090879155425, -1.099117466763120928733632));
111 1 : break;
112 1 : case 12:
113 1 : p.push_back(std::complex<DataType>(-.9084478234140682638817772, 95506365213450398415258360.0e-27));
114 1 : p.push_back(std::complex<DataType>(-.9084478234140682638817772, -95506365213450398415258360.0e-27));
115 1 : p.push_back(std::complex<DataType>(-.8802534342016826507901575, .2871779503524226723615457));
116 1 : p.push_back(std::complex<DataType>(-.8802534342016826507901575, -.2871779503524226723615457));
117 1 : p.push_back(std::complex<DataType>(-.8217296939939077285792834, .4810212115100676440620548));
118 1 : p.push_back(std::complex<DataType>(-.8217296939939077285792834, -.4810212115100676440620548));
119 1 : p.push_back(std::complex<DataType>(-.7276681615395159454547013, .6792961178764694160048987));
120 1 : p.push_back(std::complex<DataType>(-.7276681615395159454547013, -.6792961178764694160048987));
121 1 : p.push_back(std::complex<DataType>(-.5866369321861477207528215, .8863772751320727026622149));
122 1 : p.push_back(std::complex<DataType>(-.5866369321861477207528215, -.8863772751320727026622149));
123 1 : p.push_back(std::complex<DataType>(-.3679640085526312839425808, 1.114373575641546257595657));
124 1 : p.push_back(std::complex<DataType>(-.3679640085526312839425808, -1.114373575641546257595657));
125 1 : break;
126 1 : case 13:
127 1 : p.push_back(std::complex<DataType>(-.9110914665984182781070663, 0));
128 1 : p.push_back(std::complex<DataType>(-.8991314665475196220910718, .1768342956161043620980863));
129 1 : p.push_back(std::complex<DataType>(-.8991314665475196220910718, -.1768342956161043620980863));
130 1 : p.push_back(std::complex<DataType>(-.8625094198260548711573628, .3547413731172988997754038));
131 1 : p.push_back(std::complex<DataType>(-.8625094198260548711573628, -.3547413731172988997754038));
132 1 : p.push_back(std::complex<DataType>(-.7987460692470972510394686, .5350752120696801938272504));
133 1 : p.push_back(std::complex<DataType>(-.7987460692470972510394686, -.5350752120696801938272504));
134 1 : p.push_back(std::complex<DataType>(-.7026234675721275653944062, .7199611890171304131266374));
135 1 : p.push_back(std::complex<DataType>(-.7026234675721275653944062, -.7199611890171304131266374));
136 1 : p.push_back(std::complex<DataType>(-.5631559842430199266325818, .9135900338325109684927731));
137 1 : p.push_back(std::complex<DataType>(-.5631559842430199266325818, -.9135900338325109684927731));
138 1 : p.push_back(std::complex<DataType>(-.3512792323389821669401925, 1.127591548317705678613239));
139 1 : p.push_back(std::complex<DataType>(-.3512792323389821669401925, -1.127591548317705678613239));
140 1 : break;
141 1 : default:
142 1 : throw std::out_of_range("Can't create a Bessel filter with this order");
143 : }
144 89 : }
145 :
146 : template<typename DataType, typename Container>
147 54 : void create_default_bessel_coeffs(size_t order, DataType Wn, Container& coefficients_in, Container& coefficients_out)
148 : {
149 108 : EQUtilities::ZPK<DataType> zpk;
150 :
151 54 : int fs = 2;
152 54 : create_bessel_analog_coefficients(static_cast<int>(order), zpk);
153 53 : EQUtilities::populate_lp_coeffs(Wn, fs, order, zpk, coefficients_in, coefficients_out);
154 53 : }
155 :
156 : template<typename DataType, typename Container>
157 18 : void create_bp_bessel_coeffs(size_t order, DataType wc1, DataType wc2, Container& coefficients_in, Container& coefficients_out)
158 : {
159 36 : EQUtilities::ZPK<DataType> zpk;
160 :
161 18 : int fs = 2;
162 18 : create_bessel_analog_coefficients(static_cast<int>(order/2), zpk);
163 18 : EQUtilities::populate_bp_coeffs(wc1, wc2, fs, order, zpk, coefficients_in, coefficients_out);
164 18 : }
165 :
166 : template<typename DataType, typename Container>
167 18 : void create_bs_bessel_coeffs(size_t order, DataType wc1, DataType wc2, Container& coefficients_in, Container& coefficients_out)
168 : {
169 36 : EQUtilities::ZPK<DataType> zpk;
170 :
171 18 : int fs = 2;
172 18 : create_bessel_analog_coefficients(static_cast<int>(order/2), zpk);
173 18 : EQUtilities::populate_bs_coeffs(wc1, wc2, fs, order, zpk, coefficients_in, coefficients_out);
174 18 : }
175 : }
176 :
177 : namespace ATK
178 : {
179 : template <typename DataType>
180 9 : BesselLowPassCoefficients<DataType>::BesselLowPassCoefficients(gsl::index nb_channels)
181 9 : :Parent(nb_channels, nb_channels)
182 : {
183 9 : }
184 :
185 : template <typename DataType_>
186 7 : void BesselLowPassCoefficients<DataType_>::set_cut_frequency(CoeffDataType cut_frequency)
187 : {
188 7 : if(cut_frequency <= 0)
189 : {
190 1 : throw std::out_of_range("Frequency can't be negative");
191 : }
192 6 : this->cut_frequency = cut_frequency;
193 6 : setup();
194 6 : }
195 :
196 : template <typename DataType_>
197 1 : typename BesselLowPassCoefficients<DataType_>::CoeffDataType BesselLowPassCoefficients<DataType_>::get_cut_frequency() const
198 : {
199 1 : return cut_frequency;
200 : }
201 :
202 : template <typename DataType>
203 21 : void BesselLowPassCoefficients<DataType>::set_order(unsigned int order)
204 : {
205 21 : if(order == 0)
206 : {
207 1 : throw std::out_of_range("Order can't be null");
208 : }
209 20 : in_order = out_order = order;
210 20 : setup();
211 19 : }
212 :
213 : template <typename DataType>
214 1 : unsigned int BesselLowPassCoefficients<DataType>::get_order() const
215 : {
216 1 : return in_order;
217 : }
218 :
219 : template <typename DataType>
220 36 : void BesselLowPassCoefficients<DataType>::setup()
221 : {
222 36 : Parent::setup();
223 36 : coefficients_in.assign(in_order+1, 0);
224 36 : coefficients_out.assign(out_order, 0);
225 :
226 36 : BesselUtilities::create_default_bessel_coeffs(in_order, 2 * cut_frequency / input_sampling_rate, coefficients_in, coefficients_out);
227 35 : }
228 :
229 : template <typename DataType>
230 8 : BesselHighPassCoefficients<DataType>::BesselHighPassCoefficients(gsl::index nb_channels)
231 8 : :Parent(nb_channels, nb_channels)
232 : {
233 8 : }
234 :
235 : template <typename DataType_>
236 6 : void BesselHighPassCoefficients<DataType_>::set_cut_frequency(CoeffDataType cut_frequency)
237 : {
238 6 : if(cut_frequency <= 0)
239 : {
240 1 : throw std::out_of_range("Frequency can't be negative");
241 : }
242 5 : this->cut_frequency = cut_frequency;
243 5 : setup();
244 5 : }
245 :
246 : template <typename DataType_>
247 1 : typename BesselHighPassCoefficients<DataType_>::CoeffDataType BesselHighPassCoefficients<DataType_>::get_cut_frequency() const
248 : {
249 1 : return cut_frequency;
250 : }
251 :
252 : template <typename DataType>
253 6 : void BesselHighPassCoefficients<DataType>::set_order(unsigned int order)
254 : {
255 6 : if(order == 0)
256 : {
257 1 : throw std::out_of_range("Order can't be null");
258 : }
259 5 : in_order = out_order = order;
260 5 : setup();
261 5 : }
262 :
263 : template <typename DataType>
264 1 : unsigned int BesselHighPassCoefficients<DataType>::get_order() const
265 : {
266 1 : return in_order;
267 : }
268 :
269 : template <typename DataType>
270 18 : void BesselHighPassCoefficients<DataType>::setup()
271 : {
272 18 : Parent::setup();
273 18 : coefficients_in.assign(in_order+1, 0);
274 18 : coefficients_out.assign(out_order, 0);
275 :
276 18 : BesselUtilities::create_default_bessel_coeffs(in_order, (input_sampling_rate - 2 * cut_frequency) / input_sampling_rate, coefficients_in, coefficients_out);
277 41 : for(gsl::index i = in_order - 1; i >= 0; i -= 2)
278 : {
279 23 : coefficients_in[i] = - coefficients_in[i];
280 23 : coefficients_out[i] = - coefficients_out[i];
281 : }
282 18 : }
283 :
284 : template <typename DataType>
285 9 : BesselBandPassCoefficients<DataType>::BesselBandPassCoefficients(gsl::index nb_channels)
286 9 : :Parent(nb_channels, nb_channels)
287 : {
288 9 : }
289 :
290 : template <typename DataType_>
291 7 : void BesselBandPassCoefficients<DataType_>::set_cut_frequencies(std::pair<CoeffDataType, CoeffDataType> cut_frequencies)
292 : {
293 7 : if(cut_frequencies.first <= 0 || cut_frequencies.second <= 0)
294 : {
295 2 : throw std::out_of_range("Frequencies can't be negative");
296 : }
297 5 : this->cut_frequencies = cut_frequencies;
298 5 : setup();
299 5 : }
300 :
301 : template <typename DataType_>
302 3 : void BesselBandPassCoefficients<DataType_>::set_cut_frequencies(CoeffDataType f0, CoeffDataType f1)
303 : {
304 3 : set_cut_frequencies(std::make_pair(f0, f1));
305 1 : }
306 :
307 : template <typename DataType_>
308 2 : std::pair<typename BesselBandPassCoefficients<DataType_>::CoeffDataType, typename BesselBandPassCoefficients<DataType_>::CoeffDataType> BesselBandPassCoefficients<DataType_>::get_cut_frequencies() const
309 : {
310 2 : return cut_frequencies;
311 : }
312 :
313 : template <typename DataType>
314 6 : void BesselBandPassCoefficients<DataType>::set_order(unsigned int order)
315 : {
316 6 : if(order == 0)
317 : {
318 1 : throw std::out_of_range("Order can't be null");
319 : }
320 5 : in_order = out_order = 2 * order;
321 5 : setup();
322 5 : }
323 :
324 : template <typename DataType>
325 1 : unsigned int BesselBandPassCoefficients<DataType>::get_order() const
326 : {
327 1 : return in_order / 2;
328 : }
329 :
330 : template <typename DataType>
331 18 : void BesselBandPassCoefficients<DataType>::setup()
332 : {
333 18 : Parent::setup();
334 18 : coefficients_in.assign(in_order+1, 0);
335 18 : coefficients_out.assign(out_order, 0);
336 :
337 18 : BesselUtilities::create_bp_bessel_coeffs(in_order, 2 * cut_frequencies.first / input_sampling_rate, 2 * cut_frequencies.second / input_sampling_rate, coefficients_in, coefficients_out);
338 18 : }
339 :
340 : template <typename DataType>
341 9 : BesselBandStopCoefficients<DataType>::BesselBandStopCoefficients(gsl::index nb_channels)
342 9 : :Parent(nb_channels, nb_channels)
343 : {
344 9 : }
345 :
346 : template <typename DataType_>
347 7 : void BesselBandStopCoefficients<DataType_>::set_cut_frequencies(std::pair<CoeffDataType, CoeffDataType> cut_frequencies)
348 : {
349 7 : if(cut_frequencies.first <= 0 || cut_frequencies.second <= 0)
350 : {
351 2 : throw std::out_of_range("Frequencies can't be negative");
352 : }
353 5 : this->cut_frequencies = cut_frequencies;
354 5 : setup();
355 5 : }
356 :
357 : template <typename DataType_>
358 3 : void BesselBandStopCoefficients<DataType_>::set_cut_frequencies(CoeffDataType f0, CoeffDataType f1)
359 : {
360 3 : set_cut_frequencies(std::make_pair(f0, f1));
361 1 : }
362 :
363 : template <typename DataType_>
364 2 : std::pair<typename BesselBandStopCoefficients<DataType_>::CoeffDataType, typename BesselBandStopCoefficients<DataType_>::CoeffDataType> BesselBandStopCoefficients<DataType_>::get_cut_frequencies() const
365 : {
366 2 : return cut_frequencies;
367 : }
368 :
369 : template <typename DataType>
370 6 : void BesselBandStopCoefficients<DataType>::set_order(unsigned int order)
371 : {
372 6 : if(order == 0)
373 : {
374 1 : throw std::out_of_range("Order can't be null");
375 : }
376 5 : in_order = out_order = 2 * order;
377 5 : setup();
378 5 : }
379 :
380 : template <typename DataType>
381 1 : unsigned int BesselBandStopCoefficients<DataType>::get_order() const
382 : {
383 1 : return in_order / 2;
384 : }
385 :
386 : template <typename DataType>
387 18 : void BesselBandStopCoefficients<DataType>::setup()
388 : {
389 18 : Parent::setup();
390 18 : coefficients_in.assign(in_order+1, 0);
391 18 : coefficients_out.assign(out_order, 0);
392 :
393 18 : BesselUtilities::create_bs_bessel_coeffs(in_order, 2 * cut_frequencies.first / input_sampling_rate, 2 * cut_frequencies.second / input_sampling_rate, coefficients_in, coefficients_out);
394 18 : }
395 : }
|
