© The Author(s), Published by EDP Sciences, 2023

1 Introduction

The directivity properties of sound sources, in particular musical instruments and voices, have been shown to be a key factor in the perception of acoustical sources in a reverberant environment [2, 3], for the reproduction of sound by electroacoustic devices [4, 5], and for the realistic synthesis of sources in auralization and virtual acoustic frameworks [6–9]. In all these areas of ongoing research, it would be beneficial to formulate a compact way to describe, compare, and classify sound radiation patterns.

Audio descriptors are widely and successfully used for the classification of musical and sonic material [10–12]. However, they are primarily evaluated on monophonic signals and therefore do not capture the complex spatial patterns of acoustical sources such as musical instruments. In order to analyze these patterns computationally, we require a similarity measure for pairwise comparisons, from which we can derive a multidimensional spatial description. The primary aim of our study is to investigate the cross-correlation of these patterns as a similarity measure and clustering tool.

1.1 Motivation

This work was initiated during an artistic-research residency where composer Aaron Einbond sought to extend corpus-based concatenative synthesis (CBCS) with spatial features. CBCS [13] is a sound synthesis method where short sound segments, or units, are automatically selected from a large database of audio samples and then assembled (i.e. concatenated) for playback. By re-arranging units from a corpus of live- or pre-recorded sounds, the technique is capable of synthesizing rich and musically expressive materials [14].

In CBCS, a unit selection algorithm is responsible for finding the sequence of units that best matches a given target specifying the sound or phrase to be synthesized. The selection is performed according to descriptors of the units, that is, feature characteristics extracted from or attributed to the source sounds. A multi-dimensional descriptor space is populated with the sound units and a distance function guides the unit selection algorithm to find the optimal sequence of units to fit the target.

Most often, the descriptors used are acoustical parameters computationally analyzed from the audio signals [15] and expected to be representative of the acoustical structure of the sound units. This typically includes temporal, spectral, harmonic, spectro-temporal, and energetic properties of the sound events, such as fundamental frequency, loudness, spectral centroid, etc. However, these descriptors do not commonly include spatial information, which would be useful especially when considering audio signals from acoustic instruments. One initial motivation for our work is to incorporate spatial descriptors into the CBCS method as a way to extend and improve our previous experiments on spatialized CBCS [16–18]. More specifically, we wish to apply CBCS to a corpus of musical instrument recordings, taking their radiation patterns into account as supplemental descriptors for the organization of the database.

1.2 Previous work

To incorporate radiation patterns in the context of CBCS, as described above, requires a rigorous comparison, or distance function between patterns. However, radiation patterns are in essence multidimensional data, and comparing them is therefore not an obvious task.

One way to assess their resemblance is by evaluating their spatial correlation, defined as follows. A cross-correlation coefficient is a frequency-dependent scalar value that expresses the similarity between sound pressures at two given positions in space [19]. This coefficient can then be integrated over the entire discrete set of measurement sensors for two radiation functions: the result is their spatial correlation. Looking chronologically at examples of this technique, Moreau [20, 21] employs such a spatial correlation for the qualitative assessment of a directivity model of spherical microphone arrays against measured data. Pollow [22] (Sects. 3.3.1–3.3.5) uses spatial correlation for the analysis of instrumental radiation patterns from the Technical University of Berlin (TU Berlin), discussed in further detail below. He produces scatter plots of the correlation values of partial tones radiated by the instruments, similar to what we present in Section 5 using a different algorithm. Sridhar et al. [23] use spatial correlation to compare polar responses of different loudspeakers and assess their invariance over a specified frequency range.

However, spatial correlation presents disadvantages that motivated us to pursue a more flexible and efficient method based on spherical harmonics. As radiation patterns are acoustic wave fields evaluated on the surface of a sphere, they can be conveniently represented by their spherical harmonic coefficients, using a decomposition of the angular functions into the orthonormal set of spherical harmonics [24]. This representation has already proven useful for the analysis [5], reproduction [9], interpolation [25], extrapolation [26], or auralization [27] of directivity functions.

The spherical harmonic representation can also be used to evaluate the correlation of two signals on the sphere; this will be referred to as spherical correlation. Spherical correlation offers some advantages over spatial correlation as presented above: (a) it is independent of the acquisition setup, i.e. from the number and positions of measurement microphones; this potentially allows for comparisons of data from different origins. (b) As spherical harmonic expansions can be easily rotated, spherical correlation can be used to determine the angular displacement (or lag) needed to align the two signals on the sphere. This is widely used in the field of image processing in order to perform 3-D shape-matching, also known as shape registration [28–33]. (c) As spherical harmonic expansion offers the most compact representation of radiation data for a given resolution, spherical correlation can be computed efficiently.

The spherical correlation formalism is discussed in depth in the mathematical literature. The state of the art is introduced in [34], building on the foundation of the spatial Fourier transform [35]. However, while “conventional” correlation is ubiquitous for the analysis of audio signals, so far the practical use of spherical correlation in the field of acoustics remains relatively limited and challenging.

In a chronological review of acoustics applications, spherical correlation is briefly presented in [5, 36], but actually unused. It is used by Guillon [37, 38] to compute similarities between spatial frequency response surfaces (SFRS), and later to clusterize a dataset of SFRS; by Deboy and Zotter [39] to perform rotational tracking of a moving trumpet, for which they address the question of discretizing the 3-D rotation group; and by Pollow et al. [26] to measure the quality of range extrapolation of head-related transfer functions, independent of a possible gain mismatch. In a recent publication, Pezzoli et al. [40] propose a set of metrics, including the spherical correlation coefficient, to compare the sound radiation of several historical violins. Similar to our study, Hohl and Zotter [41] use spherical correlation as a similarity measure of radiation patterns of musical instruments. They examine whether different partials at the same frequency, but originating from different played pitches, exhibit similar radiation on a given instrument. This is promising work, but their short paper does not provide much detail.

1.3 Source data and methods

Our study is both motivated and facilitated by the availability of existing datasets of radiation patterns of orchestral instruments and voices that have been obtained by previous researchers. These studies use acoustical measurements, in anechoic conditions, with surrounding circular or spherical microphone arrays. Instruments are excited either by human players or, in some cases, by electromechanical devices such as an artificial mouth for brass instruments. Many authors have reported such measurements either using synchronous multichannel recordings [3, 6, 42–52] or repeating the signals and rotating the instrument [48, 49, 53–57]. In the case of a human player or singer, the reproducibility of the signals is of course a matter of concern.

Throughout this article, we use a database published by TU Berlin [50, 58] that constitutes, to date, one of the most comprehensive datasets of publicly available 3-D high-resolution synchronous measurements, with radiation patterns for 41 modern and historical orchestral instruments.

