diff --git a/src/components/pulses/constant_pulses.jl b/src/components/pulses/constant_pulses.jl index f2891a01223a6761e74b272829cd5342d50525b3..f330df07ebfa36fb964cf8f06397aab272927f26 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 8142c496fd449ab9bbab93287b1e208402bba89c..d70c0d7f33c1df4768398e24a8a8f1d4a9b30a2f 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 99dfb532da82c302a0a676f2092e647e37474e4f..2c41504a2f81301c3664627efd0e11744d4907d5 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 0148c6de329d0e418d5a8c3628b261359c76c560..a1dabc3e4d4ea992aff0f7260b591b0bb98e9f6e 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 80603f483553bef97d2896789b03f3dfc9160b80..629c348637d44a74d4361bb8f739aee9dd3387ac 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 6cd930f44cf95947d08de956986e48941f099661..9fb9d71b4e2f1fe6ee41ac068e26f4e1bc87bfa0 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 6f741aa78de9ee8db47be0e78e9e2fbcb59294f1..d82c008203c0ce3e6c4e231fb144c5c62b99c049 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.)