有效地执行一维线性插值，而无需for循环 [英] Efficiently perform 1D linear interpolation without for loops
问题描述
我正在尝试使用特定的精度在MATLAB中执行线性插值. 我想知道是否有一种有效的方法可以在MATLAB中编写线性插值函数，从而不需要for循环并且运行速度非常快?
I am trying to perform linear interpolation in MATLAB using a specific precision. I was wondering if there is an efficient way to write the linear interpolation function in MATLAB so that it doesn't require for-loops and runs very fast?
我想将传入的数据修改为特定的位宽(使用quantize()函数)，然后我还要确保所有中间操作都使用另一个位宽完成.
I want to modify the incoming data to a specific bitwidth (using the quantize() function), then I also want to make sure that all of the intermediate operations are done using another bitwidth.
现在，在给定点x的值y的情况下，我正在使用以下方程式在特定点xq处执行线性插值.为了清楚起见，我没有在下面包含Quantize()命令.
Right now I am using the following equation to perform linear interpolation at specific points xq given the values y at points x. I have not included the quantize() commands below for clarity.
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保存复杂数据.
Here is something that matches the dimensions of my input data. Also, I have to do this linear interpolation lots of times (over 200), so it needs to be very fast and comparable to interp1 in matlab. xq is a matrix of size 501x501 with x locations we wish to interpolate. x and y are both vectors of length 4096. y holds complex data.:
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
推荐答案
是的.您可以针对每个xq值进行操作，我们可以找到该值适合的间隔.我们可以使用两个 bsxfun 调用，并利用广播来做这.让我们做一个简单的例子.假设我们有以下x和xq值:
Yup. What you can do is for each xq value, we can find which interval the value fits into. We can use two bsxfun calls and take advantage of broadcasting to do this. Let's do a quick example. Let's say we had these x and xq values:
x = [1 1.5 1.7 2 2.5 2.6 2.8];
xq = [1.2 1.6 2.2];
进行线性插值时，我们的工作是找出x的每个y点所属的间隔.具体来说，您要查找索引i，以使索引j处的xq值满足:
When doing linear interpolation, our job is to figure out which interval each y point belongs to for x. Specifically, you want to find the index i, such that a value of xq at index j satisfies:
xq(j) >= x(i) && xq(j) < x(i+1)
让我们以向量化的方式分别处理布尔表达式的每个部分.对于第一部分，我将创建一个矩阵，其中每一列都告诉我x的哪些位置满足xq >= x(i)不等式.您可以通过以下方式做到这一点:
Let's do each part of the Boolean expression separately in a vectorized manner. For the first part, I would create a matrix where each column tells me which locations of x satisfy the xq >= x(i) inequality. You can do that by:
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，均>=，依此类推.
The first column is for the first value of xq, or 1.2, the next column is 1.6, and the next column is 2.2. This tells us which values of x that satisfies xq(1) > x(i). Remember, we only want to check up to the second-last element because of the i+1 condition. Therefore, for the first column, this means that xq(1) >= x(1), which is xq(1) >= 1, and that's great. For the next column, this means that both xq(2) >= x(1) and xq(2) >= x(2), which is both >= than 1 and 1.5 and so on.
现在我们需要做表达式的另一面:
Now we need to do the other side of the expression:
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).现在，您需要做的就是将它们组合在一起:
This also makes sense. For the first column this means that for xq = 1, every value of x from the second element onwards (due to the i+1) are greater than the first value of xq, or equivalently xq(1) < x(i+1). Now, all you have to do is combine these together:
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的每一个的行位置:
Therefore, this for each value of xq, the row location tells us precisely which interval each value falls into. Therefore, for xq = 1.2, this fits into interval 1 or [1,1.2]. For xq = 1.6, this fits into interval 2 or [1.5,1.7]. For xq = 2.2, this fits into interval 4, or [2,2.5]. All you have to do now is find the row locations for each of these 1s:
[vals,~] = find(ind_final);
vals现在包含需要访问x进行插值的位置.这样，您最终得到:
vals now contains where you need to access into x to do your interpolation. As such, you finally get:
yq = y(vals) + (y(vals+1) - y(vals))./(x(vals+1) - x(vals)).*(xq - x(vals));
在进行矢量化时，请注意./和.*的 wise 运算符.
Note the element-wise operators of ./ and .* as we are doing this vectorized.
我们没有考虑到的是，当您指定的值在x点的范围内外部时.我要做的是检查几个语句.首先，检查是否有任何xq值小于第一个x点.如果是的话，让我们将它们设置为第一个输出点.当您拥有的值大于最后一个点并将输出设置为最后一个点时，请对另一端执行相同的操作.因此，您将必须执行以下操作:
What we didn't take into account is what happens when you specify values that are outside the range of the x points. What I would do is have a couple statements that check for this. First, check to see if any xq values are less than the first x point. If they are, let's just set them to the first output point. Do the same thing for the other side when you have values that are greater than the last point and set the output to be the last point. Therefore, you would have to do something like:
%// 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))
上面的代码应该可以成功地进行插值并处理边界条件.
The above code should successfully do the interpolation and handle the boundary conditions.
要查看其工作原理，让我们定义一堆x和y值以及要插值的xq值:
To see this working, let's define a bunch of x and y values, as well as xq values where we wish to interpolate:
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];
使用上面的代码，并在添加了范围检查后运行上面的代码，我们得到:
Using the above, and running through the above code once I added in the range checking, we get:
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值.其他所有内容都是线性插值的.
You can see that any xq values that are less than the first x value, we just set the output to the first y value. Any values that are larger than the last x value, we set the output to the last y value. Everything else is linearly interpolated.
这篇关于有效地执行一维线性插值，而无需for循环的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
