WFS secondary source selection

[hagenw]

The source selection for WFS is described as "Only secondary sources on the surface A that can directly ’see’ the point source should contribute to the integral in order to obtain a more accurate
result" (compare Herrin 2003 or page 27 in Schultz 2016).

I'm not sure if the current implementation of the source selection is correct as it does not consider shadow effects of other secondary sources. Here is an example of a difficult case.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import sfs

xs = 0, 2, 0
# loudspeaker array
N = 20
radius = 1
# start with a circle and split it
x0_circle, n0_circle, a0_circle = sfs.array.circular(N, radius)
assert N % 2 == 0
upper_x, lower_x = np.split(x0_circle, [N//2 + 1])
upper_n, lower_n = np.split(n0_circle, [N//2 + 1])
# reverse order of upper loudspeakers
upper_x = np.flipud(upper_x)
upper_n = np.flipud(upper_n)
# reverse direction of lower loudspeakers
lower_n = -lower_n
# put all together
x0 = np.row_stack((
    lower_x - [2 * radius, 0, 0],
    upper_x,
    lower_x + [2 * radius, 0, 0],
))
n0 = np.row_stack((
    lower_n,
    upper_n,
    lower_n,
))
a = sfs.mono.drivingfunction.source_selection_point(n0, x0, xs)

fig, ax = plt.subplots()
# selected secondary sources
for x00 in x0[a,:]:
    ss = plt.Circle(x00[0:2], .05, edgecolor='k', facecolor='k')
    ax.add_artist(ss)
# not selected secondary sources
for x00 in x0[~a,:]:
    ss = plt.Circle(x00[0:2], .05, edgecolor='k', facecolor='w')
    ax.add_artist(ss)
# virtual source
vs = plt.Circle(xs[0:2], .1, edgecolor='b', alpha=.5)
ax.add_artist(vs)
# visible elements
p1 = np.append(x0[24:29,0:2], [xs[0:2]], axis=0)
p2 = np.append(x0[0:5,0:2], [xs[0:2]], axis=0)
p3 = np.append(x0[11:18,0:2], [xs[0:2]], axis=0)
for points in p1, p2, p3:
   patch = mpatches.Polygon(points, closed=True, facecolor='r', alpha=0.5)
   ax.add_patch(patch)

ax.grid()
ax.axis('equal')
plt.axis([-3.5, 3.5, -2, 3])
plt.show()

The filled black circles show the loudspeakers that were selected by the secondary source selection for a point source, the empty circles show the non-active loudspeakers. But as I understand the visible elements, only the loudspeakers highlighted by the red area should be selected.

References

Herrin, F. Martinus, T. W. Wu, and A. F. Seybert, “A New Look at the High Frequency Boundary Element and Rayleigh Integral Approximations,” SAE Technical Paper, doi:10.4271/2003-01-1451 (2003).

Schultz, F., Sound Field Synthesis for Line Source Array Applications in Large-Scale Sound Reinforcement, Ph.D. dissertation, TU Berlin, Berlin, Germany, 2016 (submitted).

[fs446]

IMO the cited "Only secondary sources on the surface A that can directly ’see’ the point source should contribute to the integral in order to obtain a more accurate result" is an oversimplification that should be avoided in future. I think it's correct to implement (2.114) and (2.115) from [Schultz 2016] for arbitrary SSDs, i.e. the black dots in above example, since this is what you get from math. However, I did not check this in detail->TBD.

Please cite [Schultz 2016] as (Berlin changed to Rostock + release date)
Schultz, F., Sound Field Synthesis for Line Source Array Applications in Large-Scale Sound Reinforcement, Ph.D. dissertation, University of Rostock, Rostock, Germany, 2016 (submitted, version 2015-12-13).

[trettberg]

(I ran into the same problem, discussed it with Frank and he pointed me here. Some thoughts:)

Let
V be the union of the listening area and its boundary,
comp(V) its complement,
x the position of a virtual point source in comp(V).

Consider the two secondary source selection strategies based on:
I visibility,
II projection on inward normal (current implementation).

I think the following holds:
I and II are equivalent for all x iff V is convex or comp(V) is convex.

If the condition is not met, it seems that either I or II may be preferable. That would depend on geometry and whether the boundary is assumed to be rigid or not.

(Sorry if that is either obvious or obviously wrong.)

[hagenw]

Good to see that someone else is also wondering about the selection.
You are right, it does not really matter for convex surfaces and the classical solution using II states convex arrays as a requirement, compare [Ahrens 2012] or page 684 in Lax and Feshbach 1947 which reads "This restricts our discussion, for the moment, to convex surfaces."

If we have non-convex surfaces, solution I seems to be the correct one, compare Appendix in Lax and Feshbach 1947 which reads "If, however, the surface is shaped so that one part of it blocks the other and casts a shadow, then only those points h; contribute to p(r, theta) which are visible from direction theta. The others must be omitted."

But I guess for a practical setup your are right that one should test what gives a good result.
I compared both cases for the secondary source setup from above in Matlab:

clear all;
conf = SFS_config;
X = [-6 6];
Y = [-6 1.5];
Z = 0;
xs = [0, 2, 0];
src = 'ps';
conf.xref = [0, -1, 0];
f = 700;
% Create convex + concave shaped secondary sources
conf.secondary_sources.number = 20; % number of one circle
conf.secondary_sources.size = 2;
x0_circle = secondary_source_positions(conf);
if isodd(conf.secondary_sources.number)
    error('The number of secondary sources on the circle has to be even.');
end
% Split circle in two parts (the upper part has 2 sources more than the lower
% one)
x0_upper = x0_circle(1:conf.secondary_sources.number/2+1,:);
x0_lower = x0_circle(conf.secondary_sources.number/2+2:end,:);
% Reverse direction of lower secondary sources
x0_lower(:,4:6) = -1 * x0_lower(:,4:6);
% Put all together
x0(:,1:3) = [
    x0_lower(:,1:3) - repmat([conf.secondary_sources.size 0 0], ...
        [size(x0_lower,1) 1])
    x0_upper(end:-1:1,1:3)
    x0_lower(:,1:3) + repmat([conf.secondary_sources.size 0 0], ...
        [size(x0_lower,1) 1])
];
x0(:,4:7) = [
    x0_lower(:,4:7)
    x0_upper(end:-1:1,4:7)
    x0_lower(:,4:7)
];
x0 = secondary_source_selection(x0,xs,src);
conf.secondary_sources.geometry = 'custom';
conf.secondary_sources.x0 = x0;
conf.plot.usenormalisation = false;
[P,~,~,~,x0] = sound_field_mono_wfs(X,Y,Z,xs,src,f,conf);
plot_sound_field(10*P,X,Y,Z,x0,conf);
[P,~,~,~,x0] = sound_field_mono_wfs(X,-4,Z,xs,src,f,conf);
conf.plot.usedb = true;
plot_sound_field(10*P,X,-4,Z,x0,conf);

clear all;
conf = SFS_config;
X = [-6 6];
Y = [-6 1.5];
Z = 0;
xs = [0, 2, 0];
src = 'ps';
conf.xref = [0, -1, 0];
f = 700;
% Create convex + concave shaped secondary sources
conf.secondary_sources.number = 20; % number of one circle
conf.secondary_sources.size = 2;
x0_circle = secondary_source_positions(conf);
% Split circle in two parts (the upper part has 2 sources more than the lower
% one)
x0_upper = x0_circle(1:conf.secondary_sources.number/2+1,:);
x0_lower = x0_circle(conf.secondary_sources.number/2+2:end,:);
% Reverse direction of lower secondary sources
x0_lower(:,4:6) = -1 * x0_lower(:,4:6);
% Put all together
x0(:,1:3) = [
    x0_lower(:,1:3) - repmat([conf.secondary_sources.size 0 0], ...
        [size(x0_lower,1) 1])
    x0_upper(end:-1:1,1:3)
    x0_lower(:,1:3) + repmat([conf.secondary_sources.size 0 0], ...
        [size(x0_lower,1) 1])
];
x0(:,4:7) = [
    x0_lower(:,4:7)
    x0_upper(end:-1:1,4:7)
    x0_lower(:,4:7)
];
x0 = secondary_source_selection(x0,xs,src);
% Remove secondary sources that are not illuminated by source
% see: https://github.com/sfstoolbox/sfs-documentation/issues/8
x0 = [x0(1:5,:); x0(9:15,:); x0(19:23,:)];
conf.secondary_sources.geometry = 'custom';
conf.secondary_sources.x0 = x0;
conf.plot.usenormalisation = false;
[P,~,~,~,x0] = sound_field_mono_wfs(X,Y,Z,xs,src,f,conf);
plot_sound_field(10*P,X,Y,Z,x0,conf);
[P,~,~,~,x0] = sound_field_mono_wfs(X,-4,Z,xs,src,f,conf);
conf.plot.usedb = true;
plot_sound_field(10*P,X,-4,Z,x0,conf);

The results are a little bit interfered by the fact that I'm using only 20 secondary source, but I did the selection of the secondary sources for the second case graphical and was hence not able to use a higher number. Do you have an idea how we could come up with a formula for the second case (which corresponds to I)?

[trettberg]

Assuming that

• All secondary sources are non-directional.
• The boundary is only notional, i.e. has no influence on sound propagation.

Then, for a non-convex volume V neither of the two approaches seems to work.
Consider your example in time domain for a listener located in the middle arc, say at (0,0.5):

• First the desired wave front arrives from the speakers that are visible to the listener.
• Later, unwanted wave fronts arrive from wrong directions, caused by the active sources in the outer arcs.

It seems that visibility of secondary sources from the listener is crucial for non-convex V to avoid back-propagation. (In addition to one of the other two selection criteria.)

Suppose we want to avoid back-propagation for all listener positions in some region R contained in V. This seems to be possible iff the convex hull of R is contained in V. In which case it is possible for all points in the convex hull.

Caveat: All this may be horribly wrong.

[hagenw]

Mh, another good point. For me it was obvious that the listening area can only start for values of y < -1, but you are right we have no mathematical definition of what the listening area is besides V, which obviously is not correct in this case.
For the implementation this is not a big problem as we normally don't check if the given listening position makes any sense at all and leave this to the user.

[trettberg]

For the case above (finite, open, non-convex array) I'm not not even sure about the definition of V.

I think with a bit of work, one could come up with a mathematical formulation of a "desirable" selection strategy for non-convex SSDs. On the other hand, I suspect that it could not be implemented without major changes. (An example: For closed arrays: Does x0 uniquely determine a polyhedron V ? If not, the current SSD definition does not work.)
To me, it seems like a lot of trouble for a somewhat exotic use case.

Maybe we should stick with the current behaviour:

• convex case: Selection works correct.
• non-convex case: Selection works as described. Not guaranteed to be meaningful.

[hagenw]

OK, for the implementation I fully agree and would also propose to simply stay with the current solution.

For the documentation, should we then do the following?

• start with a section on convex arrays and explain that in this case a projection on inward normal is used (which we can then reference from the implementation)
• in addition, state that in principal also non-convex arrays are allowed, but in this case it might be better to switch to another selection criterion, based on visibility

For the last point I'm also not sure, if we have to state that "visibility" is the correct criterion in this case (as for example stated in the Lax and Feshbach paper) or if we only say "it might be better" to switch the criterion?

[trettberg]

For the non-convex case, personally I would be very reluctant to give any recommendation, at least for the moment. Some more nitpicking:

To my understanding, the argument in the Lax and Feshbach paper does not really apply for us. They consider far-field radiation of a finite-extent source:

1. As long as all field points of interest are sufficiently far away (outside the convex hull of the source), back-propagation is not a problem at all. (At least, not in the sense we discussed above.)
2. The surface of the radiator is physical (not only notional), so it's appropriate to use only the part of the surface visible from the field point.

The second point is in contrast to our case:
Consider in your example a listener at (-10,-1) and a virtual point source at (10,1).
Each of the selection criteria we discussed so far (visibility from virtual source, surface normals, visibility from listener) will discard some secondary sources. But for this case it seems desirable that in fact all secondary sources contribute.

I think a useful criterion for the non-convex case depends on

• the SSD, including secondary source characteristics
• the desired listening region
• the admissable region for virtual sources