Interp1 or its alternatives, working for matrices without loop?

Hi All,

I have matrices X and Y, with size(x) = [30,30000], and size(y) = [30,30000]. Columnar data in x are different between different columns, and so is y. xi is interpolation vector, with size(xi) = [60, 1]. For linear interpolation, I know I can do:

for i=1: 30000

output(:,i) = interp1(x(:,i), y(:,i), xi);

end

But it is slow to go through the 30000 loops. What is worse, x and y are daily data, and I need to do it annually. Therefore, using loop is a way too slow. I found that shared scripts like ScaleTime exists, which however do not accept x with different columnar values.

Question 1:

Is there a vecterize method to do the interpolation (linear)?

Question 2:

What about vecterized interpolation using cubic or pchip?

Many thanks,

John

Answer 1

[Image Analyst]

So you're going over column by column, creating columns of 60 out of columns of 30. Assuming your xi are just exactly between the elements and exactly on the elements, you can just assume the "in between" interpolations are just the average of the two elements on either side. Thus you can just use a sliding mean filter. Use conv2(y(:,i), [.5; .5]) to give you the "in between" values, then just interleave it with the existing values. Exercise for the reader: interleave the convolution output with the original matrix. (It's like 3 lines of code or so.) Let me know if you can't figure it out and need help.

If your xi are spaced out at non-periodic oddball locations then you can't use conv2.

[John Yang]

Thanks for your reply.

I did not realize that conv2 can be used in this way. Thank you very much for refreshing me.

However, conv2 does not solve the problem in my case. This is because the interpolation, xi, is with equal interval, such that xi = 1:100, while x is not, such that x = [1 5 23 ... 100].

I figured out the problem after a while thinking. The idea is to do translation for x. More specifically a code is as follows:

for i = 1: size(x,2)

x(:,i) = x(:,i)+max(x(:,i))+1; % to keep x monotonic increasing; 1 is to avoid zero x

end

yi = interp1(x(:),y(:),repmat(xi,size(x,2),1));

yi = reshape(yi,size(x,1),size(x,2)); % the answer

In this way, we can use interp1 to deal with matrix. I am very glad, since the solution can help my other related problems. Thanks again!

John