{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# B1 Selective RF Pulse Design"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import sigpy as sp\n",
    "import sigpy.mri as mr\n",
    "import sigpy.mri.rf as rf\n",
    "import sigpy.plot as pl\n",
    "import scipy.signal as signal\n",
    "import matplotlib.pyplot as pyplot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 2e-6 # sampling period\n",
    "d1 = 0.01 # passband ripple\n",
    "d2 = 0.01 # stopband ripple\n",
    "tb = 4 # time-bandwidth product\n",
    "ptype = 'ex' # 'st', 'ex', 'inv' or 'sat'\n",
    "pbw = 0.5 # gauss, passband width\n",
    "pbc = 5 # gauss, passband center\n",
    "flip = np.pi/4 # radians, flip angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Design Excitation Pulse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "[rf_am, rf_fm] = rf.b1sel.dz_b1_rf(dt, tb, ptype, flip, pbw, pbc, d1, d2)\n",
    "T = np.size(rf_fm)*dt\n",
    "t = np.arange(-T/2,T/2,dt)*1000\n",
    "pyplot.figure()\n",
    "pyplot.plot(t,rf_fm)\n",
    "pyplot.xlabel('ms')\n",
    "pyplot.ylabel('Hz')\n",
    "pyplot.title('FM Waveform')\n",
    "pyplot.figure()\n",
    "pyplot.plot(t,rf_am)\n",
    "pyplot.xlabel('ms')\n",
    "pyplot.ylabel('a.u.')\n",
    "pyplot.title('AM Waveform')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulate the excitation pulse's Mxy profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = np.size(rf_am)\n",
    "b1 = np.arange(0, 2*pbc, 2*pbc/np.size(rf_am)*4) # b1 grid we simulate the pulse over\n",
    "b1 = np.reshape(b1, (np.size(b1),1))\n",
    "[a, b] = rf.sim.abrm_nd(2*np.pi*dt*rf_fm, b1, 2*np.pi*4258*dt*np.reshape(rf_am, (np.size(rf_am),1)))\n",
    "Mxy = -2*np.real(a*b) + 1j*np.imag(np.conj(a)**2 - b**2)\n",
    "pyplot.figure()\n",
    "pyplot.plot(b1, np.abs(Mxy))\n",
    "pyplot.xlabel('Gauss')\n",
    "pyplot.ylabel('|Mxy|')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Design gSlider excitation pulses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = 5\n",
    "flip = np.pi/2\n",
    "ptype = 'ex' # 'ex' or 'st'\n",
    "tb = 12\n",
    "d1 = 0.01\n",
    "d2 = 0.01\n",
    "pbc = 1 # gauss, passband center\n",
    "pbw = 0.25 # passband width\n",
    "dt = 2e-6 # seconds, sampling rate\n",
    "[om1, dom] = rf.b1sel.dz_b1_gslider_rf(dt, G, tb, ptype, flip, pbw, pbc, d1, d2)\n",
    "\n",
    "# plot the first pulse\n",
    "pyplot.figure()\n",
    "pyplot.plot(om1[:, 0])\n",
    "pyplot.figure()\n",
    "pyplot.plot(dom[:, 0])\n",
    "pyplot.xlabel('ms')\n",
    "pyplot.ylabel('a.u.')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulate the gSlider pulse's excitation profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = np.shape(om1)[0]\n",
    "b1 = np.arange(0, 2*pbc, 2*pbc/n*4) # b1 grid we simulate the pulse over\n",
    "b1 = np.reshape(b1, (np.size(b1),1))\n",
    "[a, b] = rf.sim.abrm_nd(2*np.pi*dt*dom[:, 0], b1, 2*np.pi*4258*dt*np.reshape(om1[:, 0], (n, 1)))\n",
    "Mxy = -2*np.real(a*b) + 1j*np.imag(np.conj(a)**2 - b**2)\n",
    "pyplot.figure()\n",
    "pyplot.plot(b1, np.abs(Mxy))\n",
    "pyplot.plot(b1, np.real(Mxy))\n",
    "pyplot.plot(b1, np.imag(Mxy))\n",
    "pyplot.xlabel('Gauss')\n",
    "pyplot.ylabel('|Mxy|')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Design Hadamard Encoding Pulses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "G = 8\n",
    "flip = np.pi/6\n",
    "ptype = 'st' # 'ex' or 'st'\n",
    "tb = 24\n",
    "d1 = 0.01\n",
    "d2 = 0.01\n",
    "pbc = 2 # gauss, passband center\n",
    "pbw = 2 # passband width\n",
    "dt = 2e-6 # seconds, sampling rate\n",
    "[om1, dom] = rf.b1sel.dz_b1_hadamard_rf(dt, G, tb, ptype, flip, pbw, pbc, d1, d2)\n",
    "\n",
    "# plot the first pulse\n",
    "pyplot.figure()\n",
    "pyplot.plot(om1[:, 0])\n",
    "pyplot.figure()\n",
    "pyplot.plot(dom[:, 0])\n",
    "pyplot.xlabel('ms')\n",
    "pyplot.ylabel('a.u.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulate the Hadamard pulse's excitation profile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = np.shape(om1)[0]\n",
    "b1 = np.arange(0, 2*pbc, 2*pbc/n*4) # b1 grid we simulate the pulse over\n",
    "b1 = np.reshape(b1, (np.size(b1),1))\n",
    "[a, b] = rf.sim.abrm_nd(2*np.pi*dt*dom[:, 5], b1, 2*np.pi*4258*dt*np.reshape(om1[:, 5], (n, 1)))\n",
    "Mxy = -2*np.real(a*b) + 1j*np.imag(np.conj(a)**2 - b**2)\n",
    "pyplot.figure()\n",
    "pyplot.plot(b1, np.abs(Mxy))\n",
    "pyplot.plot(b1, np.real(Mxy))\n",
    "pyplot.plot(b1, np.imag(Mxy))\n",
    "pyplot.xlabel('Gauss')\n",
    "pyplot.ylabel('|Mxy|')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}