In the following sections, we review the theory of spherical correlation and expand on the studies reviewed above to present applications to the analysis of directivity functions. We argue that the spherical correlation coefficient, as a simple scalar number, is a powerful tool for the comparison and similarity measurement of radiation patterns. After introducing the mathematical preliminaries (Sect. 2), we discuss practical considerations for an actual implementation, in particular with respect to the discretization of the search space for rotational matching (Sect. 3). In Section 4, we present numerical simulations to elicit basic properties of spherical correlation and to validate our implementation. We then present several visualization techniques, based on a matrix of cross-correlation values, that can be helpful to detect similarities among radiation patterns, and we propose the use of multidimensional scaling to clusterize radiation data. We apply these tools to several case studies of measured directivity data of musical instruments (Sect. 5), and we show that spherical correlation, both with and without rotational matching, can unveil interesting radiation similarities across frequencies or partials of a given instrument, or across multiple instruments.

2 Theoretical background

2.1 Cross-correlation on the sphere

Given two shapes f and g, we wish to find the find the rotation that best aligns them. Let f be a square-integrable function on the unit sphere $f\in L^2\left({\mathbb{S}}^2\right)$. In the simple case, let us assume that g is a rotated version of f, i.e. f = Λ_R(g) for some 3-D rotation R. We denote Λ_R the rotational operator defined such that $\forall \mathbf{\Omega }\in {\mathbb{S}}^2$

$$\begin{array}{cc}\mathrm{\Lambda }:& L^2\left({\mathbb{S}}^2\right)\to L^2\left({\mathbb{S}}^2\right)\\ & {\mathrm{\Lambda }}_{\mathbf{R}}\left(g\right)\left(\mathbf{\Omega }\right)=g\left({\mathbf{R}}^{-1}\left(\mathbf{\Omega }\right)\right). \end{array}$$(1)

We wish to find the rotation R. In the more general case, given the two patterns f and g, we wish to find the rotation that best aligns the two shapes on the sphere. This can be accomplished by evaluating the cross-correlation between the two functions

$$C_{\mathbf{R}}\left(f,g\right)={\int }_{{\mathbb{S}}^2}\mathrm{}f\left(\mathbf{\Omega }\right) \overline{{\mathrm{\Lambda }}_{\mathbf{R}}\left(g\right)\left(\mathbf{\Omega }\right)} d\mathbf{\Omega },$$(2)

and finding the rotation R that maximizes the above integral [34]. However, evaluating C_R(f, g) for all possible rotations in the spatial domain is a time-consuming task. Instead, we undertake it in the spatial Fourier domain. Since f and g are square-integrable on the sphere ${\mathbb{S}}^2$, we can write their Fourier expansions [35]

$$f\left(\mathbf{\Omega }\right)=\sum _{n=0}^N\mathrm{} \sum _{m=-n}^{m=+n}\mathrm{} f_n^m Y_n^m\left(\mathbf{\Omega }\right).$$(3)

$Y_n^m\left(\mathbf{\Omega }\right)$ are the spherical harmonic functions, and $f_n^m$ are the Fourier coefficients of f. Here we have further assumed that f and g are bandlimited, with their bandwidth B = (N + 1), N being the maximum order of the Fourier expansion.

A well-known property of the conventional (Euclidian) Fourier transform is that a translation in the time domain is interpreted as a phase shift in the frequency domain. For the Fourier transform on ${\mathbb{S}}^2$, this property means that the magnitudes of the Fourier coefficients are invariant under rotation: writing f_n(Ω) for the nth frequency component of f,

$$f_n\left(\mathbf{\Omega }\right)=\sum _{m=-n}^{m=+n}\mathrm{} f_n^m Y_n^m\left(\mathbf{\Omega }\right),$$(4)

the quantity |f_n(Ω)| is invariant under rotation.

Furthermore, the spherical harmonic basis functions of each order $n$ transform among themselves under rotation according to

$${\mathrm{\Lambda }}_{\mathbf{R}}\left(Y_n^m\left(\mathbf{\Omega }\right)\right)=\sum _{m\mathrm{\prime }=-n}^{m\mathrm{\prime }=+n}\mathrm{} Y_n^m\left(\mathbf{\Omega }\right)D_{mm\mathrm{\prime }}^n\left(\mathbf{R}\right),$$(5)

where $D_{mm\mathrm{\prime }}^n\left(\mathbf{R}\right)$ is the Wigner-D function (see Appendix). Equation (5) is valid for complex-valued spherical harmonics. When considering real-valued spherical harmonics (as typically used in the field of Ambisonics research), a similar result holds, however involving the Wigner-d function instead of Wigner-D; details can be found for instance in [59] (Eqs. (16)–(20)).

Now, using equations (3)–(5), it is possible to simplify equation (2) into:

$$C_{\mathbf{R}}(f,g)=\sum _{n=0}^N\mathrm{}\sum _{\left|m\right|\le n}\mathrm{}\sum _{\left|m\mathrm{\prime }\right|\le n}\mathrm{}{f}_n^{-m}\overline{g_n^{-m\mathrm{\prime }}}(-1{)}^{m-m\mathrm{\prime }}{D}_{mm\mathrm{\prime }}^n(\mathbf{R}).$$(6)

A complete demonstration can be found for instance in [34]. Equation (6) allows us efficiently to compute the cross-correlation by combining the Fourier coefficients of f and g. From this, we can now search for which rotation R maximizes C_R(f, g).

In equations (5) and (6), the Fourier coefficients are rotated by means of the explicit formulae with the Wigner-D functions (or Wigner-d for real-valued harmonics); in our practical implementation, we rather use recurrence relations as this appears to be computationally more efficient and numerically stable [60–62].

2.2 Normalized cross-correlation

The normalized cross-correlation [31, 38, 63] is simply a variant of equation (2) normalized by the energy of f and g, written

$${\mathcal{C}}_{\mathbf{R}}\left(f,g\right)=\frac{{\int }_{{\mathbb{S}}^2}\mathrm{}f\left(\mathbf{\Omega }\right)\overline{{\mathrm{\Lambda }}_{\mathbf{R}}\left(g\right)\left(\mathbf{\Omega }\right)} d\mathbf{\Omega }}{\sqrt{{\int }_{{\mathbb{S}}^2}\mathrm{}\left|f\left(\mathbf{\Omega }\right){|}^2 d\mathbf{\Omega } \cdot {\int }_{{\mathbb{S}}^2}\mathrm{}\right|g\left(\mathbf{\Omega }\right){|}^2 d\mathbf{\Omega }}}.$$(7)

The numerator of this expression has already been developed in equation (6). The denominator can be easily calculated thanks to Parseval’s identity:

$${\int }_{{\mathbb{S}}^2}\mathrm{}\left|f\right(\mathbf{\Omega }){|}^2 d\mathbf{\Omega }=\sum _{n=0}^N\mathrm{} \sum _{m=-n}^{m=+n}\mathrm{} |f_n^m{|}^2.$$(8)

As expressed above, the cross-correlation coefficient ${\mathcal{C}}_{\mathbf{R}}(f,g)$ corresponds to a cosine similarity measure. It is also possible to construct a Pearson correlation coefficient by replacing f and g with $\breve{f}$ and ğ, where $\breve{f}$ results from centering f on its spatial average. In the spherical harmonic domain, the spatial average is simply given by the 0th-order component $f_0^0$. This approach is used, for example, in [37, 38, 40]. In [64], this is also referred to as the normalized cross-covariance coefficient.

