# Lookup table interpolation of an array

Hello,

I am not yet fully familiar with Julia and am currently trying to rewrite a Matlab code. I am stumbling with the (Matlab) function interp1. Summarized the code looks like this:

A → Type: double(2048x1)
B → Type: complex-single(2048,512)
C → Type: double(1x2048)

D = interp1(A,B,C,“spline”)

Can anyone give me a tip on how I can write this in Julia, previous attempts have only led to wild error messages. Many thanks in advance.

Kind regards,

Martin

1 Like

Have you tried Home · Interpolations.jl?

2 Likes

I would also look at GitHub - jipolanco/BSplineKit.jl: A collection of B-spline tools in Julia for irregular data.

I also list several other packages here:
https://juliamath.github.io/Interpolations.jl/latest/other_packages/#Other-Interpolation-Packages

Also post your wild error messages.

3 Likes

Thanks for the answer, I tried "interpolate* from Interpolations.jl as well as from BSplineKit.jl. The code I use is

`using BSplineKit;`

with

A → typeof : Matrix{ComplexF64}, size : (2048, 1)
B → typeof : Matrix{ComplexF64}, size : (2048, 512)
C → typeof : Vector{Int64}, size : (2048,)

and in case of

`itp = interpolate(A[:,1], B[:,1], BSplineOrder(1))`

I received

`MethodError: no method matching make_knots(::Vector{ComplexF64}, ::BSplineOrder{1}, ::Nothing)`

In case of

`itp = interpolate(A, B, BSplineOrder(1))`

I received

`MethodError: no method matching interpolate(::Matrix{ComplexF64}, ::Matrix{ComplexF64}, ::BSplineOrder{1})`

I also tried

`spl = Spline1D(A, B , C , k=3, bc="nearest")`

which results in

`MethodError: no method matching Spline1D(::Matrix{ComplexF64}, ::Matrix{ComplexF64}, ::Vector{Int64}; k::Int64, bc::String)`

So, I’m not sure what is the best approach to replace the interp1 function from MatLab.

Kind regards,

Martin

Cannot test Matlab directly, but from what I understand 1D interpolation needs real inputs, i.e., would not work on complex inputs `A`.
Here is an example – assuming real inputs – using `BSplineKit`:

``````using Distributions, BSplineKit
N = 2048
A = sort(rand(Uniform(0, 10), N))  # Note: Real input vector
B1 = rand(ComplexF64, N)  # Single complex output vector
sp1 = BSplineKit.interpolate(A, B1, BSplineKit.BSplineOrder(3))  # Get interpolation function
C = rand(Uniform(0, 10), N)
sp1.(C)  # Note: Broadcast spline function over desired new inputs

# In contrast to Matlab, to handle multiple outputs, we need to be more explicit and broadcast manually
B = rand(ComplexF64, N, 512)  # Matrix of multiple outputs
sps = BSplineKit.interpolate.(Ref(A), eachcol(B), Ref(BSplineKit.BSplineOrder(3)))  # Get vector of splines interpolating each column of B
# Finally, apply each of these splines to all input values and combine results
sps .|> (f -> f.(C)) |> stack
``````
2 Likes

Thank you for the example, it works well for real data; nevertheless I need to interpolate complex data …

He demonstrated with complex data. I’m unclear what is missing.

By the way, the above methods work more similarly to MATLAB’s gridded interpolant:

That is first we input the sample points and the values at those sample points and get an object that we can call to obtain other points.

Let’s breakdown what @bertschi demonstrated above. We first create a real vector of sample points in `A` and a complex 32-bit floating point in `B1`. These are both column vectors.

``````julia> using Distributions, BSplineKit

julia> N = 2048
2048

julia> A = sort(rand(Uniform(0, 10), N));

julia> B1 = rand(ComplexF32, N)
2048-element Vector{ComplexF32}:
0.8409336f0 + 0.19596744f0im
0.9809899f0 + 0.5659446f0im
0.79642314f0 + 0.0299052f0im
0.5722098f0 + 0.89307314f0im
0.96031356f0 + 0.1327033f0im
0.71136796f0 + 0.862921f0im
0.5338061f0 + 0.76945704f0im
0.5737985f0 + 0.1366781f0im
0.27081704f0 + 0.2231226f0im
0.25618458f0 + 0.18294215f0im
⋮
0.17688131f0 + 0.44664085f0im
0.7323388f0 + 0.94128615f0im
0.35286158f0 + 0.29547954f0im
0.8608282f0 + 0.7629858f0im
0.047004223f0 + 0.70044565f0im
0.8486275f0 + 0.44268453f0im
0.8413368f0 + 0.25111085f0im
0.23773956f0 + 0.5144081f0im
0.9474749f0 + 0.6850947f0im

