无需循环即可高效执行1D线性插值-Java 学习之路

我正在尝试使用特定的精度在MATLAB中执行线性插值 . 我想知道是否有一种有效的方法在MATLAB中编写线性插值函数，这样它就不需要for-loops并且运行速度非常快？

我想将输入数据修改为特定的位宽（使用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

4 回答

3
对 . 你可以做的是每个 xq 值，我们可以找到值适合的区间 . 我们可以使用两个bsxfun呼叫并利用广播来实现这一目标 . 让's do a quick example. Let'说我们有 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)
```
让我们以矢量化的方式分别完成布尔表达式的每个部分 . 对于第一部分，我将创建一个 matrix ，其中每列告诉我 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 . 这告诉我们满足 xq(1) > x(i) 的 x 的哪些值 . 请记住，由于 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 ）都是 greater ，而不是 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 的每个值，行位置告诉我们 precisely 每个值落入哪个区间 . 因此，对于 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 运算符，因为我们正在进行向量化 .

警告

我们没有考虑的是当您指定 outside 点范围的值时会发生什么 . 我要做的是有几个声明检查这个 . 首先，检查是否有任何 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 值 . 其他所有内容都是线性插值的 .
回复于 2024-05-03T10:44:36+08:00
1
为什么不使用内置功能的matlabs来做到这一点？
```
yq = interp1(x, y, xq);
```
回复于 2024-05-03T10:44:36+08:00
1
这应该是一个评论，但我缺乏业力 . 在半新版本的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这样的数组 .
回复于 2024-05-03T10:44:36+08:00
1
由于我还不完全确定你想以哪种方式进行插值，这里只是一些关于如何降低算法复杂性的建议 . 在目前的形式中，对于 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 -loop矢量化，但是你没有把它留给你 . 你可能想看一下rayryeng的最后计算步骤 . 使用这个 histc 预处理步骤，只需将 vals 替换为 ks 即可 .

对于问题中列出的数据大小，这应该可以为您提供~30倍的加速 .
回复于 2024-05-03T10:44:36+08:00

无需循环即可高效执行1D线性插值

4 回答

警告

相关问题