From 4395c15f640096f73cf4b739f076c535fda04c40 Mon Sep 17 00:00:00 2001 From: Nara Hahn Date: Sat, 7 Jan 2017 01:34:55 +0100 Subject: [PATCH] Time-domain NFC-HOA --- doc/examples/time_domain_nfchoa.py | 50 ++++ sfs/mono/drivingfunction.py | 29 +- sfs/time/drivingfunction.py | 452 ++++++++++++++++++++++++++++- sfs/time/source.py | 3 +- sfs/util.py | 42 +++ 5 files changed, 559 insertions(+), 17 deletions(-) create mode 100644 doc/examples/time_domain_nfchoa.py diff --git a/doc/examples/time_domain_nfchoa.py b/doc/examples/time_domain_nfchoa.py new file mode 100644 index 0000000..6195b99 --- /dev/null +++ b/doc/examples/time_domain_nfchoa.py @@ -0,0 +1,50 @@ +"""Create some examples of time-domain NFC-HOA.""" + +import numpy as np +import matplotlib.pyplot as plt +import sfs +from scipy.signal import unit_impulse + +# Parameters +fs = 44100 # sampling frequency +grid = sfs.util.xyz_grid([-2, 2], [-2, 2], 0, spacing=0.005) +N = 60 # number of secondary sources +R = 1.5 # radius of circular array +array = sfs.array.circular(N, R) + +# Excitation signal +signal = unit_impulse(512), fs, 0 + +# Plane wave +max_order = None +npw = [0, -1, 0] # propagating direction +t = 0 # observation time +delay, weight, sos, phaseshift, selection, secondary_source = \ + sfs.time.drivingfunction.nfchoa_25d_plane(array.x, R, npw, fs, max_order) +d = sfs.time.drivingfunction.nfchoa_25d_driving_signals( + delay, weight, sos, phaseshift, signal) +p = sfs.time.synthesize(d, selection, array, secondary_source, + observation_time=t, grid=grid) + +plt.figure() +sfs.plot.level(p, grid) +sfs.plot.loudspeaker_2d(array.x, array.n) +sfs.plot.virtualsource_2d([0, 0], ns=npw, type='plane') +plt.savefig('impulse_pw_nfchoa_25d.png') + +# Point source +max_order = 100 +xs = [1.5, 1.5, 0] # position +t = np.linalg.norm(xs) / sfs.default.c # observation time +delay, weight, sos, phaseshift, selection, secondary_source = \ + sfs.time.drivingfunction.nfchoa_25d_point(array.x, R, xs, fs, max_order) +d = sfs.time.drivingfunction.nfchoa_25d_driving_signals( + delay, weight, sos, phaseshift, signal) +p = sfs.time.synthesize(d, selection, array, secondary_source, + observation_time=t, grid=grid) + +plt.figure() +sfs.plot.level(p, grid) +sfs.plot.loudspeaker_2d(array.x, array.n) +sfs.plot.virtualsource_2d(xs, type='point') +plt.savefig('impulse_ps_nfchoa_25d.png') diff --git a/sfs/mono/drivingfunction.py b/sfs/mono/drivingfunction.py index fd188ca..cd08f77 100644 --- a/sfs/mono/drivingfunction.py +++ b/sfs/mono/drivingfunction.py @@ -325,14 +325,16 @@ def nfchoa_2d_plane(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None): plot(d, selection, secondary_source) """ + if max_order is None: + max_order = util.max_order_circular_harmonics(len(x0)) + x0 = util.asarray_of_rows(x0) k = util.wavenumber(omega, c) n = util.normalize_vector(n) phi, _, r = util.cart2sph(*n) phi0 = util.cart2sph(*x0.T)[0] - M = _max_order_circular_harmonics(len(x0), max_order) d = 0 - for m in range(-M, M + 1): + for m in range(-max_order, max_order + 1): d += 1j**-m / hankel2(m, k * r0) * np.exp(1j * m * (phi0 - phi)) selection = util.source_selection_all(len(x0)) return -2j / (np.pi*r0) * d, selection, secondary_source_point(omega, c) @@ -361,16 +363,18 @@ def nfchoa_25d_point(omega, x0, r0, xs, max_order=None, c=None): plot(d, selection, secondary_source) """ + if max_order is None: + max_order = util.max_order_circular_harmonics(len(x0)) + x0 = util.asarray_of_rows(x0) k = util.wavenumber(omega, c) xs = util.asarray_1d(xs) phi, _, r = util.cart2sph(*xs) phi0 = util.cart2sph(*x0.T)[0] - M = _max_order_circular_harmonics(len(x0), max_order) - hr = util.spherical_hn2(range(0, M + 1), k * r) - hr0 = util.spherical_hn2(range(0, M + 1), k * r0) + hr = util.spherical_hn2(range(0, max_order + 1), k * r) + hr0 = util.spherical_hn2(range(0, max_order + 1), k * r0) d = 0 - for m in range(-M, M + 1): + for m in range(-max_order, max_order + 1): d += hr[abs(m)] / hr0[abs(m)] * np.exp(1j * m * (phi0 - phi)) selection = util.source_selection_all(len(x0)) return d / (2 * np.pi * r0), selection, secondary_source_point(omega, c) @@ -399,15 +403,17 @@ def nfchoa_25d_plane(omega, x0, r0, n=[0, 1, 0], max_order=None, c=None): plot(d, selection, secondary_source) """ + if max_order is None: + max_order = util.max_order_circular_harmonics(len(x0)) + x0 = util.asarray_of_rows(x0) k = util.wavenumber(omega, c) n = util.normalize_vector(n) phi, _, r = util.cart2sph(*n) phi0 = util.cart2sph(*x0.T)[0] - M = _max_order_circular_harmonics(len(x0), max_order) d = 0 - hn2 = util.spherical_hn2(range(0, M + 1), k * r0) - for m in range(-M, M + 1): + hn2 = util.spherical_hn2(range(0, max_order + 1), k * r0) + for m in range(-max_order, max_order + 1): d += (-1j)**abs(m) / (k * hn2[abs(m)]) * np.exp(1j * m * (phi0 - phi)) selection = util.source_selection_all(len(x0)) return 2*1j / r0 * d, selection, secondary_source_point(omega, c) @@ -782,8 +788,3 @@ def secondary_source(position, _, grid): return _source.line(omega, position, grid, c) return secondary_source - - -def _max_order_circular_harmonics(N, max_order): - """Compute order of 2D HOA.""" - return N // 2 if max_order is None else max_order diff --git a/sfs/time/drivingfunction.py b/sfs/time/drivingfunction.py index 40d0a77..2c942d5 100644 --- a/sfs/time/drivingfunction.py +++ b/sfs/time/drivingfunction.py @@ -7,8 +7,8 @@ import matplotlib.pyplot as plt import numpy as np - from scipy.signal import unit_impulse import sfs + from scipy.signal import unit_impulse # Plane wave npw = sfs.util.direction_vector(np.radians(-45)) @@ -25,7 +25,8 @@ tf = rf / sfs.default.c # time-of-arrival at origin # Impulsive excitation - signal = unit_impulse(512), 44100 + fs = 44100 + signal = unit_impulse(512), fs # Circular loudspeaker array N = 32 # number of loudspeakers @@ -43,6 +44,8 @@ def plot(d, selection, secondary_source, t=0): """ import numpy as np from numpy.core.umath_tests import inner1d # element-wise inner product +from scipy.signal import besselap, sosfilt, zpk2sos +from scipy.special import eval_legendre as legendre from .. import default from .. import util from . import source as _source @@ -357,3 +360,448 @@ def secondary_source(position, _, signal, observation_time, grid): return _source.point(position, signal, observation_time, grid, c=c) return secondary_source + + +def matchedz_zpk(s_zeros, s_poles, s_gain, fs): + """Matched-z transform of poles and zeros. + + Parameters + ---------- + s_zeros : array_like + Zeros in the Laplace domain. + s_poles : array_like + Poles in the Laplace domain. + s_gain : float + System gain in the Laplace domain. + fs : int + Sampling frequency in Hertz. + + Returns + ------- + z_zeros : numpy.ndarray + Zeros in the z-domain. + z_poles : numpy.ndarray + Poles in the z-domain. + z_gain : float + System gain in the z-domain. + + See Also + -------- + :func:`scipy.signal.bilinear_zpk` + + """ + z_zeros = np.exp(s_zeros / fs) + z_poles = np.exp(s_poles / fs) + omega = 1j * np.pi * fs + s_gain *= np.prod((omega - s_zeros) / (omega - s_poles) + * (-1 - z_poles) / (-1 - z_zeros)) + return z_zeros, z_poles, np.real(s_gain) + + +def nfchoa_25d_plane(x0, r0, npw, fs, max_order=None, c=None, + s2z=matchedz_zpk): + r"""Virtual plane wave by 2.5-dimensional NFC-HOA. + + .. math:: + + D(\phi_0, s) = + 2\e{\frac{s}{c}r_0} + \sum_{m=-M}^{M} + (-1)^m + \Big(\frac{s}{s-\frac{c}{r_0}\sigma_0}\Big)^\mu + \prod_{l=1}^{\nu} + \frac{s^2}{(s-\frac{c}{r_0}\sigma_l)^2+(\frac{c}{r_0}\omega_l)^2} + \e{\i m(\phi_0 - \phi_\text{pw})} + + The driving function is represented in the Laplace domain, + from which the recursive filters are designed. + :math:`\sigma_l + \i\omega_l` denotes the complex roots of + the reverse Bessel polynomial. + The number of second-order sections is + :math:`\nu = \big\lfloor\tfrac{|m|}{2}\big\rfloor`, + whereas the number of first-order section :math:`\mu` is either 0 or 1 + for even and odd :math:`|m|`, respectively. + + Parameters + ---------- + x0 : (N, 3) array_like + Sequence of secondary source positions. + r0 : float + Radius of the circular secondary source distribution. + npw : (3,) array_like + Unit vector (propagation direction) of plane wave. + fs : int + Sampling frequency in Hertz. + max_order : int, optional + Ambisonics order. + c : float, optional + Speed of sound in m/s. + s2z : callable, optional + Function transforming s-domain poles and zeros into z-domain, + e.g. :func:`matchedz_zpk`, :func:`scipy.signal.bilinear_zpk`. + + Returns + ------- + delay : float + Overall delay in seconds. + weight : float + Overall weight. + sos : list of numpy.ndarray + Second-order section filters :func:`scipy.signal.sosfilt`. + phaseshift : (N,) numpy.ndarray + Phase shift in radians. + + Examples + -------- + .. plot:: + :context: close-figs + + delay, weight, sos, phaseshift, selection, secondary_source = \ + sfs.time.drivingfunction.nfchoa_25d_plane(array.x, R, npw, fs) + d = sfs.time.drivingfunction.nfchoa_25d_driving_signals( + delay, weight, sos, phaseshift, signal) + plot(d, selection, secondary_source) + + """ + if max_order is None: + max_order = util.max_order_circular_harmonics(len(x0)) + if c is None: + c = default.c + + x0 = util.asarray_of_rows(x0) + npw = util.asarray_1d(npw) + phi0, _, _ = util.cart2sph(*x0.T) + phipw, _, _ = util.cart2sph(*npw) + phaseshift = phi0 - phipw + np.pi + + delay = -r0 / c + weight = 2 + sos = [] + for m in range(max_order + 1): + _, p, _ = besselap(m, norm='delay') + s_zeros = np.zeros(m) + s_poles = c / r0 * p + s_gain = 1 + z_zeros, z_poles, z_gain = s2z(s_zeros, s_poles, s_gain, fs) + sos.append(zpk2sos(z_zeros, z_poles, z_gain, pairing='nearest')) + selection = util.source_selection_all(len(x0)) + return delay, weight, sos, phaseshift, selection, secondary_source_point(c) + + +def nfchoa_25d_point(x0, r0, xs, fs, max_order=None, c=None, s2z=matchedz_zpk): + r"""Virtual Point source by 2.5-dimensional NFC-HOA. + + .. math:: + + D(\phi_0, s) = + \frac{1}{2\pi r_\text{s}} + \e{\frac{s}{c}(r_0-r_\text{s})} + \sum_{m=-M}^{M} + \Big(\frac{s-\frac{c}{r_\text{s}}\sigma_0}{s-\frac{c}{r_0}\sigma_0}\Big)^\mu + \prod_{l=1}^{\nu} + \frac{(s-\frac{c}{r_\text{s}}\sigma_l)^2-(\frac{c}{r_\text{s}}\omega_l)^2} + {(s-\frac{c}{r_0}\sigma_l)^2+(\frac{c}{r_0}\omega_l)^2} + \e{\i m(\phi_0 - \phi_\text{s})} + + The driving function is represented in the Laplace domain, + from which the recursive filters are designed. + :math:`\sigma_l + \i\omega_l` denotes the complex roots of + the reverse Bessel polynomial. + The number of second-order sections is + :math:`\nu = \big\lfloor\tfrac{|m|}{2}\big\rfloor`, + whereas the number of first-order section :math:`\mu` is either 0 or 1 + for even and odd :math:`|m|`, respectively. + + Parameters + ---------- + x0 : (N, 3) array_like + Sequence of secondary source positions. + r0 : float + Radius of the circular secondary source distribution. + xs : (3,) array_like + Virtual source position. + fs : int + Sampling frequency in Hertz. + max_order : int, optional + Ambisonics order. + c : float, optional + Speed of sound in m/s. + s2z : callable, optional + Function transforming s-domain poles and zeros into z-domain, + e.g. :func:`matchedz_zpk`, :func:`scipy.signal.bilinear_zpk`. + + Returns + ------- + delay : float + Overall delay in seconds. + weight : float + Overall weight. + sos : list of numpy.ndarray + Second-order section filters :func:`scipy.signal.sosfilt`. + phaseshift : (N,) numpy.ndarray + Phase shift in radians. + + Examples + -------- + .. plot:: + :context: close-figs + + delay, weight, sos, phaseshift, selection, secondary_source = \ + sfs.time.drivingfunction.nfchoa_25d_point(array.x, R, xs, fs) + d = sfs.time.drivingfunction.nfchoa_25d_driving_signals( + delay, weight, sos, phaseshift, signal) + plot(d, selection, secondary_source, t=ts) + + """ + if max_order is None: + max_order = util.max_order_circular_harmonics(len(x0)) + if c is None: + c = default.c + + x0 = util.asarray_of_rows(x0) + xs = util.asarray_1d(xs) + phi0, _, _ = util.cart2sph(*x0.T) + phis, _, rs = util.cart2sph(*xs) + phaseshift = phi0 - phis + + delay = (rs - r0) / c + weight = 1 / 2 / np.pi / rs + sos = [] + for m in range(max_order + 1): + _, p, _ = besselap(m, norm='delay') + s_zeros = c / rs * p + s_poles = c / r0 * p + s_gain = 1 + z_zeros, z_poles, z_gain = s2z(s_zeros, s_poles, s_gain, fs) + sos.append(zpk2sos(z_zeros, z_poles, z_gain, pairing='nearest')) + selection = util.source_selection_all(len(x0)) + return delay, weight, sos, phaseshift, selection, secondary_source_point(c) + + +def nfchoa_3d_plane(x0, r0, npw, fs, max_order=None, c=None, s2z=matchedz_zpk): + r"""Virtual plane wave by 3-dimensional NFC-HOA. + + .. math:: + + D(\phi_0, s) = + \frac{\e{\frac{s}{c}r_0}}{r_0} + \sum_{n=0}^{N} + (-1)^n (2n+1) P_{n}(\cos\Theta) + \Big(\frac{s}{s-\frac{c}{r_0}\sigma_0}\Big)^\mu + \prod_{l=1}^{\nu} + \frac{s^2}{(s-\frac{c}{r_0}\sigma_l)^2+(\frac{c}{r_0}\omega_l)^2} + + The driving function is represented in the Laplace domain, + from which the recursive filters are designed. + :math:`\sigma_l + \i\omega_l` denotes the complex roots of + the reverse Bessel polynomial. + The number of second-order sections is + :math:`\nu = \big\lfloor\tfrac{|m|}{2}\big\rfloor`, + whereas the number of first-order section :math:`\mu` is either 0 or 1 + for even and odd :math:`|m|`, respectively. + :math:`P_{n}(\cdot)` denotes the Legendre polynomial of degree :math:`n`, + and :math:`\Theta` the angle between :math:`(\theta, \phi)` + and :math:`(\theta_\text{pw}, \phi_\text{pw})`. + + Parameters + ---------- + x0 : (N, 3) array_like + Sequence of secondary source positions. + r0 : float + Radius of the spherical secondary source distribution. + npw : (3,) array_like + Unit vector (propagation direction) of plane wave. + fs : int + Sampling frequency in Hertz. + max_order : int, optional + Ambisonics order. + c : float, optional + Speed of sound in m/s. + s2z : callable, optional + Function transforming s-domain poles and zeros into z-domain, + e.g. :func:`matchedz_zpk`, :func:`scipy.signal.bilinear_zpk`. + + Returns + ------- + delay : float + Overall delay in seconds. + weight : float + Overall weight. + sos : list of numpy.ndarray + Second-order section filters :func:`scipy.signal.sosfilt`. + phaseshift : (N,) numpy.ndarray + Phase shift in radians. + + """ + if max_order is None: + max_order = util.max_order_spherical_harmonics(len(x0)) + if c is None: + c = default.c + + x0 = util.asarray_of_rows(x0) + npw = util.asarray_1d(npw) + phi0, theta0, _ = util.cart2sph(*x0.T) + phipw, thetapw, _ = util.cart2sph(*npw) + phaseshift = np.arccos(np.dot(x0 / r0, -npw)) + + delay = -r0 / c + weight = 4 * np.pi / r0 + sos = [] + for m in range(max_order + 1): + _, p, _ = besselap(m, norm='delay') + s_zeros = np.zeros(m) + s_poles = c / r0 * p + s_gain = 1 + z_zeros, z_poles, z_gain = s2z(s_zeros, s_poles, s_gain, fs) + sos.append(zpk2sos(z_zeros, z_poles, z_gain, pairing='nearest')) + selection = util.source_selection_all(len(x0)) + return delay, weight, sos, phaseshift, selection, secondary_source_point(c) + + +def nfchoa_3d_point(x0, r0, xs, fs, max_order=None, c=None, s2z=matchedz_zpk): + r"""Virtual point source by 3-dimensional NFC-HOA. + + .. math:: + + D(\phi_0, s) = + \frac{\e{\frac{s}{c}(r_0-r_\text{s})}}{4 \pi r_0 r_\text{s}} + \sum_{n=0}^{N} + (2n+1) P_{n}(\cos\Theta) + \Big(\frac{s-\frac{c}{r_\text{s}}\sigma_0}{s-\frac{c}{r_0}\sigma_0}\Big)^\mu + \prod_{l=1}^{\nu} + \frac{(s-\frac{c}{r_\text{s}}\sigma_l)^2-(\frac{c}{r_\text{s}}\omega_l)^2} + {(s-\frac{c}{r_0}\sigma_l)^2+(\frac{c}{r_0}\omega_l)^2} + + The driving function is represented in the Laplace domain, + from which the recursive filters are designed. + :math:`\sigma_l + \i\omega_l` denotes the complex roots of + the reverse Bessel polynomial. + The number of second-order sections is + :math:`\nu = \big\lfloor\tfrac{|m|}{2}\big\rfloor`, + whereas the number of first-order section :math:`\mu` is either 0 or 1 + for even and odd :math:`|m|`, respectively. + :math:`P_{n}(\cdot)` denotes the Legendre polynomial of degree :math:`n`, + and :math:`\Theta` the angle between :math:`(\theta, \phi)` + and :math:`(\theta_\text{s}, \phi_\text{s})`. + + Parameters + ---------- + x0 : (N, 3) array_like + Sequence of secondary source positions. + r0 : float + Radius of the spherial secondary source distribution. + xs : (3,) array_like + Virtual source position. + fs : int + Sampling frequency in Hertz. + max_order : int, optional + Ambisonics order. + c : float, optional + Speed of sound in m/s. + s2z : callable, optional + Function transforming s-domain poles and zeros into z-domain, + e.g. :func:`matchedz_zpk`, :func:`scipy.signal.bilinear_zpk`. + + Returns + ------- + delay : float + Overall delay in seconds. + weight : float + Overall weight. + sos : list of numpy.ndarray + Second-order section filters :func:`scipy.signal.sosfilt`. + phaseshift : (N,) numpy.ndarray + Phase shift in radians. + + """ + if max_order is None: + max_order = util.max_order_spherical_harmonics(len(x0)) + if c is None: + c = default.c + + x0 = util.asarray_of_rows(x0) + xs = util.asarray_1d(xs) + phi0, theta0, _ = util.cart2sph(*x0.T) + phis, thetas, rs = util.cart2sph(*xs) + phaseshift = np.arccos(np.dot(x0 / r0, xs / rs)) + + delay = (rs - r0) / c + weight = 1 / r0 / rs + sos = [] + for m in range(max_order + 1): + _, p, _ = besselap(m, norm='delay') + s_zeros = c / rs * p + s_poles = c / r0 * p + s_gain = 1 + z_zeros, z_poles, z_gain = s2z(s_zeros, s_poles, s_gain, fs) + sos.append(zpk2sos(z_zeros, z_poles, z_gain, pairing='nearest')) + selection = util.source_selection_all(len(x0)) + return delay, weight, sos, phaseshift, selection, secondary_source_point(c) + + +def nfchoa_25d_driving_signals(delay, weight, sos, phaseshift, signal): + """Get 2.5-dimensional NFC-HOA driving signals. + + Parameters + ---------- + delay : float + Overall delay in seconds. + weight : float + Overall weight. + sos : list of array_like + Second-order section filters :func:`scipy.signal.sosfilt`. + phaseshift : (N,) array_like + Phase shift in radians. + signal : (L,) array_like + float + Excitation signal consisting of (mono) audio data and a sampling + rate (in Hertz). A `DelayedSignal` object can also be used. + + Returns + ------- + `DelayedSignal` + A tuple containing the delayed signals (in a `numpy.ndarray` + with shape ``(L, N)``), followed by the sampling rate (in Hertz) + and a (possibly negative) time offset (in seconds). + + """ + data, fs, t_offset = util.as_delayed_signal(signal) + N = len(phaseshift) + out = np.tile(np.expand_dims(sosfilt(sos[0], data), 1), (1, N)) + for m in range(1, len(sos)): + modal_response = sosfilt(sos[m], data)[:, np.newaxis] + out += modal_response * np.cos(m * phaseshift) + return util.DelayedSignal(2 * weight * out, fs, t_offset + delay) + + +def nfchoa_3d_driving_signals(delay, weight, sos, phaseshift, signal): + """Get 3-dimensional NFC-HOA driving signals. + + Parameters + ---------- + delay : float + Overall delay in seconds. + weight : float + Overall weight. + sos : list of array_like + Second-order section filters :func:`scipy.signal.sosfilt`. + phaseshift : (N,) array_like + Phase shift in radians. + signal : (L,) array_like + float + Excitation signal consisting of (mono) audio data and a sampling + rate (in Hertz). A `DelayedSignal` object can also be used. + + Returns + ------- + `DelayedSignal` + A tuple containing the delayed signals (in a `numpy.ndarray` + with shape ``(L, N)``), followed by the sampling rate (in Hertz) + and a (possibly negative) time offset (in seconds). + + """ + data, fs, t_offset = util.as_delayed_signal(signal) + N = len(phaseshift) + out = np.tile(np.expand_dims(sosfilt(sos[0], data), 1), (1, N)) + for m in range(1, len(sos)): + modal_response = sosfilt(sos[m], data)[:, np.newaxis] + out += (2 * m + 1) * modal_response * legendre(m, np.cos(phaseshift)) + return util.DelayedSignal(weight / 4 / np.pi * out, fs, t_offset + delay) diff --git a/sfs/time/source.py b/sfs/time/source.py index cfb6898..cc0154e 100644 --- a/sfs/time/source.py +++ b/sfs/time/source.py @@ -17,7 +17,8 @@ ts = rs / sfs.default.c # time-of-arrival at origin # Impulsive excitation - signal = unit_impulse(512), 44100 + fs = 44100 + signal = unit_impulse(512), fs grid = sfs.util.xyz_grid([-2, 3], [-1, 2], 0, spacing=0.02) diff --git a/sfs/util.py b/sfs/util.py index 8f02ef7..a97bfff 100644 --- a/sfs/util.py +++ b/sfs/util.py @@ -616,3 +616,45 @@ def source_selection_focused(ns, x0, xs): def source_selection_all(N): """Select all secondary sources.""" return np.ones(N, dtype=bool) + + +def max_order_circular_harmonics(N): + r"""Maximum order of 2D/2.5D HOA. + + It returns the maximum order for which no spatial aliasing appears. + It is given on page 132 of [Ahrens2012]_ as + + .. math:: + \text{max_order} = + \begin{cases} + N/2 - 1 & \text{even}\;N \\ + (N-1)/2 & \text{odd}\;N, + \end{cases} + + which is equivalent to + + .. math:: + \text{max_order} = \big\lfloor \frac{N - 1}{2} \big\rfloor. + + Parameters + ---------- + N : int + Number of secondary sources. + + """ + return (N - 1) // 2 + + +def max_order_spherical_harmonics(N): + r"""Maximum order of 3D HOA. + + .. math:: + \text{max_order} = \lfloor \sqrt{N} \rfloor - 1. + + Parameters + ---------- + N : int + Number of secondary sources. + + """ + return int(np.sqrt(N) - 1)