julia> typeof(A)
Vector{Float64} (alias for Array{Float64, 1})

julia> typeof(B1)
Vector{ComplexF32} (alias for Array{Complex{Float32}, 1})
``````

Then we create the `SplineInterpolation object`. Here I use `BSplineOrder(4)` for a cubic spline, which is what MATLAB uses.

``````julia> sp1 = BSplineKit.interpolate(A, B1, BSplineKit.BSplineOrder(4))
SplineInterpolation containing the 2048-element Spline{ComplexF32}:
basis: 2048-element BSplineBasis of order 4, domain [0.000857447419672086, 9.992905397874818]
order: 4
knots: [0.000857447, 0.000857447, 0.000857447, 0.000857447, 0.00779442, 0.0123642, 0.0176373, 0.0214538, 0.0264004, 0.0297231  …  9.9547, 9.95925, 9.96464, 9.96861, 9.97387, 9.98188, 9.99291, 9.99291, 9.99291, 9.99291]
coefficients: ComplexF32[0.840934+0.195967im, 1.13049+1.21022im, 0.780093-1.07272im, 0.235346+1.87558im, 1.19961-0.603868im, 0.938954+1.03391im, 0.0414462+0.989143im, 2.25266-1.20595im, -0.204628+0.862867im, 0.804739-1.07816im  …  2.00787+2.41364im, -2.50337-2.01189im, 5.14302+4.93931im, -1.16833-0.828984im, 1.90401+0.996272im, -0.889589+0.651825im, 2.48836+0.234553im, 1.07138+0.456499im, -4.38962-0.590771im, 0.947475+0.685095im]
interpolation points: [0.000857447, 0.00305412, 0.00675145, 0.0135777, 0.0167635, 0.0225707, 0.0250273, 0.0316031, 0.032539, 0.0327446  …  9.94789, 9.95161, 9.95218, 9.9603, 9.96527, 9.96834, 9.97223, 9.98104, 9.99238, 9.99291]
``````

We then query the `SplineInterpolation` object at various points via column vector `C` which is also randomly generated.

``````julia> C = rand(Uniform(0,10), N)
2048-element Vector{Float64}:
9.483733019514059
5.43085430896525
8.318840782645099
5.288133871920023
4.0813642928007
6.78805933501573
6.291783894803595
5.170340091879966
5.295940125355911
6.382796375189246
⋮
8.079857717607496
4.940743863230217
0.8397228433210491
8.198655558986061
5.212823437408405
9.677186608905854
1.7191592119233634
9.625444117783816
5.039112932203945
``````

We can then query the spline as follows using the dot syntax.

``````julia> sp1.(C)
2048-element Vector{ComplexF64}:
0.6628213223809246 + 0.13845306052661216im
7.266404546511955 + 1.4918834901186397im
2.529257327289871 - 3.0934900229509243im
2.082829065843516 + 0.35194745428180896im
0.42051821379475274 + 0.11991317115395697im
-0.06217035130524348 + 0.42779375900348354im
-1.29880373722419 + 0.410014463695485im
0.7360561045406451 + 1.4965432408578285im
0.29754144965019985 + 0.6998459950323489im
0.8115314971604752 + 0.7237799271783811im
⋮
3.3251618089535615 - 5.701583525945178im
-21.331010022778777 + 73.7176687071869im
1.2751402913844252 + 1.268245785368036im
-0.06294013906448183 + 1.631632769479881im
0.766038339089081 + 1.10628834593872im
0.20352877087166413 + 0.6675186029274851im
-0.2737631672819028 + 0.2230333372408741im
0.41110098498380204 + 0.1973548701882316im
18.76530445342889 + 22.771700462031767im

julia> typeof(ans)
Vector{ComplexF64} (alias for Array{Complex{Float64}, 1})
``````

We get a complex valued answer back.

For your original question, you provided multiple columns, 512 of them. This corresponds to 512 independent splines. bertschi’s second example demonstrated how to do this, but I will spell it out in a simple to use function for you that I will save in a file called `interp1.jl`.

``````import BSplineKit

function interp1(
x::AbstractVector,
v::AbstractVector,
xq::AbstractVector,
method::String = "spline"
)
spline = interp1(x, v, method)
return spline.(xq)
end