3 Finding the optimal rotation R

A rotation in ${\mathbb{R}}^3$ can be equivalently represented by: (a) a 3 × 3 rotation matrix, (b) a set of 3 Euler angles (or alternatively 3 Tait–Bryan angles), (c) a unit quaternion, also known as a versor, or (d) an axis–angle representation (i.e. a unit vector indicating the direction of an axis of rotation, and an angle describing the magnitude of the rotation about the axis). The different properties of these formalisms – such as compactness, numerical stability, computational cost, singularity or gimbal lock, etc. – can easily be found in the literature, as well as conversion formulae from one representation to another [65–67].

We denote as SO(3) the 3-D rotation group, also known as the special orthogonal group. SO(3) contains all rotations R about the origin in 3-D Euclidean space: $\mathrm{SO}\left(3\right) = \{\mathit{R} \in {\mathbb{R}}^3|\mathit{RR}T = {\mathit{I}}_3$ and det(R) = +1}, where I₃ is the 3 × 3 identity matrix. In other words, SO(3) is the subgroup of orthogonal matrices with determinant +1.

We must explore the SO(3) space in order to determine the rotation R ∈ SO(3) that maximizes equation (6) or (7) correlation. One possibility would be to use a gradient descent technique, which has been formalized on the SO(3) rotation group in [68, 69]. However, we have decided not to use this approach, as it presents several challenges: the procedure requires approximation of the directional derivatives, with an adjustable small step size, in order to converge to the nearest local maximum. It might also be necessary to restart the algorithm from several initial values, and choosing the optimal step size or the initial value is not a trivial task. Therefore, in the scope of this paper, we instead proceed with a sampling-based exploration of the SO(3) space, as this is conceptually simpler and straightforward in its implementation.

3.1 Sampling the rotation group SO(3)

Due to the topology of SO(3), different choices of parametrization may yield distributions of samples with various properties or biases. To implement our approach, we would prefer a uniform and deterministic sampling scheme. The definition of uniform is subject to interpretation, but essentially the sampling grid should ensure both global coverage and local separation. Note that sampling the SO(3) space is not the same task as uniformly distributing samples on the sphere ${\mathbb{S}}^2$, which is a well-studied question in the field of spherical acoustics processing [70]. In contrast, the SO(3) sampling problem is examined for instance in [71–74]. While an extensive analysis is beyond the scope of this article, in the following subsections we discuss and evaluate the most relevant approaches.

3.1.1 Parametrization by Euler angles

For the sake of simplicity, we will follow [34, 37] and use regular sampling in terms of ZYZ Euler angles (α, β, γ):

$${\alpha }_j={\gamma }_j=\frac{2\pi j}{2B} \mathrm{and} {\beta }_k=\frac{\pi \left(2k+1\right)}{4B},$$(9)

with 0 ≤ j, k < 2B. According to [34], this sampling scheme is suitable for the analysis of band-limited functions with bandwidth B. In this case, the size of the search space is 2B × 2B × 2B.

Note that the null rotation (α = β = γ = 0) is not part of the sampling grid. Consequently, correlating a signal with itself will not yield α = β = γ = 0, but rather β = π/(4B). This is counter-intuitive but not problematic. To circumvent this, it has been proposed [37] first to rotate one of the patterns, f or g, by −π/(4B) around the y axis.

We should also note that the true rotation R might not be on the sampling grid. Therefore, we might only find an approximate solution $\widehat{\mathbf{R}}$ that should be close to R. The notion of distance in SO(3) will be discussed in Section 3.1.4.

3.1.2 Sampling by Halton sequences

Another sampling strategy is proposed in [74], extending and improving a probabilistic approach first proposed in [75, 76] with the use of Halton sequences in the unit cube. This approach is easy to implement, and allows the generation of sampling grids with arbitrary numbers of samples.

3.1.3 Axis-angle visualization

Any 3-D rotation can be represented by an axis-angle representation. This formalism is useful for visualizing the projective space in ${\mathbb{R}}^3$ [72]: each rotation is drawn as a vector with direction n (the axis of the rotation) and a magnitude corresponding to ϑ (the rotation angle). In Figure 1, we present several sampling grids generated via (a) Euler parametrization, (b) Halton sequences, or (c) the Hopf fibration sampling proposed in [72]. The later is another, more elaborated, deterministic strategy that produces dense and highly uniform grids. It can be observed graphically in Figure 1 that Euler sampling is not perfectly uniform, while the Hopf fibration grid divides the surface of SO(3) into regions of (apparently) equal volume. However, the Hopf fibration approach is restricted to fixed sequences of samples, i.e. it cannot generate an arbitrary number of samples. Therefore, for the remainder of this article, we will use the basic Euler sampling (Eq. (9)) for its simplicity and scalability (with respect to B).

Figure 1 Visualization of SO(3) sampling grids using the axis-angle representation. The color represents the magnitude ϑ of the rotation (from blue to red). (a) Regular sampling of Euler angles (as in Eq. (9)) with B = 4, leading to (2B)3 = 512 samples; (b) Uniform distribution using Halton sequences [74] with 512 sampling points; (c) Uniform incremental sampling using the Hopf fibration [72] with 576 sampling points.

3.1.4 Distance measure on SO(3)

Various functions for measuring the “closeness” between 3-D rotations have been proposed in the literature. The usual (angular) distance between two rotation matrices R, Q ∈ SO(3) is given by [39, 71, 77] $d_1(\mathbf{R},\mathbf{Q})=\mathrm{acos}\left(\frac{1}{2} \left(\mathrm{tr}\left(\mathbf{Q}{\mathbf{R}}^T\right)-1\right)\right)$. This metric measures the angle of rotation needed to map the transformation R to the transformation Q, or equivalently the angle of rotation associated to the transformation QR⁻¹. Alternatively, we can interpret d₁(R, Q) as a scaled Froebenius norm, since $d_2(\mathbf{R},\mathbf{Q})={\left|\right|\mathbf{R}-\mathbf{Q}\left|\right|}_{\mathrm{F}}=2 \sqrt{2} \sin \left(\frac{d_1(\mathbf{R},\mathbf{Q})}{2}\right)$, where ${\left|\right|\cdot \left|\right|}_{\mathrm{F}}$ is the Frobenius norm defined as ${\left|\right|\mathbf{A}\left|\right|}_{\mathrm{F}}=\sqrt{\mathrm{tr}\left(\mathbf{A}{\mathbf{A}}^T\right)}$. Several other distance functions have also been proposed, producing values in different ranges and of different units. Huynh [78] presents a detailed review and analysis of various metrics, and demonstrates that many of them are functionally equivalent. Huynh concludes that the following metric, based on quaternions, is both spatially and computationally more efficient:

$$\mathrm{d}\left(\mathbf{R},\mathbf{Q}\right)=1-\left|{\mathbf{q}}_{\mathbf{R}}·{\mathbf{q}}_{\mathbf{Q}}\right|,$$(10)

