From 9f4dfbc27842f080704a4f710313b56f897d5322 Mon Sep 17 00:00:00 2001
From: Michiel Cottaar <michiel.cottaar@ndcn.ox.ac.uk>
Date: Wed, 17 Apr 2024 15:23:44 +0100
Subject: [PATCH] Compute RF pulse frequency at specific times
---
src/components/pulses/constant_pulses.jl | 1 +
src/components/pulses/generic_pulses.jl | 24 +++++++++++++++++++++++-
src/components/pulses/sinc_pulses.jl | 1 +
src/containers/abstract.jl | 16 +++++++++++++---
src/containers/base_sequences.jl | 4 ++--
src/containers/building_blocks.jl | 4 ++--
test/test_sequences.jl | 4 ++++
7 files changed, 46 insertions(+), 8 deletions(-)
diff --git a/src/components/pulses/constant_pulses.jl b/src/components/pulses/constant_pulses.jl
index f2891a0..f330df0 100644
--- a/src/components/pulses/constant_pulses.jl
+++ b/src/components/pulses/constant_pulses.jl
@@ -48,6 +48,7 @@ effective_time(pulse::ConstantPulse) = duration(pulse) / 2
amplitude(pulse::ConstantPulse, time::Number) = amplitude(pulse)
phase(pulse::ConstantPulse, time::Number) = phase(pulse) + frequency(pulse) * (time - effective_time(pulse))
+frequency(pulse::ConstantPulse, time::Number) = frequency(pulse)
function make_generic(block::ConstantPulse)
d = duration(block)
diff --git a/src/components/pulses/generic_pulses.jl b/src/components/pulses/generic_pulses.jl
index 8142c49..d70c0d7 100644
--- a/src/components/pulses/generic_pulses.jl
+++ b/src/components/pulses/generic_pulses.jl
@@ -2,7 +2,7 @@ module GenericPulses
import Polynomials: fit
import ...AbstractTypes: RFPulseComponent, split_timestep
-import ....Variables: duration, amplitude, effective_time, flip_angle, make_generic, phase
+import ....Variables: duration, amplitude, effective_time, flip_angle, make_generic, phase, frequency
"""
@@ -99,6 +99,28 @@ for fn in (:amplitude, :phase)
end
+function frequency(gp::GenericPulse, time::Number)
+ i2 = findfirst(t -> t > time, gp.time)
+ if isnothing(i2)
+ @assert time ≈ fp.time[end]
+ i2 = length(time)
+ end
+ if i2 != length(time) && time ≈ fp.time[i2 + 1]
+ i2 += 1
+ end
+ if time ≈ fp.time[i2]
+ if i2 == 1
+ return (gp.phase[2] - gp.phase[1]) / (gp.time[2] - gp.time[1]) / 360
+ elseif i2 == length(time)
+ return (gp.phase[end] - gp.phase[end-1]) / (gp.time[end] - gp.time[end-1]) / 360
+ end
+ return (gp.phase[i2 + 1] - gp.phase[i2 - 1]) / (gp.time[i2 + 1] - gp.time[i2 - 1]) / 360
+ end
+ @assert !isone(i2)
+ return (gp.phase[i2] - gp.phase[i2 - 1]) / (gp.time[i2] - gp.time[i2 - 1]) / 360
+end
+
+
function get_weights(pulse::GenericPulse)
weights_double = pulse.time[2:end] - pulse.time[1:end-1]
return [weights_double[1] / 2, ((weights_double[1:end-1] + weights_double[2:end]) / 2)..., weights_double[end]/2]
diff --git a/src/components/pulses/sinc_pulses.jl b/src/components/pulses/sinc_pulses.jl
index 99dfb53..2c41504 100644
--- a/src/components/pulses/sinc_pulses.jl
+++ b/src/components/pulses/sinc_pulses.jl
@@ -92,6 +92,7 @@ effective_time(pulse::SincPulse) = N_left(pulse) * lobe_duration(pulse)
amplitude(pulse::SincPulse, time::Number) = amplitude(pulse) * normalised_function(abs((time - effective_time(pulse))) / lobe_duration(pulse), N_left(pulse), N_right(pulse); apodise=pulse.apodise)
phase(pulse::SincPulse, time::Number) = phase(pulse) + frequency(pulse) * (time - effective_time(pulse))
+frequency(pulse::SincPulse, time::Number) = frequency(pulse)
function make_generic(block::SincPulse)
normed_times = -N_left(block):0.1:N_right(block) + 1e-5
diff --git a/src/containers/abstract.jl b/src/containers/abstract.jl
index 0148c6d..a1dabc3 100644
--- a/src/containers/abstract.jl
+++ b/src/containers/abstract.jl
@@ -72,7 +72,7 @@ function gradient_strength end
"""
amplitude(sequence, time)
-Returns the RF amplitude at a particular time within the sequence.
+Returns the RF amplitude at a particular time within the sequence in kHz.
"""
function amplitude end
@@ -80,11 +80,21 @@ function amplitude end
"""
phase(sequence, time)
-Returns the RF phase at a particular time within the sequence.
-NaN if the ampltiude is zero.
+Returns the RF phase at a particular time within the sequence in degrees.
+
+NaN is returned if there is no pulse activate at that `time`.
"""
function phase end
+"""
+ frequency(sequence, time)
+
+Returns the RF frequency at a particular time within the sequence in kHz.
+
+NaN is returned if there is no pulse activate at that `time`.
+"""
+function frequency end
+
"""
iter(sequence, get_type)
diff --git a/src/containers/base_sequences.jl b/src/containers/base_sequences.jl
index 80603f4..629c348 100644
--- a/src/containers/base_sequences.jl
+++ b/src/containers/base_sequences.jl
@@ -4,7 +4,7 @@ Defines [`BaseSequence`](@ref) and [`Sequence`](@ref)
module BaseSequences
import StaticArrays: SVector
import JuMP: @constraint
-import ...Variables: get_free_variable, repetition_time, VariableType, duration, variables, VariableNotAvailable, Variables, set_simple_constraints!, TR, make_generic, gradient_strength, amplitude, phase, gradient_strength3, get_gradient, get_pulse
+import ...Variables: get_free_variable, repetition_time, VariableType, duration, variables, VariableNotAvailable, Variables, set_simple_constraints!, TR, make_generic, gradient_strength, amplitude, phase, gradient_strength3, get_gradient, get_pulse, frequency
import ...BuildSequences: global_model, global_scanner
import ...Components: EventComponent, NoGradient
import ...Scanners: Scanner, B0
@@ -114,7 +114,7 @@ function edge_times(seq::BaseSequence; tol=1e-6)
return sort(unique_res)
end
-for fn in (:gradient_strength, :amplitude, :phase, :gradient_strength3, :get_gradient, :get_pulse)
+for fn in (:gradient_strength, :amplitude, :phase, :frequency, :gradient_strength3, :get_gradient, :get_pulse)
@eval function $fn(sequence::BaseSequence, time::AbstractFloat)
(block_time, block) = sequence(time)
return $fn(block, block_time)
diff --git a/src/containers/building_blocks.jl b/src/containers/building_blocks.jl
index 6cd930f..9fb9d71 100644
--- a/src/containers/building_blocks.jl
+++ b/src/containers/building_blocks.jl
@@ -8,7 +8,7 @@ import StaticArrays: SVector
import ..Abstract: ContainerBlock, start_time, readout_times, edge_times, end_time, iter
import ...BuildSequences: global_model
import ...Components: BaseComponent, GradientWaveform, EventComponent, NoGradient, ChangingGradient, ConstantGradient, split_gradient, RFPulseComponent, ReadoutComponent, InstantGradient
-import ...Variables: qval, bmat_gradient, effective_time, get_free_variable, qval3, slew_rate, gradient_strength, amplitude, phase
+import ...Variables: qval, bmat_gradient, effective_time, get_free_variable, qval3, slew_rate, gradient_strength, amplitude, phase, frequency
import ...Variables: VariableType, duration, make_generic, get_pulse, get_readout, scanner_constraints!, get_gradient
"""
@@ -230,7 +230,7 @@ function get_pulse(bb::BaseBuildingBlock, time::Number)
return nothing
end
-for (fn, default_value) in ((:amplitude, 0.), (:phase, NaN))
+for (fn, default_value) in ((:amplitude, 0.), (:phase, NaN), (:frequency, NaN))
@eval function $fn(bb::BaseBuildingBlock, time::Number)
pulse = get_pulse(bb, time)
if isnothing(pulse)
diff --git a/test/test_sequences.jl b/test/test_sequences.jl
index 6f741aa..d82c008 100644
--- a/test/test_sequences.jl
+++ b/test/test_sequences.jl
@@ -117,6 +117,10 @@
gp = GenericPulse(pulse, 0., 1.)
@test gp.amplitude[1] ≈ 0. atol=1e-8
@test gp.amplitude[end] ≈ amplitude(pulse, 1.) rtol=1e-2
+ @test iszero(phase(pulse, 1.)) rtol=1e-2
+ @test iszero(frequency(pulse, 1.)) rtol=1e-2
+ @test isnan(phase(pulse, 1.)) rtol=1e-2
+ @test isnan(frequency(pulse, 1.)) rtol=1e-2
@test all(iszero.(gp.phase))
(pulse, t_pulse) = get_pulse(seq, 35.)
--
GitLab