A modal analysis of multichannel crosstalk cancellation systems and their relationship to amplitude panning

The single-listener multichannel crosstalk cancellation system is analysed using the singu- lar value decomposition. Analytic expressions for the singular system are derived considering a sound-ﬁeld of point sources placed in the far-ﬁeld and in free-ﬁeld conditions. The derived singular system yields theoretical insight into the beneﬁt of multichannel loudspeaker geometries over the conventional two-channel system. From the singular system the source strength solution is derived for a target far-ﬁeld virtual source. The low frequency approximation of the derived source strengths is an amplitude panning function. This panning function is a generalisation of the sine law solution that allows for head rotation and asymmetric loudspeaker geometry. When the geometry is symmetric about the median plane of the listener the low frequency sources strengths are the minimum (cid:2) 2 norm solution of the multichannel sine law. The low frequency approximation is also shown to be representative when the free-ﬁeld assumption is relaxed.


Introduction
Crosstalk cancellation (CTC) systems are designed to deliver a binaural signal at the ears of a listener using loudspeakers rather than headphones [1,2] . The binaural signal must be filtered prior to reaching the loudspeakers in such a way that the desired signals are reproduced at designated control points. Various optimisation problems have been proposed to obtain the necessary loudspeaker filters, e.g., least-mean-squares (LMS) approaches in [3,4] , a pressure matching method (PMM) approach in [5,6] , a recursive least-squares (RLMS) approach in [7] , a minimax approach in [8] , controller methods in [9,10] , and a listener-adaptive approach based on the PMM in [11] . Thus far, study of the physical mechanisms underlying the CTC system have been largely limited to the two-channel two-point (single-listener) case under ideal acoustic assumptions, e.g., [12,13] . The important conclusion drawn from previous two-channel analyses is that the performance of the CTC system is frequency and geometry dependent [13] . For example, in [13] it was shown, assuming monopole sources and a shadowless head model, that the two-channel radiation matrix is ill-conditioned at uniformly distributed frequencies. An ill-conditioned plant matrix poses a practical difficulty as the accuracy of the reproduced solution can be easily compromised [14][15][16] . Furthermore, at frequencies where the system is ill-conditioned the inverse plant matrix can have a large 2 norm [13] . This potentially large amplification can lead to practical issues such as loss of dynamic range, loss of audio quality, and in extreme 1. The analytic expressions for the singular system of the L -channel CTC system and inverse filters are derived. The derivation accounts for asymmetric loudspeaker geometry and head rotation under a far-field source assumption. The singular values in particular give rigorous evidence of the physical aspects of the system that determine the energy requirement needed in the system inversion and the system robustness. 2. A low frequency approximation of the L -channel source strength solution is derived assuming a virtual far-field monopole source. This low frequency approximation is in fact an amplitude panning function that extends the known solution to the stereo sine law to an arbitrary channel count and geometry. When geometric symmetry is enforced about the median plane the low frequency source strengths are the minimum 2 norm solution of the multichannel sine law. This result bridges two well-established sound reproduction techniques, amplitude panning and CTC, under the same inverse problem framework. 3. A comparison of the derived source strengths to source strengths calculated using freely available Neumann KU 100 head-related transfer functions (HRTFs).
The remainder of this work is outlined as follows. Section 2 establishes the necessary background theory of the multichannel CTC system. In Section 3 basic SVD theory is reviewed in the context of the CTC system in order to aid the interpretation of the following modal analysis. Section 4 defines the acoustic model and radiation matrix used in the derivation. Section 5 presents the results of the modal analysis for asymmetric geometries and discusses aspects of the mode efficiency and inverse filters, specifically for multichannel geometries and with head rotation. Section 6 presents a modal analysis of the same system from Section 5 after forcing geometric symmetry about the median plane of the listener. Having established both singular systems, Section 7 derives the source strengths when the binaural signals are the pressures at the ears due to a far-field monopole source. Therein, the low frequency approximation of the source strengths is derived and found to be a multichannel amplitude panning function. In Section 8 an evaluation of the derived source strengths is presented having relaxed the free-field condition. Finally, Section 9 provides a summary of this work and some concluding remarks. Much of the contribution of this work is based on the analytical results, the fine details of which are included in the appendices. The derivations of the singular systems for asymmetric and symmetric geometries can be found in Appendix A.1 and A.2 , respectively.