where q_R is the unit quaternion corresponding to matrix R, and · denotes the quaternion inner product (scalar inner product of two 4-D unit vectors). Equation (10) gives values in the range [0, 1], with 0 denoting that R and Q are close [73, 78]. This metric will be used for the remainder of this article.

Note that, while a distance metric in SO(3) is useful, it must be handled with care in the context of this paper: consider for example an axis-symmetric pattern f, for example a cardioid in the y direction. Rotate f around the z-axis with R₀ ≡ (α₀, β₀, γ₀) = (π, 0, 0). Rotate f around the x-axis with R₁ ≡ (α₁, β₁, γ₁) = (0, 0, π). Both scenarios result in the exact same pattern g; however, the two rotations differ, and their distance is d(R₀, R₁) = 1.

4 Numerical simulations

We run numerical simulations on simple test cases in order to elicit basic properties of spherical correlation and to validate our implementation.

4.1 Impact of SO(3) sampling

As discussed above (Sect. 3.1), the choice of the sampling grid for SO(3) might have an influence on the accuracy and efficiency of the maximization problem equation (7) correlation. For the sake of simplicity, here we consider only the regular Euler sampling of equation (9). However, we investigate the use of oversampling i.e. we use a smaller step size in discretizing the Euler angles (α, β, γ). Instead of using 2B samples for each parameter, we use $\mathcal{N}\ge 2B$. This seems relevant as our scenario involves directional functions with relatively low bandwidth (B = 5 if we consider the radiation patterns available in the TU Berlin database) compared to other authors (Kostelec [34] typically presents results with B = 128 or B = 256). With such low bandwidth (B = 5), the angular step size is very large (π/B = 36°), and therefore high angular misalignment might occur. Of course, oversampling the search grid results in higher computation time, but that is not the focus of this paper.

We run the following numerical simulation: (a) generate a random directional function f with a given bandwidth B; (b) generate a random rotation matrix [79] R ∈ SO(3) and compute g = Λ_R(f), a rotated version of f; (c) sample SO(3) and search for the rotation matrix $\widehat{\mathbf{R}}$ that maximizes the cross-correlation; (d) compute the cross-correlation after rotational matching, i.e. the cross-correlation between g and ${\mathrm{\Lambda }}_{\widehat{\mathbf{R}}}\left(f\right)$; (e) repeat the simulation for various samplings of SO(3), varying the oversampling factor $\mathcal{N}/\left(2B\right)$; (f) repeat the simulation for various bandwidths 1 ≤ B ≤ 6. For each test case, we perform 10,000 Monte-Carlo runs. The results are presented in Figure 2.

Figure 2 Impact of oversampling the SO(3) search space.

For each simulated bandwidth B, it can be observed that oversampling the search space improves accuracy, allowing for a higher cross-correlation (closer to 1), and a smaller distance (closer to 0) between the expected and estimated rotation matrices. For B = 2, the distance criteria $\mathrm{d}(\mathbf{R},\widehat{\mathbf{R}})$ appears bounded: this is due to the non-uniqueness of the solution, as several matrices $\widehat{\mathbf{R}}$ achieve rotational matching (with normalized cross-correlation values close to 1), while being “far” from the expected/simulated matrix R (see concluding remark in Sect. 3.1.4).

As a result of this simulation, we speculate that sufficiently oversampled Euler schemes yield accuracy comparable to other highly uniform grids (such as Halton sequences or Hopf fibration).

4.2 Auto-correlation of elementary patterns

We now perform a simulation to highlight some basic properties of the spherical auto-correlation, i.e. the correlation ${\mathcal{C}}_{\mathbf{R}}(f,f)$ of a signal f with a rotated copy of itself. For f we use different elementary functions (see legend of Fig. 3). For the sake of simplicity and visualization, we examine rotations only around the yaw angle: R ≡ R(α). Figure 3 presents the auto-correlation as a function of the rotational lag α.

Figure 3 Top: polar pattern of elementary functions. Bottom: auto-correlation of f as a function of the rotational lag α (yaw angle). Columns (a) to (e): $f\equiv Y_n^n\left(\mathbf{\Omega }\right)$ for 0 ≤ n ≤ 4. Column (f): Dirac delta with N = 4.

We can make a few simple observations in agreement with theoretical expectations: (1) the auto-correlation coefficient ${\mathcal{C}}_{\mathbf{R}}(f,f)$ varies in range [−1; 1]; (2) as the omnidirectional pattern $Y_0^0\left(\mathbf{\Omega }\right)$ is rotationally invariant, its auto-correlation always equals 1; (3) as all examples of f are real functions, the auto-correlation is an even function of α; (4) the auto-correlation reaches its peak at the origin (α = 0), and for any lag α we have: $\left|{\mathcal{C}}_{\mathbf{R}\left(\alpha \right)}\right(f,f\left)\right|\le {\mathcal{C}}_{\mathbf{R}\left(0\right)}(f,f)$; (5) the auto-correlation of a periodic function is, itself, periodic with the same period; (6) the Dirac distribution is an eigenfunction of the auto-correlation function, i.e. the auto-correlation of a Dirac delta is a Dirac distribution itself. All of these are well-known properties of the auto-correlation function of time signals, here pertained to spherical auto-correlation of directional signals on ${\mathbb{S}}^2$.

4.3 Comparing two patterns of different shapes

So far, we have assumed that the two patterns f and g are rotated cousins and therefore have the same bandwidth. In the most general case, we are interested in comparing two arbitrary radiation patterns, potentially with different bandwidths. Such a scenario generally yields $\left|{\mathcal{C}}_{\widehat{\mathbf{R}}}\right(f,g\left)\right|<1$, but the rotational alignment can still be effective.

As a proof of concept, we investigate the cross-correlation function of two patterns f and g having different “shapes”. For f, we choose an Nth-order Dirac delta, i.e. a maximally directional pattern in direction Ω₀:

$$\begin{array}{c}f\left(\mathbf{\Omega }\right)=\sum _{n=0}^N\mathrm{} \sum _{m=-n}^{m=+n}\mathrm{} f_n^m Y_n^m\left(\mathbf{\Omega }\right),\\ \mathrm{with} f_n^m=\frac{4\pi }{(N+1{)}^2} \overline{Y_n^m\left({\mathbf{\Omega }}_{\mathbf{0}}\right)}.\end{array}$$(11)

For g, we build a “directionally reduced” version of $f$ such as:

$$\forall \left(n,m\right), g_n^m=f_n^m w_n\left(\zeta \right),$$(12)

where ζ ∈ [0–100] is a directivity (or “aperture”) factor, and w_n(ζ) is a spherical harmonics weighting function. The precise choice of w_n(ζ) is not relevant here; we simply build a weighting function that allows us to generate a series of patterns with significantly different characteristics. We further impose that ζ = 100% produces a Dirac delta, and ζ = 0% produces an omnidirectional pattern. In other words, w_n(ζ) is used to simulate patterns with varying bandwidths (see Fig. 4). Techniques for building such weighting functions have been proposed e.g. in [80]. Finally, we rotate g with a random matrix R ∈ SO(3), repeating the simulation with 10,000 Monte-Carlo runs. One example of resulting rotational alignment is presented in Figure 5.

