0-D interpolation, piecewise constant function definition

Specific domains Data handling
1

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)
2

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.

3

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)

image

David

4

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));

Capture d’écran 2023-12-22 à 23.53.52
Capture d’écran 2023-12-22 à 23.53.52722×668 49.1 KB

1 Like
5

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

6

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