高效执行一维线性插值，无需 for 循环

Question

我正在尝试使用特定的精度在 MATLAB 中执行线性插值。我想知道是否有一种有效的方法可以在 MATLAB 中编写线性插值函数，使其不需要 for 循环并且运行速度非常快？

我想将传入的数据修改为特定的位宽(使用 quantize(( 函数(，然后我还想确保所有中间操作都使用另一个位宽完成。

现在，我使用以下等式在特定点 xq 处执行线性插值，给定点 x 处的值 y。为了清楚起见，我没有包括下面的 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

Answer 1

是的。 您可以做的是，对于每个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，从第二个元素开始的每个x值(由于i+1(都大于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));

请注意 ./ 和 .* 的元素运算符，因为我们正在进行矢量化。

---

警告

我们没有考虑到的是，当您指定的值超出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值。 其他一切都是线性插值的。

Answer 2

为什么不使用内置函数中的 matlabs 来做到这一点呢？

yq = interp1(x, y, xq);

Answer 3

这可能应该是一个评论，但我缺乏业力。在半最新版本的 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 等数组。

Answer 4

由于我仍然不完全确定您希望以哪种方式进行插值本身，因此这里只是关于如何降低算法复杂性的一些建议。在当前形式中，该算法对于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，它为你做了这件事。它可以返回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 倍的加速。