[PD] different float accuracy betw. Pd-0.52-1-msw-i386 and Pd-0.52-1-msw-amd64

Thu Feb 10 22:45:55 CET 2022

Generally, I can reproduce the output, including the funny noises on 
Pd-0.52.1 64-bit.

> IOhannes z. and Chrostof R. figured out, it depends on compiler 
> options -ffast-math and / or -fassociative-math.
That's only half the story. Generally, one shouldn't expect floating 
point computations to yield the exact same result with different 
compilers/machines. The rounding errors themselves are very small, but 
they can accumulate over long periods of time or get amplified by 
recursive algorithms.

However, in this case, the actual problem is really the implementation 
of [biquad~]. It uses direct form 2 - which requires fewer delay units, 
but is less numerically stable. 
(https://en.wikipedia.org/wiki/Digital_biquad_filter)

I've made an alternative implementation [biquad2~] that uses direct form 
1. It shows significantly less noise than the [biquad~] examples. See 
attachments.

Christof

On 10.02.2022 14:01, musil at iem.at wrote:
> ISSUE: floating-point-inaccuracy of low frequency filters with low 
> frequency signals in pd vanilla
>
> The patch "Test_HP3_butterworth_ at _20_Hz.pd" is part of the Live 
> Electronic of K.H.Stockhausens piece Mikrophonie I.
> My colleague David P. told me, that there occur hearable noise, 
> distortion, oscillations and rustling during perfoming this piece.
>
> The test-patch is a 0.1 Hz oscillator sourcing a highpass filter 
> 3.order with butterworth characteristic at 20 Hz.
>
> Try this patch with pd-0.52-1-msw-i386 and amd64 aka 32-bit and 64-bit 
> and older pd versions.
> The 32 bit version has 20 dB less noise than the 64 bit version (and 
> no oscillations).
>
> Double precision filters of iemlib ("hp3_butt_dp~.pd") have less noise 
> than single precision filters.
>
> IOhannes z. and Chrostof R. figured out, it depends on compiler 
> options -ffast-math and / or -fassociative-math.
>
> Some screenshots are added.
>
> have fun
>
> Thomas Musil
>
> _______________________________________________
> Pd-list at lists.iem.at mailing list
> UNSUBSCRIBE and account-management -> https://lists.puredata.info/listinfo/pd-list
-------------- next part --------------
A non-text attachment was scrubbed...
Name: biquad2~.zip
Type: application/zip
Size: 52051 bytes
Desc: not available
URL: <http://lists.puredata.info/pipermail/pd-list/attachments/20220210/eb35fe5a/attachment-0001.zip>
-------------- next part --------------
#N canvas 192 20 1345 798 12;
#N canvas 25 25 801 679 hp_3.order_butterworth_ at _20_Hz 0;
#X obj 83 52 inlet~;
#X obj 83 597 outlet~;
#X text 76 617 signal out;
#X text 74 33 signal in;
#X obj 336 51 inlet;
#X obj 525 38 loadbang;
#N canvas 0 50 654 536 freq_Q_transform 0;
#X obj 77 60 inlet;
#X text 79 42 freq;
#X obj 77 357 outlet;
#X obj 104 233 expr 4*atan(1)/$f1;
#X obj 104 209 samplerate~;
#X obj 77 252 * 1;
#X obj 77 277 clip 1e-20 1.57079;
#X obj 240 211 expr 2*atan(1);
#X obj 77 302 expr cos($f1)/sin($f1);
#X obj 240 233 - 1e-07;
#X text 147 336 l a;
#X text 251 254 pi/2;
#X text 120 252 *pi/SR;
#X text 12 3 fq-al : freq Q to l a transformer;
#X obj 510 57 inlet;
#X text 511 38 init;
#X obj 510 82 t b b;
#X obj 77 334 pack 0 1;
#X connect 0 0 5 0;
#X connect 3 0 5 1;
#X connect 4 0 3 0;
#X connect 5 0 6 0;
#X connect 6 0 8 0;
#X connect 7 0 9 0;
#X connect 8 0 17 0;
#X connect 9 0 6 2;
#X connect 14 0 16 0;
#X connect 16 0 4 0;
#X connect 16 1 7 0;
#X connect 17 0 2 0;
#X restore 124 168 pd freq_Q_transform;
#X text 267 168 l a;
#X obj 597 38 r pd-dsp-started;
#X obj 124 196 expr $f1*$f1+1 \; $f1*$f2;
#X text 103 229 l2;
#X text 234 232 al;
#X obj 136 241 expr $f1 \; $f2 \; 1/($f1+$f2);
#X text 99 283 l2;
#X text 185 283 al;
#X text 254 283 rcp;
#X obj 113 448 pack 0 0 0 0 0;
#X text 513 276 rcp;
#X obj 406 420 pack 0 0 0 0 0;
#X obj 427 249 expr $f1 \; 1/($f1+1);
#X text 397 276 l;
#N canvas 259 301 581 451 freq_transform 0;
#X obj 104 283 expr 4*atan(1)/$f1;
#X obj 104 209 samplerate~;
#X obj 77 302 * 1;
#X obj 77 60 inlet;
#X text 79 42 freq;
#X obj 77 327 clip 1e-20 1.57079;
#X obj 240 211 expr 2*atan(1);
#X obj 77 352 expr cos($f1)/sin($f1);
#X obj 240 283 - 1e-07;
#X text 131 389 l;
#X text 251 304 pi/2;
#X text 120 302 *pi/SR;
#X obj 77 388 outlet;
#X obj 164 131 t b b;
#X floatatom 125 242 5 0 0 0 - - - 0;
#X obj 164 62 inlet;
#X text 163 43 init;
#X connect 0 0 2 1;
#X connect 1 0 0 0;
#X connect 1 0 14 0;
#X connect 2 0 5 0;
#X connect 3 0 2 0;
#X connect 5 0 7 0;
#X connect 6 0 8 0;
#X connect 7 0 12 0;
#X connect 8 0 5 2;
#X connect 13 0 1 0;
#X connect 13 1 6 0;
#X connect 15 0 13 0;
#X restore 427 213 pd freq_transform;
#X text 553 211 l;
#X obj 525 68 t b b;
#X obj 113 299 expr 2*$f3*($f1-2) \; $f3*($f2-$f1) \; $f3*($f1-1) \;
-2*$f3*($f1-1) \; $f3*($f1-1);
#X obj 406 297 expr ($f1-1)*$f2 \; $f1*$f2 \; -$f1*$f2;
#X msg 720 38 bang;
#N canvas 169 171 504 523 stability 0;
#X obj 87 43 inlet;
#X obj 89 339 outlet;
#X msg 124 148 -1 2e-07;
#X obj 124 173 +;
#X obj 146 234 +;
#X msg 146 209 1 -2e-07;
#X obj 87 265 max -1;
#X obj 87 291 min 1, f 9;
#X obj 124 123 t b b, f 6;
#X obj 126 94 loadbang;
#X text 92 15 fb1;
#X text 92 370 fb1;
#X text 178 255 if(fb1 <= -0.999 9998) \; ..fb1 = -0.999 9998 \; else
if(fb1 >= 0.999 9998) \; ..fb1 = 0.999 9998 \;;
#X connect 0 0 6 0;
#X connect 2 0 3 0;
#X connect 3 0 6 1;
#X connect 4 0 7 1;
#X connect 5 0 4 0;
#X connect 6 0 7 0;
#X connect 7 0 1 0;
#X connect 8 0 2 0;
#X connect 8 1 5 0;
#X connect 9 0 8 0;
#X restore 383 357 pd stability check;
#X f 9;
#N canvas 169 171 1356 1002 stability 0;
#X obj 87 53 inlet;
#X obj 87 859 outlet;
#X obj 124 233 +;
#X obj 163 259 +;
#X obj 124 183 t b b, f 6;
#X obj 124 144 loadbang, f 4;
#X text 91 28 fb1;
#X text 90 887 fb1;
#X obj 407 53 inlet;
#X obj 307 859 outlet;
#X msg 444 208 -1 2e-07;
#X obj 444 233 +;
#X obj 483 294 +;
#X msg 483 269 1 -2e-07;
#X obj 407 325 max -1;
#X obj 407 351 min 1, f 11;
#X obj 444 183 t b b, f 6;
#X obj 444 143 loadbang, f 4;
#X text 411 28 fb2;
#X text 310 887 fb2;
#X msg 163 234 2 -4e-07;
#X msg 124 208 -2 4e-07;
#X obj 87 325 max -2;
#X obj 87 351 min 2, f 11;
#X obj 407 80 t f f, f 66;
#X obj 87 113 t f f, f 92;
#X obj 728 230 expr $f1*$f1 + 4*$f2;
#X obj 728 259 >= 0;
#X obj 87 445 route 0 1, f 63;
#X obj 87 391 pack 0 0 0, f 92;
#X obj 87 823 unpack 0 0, f 32;
#X obj 306 512 unpack 0 0, f 15;
#X text 990 404 if(x->filter_function_is_first_order) \; if(fb1 <=
-0.999 9998) \; ..fb1 = -0.999 9998 \; else if(fb1 >= 0.999 9998) \;
..fb1 = 0.999 9998 \; \; else \; \; discriminant = x->b1 * x->b1 +
4 * x->b2 \; \; if(fb1 <= -1.999 9996) \; ..fb1 = -1.999 9996 \; else
if(fb1 >= 1.999 9996) \; ..fb1 = 1.999 9996 \; \; if(fb2 <= -0.999
9998) \; ..fb2 = -0.999 9998 \; else if(fb2 >= 0.999 9998) \; ..fb2
= 0.999 9998 \; \; if(discriminant >= 0.0) \; ..if(0.999 9998 - fb1
- fb2 < 0.0) \; ....fb2 = 0.999 9998 - fb1 \; ..if(0.999 9998 + fb1
- fb2 < 0.0) \; ....fb2 = 0.999 9998 + fb1 \; \; once again \; \; if(discriminant
>= 0.0) \; ..if(fb2 >= 0.999 9998 - fb1) \; ....fb2 = 0.999 9998 -
fb1 \; ..if(fb2 >= 0.999 9998 + fb1) \; ....fb2 = 0.999 9998 + fb1
\;;
#X text 187 139 if(fb1 <= -1.999 9996) \; ..fb1 = -1.999 9996 \; else
if(fb1 >= 1.999 9996) \; ..fb1 = 1.999 9996 \;;
#X text 877 231 discriminant = x->b1 * x->b1 + 4 * x->b2 \; if(discriminant
>= 0.0) \; .... \;;
#X text 493 140 if(fb2 <= -0.999 9998) \; ..fb2 = -0.999 9998 \; else
if(fb2 >= 0.999 9998) \; ..fb2 = 0.999 9998 \;;
#X text 505 557 if(discriminant >= 0.0) \; .. \; ..if(fb2 >= 0.999
9998 - fb1) \; ....fb2 = 0.999 9998 - fb1 \; ..if(fb2 >= 0.999 9998
+ fb1) \; ....fb2 = 0.999 9998 + fb1 \;;
#X msg 87 418 \$3 \$1 \$2;
#X msg 306 482 \$2 \$1;
#X obj 454 605 +;
#X obj 431 575 * -1;
#X obj 306 631 min, f 21;
#X obj 443 663 +, f 8;
#X obj 306 690 min, f 20;
#X obj 306 717 pack 0 0;
#X obj 408 541 t f f f;
#X msg 306 745 \$2 \$1;
#X connect 0 0 25 0;
#X connect 2 0 22 1;
#X connect 3 0 23 1;
#X connect 4 0 21 0;
#X connect 4 1 20 0;
#X connect 5 0 4 0;
#X connect 8 0 24 0;
#X connect 10 0 11 0;
#X connect 11 0 14 1;
#X connect 12 0 15 1;
#X connect 12 0 39 1;
#X connect 12 0 42 1;
#X connect 13 0 12 0;
#X connect 14 0 15 0;
#X connect 15 0 29 1;
#X connect 16 0 10 0;
#X connect 16 1 13 0;
#X connect 17 0 16 0;
#X connect 20 0 3 0;
#X connect 21 0 2 0;
#X connect 22 0 23 0;
#X connect 23 0 29 0;
#X connect 24 0 14 0;
#X connect 24 1 26 1;
#X connect 25 0 22 0;
#X connect 25 1 26 0;
#X connect 26 0 27 0;
#X connect 27 0 29 2;
#X connect 28 0 30 0;
#X connect 28 1 38 0;
#X connect 29 0 37 0;
#X connect 30 0 1 0;
#X connect 30 1 9 0;
#X connect 31 0 41 0;
#X connect 31 1 45 0;
#X connect 37 0 28 0;
#X connect 38 0 31 0;
#X connect 39 0 41 1;
#X connect 40 0 42 0;
#X connect 41 0 43 0;
#X connect 42 0 43 1;
#X connect 43 0 44 0;
#X connect 44 0 46 0;
#X connect 45 0 44 1;
#X connect 45 1 40 0;
#X connect 45 2 39 0;
#X connect 46 0 30 0;
#X restore 92 387 pd stability check;
#X f 9;
#X text 199 10 hp3_butt~;
#X text 338 33 freq [Hz];
#X obj 123 145 f 20;
#X obj 427 186 f 20;
#X obj 83 500 biquad2~ 1.99715 -0.997155 0.998575 -1.99715 0.998575
;
#X obj 83 561 biquad2~ 0.997154 0 0.998577 -0.998577 0;
#X connect 0 0 33 0;
#X connect 4 0 31 0;
#X connect 4 0 32 0;
#X connect 5 0 23 0;
#X connect 6 0 9 0;
#X connect 8 0 23 0;
#X connect 9 0 12 0;
#X connect 9 1 12 1;
#X connect 12 0 24 0;
#X connect 12 1 24 1;
#X connect 12 2 24 2;
#X connect 16 0 33 0;
#X connect 18 0 34 0;
#X connect 19 0 25 0;
#X connect 19 1 25 1;
#X connect 21 0 19 0;
#X connect 23 0 31 0;
#X connect 23 0 32 0;
#X connect 23 1 6 1;
#X connect 23 1 21 1;
#X connect 24 0 28 0;
#X connect 24 1 28 1;
#X connect 24 2 16 2;
#X connect 24 3 16 3;
#X connect 24 4 16 4;
#X connect 25 0 27 0;
#X connect 25 1 18 2;
#X connect 25 2 18 3;
#X connect 26 0 23 0;
#X connect 27 0 18 0;
#X connect 28 0 16 0;
#X connect 28 1 16 1;
#X connect 31 0 6 0;
#X connect 32 0 21 0;
#X connect 33 0 34 0;
#X connect 34 0 1 0;
#X restore 117 323 pd hp_3.order_butterworth_ at _20_Hz;
#X obj 54 186 tgl 25 0 empty empty empty 17 7 0 10 #fce0c4 #000000
#000000 0 1;
#X msg 54 221 \; pd dsp \$1;
#X text 89 190 1.) DSP on/off;
#X obj 116 297 osc~ 0.1;
#X floatatom 286 299 8 0 0 0 - - - 0;
#X msg 286 271 20;
#X obj 43 354 vu 15 120 empty empty -1 -8 0 10 #404040 #000000 1 0
;
#X obj 108 461 dac~;
#X obj 42 326 - 100;
#X obj 117 431 *~ 0;
#X obj 399 246 vradio 15 1 0 11 empty empty empty 0 -8 0 10 #fce0c4
#000000 #000000 0;
#X text 433 243 0 dB;
#X text 419 273 -20 dB;
#X text 419 303 -40 dB;
#X text 419 333 -60 dB;
#X text 419 393 -100 dB;
#X text 419 363 -80 dB;
#X obj 399 469 dbtorms;
#X obj 399 417 * -10;
#X obj 399 444 + 100;
#X obj 399 198 loadbang;
#X msg 399 222 10;
#X floatatom 143 408 5 0 0 0 - - - 0;
#N canvas 0 50 450 250 (subpatch) 0;
#X array \$0-recorder 441000 float 2;
#X coords 0 -0.02 441000 0.02 600 140 1 0 0;
#X restore 564 383 graph;
#X obj 232 348 bng 25 250 50 0 empty empty empty 17 7 0 10 #fce0c4
#000000 #000000;
#X floatatom 305 435 5 0 0 0 - - - 0;
#X text 509 444 0 sec;
#X text 1169 445 10 sec;
#X floatatom 117 275 5 0 0 0 - - - 0;
#X obj 211 382 tabwrite~ \$0-recorder;
#X floatatom 501 375 8 0 0 0 - \$0-upper - 0;
#X floatatom 500 508 8 0 0 0 - \$0-lower - 0;
#X obj 221 495 r \$0-lower;
#X obj 273 517 r \$0-upper;
#X obj 221 594 s \$0-recorder;
#X msg 221 565 bounds 0 \$1 441000 \$2;
#X obj 221 540 pack 0 0;
#X floatatom 35 580 8 0 0 0 - - - 0;
#X obj 35 664 t f f, f 8;
#X obj 35 695 * -1;
#X obj 35 610 abs;
#X obj 35 637 clip 0 1;
#X obj 35 721 s \$0-lower;
#X obj 88 693 s \$0-upper;
#X msg 92 524 0.02;
#X msg 132 523 0.05;
#X msg 38 549 0.1;
#X msg 71 548 0.2;
#X msg 104 548 0.5;
#X text 31 499 array bounds;
#X msg 136 550 1;
#X msg 152 250 0.1;
#N canvas 0 50 450 378 clock_10_sec 0;
#X obj 132 40 inlet;
#X obj 132 172 f;
#X obj 245 174 + 0.1;
#X obj 125 351 outlet;
#X obj 133 71 t b b;
#X msg 234 116 0;
#X obj 132 131 metro 100;
#X msg 163 104 0;
#X msg 131 102 1;
#X obj 132 200 t f f f;
#X obj 213 230 >= 10;
#X obj 213 256 sel 1;
#X obj 130 236 * 10;
#X obj 131 311 / 10;
#X obj 130 263 + 0.49;
#X obj 130 287 int;
#X connect 0 0 4 0;
#X connect 1 0 9 0;
#X connect 2 0 1 1;
#X connect 4 0 8 0;
#X connect 4 1 5 0;
#X connect 5 0 1 1;
#X connect 6 0 1 0;
#X connect 7 0 6 0;
#X connect 8 0 6 0;
#X connect 9 0 12 0;
#X connect 9 1 2 0;
#X connect 9 2 10 0;
#X connect 10 0 11 0;
#X connect 11 0 7 0;
#X connect 12 0 14 0;
#X connect 13 0 3 0;
#X connect 14 0 15 0;
#X connect 15 0 13 0;
#X restore 253 409 pd clock_10_sec;
#X obj 138 498 loadbang;
#X text 159 275 Hz;
#X text 152 227 osc-freq;
#X text 280 251 filter-cutoff;
#X text 349 298 Hz;
#X text 12 86 the 32 bit version has less noise than the 64 bit version
(execution bit width), f 59;
#X text 14 7 ISSUE: floating-point-inaccuracy of low frequency filters
with low frequency signals in pd vanilla, f 100;
#X text 383 177 2.) output gain;
#X msg 342 273 9;
#X text 14 50 try this patch with pd-0.52-1 i386 and amd64 aka 32-bit
and 64-bit and older pd versions, f 69;
#X msg 52 524 0.01;
#X obj 42 298 env~ 2048;
#X floatatom 441 417 5 0 0 0 - - - 0;
#X text 261 352 3.) rec;
#X connect 0 0 10 0;
#X connect 0 0 30 0;
#X connect 1 0 2 0;
#X connect 4 0 0 0;
#X connect 5 0 0 1;
#X connect 6 0 5 0;
#X connect 9 0 7 0;
#X connect 10 0 8 0;
#X connect 10 0 8 1;
#X connect 10 0 65 0;
#X connect 11 0 19 0;
#X connect 18 0 10 1;
#X connect 18 0 23 0;
#X connect 19 0 20 0;
#X connect 19 0 66 0;
#X connect 20 0 18 0;
#X connect 21 0 22 0;
#X connect 22 0 11 0;
#X connect 25 0 53 0;
#X connect 25 0 30 0;
#X connect 29 0 4 0;
#X connect 33 0 37 0;
#X connect 34 0 37 1;
#X connect 36 0 35 0;
#X connect 37 0 36 0;
#X connect 38 0 41 0;
#X connect 39 0 40 0;
#X connect 39 1 44 0;
#X connect 40 0 43 0;
#X connect 41 0 42 0;
#X connect 42 0 39 0;
#X connect 45 0 38 0;
#X connect 46 0 38 0;
#X connect 47 0 38 0;
#X connect 48 0 38 0;
#X connect 49 0 38 0;
#X connect 51 0 38 0;
#X connect 52 0 29 0;
#X connect 53 0 26 0;
#X connect 54 0 45 0;
#X connect 62 0 5 0;
#X connect 64 0 38 0;
#X connect 65 0 9 0;
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pd0.52.1-64bit-biquad2~_0.001bound.png
Type: image/png
Size: 54586 bytes
Desc: not available
URL: <http://lists.puredata.info/pipermail/pd-list/attachments/20220210/eb35fe5a/attachment-0001.png>