执行 "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-1
和 kk
相反,如果您愿意,只需初始化 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
循环的目的相同,都是您从任何地方读取数据的地方。
我有一个问题,我需要对一些数据进行线性插值 ,因为它是从传感器获取的(从技术上讲,它是位置数据,但数据的性质并不真的很重要)。我现在正在 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-1
和 kk
相反,如果您愿意,只需初始化 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
循环的目的相同,都是您从任何地方读取数据的地方。