On 05/11, Oleg Nesterov wrote:
>
> Lets look at the code,
>
> cub = it.interpolate_cubic(d,y0,y1,y2,y3)
> with {
> x0 = x1-1;
> x1 = int(id);
> ...
> y0 = rdtable(S, wf, rid(x0, C));
> y1 = rdtable(S, wf, rid(x1, C));
>
> I see two problems:
>
> 1. It always need C=1 or "-ct 1" compiler option, even if the
> input "x" is
> strictly in [r0,r1] range. Otherwise x0=-1 when x=r0.
>
> Not good, C=1 or "-ct 1" implies performance penalty.
>
> 2. Even if C=1, we loss the prcision at r0. In this case y0 =
> y1 and this
> breaks the idea of cubic interpolation.
>
> The same for r1.
I'd like to "complete" this discussion, so let me change the subject
and add Stephane back to CC list.
So, tabulate(r0,r1).cub is very inaccurate near r0/r1 and I am wondering
if is this worth fixing.
One possible (simple) fix is that we can slightly extend [r0,r1] to
[((mid-2)*r0-r1)/(mid-3),((mid-1)*r1-2*r0)/(mid-3)] to ensure that
id(r0) = 1;
id(r1) = mid-2;
this way interpolate_cubic() always has 4 valid points for interpolation.
We can change tabulate().cub but this is not safe in that FX() is not
necessarily well defined on the extended interval.
Or we can add something like
tabulate_cub(C, FX, S, r0, r1, x) = ba.tabulate(C, FX, S, R0, R1, x).cub
with {
R0 = ((S-3)*r0 - r1) / (S-4);
R1 = ((S-2)*r1 - 2*r0) / (S-4);
};
Test-case:
import("stdfaust.lib");
FX = sin;
S = 100;
X0 = 0;
X1 = ma.PI;
tab(x) = ba.tabulate (1, FX, S, X0,X1, x).cub;
cub(x) = tabulate_cub (1, FX, S, X0,X1, x);
maxerr = abs : max ~ _;
line(n, x0,x1) = x0 + (ba.time/n) * (x1-x0);
process = line(50000, X0,X1)<: FX-tab, FX-cub : par(i,2,maxerr);
compiled with "-a plot.cpp"
$ test-plot -n 50001 | tail -1
0.00235060556 5.6214617e-07
As for .lin, I don't think it worth fixing, but we can also do, say,
tabulate_fix(C, FX, S, r0, r1, x) = environment {
val = ba.tabulate(C, FX, S, r0, r1, x).val;
lin = ba.tabulate(C, FX, S, r0, R1, x).lin
with {
R1 = ((S-1)*r1 - r0) / (S-2);
};
cub = ba.tabulate(C, FX, S, R0, R1, x).cub
with {
R0 = ((S-3)*r0 - r1) / (S-4);
R1 = ((S-2)*r1 - 2*r0) / (S-4);
};
};
What do you think?
Oleg.
_______________________________________________
Faudiostream-users mailing list
Faudiostream-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/faudiostream-users
