Custom Matrix Interpolation for every 10 numbers in each column

Dear all, I've got a matrix (Mx) 400 x 410. I need to interpolate it to obtain a 2000x2050 matrix. The final result should be almost equal to Fx (attached matrix. Please note, Fx should be 2000x2050 but the file is too big, so I cropped it to 1000x1000). Given the nature of the matrix (latitude file from a satellite scene), a linear/bilinear/cubic/etc. interpolation introduces an error which makes the interpolated file useless. I figured out what works for me, and I reproduced it in a small sample in Excel (file attached), but I don't have the competencies to write a MATLAB script that can do the same. These are the steps that I need, to obtain the desired result:

1) EDIT*: First, a horizontal interpolation for every row. Taking differences between adjacent points and using 1/5th of those difference to compute the 4 interior points. Produces a matrix (400x2050):

2) Arrange the matrix for vertical interpolation (this step I guess is not needed in MATLAB, it is just to visually explain the theory):

3) I need to compute the interpolation in blocks of 10 numbers. That would be, in the example above, A1 to A10, A11 to A20, etc. This should be done by applying the method below:

x = (A6 - A1)/5;

A1 = A1; This has to remain unchanged.

A2 = A1+x;

A3 = A2+x;

A4 = A3+x;

A5 = A4+x;

A6 = A6; This has to remain unchanged.

A7 = A6+x;

A8 = A7+x;

A9 = A8+x;

A10 = A9+x;

At this point, I need to move to the second block of 10 (i.e., A11 to A20):

x' = (A16 - A11)/5;

A11 = A11; This has to remain unchanged.

A12 = A11+x';

A13 = A12+x';

A14 = A13+x';

A15 = A14+x';

A16 = A16; This has to remain unchanged.

A17 = A16+x';

A18 = A17+x';

A19 = A18+x';

A20 = A19+x';

Then, I should move to the third block of 10 (i.e., A21 to A30)... And so on...

I should repeat the same procedure across all the columns. The final result should look like:

Your help is massively appreciated!

Answer 1

[Matt J]

A=reshape(1:40,[],2); 
A(10:end,:)=A(10:end,:)*2;
A=A(1:5:end,:); %hypothetical input

M=5; %upsampling factor


A=kron(A,[1;zeros(M-1,1)]);

N=size(A,1);
fcn=@(z) repelem(z,2*M,1);

