{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import spharpy\n",
    "from pyfar import Coordinates\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colorbar import Colorbar\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plane Wave Decomposition\n",
    "\n",
    "\n",
    "\n",
    "The sound pressure at the microphones of a spherical array resulting from a plane wave incident from $(\\theta_i, \\phi_i)$ may be written as the expansion of spherical harmonic coefficients as\n",
    "\n",
    "$$ \\mathbf{p}(kr) = \\mathbf{Y} \\mathbf{B}  \\mathbf{a}(k, \\theta, \\phi)$$\n",
    "\n",
    "with the plane wave density function\n",
    "\n",
    "$$ \\mathbf{a}(k, \\theta, \\phi) = [{Y_0^0}^*(\\theta_i, \\phi_i), {Y_1^{-1}}^*(\\theta_i, \\phi_i), \\dots, {Y_N^N}^*(\\theta_i, \\phi_i)] $$\n",
    "\n",
    "for the single plane wave with unit amplitude, and the modal strength matrix \n",
    "\n",
    "$$ \\mathbf{B} = \\mathrm{diag}([b_0(kr), b_1(kr), b_1(kr), \\dots, b_N(kr)]) $$\n",
    "\n",
    "representing the type of sphere, which may either be an open sphere (microphones mounted on open frame support) or a rigid sphere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_max = 4\n",
    "kr = np.linspace(0.5, 10, 512)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "B_open = spharpy.spherical.modal_strength(n_max, kr, arraytype='open')\n",
    "B_rigid = spharpy.spherical.modal_strength(n_max, kr, arraytype='rigid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16, 4))\n",
    "ax = plt.subplot(1, 2, 1)\n",
    "n_max = 4\n",
    "mask = spharpy.spherical.nm_to_acn(np.arange(0, n_max+1), np.zeros(n_max+1))\n",
    "ax.plot(kr, 10*np.log10(np.abs(np.diagonal(B_open, axis1=1, axis2=2)[:, mask])))\n",
    "plt.grid(True)\n",
    "ax.set_ylabel('Magnitude in dB')\n",
    "ax.set_xlabel('kr')\n",
    "plt.legend(['$b_0(kr)$', '$b_1(kr)$', '$b_2(kr)$', '$b_3(kr)$', '$b_4(kr)$'])\n",
    "ax.set_title('Open Sphere')\n",
    "ax.set_ylim((-20, 15))\n",
    "ax = plt.subplot(1, 2, 2)\n",
    "mask = spharpy.spherical.nm_to_acn(np.arange(0, n_max+1), np.zeros(n_max+1))\n",
    "ax.plot(kr, 10*np.log10(np.abs(np.diagonal(B_rigid, axis1=1, axis2=2)[:, mask])))\n",
    "plt.grid(True)\n",
    "ax.set_ylabel('Magnitude in dB')\n",
    "ax.set_xlabel('kr')\n",
    "plt.legend(['$b_0(kr)$', '$b_1(kr)$', '$b_2(kr)$', '$b_3(kr)$', '$b_4(kr)$'])\n",
    "ax.set_title('Rigid Sphere')\n",
    "ax.set_ylim((-20, 15));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kr = np.array([2.5])\n",
    "arraytype = 'rigid'\n",
    "doa_pw = Coordinates(0, 1, 0)\n",
    "plane_wave_density = spharpy.spherical.spherical_harmonic_basis(n_max, doa_pw)\n",
    "\n",
    "B = np.squeeze(spharpy.spherical.modal_strength(n_max, kr, arraytype=arraytype))\n",
    "\n",
    "p_nm = B @ plane_wave_density.T.conj()\n",
    "\n",
    "sphere = spharpy.samplings.hyperinterpolation(n_max=30)\n",
    "Y_sphere = spharpy.spherical.spherical_harmonic_basis(n_max, sphere)\n",
    "\n",
    "p_sphere =  np.squeeze(Y_sphere @ p_nm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7, 7))\n",
    "spharpy.plot.balloon(sphere, p_sphere)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decomposition of the Sound Field\n",
    "\n",
    "Before the actual plane wave decomposition, the sound pressure recorded at the microphone positions needs to be decomposed into spherical harmonic coefficients using the spherical harmonic transform\n",
    "\n",
    "$$ \\mathbf{p_{nm}} = \\mathbf{Y}^+ \\mathbf{p}(kr) $$\n",
    "\n",
    "We can then approximate the plane wave density function using the matrix vector product\n",
    "\n",
    "$$ \\mathbf{a}(k) = \\frac{4\\pi}{(N+1)^2}\\mathbf{Y_s}^H \\mathbf{B}^{-1} \\mathbf{p_{nm}}$$\n",
    "\n",
    "with the steering matrix $\\mathbf{Y_s}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "steering_directions = spharpy.samplings.gaussian(30)\n",
    "Y_steering = spharpy.spherical.spherical_harmonic_basis(n_max, steering_directions)\n",
    "plane_wave_density = 4*np.pi/(n_max+1)**2 * np.squeeze(np.real(Y_steering @ np.linalg.pinv(B) @ p_nm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 4))\n",
    "spharpy.plot.contour_map(steering_directions, np.abs(plane_wave_density), projection='mollweide')"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "9131b587d0f7e98eda94de94aac65297fc22abcfb5aba426e123fc3ffc316114"
  },
  "kernelspec": {
   "display_name": "Python 3.7.2 64-bit ('spharpy': conda)",
   "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.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}