Changing the thread title to reflect the new approach. Extract of the original thread;
- I suggested using iemlib's hi pass filter recipe to improve frequency response of [hip~] - Christof Ressi pointed to formula in http://www.arpchord.com/pdf/coeffs_first_order_filters_0p1.pdf - this formula calculates feedback coefficient k = (1 - sin(a)) / cos(a) where a = 2 * pi * fc / SR - the filter implementation is y[n] = (x[n] - x[n-1]) * (1 + k) / 2 + k * y[n-1] - following convention in d_filter.c (and pd tilde classes in general), trig functions should best be approximated - Cyrille provided libre office linear regression result for (1-sin(x))/cos(x) Thanks for the useful infos and discussion. My 'math coach' suggested using odd powers of -(x-pi/2) in an approximation polynomial for (1-sin(x))/cos(x). The best accuracy/performance balance I could get is with this 5th degree polynomial: (-(x-pi/2))*0.4908 - (x-pi/2)^3*0.04575 - (x-pi/2)^5*0.00541 Using this approximation in the filter formula, response at cutoff frequency is -3 dB with +/-0.06 dB accuracy in the required range 0 < x < pi. It can be efficiently implemented in C, analogous to an approximation Miller uses in [bp~]. So that is what I'll try next. Attached patch hip~-models.pd illustrates and compares filter recipes using vanilla objects: - current implementation, most efficient, accuracy +/- 3 dB - implementation with trig functions, least efficient, accuracy +/- 0.01 dB - implementation with approximation for trig functions, efficient, accuracy +/- 0.06 dB A note on efficiency: coefficients in [hip~] are only recalculated when cutoff frequency is changed. How important is performance for a function rarely called? I'm much aware of the motto 'never optimize early', yet I spent much time on finding a fast approximation, for several reasons: it's a nice math challenge, instructive for cases where performance matters more, and I want to respect Miller's code efficiency when proposing a change. Today pd is even deployed on embedded devices so the frugal coding approach is still relevant. After 20 years. Katja On Tue, Oct 18, 2016 at 10:28 AM, cyrille henry <c...@chnry.net> wrote: > > > Le 18/10/2016 à 00:47, katja a écrit : >> >> The filter recipe that Christof pointed to was easy to plug into the C >> code of [hip~] and works perfectly. But when looking further in >> d_filter.c I came across an approximation function 'sigbp_qcos()' used >> in the bandpass filter. It made me realize once more how passionate >> Miller is about efficiency. I'm not going to make a fool of myself by >> submitting a 'fix' using two trig functions to calculate a filter >> coefficient when a simple approximation could do the job. So that is >> what I'm now looking into, with help of a math friend: an efficient >> polynomial approximation for (1-sin(x))/cos(x). > > according to libre office linear regression, for x between 0 and Pi, > (1-sin(x))/cos(x) is about : > -0.057255x³ + 0.27018x² - 0.9157x + 0.99344 > > the calc is in attachment, if you want to tune the input source or > precision. > cheers > c
#N canvas 56 116 600 483 10; #X obj 142 161 osc~; #X msg 501 129 \; pd dsp 1; #X msg 501 176 \; pd dsp 0; #X floatatom 142 98 8 0 0 0 - - -, f 8; #X obj 142 260 env~ 8192; #X floatatom 142 297 5 0 0 0 - - -, f 5; #X text 171 161 test signal; #X text 54 447 Katja Vetter Oct. 2016; #X text 194 97 Hz; #X obj 227 260 env~ 8192; #X floatatom 227 301 5 0 0 0 - - -, f 5; #X obj 227 331 -; #X floatatom 227 365 5 0 0 0 - - -, f 5; #X obj 145 45 hsl 200 15 0 22050 0 0 empty empty frequency(coarse) 20 8 0 10 -262144 -260097 -1 19900 0; #X text 19 58 edge cases; #X obj 57 257 env~ 8192; #X floatatom 57 298 5 0 0 0 - - -, f 5; #X obj 57 328 -; #X floatatom 57 362 5 0 0 0 - - -, f 5; #X obj 359 260 env~ 8192; #X floatatom 359 301 5 0 0 0 - - -, f 5; #X obj 359 331 -; #X floatatom 359 365 5 0 0 0 - - -, f 5; #X obj 145 16 hsl 200 15 0 1000 0 0 empty empty frequency(fine) 20 8 0 10 -262144 -260097 -1 9000 0; #N canvas 110 156 413 534 hip1~ 0; #X obj 40 25 inlet~; #X obj 40 436 outlet~; #X floatatom 122 414 8 0 0 0 - - -, f 8; #X text 177 413 normalisation; #X obj 40 405 *~ 0; #X obj 40 345 rpole~; #X obj 90 56 cnv 15 300 300 empty empty empty 20 12 0 14 -233017 -66577 0; #X obj 306 126 atan; #X msg 306 99 1; #X obj 122 24 inlet; #X floatatom 306 181 8 0 0 0 - - -, f 8; #X obj 137 88 samplerate~; #X obj 122 118 /; #X obj 306 71 loadbang; #X obj 122 169 *; #X floatatom 155 119 8 0 0 0 - - -, f 8; #X floatatom 152 186 8 0 0 0 - - -, f 8; #X text 361 181 pi; #X text 161 23 cut off frequency f(0); #X obj 40 378 rzero~ 1; #X obj 306 152 * 8; #X obj 122 388 / 2; #X text 204 185 a = 2*pi*fc/SR; #X obj 122 361 + 1; #X floatatom 154 330 8 0 0 0 - - -, f 8; #X obj 122 258 t b f; #X msg 122 283 1; #X obj 122 313 -; #X text 210 328 1-a; #X text 36 474 Efficient [hip~] as implemented in vanilla pd (d_filter.c). ; #X obj 122 231 clip 0 1; #X connect 0 0 5 0; #X connect 4 0 1 0; #X connect 5 0 19 0; #X connect 7 0 20 0; #X connect 8 0 7 0; #X connect 9 0 12 0; #X connect 11 0 12 1; #X connect 11 0 15 0; #X connect 12 0 14 0; #X connect 13 0 11 0; #X connect 13 0 8 0; #X connect 14 0 16 0; #X connect 14 0 30 0; #X connect 19 0 4 0; #X connect 20 0 10 0; #X connect 20 0 14 1; #X connect 21 0 2 0; #X connect 21 0 4 1; #X connect 23 0 21 0; #X connect 25 0 26 0; #X connect 25 1 27 1; #X connect 26 0 27 0; #X connect 27 0 23 0; #X connect 27 0 24 0; #X connect 27 0 5 1; #X connect 30 0 25 0; #X restore 57 208 pd hip1~; #N canvas 547 139 433 582 hip2~ 0; #X obj 40 25 inlet~; #X obj 40 453 outlet~; #X floatatom 121 431 8 0 0 0 - - -, f 8; #X text 176 430 normalisation; #X obj 40 422 *~ 0; #X obj 40 362 rpole~; #X obj 89 52 cnv 15 300 300 empty empty empty 20 12 0 14 -233017 -66577 0; #X obj 121 287 /; #X floatatom 165 304 8 0 0 0 - - -, f 8; #X obj 306 113 atan; #X msg 306 89 1; #X obj 121 24 inlet; #X floatatom 306 168 8 0 0 0 - - -, f 8; #X obj 136 78 samplerate~; #X obj 121 105 /; #X obj 306 61 loadbang; #X obj 121 156 *; #X floatatom 154 106 8 0 0 0 - - -, f 8; #X floatatom 151 173 8 0 0 0 - - -, f 8; #X text 160 23 cut off frequency f(0); #X obj 40 395 rzero~ 1; #X obj 306 139 * 8; #X obj 186 230 cos; #X obj 152 230 sin; #X obj 121 199 t b f f; #X obj 121 262 -; #X msg 121 230 1; #X obj 121 405 / 2; #X text 220 302 k = (1-sin(a))/cos(a); #X text 203 172 a = 2*pi*fc/SR; #X obj 121 378 + 1; #X text 35 528 http://www.arpchord.com/pdf/coeffs_first_order_filters_0p1.pdf ; #X text 36 478 filter recipe adapted from Cliff Sparks \, Coefficients of Recursive Linear Time-Invariant First-Order Low-Pass and High-Pass Filters; #X floatatom 164 279 8 0 0 0 - - -, f 8; #X text 216 279 1-sin(a); #X text 357 169 2*pi; #X obj 121 324 clip -1 1; #X floatatom 135 357 8 0 0 0 - - -, f 8; #X connect 0 0 5 0; #X connect 4 0 1 0; #X connect 5 0 20 0; #X connect 7 0 8 0; #X connect 7 0 36 0; #X connect 9 0 21 0; #X connect 10 0 9 0; #X connect 11 0 14 0; #X connect 13 0 14 1; #X connect 13 0 17 0; #X connect 14 0 16 0; #X connect 15 0 13 0; #X connect 15 0 10 0; #X connect 16 0 18 0; #X connect 16 0 24 0; #X connect 20 0 4 0; #X connect 21 0 12 0; #X connect 21 0 16 1; #X connect 22 0 7 1; #X connect 23 0 25 1; #X connect 24 0 26 0; #X connect 24 1 23 0; #X connect 24 2 22 0; #X connect 25 0 7 0; #X connect 25 0 33 0; #X connect 26 0 25 0; #X connect 27 0 2 0; #X connect 27 0 4 1; #X connect 30 0 27 0; #X connect 36 0 30 0; #X connect 36 0 5 1; #X connect 36 0 37 0; #X restore 228 211 pd hip2~; #N canvas 558 130 433 538 hip3~ 0; #X obj 40 25 inlet~; #X obj 40 436 outlet~; #X floatatom 122 414 8 0 0 0 - - -, f 8; #X text 177 413 normalisation; #X obj 40 405 *~ 0; #X obj 40 345 rpole~; #X obj 90 56 cnv 15 300 300 empty empty empty 20 12 0 14 -233017 -66577 0; #X obj 306 126 atan; #X msg 306 99 1; #X obj 122 24 inlet; #X floatatom 306 181 8 0 0 0 - - -, f 8; #X obj 137 88 samplerate~; #X obj 122 118 /; #X obj 306 71 loadbang; #X obj 122 169 *; #X floatatom 155 119 8 0 0 0 - - -, f 8; #X floatatom 152 186 8 0 0 0 - - -, f 8; #X text 361 181 pi; #X text 161 23 cut off frequency f(0); #X obj 40 378 rzero~ 1; #X obj 306 152 * 8; #X obj 122 388 / 2; #X text 204 185 a = 2*pi*fc/SR; #X obj 122 361 + 1; #X floatatom 142 262 10 0 0 0 - - -, f 10; #N canvas 235 281 517 393 approx(1-sin(a))/cos(a) 0; #X obj 47 12 inlet; #X obj 47 266 outlet; #X floatatom 61 114 10 0 0 0 - - -, f 10; #X obj 47 232 clip -1 1; #X obj 245 29 loadbang; #X obj 47 141 expr ($f1*$f2)+(pow($f1 \, 3)*$f3)+(pow($f1 \, 5)*$f4) ; #X text 138 179 (x-pi/2)*a + (x-pi/2)^3*b + (x-pi/2)^5*c; #X text 196 93 a; #X text 302 93 b; #X text 402 93 c; #X floatatom 61 179 10 0 0 0 - - -, f 10; #X obj 23 45 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144 -1 -1; #X msg 146 94 -0.4908; #X msg 245 94 -0.04575; #X msg 344 95 -0.00541; #N canvas 1 90 450 300 pi/2 0; #X obj 48 86 atan; #X msg 48 59 1; #X obj 48 31 loadbang; #X obj 48 148 outlet; #X floatatom 105 129 8 0 0 0 - - -, f 8; #X obj 48 112 * 2; #X text 163 128 pi/2; #X connect 0 0 5 0; #X connect 1 0 0 0; #X connect 2 0 1 0; #X connect 5 0 3 0; #X connect 5 0 4 0; #X restore 62 46 pd pi/2; #X obj 47 75 -; #X floatatom 73 77 8 0 0 0 - - -, f 8; #X text 140 205 Approximation for (1-sin(x))/cos(x). General formula proposed by Joe Deken. Factors a b c tuned 'by hand' \, Katja Vetter Oct. 2016; #X text 45 302 Results would be slightly different in double precision \, being not exactly 1 for fc=0 and -1 for fc=SR/2. To be sure \, result -0.999999 should probably be forced to -1 and 0.999999 to 1; #X connect 0 0 16 0; #X connect 3 0 1 0; #X connect 4 0 12 0; #X connect 4 0 13 0; #X connect 4 0 14 0; #X connect 5 0 3 0; #X connect 5 0 10 0; #X connect 11 0 16 0; #X connect 12 0 5 1; #X connect 13 0 5 2; #X connect 14 0 5 3; #X connect 15 0 16 1; #X connect 15 0 17 0; #X connect 16 0 5 0; #X connect 16 0 2 0; #X restore 122 228 pd approx(1-sin(a))/cos(a); #X obj 122 297 clip -1 1; #X floatatom 140 331 8 0 0 0 - - -, f 8; #X text 195 333 feedback coefficient; #X text 39 466 Hipass filter recipe using an approximation of (1-sin(a))/cos(a) for efficient calculation of the coefficient.; #X connect 0 0 5 0; #X connect 4 0 1 0; #X connect 5 0 19 0; #X connect 7 0 20 0; #X connect 8 0 7 0; #X connect 9 0 12 0; #X connect 11 0 12 1; #X connect 11 0 15 0; #X connect 12 0 14 0; #X connect 13 0 11 0; #X connect 13 0 8 0; #X connect 14 0 16 0; #X connect 14 0 25 0; #X connect 19 0 4 0; #X connect 20 0 10 0; #X connect 20 0 14 1; #X connect 21 0 2 0; #X connect 21 0 4 1; #X connect 23 0 21 0; #X connect 25 0 24 0; #X connect 25 0 26 0; #X connect 26 0 23 0; #X connect 26 0 5 1; #X connect 26 0 27 0; #X restore 359 211 pd hip3~; #X text 62 231 cheapest; #X text 233 234 trig functions; #X text 367 233 approximation; #X text 92 362 dB gain; #X text 261 365 dB gain; #X text 394 364 dB gain; #X obj 235 77 hsl 200 15 0 22050 0 0 empty empty frequency(coarse) 20 8 0 10 -262144 -260097 -1 2300 1; #X floatatom 232 101 8 0 0 0 - - -, f 8; #X msg 22 81 -10; #X msg 68 81 600000; #X text 365 31 << cutoff and test signal; #X text 447 77 << test signal only; #X text 53 404 Test attenuation at cutoff frequency for different hipass filter recipes. Output should be -3dB at cutoff.; #X connect 0 0 4 0; #X connect 0 0 25 0; #X connect 0 0 24 0; #X connect 0 0 26 0; #X connect 3 0 25 1; #X connect 3 0 0 0; #X connect 3 0 24 1; #X connect 3 0 26 1; #X connect 4 0 5 0; #X connect 5 0 11 1; #X connect 5 0 17 1; #X connect 5 0 21 1; #X connect 9 0 10 0; #X connect 10 0 11 0; #X connect 11 0 12 0; #X connect 13 0 3 0; #X connect 15 0 16 0; #X connect 16 0 17 0; #X connect 17 0 18 0; #X connect 19 0 20 0; #X connect 20 0 21 0; #X connect 21 0 22 0; #X connect 23 0 13 0; #X connect 24 0 15 0; #X connect 25 0 9 0; #X connect 26 0 19 0; #X connect 33 0 0 0; #X connect 33 0 34 0; #X connect 35 0 3 0; #X connect 36 0 3 0;
_______________________________________________ Pd-list@lists.iem.at mailing list UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
