simplify referencing code, add reference, add difference plot, small text mods

Author: fs446

doc/examples/wfs-referencing.ipynb

@@ -6,8 +6,7 @@
    "source": [
     "# WFS Referencing Schemes\n",
     "\n",
-    "Illustrates the usage of the SFS toolbox for the simulation of different sound fields using 2.5D WFS and referencing schemes for contours that exhibit amplitude correct synthesis, cf. Ch. 4.1.3 in Gergely Firtha's doctoral thesis [A Generalized Wave Field Synthesis Framework with Application for Moving Virtual Sources](https://last.hit.bme.hu/download/firtha/PhD_thesis/firtha_phd_thesis.pdf) of 2019.\n",
-    "\n"
+    "Illustrates the usage of the SFS toolbox for the simulation of different 2.5D WFS referencing schemes for contours within the listening area that exhibit amplitude correct synthesis, cf. Ch. 4.1.3 in [[Fir19]](https://last.hit.bme.hu/download/firtha/PhD_thesis/firtha_phd_thesis.pdf) and [[FFSS17]](doi:10.1109/TASLP.2017.2689245)"
    ]
   },
   {
@@ -31,7 +30,7 @@
     "array = sfs.array.circular(N=64, R=R)\n",
     "grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02)\n",
     "\n",
-    "xs = -4, 0, 0  # point source on negative x-axis\n",
+    "xs = sfs.util.asarray_of_rows((-4, 0, 0))  # point source on negative x-axis\n",
     "normalize_gain = 4 * np.pi * np.linalg.norm(xs)\n",
     "wavelength = 1 / 4  # m\n",
     "f = sfs.default.c / wavelength  # Hz\n",
@@ -63,9 +62,10 @@
     "                                tapering_window,\n",
     "                                size=0.125, ax=axs[i])\n",
     "        axs[i].plot(xref[:, 0][selection],\n",
-    "                    xref[:, 1][selection], 'C5o', ms=4)\n",
+    "                    xref[:, 1][selection], 'o', color='dimgray', ms=4)\n",
     "        axs[i].grid(True)\n",
-    "    plt.tight_layout()"
+    "    plt.tight_layout()\n",
+    "    return p\n"
    ]
   },
   {
@@ -74,7 +74,7 @@
    "source": [
     "### Line as reference contour\n",
     "\n",
-    "The reference contour is calculated according to https://github.com/spatialaudio/wfs_chapter_hda/blob/master/python/wfs25d_circSSD.py#L91 for a virtual point source on x-axis."
+    "The reference contour is calculated according to [Fir19, eq. (52)] for a virtual point source on x-axis."
    ]
   },
   {
@@ -85,28 +85,16 @@
    "source": [
     "# reference contour is a straight line\n",
     "xref_line = 0\n",
-    "# calc reference contour xref(x0):\n",
-    "x0_tmp = array.x.T[np.newaxis, :]\n",
-    "xs_tmp = np.array(xs)[np.newaxis, :, np.newaxis]\n",
-    "x0xs = x0_tmp - xs_tmp\n",
-    "x0xs_length = np.linalg.norm(x0xs, axis=1)\n",
-    "x0xs_unit = x0xs / x0xs_length\n",
-    "n0_ = array.n.T[np.newaxis, :]\n",
-    "xref = np.zeros_like(x0_tmp)\n",
-    "# code assumes that virtual point source is on x-axis:\n",
-    "for i in range(array.x.shape[0]):\n",
-    "    cosbeta = np.dot(-n0_[0, :, i], [-1, 0, 0])  # use outward SSD normal\n",
-    "    tmp = x0xs_unit[0, :, i]\n",
-    "    tmp *= -x0xs_length[0, i]\n",
-    "    tmp *= xref_line + R * cosbeta\n",
-    "    tmp /= xs_tmp[0, 0, 0] + R * cosbeta\n",
-    "    xref[0, :, i] = x0_tmp[0, :, i] + tmp\n",
-    "xref = np.squeeze(xref).T\n",
+    "# calc reference contour xref(x0), cf. [Fir19, eq. (24), (31), (52)]\n",
+    "# this code assumes virtual point source on negative x-axis\n",
+    "cosbeta = (array.n @ [1, 0, 0]).reshape(-1, 1)\n",
+    "xref = array.x + (xs - array.x) * \\\n",
+    "    (xref_line + R * cosbeta) / (xs[0, 0] + R * cosbeta)\n",
     "\n",
     "d, selection, secondary_source = sfs.fd.wfs.point_25d(\n",
     "    omega, array.x, array.n, xs, xref=xref)\n",
-    "sound_field(d * normalize_gain, xref, selection,\n",
-    "            secondary_source, array, grid, tapering=False)"
+    "p_line = sound_field(d * normalize_gain, xref, selection,\n",
+    "                     secondary_source, array, grid, tapering=False)"
    ]
   },
   {
@@ -115,7 +103,7 @@
    "source": [
     "### Circle as reference contour\n",
     "\n",
-    "This reference contour is a circle with origin xs and a radius, such that the origin is on this circle. This contour is straightforward with some obvious vector calculus."
+    "This reference contour is a circle with origin xs and a radius, such that the origin is on this circle. This contour is obtained with straightforward vector calculus."
    ]
   },
   {
@@ -126,19 +114,31 @@
    "source": [
     "# reference contour is a circle with origin xs\n",
     "xref_dist = np.linalg.norm(xs)\n",
-    "# calc reference contour xref(x0):\n",
-    "x0_tmp = array.x.T[np.newaxis, :]\n",
-    "xs_tmp = np.array(xs)[np.newaxis, :, np.newaxis]\n",
-    "x0xs = x0_tmp - xs_tmp\n",
-    "x0xs_length = np.linalg.norm(x0xs, axis=1)\n",
-    "x0xs_unit = x0xs / x0xs_length\n",
-    "xref = x0_tmp + (xref_dist - x0xs_length) * x0xs_unit\n",
-    "xref = np.squeeze(xref).T\n",
+    "# calc reference contour xref(x0), cf. [Fir19, eq. (24), (31)]\n",
+    "x0xs = array.x - xs\n",
+    "x0xs_length = np.linalg.norm(x0xs, axis=1, keepdims=True)\n",
+    "x0xs_unit = x0xs / x0xs_length  # unit vec\n",
+    "xref = xs + xref_dist * x0xs_unit\n",
     "\n",
     "d, selection, secondary_source = sfs.fd.wfs.point_25d(\n",
     "    omega, array.x, array.n, xs, xref=xref)\n",
-    "sound_field(d * normalize_gain, xref, selection,\n",
-    "            secondary_source, array, grid, tapering=False)"
+    "p_circ = sound_field(d * normalize_gain, xref, selection,\n",
+    "                     secondary_source, array, grid, tapering=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# (complex-valued) ratio between the two soundfields tells\n",
+    "# us about the difference of both approaches\n",
+    "sfs.plot2d.level(p_line / p_circ, grid, vmax=1, vmin=-1,\n",
+    "                 cmap='seismic', colorbar_kwargs={'label': 'dB'})\n",
+    "sfs.plot2d.loudspeakers(array.x, array.n, size=0.125)\n",
+    "plt.grid(True)\n",
+    "plt.tight_layout()"
    ]
   }
  ],

