Converting MATLAB's interp1 to Python interp1d

Question

I'm converting a MATLAB code into a Python code.

The code uses the function interp1 in MATLAB. I found that the scipy function interp1d should be what I'm after, but I'm not sure. Could you tell me if the code, I implemented is correct? My Python version is 3.4.1, MATLAB version is R2013a. However, the code has been implemented around 2010].

MATLAB:

S_T = [0.0, 2.181716948, 4.363766232, 6.546480392, 8.730192373, ...
      10.91523573, 13.10194482, 15.29065504, 17.48170299, 19.67542671, ...
      21.87216588, 24.07226205, 26.27605882, 28.48390208; ...
     1.0, 1.000382662968538, 1.0020234819906781, 1.0040560245904753, ...
       1.0055690037530718, 1.0046180687475195, 1.000824223678225, ...
       0.9954866694014762, 0.9891408937764872, 0.9822543350571298, ...
       0.97480163751874, 0.9666158376141503, 0.9571711322843011, ...
       0.9460998105962408; ...
     1.0, 0.9992731388936672, 0.9995093132493109, 0.9997021748479805, ...
       0.9982835412406582, 0.9926319477117723, 0.9833685776596993, ...
       0.9730725288209638, 0.9626092685176822, 0.9525234896714959, ...
       0.9426698515488858, 0.9326788630704709, 0.9218100196936996, ...
       0.9095717918978693];

S = transpose(S_T);

dist = 0.00137;
old = 15.61;
ll = 125;
ref = 250;
start = 225;
high = 7500;
low = 2;

U = zeros(low,low,high);

for ii=1:high
    g0= start-ref*dist*ii;
    g1= g0+ll;
    if(g0 <=0.0 && g1 >= 0.0)
        temp= old/2*(1-cos(2*pi*g0/ll));
        for jj=1:low
            U(jj,jj,ii)= temp;
        end
    end
end

for ii=1:low
 S_mod(ii,1,:)=interp1(S(:,1),S(:,ii+1),U(ii,ii,:),'linear');
end

Python:

import numpy
import os
from scipy import interpolate

S = [[0.0, 2.181716948, 4.363766232, 6.546480392, 8.730192373, 10.91523573, 13.10194482, 15.29065504, \
      17.48170299, 19.67542671, 21.87216588, 24.07226205, 26.27605882, 28.48390208], \
     [1.0, 1.000382662968538, 1.0020234819906781, 1.0040560245904753, 1.0055690037530718, 1.0046180687475195, \
      1.000824223678225, 0.9954866694014762, 0.9891408937764872, 0.9822543350571298, 0.97480163751874, \
      0.9666158376141503, 0.9571711322843011, 0.9460998105962408], \
     [1.0, 0.9992731388936672, 0.9995093132493109, 0.9997021748479805, 0.9982835412406582, 0.9926319477117723, \
      0.9833685776596993, 0.9730725288209638, 0.9626092685176822, 0.9525234896714959, 0.9426698515488858, \
      0.9326788630704709, 0.9218100196936996, 0.9095717918978693]]

dist = 0.00137
old = 15.61
ll = 125
ref = 250
start = 225
high = 7500
low = 2

U = [numpy.zeros( [low, low] ) for _ in range(high)]

for ii in range(high):
    g0 = start - ref * dist * (ii+1)   
    g1 = g0 + ll
    if g0 <=0.0 and g1 >= 0.0:
        for jj in range(low):
            U[ii][jj,jj] = old / 2 * (1 - numpy.cos( 2 * numpy.pi * g0 / ll) )

S_mod = []

for jj in range(high):
    temp = []
    for ii in range(low):
        temp.append(interpolate.interp1d( S[0], S[ii+1], U[jj][ii,ii]))

    S_mod.append(temp)

Matlab used to have a clearer documentation for `interp1`. You can have a look at [my answer in this question](http://stackoverflow.com/questions/27406562/collapse-mean-data-in-matlab-with-respect-to-a-different-set-of-data/27407674#27407674) for an explanation of what `interp1` is doing (there is also a [graphical explanation](http://i.imgur.com/xQUk24x.png) in the answer). — Hoki, Jan 09 '15 at 20:13

score 3 · Accepted Answer · edited Mar 31 '17 at 04:37

3

Ok so I've solved my own problem (thanks to the explanation on the MATLAB interp1 from Alex!).

The python interp1d doesn't have query points in itself, but instead creates a function which you then use to get your new data points. Thus, it should be:

        f = interpolate.interp1d( S[0], S[ii+1])
        temp.append(f(U[jj][ii,ii]))

edited Mar 31 '17 at 04:37

Anton K

4,658
2
47
60

answered Jan 09 '15 at 20:27

Steve

614
1
10
20

Yup that's exactly it. `numpy.interpolate.interp1d` gives you a **function** based on the keypoints you use. You would then use this function and interpolate your values. – rayryeng Jan 09 '15 at 20:37

score 1 · Answer 2 · answered Jan 09 '15 at 19:53

1

There is a python library that let's you use MATLAB functions through wrappers: mlabwrap. If you don't need to change the code of the functions itself this could save you some time.

answered Jan 09 '15 at 19:53

dudenr33

1,119
1
10
26

Thanks for the help, I'll be sure to look into that for future use! (I'm more or less at the end of this conversion but I'm sure I'll have more in the future!) – Steve Jan 09 '15 at 20:32

score 1 · Answer 3 · answered Jan 09 '15 at 20:08

I don't know scipy, but I can tell you what the interp1 call in MATLAB is doing:

http://www.mathworks.com/help/matlab/ref/interp1.html

You are using the syntax:

vq = interp1(x,v,xq,method)

"Vector x contains the sample points, and v contains the corresponding values, v(x). Vector xq contains the coordinates of the query points."

So, in your code, S(:,1) contains the sample points where your grid is defined, S(:,ii+1) contains your sampled values for your 1-D function, and U(ii,ii,:) contains the query points where you want to interpolate to find new functional values between known values in your grid. You are using linear interpolation.

1-D interpolation is an extremely well defined operation, and interp1 is a relatively straightforward interface for this operation. What exactly do you not understand? Are you clear what interpolation is?

Essentially, you have a discretely defined function f[x], the first argument to interp1 is x, the second argument is f[x], and the third argument are arbitrarily defined query points Xq at which you want to find new function values f[Xq]. Since these values are not known, you have to use an interpolation method for how you will approximate f[Xq]. 'linear' means you will use a linear weighted average of the two known sampled neighbors (left and right neighbors) nearest to Xq.

Ah, ok, it was the xq that had confused me on that front. Don't know why but "query points" in the help completely threw me. What you've said makes perfect sense. You're giving it a series of X points that correspond to V points, then it interpolates between those points using whatever method. Then it gives you the V points for any X points you provide (xq). Got it. Thanks for the help. No idea why I couldn't wrap my head around something so simple. I'll see if I can apply this understanding to the Python method. — Steve, Jan 09 '15 at 20:21

Converting MATLAB's interp1 to Python interp1d

3 Answers3