B=diff(A(1:M:end,:),1,1);
B=B(1:2:end,:)/M;
C=repmat((0:2*M-1)' ,height(B),1);
B=fcn(B);



result=fcn(A(1:2*M:end,:))+B.*C;

table(A,result)

ans = 20×2 table
       A         result 
    ________    ________

     1    21     1    21
     0     0     2    22
     0     0     3    23
     0     0     4    24
     0     0     5    25
     6    26     6    26
     0     0     7    27
     0     0     8    28
     0     0     9    29
     0     0    10    30
    22    62    22    62
     0     0    24    64
     0     0    26    66
     0     0    28    68
     0     0    30    70
    32    72    32    72

[Simone A.]

Hi Matt, Thanks for your reply.

However, assuming A being a matrix (400 x 410), the output I get (result), is still a matrix of 400 x 410 (it should be a 2000 x 2050 matrix). Additionally, the values that i get are not consistent with my wanted output (see Fx.mat file attached).

Thanks for your help!

[Matt J]

However, assuming A being a matrix (400 x 410), the output I get (result), is still a matrix of 400 x 410.

That's not what I get. I get 2000 rows for sure:

A=rand(400,410);

M=5; %upsampling factor


A=kron(A,[1;zeros(M-1,1)]);

N=size(A,1);
fcn=@(z) repelem(z,2*M,1);

B=diff(A(1:M:end,:),1,1);
B=B(1:2:end,:)/M;
C=repmat((0:2*M-1)' ,height(B),1);
B=fcn(B);



result=fcn(A(1:2*M:end,:))+B.*C;

whos result

  Name           Size               Bytes  Class     Attributes

  result      2000x410            6560000  double              

[Simone A.]

Dear @Matt J,

Thanks a lot for getting back to me! I don't really know what happened before (I surely made a mistake copying the code), but now it actually works!

Would you be able to help me obtaining the final result?

Now that I've got a 2000x400 matrix, I need to interpolate it to create a 2000x2050 matrix.

This (horizontal) interpolation should be done with the following method:

Taking differences between adjacent points and using 1/5th of those difference to compute the 4 interior points.

47.1493 47.1497 47.1501 47.1505 47.1508 47.1512 47.1514 47.1517 47.1519 47.1521 47.1524 
47.1306 47.1310 47.1314 47.1318 47.1323 47.1327 47.1329 47.1332 47.1335 47.1338 47.1340 
47.1119 47.1123 47.1128 47.1132 47.1137 47.1141 47.1145 47.1148 47.1151 47.1154 47.1157 
47.0931 47.0936 47.0941 47.0946 47.0951 47.0956 47.0960 47.0963 47.0967 47.0970 47.0974 
47.0744 47.0749 47.0755 47.0760 47.0765 47.0771 47.0775 47.0779 47.0783 47.0787 47.0790 
47.0556 47.0562 47.0568 47.0574 47.0580 47.0585 47.0590 47.0594 47.0598 47.0603 47.0607 
47.0369 47.0375 47.0381 47.0388 47.0394 47.0400 47.0405 47.0410 47.0414 47.0419 47.0424 
47.0182 47.0188 47.0195 47.0202 47.0208 47.0215 47.0220 47.0225 47.0230 47.0235 47.0241 
46.9994 47.0001 47.0008 47.0015 47.0022 47.0029 47.0035 47.0041 47.0046 47.0052 47.0057 
46.9807 46.9814 46.9822 46.9829 46.9837 46.9844 46.9850 46.9856 46.9862 46.9868 46.9874 
47.0623 47.0626 47.0630 47.0633 47.0637 47.0640 47.0643 47.0645 47.0647 47.0649 47.0652

The result of this step, taking Mx (attached) as an example, should return the following:

Thanks a lot, I really appreciate your help!

[Matt J]

For the final step, you can just use interp1

A=rand(2000,410);

M=5;
n=size(A,2);

A=interp1(0:n-1,A',0:1/M:n-1/M,'linear','extrap')';
whos A

  Name         Size                 Bytes  Class     Attributes

  A         2000x2050            32800000  double              

[Simone A.]

That's great too! Thanks a lot Matt!

Answer 2

[Walter Roberson]

[Mxrows, Mxcols] = size(Mx);
output = interp1(1:Mxrows, interp1(1:Mxcols, Mx.', 1:.1:Mxcols).', 1:.1:Mxrows);

[Simone A.]

Hi Walter, thanks a lot for your help. I tried to run this code, but I get an output of 3991x4091 instead of 2000x2050. Do you know why?

[Walter Roberson]

output = interp1(1:Mxrows, interp1(1:Mxcols, Mx.', linspace(1,Mxcols,2000)).', linspace(1,Mxrows,2050));

However, this disagrees with the interpolation diagram that you had earlier posted and then removed.

x = (A6 - A1)/5; A1 = A1; This has to remain unchanged. A2 = A1+x; A3 = A2+x; A4 = A3+x; A5 = A4+x; A6 = A6;

1. So x = (original second - original first)/5
2. new first = original first point
3. new second = new first + x = original first + 1*(original second - original first)/5
4. new third = new second + x = original first + 2*(original second - original first)/5
5. new fourth = new third + x = original first + 3*(original second - original first)/5
6. new fifth = new fourth + x = original first + 4*(original second - original first)/5
7. new sixth = original second = original first + 5*(original second - original first)/5

Number of output points = 6. If you were to continue your steps for another cycle, you would add 5 more rows, getting out 11 total. Another cycle, another 5 points... The total number of output points is 1 + 5*number of cycles.... Which is never a number divisible by 5 !!

A7 = A6+x; A8 = A7+x;

You are using the same increment here as you used for the previous batch of 5 !! Are you sure that is what you want ?? This is not what you explained in those diagrams you had posted!

For example if the input were [1 21 221] then your text says you would calculate x = (21-1)/5 = 4, A1=1, A2=A1+x=1+4=5, A3=A2+x=5+4=9, A4=A3+x=9+4=13, A5=A4+13+4=17, A6=21 (copied), A7=A6+x=21+4=25, A8=A7+x=25+4=29, A9=A8+x=29+4=33, A10=A9+x=33+4=37, A11=221 (copied) . Are you sure you want to take the difference between two adjacent points, divide by 5, and then use that to predict the next 10 points?? Rather than taking differences between adjacent points and using 1/5th of those difference to compute the 4 interior points??

[Simone A.]

Dear Walter, once again, thanks for your reply. I am not sure why the diagrams went missing, however, I uploded them back.

"A7 = A6+x; A8 = A7+x;

You are using the same increment here as you used for the previous batch of 5 !! Are you sure that is what you want ??"

Yes, this is what I need indeed. I need to use the same increment from A1 to A10 (nemely x=(A6-A1)/5), keeping A1 and A6 unchanged... Then from A11 to A20 (nemely x'=(A16-A11)/5), keeping A11 and A16 unchanged.... and so on.

As you correctly stated, this could also work as:

1. So x = (original second - original first)/5
2. new first = original first point
3. new second = new first + x = original first + 1*(original second - original first)/5
4. new third = new second + x = original first + 2*(original second - original first)/5
5. new fourth = new third + x = original first + 3*(original second - original first)/5
6. new fifth = new fourth + x = original first + 4*(original second - original first)/5
7. new sixth = original second = original first + 5*(original second - original first)/5
8. new seventh= new sixth + x = original first + 6*(original second - original first)/5
9. new eigth= new seventh + x = original first + 7*(original second - original first)/5
10. new ninth= new eigth + x = original first + 8*(original second - original first)/5
11. new tenth= new nineth + x = original first + 9*(original second - original first)/5

Then, I should move to the second batch of 10 (from A11 to A20, with x' = (A16 - A11)/5;), repeating the same procedure above.

'Are you sure you want to take the difference between two adjacent points, divide by 5, and then use that to predict the next 10 points??'

Yes, as by taking differences between adjacent points and using 1/5th of those difference to compute the 4 interior points, considering the peculiarity of my data and the wanted output, I get an "erroneous interpolation".

I hope this was a bit clearer. I also replaced the diagram, hoping that they can clarify further. The Excel file is also attached providing additional support.

Thank you very much indeed.

[Walter Roberson]

If the input were [1 21 221] then could you work through the desired outputs please?

[Simone A.]

Hi Walter, if my input were [1 21 221] in [A1 A6 A11],respectively, the outputs you suggested are correct.

Just to confirm, the first batch will go from A1 to A10, with A11 being part of the second batch, running from A11 to A20.

EDIT*: I noticed that step 1 (horizontal linear interpolation, for example using interp1), won't produce the result I need. Instead, could step 1 be adjusted to interpolate the values by taking differences between adjacent points and using 1/5th of those difference to compute the 4 interior points?

I amended the original post accordingly.

Thanks a bunch!