无需 for 循环即可高效执行一维线性插值
Efficiently perform 1D linear interpolation without for loops
我正在尝试使用特定精度在 MATLAB 中执行线性插值。
我想知道是否有一种有效的方法可以在 MATLAB 中编写线性插值函数，使其不需要 for 循环并且运行速度非常快？
我想将传入的数据修改为特定的位宽（使用 quantize() 函数），然后我还想确保所有的中间操作都是使用另一个位宽完成的。
现在我正在使用以下方程在给定点 x 处的值 y 的特定点 xq 处执行线性插值。为了清楚起见，我没有在下面包含 quantize() 命令。
for j = 1:length(xq)
for k = 1:length(x)
if ((x(k) <= xq(j)) && (xq(j) < x(k+1)))
yq(j) = y(k) + (y(k+1) - y(k))/(x(k+1)-x(k)) * (xq(j)-x(k));
break;
end
end
end
这是与我输入数据的尺寸相匹配的东西。此外，我必须多次执行此线性插值（超过 200 次），因此它需要非常快并且可以与 matlab 中的 interp1 相媲美。 xq 是一个大小为 501x501 的矩阵，其中包含我们希望插值的 x 个位置。 x 和 y 都是长度为 4096 的向量。y 包含复数数据。:
x = linspace(-100.3,100.5,4096);
y = (cos(rand(4096,1))+j*sin(rand(4096,1)))*1/100000;
for i = 1:501
xq(i,:) = (200).*rand(501,1)-100;
end
这应该是评论，但我缺乏因果关系。
在半新版本的 Matlab 中，For 循环并不等于慢。我以前见过这样的情况，人们为了消除任何环路而跳来跳去，结果表现更差。
无论如何，您的循环的很大一部分似乎独立于 j。尝试按照这个思路开始。我怀疑通过删除循环你会看到性能有很大提高，但我很乐意被证明是错误的:)
C = y(1:end-1) + (y(2:end) - y(1:end-1))./(x(2:end)-x(1:end-1))';
for j = 1:length(xq)
for k = 1:length(x)-1
if ((x(k) <= xq(j)) && (xq(j) < x(k+1)))
yq(j) = C(k) * (xq(j)-x(k));
break;
end
end
end
使用示例数据，它的运行速度比您的原始代码快约 25%。另外，请注意第一行末尾的 ' 仅当 x 和 y 的大小不同时才需要，如示例中的 (1,4096) 与 (4096,1)。还要考虑预分配数组，如 xq 和 yq。
是的。您可以做的是对于每个 xq 值，我们可以找到该值适合的区间。我们可以使用两个 bsxfun 调用并利用广播来完成此操作。让我们做一个简单的例子。假设我们有这些 x 和 xq 值：
x = [1 1.5 1.7 2 2.5 2.6 2.8];
xq = [1.2 1.6 2.2];
进行线性插值时，我们的工作是找出每个y点属于x的哪个区间。具体来说，您想要找到索引 i，使得索引 j 处的 xq 的值满足：
xq(j) >= x(i) && xq(j) < x(i+1)
让我们以向量化的方式分别处理布尔表达式的每个部分。对于第一部分，我将创建一个 矩阵 ，其中每一列告诉我 x 的哪些位置满足 xq >= x(i) 不等式。您可以通过以下方式做到这一点：
ind = bsxfun(@ge, xq, x(1:end-1).');
我们得到的是：
ind =
1 1 1
0 1 1
0 0 1
0 0 1
0 0 0
0 0 0
第一列是xq的第一个值，即1.2，下一列是1.6，再下一列是2.2。这告诉我们 x 的哪些值满足 xq(1) > x(i)。请记住，由于 i+1 条件，我们只想检查倒数第二个元素。因此，对于第一列，这意味着 xq(1) >= x(1)，即 xq(1) >= 1，这很好。对于下一列，这意味着 xq(2) >= x(1) 和 xq(2) >= x(2)，即比 1 和 1.5 都 >=，依此类推。
现在我们需要做表达式的另一面：
ind2 = bsxfun(@lt, xq, x(2:end).')
ind2 =
1 0 0
1 1 0
1 1 0
1 1 1
1 1 1
1 1 1
这也是有道理的。对于第一列，这意味着对于 xq = 1，从第二个元素开始（由于 i+1）x 的每个值都 大于 xq 的第一个值，或等效的 xq(1) < x(i+1)。现在，您所要做的就是将这些组合在一起：
ind_final = ind & ind2
ind_final =
1 0 0
0 1 0
0 0 0
0 0 1
0 0 0
0 0 0
因此，对于xq的每个值，行位置精确地告诉我们每个值属于哪个区间。因此，对于 xq = 1.2，这符合区间 1 或 [1,1.2]。对于 xq = 1.6，这适合区间 2 或 [1.5,1.7]。对于 xq = 2.2，这符合区间 4，或 [2,2.5]。您现在要做的就是找到每个 1 的行位置：
[vals,~] = find(ind_final);
vals 现在包含您需要访问 x 以进行插值的位置。这样，你最终得到：
yq = y(vals) + (y(vals+1) - y(vals))./(x(vals+1) - x(vals)).*(xq - x(vals));
注意 element-wise 运算符 ./ 和 .* 因为我们正在做这个向量化。
警告
我们没有考虑到当您指定的值 超出 x 点范围时会发生什么。我会做的是有几个声明来检查这一点。首先，检查是否有任何 xq 值小于第一个 x 点。如果是，让我们将它们设置为第一个输出点。当您的值大于最后一个点并将输出设置为最后一个点时，对另一侧执行相同的操作。因此，您必须执行以下操作：
%// Allocate output array first
yq = zeros(1, numel(xq));
%// Any xq values that are less than the first value, set this to the first
yq(xq < x(1)) = y(1);
%// Any xq values that are less than the last value, set this to the last value
yq(xq >= x(end)) = y(end);
%// Get all of the other values that are not within the above ranges
ind_vals = (xq >= x(1)) & (xq < x(end));
xq = xq(ind_vals);
%// Do our work
ind = bsxfun(@ge, xq, x(1:end-1).');
ind2 = bsxfun(@lt, xq, x(2:end).');
ind_final = ind & ind2;
[vals,~] = find(ind_final);
%// Place output values in the right spots, escaping those values outside the ranges
yq(ind_vals) = y(vals) + (y(vals+1) - y(vals))./(x(vals+1) - x(vals)).*(xq - x(vals))
以上代码应该可以成功进行插值并处理边界条件。
为了看到这个效果，让我们定义一堆 x 和 y 值，以及我们希望插入的 xq 值：
x = [1 1.5 1.7 2 2.5 2.6 2.8];
y = [2 3 4 5 5.5 6.6 7.7];
xq = [0.9 1 1.1 1.2 1.8 2.2 2.5 2.75 2.8 2.85];
使用上面的代码，并通过上面的代码 运行 添加范围检查后，我们得到：
yq =
2.0000 2.0000 2.2000 2.4000 4.3333 5.2000 5.5000 7.4250 7.7000 7.7000
您可以看到任何小于第一个 x 值的 xq 值，我们只是将输出设置为第一个 y 值。任何大于最后一个 x 值的值，我们将输出设置为最后一个 y 值。其他一切都是线性插值的。
为什么不使用 matlabs 内置函数来做这个？
yq = interp1(x, y, xq);
由于我仍然不完全确定您想以哪种方式进行插值本身，这里只是关于如何降低算法复杂性的一些建议。
在其当前形式中，算法对于 m=length(x) 和 n=length(xq) 具有复杂性 O(m*n)。可以将其减少到 O(m*log(m) + n*log(n)).
在低级语言中，您将首先对向量 xq 和 x 进行排序（并重新排列 y 的相应值），然后遍历向量 xq.通过这种方式搜索正确的区间索引 k 比从 k=1 到 length(x) 的搜索更有效。相反，您可以从最后找到的值 k.
开始搜索
不过，MATLAB 确实提供了一个方便的函数 histc，它可以为您完成这项工作。它可以 return 值 xq 所在区间的 k 索引。尽管如此，您仍然需要对 x 和 y 进行排序。
%%// Find the intervals the values lie in.
[x,I] = sort(x);
y = y(I);
[~, ks] = histc(xq, x);
%%// Do the interpolation.
for j = 1:length(xq)
k = ks(j);
yq(j) = y(k) + (y(k+1) - y(k))/(x(k+1)-x(k)) * (xq(j)-x(k));
end
您也许可以将 for 循环矢量化，但由于您没有指定实际执行量化步骤的方式，我将把这个留给您。不过，您可能想看看 rayryeng 的最后一个计算步骤。使用这个 histc 预处理步骤，只需将 vals 换成 ks.
对于问题中列出的数据大小，这应该会给您带来约 30 倍的加速。
我正在尝试使用特定精度在 MATLAB 中执行线性插值。 我想知道是否有一种有效的方法可以在 MATLAB 中编写线性插值函数，使其不需要 for 循环并且运行速度非常快？
我想将传入的数据修改为特定的位宽（使用 quantize() 函数），然后我还想确保所有的中间操作都是使用另一个位宽完成的。
现在我正在使用以下方程在给定点 x 处的值 y 的特定点 xq 处执行线性插值。为了清楚起见，我没有在下面包含 quantize() 命令。
for j = 1:length(xq)
for k = 1:length(x)
if ((x(k) <= xq(j)) && (xq(j) < x(k+1)))
yq(j) = y(k) + (y(k+1) - y(k))/(x(k+1)-x(k)) * (xq(j)-x(k));
break;
end
end
end
这是与我输入数据的尺寸相匹配的东西。此外，我必须多次执行此线性插值（超过 200 次），因此它需要非常快并且可以与 matlab 中的 interp1 相媲美。 xq 是一个大小为 501x501 的矩阵，其中包含我们希望插值的 x 个位置。 x 和 y 都是长度为 4096 的向量。y 包含复数数据。:
x = linspace(-100.3,100.5,4096);
y = (cos(rand(4096,1))+j*sin(rand(4096,1)))*1/100000;
for i = 1:501
xq(i,:) = (200).*rand(501,1)-100;
end
这应该是评论，但我缺乏因果关系。 在半新版本的 Matlab 中，For 循环并不等于慢。我以前见过这样的情况，人们为了消除任何环路而跳来跳去，结果表现更差。
无论如何，您的循环的很大一部分似乎独立于 j。尝试按照这个思路开始。我怀疑通过删除循环你会看到性能有很大提高，但我很乐意被证明是错误的:)
C = y(1:end-1) + (y(2:end) - y(1:end-1))./(x(2:end)-x(1:end-1))';
for j = 1:length(xq)
for k = 1:length(x)-1
if ((x(k) <= xq(j)) && (xq(j) < x(k+1)))
yq(j) = C(k) * (xq(j)-x(k));
break;
end
end
end
使用示例数据，它的运行速度比您的原始代码快约 25%。另外，请注意第一行末尾的 ' 仅当 x 和 y 的大小不同时才需要，如示例中的 (1,4096) 与 (4096,1)。还要考虑预分配数组，如 xq 和 yq。
是的。您可以做的是对于每个 xq 值，我们可以找到该值适合的区间。我们可以使用两个 bsxfun 调用并利用广播来完成此操作。让我们做一个简单的例子。假设我们有这些 x 和 xq 值：
x = [1 1.5 1.7 2 2.5 2.6 2.8];
xq = [1.2 1.6 2.2];
进行线性插值时，我们的工作是找出每个y点属于x的哪个区间。具体来说，您想要找到索引 i，使得索引 j 处的 xq 的值满足：
xq(j) >= x(i) && xq(j) < x(i+1)
让我们以向量化的方式分别处理布尔表达式的每个部分。对于第一部分，我将创建一个 矩阵 ，其中每一列告诉我 x 的哪些位置满足 xq >= x(i) 不等式。您可以通过以下方式做到这一点：
ind = bsxfun(@ge, xq, x(1:end-1).');
我们得到的是：
ind =
1 1 1
0 1 1
0 0 1
0 0 1
0 0 0
0 0 0
第一列是xq的第一个值，即1.2，下一列是1.6，再下一列是2.2。这告诉我们 x 的哪些值满足 xq(1) > x(i)。请记住，由于 i+1 条件，我们只想检查倒数第二个元素。因此，对于第一列，这意味着 xq(1) >= x(1)，即 xq(1) >= 1，这很好。对于下一列，这意味着 xq(2) >= x(1) 和 xq(2) >= x(2)，即比 1 和 1.5 都 >=，依此类推。
现在我们需要做表达式的另一面：
ind2 = bsxfun(@lt, xq, x(2:end).')
ind2 =
1 0 0
1 1 0
1 1 0
1 1 1
1 1 1
1 1 1
这也是有道理的。对于第一列，这意味着对于 xq = 1，从第二个元素开始（由于 i+1）x 的每个值都 大于 xq 的第一个值，或等效的 xq(1) < x(i+1)。现在，您所要做的就是将这些组合在一起：
ind_final = ind & ind2
ind_final =
1 0 0
0 1 0
0 0 0
0 0 1
0 0 0
0 0 0
因此，对于xq的每个值，行位置精确地告诉我们每个值属于哪个区间。因此，对于 xq = 1.2，这符合区间 1 或 [1,1.2]。对于 xq = 1.6，这适合区间 2 或 [1.5,1.7]。对于 xq = 2.2，这符合区间 4，或 [2,2.5]。您现在要做的就是找到每个 1 的行位置：
[vals,~] = find(ind_final);
vals 现在包含您需要访问 x 以进行插值的位置。这样，你最终得到：
yq = y(vals) + (y(vals+1) - y(vals))./(x(vals+1) - x(vals)).*(xq - x(vals));
注意 element-wise 运算符 ./ 和 .* 因为我们正在做这个向量化。
警告
我们没有考虑到当您指定的值 超出 x 点范围时会发生什么。我会做的是有几个声明来检查这一点。首先，检查是否有任何 xq 值小于第一个 x 点。如果是，让我们将它们设置为第一个输出点。当您的值大于最后一个点并将输出设置为最后一个点时，对另一侧执行相同的操作。因此，您必须执行以下操作：
%// Allocate output array first
yq = zeros(1, numel(xq));
%// Any xq values that are less than the first value, set this to the first
yq(xq < x(1)) = y(1);
%// Any xq values that are less than the last value, set this to the last value
yq(xq >= x(end)) = y(end);
%// Get all of the other values that are not within the above ranges
ind_vals = (xq >= x(1)) & (xq < x(end));
xq = xq(ind_vals);
%// Do our work
ind = bsxfun(@ge, xq, x(1:end-1).');
ind2 = bsxfun(@lt, xq, x(2:end).');
ind_final = ind & ind2;
[vals,~] = find(ind_final);
%// Place output values in the right spots, escaping those values outside the ranges
yq(ind_vals) = y(vals) + (y(vals+1) - y(vals))./(x(vals+1) - x(vals)).*(xq - x(vals))
以上代码应该可以成功进行插值并处理边界条件。
为了看到这个效果，让我们定义一堆 x 和 y 值，以及我们希望插入的 xq 值：
x = [1 1.5 1.7 2 2.5 2.6 2.8];
y = [2 3 4 5 5.5 6.6 7.7];
xq = [0.9 1 1.1 1.2 1.8 2.2 2.5 2.75 2.8 2.85];
使用上面的代码，并通过上面的代码 运行 添加范围检查后，我们得到：
yq =
2.0000 2.0000 2.2000 2.4000 4.3333 5.2000 5.5000 7.4250 7.7000 7.7000
您可以看到任何小于第一个 x 值的 xq 值，我们只是将输出设置为第一个 y 值。任何大于最后一个 x 值的值，我们将输出设置为最后一个 y 值。其他一切都是线性插值的。
为什么不使用 matlabs 内置函数来做这个？
yq = interp1(x, y, xq);
由于我仍然不完全确定您想以哪种方式进行插值本身，这里只是关于如何降低算法复杂性的一些建议。
在其当前形式中，算法对于 m=length(x) 和 n=length(xq) 具有复杂性 O(m*n)。可以将其减少到 O(m*log(m) + n*log(n)).
在低级语言中，您将首先对向量 xq 和 x 进行排序（并重新排列 y 的相应值），然后遍历向量 xq.通过这种方式搜索正确的区间索引 k 比从 k=1 到 length(x) 的搜索更有效。相反，您可以从最后找到的值 k.
不过，MATLAB 确实提供了一个方便的函数 histc，它可以为您完成这项工作。它可以 return 值 xq 所在区间的 k 索引。尽管如此，您仍然需要对 x 和 y 进行排序。
%%// Find the intervals the values lie in.
[x,I] = sort(x);
y = y(I);
[~, ks] = histc(xq, x);
%%// Do the interpolation.
for j = 1:length(xq)
k = ks(j);
yq(j) = y(k) + (y(k+1) - y(k))/(x(k+1)-x(k)) * (xq(j)-x(k));
end
您也许可以将 for 循环矢量化，但由于您没有指定实际执行量化步骤的方式，我将把这个留给您。不过，您可能想看看 rayryeng 的最后一个计算步骤。使用这个 histc 预处理步骤，只需将 vals 换成 ks.
对于问题中列出的数据大小，这应该会给您带来约 30 倍的加速。