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
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.
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
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}:
...
Hello,
I think now works. Thank you very much.
Kind regards,
/Martin