dear all,

I’d like simply to interpolate data following piecewise constant method, as illustrated in code below. I looked into available interpolation functions but I didn’t found proper configuration to do this. I would like extrapolation option “edgevalue” from interp1() but for interpolation… more or less the same behaviour as in plot2d2() graphic function but as function defintion.

The code is doing what I want but the loop makes it not performant when dealing with huge datasets.

Do you have any suggestion to make it clever and faster ?

Thanks

David

```
I=1:10
inom= [1 2 5 7 8 10]
innom= setdiff(I,inom)
iref= innom
for k=1:length(innom), iref(k)=inom(max(find(innom(k)>inom))) ; end
I(innom)=I(iref)
```

Hello David,

A piecewise constant function f on [a,b[ is defined e.g. by a sequence (x_i)_{i=1..n} where a=x_1<x_2<\dots<x_n=b and values (y_i)_{i=1..n-1} such that \forall x\in[x_i,x_{i+1}[, f(x)=y_i. Where are defined x and y in your code?

If you have huge sequences a built-in function should be preferable: you can use `interp1`

in nearest neighbor mode but in that case you should define abscicae values (u_i)_{i=0..n} such that (u_i+u_{i-1})/2=x_{i}, i=1\dots n, to synchronize the discontinuity (the nearest neighbor interpolant value changes at the middle of intervals).

Since the choice is not unique you have to define arbitrarily u_0, for example u_0=x_1-(x_2-x_1)/2 (it places u_1 right in the middle of [x_1,x_2]).

Then you have u_{i}=2x_{i}-u_{i-1}, i=1\dots n, and you can feed `interp1`

with (u_i)_{i=0..n} and (y_i)_{i=0..n}, where y_0 is arbitary, i.e. y_0=y_1.

EDIT : the method is not general, i.e. there are cases when there does not exist a nearest neighbor equivalent of an arbitrary piecewise C0 function…

S.

Hello Stephane,

the previous code was very specific to current problem,not general at all, I was focused on the index and somewhat confusing myself.

I tried to generalize a little bit but still far from being robust ! :

```
function yp= interp0(x,y,xp)
x_=xp
for k=1:length(x_), x_(k)=max(find(xp(k)>=x)) ; end
yp=y(x_)
endfunction
x=[1 2 5 7 8 10]
y= [1 2 5 7 8 10]+1.6
xp= [1:10] +0.1
yp= interp0(x,y,xp)
plot2d2(x,y)
plot2d(xp,yp,-2)
```

David

Sorry for my first complicated/not working solution. The following one is much clever…

```
function yp= interp0(x,y,xp)
i = floor(interp1(x,1:length(x),xp,"linear","edgevalue"));
yp = y(i);
endfunction
x = [1 2 5 7 8 10]
y = [1 2 5 7 8 10]+1.6
xp = linspace(-1,12,10000);
plot(xp,interp0(x,y,xp));
```

1 Like

Hi Stephane,

thank you for this simple and fast solution !

May this implementation of piecewise constant interpolation feature be included in future release of Scilab as an extra option of `interp1()`

?

best wishes,

David

Hi ,

looking for highest speed, I experimented successfully to use native function `linear_interpn`

instead of `interp1`

macro, significantly faster on large dataset.

```
function yp= interp0(x,y,xp)
i = floor(linear_interpn(xp,x,1:length(x),"C0"))
yp = y(i)
endfunction
```

So the suggestion is to add the piecewise constant interplotation feature rather to `linear_interpn`

.

best,

David