Multichannel crosstalk cancellation system
The L -channel ( L ≥ 2) two-point crosstalk cancellation system can be described in the frequency domain by where p (ω) ∈ C 2 are the reproduced binaural signals at the ears (or "control points"), G (ω) ∈ C 2 ×L is the radiation matrix of electro-acoustic transfer functions between the L sources and the two control points, q (ω) ∈ C L are the unknown source strengths necessary to drive the loudspeakers in order to recreate the pressures p ( ω), and ω ∈ R is radian frequency. Fig. 1 shows a top-down view of an example L -channel two-point CTC geometry. The desired outcome is that the reproduced pressures at the control points equal an exact delayed copy of the desired binaural signal d (ω) ∈ C 2 such that where t is a modelling delay in seconds necessary to ensure a causal solution, e.g., [3] . In order to achieve the equality stated in Eq. (2) the source strengths are calculated as where H (ω) ∈ C L ×2 is the matrix of filters such that G (ω) H (ω) = I e −j ω t (assuming no regularisation is applied). A particular solution, but not the only solution, is given by where (·) + denotes the pseudoinverse, e.g., [19] . When L > 2 the system is underdetermined and the pseudoinverse will yield a unique and exact solution, in theory, that minimises the Frobenius norm of the inverse filters, i.e., min H ( ω ) F . The filters given by Eq. (4) are thus referred to as being "minimum norm" [19] .
Note that the delay e −j ω t and the explicit dependence on ω will be omitted in the following analysis for the sake of brevity although their presence is implied.

Singular value decomposition
In this section the SVD is introduced in the context of the CTC system along with the concepts of radiation matrix modes and the associated mode efficiencies.

Singular value decomposition of G
The singular value decomposition of the radiation matrix G is a matrix factorisation of the form [19] where U ∈ C 2 ×2 is a unitary matrix of left singular vectors, V ∈ C L ×L is a unitary matrix of right singular vectors, and ∈ R 2 ×L ≥0 is a rectangular diagonal matrix containing the singular values along the main diagonal such that where all other entries off the main diagonal are zero, ( · ) T denotes the matrix transpose and ( · ) H denotes the Hermitian transpose. The chosen ordering of the singular values is by descending magnitude, i.e., σ 1 ≥ σ 2 , unless otherwise stated. The pseudoinverse of G can be found by rearranging Eq. (5) to where + ∈ R L ×2 ≥0 is a rectangular diagonal matrix whose elements along the main diagonal are the inverse singular values of G such that where all elements off the main diagonal are zero. The 2 norm of both the radiation matrix and its pseudoinverse are given by the largest singular value of each matrix, respectively, i.e., [19] Additionally, the condition number of both matrices is [19] Thus, if the singular values are known in terms of the physical parameters of the system, one has obtained some knowledge of how those physical parameters influence the system condition number and the 2 norm of the inverse filters.

Radiation matrix modes
Given G , the left singular vectors U = u 1 u 2 are considered the "field pressure" modes while the right singular vec- . v L are considered the "source strength" modes, e.g., [22,32] . The field pressure and source strength modes form an orthonormal basis for C 2 and C L , respectively. In this context the modes form a basis for the spaces of possible reproduced field pressures and source strengths. The individual modes and the radiation matrix are related by where i = 1 , . . . , rank ( G ) . The right singular vectors with indices rank( G ) < i ≤ L span the null space of G [32] . Conceptually these source strengths modes are not "seen" at the points in space where the field pressure modes are located [32] . As such, these unseen modes are not utilised. For this reason, and for compactness, the following only considers the source strength modes v 1 and v 2 explicitly.

Sound-field model
In this section a sound-field model is introduced that allows for some physical interpretation of the modal decomposition. The coordinate system and the radiation matrix are defined.

Interaural-polar coordinate system
A variation of the interaural-polar coordinate system, e.g., [33] , as shown in Fig. 2 , is used throughout. The origin of the coordinate system is placed at the centre of the listener's head. The conversion between Cartesian and interaural-polar coordinates is given by x = r cos γ cos ψ, (13) z = r cos γ sin ψ. (15) Note that γ ∈ [0, 2 π ] and ψ ∈ − π 2 , π 2 (this convention swaps the domains of γ and ψ compared to [33] ). The forward facing position is when the median plane of the listener is the xz plane.

