执行 "online" 线性插值

Performing an "online" linear interpolation

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

我的实现最初似乎还不错,但是当根据 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 样本之间有关,但我不确定。

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

你循环 kk 而不是 cc 的逻辑让我很困惑。我会在 cc 上写一个简单的 for 循环,这是你想要插值的点。对于这些点中的每一个,如有必要,提前 kk,使得 qp(cc) 介于 x_data(kk)x_data(kk+1) 之间(您可以使用 kk-1kk 相反,如果您愿意,只需初始化 kk=2 以确保 kk-1 存在,我发现从 kk=1 开始更直观)。

为了简化这里的逻辑,我将 qp 中的值限制在 x_data 的范围内,这样我们就不需要测试来确保 [=22] =] 存在,也不存在 x_data(1)<pq(cc)。如果您愿意,可以添加这些测试。

这是我的代码:

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 循环的目的相同,都是您从任何地方读取数据的地方。