Excellent!
+1
Thanks very much!
- Julius
On Fri, Nov 25, 2022 at 12:48 PM Oleg Nesterov <o...@redhat.com> wrote:
> Hello,
>
> Yesterday I played with ba.tabulate() and I think we can have a better
> helper which approximates an unary function using chebyshev polynomials
> more accurately:
>
> chebtab(CK, FX, NX,CD, X0,X1, x) = y with {
> ck(0) = _;
> ck(1) = max(0) : min(NX-1);
>
> NC = CD + 1;
> DX = (X1 - X0) / NX;
> nx = (x - X0) / DX : int : ck(CK);
> xc(n) = X0 + DX * (n + 1/2);
> ch(i) = chtab(NC * nx + i);
>
> y = (x - xc(nx)) * 2/DX <: sum(i,NC, ch(i) *
> ma.chebychev(i));
>
> mapfx(nx, x) = FX(xc(nx) + DX/2 * x);
>
> gench(nx, nc) = (1+(nc!=0))/NC * sum(k,NC,
> mapfx(nx, cos(ma.PI*(k+1/2)/NC)) *
> cos(ma.PI*nc*(k+1/2)/NC));
>
> chtab = rdtable(NX*NC, (ba.time <: int(/(NC)), %(NC) :
> gench));
> };
>
> Usage:
>
> _ : chebtab(CK, FX, NX,CD, X0,X1) : _
>
> Where
> CK: force the index in [X0, X1] range
>
> FX: unary function
>
> NX: number of segments
> CD: degree of the chebyshev polynomials
>
> X0: range minimal value
> X1: range maximal value
>
> chebtab() splits the interval [X0,X1] into NX segments, and for each
> segment
> it calculates CD+1 coefficients for ma.chebychev(0) .. ma.chebychev(CD).
> IOW,
> it uses rdtable as 2-dimensional array [NX][CD+1], so its size is NX *
> (CD+1).
>
> Lets compare chebtab() with ba.tabulate().cub, test-case:
>
> CK = 0;
> FX = sin;
> NX = 50;
> CD = 3;
> X0 = 0;
> X1 = ma.PI;
>
> // Both ch() and tb() create rdtable of size = 200, both use
> // cubic polynomials, so this comparison is more or less fair.
> ch(x) = chebtab(CK, FX, NX, CD, X0,X1, x);
> tb(x) = ba.tabulate(CK, FX, NX*(CD+1), X0,X1, x).cub;
>
> maxerr = abs : max ~ _;
> line(n, x0,x1) = x0 + (ba.time%n)/n * (x1-x0);
> process = line(50000, X0,X1) <: FX-tb,FX-ch : par(i,2,maxerr);
>
> Now,
>
> $ cheb-plot -n 50000 | tail -1
>
> 0.00349552557 5.07155429e-09
>
> as you can see chebtab() is much more accurate. And according to my quick
> testing it is even a little bit faster but I won't insist on that.
>
> Just in case... tabulate().cub can be improved, it.interpolate_cubic() is
> simply wrong near X0 or X1, but chebtab() will be more accurate anyway.
>
> Do you think it can live in basics.lib ?
>
> 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