Source model and radiation matrix
Throughout it is assumed that each loudspeaker is an acoustic monopole in free-field and that the pressure is measured at the control points by point source microphones. The transfer function from the l th monopole source positioned at s l ∈ R 3 to the m th control point positioned at where k = 2 π f /c is the wavenumber ( f is frequency in Hz, c is the speed of sound (m · s −1 )), R ml = s l − x m 2 is the distance in metres from the l th source to the m th control point and l = 1 , 2 , . . . , L and m = 1 , 2 . This model assumes that the physical acoustics of the listener are not accounted for, i.e., a shadowless head model. The lack of a real head in the model is a suitable assumption for the low frequency analysis, e.g., [18,31] . The sources are assumed to be positioned sufficiently far away from the listener such that s l 2 x m 2 . A far-field approximation is therefore introduced, e.g., [34] , such that the transfer functions are approximately where ˆ n l ∈ R 3 is the unit vector pointing to the l th point source and R l = s l 2 is the distance from the l th point source to the origin. Furthermore, it is assumed that the control points are diametrically opposed such that x 1 = −x 2 . With these assumptions the elements of the radiation matrix are defined compactly as and where δ l = e −j kR l /R l is the l th propagation delay and attenuation and φ l = k T l x 1 is the dot product of the l th wavevector and the left ear position vector, where the l th wavevector is k l = k ˆ n l . This notation makes it clear that the elements of the proposed radiation matrix are delayed and attenuated plane waves. The radiation matrix is now defined as where is the diagonal matrix of propagation delays and attenuations.
With the established physical model the radiation matrix SVD is presented in the following section.

Modal analysis: asymmetric geometry
This section presents the SVD of G (as defined in Section 4.2 ), i.e., for an asymmetric geometry. Analytic expressions are presented for the field pressure and source strength modes as well as the associated singular values.

Singular system: asymmetric geometry
As derived in Appendix A.1 , the SVD of G for asymmetric geometries is and and ( · ) * denotes complex conjugation. Eqs. (22) to (23) indicate that the field pressure modes no longer have the strict in-phase and out-of-phase form as previously seen in the two-and three-channel analyses [13] and [17] , respectively. The asymmetric field pressure modes have been altered by a phase shift of ࢬ α that is present in the second component of each mode. This additional phase shift is a consequence of the asymmetric geometry and does not appear when the geometry is symmetric (see Section 6.2 ). Overall, there is a clear sum and difference trend in the singular system. For example, when applying the inverse filters H to an arbitrary binaural input, i.e., Hd , the matrix U forms a vector of sum and differences of the two components of the binaural signal with one of the binaural signal components having undergone a phase shift by ࢬ α, e.g., A similar sum and difference structure appears in V whose elements are the sum and difference of the source-to-left and -right ear transfer functions (with a phase shift of −∠ α applied to the latter). Additionally, the singular values are distinguished from each other by a sum or difference with | α|. The natural ordering of the singular values is apparent from Eqs. (24) to (25) . Both A and | α| are always real and positive (distances ≤ 0 are not allowed) and thus σ 1 ≥ σ 2 must always be true. From Eq. (11) and Eqs. (24) to (25) the condition number is Inspection of Eq. (31) indicates a condition number of one occurs when | α| = 0 , i.e., the singular values are equal.

