test_verse_flatten.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VERSE for B1+ Selective RF Pulses\n",
    "This notebook was inspired by a recent research project! We have experimented with a class of B1+-selective RF pulses as a part of B0-gradient free imaging systems. These systems often use low-cost amplifiers, and the particular class of amplifier we were investigating was only able to produce constant-amplitude waveforms. This notebook designs a B1+-selective pulse by the Bloch-Siegert localized method, then uses VERSE to flatten the amplitude of the pulse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pulpy as pp\n",
    "import pulpy.rf as rf\n",
    "import pulpy.grad as grad\n",
    "import numpy as np\n",
    "\n",
    "from sigpy import plot as pl\n",
    "from scipy.io import loadmat, savemat \n",
    "from scipy.interpolate import CubicSpline\n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# design a bssel pulse \n",
    "dtb1sel = 1e-6\n",
    "bsrf, rfp, rw = rf.b1sel.dz_bssel_rf(dt=dtb1sel, bs_offset=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rfb1sel = bsrf + rfp\n",
    "plt.plot(np.abs(rfb1sel.flatten()))\n",
    "plt.plot(np.real(rfb1sel.flatten()))\n",
    "plt.plot(np.imag(rfb1sel.flatten()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rfflat = pp.verse.flatten_bsse(rfb1sel, dtb1sel)\n",
    "plt.plot(np.abs(rfflat.flatten()))\n",
    "plt.plot(np.real(rfflat.flatten()))\n",
    "plt.plot(np.imag(rfflat.flatten()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rfSimCplx(rf, gamgdt, xx):\n",
    "\n",
    "    Ns = np.size(xx) # Ns: # of spatial locs\n",
    "    Nt = np.size(gamgdt) # Nt: # time points\n",
    "\n",
    "    a = np.ones((Ns,))\n",
    "    b = np.zeros((Ns,))\n",
    "\n",
    "    Qc = np.cos(np.pi / 4)\n",
    "    Qs = np.sin(np.pi / 4)\n",
    "\n",
    "    for ii in range(Nt):\n",
    "\n",
    "        # apply phase accural\n",
    "        z = np.exp(-1j * (xx * gamgdt[ii]))\n",
    "        b = b * z\n",
    "\n",
    "        # apply rf\n",
    "        C = np.cos(np.abs(rf[ii]) / 2)\n",
    "        S = 1j * np.exp(1j * np.angle(rf[ii])) * np.sin(np.abs(rf[ii]) / 2)\n",
    "        at = a * C - b * np.conj(S)\n",
    "        bt = a * S + b * C\n",
    "        a = at\n",
    "        b = bt\n",
    "\n",
    "    z = np.exp(1j / 2 * (xx * np.sum(gamgdt)))\n",
    "    a = a * z\n",
    "    b = b * z\n",
    "\n",
    "    return a, b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 128 # number of frequency and b1 points to simulate \n",
    "\n",
    "xx = dtb1sel * 2 * np.pi * np.arange(-1000, 1000, 2000 / N).reshape((N,))\n",
    "xx = np.array([0])\n",
    "omOffRes = np.ones((rfb1sel.size,))\n",
    "\n",
    "MxyOrig = np.zeros((N, N), dtype=complex)\n",
    "\n",
    "b1 = np.arange(0, 2, 2 / 128) # 0-40 Gauss\n",
    "\n",
    "# run simulation, looping over b1 amplitudes\n",
    "for i in range(N):\n",
    "    a, b = rfSimCplx(dtb1sel * 2 * np.pi * 4258 * b1[i] * rfb1sel.flatten(), omOffRes, xx)\n",
    "    MxyOrig[:, i] = 2 * np.conj(a) * b \n",
    "\n",
    "#pl.ImagePlot(Mxy)\n",
    "\n",
    "#%%\n",
    "omOffRes = np.ones((rfflat.size,))\n",
    "MxyFlat = np.zeros((N, N), dtype=complex)\n",
    "\n",
    "# run simulation, looping over b1 amplitudes\n",
    "for i in range(N):\n",
    "    a, b = rfSimCplx(dtb1sel * 2 * np.pi * 4258 * b1[i] * rfflat.flatten(), omOffRes, xx)\n",
    "    MxyFlat[:, i] = 2 * np.conj(a) * b \n",
    "\n",
    "#pl.ImagePlot(Mxy)\n",
    "\n",
    "#%%\n",
    "omOffRes = np.ones((rfb1sel.size,))\n",
    "MxyClipped = np.zeros((N, N), dtype=complex)\n",
    "\n",
    "# run simulation, looping over b1 amplitudes\n",
    "for i in range(N):\n",
    "    a, b = rfSimCplx(dtb1sel * 2 * np.pi * 4258 * b1[i] * rfb1sel.flatten() / np.abs(rfb1sel.flatten()), omOffRes, xx)\n",
    "    MxyClipped[:, i] = 2 * np.conj(a) * b \n",
    "\n",
    "#pl.ImagePlot(Mxy)\n",
    "\n",
    "#%% \n",
    "fig, axs = plt.subplots(3, 1, layout='constrained')\n",
    "axs[0].plot(dtb1sel * 1000 * np.arange(rfb1sel.size), np.abs(rfb1sel.flatten()), label='Magnitude')\n",
    "axs[0].plot(dtb1sel * 1000 * np.arange(rfb1sel.size), np.real(rfb1sel.flatten()), label='Real')\n",
    "axs[0].plot(dtb1sel * 1000 * np.arange(rfb1sel.size), np.imag(rfb1sel.flatten()), label='Imag')\n",
    "axs[0].set_title('Original')\n",
    "axs[0].legend()\n",
    "axs[0].set_xlabel('ms')\n",
    "axs[1].plot(dtb1sel * 1000 * np.arange(rfflat.size), np.abs(rfflat), label='Magnitude')\n",
    "axs[1].plot(dtb1sel * 1000 * np.arange(rfflat.size), np.real(rfflat), label='Real')\n",
    "axs[1].plot(dtb1sel * 1000 * np.arange(rfflat.size), np.imag(rfflat), label='Imag')\n",
    "axs[1].set_title('Flattened')\n",
    "axs[1].legend()\n",
    "axs[1].set_xlabel('ms')\n",
    "axs[2].plot(b1 * 100, np.abs(MxyOrig[0,:]), label='Original')\n",
    "axs[2].plot(b1 * 100, np.abs(MxyFlat[0,:]), label='Flattened')\n",
    "axs[2].plot(b1 * 100, np.abs(MxyClipped[0,:]), label='Scaled')\n",
    "axs[2].legend()\n",
    "axs[2].set_xlabel('uT')\n",
    "axs[2].set_ylabel('|Mxy|')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.19"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}