Thanks for the detective work!
I like the idea of setting the pole protection margin according to the
working precision.
The boldest choice would be 1.0 - machineEpsilon.
I think it was Stéphane who gave us ma.EPSILON to use in this kind of
situation (it arose fairly recently for fi.tau2pole's divide-by-zero check).
In filters.lib, I just now changed (in faustlibraries / master)
smax = 0.9999;
to
smax = 1.0-ma.EPSILON;
and now it behaves more on par with the others:
> faust2plot tfilterforms2.dsp ; tfilterforms2 -n 100000 | tail -16
tfilterforms2;
0.00228387304 0.000913831813 0.00317438412; ...
0.00228387304 0.000913831813 0.00317438412; ...
0.00228387304 0.000913831813 0.00317438412; ...
0.00228387304 0.000913831813 0.00317438412; ...
0.00228387304 0.000913831813 0.00317438412; ...
0.00228387304 0.000913831813 0.00317438412; ...
0.00228387304 0.000913831813 0.00317438412; ...
0.00228387304 0.000913831813 0.00317438412; ...
...
> faust2plot -double tfilterforms2.dsp ; tfilterforms2 -n 100000 | tail -16
tfilterforms2;
5.1756238184097469e-12 3.7870455898468107e-12 8.7758306858103949e-12; ...
5.1756238184097469e-12 3.7870455898468107e-12 8.7758306858103949e-12; ...
5.1756238184097469e-12 3.7870455898468107e-12 8.7758306858103949e-12; ...
5.1756238184097469e-12 3.7870455898468107e-12 8.7758306858103949e-12; ...
5.1756238184097469e-12 3.7870455898468107e-12 8.7758306858103949e-12; ...
5.1756238184097469e-12 3.7870455898468107e-12 8.7758306858103949e-12; ...
5.1756238184097469e-12 3.7870455898468107e-12 8.7758306858103949e-12; ...
5.1756238184097469e-12 3.7870455898468107e-12 8.7758306858103949e-12; ...
...
> faust2plot -quad tfilterforms2.dsp ; tfilterforms2 -n 100000 | tail -16
tfilterforms2;
3.71743109e-15 1.88456022e-15 5.46351836e-15; ...
3.71743109e-15 1.88456022e-15 5.46351836e-15; ...
3.71743109e-15 1.88456022e-15 5.46351836e-15; ...
3.71743109e-15 1.88456022e-15 5.46351836e-15; ...
3.71743109e-15 1.88456022e-15 5.46351836e-15; ...
3.71743109e-15 1.88456022e-15 5.46351836e-15; ...
3.71743109e-15 1.88456022e-15 5.46351836e-15; ...
3.71743109e-15 1.88456022e-15 5.46351836e-15; ...
...
The question remains as to which filter form is more accurate.
Another thing worth mentioning, by the way, is that tf2snp can be
modulated, even using white noise for its coefficients, without affecting
signal energy.
That's the usual and original principal benefit of going to the normalized
ladder form - "power-invariant modulatability".
The superior numerical robustness for closely spaced poles was observed
later and documented in the literature by Gray and Markel, as I recall.
- Julius
On Sun, Dec 20, 2020 at 12:07 PM Oleg Nesterov <o...@redhat.com> wrote:
> On 12/20, Dario Sanfilippo wrote:
> >
> > > --- a/filters.lib
> > > +++ b/filters.lib
> > > @@ -1004,7 +1004,7 @@ declare tf2np copyright "Copyright (C) 2003-2019
> by
> > > Julius O. Smith III <jos@ccr
> > > declare tf2np license "MIT-style STK-4.3 license";
> > > tf2np(b0,b1,b2,a1,a2) = allpassnnlt(M,sv) : sum(i,M+1,*(tghr(i)))
> > > with {
> > > - smax = 0.9999; // maximum reflection-coefficient magnitude allowed
> > > + smax = 0.999999999; // maximum reflection-coefficient magnitude
> allowed
> > > s2 = max(-smax, min(smax,a2)); // Project both
> reflection-coefficients
> > > s1 = max(-smax, min(smax,a1/(1+a2))); // into the defined
> > > stability-region.
> > > sv = (s1,s2); // vector of sin(theta) reflection coefficients
> > >
> > >
> > If I'm not wrong, anything above 0.9999999 would be rounded to 1 in
> single
> > precision, right?
>
> Quite possibly, I didn't even bother to check.
>
> In case it was not clear, I didn't try to propose a fix, I just tried to
> identify where does the problem come from.
>
> > Would it be possible to choose different constants based on different
> > options given to the compiler?
>
> Yes, perhaps we should use singleprecision/doubleprecision I dunno. (Can't
> resist I think this feature was a mistake but this is offtopic ;)
>
> Even if we forget about single precision, I simply do not know how much
> we can enlarge this limit. The 0.999999999 value I used is just the "random
> number closer to 1".
>
> > If not, a philosophical question (not really) for these situations might
> > be: should we prioritise single precision or double precision
> > performance/stability?
>
> Good question! please inform me when you know the answer? ;)
>
> Oleg.
>
>
--
"Anybody who knows all about nothing knows everything" -- Leonard Susskind
_______________________________________________
Faudiostream-users mailing list
Faudiostream-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/faudiostream-users