doc/references.bib

@@ -96,3 +96,10 @@ @phdthesis{Schultz2016
     year = {2016},
     doi = {10.18453/rosdok_id00001765}
 }
+@phdthesis{Firtha2019,
+    author = {Firtha, G.},
+    title = {{A Generalized Wave Field Synthesis Framework with Application
+        for Moving Virtual Sources}},
+    school = {Budapest University of Technology and Economics},
+    year = {2019}
+}

sfs/fd/wfs.py

@@ -205,7 +205,7 @@ def point_25d(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None):
 
     is implemented.
     The theoretical link of `point_25d()` and `point_25d_legacy()` was
-    introduced as *unified WFS framework* in :cite:`Firtha2017`.
+    introduced as *unified WFS framework* in :cite:`Firtha2017`, :cite:`Firtha2019`.
 
     Examples
     --------
@@ -292,7 +292,7 @@ def point_25d_legacy(omega, x0, n0, xs, xref=[0, 0, 0], c=None, omalias=None):
             \e{-\i\wc |\x_0-\x_\text{s}|}
 
     The theoretical link of `point_25d()` and `point_25d_legacy()` was
-    introduced as *unified WFS framework* in :cite:`Firtha2017`.
+    introduced as *unified WFS framework* in :cite:`Firtha2017`, :cite:`Firtha2019`.
     Also cf. Eq. (2.145)-(2.147) :cite:`Schultz2016`.
 
     Examples