Interesting test!
Will have to drill down on this...
On Tue, Dec 22, 2020 at 11:43 AM Oleg Nesterov <o...@redhat.com> wrote:
> Julius, et al,
>
> sorry for delay, I was busy.
>
> And let me apologize in advance, most probably I won't be responsive
> till next week.
>
> On 12/20, Julius Smith wrote:
> >
> > In filters.lib, I just now changed (in faustlibraries / master)
> > smax = 0.9999;
> > to
> > smax = 1.0-ma.EPSILON;
>
> OK, I ran the test below with this change applied,
>
> > The question remains as to which filter form is more accurate.
>
> Hmm. QUITE POSSIBLY I am wrong, I'll recheck. But see below, it seems
> that svf is more accurate, at least with -single option.
>
> > Another thing worth mentioning, by the way, is that tf2snp can be
> > modulated, even using white noise for its coefficients, without affecting
> > signal energy.
>
> I'll write another email about this when I have time...
>
>
> -------------------------------------------------------------------------------
>
> test:
> import("stdfaust.lib");
>
> F = 1;
> w1 = 2*ma.PI*F;
> a1s = -2.0*cos((ma.PI)*-1.0 + ma.PI/4.0);
>
> hp_df = fi.tf2s (1,0,0,a1s,1,w1);
> hp_np = fi.tf2snp (1,0,0,a1s,1,w1);
> hp_svf = fi.svf.hp(F, 1/sqrt(2));
>
> process = 1-1' <: hp_df, hp_np, hp_svf;
>
> Compile with "-single -a plot.cpp", save the output to file:
>
> $ ./test-plot -n 50000 >/tmp/ir
>
> Now, lets plot the difference with the theoretical impulse response
> calculated by maxima:
>
> F0 : 1 $
> SR : 44100 $
>
> H(s,Q) := s^2 / (s^2 + s/Q + 1) $
>
> subst(s = G * (z-1)/(z+1), H(s,Q)) $
> define(TF(z), facsum(ratsimp(%),z)) $
>
> res2(tf, z) := block([t,d, p1,p2, r1,r2],
> t: facsum(ratsimp(tf), z), d: denom(t),
> solve(d, z), [p1,p2] : map(rhs, %%),
> t: num(t) / (coeff(d,z,2) * (z-p1) * (z-p2)),
>
> t * (z-p1)/z, ev(%%, z=p1), r1: map(factor, %%),
> t * (z-p2)/z, ev(%%, z=p2), r2: map(factor, %%),
>
> [r1,r2, p1,p2]
> ) $
>
> [r1,r2, p1,p2]: res2(TF(z),z) $
> /* check, should print zero */
> ratsimp(r1/(1-p1/z) + r2/(1-p2/z) + TF(0) - TF(z));
>
> G: 1/tan(%pi*F0/SR) $ Q: 1/sqrt(2) $
> /* check, should print zero */
> r1-conjugate(r2), eval, bfloat, polarform;
> p1-conjugate(p2), eval, bfloat, polarform;
>
> [i,r,p]: ev([TF(0),r1,p1], eval, bfloat, polarform) $
>
> /*
> declare(n, integer) $
> assume(n >= 0) $
> declare([rm,ra,pm,pa],real) $
> rm*exp(%i*ra) * (pm*exp(%i*pa))^n +
> conjugate(rm*exp(%i*ra)) * conjugate(pm*exp(%i*pa))^n $
> demoivre(%)$ ratsimp(%);
>
> */
>
> define(IR(n),
> charfun(n<1) * i +
> 2 * cabs(r) * cabs(p)^n * cos(carg(r) + carg(p) * n)
> );
>
> /* OK, lets print the difference */
> m: transpose(read_matrix("/tmp/ir")) $
> I: makelist(IR(n),n,0,length(m[3])-1) $
> e: m - matrix(I,I,I) $
> plot2d([[discrete, e[1]], [discrete, e[2]], [discrete, e[3]]],
> [yx_ratio, 0], grid2d) $
>
> See the result:
>
> http://people.redhat.com/onestero/svf/ir-diff.png
>
> 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
