diff --git a/CompasToolkit.jl/examples/signal.jl b/CompasToolkit.jl/examples/signal.jl
index 9347925..65ac5dd 100644
--- a/CompasToolkit.jl/examples/signal.jl
+++ b/CompasToolkit.jl/examples/signal.jl
@@ -79,15 +79,14 @@ signal_ref = simulate(CUDALibs(), pssfp_ref, parameters_ref, trajectory_ref, coi
 signal_ref = collect(signal_ref)
 signal_ref = reshape(signal_ref, ns, nr)
 
-signal = zeros(ComplexF32, ns, nr, ncoils)
-CompasToolkit.magnetization_to_signal(
+signal = CompasToolkit.magnetization_to_signal(
     context,
-    signal,
     echos,
     parameters,
     trajectory,
     coil_sensitivities)
 
+signal = collect(signal)
 
 for c in 1:ncoils
     expected = map(x -> x[c], signal_ref)
diff --git a/CompasToolkit.jl/src/CompasToolkit.jl b/CompasToolkit.jl/src/CompasToolkit.jl
index 81f9873..a0206e3 100644
--- a/CompasToolkit.jl/src/CompasToolkit.jl
+++ b/CompasToolkit.jl/src/CompasToolkit.jl
@@ -34,11 +34,12 @@ make_context(device::Integer)::Context = Context(device)
 make_context() = make_context(0)
 
 mutable struct CompasArray{T, N} <: AbstractArray{T, N}
+    context::Context
     ptr::Ptr{Cvoid}
     sizes::Dims{N}
 
-    function CompasArray{T, N}(ptr::Ptr{Cvoid}, sizes::Dims{N}) where {T, N}
-        obj = new(ptr, sizes)
+    function CompasArray{T, N}(context::Context, ptr::Ptr{Cvoid}, sizes::Dims{N}) where {T, N}
+        obj = new(context, ptr, sizes)
         destroy = (obj) -> @ccall LIBRARY.compas_destroy_array(ptr::Ptr{Cvoid})::Cvoid
         finalizer(destroy, obj)
     end
@@ -50,23 +51,47 @@ Base.getindex(array::CompasArray, i) = throw(ArgumentError("cannot index into a
 function make_array(context::Context, input::Array{Float32, N})::CompasArray{Float32, N} where {N}
     sizes::Vector{Int64} = [reverse(size(input))...]
 
-    @ccall LIBRARY.compas_make_array_float(
+    ptr = @ccall LIBRARY.compas_make_array_float(
         context.ptr::Ptr{Cvoid},
         pointer(input)::Ptr{Float32},
         N::Int32,
         pointer(sizes)::Ptr{Int64}
+    )::Ptr{Cvoid}
+
+    return CompasArray{Float32, N}(context, ptr, Dims{N}(sizes))
+end
+
+function collect(input::CompasArray{ComplexF32, N}) where {N}
+    result = Array{ComplexF32, N}(undef, input.sizes...)
+    @ccall LIBRARY.compas_read_array_complex(
+        input.context::Ptr{Cvoid},
+        input.ptr::Ptr{Cvoid},
+        pointer(result)::Ptr{ComplexF32}
     )::Cvoid
+    return result
+end
+
+function collect(input::CompasArray{Float32, N}) where {N}
+    result = Array{Float32, N}(undef, input.sizes...)
+    @ccall LIBRARY.compas_read_array_float(
+        input.context::Ptr{Cvoid},
+        input.ptr::Ptr{Cvoid},
+        pointer(result)::Ptr{Float32}
+    )::Cvoid
+    return result
 end
 
 function make_array(context::Context, input::Array{ComplexF32, N})::CompasArray{ComplexF32, N} where {N}
     sizes::Vector{Int64} = [reverse(size(input))...]
 
-    @ccall LIBRARY.compas_make_array_complex(
+    ptr = @ccall LIBRARY.compas_make_array_complex(
         context.ptr::Ptr{Cvoid},
         pointer(input)::Ptr{ComplexF32},
         N::Int32,
         pointer(sizes)::Ptr{Int64}
-    )::Cvoid
+    )::Ptr{Cvoid}
+
+    return CompasArray{ComplexF32, N}(context, ptr, Dims{N}(sizes))
 end
 
 function assert_size(input::AbstractArray, expected::Dims{N}) where {N}
@@ -76,7 +101,7 @@ function assert_size(input::AbstractArray, expected::Dims{N}) where {N}
     end
 end
 
-function convert_array(
+function convert_array_old(
     ty::Type{T},
     dims::Dims{N},
     input::Array{T,N},
@@ -85,7 +110,7 @@ function convert_array(
     return input
 end
 
-function convert_array(
+function convert_array_old(
     ty::Type{T},
     dims::Dims{N},
     input::AbstractArray,
@@ -94,10 +119,20 @@ function convert_array(
     return convert(Array{T,N}, input)
 end
 
-function convert_array(ty::Type{T}, dims::Dims{N}, input::Number)::Array{T,N} where {T,N}
+function convert_array_old(ty::Type{T}, dims::Dims{N}, input::Number)::Array{T,N} where {T,N}
     return fill(convert(ty, input), dims)
 end
 
+function convert_array(
+    context::Context,
+    ty::Type{T},
+    dims::Dims{N},
+    input::Array{T,N},
+)::CompasArray{T,N} where {T,N}
+    assert_size(input, dims)
+    return make_array(context, input)
+end
+
 function unsafe_destroy_object!(obj)
     @ccall LIBRARY.compas_destroy(obj.ptr::Ptr{Cvoid})::Cvoid
 end
@@ -117,7 +152,7 @@ mutable struct CartesianTrajectory <: Trajectory
         k_start::AbstractVector,
         delta_k::Number,
     )
-        k_start = convert_array(ComplexF32, (nreadouts,), k_start)
+        k_start = convert_array_old(ComplexF32, (nreadouts,), k_start)
 
         ptr = @ccall LIBRARY.compas_make_cartesian_trajectory(
             pointer(context)::Ptr{Cvoid},
@@ -146,8 +181,8 @@ mutable struct SpiralTrajectory <: Trajectory
         k_start::AbstractVector,
         delta_k::AbstractVector,
     )
-        k_start = convert_array(ComplexF32, (nreadouts,), k_start)
-        delta_k = convert_array(ComplexF32, (nreadouts,), delta_k)
+        k_start = convert_array_old(ComplexF32, (nreadouts,), k_start)
+        delta_k = convert_array_old(ComplexF32, (nreadouts,), delta_k)
 
         ptr = @ccall LIBRARY.compas_make_spiral_trajectory(
             pointer(context)::Ptr{Cvoid},
@@ -181,15 +216,15 @@ mutable struct TissueParameters
         y::AbstractVector,
         z::AbstractVector,
     )
-        T1 = convert_array(Float32, (nvoxels,), T1)
-        T2 = convert_array(Float32, (nvoxels,), T2)
-        B1 = convert_array(Float32, (nvoxels,), B1)
-        B0 = convert_array(Float32, (nvoxels,), B0)
-        rho_x = convert_array(Float32, (nvoxels,), rho_x)
-        rho_y = convert_array(Float32, (nvoxels,), rho_y)
-        x = convert_array(Float32, (nvoxels,), x)
-        y = convert_array(Float32, (nvoxels,), y)
-        z = convert_array(Float32, (nvoxels,), z)
+        T1 = convert_array_old(Float32, (nvoxels,), T1)
+        T2 = convert_array_old(Float32, (nvoxels,), T2)
+        B1 = convert_array_old(Float32, (nvoxels,), B1)
+        B0 = convert_array_old(Float32, (nvoxels,), B0)
+        rho_x = convert_array_old(Float32, (nvoxels,), rho_x)
+        rho_y = convert_array_old(Float32, (nvoxels,), rho_y)
+        x = convert_array_old(Float32, (nvoxels,), x)
+        y = convert_array_old(Float32, (nvoxels,), y)
+        z = convert_array_old(Float32, (nvoxels,), z)
 
 
         ptr = @ccall LIBRARY.compas_make_tissue_parameters(
@@ -217,8 +252,12 @@ mutable struct TissueParameters
 end
 
 mutable struct FispSequence
-    ptr::Ptr{Cvoid}
-    nreadouts::Int32
+    RF_train::CompasArray{ComplexF32, 1}
+    slice_profiles::CompasArray{ComplexF32, 2}
+    TR::Float32
+    TE::Float32
+    max_state::Int32
+    TI::Float32
 
     function FispSequence(
         context::Context,
@@ -232,29 +271,22 @@ mutable struct FispSequence
         nreadouts = size(RF_train, 1)
         nslices = size(slice_profiles, 2)
 
-        RF_train = convert_array(ComplexF32, (nreadouts,), RF_train)
-        slice_profiles = convert_array(ComplexF32, (nreadouts, nslices), slice_profiles)
+        RF_train = convert_array_old(ComplexF32, (nreadouts,), RF_train)
+        slice_profiles = convert_array_old(ComplexF32, (nreadouts, nslices), slice_profiles)
 
-        ptr = @ccall LIBRARY.compas_make_fisp_sequence(
-            pointer(context)::Ptr{Cvoid},
-            nreadouts::Int32,
-            nslices::Int32,
-            pointer(RF_train)::Ptr{ComplexF32},
-            pointer(slice_profiles)::Ptr{ComplexF32},
-            TR::Float32,
-            TE::Float32,
-            max_state::Int32,
-            TI::Float32,
-        )::Ptr{Cvoid}
-
-        obj = new(ptr, nreadouts)
-        finalizer(unsafe_destroy_object!, obj)
+        return new(RF_train, slice_profiles, TR, TE, max_state, TI)
     end
 end
 
 mutable struct pSSFPSequence
-    ptr::Ptr{Cvoid}
-    nreadouts::Int32
+    RF_train::CompasArray{ComplexF32, 1}
+    TR::Float32
+    nRF::Int32
+    nTR::Int32
+    gamma_dt_GRz::CompasArray{ComplexF32, 1}
+    dt::NTuple{3, Float32}
+    gamma_dt_GRz::NTuple{3, Float32}
+    z::CompasArray{Float32, 1}
 
     function pSSFPSequence(
         context::Context,
@@ -269,29 +301,11 @@ mutable struct pSSFPSequence
         nRF = size(gamma_dt_RF, 1)
         nslices = size(z, 1)
 
-        RF_train = convert_array(ComplexF32, (nreadouts,), RF_train)
-        gamma_dt_RF = convert_array(ComplexF32, (nRF,), gamma_dt_RF)
-        z = convert_array(Float32, (nslices,), z)
+        RF_train = convert_array_old(ComplexF32, (nreadouts,), RF_train)
+        gamma_dt_RF = convert_array_old(ComplexF32, (nRF,), gamma_dt_RF)
+        z = convert_array_old(Float32, (nslices,), z)
 
-        ptr = @ccall LIBRARY.compas_make_pssfp_sequence(
-            pointer(context)::Ptr{Cvoid},
-            nRF::Int32,
-            nreadouts::Int32,
-            nslices::Int32,
-            pointer(RF_train)::Ptr{ComplexF32},
-            TR::Float32,
-            pointer(gamma_dt_RF)::Ptr{ComplexF32},
-            dt[1]::Float32,
-            dt[2]::Float32,
-            dt[3]::Float32,
-            gamma_dt_GRz[1]::Float32,
-            gamma_dt_GRz[2]::Float32,
-            gamma_dt_GRz[3]::Float32,
-            pointer(z)::Ptr{Float32},
-        )::Ptr{Cvoid}
-
-        obj = new(ptr, nreadouts)
-        finalizer(unsafe_destroy_object!, obj)
+        return new(RF_train, TR, nRF, nTR, gamma_dt_RF, dt, gamma_dt_GRz, z)
     end
 end
 
@@ -303,13 +317,18 @@ function simulate_magnetization(
 )
     nvoxels::Int64 = parameters.nvoxels
     nreadouts::Int64 = sequence.nreadouts
-    echos = convert_array(ComplexF32, (nvoxels, nreadouts), echos)
+    echos = convert_array_old(ComplexF32, (nvoxels, nreadouts), echos)
 
     @ccall LIBRARY.compas_simulate_magnetization_fisp(
         pointer(context)::Ptr{Cvoid},
         pointer(echos)::Ptr{ComplexF32},
         parameters.ptr::Ptr{Cvoid},
-        sequence.ptr::Ptr{Cvoid},
+        sequence.RF_train.ptr::Ptr{Cvoid},
+        sequence.slice_profiles.ptr::Ptr{Cvoid},
+        sequence.TR,
+        sequence.TE,
+        sequence.max_state,
+        sequence.TI
     )::Cvoid
 
     return echos
@@ -323,13 +342,21 @@ function simulate_magnetization(
 )
     nvoxels::Int64 = parameters.nvoxels
     nreadouts::Int64 = sequence.nreadouts
-    echos = convert_array(ComplexF32, (nvoxels, nreadouts), echos)
+    echos = convert_array_old(ComplexF32, (nvoxels, nreadouts), echos)
 
     @ccall LIBRARY.compas_simulate_magnetization_pssfp(
         pointer(context)::Ptr{Cvoid},
         pointer(echos)::Ptr{ComplexF32},
         parameters.ptr::Ptr{Cvoid},
-        sequence.ptr::Ptr{Cvoid},
+        sequence.RF_train.ptr::Ptr{Cvoid},
+        sequence.TR::Float32,
+        sequence.gamma_dt_RF.ptr::Ptr{Cvoid},
+        sequence.dt[0]::Float32,
+        sequence.dt[1]::Float32,
+        sequence.dt[2]::Float32,
+        sequence.gamma_dt_GRz[0]::Float32,
+        sequence.gamma_dt_GRz[1]::Float32,
+        sequence.gamma_dt_GRz[2]::Float32
     )::Cvoid
 
     return echos
@@ -337,7 +364,6 @@ end
 
 function magnetization_to_signal(
     context::Context,
-    signal::AbstractArray{<:Any,3},
     echos::AbstractMatrix,
     parameters::TissueParameters,
     trajectory::Trajectory,
@@ -348,21 +374,19 @@ function magnetization_to_signal(
     samples_per_readout::Int64 = trajectory.samples_per_readout
     nvoxels::Int64 = parameters.nvoxels
 
-    signal = convert_array(ComplexF32, (samples_per_readout, nreadouts, ncoils), signal)
-    echos = convert_array(ComplexF32, (nvoxels, nreadouts), echos)
-    coils = convert_array(Float32, (nvoxels, ncoils), coils)
+    echos = convert_array(context, ComplexF32, (nvoxels, nreadouts), echos)
+    coils = convert_array(context, Float32, (nvoxels, ncoils), coils)
 
-    @ccall LIBRARY.compas_magnetization_to_signal(
+    signal_ptr = @ccall LIBRARY.compas_magnetization_to_signal(
         pointer(context)::Ptr{Cvoid},
         ncoils::Int32,
-        pointer(signal)::Ptr{ComplexF32},
-        pointer(echos)::Ptr{ComplexF32},
+        pointer(echos)::Ptr{Cvoid},
         parameters.ptr::Ptr{Cvoid},
         trajectory.ptr::Ptr{Cvoid},
         pointer(coils)::Ptr{Float32}
-    )::Cvoid
+    )::Ptr{Cvoid}
 
-    return signal
+    return CompasArray{ComplexF32, 3}(context, signal_ptr, (ncoils,  nreadouts, samples_per_readout))
 end
 
 Base.pointer(c::Context) = c.ptr
diff --git a/julia-bindings/src/bindings.cpp b/julia-bindings/src/bindings.cpp
index ea63d3a..73f5fac 100644
--- a/julia-bindings/src/bindings.cpp
+++ b/julia-bindings/src/bindings.cpp
@@ -62,6 +62,15 @@ extern "C" const kmm::ArrayBase* compas_make_array_float(
     });
 }
 
+extern "C" void compas_read_array_float(
+    const compas::CudaContext* context,
+    const kmm::ArrayBase* input_array,
+    float* dest_ptr) {
+    catch_exceptions([&]() -> kmm::ArrayBase* {
+        input_array->block()->read(dest_ptr, input_array->size() * sizeof(float));
+    });
+}
+
 extern "C" const kmm::ArrayBase* compas_make_array_complex(
     const compas::CudaContext* context,
     const compas::cfloat* data_ptr,
@@ -82,6 +91,15 @@ extern "C" const kmm::ArrayBase* compas_make_array_complex(
     });
 }
 
+extern "C" void compas_read_array_complex(
+    const compas::CudaContext* context,
+    const kmm::ArrayBase* input_array,
+    compas::cfloat* dest_ptr) {
+    catch_exceptions([&]() -> kmm::ArrayBase* {
+        input_array->block()->read(dest_ptr, input_array->size() * sizeof(compas::cfloat));
+    });
+}
+
 extern "C" void compas_destroy_array(const kmm::ArrayBase* array) {
     return catch_exceptions([&] { delete array; });
 }
@@ -156,76 +174,38 @@ extern "C" const compas::TissueParameters* compas_make_tissue_parameters(
     });
 }
 
-extern "C" const compas::pSSFPSequence* compas_make_pssfp_sequence(
-    const compas::CudaContext* context,
-    int nRF,
-    int nreadouts,
-    int nslices,
-    const compas::cfloat* RF_train,
-    float TR,
-    const compas::cfloat* gamma_dt_RF,
-    float dt_ex,
-    float dt_inv,
-    float dt_pr,
-    float gamma_dt_GRz_ex,
-    float gamma_dt_GRz_inv,
-    float gamma_dt_GRz_pr,
-    const float* z) {
-    return catch_exceptions([&] {
-        auto seq = compas::make_pssfp_sequence(
-            *context,
-            make_view(RF_train, nreadouts),
-            TR,
-            make_view(gamma_dt_RF, nRF),
-            {dt_ex, dt_inv, dt_pr},
-            {gamma_dt_GRz_ex, gamma_dt_GRz_inv, gamma_dt_GRz_pr},
-            make_view(z, nslices));
-
-        return new compas::pSSFPSequence(seq);
-    });
-}
-
-extern "C" const compas::FISPSequence* compas_make_fisp_sequence(
-    const compas::CudaContext* context,
-    int nreadouts,
-    int nslices,
-    const compas::cfloat* RF_train,
-    const compas::cfloat* slice_profiles,
-    float TR,
-    float TE,
-    int max_state,
-    float TI) {
-    return catch_exceptions([&] {
-        auto seq = compas::make_fisp_sequence(
-            *context,
-            make_view(RF_train, nreadouts),
-            make_view(slice_profiles, nslices, nreadouts),
-            TR,
-            TE,
-            max_state,
-            TI);
-
-        return new compas::FISPSequence(seq);
-    });
-}
-
 extern "C" void compas_simulate_magnetization_fisp(
     const compas::CudaContext* context,
     compas::cfloat* echos_ptr,
     const compas::TissueParameters* parameters,
-    const compas::FISPSequence* sequence) {
+    compas::CudaArray<compas::cfloat>* RF_train,
+    compas::CudaArray<compas::cfloat, 2>* sliceprofiles,
+    float TR,
+    float TE,
+    int max_state,
+    float TI
+    ) {
     return catch_exceptions([&] {
-        int nreadouts = sequence->RF_train.size();
+        int nreadouts = RF_train->size();
         int nvoxels = parameters->nvoxels;
 
         auto echos = make_view(echos_ptr, nreadouts, nvoxels);
         auto d_echos = compas::CudaArray<compas::cfloat, 2>(nreadouts, nvoxels);
 
+        auto sequence = compas::FISPSequence {
+                *RF_train,
+                *sliceprofiles,
+                TR,
+                TE,
+                max_state,
+                TI
+        };
+
         context->submit_device(
             compas::simulate_magnetization_fisp,
             write(d_echos),
             *parameters,
-            *sequence);
+            sequence);
 
         d_echos.read(echos);
     });
@@ -235,29 +215,46 @@ extern "C" void compas_simulate_magnetization_pssfp(
     const compas::CudaContext* context,
     compas::cfloat* echos_ptr,
     const compas::TissueParameters* parameters,
-    const compas::pSSFPSequence* sequence) {
+    const compas::CudaArray<compas::cfloat>* RF_train,
+    float TR,
+    const compas::CudaArray<compas::cfloat>* gamma_dt_RF,
+    float dt_ex,
+    float dt_inv,
+    float dt_pr,
+    float gamma_dt_GRz_ex,
+    float gamma_dt_GRz_inv,
+    float gamma_dt_GRz_pr,
+    const compas::CudaArray<float>* z) {
     return catch_exceptions([&] {
-        int nreadouts = sequence->RF_train.size();
+        int nreadouts = RF_train->size();
         int nvoxels = parameters->nvoxels;
 
         auto echos = make_view(echos_ptr, nreadouts, nvoxels);
         auto d_echos = context->allocate(echos);
 
+        auto sequence = compas::pSSFPSequence {
+                *RF_train,
+                TR,
+                *gamma_dt_RF,
+                {dt_ex, dt_inv, dt_pr},
+                {gamma_dt_GRz_ex, gamma_dt_GRz_inv, gamma_dt_GRz_pr},
+                *z
+        };
+
         context->submit_device(
             compas::simulate_magnetization_pssfp,
             write(d_echos),
             *parameters,
-            *sequence);
+            sequence);
 
         d_echos.read(echos);
     });
 }
 
-extern "C" void compas_magnetization_to_signal(
+extern "C" kmm::ArrayBase* compas_magnetization_to_signal(
     const compas::CudaContext* context,
     int ncoils,
-    compas::cfloat* signal_ptr,
-    const compas::cfloat* echos_ptr,
+    const compas::CudaArray<compas::cfloat, 2>* echos,
     const compas::TissueParameters* parameters,
     const compas::Trajectory* trajectory,
     const float* coils_ptr) {
@@ -266,16 +263,13 @@ extern "C" void compas_magnetization_to_signal(
         int samples_per_readout = trajectory->samples_per_readout;
         int nvoxels = parameters->nvoxels;
 
-        auto signal = make_view(signal_ptr, ncoils, nreadouts, samples_per_readout);
-        auto echos = make_view(echos_ptr, nreadouts, nvoxels);
         auto coils = make_view(coils_ptr, ncoils, nvoxels);
 
-        auto d_echos = context->allocate(echos);
         auto d_coils = context->allocate(coils);
 
-        auto d_signal =
-            compas::magnetization_to_signal(*context, d_echos, *parameters, *trajectory, d_coils);
+        auto signal =
+            compas::magnetization_to_signal(*context, *echos, *parameters, *trajectory, d_coils);
 
-        d_signal.read(signal);
+        return new compas::CudaArray<compas::cfloat, 3>(signal);
     });
 }
\ No newline at end of file