#!/usr/bin/env python
# coding: utf-8

# # Audio Reproduction with Loudspeakers
# 
# 
# [return to main page](index.ipynb)
# 
# For the analysis of existing and the development of new reproduction methods,
# simulations are a very helpful tool.
# In the following exercises, we will simulate some sound fields by means of the
# [Sound Field Synthesis (SFS) Toolbox for Python](http://sfs.readthedocs.org/).
# 
# These simulations are assuming *free-field* conditions, i.e. the simulated
# loudspeakers are not located in a conventional room but in an infinitely large
# volume of air.
# In addition, the loudspeakers are modeled as idealized point sources which radiate
# uniformly in all directions and for all frequencies.

# ## Preparations
# 
# If it's not installed already, you can install the SFS module with:
# 
#     python3 -m pip install sfs --user
# 
# If you have only Python 3 installed on your system, you may have to use `python` instead of `python3`.
# 
# Afterwards, you should re-start any running IPython kernels.
# 
# Once installed, we can import it into our Python session:

# In[ ]:


import sfs


# While we're at it, let's also import some other stuff and enable plotting:

# In[ ]:


import matplotlib.pyplot as plt
import numpy as np


# We will simulate the sound pressure at different points in space.
# To specify the area we are interested in, we create a *grid*:

# In[ ]:


grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.02)


# 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.
# 
# *Exercise:* What does the third argument mean in our case?
# How many dimensions does our grid have?

# In[ ]:





# ## Sound Sources
# 
# Before we start analyzing loudspeaker systems, let's see what kinds of sound
# sources we can simulate.

# ### Point Source
# 
# 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.

# In[ ]:


x0 = 0, 1.5, 0
f = 1000  # Hz
omega = 2 * np.pi * f  #rad/s


# In[ ]:


p_point = sfs.fd.source.point(omega, x0, grid)
sfs.plot2d.amplitude(p_point, grid)
plt.title("Point Source at {} m".format(x0))


# The amplitude of the sound field is a bit weak ...
# 
# *Exercise:* Multiply the sound pressure field by a scaling factor of $4\pi$ to get an appropriate amplitude.

# In[ ]:


scaling_factor_point_source = 4 * np.pi


# In[ ]:





# *Exercise:* Try different source positions and different frequencies.

# In[ ]:





# ### Line Source
# 
# 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.

# In[ ]:


x0 = 0, 1.5
f = 1000  # Hz
omega = 2 * np.pi * f  # rad/s


# In[ ]:


p_line = sfs.fd.source.line(omega, x0, grid)
sfs.plot2d.amplitude(p_line, grid)
plt.title("Line Source at {} m".format(x0[:2]));


# Again, the amplitude is a bit weak, let's scale it up!
# This time, the scaling factor is a bit more involved:

# In[ ]:


scaling_factor_line_source = np.sqrt(8 * np.pi * omega / sfs.default.c) * np.exp(1j * np.pi / 4)


# BTW, you can get (and set) the speed of sound currently used by the SFS toolbox via the variable `sfs.default.c`.
# 
# 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.
# 
# *Exercise:* Scale the sound field by the given factor.

# In[ ]:





# *Exercise:* Again, try different source positions and different frequencies.

# In[ ]:





# *Exercise:* What's the difference between the sound fields of a point source and a line source?

# In[ ]:





# ### Plane Wave
# 
# 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.

# In[ ]:


x0 = 0, 1.5, 0
n0 = 0, -1, 0
f = 1000  # Hz
omega = 2 * np.pi * f  # rad/s


# In[ ]:


p_plane = sfs.fd.source.plane(omega, x0, n0, grid)
sfs.plot2d.amplitude(p_plane, grid);
plt.title("Plane wave with $n_0 = {}$".format(n0));


# This time, we don't need to scale the sound field.
# 
# *Exercise:* How can you see that the plane wave in the plot travels down and not up?

# In[ ]:





# *Exercise:* Try different propagation angles and different frequencies.

# In[ ]:





# *Exercise:* Compared to point source and line source, how does the level of the plane wave decay over distance?

# In[ ]:





# ## Two-channel Stereophony
# 
# As a first reproduction method, we'll have a look at *two-channel stereophony*.

# In[ ]:


def plot_stereo(f, weights=[1, 1]):
    """Plot a monochromatic stereo sound field.

    f: frequency in Hz
    weights: pair of weighting factors for the loudspeakers (can be real or complex)

    """
    omega = 2 * np.pi * f  # rad/s
    scaling_factor_point_source = 4 * np.pi
    weight_l, weight_r = np.asarray(weights) * scaling_factor_point_source

    r, phi = np.sqrt(3), 30
    x1 = -r * np.sin(phi*np.pi/180), r * np.cos(phi*np.pi/180), 0
    x2 = +r * np.sin(phi*np.pi/180), r * np.cos(phi*np.pi/180), 0

    p_1 = weight_l * sfs.fd.source.point(omega, x1, grid)
    p_2 = weight_r * sfs.fd.source.point(omega, x2, grid)

    sfs.plot2d.amplitude(p_1 + p_2, grid)


# In[ ]:


plot_stereo(1000)


# *Exercise:* Try different frequencies.

# In[ ]:





# *Exercise:* Compare the sound fields with the sound field of a single point source.
# What effects can be seen in the stereo sound field?

# In[ ]:





# *Exercise:* Where is the *phantom source*?

# In[ ]:





# *Exercise:* Move the phantom source using *intensity stereophony*.
# Try level differences of, say, 6 dB and 12 dB.

# In[ ]:





# *Exercise:* Try to move the phantom source using *time-of-arrival stereophony* (using phase differences between complex weighting factors).

# In[ ]:





# ## A Linear Loudspeaker Array
# 
# For the following exercises, we need a [loudspeaker array](http://sfs.readthedocs.org/en/latest/#module-sfs.array).

# In[ ]:


x0, n0, a0 = sfs.array.linear(20, 0.15, center=[0, 1, 0],
                              orientation=[0, -1, 0])


# ## Wave Field Synthesis, Point Source
# 
# Let's write a little function to plot the sound field of a virtual point source reproduced by WFS:

# In[ ]:


def plot_wfs_point_source(f, xs):
    """Plot a point source using Wave Field Synthesis.
    
    f: frequency in Hz
    xs: position vector of the virtual source
    
    """
    omega = 2 * np.pi * f  # rad/s
    d, selection, secondary_source = sfs.fd.wfs.point_25d(omega, x0, n0, xs)
    normalize_gain = 4 * np.pi
    
    twin = sfs.tapering.tukey(selection, alpha=.3)
    p = sfs.fd.synthesize(d, twin, [x0, n0, a0], secondary_source, grid=grid)
    sfs.plot2d.amplitude(normalize_gain * p, grid)
    sfs.plot2d.loudspeakers(x0, n0, selection * a0, size=0.135)


# In[ ]:


plot_wfs_point_source(1000, [0, 1.5, 0])


# *Exercise:* Compare the plot with the sound field of an ideal point source.

# In[ ]:





# *Exercise:* Increase the frequency to 2000 Hz.
# Suddenly, the sound field doesn't look like that of a point source anymore.
# What's the reason for the differences?

# In[ ]:





# *Exercise:* Try to empirically find the frequency where those artifacts start to appear.

# In[ ]:





# *Exercise:* Increase the loudspeaker distance from 15 to 20 cm.
# At which frequency do the artifacts show up now?

# In[ ]:





# *Exercise:* What does *2.5D* mean?

# In[ ]:





# ## WFS, Plane Wave
# 
# Now let's create a plane wave instead of a point source ...
# 
# Don't forget to set the loudspeaker distance to 15 cm again:

# In[ ]:


x0, n0, a0 = sfs.array.linear(20, 0.15, center=[0, 1, 0],
                              orientation=[0, -1, 0])


# In[ ]:


def plot_wfs_plane_wave(f, npw):
    """Plot a plane wave using Wave Field Synthesis.
    
    f: frequency in Hz
    npw: vector with the propagation direction of the virtual source
    
    """
    omega = 2 * np.pi * f  # rad/s
    
    d, selection, secondary_source = sfs.fd.wfs.plane_25d(omega, x0, n0, npw)
    
    twin = sfs.tapering.tukey(selection, alpha=.3)
    p = sfs.fd.synthesize(d, twin, [x0, n0, a0], secondary_source, grid=grid)
    sfs.plot2d.amplitude(p, grid)
    sfs.plot2d.loudspeakers(x0, n0, selection * a0, size=0.135)