Figure 4 Polar representation of angular pattern g(Ω) for N = 4 and for various values of the directivity factor ζ. The radial scale is linear.

Figure 5 Example of rotational alignment of two patterns with different bandwidths. In magenta, the original pattern f (Dirac delta with N = 4); in black, the rotated and directionally reduced pattern g (with ζ = 40%); in red, pattern f after rotational alignment with the estimated matrix $\widehat{\mathbf{R}}$. We observe that the rotation matching is still effective, although the main lobes of the patterns are not perfectly aligned ($\mathrm{d}(\mathbf{R},\widehat{\mathbf{R}})\approx 10^{-1}$). The correlation coefficient is ${\mathcal{C}}_{\mathbf{R}}(f,g)\approx 0.5$.

In Figure 6, we present the cross-correlation of the functions f and g for varying values of ζ. As expected, we observe that the spherical correlation coefficient decreases when ζ tends towards 0, i.e. when g becomes more omnidirectional. When ζ = 0%, the correlation value is not 0, but the curve is bounded by 0.2. Indeed, when g reduces to the omnidirectional pattern, it is easy to show that ∀f:

$${\mathcal{C}}_{\mathbf{R}}\left(f,g\right)=\frac{f_0^0}{\sum _{n=0}^N\mathrm{} \sum _{m=-n}^{m=+n}\mathrm{} |f_n^m{|}^2}.$$(13)

Figure 6 Numerical simulation of correlation between Dirac delta f and directionally reduced g, as a function of the directivity factor ζ, for N = 4. Top: Normalized cross-correlation (after rotational matching) between f and g. Bottom: Distance $\mathrm{d}(\mathbf{R},\widehat{\mathbf{R}})$ between expected and estimated rotation matrices.

When f is a Dirac delta, this further simplifies to ${\mathcal{C}}_{\mathbf{R}}(f,g)=\frac{1}{B}$, which equals 0.2 in our example. The actual shape of the ${\mathcal{C}}_{\mathbf{R}}(f,g)$ curve (Fig. 6) depends on the chosen weighting functions w_n(ζ). Similarly, the distance between the estimated and expected rotation matrix significantly increases as ζ gets smaller.

These results suggest that the spherical correlation coefficient can be a relevant indicator of similarity between patterns of different shapes/bandwidth. As with any correlation measure, this essentially allows for qualitative interpretation; a quantitative analysis of the correlation coefficient depends on the context and purposes.

5 Application to measured data

We now apply the approaches presented in the previous sections to measured data. As mentioned in the introduction, our source is the acoustic instrumental radiation database, containing 41 modern and historical orchestral instruments, made available by TU Berlin [50, 58]. This database contains single-tone recordings at two dynamic levels (pp and ff). For our processing, we use the ff data due to its better signal-to-noise ratio.

In the TU Berlin database, radiation patterns were measured by a surrounding spherical array of 32 microphones with a radius of 2.1 m. The database is published with spherical harmonic coefficients calculated for the first 10 partials of each played note, and third-octave band-averaged patterns are also provided [81]. Methods to obtain the spherical harmonic coefficients from such measurements are not discussed here; readers can refer e.g. to [5, 9, 24, 36, 54]. In our work, we directly exploit the precomputed spherical harmonic coefficients for the order N = 4 (25 coefficients) as provided by the TU database authors. An acoustic source centering algorithm [82, 83] has also been applied, by the TU researchers, to these spherical harmonic coefficients (below 1 kHz) in order to align the acoustic center of the sound source to the geometrical center of the microphone array, and to account for the resulting phase shifts at the microphone positions. This source re-alignment procedure is known to be particularly important for directivity functions that model the complex sound pressure (as is the case in the TU database), with significant impact on the spherical harmonic coefficients; when only the absolute values (sound pressure levels) are considered, the impact of displaced sources might be less severe [84, 85]. The TU database includes both uncentered and centered data, and we have used the latter to minimize the influence on our calculations of changes to instrument alignment. The practical setup used in the collection of the TU data exhibits a spatial aliasing frequency of approximately 1.1 kHz [58, 82]. Observations made above this frequency should therefore be interpreted cautiously.

Finally, let us emphasize that our aim in this section is not to provide a thorough analysis of the TU database, but rather to exemplify the qualitative results that can potentially be obtained through spherical correlation. More in-depth discussions about the TU data can be found e.g. in [22, 58, 86–88].

5.1 Similarities over frequencies

In Figure 7, we propose scatter plots to visualize the matrix of cross-correlation values, i.e. the normalized cross-correlation of all pairs of partial tones of a given instrument, for four different instruments. The color code depicts the magnitude of the correlation coefficient (its absolute value), from blue (low correlation) to red (high correlation). The x and y axes represent the frequency on a logarithmic scale. The size of the scatter dots is slightly adjusted with frequency in order to avoid overlapping dots and to improve readability. This visualization is similar to what is proposed in [22, 89]. Such a diagram allows the compact combination of a significant amount of information: the number of dots is equal to (N_p × N_t)² where N_t is the number of played tones and N_p the number of partials for each tone. In the TU dataset, N_p = 10. As the matrix of cross-correlation values is persymmetric, so is the scatter diagram. The cross-correlation values presented here have been calculated without rotational alignment. When applying rotational alignment, the scatter diagrams (not shown here) do not exhibit significant differences, suggesting the musicians did not rotate dramatically during the recording session.

Figure 7 Matrix of cross-correlation values for all recorded tones, including all partials, organized by frequency. Results are presented for clarinet in B♭, cello, tuba, and double bass.

Inspecting these diagrams, as well as those of other instruments not displayed here for brevity, we can draw several general observations. The results are consistent with the categorization of the TU database into three groups of instruments, as proposed by Shabtai et al. [58]: those with one expected radiation point, such as brass instruments; those with several expected radiation points, such as woodwind instruments, with sound radiated by the bell, the fingering holes, and the mouthpiece; and those with a full body radiating sound, such as string instruments. For brass instruments (tuba in Fig. 7), the cross-correlation values vary slowly with frequency, and the diagonal in the matrix has a wide “spread”. This suggests that the radiation patterns evolve smoothly over frequency. In other words, nearby frequencies produce similar directivity patterns. For woodwind instruments (clarinet in Fig. 7), there is greater correlation of all partials for low frequencies, approximately below 300 Hz. With increasing frequencies, the correlation values exhibit abrupt changes, manifested by vertical and horizontal “stripes” and “clusters” in the depicted matrix containing the cross-correlation values. These discontinuities are likely related to changes of fingering, for example using the register key. For string instruments (cello and double bass in Fig. 7), the cross-correlation is relatively large in the low frequency region (for wavelengths larger than the dimensions of the instrument), but drops sharply at higher frequencies, suggesting greater variations in their radiation patterns. For further interpretations see e.g. [22, 46, 89].

