new HSn desginer!

docs/mri_rf.rst

@@ -11,6 +11,7 @@ Adiabatic Pulse Design Functions
     :nosignatures:
 
     sigpy.mri.rf.adiabatic.bir4
+    sigpy.mri.rf.adiabatic.hypsec_n
     sigpy.mri.rf.adiabatic.hypsec
     sigpy.mri.rf.adiabatic.wurst
     sigpy.mri.rf.adiabatic.goia_wurst

sigpy/mri/rf/adiabatic.py

@@ -4,7 +4,39 @@
 """
 import numpy as np
 
-__all__ = ['bir4', 'hypsec', 'wurst', 'goia_wurst', 'bloch_siegert_fm']
+__all__ = ['bir4', 'hypsec', 'hypsec_n', 'wurst', 'goia_wurst',
+           'adiabatic_bs_fm']
+
+def hypsec_n(n=500, beta=8, a_max=12000, pwr=8):
+    r"""Design a HS-n hyperbolic secant adiabatic pulse.
+
+    Args:
+        n (int): number of samples in the final pulse.
+        beta (float): AM waveform parameter (rad/ms).
+        a_max (float): maximum values of the frequency sweep (Hz).
+        pwr (int): power to exponentiate time by within the
+            hyperbolic sinc.
+
+    Returns:
+        2-element tuple containing
+
+        - **am_sechn** (*array*): AM waveform.
+        - **fm_sechn** (*array*): FM waveform (Hz/ms).
+
+    References:
+        Tannus, A. and Garwood, M. 'Improved performance
+        of frequency-swept pulses using offset-independent
+        adiabaticity'. J. Magn. Reson. A 120, 133-137 (1996).
+     """
+
+    t = 2 * np.arange(-n // 2, n // 2) / n
+    dt = t[1]-t[0]
+
+    am_sechn = np.cosh(beta * t ** pwr) ** -1
+    f2_sechn = np.cumsum(am_sechn ** 2) * dt
+    fm_sechn = -2 * a_max * (f2_sechn / np.max(f2_sechn) - 1 / 2)
+
+    return am_sechn, fm_sechn
 
 
 def bir4(n, beta, kappa, theta, dw0):

tests/mri/rf/test_adiabatic.py

@@ -10,6 +10,37 @@
 
 class TestAdiabatic(unittest.TestCase):
 
+    def test_hypsec_n(self):
+        # test an HSn adiabatic inversion pulse
+        n = 500
+        beta = 8  # rad/ms
+        a_max = 12000  # Hz
+        pwr = 8
+        dur = 0.012  # s
+
+        [am_sechn, fm_sechn] = rf.adiabatic.hypsec_n(n, beta, a_max, pwr)
+
+        # check relatively homogeneous over range of B1 values
+        b1 = np.arange(0.2, 0.8, 0.1)
+        b1 = np.reshape(b1, (np.size(b1), 1))
+
+        a = np.zeros(np.shape(b1), dtype='complex')
+        b = np.zeros(np.shape(b1), dtype='complex')
+        for ii in range(0, np.size(b1)):
+            [a[ii], b[ii]] = rf.sim.abrm_nd(
+                2 * np.pi * (dur / n) * 4258 * b1[ii] * am_sechn, np.ones(1),
+                dur / n * np.reshape(fm_sechn, (np.size(fm_sechn), 1)))
+        mz = 1 - 2 * np.abs(b) ** 2
+
+        test = np.ones(mz.shape) * -1  # magnetization value we expect
+
+        # test we get our inversion profile
+        npt.assert_array_almost_equal(mz, test, 2)
+        # test that the AM is scaled to 1
+        npt.assert_almost_equal(np.max(am_sechn), 1, 3)
+        # test that the FM is scaled to A
+        npt.assert_almost_equal(np.max(fm_sechn), a_max, 3)
+
     def test_bir4(self):
         # test an excitation bir4 pulse
         n = 1176