function interp1(
x::AbstractVector,
v::AbstractArray,
xq::AbstractArray,
method::String = "spline"
)
splines = interp1(x, v, method)
results = map(splines) do spline
spline.(xq)
end
return stack(results)
end

function interp1(
x::AbstractVector,
v::AbstractVector,
method::String = "spline"
)
order = if method == "spline" || method == "cubic"
BSplineKit.BSplineOrder(4)
elseif method == "quadratic"
BSplineKit.BSplineOrder(3)
elseif method == "linear"
BSplineKit.BSplineOrder(2)
elseif method == "constant"
BSplineKit.BSplineOrder(1)
end

return BSplineKit.interpolate(x, v, order)
end

function interp1(
x::AbstractVector,
v::AbstractArray,
method::String = "spline"
)
v_matrix = reshape(v, size(v, 1), prod(size(v)[2:end]))
splines = map(eachcol(v_matrix)) do col
interp1(x, col, method)
end
return splines
end

function interp1(
x::AbstractVector,
v::AbstractArray,
method::String = "spline"
)
v_matrix = reshape(v, size(v, 1), prod(size(v)[2:end]))
splines = map(eachcol(v_matrix)) do col
interp1(x, col, method)
end
return reshape(splines, size(v)[2:end])
end

function interp1(
x::AbstractArray,
v::AbstractArray,
method::String = "spline"
)
if !all(size(x) .== size(v))
error("The size of x and v must the be same. Got \$(size(x)) and \$(size(v))")
end
x_matrix = reshape(v, size(x, 1), prod(size(x)[2:end]))
v_matrix = reshape(v, size(v, 1), prod(size(v)[2:end]))
splines = map(eachcol(x_matrix), eachcol(v_matrix)) do x_col, v_col
interp1(x_col, v_col, method)
end
return reshape(splines, size(v)[2:end])
end
``````

Now we can use this as follows.

``````julia> using BSplineKit, Distributions

julia> include("interp1.jl")
interp1 (generic function with 10 methods)

julia> N = 2048; M = 512;

julia> A = sort(rand(Uniform(0,10), N)); size(A)
(2048,)

julia> B = rand(ComplexF32, N, M); size(B)
(2048, 512)

julia> C = rand(Uniform(0,10), N); size(C)
(2048,)

julia> D = interp1(A, B, C, "spline")
2048×512 Matrix{ComplexF64}:
-3.56976+0.601276im   -0.219194-0.447994im    …    -1.13067+0.473447im
-6.33133-0.066161im    0.965402+2.17368im           0.11428+7.66711im
0.265947+0.564525im   0.0940509+0.639697im        0.0618999+0.717675im
-0.525911-0.268237im    0.030246+0.119171im         0.433661+0.519061im
0.5511+0.237039im    0.589046+0.306925im       -0.0207684+1.00121im
0.601093+0.613032im    0.681022+0.923576im    …    0.622523+0.701978im
-0.310867+2.16622im     0.057356-0.517554im          2.94041+1.27067im
0.513746+0.668517im   -0.145421+0.247885im         0.966718-0.471837im
1.90147+0.734572im      0.1871+1.62462im           3.60835-1.48335im
0.747279+0.0694702im    1.03762+0.918896im         0.798624+0.0975493im
⋮                                     ⋱
-1.75224-0.431601im     1.33743-2.43919im           0.67119+2.0219im
0.26693+0.842739im   0.0213721+0.180948im    …    0.147373+0.755189im
2.36164+0.665174im    0.710751+0.182283im          1.62779+1.1795im
4.38767+3.99549im     -2.68346-0.922896im         -1.24872+0.693219im
0.691161+0.844551im    0.664908+0.00167152im       0.635339+0.898001im
0.552734+0.512182im     0.31708+0.480791im         0.807158+0.48275im
0.570123+0.512312im    0.215266+0.285696im    …    0.640996+0.0682268im
0.322452-0.137867im   -0.147389+0.277156im       0.00488787+0.61134im
-0.204593+0.740612im     0.89808-0.270704im          1.12213-0.0284527im
``````

Also you gave `C` as a row vector for some reason. This will also work. What I’m not entirely sure about is the expected size of the output `D` for your input.

``````julia> C2 = transpose(C)
1×2048 transpose(::Vector{Float64}) with eltype Float64:
9.15085  0.0196092  5.6302  3.97698  …  9.05284  8.40113  5.87534  7.91804

julia> D2 = interp1(A, B, C2, "spline")
1×2048×512 Array{ComplexF64, 3}:
...
``````
1 Like

Hello,

I think now works. Thank you very much.

Kind regards,
/Martin