Mode efficiency
The mode efficiencies are synonymous with the singular values, e.g., [22] , and will be defined here as the squared singular values which are the eigenvalues of GG H , the so-called Gram matrix, e.g., [19] . The specific form of the singular values given by Eqs. (24) to (25) gives insight into the physical behaviour of the multichannel CTC system. For example, consider that λ −1 2 = H 2 2 and that Additionally, consider that | g H 1 g 1 | and | g H 1 g 2 | represent the beamforming gain, e.g., [35] , measured at the first and second control points, respectively, as a result of filtering with the focusing filters g H 1 , e.g., [36] . In this case the beamforming gain where g m is the set of transfer functions from the loudspeakers to one of the two control points x m , where m = 1 , 2 , and g ( x ) is the set of transfer functions from the loudspeakers to an arbitrary point in space x ∈ R 3 . That is to say, g 1 ≡ g ( x 1 ) and The beamforming gain as it has been defined is not to be confused with a far-field directivity pattern. This form of the eigenvalues is rooted in the fact that the pseudoinverse itself represents a two-part operation, one of which is a focusing stage wherein two sets of focusing filters are formed to create the desired pressures at the "focus points", and the other being the inverse Gram which introduces the actual crosstalk cancellation wherein the null is created at the "null" points (also with normalisation at the focus points). For example, if focusing at x 1 the null point would be x 2 (the opposite ear), and vice versa. This two stage process is made clear from the pseudoinverse written in the form where the focusing filters (each column of G H ) ensure the desired pressure is reproduced at each focus point and, on the other hand, the inverse Gram matrix ensures that the null is created at the opposite ear simultaneously along with normalisation at the focus points. With this in mind, Eq. (34) makes a profound statement about the behaviour of the CTC system, in that, the system is ill-conditioned when the magnitude of pressure at the null point (as a result of the focusing operation) is near or equal to the beamforming gain at the focus point. Put another way, the inversion in the CTC system takes more energy if larger pressure exists at the null point as a result of the focusing filters. In line with this interpretation, the system is well-conditioned when | g H 2 g 1 | ≈ 0 (see Eq. (31) ). Fig. 3 a and b illustrate this concept further by showing the beamforming gain, measured along an arc of points one metre from where the centre speaker of the array would be placed (which was one metre from the centre of the forward facing listener's head) overlaid with the two-channel array geometry for two different frequencies, f ≈ 967 Hz and f ≈ 1963 Hz, respectively. In Fig. 3 a, the beamforming gain at the null point is nearly zero and thus the condition number is near one. In Fig. 3 b, on the other hand, a grating lobe appears in the direction of the null point and thus the condition number is substantially higher. At low frequencies the array is no longer directive and in general the system always has a condition number greater than one. The lack of directivity at low frequencies is shown in Fig. 3 c for the same two-channel geometry when f ≈ 117 Hz. The physical factors that determine the mode efficiency magnitudes can be explored in more detail from the expansions of Eqs. (24) to (25) , which are where φ lk = 2(φ l − φ k ) is the scaled phase difference at the first control point due the l th and k th source. The magnitude of the mode efficiencies varies with the frequency-dependent phase difference φ l k and generally increases with increasing number of sources (see the denominators of Eq. (59) for an example of the latter observation). Inspection suggests that the magnitude of λ 1 is maximised when the phase differences are as close to zero as possible. On the other hand, λ 2 is maximised when the phase differences are as close to ± π /2 as possible. When one mode efficiency is maximised the other is minimised. It is expected that at low frequencies, when the phases are similar, λ 1 will reach a maximum while λ 2 will reach a minimum. Inspection of Eq. (37) also suggests that the two-channel system is prone to a singularity when the sources are equidistant and the head is rotated 90 • such that φ 12 = 0 , resulting in λ 2 = 0 → H 2 = ∞ . This instability can be mitigated by the addition of one or more channels, e.g., addition of a centre source or another pair. The effect of head rotation on the mode efficiencies can be understood in light of Eq. (34) . The placement of the null control point relative to the beamforming gain will dictate the system condition number regardless of the head rotation angle. Fig. 4 a to c show, for L = { 3 , 6 , 12 } symmetric geometries with uniform angular spacing and angular span γ = 60 • (the span is formed by the outermost sources relative to the median plane), respectively, the system condition number κ( G ) as a function of frequency in Hz and head rotation angle in degrees (note that there is no need to consider head rotations . 3. Two-channel beamforming gain | g H 1 g (x ) | measured along an arc one metre from the centre speaker location for frequencies: larger than 90 • due to the symmetry of this acoustic model). The condition number has been calculated by substituting Eqs. (36) to (37) into Eq. (11) . Also note that condition numbers larger than 40 dB appear as 40 dB. These surfaces indicate that the system robustness increases, i.e., κ( G ) decreases as L increases above approximately 1 kHz. This finding can be understood physically from Fig. 5 which shows that in the move from three ( Fig. 5 a) to twelve channels ( Fig. 5 b) the array becomes more directive at higher frequencies and thus | g H 1 g 2 | becomes smaller, thereby improving the system condition number. As a result of the increased directivity the system requires less effort to cancel the energy at the null point. Furthermore, due to the fixed span, the low frequency performance of the system, i.e., below 1 kHz, is largely unaffected by the increase in loudspeaker density, as shown in Fig. 4 .

Inverse filters
From Eq. (7) each column of the inverse filters can be expressed as a linear combination of the source strength modes such that Substitution of the singular system derived in Section 5.1 into Eq. (38) , along with some minor rearrangement of terms, yields Note that λ 1 λ 2 = det (GG H ) and λ 1 + λ 2 = tr (GG H ) , where det( · ) is the determinant of a matrix and tr( · ) is the trace of a matrix. From Eqs. (10) and (25) it can be said that Thus the 2 norm of the inverse filters is minimised per frequency when the beamforming gain at the null point is 0 (see Section 5.2 ). This state also corresponds to the special condition of equal mode efficiency, i.e., λ 1 = λ 2 = A .

Modal analysis: symmetric geometry
This section builds upon the asymmetric geometry singular system analysis in Section 5 by considering the special case of loudspeaker and listener geometry symmetric about the median plane.

Radiation matrix: symmetric geometry
For the symmetric geometry it is assumed that the transfer functions are mirrored about the median plane (see Section 4.1 ). Odd numbers of loudspeakers have an additional centre source (placed directly in front of the listener) and no restriction is put on the height of the loudspeakers relative to the listener. The loudspeaker arrangement itself is symmetric about the median plane and the listener is in the centre forward facing position. Recall that the control points are also assumed to be diametrically opposed. Thus the ear position vectors simplify to where a = x 1 2 = x 2 2 is the head radius in metres. This means that φ l = ka sin γ l under these symmetric conditions, where γ l is the l th loudspeaker angle from the median plane. With the stated assumptions the radiation matrices for even and odd L are

Singular system: symmetric geometry
The detailed derivation of the SVD of G under geometric symmetry is given in Appendix A.2 . For both even and odd L the field pressure modes are For even L the singular values and source strength modes are and for odd L the singular values and source strength modes are σ + ,odd = 2 Note that the subscripts + and -have been introduced to distinguish the symmetric singular system from the asymmetric.
This change in notation also implies that there is no logical ordering by magnitude in this derivation, i.e., σ + is not necessarily larger than σ − , and vice versa (see Appendix A.2 as to why this is the case). Previous studies, e.g., [13,22] , reported the same real-valued field pressure modes, i.e., Eqs. (45) to (46) , for the symmetric two-channel two-point system also using a monopole in free-field model. This simple form led to naming the modes associated with σ + and σ − as "in-phase" and "out-of-phase", respectively [13,22] . The results here show that this simple form is maintained regardless of the number of sources as long as geometric symmetry is enforced. Furthermore, the source strength modes simplify, in comparison to Eqs. (26) and (27) , to only trigonometric functions while exhibiting a distinct in-phase and out-of-phase structure. This structure was previously observed by the authors in [17] for the symmetric three-channel geometry derived using a plane wave model. The results here show the same structure is repeated (with the inclusion of a propagation attenuation and delay matrix) and extended, i.e., the number of elements of the modes is increased by two identical terms (differing by a sign in v −,e v en/odd ) for every additional symmetric two-channel pair. The primary distinction between L -channel even and odd source strength modes is also maintained, e.g., [17] , in that only v + ,odd utilises the centre source.
As reported previously in [17,37] , the theoretical advantage of the three-channel versus two-channel system is the increase in in-phase mode efficiency attributed to the additional centre source. It was shown that in the move from two to three channels half of the symmetric two-channel singular value minima as functions of frequency increase substantially, thereby improving the global in-phase mode efficiency as a function of frequency. This improvement in mode efficiency resulted in a lower inverse filter amplification factor H 2 and improved system condition number κ( G ) for frequencies at which the system was previously ill-conditioned due to the in-phase mode inefficiency [17] . However, the behaviour for L > 3 was then unknown. From the results derived here it is evident that this benefit persists in general for symmetric L > 3, where the increase in in-phase mode efficiency is a direct consequence of the added contribution of the centre source. Under these theoretical assumptions odd multichannel geometries are preferred for this reason and it is the authors' experience that a centre source contributes notable perceived stability and clarity in practical CTC systems.
It was also found previously in [17] that the centre source plays no role in the symmetric three-channel out-of-phase mode efficiency. Examination of Eq. (52) shows that this remains the case for symmetric odd L > 3.

Source strengths
Up to this point in the analysis the inverse of the plant matrix has been examined for an arbitrary target binaural signal. The following restricts the analysis by assuming a target virtual source. The primary focus of this section is to arrive at a low frequency approximation of the derived source strengths. It is shown that this low frequency approximation is a frequencyindependent real-valued "amplitude panning" function. This panning function inherently compensates for head rotation and asymmetric loudspeaker layouts. In the symmetric case it is shown that the derived low frequency solution is the minimum

Source strength definition
From Eqs. (3) and (7) q can be expressed as a linear combination of the source strengths modes such that Eq. (55) indicates that the frequency domain amplification factor on each source strength mode is determined by an associated complex mode weight u H m d /σ m , where m = 1 , 2 . The magnitudes of the inner products largely affect the resulting source strength energy requirement [17] . This energy requirement is given by the squared 2 norm of q , which is Eq. (56) indicates that when a particular mode approaches a minimum in efficiency the potential for division by a small singular value can result in a large magnitude. However, if the magnitude of the Hermitian inner product is comparatively small or smaller the contribution to the overall source strength energy requirement can be low or tolerable, despite low mode efficiency. An intuitive takeaway from this result is that if a mode radiates inefficiently the potentially high energy requirement may be lessened or mitigated entirely as the binaural signal may lack components in the "direction" of the inefficient field pressure modes.

A target virtual source
A specific desired binaural signal is now introduced. Let the desired binaural signal equal the pressure at the control points due to a monopole source in the far-field and free-field conditions such that d = δ s e j φs e −j φs , where δ s = e −j kR s /R s is the source propagation delay and attenuation factor, φ s = k T s x 1 , and k s = k ˆ n s (where ˆ n s is the unit vector in the direction of the virtual source).

Source strengths: asymmetric geometry
Substituting Eq. (57) and the singular system derived in Section 5.1 into Eq. (55) yields the source strengths While Eq. (58) could be explored further, the focus here is to reach a low frequency approximation. Before deriving the low frequency approximation the geometry is simplified to equidistant sources, including the virtual source. This causes D * → δ * and δ s → δ, where δ = e −j kR /R ( R being the common distance in metres), and thus Eq. (58) becomes q = e j φs +e −j ( φs −∠ α) where λ 1 and λ 2 have been expanded and simplified under the equidistant assumption (see A.1 for the expansion before this assumption). Note that the source strengths are no longer dependent on distance.

Source strengths: symmetric geometry and equidistant sources
The source strengths for the symmetric system are found in the same manner, however, this derivation uses the results derived in Section 6.2 and substitutes those into Eq. (55) . This substitution yields q = cos ( φ s ) and Recall that here φ l = ka sin γ l , where l = 1 , 2 , . . . , L/ 2 or (L − 1) / 2 , and φ s = ka sin γ s , where γ s is the target virtual source interaural-polar azimuth angle. In this special case q is real-valued. Each individual source strength is the resultant of components corresponding to the in-phase and out-of-phase components of the l th two-channel pair (of the larger array) scaled by the associated complex mode weight. This multichannel result is a direct extension of the three-channel source strengths derived by the authors in [17] . Furthermore, Eq. (60) indicates that for source phases φ s = ±nπ / 2 , n = 1 , 3 , 5 , · · · the inphase component will go to zero, while for source phases φ s = ±nπ / 2 , n = 0 , 2 , 4 , · · · the out-of-phase component will go to zero. This result was previously observed in [17] for the three-channel system where it was expressed in terms of the virtual source angle.

Low frequency approximation: a generalised amplitude panning function
The low frequency approximation of the sources strengths given by Eq. (59) is where s is the angle between ˆ n s and x 1 and l is the angle between ˆ n l and x 1 , in radians. The cosine of the angles are calculated from cos s = ˆ n T s x 1 /a and cos l = ˆ n T l x 1 /a . The vectors ˆ n s , ˆ n l , and x 1 are not restricted to the horizontal plane here and thus full 3D distributions of sources and all possible head rotations are considered in Eq. (63) . Note that the control points are still assumed to be diametrically opposed.
Eq. (63) is a frequency-independent real-valued amplitude panning function where q l is the l th loudspeaker gain. This panning function is an extension of the solution to the multichannel sine law, e.g., [38,39] . The extended capabilities are compensation for head rotation and asymmetric loudspeaker placement. The result given by Eq. (63) extends the authors' previous work in [18] from the two-channel to multichannel case. In order to fully utilise such a panning function the listener's head position must be tracked, as is done for example with compensated amplitude panning (CAP) [40] . If the head position is not tracked the reproduced low-frequency interaural-time-difference (ITD) will be incorrect and the perceived source image will deviate from its intended target position [40] .
The assertion that Eq. (63) is an extension of the solution to the multichannel sine law has so far been unsupported and it is not obvious from inspection that they are related. The following section will show that the derived low frequency source strengths are a solution to the multichannel sine law under geometric symmetry.

Low frequency approximation: minimum 2 norm solution to the multichannel sine law
Starting from Eq. (60) , the low frequency approximation of the source strengths under geometric symmetry is and 1/ L is q + in the low frequency limit. The relationship to the stereo sine law solution, with the additional constraint that the loudspeakers gains sum to unity, can be verified considering Eqs. (64) to (65) . For L = 2 the components of Eq. (64) are which is the solution to the stereo sine law, e.g., [39] , restated here as where q 1 ∈ R and q 2 ∈ R are the sine law loudspeaker gains with the constraint with the constraint given by Eq. (68) . As with the two-channel case, the assertion that Eq. (64) is a solution to Eq. (69) can be verified by substitution. Furthermore, the elements of Eq. (64) sum to unity in the multichannel case. The assertion that Eq. (64) is the minimum 2 norm solution of the multichannel sine law given by Eq. (69) with the constraint given by Eq. (68) stems from the fact that the inverse-filters from which Eq. (64) is derived are minimum norm. To support this assertion Eq. (64) is compared to numerically calculated gains. The numerically calculated gains are found using CVX, e.g., [41,42] , a software package for specifying and solving convex optimisation programs, to solve the following minimisation problem minimise || q || 2 (70) subject to sin γ S = L l=1 q l sin γ l L l=1 q l L l=1 q l = 1 , which is posed to give the minimum 2 norm solution to (69) .  ( Fig. 6 a and b, respectively). The outermost loudspeakers for each geometry are positioned at 30 • and 330 • in the horizontal plane. The results show that the derived low frequency source strengths are in fact a one-to-one match to the numerical solution of Eq. (70) . It can be concluded that, under the stated assumptions, the derived L -channel CTC source strengths in the low frequency limit are asymptotically equivalent to the minimum 2 norm solution to the multichannel sine law when the gains sum to unity. It is rather remarkable that, despite the known low frequency out-of-phase mode inefficiency which leads to an inversely proportional inverse filter amplification factor, the low frequency source strengths derived here are bounded and rather well behaved when assuming the outermost loudspeaker angular span is not very small. This result means the virtual source can be reproduced with relatively little amplification and can be seen as an example of when the magnitude of the out-of-phase Hermitian inner product in the numerator of Eq. (55) is comparable in magnitude to the associated singular value in the denominator.
On the other hand, it can be seen in Fig. 6 b when L > 2 that the solution is not sparse given that a majority of the gains are non-zero for a given virtual source angle. This property of the low frequency source strengths makes the resulting reproduction highly dependent on the listener position in practical settings as the wavefronts from each loudspeaker are assumed to be exactly time aligned. For example, if the listener deviates from the centre position, or if there are any misaligned loudspeakers, the wavefronts from the nearest loudspeaker will typically reach the listener first and the listener may experience the so-called precedence effect, e.g., [39] , wherein the perceived location of the panned sound source collapses to the nearest loudspeaker, or grouping of nearby loudspeakers.

Comparison to measured source strengths
In this section the results of simulations are presented that compare the derived symmetric source strength solution to source strengths calculated using freely available Neumann KU 100 HRTFs (as opposed to an analytic model of the transfer functions) so as to examine the accuracy of the derived results when the free-field assumption is relaxed.

Simulation setup
The geometry used in the simulation was a five-channel symmetric geometry with uniform angular spacing and outermost angular span of γ = 60 °, constrained to the horizontal plane. This model geometry is shown in Fig. 7 . The measured source strengths ˜ q were calculated according to Eq. (3) with a measured radiation matrix ˜ G comprised of HRTFs corresponding to sources at the angles of the specified five-channel geometry. The HRTFs were originally measured with a Neumann KU 100 mannequin head in an anechoic chamber at the Cologne University of Applied Sciences [43] . These HRTFs were chosen specifically because of the low frequency extension that was applied which allows for the low frequency limit behaviour to be studied below the lowest reproduction frequency of the measurement loudspeaker [43] . The reported distance from the loudspeakers to the centre of the head was 1.5 metres [43] . The desired measured binaural signals ˜ d were also created using the KU 100 HRTFs corresponding to sources at incident angles γ s ∈ [0 • , 359 • ] (with one degree angular resolution), also at a distance of 1.5 metres. For virtual source angles 30 • < γ s < 330 • the source strengths are designed to render a virtual source beyond the physical extent of the loudspeakers. On the other hand, the derived source strengths were calculated using Eq. (60) for the same set of angles γ s . The head radius parameter for the derived filters was chosen as a = 0 . 09 metres, which is approximately the KU 100 head radius. The speed of sound was assumed to be c = 343 m · s −1 .  Fig. 8 c. The results are consistent with the analysis of the three-channel CTC system conducted by the authors in [17] , in that, the magnitude response of the measured solution is comparable to the magnitude response of the analytic solution up to approximately 1 kHz. Furthermore, the measured solution is mostly frequency-independent in this region as predicted by the low frequency approximation given by Eq. (64) . As the measured solution is complex and the analytic solution is real, it is useful to compare the real and imaginary parts separately to learn more about the exact behaviour of the measured solution. Fig. 8 c shows the absolute value of the real part of the measured solution which can be visually compared to the magnitude of the analytic solution in Fig. 8 a. In contrast to the measured magnitude surface in Fig. 8 b the magnitude of the real part of the measured solution displays much more similarity to the analytic solution across the entire frequency range, e.g., the position and number of nodal lines in the high frequencies is comparable. However, there are clearly more complicated phenomena occurring in the measured source strengths at higher frequencies that are not represented by the analytic solution, e.g., variations in curvature of the nodal lines. All of these phenomena observed in the measured solution are attributed to the complicated scattering at high frequencies, head shadowing, and pinna effects that are not taken into account in the shadowless head model. On the other hand, while the magnitude of the real part of the measured solution is comparable to the analytic solution, it is expected that the primary differences between the solutions are due to the imaginary component of the measured solution. However, it is expected that this difference exists for higher frequencies where the imaginary component of the measured solution is expected to have larger magnitude, i.e., above 1 kHz. Fig. 9 shows that the magnitude of the difference between the analytic and measured source strengths ( Fig. 9 a) is indeed similar to the magnitude of the imaginary component of the measured source strengths ( Fig. 9 b), reinforcing this hypothesis. However, the two results are clearly not identical which can be attributed to the differences between the real parts of each solution previously noted.

Simulation results
In summary, the comparison between the analytic and measured source strengths indicates that the analytic solution represents the real part of the measured solution reasonably well, best so at frequencies below 1 kHz. Above 1 kHz high Fig. 9. a) Magnitude of the difference between the analytic and measured source strengths; b) Magnitude of the imaginary part of the measured source strengths. frequency phenomena appear in the real part of the measured solution that are not visible in the derived solution. Additionally, a considerable amount of the difference between solutions above 1 kHz was attributed to the imaginary component of the measured solution. This result is to be expected as the measured solution is not exactly symmetric in practice and includes additional acoustic effects that cause an imaginary component to be present.

