Performing an "online" linear interpolation

Question

I have a problem where I need to do a linear interpolation on some data as it is acquired from a sensor (it's technically position data, but the nature of the data doesn't really matter). I'm doing this now in matlab, but since I will eventually migrate this code to other languages, I want to keep the code as simple as possible and not use any complicated matlab-specific/built-in functions.

My implementation initially seems OK, but when checking my work against matlab's built-in interp1 function, it seems my implementation isn't perfect, and I have no idea why. Below is the code I'm using on a dataset already fully collected, but as I loop through the data, I act as if I only have the current sample and the previous sample, which mirrors the problem I will eventually face.

%make some dummy data
np = 109; %number of data points for x and y
x_data = linspace(3,98,np) + (normrnd(0.4,0.2,[1,np]));
y_data = normrnd(2.5, 1.5, [1,np]);

 %define the query points the data will be interpolated over
qp = [1:100];

kk=2; %indexes through the data
cc = 1; %indexes through the query points

qpi = qp(cc); %qpi is the current query point in the loop
y_interp = qp*nan; %this will hold our solution

while kk<=length(x_data)
    kk = kk+1; %update the data counter

    %perform online interpolation
    if cc<length(qp)-1
        if qpi>=y_data(kk-1) %the query point, of course, has to be in-between the current value and the next value of x_data
            y_interp(cc) = myInterp(x_data(kk-1), x_data(kk), y_data(kk-1), y_data(kk), qpi);
        end

        if qpi>x_data(kk), %if the current query point is already larger than the current sample, update the sample
            kk = kk+1;
        else %otherwise, update the query point to ensure its in between the samples for the next iteration
            cc = cc + 1;
            qpi = qp(cc);

            %It is possible that if the change in x_data is greater than the resolution of the query
            %points, an update like the above wont work. In this case, we must lag the data
            if qpi<x_data(kk),
                kk=kk-1;
            end
        end
    end
end

%get the correct interpolation
y_interp_correct = interp1(x_data, y_data, qp);

%plot both solutions to show the difference
figure;
plot(y_interp,'displayname','manual-solution'); hold on;
plot(y_interp_correct,'k--','displayname','matlab solution');
leg1 = legend('show');
set(leg1,'Location','Best');
ylabel('interpolated points');
xlabel('query points');

Note that the "myInterp" function is as follows:

function yi = myInterp(x1, x2, y1, y2, qp)
%linearly interpolate the function value y(x) over the query point qp
yi = y1 + (qp-x1) * ( (y2-y1)/(x2-x1) );
end

And here is the plot showing that my implementation isn't correct :-(

Can anyone help me find where the mistake is? And why? I suspect it has something to do with ensuring that the query point is in-between the previous and current x-samples, but I'm not sure.

Answer 1

The problem in your code is that you at times call myInterp with a value of qpi that is outside of the bounds x_data(kk-1) and x_data(kk). This leads to invalid extrapolation results.

Your logic of looping over kk rather than cc is very confusing to me. I would write a simple for loop over cc, which are the points at which you want to interpolate. For each of these points, advance kk, if necessary, such that qp(cc) is in between x_data(kk) and x_data(kk+1) (you can use kk-1 and kk instead if you prefer, just initialize kk=2 to ensure that kk-1 exists, I just find starting at kk=1 more intuitive).

To simplify the logic here, I'm limiting the values in qp to be inside the limits of x_data, so that we don't need to test to ensure that x_data(kk+1) exists, nor that x_data(1)<pq(cc). You can add those tests in if you wish.

Here's my code:

qp = [ceil(x_data(1)+0.1):floor(x_data(end)-0.1)];
y_interp = qp*nan; % this will hold our solution

kk=1; % indexes through the data
for cc=1:numel(qp)
   % advance kk to where we can interpolate
   % (this loop is guaranteed to not index out of bounds because x_data(end)>qp(end),
   %    but needs to be adjusted if this is not ensured prior to the loop)
   while x_data(kk+1) < qp(cc)
      kk = kk + 1;
   end
   % perform online interpolation
   y_interp(cc) = myInterp(x_data(kk), x_data(kk+1), y_data(kk), y_data(kk+1), qp(cc));
end

As you can see, the logic is a lot simpler this way. The result is identical to y_interp_correct. The inner while x_data... loop serves the same purpose as your outer while loop, and would be the place where you read your data from wherever it's coming from.