Multidimensional interpolation

Hello,

I have the following multidimensional interpolation that I am trying to do. I think it is subtly different from interpn that Matlab has as an intrinsic function.

I have two series of arrays, A and B.

Series A corresponds to a combination of variables that were used as inputs to a calculation (one who's function is unknown and is unimportant):

[Wa,Xa,Ya,Za] = ndgrid(1:1:5,30:30:360,2:2:20,1:1:10); You'll notice that this corresponds to 6000 calculations as length(Wa(:)) = 6000.

Series B has exactly the same size and dimensions as A except the values in B are the results of the calculation. [Wb,Xb,Yb,Zb] = ngrid(1:1:5,[some numbers],[some numbers],[some numbers])

What I want to do is:

1) For an arbitary value of W,X,Y,Z, work out what its position (or address) would be in arrays for series A. Even if it does not correspond to an exact elemnt address (i.e. it lies slightly beteen elements).

2) Using the address obtained in (1), or addresses of the elements surrounding the position in (1), work out what the corresponding values of W,X,Y,Z woud be for that position in the arrays in series B.

Any suggestions or examples of this kind of calculation? The idea is that, having done this series of calculations in A (which takes a long time), I don't have to repeat them for new data - I can just interpolate them from existing results in B.

Best regards,

Ben

Answer 1

[Walter Roberson]

For any A one coordinate, subtract the minimum coordinate value for that dimension, divide by the A step size, and add 1 to get the corresponding index in B. Which you would convert to a B coordinate by subtracting one, multiplying by the B step size, and adding the B minimum.

If minXa = min(30:30:360) (i.e., 30), and stepXa is the increment (i.e., 30), and minXb = min([some numbers #1]), and stepXb is the increment in B for X, then

(((Xa - minXa)/stepXa + 1) - 1) * stepXb + minXb

simplifies to

   (Xa - minXa) / stepXa * stepXb + minXb

If you prefer, this can be expressed as

tsX = stepXb / stepXa;
tincrX = minXb - minXa * tsX;

then

Xa * tsX + tincrX;

to within round-off error.

[Ben Williams]

Hi,

Thanks for the info! I should probably highight that I am specifically interested in cases where the new data for W,X,Y,Z does NOT correspond to an exact address - i.e. there will be some linear interpolation required between the closest surrounding elements. But I'll have a play with the code you suggest. Thanks, Ben

[Walter Roberson]

For linear interpolation, you would have

floor((Xa - minXa) / stepXa) * stepXb + minXb

and

ceil((Xa - minXa) / stepXa) * stepXb + minXb

as the nearest neighbours, with the two being exactly the same if (Xa - minXa) / stepXa comes out as an exact integer.

You can work from those two together with the unconstrained value (Xa - minXa) / stepXa * stepXb + minXb to do the linear interpolation.

I believe that if you look at these equations more closely, you will find that (1 - mod(Xa-minXa, stepXa) / stepXa) gives you the weighting factor you need.