# -*- coding: utf-8 -*-
"""Multiband RF Pulse Design functions.
"""
import numpy as np
__all__ = ["mb_phs_tab", "mb_rf", "dz_pins"]
from sigpy.mri.rf import slr as slr
from sigpy.mri.rf.trajgrad import trap_grad
[docs]def mb_rf(pulse_in, n_bands=3, band_sep=5, phs_0_pt="None"):
r"""Multiband an input RF pulse.
Args:
pulse_in (array): samples of single-band RF pulse.
n_bands (int): number of bands.
band_sep (float): normalized slice separation.
phs_0_pt (string): set of phases to use. Can be 'phs_mod' (Wong),
'amp_mod' (Malik), 'quad_mod' (Grissom), or 'None'
band_sep = slice_sep/slice_thick*tb, where tb is time-bandwidth product
of the single-band pulse
Returns:
array: multibanded pulse out
References:
Wong, E. (2012). 'Optimized Phase Schedules for Minimizing Peak RF
Power in Simultaneous Multi-Slice RF Excitation Pulses'. Proc. Intl.
Soc. Mag. Reson. Med., 20 p. 2209.
Malik, S. J., Price, A. N., and Hajnal, J. V. (2015). 'Optimized
Amplitude Modulated Multi-Band RF pulses'. Proc. Intl. Soc. Mag.
Reson. Med., 23 p. 2398.
"""
if phs_0_pt != "None":
phs = mb_phs_tab(n_bands, phs_0_pt)
else:
phs = np.zeros(n_bands)
# build multiband modulation function
n = np.size(pulse_in)
b = np.zeros(n, dtype="complex")
for ii in range(0, n_bands):
b += np.exp(
1j
* 2
* np.pi
/ n
* band_sep
* np.arange(-n / 2, n / 2, 1)
* (ii - (n_bands - 1) / 2)
) * np.exp(1j * phs[ii])
pulse_out = b * pulse_in
return pulse_out
def mb_phs_tab(n_bands, phs_type="phs_mod"):
# Return phases to minimize peak b1 amplitude of an MB pulse
if phs_type == "phs_mod":
if n_bands < 3 or n_bands > 16:
raise Exception("Wongs phases valid for 2 < nBands < 17.")
# Eric Wong's phases: From E C Wong, ISMRM 2012, p. 2209
p = np.zeros((14, 16))
p[0, 1:3] = np.array([0.73, 4.602])
p[1, 1:4] = np.array([3.875, 5.94, 6.197])
p[2, 1:5] = np.array([3.778, 5.335, 0.872, 0.471])
p[3, 1:6] = np.array([2.005, 1.674, 5.012, 5.736, 4.123])
p[4, 1:7] = np.array([3.002, 5.998, 5.909, 2.624, 2.528, 2.440])
p[5, 1:8] = np.array([1.036, 3.414, 3.778, 3.215, 1.756, 4.555, 2.467])
p[6, 1:9] = np.array(
[1.250, 1.783, 3.558, 0.739, 3.319, 1.296, 0.521, 5.332]
)
p[7, 1:10] = np.array(
[4.418, 2.360, 0.677, 2.253, 3.472, 3.040, 3.974, 1.192, 2.510]
)
p[8, 1:11] = np.array(
[
5.041,
4.285,
3.001,
5.765,
4.295,
0.056,
4.213,
6.040,
1.078,
2.759,
]
)
p[9, 1:12] = np.array(
[
2.755,
5.491,
4.447,
0.231,
2.499,
3.539,
2.931,
2.759,
5.376,
4.554,
3.479,
]
)
p[10, 1:13] = np.array(
[
0.603,
0.009,
4.179,
4.361,
4.837,
0.816,
5.995,
4.150,
0.417,
1.520,
4.517,
1.729,
]
)
p[11, 1:14] = np.array(
[
3.997,
0.830,
5.712,
3.838,
0.084,
1.685,
5.328,
0.237,
0.506,
1.356,
4.025,
4.483,
4.084,
]
)
p[12, 1:15] = np.array(
[
4.126,
2.266,
0.957,
4.603,
0.815,
3.475,
0.977,
1.449,
1.192,
0.148,
0.939,
2.531,
3.612,
4.801,
]
)
p[13, 1:16] = np.array(
[
4.359,
3.510,
4.410,
1.750,
3.357,
2.061,
5.948,
3.000,
2.822,
0.627,
2.768,
3.875,
4.173,
4.224,
5.941,
]
)
out = p[n_bands - 3, 0:n_bands]
elif phs_type == "amp_mod":
# Malik's Hermitian phases: From S J Malik, ISMRM 2015, p. 2398
if n_bands < 4 or n_bands > 12:
raise Exception("Maliks phases valid for 3 < nBands < 13.")
p = np.zeros((9, 12))
p[0, 0:4] = np.array([0, np.pi, np.pi, 0])
p[1, 0:5] = np.array([0, 0, np.pi, 0, 0])
p[2, 0:6] = np.array([1.691, 2.812, 1.157, -1.157, -2.812, -1.691])
p[3, 0:7] = np.array([2.582, -0.562, 0.102, 0, -0.102, 0.562, -2.582])
p[4, 0:8] = np.array(
[2.112, 0.220, 1.464, 1.992, -1.992, -1.464, -0.220, -2.112]
)
p[5, 0:9] = np.array(
[0.479, -2.667, -0.646, -0.419, 0, 0.419, 0.646, 2.667, -0.479]
)
p[6, 0:10] = np.array(
[
1.683,
-2.395,
2.913,
0.304,
0.737,
-0.737,
-0.304,
-2.913,
2.395,
-1.683,
]
)
p[7, 0:11] = np.array(
[
1.405,
0.887,
-1.854,
0.070,
-1.494,
0,
1.494,
-0.070,
1.854,
-0.887,
-1.405,
]
)
p[8, 0:12] = np.array(
[
1.729,
0.444,
0.722,
2.190,
-2.196,
0.984,
-0.984,
2.196,
-2.190,
-0.722,
-0.444,
-1.729,
]
)
out = p[n_bands - 4, 0:n_bands]
elif phs_type == "quad_mod":
# Grissom's quadratic phases (unpublished)
k = 3.4 / n_bands # quadratic phase coefficient
out = k * (np.arange(0, n_bands, 1) - (n_bands - 1) / 2) ** 2
else:
raise Exception('phase type ("{}") not recognized.'.format(phs_type))
return out
[docs]def dz_pins(
tb,
sl_sep,
sl_thick,
g_max,
g_slew,
dt,
b1_max=0.18,
ptype="ex",
ftype="ls",
d1=0.01,
d2=0.01,
gambar=4258,
):
r"""PINS multiband pulse design.
Args:
tb (float): time-bandwidth product.
sl_sep (float): slice separation in cm.
sl_thick (float): slice thickness in cm.
g_max (float): max gradient amplitude in gauss/cm
g_slew (float): max gradient sliew in gauss/cm/s
dt (float): RF + gradient dwell time in s.
b1_max (float): Maximum RF amplitude
ptype (string): pulse type, 'st' (small-tip excitation), 'ex' (pi/2
excitation pulse), 'se' (spin-echo pulse), 'inv' (inversion), or
'sat' (pi/2 saturation pulse).
ftype (string): type of filter to use in pulse design
d1 (float): passband ripple level in :math:'M_0^{-1}'.
d2 (float): stopband ripple level in :math:'M_0^{-1}'.
gambar (float): Appropriate gyromagnetic ratio in Hz/gauss.
Returns:
2-element tuple containing:
- **rf** (*array*): RF Pulse in Gauss
- **g** (*array*): Gradient waveform in Gauss/cm
References:
Norris, D.G. and Koopmans, P.J. and Boyacioglu, R and Barth, M (2011).
'Power independent of number of slices (PINS) radiofrequency Pulses
for low-power simultaneous multislice excitation'.
Magn. Reson. Med., 66(5):1234-1240.
"""
kz_width = tb / sl_thick # 1/cm, width in k-space we must go
# calcualte number of subpulses (odd)
n_pulses = int(2 * np.floor(np.ceil(kz_width / (1 / sl_sep)) / 2))
# call SLR to get envelope
rf_soft = slr.dzrf(n_pulses, tb, ptype, ftype, d1, d2)
# design the blip trapezoid
g_area = 1 / sl_sep / gambar
[gz_blip, _] = trap_grad(g_area, g_max, g_slew, dt)
# Calculate the block/hard RF pulse width based on
b1_scaled = 2 * np.pi * gambar * b1_max * dt
hpw = int(np.ceil(np.max(np.abs(rf_soft)) / b1_scaled))
# interleave RF subpusles with gradient subpulses to form full pulses
rf = np.kron(
rf_soft[:-1],
np.concatenate((np.ones(hpw), np.zeros((np.size(gz_blip))))),
)
rf = np.concatenate((rf, rf_soft[-1] * np.ones(hpw)))
rf = rf / (np.sum(rf) * 2 * np.pi * gambar * dt) * np.sum(rf_soft)
g = np.concatenate([np.zeros(hpw), np.squeeze(gz_blip)])
g = np.tile(g, n_pulses - 1)
g = np.concatenate((g, np.zeros(hpw)))
return rf, g