Add hypersec_n pulse type and overdrive factor to make_adiabatic_pulse and include inversion pulses in rf_times

src/pypulseq/Sequence/sequence.py

@@ -293,7 +293,7 @@ def calculate_kspace(
 
         total_duration = sum(self.block_durations.values())
 
-        t_excitation, fp_excitation, t_refocusing, _ = self.rf_times()
+        t_excitation, fp_excitation, t_refocusing, _, _, _ = self.rf_times()
         t_adc, _ = self.adc_times()
 
         # Convert data to piecewise polynomials
@@ -1239,7 +1239,7 @@ def rf_times(
         self, time_range: Union[List[float], None] = None
     ) -> Tuple[List[float], np.ndarray, List[float], np.ndarray, np.ndarray]:
         """
-        Return time points of excitations and refocusings.
+        Return time points of excitations, refocusings and inversions.
 
         Returns
         -------
@@ -1250,13 +1250,19 @@ def rf_times(
         t_refocusing : List[float]
             Contains time moments of the refocusing RF pulses
         fp_refocusing : np.ndarray
-            Contains frequency and phase offsets of the excitation RF pulses
+            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:
@@ -1292,6 +1298,9 @@ def rf_times(
                 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 += self.block_durations[block_counter]
 
@@ -1305,7 +1314,19 @@ def rf_times(
         else:
             fp_refocusing = np.empty((2, 0))
 
-        return t_excitation, fp_excitation, t_refocusing, fp_refocusing
+        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,
+        )
 
     def set_block(self, block_index: int, *args: SimpleNamespace) -> None:
         """
@@ -1578,7 +1599,7 @@ def waveforms_and_times(
             Contains frequency and phase offsets of each ADC object (not samples).
         """
         wave_data = self.waveforms(append_RF=append_RF, time_range=time_range)
-        t_excitation, fp_excitation, t_refocusing, fp_refocusing = self.rf_times(time_range=time_range)
+        t_excitation, fp_excitation, t_refocusing, fp_refocusing, _, _ = self.rf_times(time_range=time_range)
         t_adc, fp_adc = self.adc_times(time_range=time_range)
 
         # Join times, frequency and phases of RF pulses for compatibility with previous implementation

src/pypulseq/make_adiabatic_pulse.py

@@ -24,8 +24,10 @@ def make_adiabatic_pulse(
     freq_offset: float = 0,
     max_grad: Union[float, None] = None,
     max_slew: Union[float, None] = None,
-    n_fac: int = 40,
     mu: float = 4.9,
+    n_fac: int = 40,
+    order: float = 1.0,
+    overdrive: float = 1.0,
     phase_offset: float = 0,
     return_gz: bool = False,
     slice_thickness: float = 0,
@@ -62,6 +64,22 @@ def make_adiabatic_pulse(
             inversion of a two-level system by phase-modulated pulses'.
             Phys. Rev. A., 32:3435-3447.
 
+    hypsec_n(n=512, beta=400, mu=9.8, order=4, dur=0.012)
+        Design a generalized hyperbolic secant (HSN) adiabatic pulse. `mu` * `beta` becomes the amplitude of the frequency sweep.
+
+        Args:
+            - n (int): number of samples (should be a multiple of 4).
+            - beta (float): AM waveform parameter.
+            - mu (float): a constant, determines amplitude of frequency sweep.
+            - order (float): order of the hyperbolic secant pulse.
+            - dur (float): pulse time (s).
+
+    Returns
+    -------
+            2-element tuple containing
+            - **a** (*array*): AM waveform.
+            - **om** (*array*): FM waveform (radians/s).
+
     wurst(n=512, n_fac=40, bw=40000.0, dur=0.002)
         Design a WURST (wideband, uniform rate, smooth truncation) adiabatic inversion pulse
 
@@ -106,8 +124,12 @@ def make_adiabatic_pulse(
         Constant determining amplitude of frequency sweep.
     n_fac : int, default=40
         Power to exponentiate to within AM term. ~20 or greater is typical.
+    overdrive : float, default=1
+        Overdrive factor.
     phase_offset : float, default=0
         Phase offset.
+    order: float, default=1.0,
+        Order of the generalized hyperbolic secant (HSN) adiabatic pulse.
     return_gz : bool, default=False
         Boolean flag to indicate if the slice-selective gradient has to be returned.
     slice_thickness : float, default=0
@@ -138,7 +160,7 @@ def make_adiabatic_pulse(
     if return_gz and slice_thickness <= 0:
         raise ValueError('Slice thickness must be provided')
 
-    valid_pulse_types = ['hypsec', 'wurst']
+    valid_pulse_types = ['hypsec', 'hypsec_n', 'wurst']
     if (not pulse_type) or (pulse_type not in valid_pulse_types):
         raise ValueError(f'Invalid type parameter. Must be one of {valid_pulse_types}.Passed: {pulse_type}')
     valid_rf_use_labels = get_supported_rf_uses()
@@ -155,6 +177,8 @@ def make_adiabatic_pulse(
 
     if pulse_type == 'hypsec':
         amp_mod, freq_mod = _hypsec(n=n_samples, beta=beta, mu=mu, dur=duration)
+    elif pulse_type == 'hypsec_n':
+        amp_mod, freq_mod = _hypsec_n(n=n_samples, beta=beta, mu=mu, order=order, dur=duration)
     elif pulse_type == 'wurst':
         amp_mod, freq_mod = _wurst(n=n_samples, n_fac=n_fac, bw=bandwidth, dur=duration)
 
@@ -189,7 +213,7 @@ def make_adiabatic_pulse(
     amp = np.sqrt(rate_of_freq_change * adiabaticity) / (2 * np.pi * amp_at_zero_freq)
 
     # Create the modulated signal
-    signal = amp * amp_mod * np.exp(1j * phase_mod)
+    signal = amp * amp_mod * np.exp(1j * phase_mod) * overdrive
 
     # Adjust the number of samples if needed
     if n_samples != n_raw:
@@ -224,7 +248,7 @@ def make_adiabatic_pulse(
         max_grad_slice_select = max_grad
         max_slew_slice_select = max_slew
 
-        if pulse_type == 'hypsec':
+        if pulse_type in ['hypsec', 'hypsec_n']:
             bandwidth = mu * beta / np.pi
         elif pulse_type == 'wurst':
             bandwidth = bandwidth
@@ -504,3 +528,41 @@ def _bloch_siegert_fm(
     om = np.concatenate((om, om[::-1]))
 
     return om
+
+
+def _hypsec_n(
+    n: int = 512,
+    beta: float = 800.0,
+    mu: float = 4.9,
+    order: float = 1.0,
+    dur: float = 0.012,
+):
+    r"""
+    Design a generalized hyperbolic secant (HSN) adiabatic pulse.
+
+    Extends the traditional hyperbolic secant pulse to allow control over the shape
+    via the `order` 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.
+        order (float): Order of the hyperbolic secant function. Default is 1 (standard hyperbolic secant).
+        dur (float): Pulse duration (s).
+
+    Returns
+    -------
+        tuple:
+        - **a** (numpy.ndarray): AM waveform (envelope).
+        - **om** (numpy.ndarray): FM waveform (frequency modulation in radians/s).
+
+    References
+    ----------
+
+    """
+    t = np.arange(-n // 2, n // 2) / n * dur
+
+    a = np.cosh((beta * t) ** order) ** -1
+    om = -2 * mu * beta * (np.cumsum(a**2) / np.sum(a**2) - 0.5)
+
+    return a, om

tests/test_make_adiabatic_pulse.py

@@ -13,7 +13,7 @@
 
 def test_pulse_select():
     valid_rf_use_labels = get_supported_rf_uses()
-    valid_pulse_types = ('hypsec', 'wurst')
+    valid_pulse_types = ('hypsec', 'hypsec_n', 'wurst')
 
     # Check all use and valid pulse combinations return a sensible object
     # with default parameters.
@@ -65,6 +65,12 @@ def test_hypsec_options():
     assert np.isclose(pobj.shape_dur, 0.05)
 
 
+def test_hypsec_n_options():
+    pobj = make_adiabatic_pulse(pulse_type='hypsec_n', beta=400, mu=9.8, order=4, duration=0.01)
+
+    assert np.isclose(pobj.shape_dur, 0.01)
+
+
 def test_wurst_options():
     pobj = make_adiabatic_pulse(pulse_type='wurst', n_fac=25, bandwidth=30000, duration=0.05)