5.2 Similarities over partials

Following an idea from Hohl and Zotter [41, 46], we now plot the matrix of cross-correlations organized by partials rather than by frequencies. In other words, we generate a 2-D color-coded diagram, similar to the one presented in the previous section, but with x and y axes now organized by partial tones instead of frequencies. The matrix is organized in blocks each containing a chromatic scale, with each successive block representing a different partial for each played note. Black dashed lines separate each partial index, for ease of readability. This sorting makes the comparison of cross-correlation between partials more evident. With this matrix arrangement, partials with matching frequency are located on the secondary diagonals. Figure 8 illustrates examples of results for the same four instruments as in Figure 7.

Figure 8 Matrix of cross-correlation values (for all partial tones of one chromatic scale), organized by partial tones.

For the tuba, it can be observed that most partials at the same frequency exhibit strongly correlated radiation, regardless of the played pitch from which they originate. This is in line with the observations made in Section 5.1. Similar results can be observed for other brass instruments (not shown here). For the clarinet, however, the diagram is less regular. Along the diagonal, a checkered pattern appears, revealing that partial tones with similar or nearby frequencies radiate differently. More precisely, pitches from C to G# are strongly correlated, but differ significantly from pitches A to B. The abrupt change between these two blocks corresponds to the transition between the first and second registers of the B♭ clarinet. This suggests that the addition of the register key, and consequently discontinuous change in fingerings, strongly impacts the radiation pattern. For string instruments such as the cello and double bass, except for a few narrow diagonal red lines, the overall similarities between partials are low. Each partial seems to have its own radiation pattern, uncorrelated with the other partials, suggesting that the directivities of these instruments are rather complex. The diagrams presented in Figure 8 are qualitatively comparable to ones obtained by Hohl and Zotter with another database of instrument recordings. Further analysis can be found in [41, 46].

5.3 Visualization with multidimensional scaling

Another motivation for this work is to classify the TU Berlin database for subsequent manipulation with corpus-based synthesis techniques: by organizing instrumental samples according to their directional characteristics, the resulting low-dimensional representation can be used, along with other audio descriptors, to navigate the sample corpus. With an approach to pairwise spherical correlation in place, we can now examine applications to classify and navigate larger collections of 3-D radiation patterns by using multidimensional scaling (MDS) [90]. MDS is used to translate information about the pairwise distances (or dissimilarities) among a set of p objects into a configuration of p points mapped onto an abstract Cartesian space. MDS is therefore a means of visualizing the level of similarity of samples within a dataset.

We apply MDS analysis, using a dimensionality of 2 for simplicity of visualization and interpretation, to measured radiation patterns of different instruments for the fundamental frequency of various played pitches. Examples of MDS for two different pitches are presented as scatter plots in Figure 9. The color code corresponds to the three categories of instruments cited above [58]: in red, instruments with one expected radiation point such as brass; in blue, instruments with several expected radiation points, such as woodwinds; and in black, instruments with a complex radiation pattern, such as strings. The MDS plots for these two pitches, as well as for others not shown for brevity, show groupings among instruments of these three categories. Further sub-classes of instruments of similar construction appear to cluster together: for example, the cylindrical-bore brass instruments in the trombone family, which are separated from conical-bore French horn and tuba. Members of the single-reed clarinet and saxophone families cluster together, as do members of the double-reed bassoon and dulcian family. As expected, historical instruments are positioned near their modern counterparts. Therefore we observe that MDS based on spherical correlation allows us efficiently to segregate and organize the instruments in accordance with their predicted categories.

Figure 9 2D MDS map of radiation patterns of the pitches B3 (a) and E3 (b) for various brass (red), woodwind (blue), and string instruments (black).

The examples presented in Figure 9 were evaluated with pitches performed at dynamic ff. When both dynamic levels pp and ff are included in the MDS analysis, the corresponding points (not shown for clarity) are mapped closely in the MDS space, indicating that the played dynamic level does not substantially affect the clustering operation.

Rotational matching has not been applied to these examples in order better to visualize the distinctions between instrument categories – as, by definition, rotational alignment always yields higher correlations. We examine the effects of rotational matching in the following section.

5.4 Exploring rotational matches

We now show that rotational matching can be a convenient tool to explore the database further and search for similarities, regardless of the orientation of the instruments. The following analyses refer to the radiation coefficients averaged over third-octave bands [58]. While both pp and ff dynamics have been analyzed, the results are presented only for ff data for the sake of readability.

5.4.1 Across instruments

In Figure 10, we examine the cross-correlation between all instruments in a single frequency band, before and after applying rotational matching. The axes are organized according to the expected three main categories of radiation characteristics, separated by black lines (“brass-like”, “woodwind-like”, and “string-like”, as discussed in Sect. 5.3). Within each category, instruments are presented in alphabetical order without presupposing any finer categorization. White spaces indicate instruments for which no data is present in the chosen frequency band (timpani, flutes, and oboes in Fig. 10). As the matrix of cross-correlation values is persymmetric, only half of it is shown: instead, in the upper and lower triangles we present the correlation values with and without rotational alignment, respectively.

Figure 10 Matrix of cross-correlation values between instruments for the third-octave band centered at 198 Hz. Lower triangle: without rotational matching; Upper triangle: after rotational matching.

We can make a number of observations, ranging from more general to more specific:

1. The upper triangle has higher values than the lower triangle, by definition of rotational matching.

2. In the upper triangle, nearly all correlation values are greater than 0.5, indicating some degree of similarities across most instruments. This is somewhat expected: in this relatively low frequency band, most instruments have energy concentrated in the lower spherical harmonic orders, and they behave roughly like first-order functions, with cardioid- or dipole-like patterns. The concentration of energy within the lower orders is also, partly, a consequence of the acoustic centering process that has been applied [82].

3. The three “diagonal blocks”, representing correlation among instruments of the same category, show overall higher correlation than the other blocks, comparing instruments of unlike categories. This is observed both before and after rotational matching, but with several nuances: for example, the “brass-like” instruments show a particularly large increase in correlation after rotational matching, which is consistent with the assumption that they each have a single point of radiation but oriented in different directions. To a lesser extent, the “string-like” instruments show a significant improvement in correlation following rotational matching, which again fits the observation that several of these instruments feature a similar construction but oriented in different directions (for example violin, viola, cello, and guitar). In these cases, therefore, the similarities among instruments are most clearly revealed when they are rotated.

