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.
_______________________________________________
Faudiostream-users mailing list
Faudiostream-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/faudiostream-users
