<reponame>kmjohnson3/sigpy
# -*- coding: utf-8 -*-
"""Multiband RF Pulse Design functions.
"""
import numpy as np

__all__ = ['mb_phs_tab', 'mb_rf', 'dz_pins']

from sigpy.mri.rf.trajgrad import trap_grad
from sigpy.mri.rf import slr as slr


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:
         <NAME>. (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.
         <NAME>., <NAME>., and <NAME>. (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.')

        # <NAME>'s phases: From <NAME>, 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 <NAME>, 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


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:
        <NAME>. and <NAME>. and <NAME> and <NAME> (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