4. A few instruments visually emerge as “outliers”, exhibiting a line (vertical or horizontal) with strikingly low correlation values, both with and without rotational matching. In particular, this is the case for the double action harp, pedal timpani, French horn, and to a lesser extent dulcian. Indeed, it appears that these instruments have significant energy in higher orders (N > 1), and consequently their radiation patterns are dissimilar to other instruments. We can propose several hypothetical explanations: (a) These instruments are also outliers in terms of their construction and performance technique, with few obvious correlates in the database: the timpani and harp are percussion instruments that are only provisionally grouped with “strings” due to their full-body radiation. The horn, while technically a brass instrument, has historically been acknowledged as an exceptional member of that group due to its timbre and directivity [2]: it is the only instrument whose bell is oriented behind the performer, with the performer’s hand mediating the output. Its directional characteristic is therefore expected to be more complex (of higher order) than other brass instruments performed with open bells. While the dulcian may be expected to correlate with other members of the bassoon family, it is an outlier as the oldest instrument of the database by nearly a century (ca. 1600) [50] and we hypothesize that its unusual directivity is consistent with historical trends favoring more focused directivity as instruments have modernized. (b) Consistent with these observations, the presence of higher order components may be caused by the contribution of large, spatially extended, excitation sources (i.e. the absence of a unique natural acoustic center) and/or by off-center sound sources and aliasing errors; this seems especially plausible for the harp and timpani, as they are among the largest instruments in the database. We hypothesize that the performance of the centering algorithm might be degraded for these instruments (in spite of the relatively low frequency). (c) When measuring harp and timpani, “a relatively small-area floor was used” to support the instrument [83]; for the wavelength of interest, this might have contributed to undesired reflections or scattering (based on the available pictures, this floor was covered with absorptive materials for the timpani only). (d) Finally, the timpani are unique in the database as only one played pitch was recorded for each (historical and modern); this could explain differences from the other instruments, due to relatively fewer data for the third-octave averaging.

5. As already noted in Section 5.3, modern instruments and their historical counterparts usually exhibit high correlation values, even without rotational matching. The best rotational alignment is typically achieved with a relatively small angular displacement (not shown here for brevity), suggesting that the difference between modern and historical instruments might be essentially due to slight changes in the performer’s orientation and instrument construction. It is, however, not obvious why modern and historical viola exhibit low correlation values in this particular frequency band.

6. Double bass exhibits only moderate correlation with other members of the string quartet; this bias is possibly an effect of its large size relative to the frequency band being considered: we have noted in Section 5.1 that, for string instruments, partial tones within a small frequency interval show low correlation values for wavelengths smaller than the size of the instrument; consequently, the averaged directivity pattern of Double bass in the frequency band centered at 198 Hz is questionable.

7. While it is surprising there is not higher correlation between French horn and natural horn, the latter may be anomalous because relatively few performed pitches in the natural harmonic scale leave sparse data for third-octave averaging in this register (similar to timpani as mentioned above).

8. We observe high correlation between several unexpected pairs of instruments before rotational matching, such as trumpet and historical violin, tuba and historical violin, or trumpet and clarinet. There is not space to analyze every example in this paper; however, spherical correlation and the proposed visualization strategies appear as promising tools for further research. For brevity we also have not systematically discussed here the specific angle of rotation that produces the optimal rotational alignment; again, this would be a fruitful topic for further study.

5.4.2 Across frequencies

As mentioned in observation 3, rotational alignment is especially relevant for brass, so we will focus on the tenor trombone as one particular representative. Still considering its data in the 198 Hz-centered band, we present in Figure 11 its cross-correlation with all other instruments in all available frequency bands.

Figure 11 Matrix of cross-correlation values between tenor trombone in the third-octave band centered at 198 Hz, and all available instruments. (a) Without rotational matching; (b) After rotational matching.

Before rotational matching, we observe high correlation values with the other trombones at the same frequency (historical tenor trombone, modern and historical bass trombone), but also with the bass trombone in nearby frequency bands (99 Hz and 157 Hz). This expected affinity between all trombones in nearby frequencies is further increased after applying rotational alignment. There is also some correlation with the trumpet in the 314 Hz band, where the shift to higher frequencies could be reasonably explained by the trumpet’s smaller dimensions.

Rotational matching additionally reveals high similarities with certain woodwind instruments, such as the bass clarinet and bassoons, in slightly lower frequency regions. In Figure 11b, the highest correlation values (≥0.95) are observed for: bass trombone (99 Hz and 198 Hz), bassoon (79 Hz), contrabassoon (63 Hz), bass clarinet (99 Hz), and historical bass trombone (99 Hz).

The three instruments that “benefit” the most from rotational alignment are contrabassoon, bass clarinet, and historical bass trombone. One hypothesis is that these instruments have a relatively similar size, but are performed with their bells oriented in different directions: the trombone points forwards in front of the performer, while the bassoon and contrabassoon bells point upwards, and the bass clarinet bell is located at the bottom of the instrument but pointing upwards. This is partially confirmed by Shabtai et al. [58] (Figs. 8a and 9a) in their measurements of acoustic source centers, where they observe that both trombone and bassoon have strong excitation sources at their bells, but with different spatial orientations. Therefore it is plausible that a suitable rotation would produce a high correlation between the radiation patterns of these instruments.

More surprisingly, there is also noticeable correlation, before and after rotational matching, with a few strings instruments (violas, basses, historical cello, and acoustic guitar) in certain frequency bands. As in Section 5.4.2, observation 2, this could be due to the relatively low-order radiation patterns of all instruments in low frequency bands.

In summary, beyond the expected correlations with other brass, Figure 11 reveals that rotational matching can capture less obvious similarities with different instrument categories, such as some woodwinds and strings, in different frequency regions. Nevertheless, a much more systematic study would be required to confidently relate these observations to the instruments’ construction and performance technique.

Now focusing on one of the above examples, in Figure 12 we graphically represent the radiation patterns for tenor trombone and contrabassoon in the respective frequency bands centered at 198 Hz and 63 Hz, where rotational alignment induced a large increase in cross-correlation. As the normalized cross-correlation is used here, it is possible to detect similarities regardless of differences in overall energies of the two instruments, as already noted by Pollow [22] (Sect. 2.3.7). The contrabassoon radiation pattern shown in Figure 12c is essentially a vertical dipole, consistent with the vertical orientation of the instrument, as compared to the horizontal orientation of the trombone. However, after rotating the contrabassoon radiation pattern, in Figure 12f we can observe that the 2-D polar patterns align nicely. This is not always the case: indeed, maximizing the 3-D shape alignment does not necessarily imply strong similarity on the horizontal plane.

Figure 12 Example of 3-D rotational matching. (a) and (b): tenor trombone for the third-octave band centered at 198 Hz. (c) and (d): contrabassoon for the third-octave band centered at 63 Hz. (e): contrabassoon after rotational matching. (f) black curve: tenor trombone; similar to (b). (f) solid red curve: contrabassoon after rotational matching. (f) dashed red curve: similar to solid curve, scaled to match the overall energy of the tenor trombone. Left: visualization of sound pressure (in dB) on the sphere. Right: restriction to the horizontal plane (radial scale is in dB).

Beyond acoustical analysis, such high correlations suggest clear applications from an artistic perspective in the case of CBCS, as discussed in the introduction. By concatenating instrumental samples according to spherical correlation, smooth transitions can be made between radiation characteristics of disparate instruments, with or without rotational alignment, suggesting a novel approach to corpus-based spatial sound synthesis.

6 Conclusion