Conclusion
The singular value decomposition of the L -channel single-listener CTC system was derived for asymmetric and symmetric geometries assuming far-field monopole sources in free-field and a shadowless head. From analysing the singular values it was concluded that the system condition number and inverse filter energy requirement are determined by the magnitude of the beamforming gain at the null point after focusing filters have been applied. It was shown that the magnitude of the pressure at the null point after focusing is largely determined by the phase differences of the transfer functions at the ears due to each source and between sources. This realisation provided a rigorous explanation as to why increased channel count can increase the mode efficiencies and improve system robustness at mid-low to high frequencies. The physical reason is that the focusing aspect of the CTC system becomes increasingly directional with increasing channel density and therefore the system needs to cancel less pressure at the null point. As such, the effect of head rotation can be generally understood by just considering where the control points lie in the field relative to the beamforming gain. Thus, it can be concluded that multichannel arrays also benefit CTC systems that would involve dynamic head rotation or changing listener position for the same fundamental physical reasons the symmetric system benefits.
Furthermore, it was shown that when the system geometry is symmetric the field pressure modes remain unchanged for any number of sources. The simple form of the derived symmetric singular system extends the previously established concept of in-phase and out-of-phase modes from two-and three-channel versions to the multichannel symmetric geometry. It was shown for any number of odd symmetric channels that the centre source increases the in-phase mode efficiency whilst having no effect on the out-of-phase mode efficiency.
Derivation of the singular systems subsequently allowed specific analysis of a source strength solution for a target virtual point source. The low frequency source strength approximation was derived and shown to be a frequency-independent real-valued amplitude panning function. This panning function naturally extends the solution to the stereo sine law to asymmetric multichannel geometries. It was also shown that when the system geometry is symmetric about the median plane the low frequency source strength approximation is the minimum 2 norm solution of the multichannel sine law with the constraint that the gains sum to unity. This result rigorously connected the multichannel CTC system and amplitude panning under the same inverse problem framework. Furthermore, this simplification to amplitude panning showed that the loudspeaker energy requirement can be tolerable despite low out-of-phase mode efficiency at low frequencies. This result is due to the fact that the target signal has little out-of-phase mode contribution at low frequencies.
Finally, a comparison study was conducted between the analytic symmetric source strengths and source strengths calculated from freely available Neumann KU 100 HRTFs (which were chosen because of their extended low frequency response, see [43] ). The simulation was run for a model five-channel geometry in the horizontal plane. It was shown that the magnitude of the analytic and measured source strengths were comparable up to 1 kHz, above which the derived results began to deviate from the measured. Further observation concluded that the derived solution was comparable to the real part of the measured solution across the audible bandwidth. Specific deviations between solutions were noted in the high frequencies that were attributed to the complex acoustic effects not present in the shadowless head model underlying the derived source strengths. It was also observed that the measured solution has an increasing imaginary component above 1 kHz which is not present in the derived symmetric source strengths. The results from the comparison simulation confirm that the chosen acoustic model does not account well for high frequency head scattering phenomena, however, the perceptual limitations of the CTC system based on the derived inverse filters have not been investigated.

Declaration of Competing Interest
None. The derived singular values and left singular vectors are substituted directly into Eq. (A.3) to obtain the right singular vectors Eqs. (26) to (27) .

A2. Symmetric geometry
The SVD of G under geometric symmetry, i.e., Eqs. (43) to (44) , is derived following the same steps used to derive the asymmetric geometry singular system given in Appendix A.1 . While the form of the two derivations is nearly identical, there is a key distinction. Under geometric symmetry the Gram matrix takes the distinct real form ( 1 ± cos ( 2 φ l ) ) , L is even The singular values under symmetric conditions are then σ ±,sym = λ ±,sym (note that the double angle formula has been applied to obtain the forms Eqs. (47) to (48) and Eqs. (51) to (52) ). The same steps taken to obtain the left singular vectors detailed in Appendix A.1 are repeated. It is found that the left singular vectors must be chosen such that u 11 = u 21