diff --git a/CHANGELOG.md b/CHANGELOG.md
index 832fb18e14346886f769e1ee668a46e51d7b4324..02239c641fa6ea6837ff3e9d21b89f56466dbb4b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,7 +7,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 - Finite RF pulses (`RFPulse`) to more realistically model the effect of these pulses.
 - Hierchical properties: MRI relaxation parameters and collision parameters (e.g., permeability) can now set at both the global level and for each obstruction.
-- `run_tests.sh` function which will run the tests and produce a coverage that can be visualised in VS Code by the Coverage Gutters plugin
+- Sequences can now be read from the pulseq format (http://pulseq.github.io/) using `read_pulseq`.
+- `run_tests.sh` function which will run the tests and produce a coverage that can be visualised in VS Code by the Coverage Gutters plugin.
 ### Changed
 - Made units consistent with angles (i.e., phases and flip angles) in degrees, off-resonance fields in kHz, and gradients in kHz/um.
 - `time` function has been renamed `get_time`, so as not to conflict with `Base.time`
diff --git a/Project.toml b/Project.toml
index e59e3e72084660b4a9e5b42ef2b16a5d4f97b1bb..7f4456b0df32f7682e6fc242bd501a56a0855ef6 100644
--- a/Project.toml
+++ b/Project.toml
@@ -7,6 +7,7 @@ version = "0.4.0"
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
diff --git a/docs/src/sequence.md b/docs/src/sequence.md
index f29795f7add334c0be63cd5f12939d2d846074af..eaf80525e9e751910ece3991bf9e51568abcad12 100644
--- a/docs/src/sequence.md
+++ b/docs/src/sequence.md
@@ -1,4 +1,8 @@
 # [MR sequences](@id sequence)
+## Pulseseq sequences
+The recommended way to define sequences is to use pulseq sequence format (http://pulseq.github.io/).
+Pulseq sequences can be generated using matlab (http://pulseq.github.io/) or python (https://pypulseq.readthedocs.io/en/master/).
+They can be loaded into MCMRSimulator using [`read_pulseq`](@ref).
 ## Built-in MR sequences
 ### Diffusion-weighted MRI
 A pulsed-gradient spin-echo can be created using [`dwi`](@ref)
diff --git a/src/MCMRSimulator.jl b/src/MCMRSimulator.jl
index bd582d8317df73cef446a09ec42e03e20b5a6718..b62d53e05d7b98f9da251cdd008e55017f716eff 100644
--- a/src/MCMRSimulator.jl
+++ b/src/MCMRSimulator.jl
@@ -25,6 +25,7 @@ import Accessors: @set
 import Roots
 import PlyIO
 import Optim
+import DataStructures: OrderedDict
 
 
 include("constants.jl")
@@ -40,7 +41,7 @@ include("plot/plot.jl")
 
 export Sequence, InstantRFPulse, Readout, InstantGradient
 export MRGradients, gradient, rotate_bvec
-export dwi
+export dwi, read_pulseseq
 export Scanner, B0, max_gradient, max_slew_rate
 export Siemens_Connectom, Siemens_Prisma, Siemens_Terra
 export longitudinal, transverse, phase, vector, Spin, Snapshot, SpinOrientation, isinside, off_resonance
diff --git a/src/sequence/gradients.jl b/src/sequence/gradients.jl
index 91e92606a6aa84f91d59f4dcb45798a2740c7851..96ac028dc210782bd8b706a8eb7c417a64caccb0 100644
--- a/src/sequence/gradients.jl
+++ b/src/sequence/gradients.jl
@@ -49,11 +49,11 @@ function MRGradients(times::AbstractVector{<:Number}, amplitudes::AbstractVector
 end
 
 add_TR(g::MRGradients, TR::Number) = MRGradients(add_TR(g.Gx, Tr),add_TR(g.Gy, TR), add_TR(g.Gz, TR), g.origin)
-control_points(g::MRGradients) = [
+control_points(g::MRGradients) = sort(unique([
     control_points(g.Gx)...,
     control_points(g.Gy)...,
     control_points(g.Gz)...,
-]
+]))
 
 
 """
diff --git a/src/sequence/pulseseq.jl b/src/sequence/pulseseq.jl
new file mode 100644
index 0000000000000000000000000000000000000000..c389f64ee1d05f43bca57f039e05aa94106a1d4c
--- /dev/null
+++ b/src/sequence/pulseseq.jl
@@ -0,0 +1,273 @@
+struct PulseSeqSection
+    title :: String
+    content :: Vector{String}
+end
+
+struct CompressedPulseSeqShape
+    id :: Int
+    num :: Int
+    samples :: Vector{Float}
+end
+
+
+"""
+    Shape(shape::CompressedPulseSeqShape, step_size)
+
+Create a `Shape` object based on the pulseseq Shape format presuming each sample is `step_size` milliseconds apart.
+"""
+function Shape(pulseseq::CompressedPulseSeqShape)
+    compressed = length(pulseseq.samples) != pulseseq.num
+    if !compressed
+        times = range(0, 1, length=pulseseq.num)
+        return Shape(times, amplitudes)
+    end
+    times = [zero(Float)]
+    amplitudes = [pulseseq.samples[1]]
+    repeating = false
+    time_norm = 1 / (pulseseq.num - 1)
+    prev_sample = pulseseq.samples[1]
+    prev_applied = true
+    for sample in pulseseq.samples[2:end]
+        if repeating
+            nrepeats = Int(sample) + 2
+            if prev_applied
+                nrepeats -= 1
+            end
+            push!(times, times[end] + nrepeats * time_norm)
+            push!(amplitudes, amplitudes[end] + nrepeats * prev_sample)
+            repeating = false
+            prev_applied = true
+        elseif sample == prev_sample
+            repeating = true
+        else
+            if !prev_applied
+                push!(times, times[end] + time_norm)
+                push!(amplitudes, amplitudes[end] + prev_sample)
+            end
+            prev_sample = sample
+            prev_applied = false
+        end
+    end
+    if !prev_applied
+        push!(times, times[end] + time_norm)
+        push!(amplitudes, amplitudes[end] + prev_sample)
+    end
+    return Shape(times, amplitudes)
+end
+
+"""
+    read_pulseseq(filename; scanner=Scanner(B0=B0), B0=3., TR=<sequence duration>)
+
+Reads a sequence from a pulseq file (http://pulseq.github.io/).
+Pulseq files can be produced using matlab (http://pulseq.github.io/) or python (https://pypulseq.readthedocs.io/en/master/).
+"""
+function read_pulseseq(filename; kwargs...)
+    keywords = open(read_pulseseq_sections, filename)
+    build_sequence(; kwargs..., keywords...)
+end
+
+function read_pulseseq_sections(io::IO)
+    sections = PulseSeqSection[]
+    version = nothing
+    title = ""
+    for line in readlines(io)
+        line = strip(line)
+        if length(line) == 0 || line[1] == '#'
+            continue  # ignore comments
+        end
+        if line[1] == '[' && line[end] == ']'
+            if title == "VERSION"
+                version = parse_pulseq_section(sections[end]).second
+            end
+            # new section starts
+            title = line[2:end-1]
+            push!(sections, PulseSeqSection(title, String[]))
+        elseif length(sections) > 0
+            push!(sections[end].content, line)
+        else
+            error("Content found in pulseseq file before first section")
+        end
+    end
+    Dict(parse_pulseq_section.(sections, version)...)
+end
+
+function parse_pulseq_ordered_dict(strings::Vector{<:AbstractString}, names, dtypes)
+    as_dict = OrderedDict{Int, NamedTuple{Tuple(names), Tuple{dtypes...}}}()
+    for line in strings
+        parts = split(line)
+        @assert length(parts) == length(names)
+        values = parse.(dtypes, split(line))
+        @assert names[1] == :id
+        as_dict[values[1]] = (; zip(names, values)...)
+    end
+    return as_dict
+end
+
+function parse_pulseq_properties(strings::Vector{<:AbstractString})
+    result = Dict{String, Any}()
+    for s in strings
+        (name, value) = split(s, limit=2)
+        result[name] = value
+    end
+    return result
+end
+
+section_headers = Dict(
+    "BLOCKS" => ([:id, :duration, :rf, :gx, :gy, :gz, :adc, :ext], fill(Int, 8)),
+    ("BLOCKS", v"1.3.1") => ([:id, :delay, :rf, :gx, :gy, :gz, :adc, :ext], fill(Int, 8)),
+    ("DELAYS", v"1.3.1") => ([:id, :delay], [Int, Int]),
+    "RF" => ([:id, :amp, :mag_id, :phase_id, :delay, :freq, :phase], [Int, Float, Int, Int, Int, Float, Float]),
+    "Gradients" => ([:id, :amp, :shape_id, :delay], [Int, Float, Int, Int]),
+    "TRAP" => ([:id, :amp, :rise, :flat, :fall, :delay], [Int, Float, Int, Int, Int, Int]),
+    "ADC" => ([:id, :num, :dwell, :delay, :freq, :phase], [Int, Int, Float, Int, Float, Float]),
+)
+
+function parse_pulseq_section(section::PulseSeqSection, version=nothing)
+    if section.title == "VERSION"
+        props = parse_pulseq_properties(section.content)
+        result = VersionNumber(
+            parse(Int, props["major"]),
+            parse(Int, props["minor"]),
+            parse(Int, props["revision"]),
+        )
+    elseif section.title == "DEFINITIONS"
+        result = parse_pulseq_properties(section.content)
+    elseif (section.title, version) in keys(section_headers)
+        result = parse_pulseq_ordered_dict(section.content, section_headers[(section.title, version)]...)
+    elseif section.title in keys(section_headers)
+        result = parse_pulseq_ordered_dict(section.content, section_headers[section.title]...)
+    elseif section.title == "EXTENSION"
+        println("Ignoring all extensions in pulseseq")
+    elseif section.title == "SHAPES"
+        current_id = -1
+        shapes = CompressedPulseSeqShape[]
+        for line in section.content
+            if length(line) > 8 && lowercase(line[1:8]) == "shape_id"
+                current_id = parse(Int, line[9:end])
+                continue
+            end
+            for text in ("num_uncompressed", "num_samples")
+                if startswith(lowercase(line), text)
+                    @assert current_id != -1
+                    push!(shapes, CompressedPulseSeqShape(current_id, parse(Int, line[length(text)+1:end]), Float[]))
+                    current_id = -1
+                    break
+                end
+            end
+            if !startswith(lowercase(line), "num")
+                @assert current_id == -1
+                push!(shapes[end].samples, parse(Float, line))
+            end
+        end
+        result = Dict(
+            [s.id => (s.num, Shape(s)) for s in shapes]...
+        )
+    elseif section.title in ["SIGNATURE"]
+        # silently ignore these sections
+    else
+        error("Unrecognised pulseq section: $(section.title)")
+    end
+    return Symbol(lowercase(section.title)) => result
+end
+
+
+function build_sequence(; scanner=nothing, B0=3., TR=nothing, definitions, version, blocks, rf=nothing, gradients=nothing, trap=nothing, delays=nothing, shapes=nothing, adc=nothing)
+    if isnothing(scanner)
+        scanner = Scanner(B0=B0)
+    end
+    if version == v"1.4.0"
+        # load raster times (converting seconds to milliseconds)
+        convert = key -> parse(Float, definitions[key]) * Float(1e3)
+        gradient_raster = convert("GradientRasterTime")
+        rf_raster = convert("RadiofrequencyRasterTime")
+        adc_raster = convert("AdcRasterTime")
+        block_duration_raster = convert("BlockDurationRasterTime")
+    elseif version == v"1.3.1"
+        gradient_raster = rf_raster = block_duration_raster = Float(1e-3) # 1 microsecond as default raster
+        adc_raster = Float(1e-6) # ADC dwell time is in ns by default
+    else
+        error("Can only load pulseseq files with versions v1.3.1 and v1.4.0, not $(version)")
+    end
+
+    gradient_control_points = (
+        # times and amplitudes of gradient control points
+        ([zero(Float)], [zero(Float)]),
+        ([zero(Float)], [zero(Float)]),
+        ([zero(Float)], [zero(Float)]),
+    )
+    
+    block_start_time = zero(Float)
+    rf_pulses = RFPulse[]
+    readouts = Readout[]
+    for block in values(blocks)
+        block_duration = zero(Float)
+        if !iszero(block.rf)
+            proc = rf[block.rf]
+            start_time = block_start_time + proc.delay * 1e-3
+            @assert version == v"1.3.1" || iszero(proc.time_id)
+            (num, mag_shape) = shapes[proc.mag_id]
+            block_duration = max(num * rf_raster + proc.delay * 1e-3, block_duration)
+            mag = ConcreteShape(start_time, start_time + num * rf_raster, proc.amp * 1e-3, 0, mag_shape)
+            if iszero(proc.phase_id)
+                phase = ConcreteShape([start_time, start_time + num * rf_raster], [0, 0])
+            else
+                (num, phase_shape) = shapes[proc.phase_id]
+                block_duration = max(num * rf_raster + proc.delay * 1e-3, block_duration)
+                phase = ConcreteShape(start_time, start_time + num * rf_raster, proc.freq * 1e-3 * 360, rad2deg(proc.phase), phase_shape)
+            end
+            push!(rf_pulses, RFPulse(mag, phase, sqrt(proc.amp^2 + proc.freq^2)))
+        end
+        if !iszero(block.adc)
+            proc = adc[block.adc]
+            mid_readout = block_start_time + proc.delay * 1e-3 + proc.dwell * proc.num * adc_raster / 2
+            push!(readouts, Readout(mid_readout))
+            block_duration = max(proc.delay * 1e-3 + proc.dwell * proc.num * adc_raster, block_duration)
+        end
+        for ((times_grad, amps_grad), symbol_grad) in zip(gradient_control_points, [:gx, :gy, :gz])
+            grad_id = getfield(block, symbol_grad)
+            if iszero(grad_id)
+                continue
+            end
+            if grad_id in keys(gradients)
+                proc = gradients[grad_id]
+                @assert version == v"1.3.1" || iszero(proc.time_id)
+                start_time = block_start_time + proc.delay * 1e-3
+                (num, grad_shape) = shapes[proc.shape_id]
+                append!(times_grad, grad_shape.times .* (num * gradient_raster) .+ start_time)
+                append!(amps_grad, grad_shape.amplitudes .* (proc.amp * 1e3))
+                block_duration = max(proc.delay * 1e-3 + num * gradient_raster, block_duration)
+            elseif grad_id in keys(trap)
+                proc = trap[grad_id]
+                @assert version == v"1.3.1" || iszero(proc.time_id)
+                start_time = block_start_time + proc.delay * 1e-3
+                times = (cumsum([0, proc.rise, proc.flat, proc.fall]) .* 1e-3) .+ start_time
+                amplitudes = [0, proc.amp * 1e3, proc.amp * 1e3, 0]
+                append!(times_grad, times)
+                append!(amps_grad, amplitudes)
+                block_duration = max((proc.delay + proc.rise + proc.flat + proc.fall) * 1e-3, block_duration)
+            else
+                error("Gradient ID $grad_id not found in either of [GRADIENTS] or [TRAP] sections")
+            end
+        end
+        if version == v"1.3.1"
+            if !iszero(block.delay)
+                duration = max(block_duration, delays[block.delay].delay * 1e-3)
+            else
+                duration = block_duration
+            end
+        else
+            duration = block.duration * block_duration_raster
+            @assert duration >= block_duration
+        end
+        block_start_time += duration
+    end
+    if isnothing(TR)
+        TR = get(definitions, "TotalDuration", block_start_time)
+    end
+    for (times, amplitudes) in gradient_control_points
+        push!(times, TR)
+        push!(amplitudes, 0)
+    end
+    gradients = MRGradients([ConcreteShape(t, a) for (t, a) in gradient_control_points]..., zero(PosVector))
+    Sequence(; scanner=scanner, TR=TR, pulses=[rf_pulses..., readouts...], gradients=gradients)
+end
\ No newline at end of file
diff --git a/src/sequence/radio_frequency.jl b/src/sequence/radio_frequency.jl
index 3cd879cc0ca43feef8f48e6eda696d16b0fb5f7f..064561134ce7fa7c3c58425915d8fa672c25c005 100644
--- a/src/sequence/radio_frequency.jl
+++ b/src/sequence/radio_frequency.jl
@@ -31,7 +31,7 @@ function RFPulse(times::AbstractVector, amplitudes::AbstractVector, phases=nothi
 end
 
 start_time(pulse::RFPulse) = min(start_time(pulse.amplitude), start_time(pulse.phase))
-end_time(pulse::RFPulse) = min(end_time(pulse.amplitude), end_time(pulse.phase))
+end_time(pulse::RFPulse) = max(end_time(pulse.amplitude), end_time(pulse.phase))
 
 """
     amplitude(rf_pulse, t1[, t2])
diff --git a/src/sequence/sequence.jl b/src/sequence/sequence.jl
index 6fc7a715a16b71ea4f44f9e917c99eff934f8d3a..13a1849ddf60fd2f336d5e3ffeb651fb5d8dc853 100644
--- a/src/sequence/sequence.jl
+++ b/src/sequence/sequence.jl
@@ -4,7 +4,7 @@ include("gradients.jl")
 include("radio_frequency.jl")
 
 """
-    Sequence(;TR, gradients=nothing, pulses=nothing, B0=3., interplate_gradients=:step)
+    Sequence(;TR, gradients=nothing, pulses=nothing, scanner=Scanner(B0), B0=3., interplate_gradients=:step)
 
 An MR sequence represented by a series of pulses repeated with a given repetition time (`TR`).
 
@@ -135,4 +135,5 @@ gradient(seq::Sequence, time::Number) = gradient(seq.gradient, mod(time, seq.TR)
 gradient(seq::Sequence, time1::Number, time2::Number) = gradient(seq.gradient, mod(time1, seq.TR), iszero(mod(time2, seq.TR)) ? seq.TR : mod(time2, seq.TR))
 
 include("timestep.jl")
-include("diffusion.jl")
\ No newline at end of file
+include("diffusion.jl")
+include("pulseseq.jl")
\ No newline at end of file
diff --git a/src/sequence/shape.jl b/src/sequence/shape.jl
index bbd77702776f32151db759f8ecd32b7fb2107dfe..15003cd593430a6e55080081a38a42fa3ae94cac 100644
--- a/src/sequence/shape.jl
+++ b/src/sequence/shape.jl
@@ -135,8 +135,9 @@ struct ConcreteShape{N}
     t0::Float
     t1::Float
     max_amplitude::Float
+    shift::Float
     shape::Shape{N}
-    ConcreteShape(t0::Number, t1::Number, max_amplitude::Number, shape::Shape{N}) where {N} = new{N}(Float(t0), Float(t1), Float(max_amplitude), shape)
+    ConcreteShape(t0::Number, t1::Number, max_amplitude::Number, shift::Number, shape::Shape{N}) where {N} = new{N}(Float(t0), Float(t1), Float(max_amplitude), Float(shift), shape)
 end
 
 """
@@ -156,7 +157,7 @@ function ConcreteShape(times::AbstractVector{<:Number}, amplitudes::AbstractVect
     t1 = times[end]
     max_amplitude = maximum(abs.(amplitudes))
     ConcreteShape(
-        t0, t1, max_amplitude,
+        t0, t1, max_amplitude, 0,
         Shape(
             (times .- t0) ./ (t1 - t0),
             iszero(max_amplitude) ? amplitudes : (amplitudes ./ max_amplitude)
@@ -172,14 +173,14 @@ Produces a ConcreteShape with only zero amplitude.
 ConcreteShape() = ConcreteShape(Float[], Float[])
 
 convert_time(concrete::ConcreteShape, time::Number) = (time - concrete.t0) / (concrete.t1 - concrete.t0)
-amplitude(concrete::ConcreteShape, time::Number) = amplitude(concrete.shape, convert_time(concrete, time)) * concrete.max_amplitude
+amplitude(concrete::ConcreteShape, time::Number) = amplitude(concrete.shape, convert_time(concrete, time)) * concrete.max_amplitude + concrete.shift
 amplitude_derivative(concrete::ConcreteShape, time::Number) = amplitude_derivative(concrete.shape, convert_time(concrete, time)) * concrete.max_amplitude / (concrete.t1 - concrete.t0)
-amplitude_integral(concrete::ConcreteShape, t0::Number, t1::Number) = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time(concrete, t1)) * concrete.max_amplitude * (concrete.t1 - concrete.t0)
-amplitude_integral(concrete::ConcreteShape, t0::T, times::AbstractVector{T}) where {T<:Number} = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time.(concrete, times)) .* (concrete.max_amplitude * (concrete.t1 - concrete.t0))
-amplitude_integral(concrete::ConcreteShape, t0::T, times::AbstractRange{T}) where {T<:Number} = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time(concrete, times)) .* (concrete.max_amplitude * (concrete.t1 - concrete.t0))
+amplitude_integral(concrete::ConcreteShape, t0::Number, t1::Number) = (amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time(concrete, t1)) * concrete.max_amplitude + concrete.shift) * (concrete.t1 - concrete.t0)
+amplitude_integral(concrete::ConcreteShape, t0::T, times::AbstractVector{T}) where {T<:Number} = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time.(concrete, times)) .* (concrete.max_amplitude * (concrete.t1 - concrete.t0)) .+ (concrete.shift * (concrete.t1 - concrete.t0))
+amplitude_integral(concrete::ConcreteShape, t0::T, times::AbstractRange{T}) where {T<:Number} = amplitude_integral(concrete.shape, convert_time(concrete, t0), convert_time(concrete, times)) .* (concrete.max_amplitude * (concrete.t1 - concrete.t0)) .+ (concrete.shift * (concrete.t1 - concrete.t0))
 control_points(concrete::ConcreteShape) = control_points(concrete.shape) .* (concrete.t1 - concrete.t0) .+ concrete.t0
 
-add_TR(concrete::ConcreteShape, TR::Number) = ConcreteShape(concrete.t0 + TR, concrete.t1 + TR, concrete.max_amplitude, concrete.shape)
+add_TR(concrete::ConcreteShape, TR::Number) = ConcreteShape(concrete.t0 + TR, concrete.t1 + TR, concrete.max_amplitude, concrete.shift, concrete.shape)
 
 start_time(concrete::ConcreteShape) = concrete.t0
 end_time(concrete::ConcreteShape) = concrete.t1
diff --git a/test/example_pulseseq/fiddisp.seq b/test/example_pulseseq/fiddisp.seq
new file mode 100644
index 0000000000000000000000000000000000000000..bfab4c79c7acfeccf22ba8d7ee60b7cb15c38e07
--- /dev/null
+++ b/test/example_pulseseq/fiddisp.seq
@@ -0,0 +1,55 @@
+# Pulseq sequence file
+# Created by MATLAB mr toolbox
+
+[VERSION]
+major 1
+minor 3
+revision 1
+
+[DEFINITIONS]
+Name fid 
+
+# Format of blocks:
+#  #  D RF  GX  GY  GZ ADC EXT
+[BLOCKS]
+1  0  1   0   0   0  0  0
+2  1  0   0   0   0  0  0
+3  0  0   0   0   0  1  0
+
+# Format of RF events:
+# id amplitude mag_id phase_id delay freq phase
+# ..        Hz   ....     ....    us   Hz   rad
+[RF]
+1         2500 1 2 100 0 0
+
+# Format of ADC events:
+# id num dwell delay freq phase
+# ..  ..    ns    us   Hz   rad
+[ADC]
+1 1024 312500 20 0 0
+
+# Format of delays:
+# id delay (us)
+[DELAYS]
+1 5000
+
+# Sequence Shapes
+[SHAPES]
+
+shape_id 1
+num_samples 120
+1
+0
+0
+97
+-1
+0
+0
+17
+
+shape_id 2
+num_samples 120
+0
+0
+118
+
diff --git a/test/runtests.jl b/test/runtests.jl
index 5aee1ea9500e89868a88db9dddd3ce881dbf29c5..b3b1f00d037b8d609a055e3770c925124e5b931b 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -17,7 +17,8 @@ all_tests = [
     "transfer",
     "permeability",
     "radio_frequency",
-    "hierarchical_mri"
+    "hierarchical_mri",
+    "pulseseq",
 ]
 
 if length(ARGS) == 0
diff --git a/test/test_pulseseq.jl b/test/test_pulseseq.jl
new file mode 100644
index 0000000000000000000000000000000000000000..823c8be95dcd6e787797d5e4e78792b275a0188a
--- /dev/null
+++ b/test/test_pulseseq.jl
@@ -0,0 +1,38 @@
+@testset "Test pulseseq input" begin
+    directory = joinpath(pwd(), "example_pulseseq")
+    @testset "read fiddisp.seq" begin
+        keywords = open(mr.read_pulseseq_sections, joinpath(directory, "fiddisp.seq"))
+        seq = mr.build_sequence(;keywords...);
+
+        # starting with RF pulse
+        @test length(seq.pulses) == 1
+        pulse = seq.pulses[1]
+        @show pulse
+        @test mr.start_time(pulse) == 0.1  # start after 100 microseconds delay
+        @test mr.end_time(pulse) == 0.22  # pulse lasts 120 microseconds
+        for cshape in (pulse.amplitude, pulse.phase)
+            @test mr.start_time(cshape) == 0.1
+            @test mr.end_time(cshape) == 0.22 
+        end
+        for (time, ampl) in [
+            (0.09, 0.),
+            (0.11, 2.5),
+            (0.19, 2.5),
+            (0.21, 0.),
+            (0.23, 0.),
+        ]
+            @test iszero(mr.phase(pulse, time))
+            @test mr.amplitude(pulse, time) ≈ ampl
+        end
+
+        # no gradients
+        @test mr.control_points(seq.gradient) == [0, seq.TR]
+
+        # Single readout
+        @test length(seq.readout_times) == 1
+        rt = seq.readout_times[1]
+        @assert rt ≈ 0.22 + 5 + 0.02 + (1024 * 0.3125) / 2
+
+        @assert seq.TR ≈ 0.22 + 5 + 0.02 + (1024 * 0.3125)
+    end
+end
\ No newline at end of file