# In[ ]:


plot_wfs_plane_wave(1000, [0, -1, 0])


# *Exercise:* Compare the plot with the sound field of an ideal plane wave.
# Try different propagation angles.

# In[ ]:





# *Exercise:* Change the frequency to 2000 Hz (and try some other frequencies, too).

# In[ ]:





# *Exercise:* Again, change the number (and spacing) of loudspeakers and note what's changing in the sound field.

# In[ ]:





# ## A Circular Loudspeaker Array
# 
# That's easy, just have a look at the [docs](http://sfs.readthedocs.org/#sfs.array.circular).

# In[ ]:


x0, n0, a0 = sfs.array.circular(40, 1)


# ## Again, WFS Point Source

# In[ ]:


plot_wfs_point_source(1000, [0, 1.5, 0])


# *Exercise:* Same as above: different frequencies, different distances between sources ...

# In[ ]:





# ## Again, WFS Plane Wave

# In[ ]:


plot_wfs_plane_wave(1000, [0, -1, 0])


# *Exercise:* Same as always ...

# In[ ]:





# ## Higher-Order Ambisonics
# 
# Wave Field Synthesis isn't the only sound field synthesis technique ...
# 
# Let's try nearfield-corrected higher-order Ambisonics (NFC-HOA)!
# 
# Note: NFC-HOA cannot be used with linear arrays, it only works with circular and spherical arrays.

# In[ ]:


def plot_hoa_plane_wave(f, npw, R):
    """Plot a plane wave using Higher-Order Ambisonics.
    
    f: frequency in Hz
    npw: vector with the propagation direction of the virtual source
    R: radius of the loudspeaker array in metres
    
    """
    omega = 2 * np.pi * f  # rad/s
       
    d, selection, secondary_source = sfs.fd.nfchoa.plane_25d(omega, x0, R, npw)
    p = sfs.fd.synthesize(d, selection, [x0, n0, a0], secondary_source, grid=grid)
    sfs.plot2d.amplitude(p, grid)
    sfs.plot2d.loudspeakers(x0, n0, selection * a0, size=0.15)


# In[ ]:


plot_hoa_plane_wave(1000, [0, -1, 0], 1)


# *Exercise:* Same ...

# In[ ]:





# *Exercise:* WFS vs. NFC-HOA: How do the artifacts for higher frequencies differ?

# In[ ]:





# Finally, let's try to reproduce a virtual point source with NFC-HOA:

# In[ ]:


def plot_hoa_point_source(f, xs, R):
    """Plot a point source using Higher-Order Ambisonics.
    
    f: frequency in Hz
    xs: position vector of the virtual source
    R: radius of the loudspeaker array in metres
    
    """
    omega = 2 * np.pi * f  # rad/s
    
    d, selection, secondary_source = sfs.fd.nfchoa.point_25d(omega, x0, R, xs)
    p = sfs.fd.synthesize(d, selection, [x0, n0, a0], secondary_source, grid=grid)
    p *= 4 * np.pi  # ad hoc scaling factor
    sfs.plot2d.amplitude(p, grid)
    sfs.plot2d.loudspeakers(x0, n0, selection * a0, size=0.135)


# In[ ]:


x0, n0, a0 = sfs.array.circular(40, 1)


# In[ ]:


plot_hoa_point_source(1000, [0, 1.5, 0], 1)


# *Exercise:* As always ...

# In[ ]:





# ## Solutions
# 
# If you had problems solving some of the exercises, don't despair!
# Have a look at the [example solutions](reproduction-solutions.ipynb).
# 
# <p xmlns:dct="http://purl.org/dc/terms/">
#   <a rel="license"
#      href="http://creativecommons.org/publicdomain/zero/1.0/">
#     <img src="http://i.creativecommons.org/p/zero/1.0/88x31.png" style="border-style: none;" alt="CC0" />
#   </a>
#   <br />
#   To the extent possible under law,
#   <span rel="dct:publisher" resource="[_:publisher]">the person who associated CC0</span>
#   with this work has waived all copyright and related or neighboring
#   rights to this work.
# </p>