Extensions to inversion pulse modules (adiabatic HSN pulse, overdrive, get inversion time)

[markus-zimmermann]

We are currently rewriting our QRAGE sequence (https://onlinelibrary.wiley.com/doi/full/10.1002/mrm.30272) in PyPulseq to make it publicly accessible. As part of this, we had to make 3 minor modifications to the PyPulseq code and would like to discuss the possibility to upstream those changes.

First, we added hypsecn as a new option to make_adiabatic_pulse. The implementation is based PyPulseq's _hypsec and on mikgroup/sigpy#81.

def _hypsec_n(
    n: int = 512,
    beta: float = 800.0,
    mu: float = 4.9,
    dur: float = 0.012,
    pwr: float = 1.0,
):
    r"""
    Design a generalized hyperbolic secant (HSN) adiabatic pulse.

    Extends the traditional hyperbolic secant pulse to allow control over the shape
    via the `pwr` parameter, which adjusts the exponent of the hyperbolic argument.

    Args:
        n (int): Number of samples (should be a multiple of 4).
        beta (float): AM waveform parameter (steepness of the envelope).
        mu (float): Constant determining amplitude of frequency sweep.
        dur (float): Pulse duration (s).
        pwr (float): Order of the hyperbolic secant function. Default is 1 (standard hyperbolic secant).

    Returns
    -------
        tuple:
        - **a** (numpy.ndarray): AM waveform (envelope).
        - **om** (numpy.ndarray): FM waveform (frequency modulation in radians/s).
    """

    t = np.arange(-n // 2, n // 2) / n * dur

    a = np.cosh((beta * t) ** pwr) ** -1
    om = -2 * mu * beta * (np.cumsum(a**2) / np.sum(a**2) - 0.5)

    return a, om

We have tested this code on our 7T TERRA scanner with amp_mod, freq_mod = _hypsec_n(n=512, beta=400, mu=11.8, dur=0.012, pwr=4) and it worked very well. Also, for the default values, it reproduces _hypsec exactly.

Secondly, we added an overdrive option to make_adiabatic_pulse, that helped us with poor inversion efficiency, particularly in the cerebellum and the temporal lobes.

def make_adiabatic_pulse(
    ...
    overdrive: float = 1.0,
    ...
) -> Union[
    SimpleNamespace,
    Tuple[SimpleNamespace, SimpleNamespace, SimpleNamespace],
]:

    ...
    # Create the modulated signal
    signal = amp * amp_mod * np.exp(1j * phase_mod) * overdrive
    ....

And finally, because we want to do qMRI, we wanted to automatically calculate the inversion time of our sequence, but rf_times only returns the RF and refocusing pulses, but not the inversion pulses, so extended it, to also return those.

def rf_times(
    seq,
    time_range: Union[List[float], None] = None,
) -> Tuple[List[float], np.ndarray, List[float], np.ndarray, np.ndarray]:
    """
    Return time points of excitations, refocusings and inversions.

    Returns
    -------
    t_excitation : List[float]
        Contains time moments of the excitation RF pulses
    fp_excitation : np.ndarray
        Contains frequency and phase offsets of the excitation RF pulses
    t_refocusing : List[float]
        Contains time moments of the refocusing RF pulses
    fp_refocusing : np.ndarray
        Contains frequency and phase offsets of the refocusing RF pulses
    t_inversion : List[float]
        Contains time moments of the inversion RF pulses
    fp_inversion : np.ndarray
        Contains frequency and phase offsets of the inversion RF pulses
    """
    # Collect RF timing data
    t_excitation = []
    fp_excitation = []
    t_refocusing = []
    fp_refocusing = []
    t_inversion = []
    fp_inversion = []

    curr_dur = 0
    if time_range is None:
        blocks = seq.block_events
    else:
        if len(time_range) != 2:
            raise ValueError("Time range must be list of two elements")
        if time_range[0] > time_range[1]:
            raise ValueError("End time of time_range must be after begin time")

        # Calculate end times of each block
        bd = np.array(list(seq.block_durations.values()))
        t = np.cumsum(bd)
        # Search block end times for start of time range
        begin_block = np.searchsorted(t, time_range[0])
        # Search block begin times for end of time range
        end_block = np.searchsorted(t - bd, time_range[1], side="right")
        blocks = list(seq.block_durations.keys())[begin_block:end_block]
        curr_dur = t[begin_block] - bd[begin_block]

    for block_counter in blocks:
        block = seq.get_block(block_counter)

        if block.rf is not None:
            rf = block.rf
            t = rf.delay + calc_rf_center(rf)[0]
            if not hasattr(rf, "use") or block.rf.use in [
                "excitation",
                "undefined",
            ]:
                t_excitation.append(curr_dur + t)
                fp_excitation.append([block.rf.freq_offset, block.rf.phase_offset])
            elif block.rf.use == "refocusing":
                t_refocusing.append(curr_dur + t)
                fp_refocusing.append([block.rf.freq_offset, block.rf.phase_offset])
            elif block.rf.use == "inversion":
                t_inversion.append(curr_dur + t)
                fp_inversion.append([block.rf.freq_offset, block.rf.phase_offset])

        curr_dur += seq.block_durations[block_counter]

    if len(t_excitation) != 0:
        fp_excitation = np.array(fp_excitation).T
    else:
        fp_excitation = np.empty((2, 0))

    if len(t_refocusing) != 0:
        fp_refocusing = np.array(fp_refocusing).T
    else:
        fp_refocusing = np.empty((2, 0))
    
    if len(t_inversion) != 0:
        fp_inversion = np.array(fp_inversion).T
    else:
        fp_inversion = np.empty((2, 0))

    return t_excitation, fp_excitation, t_refocusing, fp_refocusing, t_inversion, fp_inversion

If there is any interest in upstreaming these changes, I'd be happy to submit a pull request.