1 #include "FilterCoefficientMaker.h"
2 #include "SurgeStorage.h"
3 #include <vt_dsp/basic_dsp.h>
4 
5 #include "filters/VintageLadders.h"
6 #include "filters/Obxd.h"
7 #include "filters/K35.h"
8 #include "filters/DiodeLadder.h"
9 #include "filters/NonlinearFeedback.h"
10 #include "filters/NonlinearStates.h"
11 
12 #include "DebugHelpers.h"
13 
14 using namespace std;
15 
16 const float smooth = 0.2f;
17 
FilterCoefficientMaker()18 FilterCoefficientMaker::FilterCoefficientMaker() { Reset(); }
19 
MakeCoeffs(float Freq,float Reso,int Type,int SubType,SurgeStorage * storageI,bool tuningAdjusted)20 void FilterCoefficientMaker::MakeCoeffs(float Freq, float Reso, int Type, int SubType,
21                                         SurgeStorage *storageI, bool tuningAdjusted)
22 {
23     storage = storageI;
24     if (storage)
25     {
26         if (tuningAdjusted && storage->tuningApplicationMode == SurgeStorage::RETUNE_ALL)
27         {
28             /*
29              * Modulations are not remapped and tuning is in efffect; remap the note
30              */
31             auto idx = (int)floor(Freq + 69);
32             float frac = (Freq + 69) - idx; // frac is 0 means use idx; frac is 1 means use idx+1
33 
34             float b0 = storage->currentTuning.logScaledFrequencyForMidiNote(idx) * 12;
35             float b1 = storage->currentTuning.logScaledFrequencyForMidiNote(idx + 1) * 12;
36 
37             auto q = (1.f - frac) * b0 + frac * b1;
38 
39             Freq = q - 69;
40         }
41     }
42     // Force compiler to error out if I miss one
43     fu_type fType = (fu_type)Type;
44 
45     switch (fType)
46     {
47     case fut_lp12:
48         if (SubType == st_SVF)
49             Coeff_SVF(Freq, Reso, false);
50         else
51             Coeff_LP12(Freq, Reso, SubType);
52         break;
53     case fut_hp12:
54         if (SubType == st_SVF)
55             Coeff_SVF(Freq, Reso, false);
56         else
57             Coeff_HP12(Freq, Reso, SubType);
58         break;
59     case fut_bp12:
60         if (SubType == st_SVF)
61             Coeff_SVF(Freq, Reso, false);
62         else
63             Coeff_BP12(Freq, Reso, SubType);
64     case fut_bp24:
65         if (SubType == st_SVF)
66             Coeff_SVF(Freq, Reso, false); // WHY FALSE? It was this way before #3006 tho
67         else
68             Coeff_BP24(Freq, Reso, SubType);
69         break;
70     case fut_notch12:
71     case fut_notch24:
72         Coeff_Notch(Freq, Reso, SubType);
73         break;
74     case fut_apf:
75         Coeff_APF(Freq, Reso, SubType);
76         break;
77     case fut_lp24:
78         if (SubType == st_SVF)
79             Coeff_SVF(Freq, Reso, true);
80         else
81             Coeff_LP24(Freq, Reso, SubType);
82         break;
83     case fut_hp24:
84         if (SubType == st_SVF)
85             Coeff_SVF(Freq, Reso, true);
86         else
87             Coeff_HP24(Freq, Reso, SubType);
88         break;
89     case fut_lpmoog:
90         Coeff_LP4L(Freq, Reso, SubType);
91         break;
92     case fut_comb_pos:
93         Coeff_COMB(Freq, Reso, SubType);
94         break;
95     case fut_comb_neg:
96         Coeff_COMB(Freq, Reso, SubType + 2); // -ve feedback is the next 2 subtypes
97         break;
98     case fut_SNH:
99         Coeff_SNH(Freq, Reso, SubType);
100         break;
101     case fut_vintageladder:
102         switch (SubType)
103         {
104         case 0:
105         case 1:
106             VintageLadder::RK::makeCoefficients(this, Freq, Reso, SubType == 1, storageI);
107             break;
108         case 2:
109         case 3:
110             VintageLadder::Huov::makeCoefficients(this, Freq, Reso, SubType == 3, storageI);
111             break;
112         default:
113             // SOFTWARE ERROR
114             break;
115         }
116         break;
117         /* When we split OBXD we went from one filter with 8 settings (L, B, H, N, L+, B+, H+, N+)
118          * to two subtypes (regular and +). So what used to be Highpass (value 2 and 6) is now 1 and
119          * 2 on a new type. But don't rewrite the filter for now. Just reconstruct the subtypes
120          * */
121     case fut_obxd_2pole_lp:
122         ObxdFilter::makeCoefficients(this, ObxdFilter::TWO_POLE, Freq, Reso, SubType * 4, storageI);
123         break;
124     case fut_obxd_2pole_bp:
125         ObxdFilter::makeCoefficients(this, ObxdFilter::TWO_POLE, Freq, Reso, SubType * 4 + 1,
126                                      storageI);
127         break;
128     case fut_obxd_2pole_hp:
129         ObxdFilter::makeCoefficients(this, ObxdFilter::TWO_POLE, Freq, Reso, SubType * 4 + 2,
130                                      storageI);
131         break;
132     case fut_obxd_2pole_n:
133         ObxdFilter::makeCoefficients(this, ObxdFilter::TWO_POLE, Freq, Reso, SubType * 4 + 3,
134                                      storageI);
135         break;
136     case fut_obxd_4pole:
137         ObxdFilter::makeCoefficients(this, ObxdFilter::FOUR_POLE, Freq, Reso, SubType, storageI);
138         break;
139     case fut_k35_lp:
140         K35Filter::makeCoefficients(this, Freq, Reso, true, fut_k35_saturations[SubType], storageI);
141         break;
142     case fut_k35_hp:
143         K35Filter::makeCoefficients(this, Freq, Reso, false, fut_k35_saturations[SubType],
144                                     storageI);
145         break;
146     case fut_diode:
147         DiodeLadderFilter::makeCoefficients(this, Freq, Reso, storageI);
148         break;
149     case fut_cutoffwarp_lp:
150     case fut_cutoffwarp_hp:
151     case fut_cutoffwarp_n:
152     case fut_cutoffwarp_bp:
153     case fut_cutoffwarp_ap:
154         NonlinearFeedbackFilter::makeCoefficients(this, Freq, Reso, Type, SubType, storageI);
155         break;
156     case fut_resonancewarp_lp:
157     case fut_resonancewarp_hp:
158     case fut_resonancewarp_n:
159     case fut_resonancewarp_bp:
160     case fut_resonancewarp_ap:
161         NonlinearStatesFilter::makeCoefficients(this, Freq, Reso, Type, storageI);
162         break;
163 
164     case n_fu_types:
165         // This should really be an error condition of course
166     case fut_none:
167         break;
168     };
169 }
170 
clipscale(float freq,int subtype)171 float clipscale(float freq, int subtype)
172 {
173     switch (subtype)
174     {
175     case st_Rough:
176         return (1.0f / 64.0f) * db_to_linear(freq * 0.55f);
177     case st_Smooth:
178         return (1.0f / 1024.0f);
179     };
180     return 0;
181 }
182 
183 #define boundfreq(freq)                                                                            \
184     {                                                                                              \
185         if (freq > 75)                                                                             \
186             freq = 75;                                                                             \
187         if (freq < -55)                                                                            \
188             freq = -55;                                                                            \
189     }
190 
Map2PoleResonance(double reso,double freq,int subtype)191 double Map2PoleResonance(double reso, double freq, int subtype)
192 {
193     switch (subtype)
194     {
195     case st_Medium:
196         reso *= max(0.0, 1.0 - max(0.0, (freq - 58) * 0.05));
197         return (0.99 - 1.0 * limit_range((double)(1 - (1 - reso) * (1 - reso)), 0.0, 1.0));
198     case st_Rough:
199         reso *= max(0.0, 1.0 - max(0.0, (freq - 58) * 0.05));
200         return (1.0 - 1.05 * limit_range((double)(1 - (1 - reso) * (1 - reso)), 0.001, 1.0));
201     default:
202     case st_Smooth:
203         return (2.5 - 2.45 * limit_range((double)(1 - (1 - reso) * (1 - reso)), 0.0, 1.0));
204     }
205     return 0;
206 }
207 
Map2PoleResonance_noboost(double reso,double freq,int subtype)208 double Map2PoleResonance_noboost(double reso, double freq, int subtype)
209 {
210     if (subtype == st_Rough)
211         return (1.0 - 0.99 * limit_range((double)(1 - (1 - reso) * (1 - reso)), 0.001, 1.0));
212     else
213         return (0.99 - 0.98 * limit_range((double)(1 - (1 - reso) * (1 - reso)), 0.0, 1.0));
214 }
215 
Map4PoleResonance(double reso,double freq,int subtype)216 double Map4PoleResonance(double reso, double freq, int subtype)
217 {
218     switch (subtype)
219     {
220     case st_Medium:
221         reso *= max(0.0, 1.0 - max(0.0, (freq - 58) * 0.05));
222         return 0.99 - 0.9949 * limit_range((double)reso, 0.0, 1.0);
223     case st_Rough:
224         reso *= max(0.0, 1.0 - max(0.0, (freq - 58) * 0.05));
225         return (1.0 - 1.05 * limit_range((double)reso, 0.001, 1.0));
226     default:
227     case st_Smooth:
228         return (2.5 - 2.3 * limit_range((double)reso, 0.0, 1.0));
229     }
230 }
231 
resoscale(double reso,int subtype)232 double resoscale(double reso, int subtype)
233 {
234     switch (subtype)
235     {
236     case st_Medium:
237         return (1.0 - 0.75 * reso * reso);
238     case st_Rough:
239         return (1.0 - 0.5 * reso * reso);
240     case st_Smooth:
241         return (1.0 - 0.25 * reso * reso);
242     }
243 
244     return 1.0;
245 }
246 
resoscale4Pole(double reso,int subtype)247 double resoscale4Pole(double reso, int subtype)
248 {
249     switch (subtype)
250     {
251     case st_Medium:
252         return (1.0 - 0.75 * reso);
253     case st_Rough:
254         return (1.0 - 0.5 * reso * reso);
255     case st_Smooth:
256         return (1.0 - 0.5 * reso);
257     }
258 
259     return 1.0;
260 }
261 
Coeff_SVF(float Freq,float Reso,bool FourPole)262 void FilterCoefficientMaker::Coeff_SVF(float Freq, float Reso, bool FourPole)
263 {
264     double f = 440.f * storage->note_to_pitch_ignoring_tuning(Freq);
265     double F1 = 2.0 * sin(M_PI * min(0.11, f * (0.25 * samplerate_inv))); // 4x oversampling
266 
267     Reso = sqrt(limit_range(Reso, 0.f, 1.f));
268 
269     double overshoot = FourPole ? 0.1 : 0.15;
270     double Q1 = 2.0 - Reso * (2.0 + overshoot) + F1 * F1 * overshoot * 0.9;
271     Q1 = min(Q1, min(2.00, 2.00 - 1.52 * F1));
272 
273     double ClipDamp = 0.1 * Reso * F1;
274 
275     const double a = 0.65;
276     double Gain = 1 - a * Reso;
277 
278     float c[n_cm_coeffs];
279     memset(c, 0, sizeof(float) * n_cm_coeffs);
280     c[0] = F1;
281     c[1] = Q1;
282     c[2] = ClipDamp;
283     c[3] = Gain;
284     FromDirect(c);
285 }
286 
Coeff_LP12(float freq,float reso,int subtype)287 void FilterCoefficientMaker::Coeff_LP12(float freq, float reso, int subtype)
288 {
289     float cosi, sinu;
290     float gain = resoscale(reso, subtype);
291 
292     boundfreq(freq) storage->note_to_omega_ignoring_tuning(freq, sinu, cosi);
293 
294     double alpha = sinu * Map2PoleResonance(reso, freq, subtype);
295 
296     if (subtype != st_Smooth)
297     {
298         alpha = min(alpha, sqrt(1.0 - cosi * cosi) - 0.0001);
299     }
300 
301     double a0 = 1 + alpha, a0inv = 1 / a0, a1 = -2 * cosi, a2 = 1 - alpha, b0 = (1 - cosi) * 0.5,
302            b1 = 1 - cosi, b2 = (1 - cosi) * 0.5;
303 
304     // double sq = a1*a1 - 4.0*a2;
305     // must be negative
306     // dvs a2 > 0.25*a1*a1
307     // 1-alpha / (1+alpha) > 0.25*(-2*cosi/(1+alpha))^2
308     // 1-alpha > 0.25*(-2*cosi)^2 / (1+alpha)
309     // (1-alpha)(1+alpha) > 0.25*(-2*cosi)^2
310     // (1-alpha)(1+alpha) > (cosi)^2
311     // >> alpha = +- sqrt(1 - cosi^2)
312 
313     // ar = 0.5 * -a1;
314     // ai = 0.5 * sqrt(-sq);
315     // ar = -(K*K - 1)/(1+alpha)
316     // ai =
317 
318     if (subtype == st_Smooth)
319         ToNormalizedLattice(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain,
320                             clipscale(freq, subtype));
321     else
322         ToCoupledForm(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain, clipscale(freq, subtype));
323 }
324 
Coeff_LP24(float freq,float reso,int subtype)325 void FilterCoefficientMaker::Coeff_LP24(float freq, float reso, int subtype)
326 {
327     float cosi, sinu;
328     float gain = resoscale(reso, subtype);
329 
330     boundfreq(freq) storage->note_to_omega_ignoring_tuning(freq, sinu, cosi);
331 
332     double Q2inv = Map4PoleResonance((double)reso, (double)freq, subtype);
333     double alpha = sinu * Q2inv;
334 
335     if (subtype != st_Smooth)
336     {
337         alpha = min(alpha, sqrt(1.0 - cosi * cosi) - 0.0001);
338     }
339 
340     double a0 = 1 + alpha, a0inv = 1 / a0, a1 = -2 * cosi, a2 = 1 - alpha, b0 = (1 - cosi) * 0.5,
341            b1 = 1 - cosi, b2 = (1 - cosi) * 0.5;
342 
343     if (subtype == st_Smooth)
344         ToNormalizedLattice(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain,
345                             clipscale(freq, subtype));
346     else
347         ToCoupledForm(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain, clipscale(freq, subtype));
348 }
349 
Coeff_HP12(float freq,float reso,int subtype)350 void FilterCoefficientMaker::Coeff_HP12(float freq, float reso, int subtype)
351 {
352     float cosi, sinu;
353     float gain = resoscale(reso, subtype);
354 
355     boundfreq(freq) storage->note_to_omega_ignoring_tuning(freq, sinu, cosi);
356 
357     double Q2inv = Map2PoleResonance(reso, freq, subtype);
358     double alpha = sinu * Q2inv;
359 
360     if (subtype != 0)
361     {
362         alpha = min(alpha, sqrt(1.0 - cosi * cosi) - 0.0001);
363     }
364 
365     double a0 = 1 + alpha, a0inv = 1 / a0, a1 = -2 * cosi, a2 = 1 - alpha, b0 = (1 + cosi) * 0.5,
366            b1 = -(1 + cosi), b2 = (1 + cosi) * 0.5;
367 
368     if (subtype == st_Smooth)
369         ToNormalizedLattice(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain,
370                             clipscale(freq, subtype));
371     else
372         ToCoupledForm(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain, clipscale(freq, subtype));
373 }
374 
Coeff_HP24(float freq,float reso,int subtype)375 void FilterCoefficientMaker::Coeff_HP24(float freq, float reso, int subtype)
376 {
377     float cosi, sinu;
378     float gain = resoscale(reso, subtype);
379 
380     boundfreq(freq) storage->note_to_omega_ignoring_tuning(freq, sinu, cosi);
381 
382     double Q2inv = Map4PoleResonance((double)reso, (double)freq, subtype);
383     double alpha = sinu * Q2inv;
384 
385     if (subtype != 0)
386     {
387         alpha = min(alpha, sqrt(1.0 - cosi * cosi) - 0.0001);
388     }
389 
390     double a0 = 1 + alpha, a0inv = 1 / a0, a1 = -2 * cosi, a2 = 1 - alpha, b0 = (1 + cosi) * 0.5,
391            b1 = -(1 + cosi), b2 = (1 + cosi) * 0.5;
392 
393     if (subtype == st_Smooth)
394         ToNormalizedLattice(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain,
395                             clipscale(freq, subtype));
396     else
397         ToCoupledForm(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain, clipscale(freq, subtype));
398 }
399 
Coeff_BP12(float freq,float reso,int subtype)400 void FilterCoefficientMaker::Coeff_BP12(float freq, float reso, int subtype)
401 {
402     float cosi, sinu;
403     float gain = resoscale(reso, subtype);
404 
405     if (subtype == st_Rough)
406     {
407         gain *= 2.f;
408     }
409 
410     boundfreq(freq) storage->note_to_omega_ignoring_tuning(freq, sinu, cosi);
411 
412     double Q2inv = Map2PoleResonance(reso, freq, subtype);
413     double Q = 0.5 / Q2inv;
414     double alpha = sinu * Q2inv;
415 
416     if (subtype != 0)
417     {
418         alpha = min(alpha, sqrt(1.0 - cosi * cosi) - 0.0001);
419     }
420 
421     double a0 = 1 + alpha, a0inv = 1 / a0, a1 = -2 * cosi, a2 = 1 - alpha, b0, b1, b2;
422 
423     {
424         b0 = Q * alpha;
425         b1 = 0;
426         b2 = -Q * alpha;
427     }
428 
429     if (subtype == st_Smooth)
430         ToNormalizedLattice(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain,
431                             clipscale(freq, subtype));
432     else
433         ToCoupledForm(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain, clipscale(freq, subtype));
434 }
435 
Coeff_BP24(float freq,float reso,int subtype)436 void FilterCoefficientMaker::Coeff_BP24(float freq, float reso, int subtype)
437 {
438     float cosi, sinu;
439     float gain = resoscale(reso, subtype);
440 
441     if (subtype == st_Rough)
442     {
443         gain *= 2.f;
444     }
445 
446     boundfreq(freq) storage->note_to_omega_ignoring_tuning(freq, sinu, cosi);
447 
448     double Q2inv = Map4PoleResonance(reso, freq, subtype);
449     double Q = 0.5 / Q2inv;
450     double alpha = sinu * Q2inv;
451 
452     if (subtype != 0)
453     {
454         alpha = min(alpha, sqrt(1.0 - cosi * cosi) - 0.0001);
455     }
456 
457     double a0 = 1 + alpha, a0inv = 1 / a0, a1 = -2 * cosi, a2 = 1 - alpha, b0, b1, b2;
458 
459     {
460         b0 = Q * alpha;
461         b1 = 0;
462         b2 = -Q * alpha;
463     }
464 
465     if (subtype == st_Smooth)
466         ToNormalizedLattice(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain,
467                             clipscale(freq, subtype));
468     else
469         ToCoupledForm(a0inv, a1, a2, b0 * gain, b1 * gain, b2 * gain, clipscale(freq, subtype));
470 }
471 
Coeff_Notch(float Freq,float Reso,int SubType)472 void FilterCoefficientMaker::Coeff_Notch(float Freq, float Reso, int SubType)
473 {
474     float cosi, sinu;
475     double Q2inv;
476 
477     boundfreq(Freq) storage->note_to_omega_ignoring_tuning(Freq, sinu, cosi);
478 
479     if (SubType == st_NotchMild)
480     {
481         Q2inv = (1.00 - 0.99 * limit_range((double)(1 - (1 - Reso) * (1 - Reso)), 0.0, 1.0));
482     }
483     else
484     {
485         Q2inv = (2.5 - 2.49 * limit_range((double)(1 - (1 - Reso) * (1 - Reso)), 0.0, 1.0));
486     }
487 
488     double alpha = sinu * Q2inv;
489 
490     double a0 = 1 + alpha, a0inv = 1 / a0, a1 = -2 * cosi, a2 = 1 - alpha, b0 = 1, b1 = -2 * cosi,
491            b2 = 1;
492 
493     ToNormalizedLattice(a0inv, a1, a2, b0, b1, b2, 0.005);
494 }
495 
Coeff_APF(float Freq,float Reso,int SubType)496 void FilterCoefficientMaker::Coeff_APF(float Freq, float Reso, int SubType)
497 {
498     float cosi, sinu;
499     double Q2inv;
500 
501     boundfreq(Freq) storage->note_to_omega_ignoring_tuning(Freq, sinu, cosi);
502 
503     Q2inv = (2.5 - 2.49 * limit_range((double)(1 - (1 - Reso) * (1 - Reso)), 0.0, 1.0));
504 
505     double alpha = sinu * Q2inv;
506 
507     double a0 = 1 + alpha, a0inv = 1 / a0, a1 = -2 * cosi, a2 = 1 - alpha, b0 = 1 - alpha,
508            b1 = -2 * cosi, b2 = 1 + alpha;
509 
510     ToNormalizedLattice(a0inv, a1, a2, b0, b1, b2, 0.005);
511 }
512 
Coeff_LP4L(float freq,float reso,int subtype)513 void FilterCoefficientMaker::Coeff_LP4L(float freq, float reso, int subtype)
514 {
515     double gg = limit_range(
516         ((double)440 * storage->note_to_pitch_ignoring_tuning(freq) * dsamplerate_os_inv), 0.0,
517         0.187); // gg
518 
519     float t_b1 = 1.f - exp(-2 * M_PI * gg);
520     float q = min(2.15f * limit_range(reso, 0.f, 1.f), 0.5f / (t_b1 * t_b1 * t_b1 * t_b1));
521 
522     float c[n_cm_coeffs];
523     memset(c, 0, sizeof(float) * n_cm_coeffs);
524     c[0] = (3.f / (3.f - q));
525     c[1] = t_b1;
526     c[2] = q;
527 
528     FromDirect(c);
529 }
530 
Coeff_COMB(float freq,float reso,int isubtype)531 void FilterCoefficientMaker::Coeff_COMB(float freq, float reso, int isubtype)
532 {
533     int subtype = isubtype & QFUSubtypeMasks::UNMASK_SUBTYPE;
534     bool extended = isubtype & QFUSubtypeMasks::EXTENDED_COMB;
535 
536     int comb_length = extended ? MAX_FB_COMB_EXTENDED : MAX_FB_COMB;
537 
538     float dtime = (1.f / 440.f) * storage->note_to_pitch_inv_ignoring_tuning(freq);
539     dtime = dtime * dsamplerate_os;
540 
541     // See comment in SurgeStorage and issue #3248
542     if (!storage->_patch->correctlyTuneCombFilter)
543     {
544         dtime -= FIRoffset;
545     }
546 
547     dtime = limit_range(dtime, (float)FIRipol_N, (float)comb_length - FIRipol_N);
548     if (extended)
549     {
550         // extended use is not from the filter bank so allow greater feedback range
551         reso = limit_range(reso, -2.f, 2.f);
552     }
553     else
554     {
555         reso = ((subtype & 2) ? -1.0f : 1.0f) * limit_range(reso, 0.f, 1.f);
556     }
557 
558     float c[n_cm_coeffs];
559     memset(c, 0, sizeof(float) * n_cm_coeffs);
560     c[0] = dtime;
561     c[1] = reso;
562     c[2] = (subtype & 1) ? 0.0f : 0.5f; // combmix
563     c[3] = 1.f - c[2];
564     FromDirect(c);
565 }
566 
Coeff_SNH(float freq,float reso,int subtype)567 void FilterCoefficientMaker::Coeff_SNH(float freq, float reso, int subtype)
568 {
569     float dtime = (1.f / 440.f) * storage->note_to_pitch_ignoring_tuning(-freq) * dsamplerate_os;
570     double v1 = 1.0 / dtime;
571 
572     float c[n_cm_coeffs];
573     memset(c, 0, sizeof(float) * n_cm_coeffs);
574     c[0] = v1;
575     c[1] = reso;
576     FromDirect(c);
577 }
578 
FromDirect(float N[n_cm_coeffs])579 void FilterCoefficientMaker::FromDirect(float N[n_cm_coeffs])
580 {
581     if (FirstRun)
582     {
583         memset(dC, 0, sizeof(float) * n_cm_coeffs);
584         memcpy(C, N, sizeof(float) * n_cm_coeffs);
585         memcpy(tC, N, sizeof(float) * n_cm_coeffs);
586 
587         FirstRun = false;
588     }
589     else
590     {
591         for (int i = 0; i < n_cm_coeffs; i++)
592         {
593             tC[i] = (1.f - smooth) * tC[i] + smooth * N[i];
594             dC[i] = (tC[i] - C[i]) * BLOCK_SIZE_OS_INV;
595         }
596     }
597 }
598 
ToNormalizedLattice(double a0inv,double a1,double a2,double b0,double b1,double b2,double g)599 void FilterCoefficientMaker::ToNormalizedLattice(double a0inv, double a1, double a2, double b0,
600                                                  double b1, double b2, double g)
601 {
602     b0 *= a0inv;
603     b1 *= a0inv;
604     b2 *= a0inv;
605     a1 *= a0inv;
606     a2 *= a0inv;
607 
608     float N[n_cm_coeffs];
609     memset(N, 0, sizeof(float) * n_cm_coeffs);
610 
611     double k1 = a1 / (1.0 + a2);
612     double k2 = a2;
613 
614     double q1 = 1.0 - k1 * k1;
615     double q2 = 1.0 - k2 * k2;
616     q1 = sqrt(fabs(q1));
617     q2 = sqrt(fabs(q2));
618 
619     double v3 = b2;
620     double v2 = (b1 - a1 * v3) / q2;
621     double v1 = (b0 - k1 * v2 * q2 - k2 * v3) / (q1 * q2);
622 
623     N[0] = k1;
624     N[1] = k2;
625     N[2] = q1;
626     N[3] = q2;
627     N[4] = v1;
628     N[5] = v2;
629     N[6] = v3;
630     N[7] = g;
631 
632     FromDirect(N);
633 }
ToCoupledForm(double a0inv,double a1,double a2,double b0,double b1,double b2,double g)634 void FilterCoefficientMaker::ToCoupledForm(double a0inv, double a1, double a2, double b0, double b1,
635                                            double b2, double g)
636 {
637     b0 *= a0inv;
638     b1 *= a0inv;
639     b2 *= a0inv;
640     a1 *= a0inv;
641     a2 *= a0inv;
642 
643     float N[n_cm_coeffs];
644     memset(N, 0, sizeof(float) * n_cm_coeffs);
645 
646     double ar, ai;
647     double sq = a1 * a1 - 4.0 * a2;
648 
649     ar = 0.5 * -a1;
650     sq = min(0.0, sq);
651     ai = 0.5 * sqrt(-sq);
652     ai = max(ai, 8.0 * 1.192092896e-07F);
653 
654     double bb1 = b1 - a1 * b0;
655     double bb2 = b2 - a2 * b0;
656 
657     double d = b0;
658     double c1 = bb1;
659     double c2 = (bb1 * ar + bb2) / ai;
660 
661     double scaler = 1.f; // 0.01 + 0.99*sqrt(c1*c1);
662 
663     N[0] = ar;
664     N[1] = ai;
665     N[2] = 1.f; // scaler;
666     N[4] = c1;  // /scaler;
667     N[5] = c2;  // /scaler;
668     N[6] = d;
669     N[7] = g;
670 
671     FromDirect(N);
672 }
673 
Reset()674 void FilterCoefficientMaker::Reset()
675 {
676     FirstRun = true;
677     memset(C, 0, sizeof(float) * n_cm_coeffs);
678     memset(dC, 0, sizeof(float) * n_cm_coeffs);
679     memset(tC, 0, sizeof(float) * n_cm_coeffs);
680 
681     storage = nullptr;
682 }
683