{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Audio Reproduction with Loudspeakers\n", "\n", "\n", "[return to main page](index.ipynb)\n", "\n", "For the analysis of existing and the development of new reproduction methods,\n", "simulations are a very helpful tool.\n", "In the following exercises, we will simulate some sound fields by means of the\n", "[Sound Field Synthesis (SFS) Toolbox for Python](http://sfs.readthedocs.org/).\n", "\n", "These simulations are assuming *free-field* conditions, i.e. the simulated\n", "loudspeakers are not located in a conventional room but in an infinitely large\n", "volume of air.\n", "In addition, the loudspeakers are modeled as idealized point sources which radiate\n", "uniformly in all directions and for all frequencies." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "## Preparations\n", "\n", "If it's not installed already, you can install the SFS module with:\n", "\n", " python3 -m pip install sfs --user\n", "\n", "If you have only Python 3 installed on your system, you may have to use `python` instead of `python3`.\n", "\n", "Afterwards, you should re-start any running IPython kernels.\n", "\n", "Once installed, we can import it into our Python session:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sfs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While we're at it, let's also import some other stuff and enable plotting:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will simulate the sound pressure at different points in space.\n", "To specify the area we are interested in, we create a *grid*:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Have a look at the [documentation](https://sfs-python.readthedocs.io/en/0.6.2/sfs.util.html?highlight=grid#sfs.util.xyz_grid) to find out what the function parameters mean.\n", "\n", "*Exercise:* What does the third argument mean in our case?\n", "How many dimensions does our grid have?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Sound Sources\n", "\n", "Before we start analyzing loudspeaker systems, let's see what kinds of sound\n", "sources we can simulate." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Point Source\n", "\n", "Let's plot a [point source](http://sfs.readthedocs.org/#sfs.mono.source.point) at the position $(0, 1.5, 0)$ metres with a frequency of 1000 Hertz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = 0, 1.5, 0\n", "f = 1000 # Hz\n", "omega = 2 * np.pi * f #rad/s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p_point = sfs.fd.source.point(omega, x0, grid)\n", "sfs.plot2d.amplitude(p_point, grid)\n", "plt.title(\"Point Source at {} m\".format(x0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The amplitude of the sound field is a bit weak ...\n", "\n", "*Exercise:* Multiply the sound pressure field by a scaling factor of $4\\pi$ to get an appropriate amplitude." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scaling_factor_point_source = 4 * np.pi" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Try different source positions and different frequencies." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Line Source\n", "\n", "Let's plot a [line source](https://sfs-python.readthedocs.io/en/0.6.2/sfs.fd.source.html#sfs.fd.source.line) (parallel to the z-axis) at the position $(0, 1.5)$ metres with a frequency of 1000 Hertz." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = 0, 1.5\n", "f = 1000 # Hz\n", "omega = 2 * np.pi * f # rad/s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p_line = sfs.fd.source.line(omega, x0, grid)\n", "sfs.plot2d.amplitude(p_line, grid)\n", "plt.title(\"Line Source at {} m\".format(x0[:2]));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Again, the amplitude is a bit weak, let's scale it up!\n", "This time, the scaling factor is a bit more involved:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "scaling_factor_line_source = np.sqrt(8 * np.pi * omega / sfs.default.c) * np.exp(1j * np.pi / 4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "BTW, you can get (and set) the speed of sound currently used by the SFS toolbox via the variable `sfs.default.c`.\n", "\n", "Don't worry too much about this strange scaling factor, just multiply the sound field of the line source with it and then you're done.\n", "\n", "*Exercise:* Scale the sound field by the given factor." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Again, try different source positions and different frequencies." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* What's the difference between the sound fields of a point source and a line source?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plane Wave\n", "\n", "Let's plot a [plane wave](http://sfs.readthedocs.org/#sfs.mono.source.plane) with a frequency of 1000 Hertz which propagates in the direction of the negative y-axis." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0 = 0, 1.5, 0\n", "n0 = 0, -1, 0\n", "f = 1000 # Hz\n", "omega = 2 * np.pi * f # rad/s" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "p_plane = sfs.fd.source.plane(omega, x0, n0, grid)\n", "sfs.plot2d.amplitude(p_plane, grid);\n", "plt.title(\"Plane wave with $n_0 = {}$\".format(n0));" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This time, we don't need to scale the sound field.\n", "\n", "*Exercise:* How can you see that the plane wave in the plot travels down and not up?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Try different propagation angles and different frequencies." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Compared to point source and line source, how does the level of the plane wave decay over distance?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Two-channel Stereophony\n", "\n", "As a first reproduction method, we'll have a look at *two-channel stereophony*." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_stereo(f, weights=[1, 1]):\n", " \"\"\"Plot a monochromatic stereo sound field.\n", "\n", " f: frequency in Hz\n", " weights: pair of weighting factors for the loudspeakers (can be real or complex)\n", "\n", " \"\"\"\n", " omega = 2 * np.pi * f # rad/s\n", " scaling_factor_point_source = 4 * np.pi\n", " weight_l, weight_r = np.asarray(weights) * scaling_factor_point_source\n", "\n", " r, phi = np.sqrt(3), 30\n", " x1 = -r * np.sin(phi*np.pi/180), r * np.cos(phi*np.pi/180), 0\n", " x2 = +r * np.sin(phi*np.pi/180), r * np.cos(phi*np.pi/180), 0\n", "\n", " p_1 = weight_l * sfs.fd.source.point(omega, x1, grid)\n", " p_2 = weight_r * sfs.fd.source.point(omega, x2, grid)\n", "\n", " sfs.plot2d.amplitude(p_1 + p_2, grid)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_stereo(1000)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Try different frequencies." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Compare the sound fields with the sound field of a single point source.\n", "What effects can be seen in the stereo sound field?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Where is the *phantom source*?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Move the phantom source using *intensity stereophony*.\n", "Try level differences of, say, 6 dB and 12 dB." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Try to move the phantom source using *time-of-arrival stereophony* (using phase differences between complex weighting factors)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A Linear Loudspeaker Array\n", "\n", "For the following exercises, we need a [loudspeaker array](http://sfs.readthedocs.org/en/latest/#module-sfs.array)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0, n0, a0 = sfs.array.linear(20, 0.15, center=[0, 1, 0],\n", " orientation=[0, -1, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Wave Field Synthesis, Point Source\n", "\n", "Let's write a little function to plot the sound field of a virtual point source reproduced by WFS:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_wfs_point_source(f, xs):\n", " \"\"\"Plot a point source using Wave Field Synthesis.\n", " \n", " f: frequency in Hz\n", " xs: position vector of the virtual source\n", " \n", " \"\"\"\n", " omega = 2 * np.pi * f # rad/s\n", " d, selection, secondary_source = sfs.fd.wfs.point_25d(omega, x0, n0, xs)\n", " normalize_gain = 4 * np.pi\n", " \n", " twin = sfs.tapering.tukey(selection, alpha=.3)\n", " p = sfs.fd.synthesize(d, twin, [x0, n0, a0], secondary_source, grid=grid)\n", " sfs.plot2d.amplitude(normalize_gain * p, grid)\n", " sfs.plot2d.loudspeakers(x0, n0, selection * a0, size=0.135)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_wfs_point_source(1000, [0, 1.5, 0])" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "*Exercise:* Compare the plot with the sound field of an ideal point source." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Increase the frequency to 2000 Hz.\n", "Suddenly, the sound field doesn't look like that of a point source anymore.\n", "What's the reason for the differences?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Try to empirically find the frequency where those artifacts start to appear." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Increase the loudspeaker distance from 15 to 20 cm.\n", "At which frequency do the artifacts show up now?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* What does *2.5D* mean?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## WFS, Plane Wave\n", "\n", "Now let's create a plane wave instead of a point source ...\n", "\n", "Don't forget to set the loudspeaker distance to 15 cm again:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0, n0, a0 = sfs.array.linear(20, 0.15, center=[0, 1, 0],\n", " orientation=[0, -1, 0])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_wfs_plane_wave(f, npw):\n", " \"\"\"Plot a plane wave using Wave Field Synthesis.\n", " \n", " f: frequency in Hz\n", " npw: vector with the propagation direction of the virtual source\n", " \n", " \"\"\"\n", " omega = 2 * np.pi * f # rad/s\n", " \n", " d, selection, secondary_source = sfs.fd.wfs.plane_25d(omega, x0, n0, npw)\n", " \n", " twin = sfs.tapering.tukey(selection, alpha=.3)\n", " p = sfs.fd.synthesize(d, twin, [x0, n0, a0], secondary_source, grid=grid)\n", " sfs.plot2d.amplitude(p, grid)\n", " sfs.plot2d.loudspeakers(x0, n0, selection * a0, size=0.135)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_wfs_plane_wave(1000, [0, -1, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Compare the plot with the sound field of an ideal plane wave.\n", "Try different propagation angles." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Change the frequency to 2000 Hz (and try some other frequencies, too)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Again, change the number (and spacing) of loudspeakers and note what's changing in the sound field." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## A Circular Loudspeaker Array\n", "\n", "That's easy, just have a look at the [docs](http://sfs.readthedocs.org/#sfs.array.circular)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0, n0, a0 = sfs.array.circular(40, 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Again, WFS Point Source" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_wfs_point_source(1000, [0, 1.5, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Same as above: different frequencies, different distances between sources ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Again, WFS Plane Wave" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_wfs_plane_wave(1000, [0, -1, 0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Same as always ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Higher-Order Ambisonics\n", "\n", "Wave Field Synthesis isn't the only sound field synthesis technique ...\n", "\n", "Let's try nearfield-corrected higher-order Ambisonics (NFC-HOA)!\n", "\n", "Note: NFC-HOA cannot be used with linear arrays, it only works with circular and spherical arrays." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_hoa_plane_wave(f, npw, R):\n", " \"\"\"Plot a plane wave using Higher-Order Ambisonics.\n", " \n", " f: frequency in Hz\n", " npw: vector with the propagation direction of the virtual source\n", " R: radius of the loudspeaker array in metres\n", " \n", " \"\"\"\n", " omega = 2 * np.pi * f # rad/s\n", " \n", " d, selection, secondary_source = sfs.fd.nfchoa.plane_25d(omega, x0, R, npw)\n", " p = sfs.fd.synthesize(d, selection, [x0, n0, a0], secondary_source, grid=grid)\n", " sfs.plot2d.amplitude(p, grid)\n", " sfs.plot2d.loudspeakers(x0, n0, selection * a0, size=0.15)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_hoa_plane_wave(1000, [0, -1, 0], 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* Same ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* WFS vs. NFC-HOA: How do the artifacts for higher frequencies differ?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, let's try to reproduce a virtual point source with NFC-HOA:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def plot_hoa_point_source(f, xs, R):\n", " \"\"\"Plot a point source using Higher-Order Ambisonics.\n", " \n", " f: frequency in Hz\n", " xs: position vector of the virtual source\n", " R: radius of the loudspeaker array in metres\n", " \n", " \"\"\"\n", " omega = 2 * np.pi * f # rad/s\n", " \n", " d, selection, secondary_source = sfs.fd.nfchoa.point_25d(omega, x0, R, xs)\n", " p = sfs.fd.synthesize(d, selection, [x0, n0, a0], secondary_source, grid=grid)\n", " p *= 4 * np.pi # ad hoc scaling factor\n", " sfs.plot2d.amplitude(p, grid)\n", " sfs.plot2d.loudspeakers(x0, n0, selection * a0, size=0.135)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "x0, n0, a0 = sfs.array.circular(40, 1)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plot_hoa_point_source(1000, [0, 1.5, 0], 1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "*Exercise:* As always ..." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solutions\n", "\n", "If you had problems solving some of the exercises, don't despair!\n", "Have a look at the [example solutions](reproduction-solutions.ipynb).\n", "\n", "

\n", " \n", " \"CC0\"\n", " \n", "
\n", " To the extent possible under law,\n", " the person who associated CC0\n", " with this work has waived all copyright and related or neighboring\n", " rights to this work.\n", "

" ] } ], "metadata": { "kernelspec": { "display_name": "mycomac", "language": "python", "name": "mycomac" }, "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.12" } }, "nbformat": 4, "nbformat_minor": 1 }