执行“在线”线性插值

Question

我有一个问题，我需要对从传感器获取的一些数据进行线性插值（从技术上讲，它是位置数据，但数据的性质并不重要） .我现在在 matlab 中执行此操作，但由于我最终会将此代码迁移到其他语言，因此我希望使代码尽可能简单，并且不使用任何复杂的 matlab 特定/内置函数。

我的实现最初看起来还不错，但是当根据 matlab 的内置 interp1 函数检查我的工作时，我的实现似乎并不完美，我不知道为什么。下面是我在已经完全收集的数据集上使用的代码，但是当我遍历数据时，我的行为就好像我只有当前样本和前一个样本，这反映了我最终将面临的问题。

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

注意“myInterp”函数如下：

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

这是显示我的实现不正确的情节:-(

谁能帮我找出错误在哪里？为什么？我怀疑这与确保查询点位于先前和当前的 x 样本之间有关，但我不确定。

Answer 1

您的代码中的问题是您有时调用myInterp 时使用的qpi 值超出了x_data(kk-1) 和x_data(kk) 的界限。这会导致外推结果无效。

您循环遍历kk 而不是cc 的逻辑让我非常困惑。我会在cc 上编写一个简单的for 循环，这是您要插入的点。对于这些点中的每一个，如有必要，提前kk，使qp(cc)介于x_data(kk)和x_data(kk+1)之间（如果您愿意，可以使用kk-1和kk代替，只需初始化kk=2为了确保kk-1 存在，我发现从kk=1 开始更直观）。

为了简化这里的逻辑，我将qp 中的值限制在x_data 的范围内，这样我们就不需要测试以确保x_data(kk+1) 存在，也不需要确保@987654341 @。如果您愿意，可以添加这些测试。

这是我的代码：

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

如您所见，这种方式的逻辑要简单得多。结果与y_interp_correct 相同。内部while x_data... 循环与外部while 循环的用途相同，是您从数据来源读取数据的地方。