Linear Interpolation code without using interp1

I need to write a code that does linear interpolation. (without using interp1) My function has 3 inputs (n,x,n2). Vector n contains the sample points, and x contains the corresponding values, x(n). Vector n2 contains the coordinates of the query points. My code is this:

function ynew = interpolate(n,x,n2)
for i = 1:length(n2)
    for j = 1:length(n)-1
        if n2(i) >= n(j) & n2(i) <= n(j+1)
            ynew(i) = x(j) + (x(j+1)-x(j))/(n(j+1)-n(j)) * (n2(i)-n(j));
        end
    end
end
end

The problem is that it is much slower than interp1(n,x,n2,'linear'). How can I write a code that works faster?

[Matt J]

What is the purpose of rewriting interp1? If it is a homework exercise, I don't think speed is supposed to matter.

Answer 1

[Torsten]

How can I write a code that works faster?

Allocate ynew = zeros(size(n2)) at the beginning of your code.

Check whether n is monotonically ordered.

Check whether min(n) <= n2 <= max(n).

After you have implemented these three points and you want to speed up your code, edit and understand interp1.

(Means: you shouldn't expect your code to be as fast as a professional code. What you've done on your own is not that bad).

[Radu-Andrei]

Thank you for the suggetions, Torsten. I've tried to look into interp1 but it is a bit overwhelming because there are so many built in functions. At some point interp1 executes this line:

VqLite = matlab.internal.math.interp1(X,V,method,method,Xqcol);

So interp1 uses itself? I don't see the logic.

[Torsten]

Most probably, interp1.m is only a driver function for the user. The real algorithmic implementation of the interpolation algorithm is hidden in this internal function.

Answer 2

[Matt J]

n=1:8;
x=rand(size(n));
n2=[1.8, 3.2, 6.7];

interp1(n,x,n2)

ans = 1×3

0.1287 0.9405 0.3349

interpolate(n,x,n2)

ans = 1×3

0.1287 0.9405 0.3349

function ynew = interpolate(n,x,n2)

  I=discretize(n2,n);

  w=(n2-n(I))./(n(I+1)-n(I));
  
  ynew= x(I+1).*w + x(I).*(1-w);

end

Answer 3

[Matt J]

function ynew = interpolate(n,x,n2)  %Doesn't use interp1

  F=griddedInterpolant(n,x);
  ynew=F(n2);

end