In this paper, we discussed the use of cross-correlation on the sphere as a similarity measure for the classification of 3-D radiation patterns. We showed that this tool can facilitate thorough analysis of directivity pattern similarities across partials of one instrument or between different instruments. We have presented several visualization tools that allow the compact organization and examination of multidimensional data, revealing or confirming the radiation behavior or categorization of instruments. Multidimensional scaling, in particular, appears to be a powerful approach to clusterize instruments based on their radiation characteristics, with potential applications to the creative exploration of the corpus. We have also discussed some challenges regarding the implementation of rotational alignment, and demonstrated that rotational matching can be useful to detect similarities among different categories of instruments, or different frequency bands, by mitigating the effects of differing instrument construction and orientation.

As a consequence, spherical correlation presents promising possibilities for human or machine navigation of a database of radiation patterns. We have already applied CBCS to computer improvisation based on audio features, in particular using timbral descriptors to structure sequences of complex sounds with an audio oracle algorithm [91]. These CBCS sequences can be synthesized spatially using spherical harmonic coefficients derived from the TU database [18] and projected with a compact spherical loudspeaker array such as the IKO [92] to approximate the dynamic radiation patterns of acoustic instruments, as implemented in Einbond’s composition Prestidigitation for percussion and 3-D electronics in 2022. A further step will be to use an MDS visualization of the TU database to train the computer improvisation agent directly on spherical correlation distances themselves, allowing for direct learning and continuation of spatial gestures.

Future work will investigate other spatial descriptors for the classification of the 3-D database of orchestral instruments as well as other metrics recently proposed in [40]. Finally, we will examine whether spherical correlation can be useful to detect or correct possible rotational misalignment of a human performer across multiple measurements, thereby improving reproducibility.

Acknowledgments

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101019164 (ERC MusAI – Music and Artificial Intelligence: Building Critical Interdisciplinary Studies).

Data availability statement

The Matlab code associated with this article is available in https://github.com/tcarpent/shxcorr

Appendix

A.1 Wigner-D function

Considering Euler angles with the ZYZ convention, the rotation matrix R can be expressed as a rotation of α about the z axis, followed by a rotation β about the y axis, and finally a rotation γ about the z axis:

$$R_{\mathrm{zyz}}\left(\alpha ,\beta ,\gamma \right)=R_z\left(\alpha \right) R_y\left(\beta \right) R_z\left(\gamma \right)$$(14)

with:

$$R_z\left(\alpha \right)=\left(\begin{array}{lll}\cos \alpha & -\sin \alpha & 0\\ \sin \alpha & \cos \alpha & 0\\ 0& 0& 1\end{array}\right),$$(15)

$$R_y\left(\beta \right)=\left(\begin{array}{lll}\cos \beta & 0& \sin \beta \\ 0& 1& 0\\ -\sin \beta & 0& \cos \beta \end{array}\right),$$(16)

$$R_z\left(\gamma \right)=\left(\begin{array}{lll}\cos \gamma & -\sin \gamma & 0\\ \sin \gamma & \cos \gamma & 0\\ 0& 0& 1\end{array}\right),$$(17)

and 0 ≤ α, γ < 2π and 0 ≤ β ≤ π. With this convention, the Wigner-D function $D_{mm\mathrm{\prime }}^n\left(R\right)$, required to rotate complex-valued spherical harmonics, is written

$$D_{mm\mathrm{\prime }}^n\left(\alpha ,\beta ,\gamma \right)=e^{-\mathrm{im\alpha }} d_{mm\mathrm{\prime }}^n\left(\beta \right) e^{-im\mathrm{\prime }\gamma },$$(18)

where $d_{m{m}^{\text{'}}}^n\left(\beta \right)$ is the Wigner-d function [34, 93].

References

All Figures

Figure 1 Visualization of SO(3) sampling grids using the axis-angle representation. The color represents the magnitude ϑ of the rotation (from blue to red). (a) Regular sampling of Euler angles (as in Eq. (9)) with B = 4, leading to (2B)3 = 512 samples; (b) Uniform distribution using Halton sequences [74] with 512 sampling points; (c) Uniform incremental sampling using the Hopf fibration [72] with 576 sampling points.

In the text

	Figure 2 Impact of oversampling the SO(3) search space.
In the text

	Figure 3 Top: polar pattern of elementary functions. Bottom: auto-correlation of f as a function of the rotational lag α (yaw angle). Columns (a) to (e): $f\equiv Y_n^n\left(\mathbf{\Omega }\right)$ for 0 ≤ n ≤ 4. Column (f): Dirac delta with N = 4.
In the text

	Figure 4 Polar representation of angular pattern g(Ω) for N = 4 and for various values of the directivity factor ζ. The radial scale is linear.
In the text

Figure 5 Example of rotational alignment of two patterns with different bandwidths. In magenta, the original pattern f (Dirac delta with N = 4); in black, the rotated and directionally reduced pattern g (with ζ = 40%); in red, pattern f after rotational alignment with the estimated matrix $\widehat{\mathbf{R}}$. We observe that the rotation matching is still effective, although the main lobes of the patterns are not perfectly aligned ($\mathrm{d}(\mathbf{R},\widehat{\mathbf{R}})\approx 10^{-1}$). The correlation coefficient is ${\mathcal{C}}_{\mathbf{R}}(f,g)\approx 0.5$.

In the text

	Figure 6 Numerical simulation of correlation between Dirac delta f and directionally reduced g, as a function of the directivity factor ζ, for N = 4. Top: Normalized cross-correlation (after rotational matching) between f and g. Bottom: Distance $\mathrm{d}(\mathbf{R},\widehat{\mathbf{R}})$ between expected and estimated rotation matrices.
In the text

	Figure 7 Matrix of cross-correlation values for all recorded tones, including all partials, organized by frequency. Results are presented for clarinet in B♭, cello, tuba, and double bass.
In the text

	Figure 8 Matrix of cross-correlation values (for all partial tones of one chromatic scale), organized by partial tones.
In the text

	Figure 9 2D MDS map of radiation patterns of the pitches B3 (a) and E3 (b) for various brass (red), woodwind (blue), and string instruments (black).
In the text

	Figure 10 Matrix of cross-correlation values between instruments for the third-octave band centered at 198 Hz. Lower triangle: without rotational matching; Upper triangle: after rotational matching.
In the text

	Figure 11 Matrix of cross-correlation values between tenor trombone in the third-octave band centered at 198 Hz, and all available instruments. (a) Without rotational matching; (b) After rotational matching.
In the text

Figure 12 Example of 3-D rotational matching. (a) and (b): tenor trombone for the third-octave band centered at 198 Hz. (c) and (d): contrabassoon for the third-octave band centered at 63 Hz. (e): contrabassoon after rotational matching. (f) black curve: tenor trombone; similar to (b). (f) solid red curve: contrabassoon after rotational matching. (f) dashed red curve: similar to solid curve, scaled to match the overall energy of the tenor trombone. Left: visualization of sound pressure (in dB) on the sphere. Right: restriction to the horizontal plane (radial scale is in dB).

In the text