Compute RF pulse frequency at specific times

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)

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]

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

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)
 

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)

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)

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.)