SYSTEM AND METHOD FOR ESTIMATING DIRECTION OF ARRIVAL AND DELAYS OF EARLY ROOM REFLECTIONS

Abstract

Aspects of embodiments pertain to systems for estimating direction of arrival and delays of reflected and non-reflected sources in a room received by a plurality of sound capturing devices in the room, comprising obtaining sound data descriptive of received sources. Embodiments may further comprise determining, based on the sound data, a correlation value between each two audio channels of a plurality of audio channels to obtain, for a same frequency of the plurality of frequencies, a set of correlation values; and to perform a transform on each set of frequency-related correlation values to obtain a plurality of subsets of delay representations of a 3D Matrix. Each set of time delay representations may be analyzed or processed to extract information about the sources.

Description

BACKGROUND

Estimation of the direction of arrival (DOA) and delay of room reflections from reverberant sound may be useful for a wide range of applications such as speech enhancement, dereverberation, source separation, optimal beamforming and room geometry inference.

Early reflections have a key role in sound perception, as they can improve speech intelligibility and listener envelopment. They are also related to the impression of source width, loudness and distance. Therefore, exploitation of the early reflections can be beneficial in parametric spatial audio methods and spatial audio coding.

Existing methods for the estimation of the parameters of early reflections can be categorized as blind and non-blind. Non-blind methods, operate on room impulse response signals, or, alternatively, assume that an anechoic recording of the sound source is available. Blind methods operate on microphone signals directly, and assume that no other information is available, as is often the case in many real world applications.

Such methods employ spatial filtering, subspace methods (e.g., MUSIC or ESPRIT), sparse recovery methods which attempt to find the smallest subset of sources that explains the measured signals. Another type of methods that can localize correlated sources is based on modeling the source signals as deterministic unknowns.

The description above is presented as a general overview of related art in this field and should not be construed as an admission that any of the information it contains constitutes prior art against the present patent application.

BRIEF DESCRIPTION OF THE FIGURES

The figures illustrate generally, by way of example, but not by way of limitation, various embodiments discussed in the present document.

For simplicity and clarity of illustration, elements shown in the figures have not necessarily been drawn to scale. For example, the dimensions of some of the elements may be exaggerated relative to other elements for clarity of presentation. Furthermore, reference numerals may be repeated among the figures to indicate corresponding or analogous elements. References to previously presented elements are implied without necessarily further citing the drawing or description in which they appear. The figures are listed below.

FIG. 1A shows entries of the unprocessed mixing matrix [M(f)]_k,k′ for different values off, with the real and imaginary part of the entries added to the k and k′ axes for the purpose of illustrating the complex function, according to some embodiments.

FIG. 1B shows entries of the phase aligned mixing matrix for different values of τ, with the absolute value of the entries added to the k axis for the purpose of illustration, according to some embodiments.

FIG. 1C shows the absolute values of the entries of the mixing matrix [M(f)]k,k′, according to some embodiments.

FIGS. 1D-1G shows the absolute values of the entries of the phase aligned mixing matrix [M(τ, 0)]_k,k′ for different values of t, according to some embodiments.

FIG. 2 is a schematic block diagram illustration of a method for estimating direction of arrival and delays of early room reflections, according to some embodiments.

FIG. 3A shows matching criterion threshold p as a function of delay τ (Equation 25), according to some embodiments.

FIG. 3B shows the angle between {circumflex over (Ω)}′(τ) (Equation (26)) and true direction of the direct sound, as a function of τ, and in which the red vertical lines correspond to the true delays of the reflections, according to some embodiment.

FIG. 4 schematically illustrates the process of Direction of Arrival (DOA) estimation, according to some embodiments.

FIG. 5 depicts local DOA estimates and delay estimates obtained by PHALCOR, according to some embodiments.

FIG. 6 depicts local DOA estimates obtained by MUSIC.

FIG. 7 illustrates probability of detection (PD) and the false positive rate (FPR) of PHALCOR and the reference method.

FIG. 8 shows a Histogram of reflection amplitudes for detected reflections and undetected reflections, according to some embodiments.

FIG. 9 shows a Histogram of reflection temporal densities for detected reflections and undetected reflections, according to some embodiments.

FIG. 10 is a schematic illustrations of a system configured to estimate direction of arrival and delays of early room reflections, according to some embodiments.

FIG. 11 is a schematic flowchart of a method for estimating direction of arrival and delays of early room reflections, according to some embodiments.

DETAILED DESCRIPTION

Aspects of disclosed embodiments pertain to systems and/or methods for blind estimation of direction of arrival and delays of reflected and non-reflected sound waves or sources originating from, for example, a single sound wave output or speaker, or multiple speakers, located in a room. It is noted that sources may contain direct sound waves and/or reflected sound waves. The systems and/or methods disclosed herein may have be employed in conjunction with a variety of applications including, for example, audio signal processing, improving audio output, improving perceived sound quality, improving Virtual or Augmented Reality experience, estimating position of walls in a room based on reflections, etc.

The term “source” as herein may also pertain to source signals or source data descriptive of the source.

Early reflections are delayed and attenuated copies of the direct sound (also: non-reflected source). Consequently, a narrowband correlation between a non-reflected source and its reflections has a phase that is linear or substantially in frequency. The property of linearity or substantial linearity can be utilized to construct a transform that can separate reflections having different delays. The transform also enables differentiating between multiple reflections having identical or similar directions of arrival (DOAs). In some examples, the transform employed for implementing part of the method may herein be referred to as Phase Aligned Correlation (PHALCOR). The expression “early reflection” may for example pertain to the 1^st room reflection, to the first 2 or the first 3 reflections, or the at least first 10 reflections, or the at least first 20 room reflections or more, or any number of initial reflections up to a certain threshold of initial reflections. The number of early reflections that may be characterized with the method presented herein may, for example, up to 100.

The term “Phase Aligned” is chosen for the reasons set forth herein. PHALCOR is based on a signal model in which the reflection signals are explicitly modeled as delayed and scaled copies of the direct sound. Before the phase alignment transformation, at each frequency, the contribution of each source to the correlation matrix is the product of a matrix with a complex number. The said matrix is not dependent on frequency. However, the complex number does depend on frequency, and specifically at each frequency, the phase (aka argument) of the number is different.

The phase alignment transform can “align” all these complex number to have the same phase. If this alignment is successful, the summation of all of the aligned correlation matrices (of different frequencies) is coherent and the source becomes more dominant. When a selected transform parameter τ matches with a delay of a certain source, the contribution of that reflection is significant. The alignment may be considered to be successful if the parameter τ of the transformation is approximately equal to the delay of the source.

The transform can separate reflections with different delays, enabling the detection and localization of reflections with similar DOAs. The DOAs and delays of the early reflections can be estimated by separately analyzing (e.g., processing) the left and right singular vectors of the transformed matrices using sparse recovery techniques.

The method may be applied on sound data that are descriptive of electronic signals which are representative of received sound waves. Sources may contain a plurality of frequencies, received via audio channels. Every audio channel may be considered to contain a mixture of all reflected sources and the direct source. It is noted that received sound waves may also encompass the meaning of the term “simulated or virtual sound waves”. In some examples, the method may include capturing reflected and non-reflected sound waves, by converting the reflected and non-reflected sound waves incident onto a plurality of microphones of one or more microphone arrays, into corresponding electronic signals. The reflected and non-reflected sound waves may be created by generating sound waves which are output by a single sound output located in the room.

In some embodiments, the method may include determining, based on the sound data, a correlation value between at least two audio channels or between each two audio channels of a plurality of audio channels to obtain, for a same frequency of a plurality of frequencies, a set of frequency-related correlation values. The set of frequency-related correlation values are represented in a first multidimensional array, which may be referred to as a volumetric or 3D spatial correlation matrix or 3D SCM. The 3D SCM may be described as being composed of a plurality of 2D SCMs for a respective plurality of frequencies. For example, each such 2D SCM relates to a spatial correlation of received sound waves of a same frequency.

In some embodiments, the method may further include performing a transpose of the frequency-related correlation values of the 3D SCM from the frequency domain to the time domain. This may be accomplished by performing an inverse weighted Fourier transform of the 3D SCM (i.e., on the values represented by the 3D SCM) to obtain a plurality of subsets (“slices”) of time delay values, represented by a second multidimensional array, which may be referred to as a 3D spatial time-delay representation matrix or 3D TDM. The 3D TDM contains a plurality of 2D TDMs relating to a respective plurality of time delays. For example, each one of such 2D TDM of the 3D TDM relates to a different time delay between two received sound waves.

In some embodiments, the method may further include analyzing each 2D TDM of time delay representations to extract information about the non-reflected sound wave and early reflected sound waves. In some embodiments, the method may include analyzing each 2D TDM of time delay representations to obtain or construct a plurality of sets of values representing a linear combination of a plurality of reflections having a same or substantially same delay. In some examples, the analyzing of the plurality of 2D TDMs of the 3D TDM may be performed through a method called Singular Value Decomposition (SVD), or any other low-rank approximation method.

In some embodiments, the method may further include analyzing the different linear combinations of sources of each set or vector, to differentiate between the source and at least one early reflection for a given delay may for example be performed through sparse recovery techniques such as, for example, Orthogonal Matching Pursuit (OMP). The column entries of the dictionary matrix are steering vectors, which are known if the geometry of the array is known, or from a calibration step.

Differentiating between at least two reflections at a given delay relative to the direct source allows determining Direction of Arrival (DOA) estimations for each one of such reflection. DOAs estimates may be determined, for example, by matching the obtained differentiated reflections and their associated delays with previously determined steering or reference vectors. A match is considered to be found if a certain matching threshold criterion is met, as is outlined herein below in more detail.

It is noted that the term “match”, as well as grammatical variations thereof also encompasses the meaning of the term “substantial match”.

In some embodiments, the method may further include performing clustering on the reflection data information for identifying outliers and clusters.

A. Notation

The notation used herein is presented briefly in this section. Lowercase boldface letters denote vectors, and uppercase boldface letters denote matrices. The k,l entry of a matrix A is denoted by [A]_k,l. The complex conjugate, transpose, and conjugate transpose are denoted by (⋅)*, (⋅)^T and (⋅)^H, respectively. The Euclidean norm of a vector is denoted by ⋅. The outer-product of two vectors a and b is the matrix ab^H. The imaginary unit is denoted by i.

S² denotes the unit sphere in R³. The symbol Ω∈S² represents a direction (DOA) in 3D-space, i.e., a pair of azimuth-elevation angles.

∠/(≠,Ω′)≙arccos(Ω·Ω′) is the angle between directions Ω and Ω′. “⋅” is the dot product in R³? ≠₀ represents the DOA of the non-reflected source.

It is noted that the delays and DOAs are determined with respect to the source signal.

B. Sound Field Representation Using Plane Wave Amplitude Density

Consider a sound field composed of M plane waves with amplitudes a₁(f), . . . , a_M(f) at frequency f, and directions Ω₁, . . . , Ω_M. The sound pressure p at any point in space x∈R³ can be formulated as follows:

$\begin{array}{cc}p\left(f,x\right)=\sum _{m=1}^M{a}_m\left(f\right)e\text{?}& \left(1\right)\end{array}$ $\text{?}\text{indicates text missing or illegible when filed}$

where k=2πf/c is the wave-number, and c is the speed of sound. When the sound field is composed of a continuum of plane waves, the summation is replaced by an integral over the entire sphere, and the amplitudes are replaced by the plane wave amplitude density (PWAD) a(f,Ω):

$p\left(f,x\right)={\int }_{{𝕊}^2}a\left(f,\Omega \right)e\text{?}d\Omega$ $\text{?}\text{indicates text missing or illegible when filed}$

For a fixed frequency, the PWAD is a function on the unit sphere. As such, it is possible to describe it by its spherical Fourier transform (SFT) coefficients (B. Rafaely, Fundamentals of Spherical Array Processing. Berlin, Germany: Springer, 2015, vol. 8. (Rafaely 2015)).

$\begin{array}{cc}{a}_{n,m}\left(f\right)\stackrel{△}{=}{\int }_{{𝕊}^2}a\left(f,\Omega \right)\left[Y_n^m\left(\Omega \right)\right]*d\Omega & \left(2\right)\end{array}$

where is the order-n and degree-m spherical harmonic. The SFT of the PWAD can be used to represent the sound pressure as follows (Rafaely 2015):

$\begin{array}{cc}p\left(f,x\right)=\sum _{n=0}^{\infty }\sum _{m=-n}^{n}4\pi i^n{j}_n\left(\mathrm{kr}\right)Y_n^m\left(\Omega \right)a_{n,m}\left(f\right)& \left(3\right)\end{array}$

where r=∥x∥, Ω=x/r, and j_n is the n'th order spherical Bessel function of the first kind. Equation (3) may be approximated by truncating the infinite sum to order N=┌kr┐ (Rafaely 2015).

k and k′ are two 2 indices in the summation. k is the index of the first source in a pair, k′ is the index of the second source in a pair.

A microphone array can be used to estimate the coefficients of the SFT of the PWAD with order less than or equal to N, by inverting the (truncated) linear transformation (3), a process known as plane wave decomposition. The existence and stability of the inverse transform depend on the frequency f and physical properties of the (typically spherical) microphone array. Further details can be found in Rafaely 2015.

The resulting signals are stacked in a vector of length (N+1)² as follows:

$a_{nm}\stackrel{\Delta }{=}{\left[a_{0,0},a_{1-1},a_{1,0},a_{1,1},\dots ,a_{N,N}\right]}^T$

Embodiments may herein be described in terms of the SFT of the PWAD. Processing and analysis in this domain may offer several advantages. First, the PWAD provides a description of the sound field that is independent of the microphone array. Second, the steering vectors, i.e., the response to a plane-wave from a given direction, are frequency independent. A steering or reference vector y(Ω) may for example be defined by the following expression:

$\begin{array}{cc}y\left(\Omega \right)\stackrel{\Delta }{=}{\frac{\sqrt{4\pi }}{N+1}\left[Y_0^0\left(\Omega \right),Y_1^{-1}\left(\Omega \right),\dots ,Y_N^N\left(\Omega \right)\right]}^H& \left(4\right)\end{array}$

In some examples, a constant

$\frac{\sqrt{4\pi }}{N+1}$

chosen for convenience such that ∥y(Ω)∥=1.

Signal Processing System Model

This section presents an example model for describing and implementing embodiments. Considering for example a sound field in a room and which is comprised of a single source, with a frequency domain signal s(f), and a DOA Ω0, relative to a measurement point in the room. As the sound wave output by the speaker propagates in the room, it is reflected from the room boundaries (e.g., walls, ceiling, floor). The k'th reflection may for example be modeled as a separate source with DOA Ω_k and signal s_k(f), which is a delayed and scaled copy of the source signal (see also: J. B. Allen and D. A. Berkley, “Image method for efficiently simulating small-room acoustics,” J. Acoustical Soc. Amer., vol. 65, no. 4 pp. 943-950, 1979 (Allen et al. 1979)), and may be expressed for instance, by the following mathematical expression:

$\begin{array}{cc}{s}_k\left(f\right)={\alpha }_k{e}^{-i2\pi f\tau k}s\left(f\right)& \left(5\right)\end{array}$

where τ_k is the delay relative to the direct sound, and αk is the scaling factor. τ₀ and a₀ are accordingly normalized to 0 and 1, respectively. In some examples it may be assumed that the delays are sorted such that τk−1≤τk.

Furthermore, a_nm(f) denotes the vector of the SFT coefficients of the PWAD, up to order N, as a function of frequency. Assuming the sources are in the far field, a_nm(f) can for example described by the following model:

$\begin{array}{cc}{a}_{nm}\left(f\right)=\mathrm{Ys}\left(f\right)+n\left(f\right)& \left(6\right)\end{array}$ $\mathrm{where}:$ $\begin{array}{cc}s\left(f\right)\stackrel{\Delta }{=}{\left[s_0\left(f\right),\dots ,s_K\left(f\right)\right]}^{\top }& \left(7\right)\end{array}$ $\begin{array}{cc}Y\left[y\left({\Omega }_0\right),\dots ,y\left({\Omega }_K\right)\right]& \left(8\right)\end{array}$

n(f) represents noise and late reverberation terms, and K is the number of early reflections. Note that although both the model and the proposed method can be generalized to include nearfield sources, this generalization may require near-field steering vectors and some information or estimation of source output distances, which could relate to source delays. Nevertheless, the near-field steering become useful for sources very close to the array (E. Fisher and B. Rafaely, “Near-field spherical microphone array processing with radial filtering,” IEEE Trans. Audio, Speech, Lang. Process., vol. 19, no. 2, pp. 256-265, February 2011 (Fisher 2011))) in which case room reflection may be considered negligible. Therefore, embodiments are herein described with the far-field assumption, which is valid for comparatively compact microphone arrays.

Any assumption described herein shall not be construed in a limiting manner. Hence, embodiments may also encompass models where these assumption may not be applied. For example, embodiments may encompass models where far-field assumptions do not apply.

In the discussion that follows R(f) denote the spatial correlation matrix (SCM) at frequency f:

$\begin{array}{cc}R\left(f\right)\stackrel{\Delta }{=}𝔼\left[a_{nm}\left(f\right){a_{nm}\left(f\right)}^H\right]& \left(9\right)\end{array}$

Substituting Eq. (6) into Eq. (9), and assuming n(f) and s(f) are uncorrelated, yields:

$\begin{array}{cc}R\left(f\right)=\mathrm{YM}\left(f\right)Y^H+N\left(f\right))& \left(10\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}N\left(f\right)\stackrel{\Delta }{=}𝔼\left[n\left(f\right){n\left(f\right)}^H\right]& \left(11\right)\end{array}$ $\begin{array}{cc}M\left(f\right)\stackrel{\Delta }{=}𝔼\left[s\left(f\right){s\left(f\right)}^H\right]& \left(12\right)\end{array}$

(Phase Alignment) Transformation of the SCM

In some embodiments, the method may be considered to include applying a transform or transformation on the SCM. The transform may be considered to be based on a “phase alignment” transformation, e.g., as defined herein.

A. Considerations for Applying the Transformation

Equation (10) can be rewritten as:

$\begin{array}{cc}R\left(f\right)=\sum _{k=0}^K\sum _{k^{\prime }=0}^K{\left[M\left(f\right)\right]}_{k,k^{\prime }}y\left({\Omega }_k\right){y\left({\Omega }_{k^{\prime }}\right)}^H+N\left(f\right)& \left(13\right)\end{array}$

When neglecting N, the matrix R is a mixture (also: linear combination) of the outer products of every pair of sources' steering vectors. The mixing coefficients are the entries of M, and therefore it is henceforth referred to as the mixing matrix. Note that the mixing coefficients are frequency dependent, but the steering vectors are not. The k,k′ entry of matrix M(f) is the correlation value between source k and source k′.

Supposing that the k,k′ entry of M(f) has a dominant magnitude, relative to all other entries. This leads to the following, for some c∈C:

$\begin{array}{cc}R\left(f\right)\approx \mathrm{cy}\left({\Omega }_k\right){y\left({\Omega }_{k^{\prime }}\right)}^H+N\left(f\right)& \left(14\right)\end{array}$

Estimating Ω_k and Ω_k′ for a case where dominant entries are likely, is easier than in the general case. However, according to (12), there may not be a dominant entry in M. Not only is the matrix M Hermitian, but also the magnitudes of its entries are products of amplitudes between pairs of two sources. Assuming the amplitude of the direct sound may be dominant, the 0, 0'th entry that corresponds to the direct sound only may indeed be dominant. Clearly, this is not useful for localizing the early reflections. The method presented herein enhance specific entries in M, allowing localizing “non-dominant” entries or reflections.

B. Phase Aligned Spatial Correlation

The following matrix is defined, herein referred to as “Phase-Aligned SCM”:

$\begin{array}{cc}\stackrel{\_}{R}\left(\tau ,f\right)\stackrel{\Delta }{=}\sum _{j=0}^{J_f-1}{w}_{j}R\left(f+j\Delta f\right)e^{i2\pi \mathrm{rj}\Delta f}& \left(15\right)\end{array}$

where τ≥0, Δf is the frequency resolution, and J_f is an integer parameter representing the overall number of frequency points. are non-negative weights. Note that when τ=0 and w_j=0, then R is identical to the SCM that may be obtained by frequency smoothing. The matrices N(τ, f) and M(τ, f) are similarly defined by replacing R in Equation (15) with N and M, respectively, such that:

$\begin{array}{cc}\stackrel{\_}{R}\left(\tau ,f\right)=\sum _{k=0}^K\sum _{k^{\prime }=0}^K{\left[\stackrel{\_}{M}\left(\tau ,f\right)\right]}_{k,k^{\prime }}y\left({\Omega }_k\right){y\left({\Omega }_{k^{\prime }}\right)}^H+\stackrel{\_}{N}\left(\tau ,f\right)& \left(16\right)\end{array}$

Similarly to Eq. (13), Eq. (16) presents R as a mixture (also: linear combination) of the outer products of every pair of the sources' steering vectors. This property of R form a basis for the method described herein.

Similarly to Eq. (13), Eq. (16) presents R as a mixture of outer products of pairs of sources' steering vectors. It is proven that for a fixed f, the k, k′ entry of M(τ,f) is maximized for τ=τ_k−τ_k′.

First, a explicit expression for the absolute value of the entries of M for an arbitrary τ is derived:

$\begin{array}{cc}\begin{array}{c}❘{\left[\stackrel{\_}{M}\left(\tau ,f\right)\right]}_{k,k^{\prime }}❘=❘\sum _{j=0}^{J_f-1}{\left[M\left(f_j\right)\right]}_{k,k^{\prime }}{w}_j{e}^{i2\pi \tau j\Delta f}❘\\ =❘\sum _{j=0}^{J_f-1}𝔼\left[s_k\left(f_j\right)s_{k^{\prime }}\left(f_j\right)\right]w_j{e}^{i2\pi \tau j\Delta f}❘\\ =❘\sum _{j=0}^{J_f-1}{\alpha }_k{\alpha }_{k^{\prime }}^{*}{\sigma }_s^2\left(f_j\right)w_j{e}^{i2\pi \tau j\Delta f\left(\tau -\left({\tau }_k-{\tau }_{k^{\prime }}\right)\right)}❘\end{array}& \left(17\right)\end{array}$

TABLE 1
Parameters for the Numerical Example:
	(ms)
0	0	1
1	2	0.7
2	6	0.5
where fj = f + jΔf and σs2(fj)   [|s(fj)|2]           (18)

The second equality in Eq. (17) is due to the definition of M(f) in Eq. (12), while the third equality is due to Eq. (5). Considering the triangle-inequality:

$\begin{array}{cc}\begin{array}{c}❘{\left[\stackrel{\_}{M}\left(\tau ,f\right)\right]}_{k,k^{\prime }}❘\le \sum _{j=0}^{J_f-1}❘{\alpha }_k{\alpha }_{k^{\prime }}^{*}{\sigma }_s^2\left(f_j\right)w_j❘\\ =❘{\left[\stackrel{\_}{M}\left({\tau }_k-{\tau }_{k^{\prime }},f\right)\right]}_{k,k^{\prime }}❘\end{array}& \left(19\right)\end{array}$

which is true since σ_s² and w_j are non-negative.

Along with Equation (16) this result implies that this result implies that among all possible delays T, it is the delay between two sources that maximizes the contribution of the outer product of their steering vectors to R (τ,f). This observation is at the core of our method. To better understand its implications, consider the following special case.

C. Special Case: White Source Signal

In this subsection the source signal is assumed to be white, such that such that σ_s²(f) is constant in f. In examples in which the weights w_j are all set to 1, equation (17) can thus be further simplified

$\begin{array}{cc}\begin{array}{c}❘{\left[\stackrel{\_}{M}\left(\tau ,f\right)\right]}_{k,k^{\prime }}❘={\sigma }_s^2❘{\alpha }_k{\alpha }_{k^{\prime }}^{*}\sum _{j=0}^{J_f-1}{e}^{i2\pi j\Delta j\left(\tau -\left({\tau }_k-{\tau }_{k^{\prime }}\right)\right)}❘\\ ={\sigma }_s^2❘{\alpha }_k{\alpha }_{k^{\prime }}^{*}{𝒟}_{J_f}\left(\Delta f\left(\tau -\left({\tau }_k-{\tau }_{k^{\prime }}\right)\right)\right)❘\end{array}& \left(20\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{𝒟}_n\left(x\right)\stackrel{\Delta }{=}\{\begin{array}{cc}n& x\in \mathbb{Z}\\ \frac{\sin \left(\pi \mathrm{nx}\right)}{\sin \left(\pi x\right)}& x\notin \mathbb{Z}\end{array}.& \left(21\right)\end{array}$

Dn(x) often arises in Fourier analysis, and is related to the Dirichlet kernel. It has a sinc-like behavior, with a main lobe centered around x=0, and a null-to-null width of 2/n. Correspondingly, |[M(τ,f)]_k,k′| has a main lobe centered at τ=τ_k−τ_k′, and a width of 2(J_fΔf)−1. Therefore, J_f determines the temporal resolution, which affects the ability to separate reflections with different delays. This result can be used as a guideline for choosing J_f. Note that the same analysis is valid for non white signals, if the weights satisfy

$w_j\propto 1/{\sigma }_s^2\left(f_j\right)$

With reference to FIGS. 1A-1G, a numerical example is now discussed for K=2 reflections and the corresponding parameter values, as listed in Table 1.

FIG. 1A presents the entries of the untransformed mixing matrix M(f). Since the signal is white, the k, k′ entry is a complex sinusoidal of the form

${\sigma }_s^2{\alpha }_k{\alpha }_{k^{\prime }}^{*}{e}^{i2\pi f\left({\tau }_k-{\tau }_k^{\prime }\right)}$

(see Eqs. (5) and (12)). Its real and imaginary parts are added to the k and k′ axes, respectively, for the purpose of illustration of the complex function. This illustration demonstrates that as the delay between two sources is decreased, the period of the corresponding entry, as a function of f, increases. At the extreme, the delay between a source and itself is zero, and so the diagonal entries are constant in frequency.

FIG. 1B presents the entries of the phase aligned mixing matrix M(τ,f), where Δf J_f=2000 Hz. In this plot, only the absolute value of the matrix entries is shown, and is added to the k-axis for the purpose of illustration.

FIG. 1C presents the absolute value of the entries of M(f). Note that the absolute value does not depend on f in this case. It is evident that M(f) is not sparse.

FIGS. 1D to 1G show several cross sections of FIG. 1B, i.e., the absolute values of the entries of M(τ,0) for selected values of r. When r is equal to a delay between two sources, which is the case in FIGS. 1E-1G, M(τ,0) is approximately sparse, and the most dominant entry is the one corresponding to the two sources. For other values of τ, all entries of M(τ,f) are relatively small. For τ=0, M(τ,f) is approximately a diagonal matrix. In the FIGS. 1C to 1G the entries are normalized such that the maximum value is 1.

Since in that case the off diagonal entries can be interpreted as correlations between different sources, this demonstrates the fact that optional frequency smoothing (τ=0) performs decorrelation of the sources. If the reflections had the same delay, M(0,0) would contain dominant off diagonal entries, and frequency smoothing would have failed to decorrelate the sources. Furthermore, while for τ=0 all 3 sources are dominant simultaneously, for values of r that correspond to delays between sources, only a subset of the sources are dominant.

D. Signal-Informed Weights Selection

The analysis presented in the previous subsection shows that if the weights {w_j}j are inversely proportional to a_s²(f_j), the phase alignment transform can effectively separate reflections with different delays. Assuming w_j=c/σ_s²(f_j) (where c is some positive constant) is equivalent to assuming w_j=1 and σ_s²(f_j)=c for all j, since their product is all that matters.

As σ_s²(f) is usually unknown, it must be estimated from the data. A very coarse, yet simple, estimate is given by the trace of R(f). By neglecting N in Eq. (13) and substituting Eqs. (12), (5) and (18), we get:

$\begin{array}{cc}\begin{array}{c}\mathrm{tr}\left(R\left(f\right)\right)\approx \sum _{k,k^{\prime }}{\left[M\left(f\right)\right]}_{k,k^{\prime }}\mathrm{tr}\left(y\left({\Omega }_k\right){y\left({\Omega }_{k^{\prime }}\right)}^H\right)\\ =\sum _{k,k^{\prime }}𝔼\left[{𝓈}_k\left(f\right){𝓈}_k^{*},\left(f\right)\right]{y\left({\Omega }_{k^{\prime }}\right)}^{B}y\left({\Omega }_k\right)\\ ={\sigma }_s^2\left(f\right)\sum _{k,k^{\prime }}{e}^{-i2\pi f\left({\tau }_k-{\tau }_{k^{\prime }}\right)}{\alpha }_k{\alpha }_k^{*},{y\left({\Omega }_{k^{\prime }}\right)}^{H}y\left({\Omega }_k\right)\\ ={\sigma }_s^2\left(f\right)\sum _k{❘{\alpha }_k❘}^2\left(1+b_k\left(f\right)\right),\end{array}& \left(22\right)\end{array}$ $\mathrm{where}$ $\begin{array}{cc}{b}_k\left(f\right)\stackrel{\Delta }{=}2\Re \left(\sum _{k^{\prime }>k}\frac{{\alpha }_k}{{\alpha }_k}{e}^{i2\pi f\left({\tau }_k-{\tau }_k^{\prime }\right)}{y\left({\Omega }_k\right)}^{H}y\left({\Omega }_{k^{\prime }}\right)\right)& \left(23\right)\end{array}$

and returns the real part of a complex scalar. We argue that b_k is typically small in comparison to 1, since usually the amplitudes decay rapidly. Furthermore, when two reflections have similar amplitudes, it is usually the case that they have very different DOAs, which implies that the inner product of their steering vectors is small (Rafaely 2015).

In some examples, the weights could have another role. Eq. (20) suggests that even if the weights are inversely proportional to σ_s², reflections with delays other than r may still be dominant in M_(τ,f), as D_n have strong side lobes. The side lobe levels can be reduced by introducing a window function, at the expense of increasing the width of the main lobe (A. V. Oppenheim, Discrete-Time Signal Processing. New York, NY, USA: Pearson Education India, 1999 (Oppenheim 1999))

E. Rank 1 Approximation of Phase Aligned SCM

In the following, the dependence on the frequency f is omitted for brevity. It is important to note that there is no direct access to M(τ); it is only indirectly observed through R(τ) as given by Eq. (16). When the proposed transformation succeeds in enhancing a single entry in M(τ), R(τ) will be dominated by a single outer product of steering vectors y(Ω)y(Ω′)^H. As an outer product is a rank 1 matrix, it is expected that the rank 1 approximation of R(τ) would perform denoising, i.e., reduce the contribution of N(τ).

The optimal rank 1 approximation (in the least squares sense) of (τ) is denoted by R(τ),and is given by truncating its singular value decomposition (SVD):

${\stackrel{\_}{R}}_1\left(\tau \right)={\sigma }_{\tau }{u}_{\tau }{v}_{\tau }^H,$

where σ_τ denotes the first (largest) singular value of R(τ), and u_τ and v_τ denote corresponding left and right singular vectors, respectively. Besides performing denoising, the SVD also separates the two steering vectors, as u_τ and v_τ are approximately equal (up to phase) to y(Ω) and y(Ω′), respectively. If there are several reflections with the same delay τ, R(τ), is still approximately of rank 1 since the dominant entries in M(τ) all appear at the same column. However, in that case u_τ is not a single steering vector, but rather a linear combination of the steering vectors of the reflections with delay τ.

V. Example Method Implementation for Estimating DOAs and Delays of Early Reflection

The method may be based on the analysis of the first singular vectors of the phase aligned SCM R(τ, f) that are presented herein.

In some embodiments, the method may first include performing plane wave decomposition on captured microphone signals (block 2110).

In some examples, analysis of SCM R(τ, f) may require that the plane wave decomposition signals to be in the frequency domain.

In some other examples, the plane decomposition signals may be approximated using, for example, the short time Fourier transform (STFT) (block 2120). STFT may enable localized analysis in both time and frequency. In some examples, it may be assumed that the window length of the STFT is sufficiently larger than τ_κ, such that the multiplicative transfer function (MTF) approximation in the STFT is applicable for Eq. (5) (Y. Avargel and I. Cohen, “On multiplicative transfer function approximation in the short-time Fourier transform domain,” IEEE Signal Process. Lett., vol. 14, no. 5, pp. 337-340, May 2007 (Avargel 2007))

In the following, r is the parameter of the phase alignment transform as in Eq. (15), and should not be confused with the time index of the STFT.

SCM Calculation:

In some embodiments, the method may include, for each time-frequency bin, estimating R (also: “3D SCM”) by replacing the expectation in Eq. (9) with averaging across J_t adjacent bins in time (block 2130). SCM calculation may, in some embodiments, involve time smoothing. Blocks 2110, 2120 and 2130 may herein be referred to as a “pre-processing” block 2100.

Frequency-to-time transform

In some embodiments, R(τ,f) may then be calculated based on R for τ=0, Δτ, . . . , (J_τ−1)Δτ using Eq. (15) (block 2210), which can be defined as performing a frequency-to-time transform. The obtained set of values R(τ,f) represents a plurality of subsets of delay representations, each set relating to a different time delay between two received sound waves. Δτ dictates the delay estimation resolution, while J_τ dictates the maximal detectable delay of a reflection. At some instances, when the discussion herein pertains to performing analysis on a single frequency, notation of the argument “f” may be omitted for brevity.

Examples for selecting these parameters is discussed further below in more detail, as well as a method to efficiently calculate R using the fast Fourier transform (FFT). The weights may for example be set using the method described in the Section titled “Signal-Informed Weights Selection”:

$\begin{array}{cc}{w}_j=\frac{W_j}{\mathrm{tr}\left(R\left(f_j\right)\right)}& \left(24\right)\end{array}$

where W_j is the j'th sample of a Kaiser window of length J_f. In some examples, the β parameter may be set to, for example, 3 (Oppenheim 1999). In some other example implementations, the β parameter may be set to values ranging from 2 to 7.

Day Detection:

The method may further include detecting delay values T that are equal to (or substantially equal or match with) reflection delays (block 2230).

The detection of such values may for example be accomplished by thresholding the following signal:

$\begin{array}{cc}\rho \left(\tau \right)=\underset{{\Omega }^{\prime }\in {𝕊}^2}{\max }❘{y\left({\Omega }^{\prime }\right)}^H{v}_{\tau }❘& \left(25\right)\end{array}$

where v_τ is a first right singular vector of R(r) obtained through SVD (block 2220). Equation 25 may herein be referred to as defining a steering matching criterion. By the Cauchy-Schwartz inequality, since both v_τ and y(Ω) are unit vectors, ρ(τ)∈[0,1] and is equal to 1 if and only if v_τ is equal (up to phase) to a steering vector. In some examples, the threshold p(τ) for associating a reflection with a steering vector may for example be set (e.g., empirically) to ρ_min=0.9. In some other example implementations, the threshold ρ(τ) may range from 0.75 to 0.99.

By {circumflex over (Ω)}′(τ) the direction is denoted that attains the maximum in Eq. (25):

$\begin{array}{cc}{\hat{\Omega }}^{\prime }\left(\tau \right)=\underset{{\Omega }^{\prime }\in {𝕊}^2}{\arg \max }❘{y\left({\Omega }^{\prime }\right)}^H{v}_{\tau }❘.& \left(26\right)\end{array}$

When τ is equal to a reflection's delay, Ω′(τ) is an estimate of Ω₀ (the DOA of the direct sound). Note that when τ is equal to a delay between two reflections (and not a delay between a reflection and the direct sound), ρ(τ) may be high as well, leading to false detections. However, by employing cluster analysis for example, such detections may be distinguishable from valid ones, as {circumflex over (Ω)}′(τ) will be different from Ω₀.

3) DOA Estimation: The method may further include estimating the DOAs of the reflections, for example, by employing a sparse recovery technique (block 2240). Estimating the DOAs may be performed separately for every τ selected in the previous step. Let u_τ denote a first left singular vector of R(τ). u_τ is approximately equal to a linear combination of the steering vectors of the reflections with a delay close to τ. If there is only a single such reflection, then its DOA can be estimated using a similar method to that of the direct sound (Eq. (26)). However, in practice there might be several reflections with similar delays. Since their number is expected to be quite small, the method may include utilizing sparse recovery to estimate the DOAs. Specifically, we aim to solve the following problem: Find the smallest set of directions {circumflex over (Ω)}₁, . . . {circumflex over (Ω)}_s and coefficients x₁, . . . x_s, such that

$\begin{array}{cc}{\sum _{𝓈=1}^{𝒮}{x}_{𝓈}y\left({\hat{\Omega }}_{𝓈}\right)-u_{\tau }}^2\le {ϵ}_u& \left(27\right)\end{array}$

where ∈_u∈(0,1) may be predefined threshold, set experimentally for example, to 0.4. In some examples, the parameter value may range from 0.2 to 0.8. The parameter ∈_u∈(0,1) represents the error between linear combinations of the steering vectors.

In the context of sparse recovery, the set of vectors {y(Ω):Ω∈S²} is known as the dictionary, and its elements are known as atoms. The optimization problem is computationally intractable and cannot be exactly solved in practice.

In some embodiments, the method may include applying the orthogonal matching pursuit (OMP) algorithm (T. T. Cai and L. Wang, “Orthogonal matching pursuit for sparse signal recovery with noise,” IEEE Trans. Inf. Theory, vol. 57, no. 7, pp. 4680-4688, July 2011 (Cal 2011)) or any other algorithm that may be extendable to infinite dictionaries.

In some examples, it may be assumed that early reflections with similar delays may have very different DOAs as they originate from different walls. If the angle between the DOAs is larger than n/N, the corresponding steering vectors are approximately orthogonal (Rafaely 2015). The sparse recovery technique may be selected such that the projection step only removes the contribution of steering vectors of DOAs the have already been found, without affecting the rest.

In some embodiments, the sparse recovery (e.g., OPM) algorithm is applied on u_τ for every detected τ. In some example implementations, values of τ where the resulting s is larger than s_max are discarded. In some examples, the value of s_max may be (e.g., empirically) set to 3, which may be the smallest set of directions meeting criterion of Equation 26. In some other example implementations, the smallest set of directions may contain 2, 3, 4, 5, 6, 7, or 8 directions.

As a result of finding the smallest set of reflections meeting the requirements defined in Equation 26, a plurality of tuples can defined, where each tuple pertains to a reflection delay estimation relative to direct sound, and a corresponding DOA estimation of the same reflection relative to direct sound, and the DOA of the direct sound. Hence, in some embodiments, each estimate is a triplet of the form ({circumflex over (r)}, {circumflex over (Ω)}, {circumflex over (Ω)}′) corresponding to the delay of a reflection, its DOA and the DOA of direct sound.

Cluster Analysis

In some embodiments, a triplet ({circumflex over (τ)}, {circumflex over (Ω)}, {circumflex over (Ω)}′) is obtained in a process of extracting direction and reflection pairs (block 2310). Clustering (block 2320) may be performed on the triplet, for example, to identify and remove outliers, for obtaining global estimates for the DOAs and delays of the early reflections (output 2400).

In some embodiments, the detection and removal of outliers may be performed by discarding reflection estimates where and angle between {circumflex over (Ω)}′ and Ω₀ is larger than a predefined outlier threshold angle. In some examples, the outlier threshold angle may be set to 10 degrees. In some other examples, the outlier threshold angle may range from 5 to 30 degrees.

In general, the DOA of the direct sound Ω₀ may not be known in advance. However, Ω₀ may be derivable by selecting the peak in the histogram of {circumflex over (Ω)}′. By discarding reflection estimates from ({circumflex over (τ)}, {circumflex over (Ω)}, {circumflex over (Ω)}′) an updated or remaining set of tuples ({circumflex over (τ)}, {circumflex over (Ω)}, {circumflex over (Ω)}′) is obtained.

In some embodiments, cluster analysis may be performed on the remaining set of tuples ({circumflex over (τ)}, {circumflex over (Ω)}, {circumflex over (Ω)}′), for example, by employing the DBSCAN algorithm (M. Ester, H.-P. Kriegel, J. Sander, and X. Xu, “A density-based algorithm for discovering clusters in large spatial databases with noise,” K dd, vol. 96, no. 34, pp. 226-231, 1996 (Ester 1996)) or any other algorithm that, for example, may not require an initial estimate of the number of clusters. Hence, in some example implementations, a clustering algorithm may be selected where the number of reflections is automatically estimated by the algorithm. Furthermore, in some examples, the clustering algorithm may be selected such that assigning an outlier to a specific cluster may not be required. Hence, in some examples, the DBSCAN cluster algorithm may be selected (Ester 1996).

However, in some embodiments, other clustering algorithms may be selected which may, for example, provide improved performance compared to DBSCAN with respect to other performance parameter values. For example, cluster algorithms may be selected from a computational viewpoint, as the complexity of DBSCAN in the present implementation is quadratic in the number of data points.

In the DBSCAN algorithm, two points are defined as neighbors of a cluster if the distance between them is less then distance threshold E. A core point of a cluster may be defined as a point having a number of neighbors that is equal or greater than MINPTS. A noise point is a non-core point, for which none of its neighbors are core points. The algorithm iterates over all the points in the dataset, and assigns two points to the same cluster if one of them is identified a core point. Noise points are not assigned to any cluster.

In some examples, the following metrics may be used:

$\begin{array}{cc}d\left(\left({\tau }_{\alpha },{\Omega }_{\alpha }\right),\left({\tau }_b,{\Omega }_b\right)\right)=\sqrt{{\left(\frac{{\mathrm{\angle \Omega }}_{\alpha },{\Omega }_b)}{\mathrm{\gamma \Omega }}\right)}^2+{\left(\frac{{\tau }_{\alpha }-{\tau }_b}{{\gamma }_{\tau }}\right)}^2}& \left(28\right)\end{array}$

where γ_Ω and γ_τ are normalization constants, set for example empirically to, e.g., 15 degrees and 500 microseconds, respectively. As the metric is normalized, the parameter value is simply set to 1. MINPTS may for example be empirically set to, e.g., 10 percent, of the number of neighbors of the point that has the largest number of neighbors.

In some embodiments, global delays and DOE estimates are calculated for each cluster, for example, by determining, for each cluster, an average or median of the local delays and DOE estimates.

The local delay estimates may be confined to a grid of resolution Δτ. However, the global delay estimates, being averages of local delay estimates, are not. The fact that each DOA estimate has an associated delay, enables the separation of clusters from one another based on the delay, even if they have similar DOAs.

Practical Considerations

Avoiding Redundant Processing

The information captured in R contains contributions from a rectangular region in the time-frequency domain of size J_t×J_f. Therefore, it can be expected that the process of frequency-to-time transformation (block 2210), analyzing the plurality of 2D TDMs of the 3D TDM resulting from the frequency-to-time transformation (block 2220) for identifying reflection delay values t that are equal to (or substantially equal or match with) a reflection delay (block 2230) of a steering vector, and further for identifying matching DOE with respect to the identified delays (block 2240), would be similar when applied to different regions having a large overlap in frequency. Hence, to reduce computation time, in some embodiments, the overlap between the time-frequency regions in blocks 2210, 2220, 2230 and 2240 can be reduced. For example, in the simulation study discussed below, an overlap of, e.g., 87.5% in frequency may be selected, such that each frequency is processed 8 times. In some examples, blocks 2210, 2220, 2230 and 2240 may herein be referred to as Phalcor—Part 1.

Acceleration Using the (Inverse) FFT

As mentioned herein, method comprises performing a frequency-to-time transformation on the correlation matrix to obtain R(τ) for a grid of values of τ. If (Δτ·Δf)⁻¹∈N, then the sequence

$\begin{array}{cc}{\left(\stackrel{\_}{R}\left(j\Delta \tau \right)\right)}_{j=0}^{J_{\tau }-1}& \left(29\right)\end{array}$

is equal (up to scaling and appropriate zero-padding) to the first J_τ terms of the inverse discrete Fourier transform (taken entry wise) of the sequence

$\begin{array}{cc}{\left(w_{j}R\left(f+j\Delta f\right)\right)}_{j=0}^{J_{\tau }-1}& \left(30\right)\end{array}$

so R(τ) can be calculated can be calculated efficiently for the grid of delays using the FFT. A further reduction in the computation time can be achieved using the following identity, obtained directly from Equation (15) and from the fact that R is Hermitian:

$\begin{array}{cc}\stackrel{\_}{R}\left(\tau \right)={\stackrel{\_}{R}\left(\Delta f^{-1}-\tau \right)}^H& \left(31\right)\end{array}$

Thus, it may in some embodiments sufficient to perform the FFT on only the upper triangular entries of R.

Periodicity of Phase Aligned SCM

It is apparent from Equation (15) that R(γ) is periodic with respect to r, with period Δf⁻¹. This periodicity does not introduce ambiguity in the delay estimation for the following reason:

When the STFT window size is T, the frequency resolution of the STFT Δf satisfies Δf≤1/T. Therefore, reflections with delay τ larger than Δf⁻¹ necessarily do not satisfy the MTF approximation criteria, since τ>T. This analysis also shows that to avoid unnecessary calculations, J_τ may in some embodiments be chosen such that J_τΔτ<Δf⁻¹/2.

Selecting the Parameters of the Phase Alignment Transform

The calculation of f may require the selection of three parameters: J_f, Δτ and J_τ. The number of frequency bins J_f, may be chosen such that the temporal resolution (given by (ΔfJ_f)⁻¹) is sufficient. For example, if J_fΔf=2000 Hz, then the phase alignment transform can separate two reflections if their delays are spaced by more than 0.5 ms. However, J_f may not be set arbitrarily high. First, the frequency independence of steering vectors is in practice limited to a given band, depending on the geometry of the microphone array. Second, some model assumptions may only be valid for bands of limited width. For example, the linear phase assumption in Equation (5) may, in practice, only hold within a local frequency region.

Once J_f has been determined, Δτ, the delay estimation resolution, may for example be set as follows:

$\begin{array}{cc}\Delta \tau =\frac{1}{M\Delta f}& \left(32\right)\end{array}$

where M is an integer that satisfies M≥J_f. This choice guarantees that Δτ≤(J_fΔf)⁻¹, and also that (Δτ·Δf)⁻¹∈N, so the FFT can be used to calculate R. Increasing M beyond J_f might increase the resolution of delay estimation; however, it would also increase the computation time of the algorithm.

In some embodiments, Jτ, the number of grid points over τ, may be chosen such that (Jτ−1)Δτ, the maximal detectable delay, is sufficiently small relative to T, the window length of the STFT such that the MTF criteria for Equation (5) holds. In some examples, (Jτ−1)Δτ≈T/10 may be considered sufficient.

Maximizing Over the Sphere

Both Equation 25 and the OMP algorithm may require maximizing functions of the form f(Ω)=|y(Ω)^Hx| over the sphere. This may be considered is equivalent to finding the maximum of a signal on the sphere whose SFT is given by x.

In some embodiments, the Newton's method may be employed for performing this maximization, with initialization obtained by sampling the sphere, for example, with a nearly uniform grid of 900 directions (J. Fliege and U. Maier, “A two-stage approach for computing cubature formulae for the sphere,” in Mathematik 139 T, Universitat Dortmund, Fachbereich Mathematik, Universitat Dortmund, 44221. Citeseer, 1996. (Fliege 1996)). Clearly, any other suitable number of grid directions may be chosen.

Computational Complexity

In some embodiments, Part 1 of PHALCOR may be independently for every selected time-frequency region, and therefore the total computation time grows linearly with the duration of the input signal. As the phase alignment transform can be calculated comparatively efficiently using the FFT, the main bottlenecks here are the SVD, delay detection and DOA estimation.

While calculating the SVD of an SCM is common in many localization methods, it is usually calculated once for every selected time-frequency region. In PHALCOR, however, it is calculated J_τ (the size of the τ-grid) times for every selected time-frequency region. As J_τ controls the maximal detectable delay, there may be a trade off between the maximal detectable delay and the computational complexity. By decreasing Jτ and increasing Δτ (τ-grid resolution), one can decrease run time without changing the maximal detectable delay. However, increasing Δτ comes at the cost of poor delay resolution.

The calculation of ρ and {circumflex over (Ω)}′(Equations (25) and (26)) may be comparatively computationally expensive as it requires a global maximum search over the sphere. For example, similar to SVD decomposition (block 2220), in PHALCOR this calculation may be required for every τ on the grid. Similarly, since the OMP algorithm is applied for every detected delay, the computational complexity of the DOA estimation step also increases with the number of delays. Accordingly, separating reflections of different delays may require the processing for each delay.

Comparison to Other Methods

Music and Frequency Smoothing

When τ=0 and w_j=1, R(τ) is a frequency-smoothed SCM (as used for example in D. Khaykin and B. Rafaely, “Coherent signals direction-of-arrival estimation using a spherical microphone array: Frequency smoothing approach,” in Proc. IEEE Workshop Appl. Signal Process. Audio Acoust., 2009, pp. 221-224 (Khaykin 2009)

Frequency smoothing is a common procedure in source localization in the presence of reverberation, as it can decorrelate signals, which is necessary for subspace methods such as MUSIC. Furthermore, R(τ) is positive semi-definite, and therefore its singular vectors are also eigenvectors, so {circumflex over (Ω)}′(0) is the estimate obtained by MUSIC if the signal subspace dimension is set to 1, and ρ(0) is equivalent to the MUSIC spectrum magnitude at this direction. While the frequency smoothing goal is to decorrelate the sources, PHALCOR actually utilizes this correlation to enhance specific reflections and to derive information about these reflections.

L1-SVD

L1-SVD is used for source localization that can address correlated sources (D. Malioutov, M. Cetin, and A. S. Willsky, “A sparse signal reconstruction perspective for source localization with sensor arrays,” IEEE Trans. Signal Process., vol. 53, no. 8, pp. 3010-3022, August 2005 (Malioutov 2005))

L1-SVD is based on the observation that the first eigenvectors of the SCM are linear combinations of steering vectors. The DOAs are estimated by decomposing the eigenvectors of the SCM into a sparse combination of steering vectors. This is similar to our method, which decomposes a first left singular vector of the phase aligned SCM to a sparse linear combination of steering vectors. In general, the performance of sparse recovery methods improves as the vectors are sparser. While in L1-SVD the sparsity is determined by the total number of reflections, in the PHALCOR employed in some embodiments, the sparsity is determined by the number of sources at a specific delay, which is significantly lower. Furthermore, like MUSIC, in L1-SVD the number of detectable sources is limited by the number of input channels ((N+1)² in our case). In PHALCOR, on the other hand, it is possible to detect more reflections than input channels, as each delay is processed independently.

Generalized Cross Correlation

The relations of PHALCOR to MUSIC and L1-SVD is related only to DOA estimation. However, PHALCOR is also related to delay estimation methods that are based on generalized cross correlation analysis. (J. Hassab and R. Boucher, “Optimum estimation of time delay by a generalized correlator,” IEEE Trans. Acoust., Speech, Signal Process., vol. A SSP-27, no. 4, pp. 373-380, August 1979 (Hassab 1979))

It can be shown that the entries of {circumflex over (R)}(τ) contain a generalized cross correlation between each pair of input channels, at lag r. Although similar, there are some important distinctions between the two methods. While cross correlation analysis is typically used to estimate the delay between two signals that are observed directly, in embodiments, PHALCOR aims to estimate the delay between multiple signals that are observed indirectly where each input channel is a linear combination of the delayed signals, given by the unknown steering matrix, which is estimated as well.

TABLE 2
Comparison of current methods and the way in which these are
overcome by employing embodiments of the method and system
Known		Example advantages of the system and
Approach	Limitations of known approach	method, according to some embodiments
Spatial	Practical arrays have restricted	R(f) is transformed to R(τ, f), making M(τ, f)
Resolution	spatial resolution, imposing a limit on	sparse, and enhancing one or a few
	the spatial separability of reflections	reflections relative to the rest, therefore
		overcoming the spatial resolution limits.
Coherent	Room reflections are coherent, see	The proposed approach exploits coherency
sources	Equation 5, and method like MUSIC	of reflection in an explicit signal model, and
	may fail, or may require frequency	does not require de-correlation of source,
	smoothing (H. Wang and M. Kaveh,	therefore overcoming this limitation
	“Coherent signal-subspace
	processing for the detection and
	estimation of angles of arrival of
	multiple wide-band sources,” IEEE
	Trans. Acoust., Speech, Signal
	Process., vol. A SS33, no. 4, pp.
	823-831, August 1985).
	The latter process may only
	decorrelate reflections with different
	delays
Reflections	Reflections with similar DOA and	The phase alignment transformation leads
with similar	different delay may not be resolved	to sparse M(τ, f) with reflections of different
DOA	by methods that rely on spatial	delay τ contributing to different sparse
	separation, such as beamforming and	matrices M(τ, f), facilitating estimation of
	MUSIC	reflections with the same DOA but different
		delay
More	Current methods may require more	PHALCOR separates reflections into groups
reflections than	microphones than reflections	with similar delay. As each group is
microphones		processed separately, the number of
		microphones limit the reflections per group,
		and not the number of reflections in total.

SIMULATION EXAMPLES

The following section provides further example implementations in conjunction with performed simulations. It is noted that the simulation configurations described herein may be analogously implemented in a real-world scenario.

TABLE 3
lists room parameters used in the simulation study discussed herein
				Average	Average
				Distance	number of
				between	reflections
	Dimen-	Reverber-		Source and	with delay
	sions	ation	Critical	Array	smaller than
Room	(m)	Time (s)	(m)	(m)	20 ms
Small	5 × 4 × 2.5	0.6	0.5	1.7	52
Medium	7 × 5 × 3	0.8	0.7	2.5	31
Large	10 × 7 × 3.5	1.1	1	3.8	21

Simulation Setup

An example simulation study is herein presented, demonstrating the performance of the PHALCOR method, according to some embodiments. First, a detailed analysis of the different steps of the algorithm is presented on a specific test case. Next, a Monte Carlo analysis is presented, demonstrating the robustness of PHALCOR.

The setup of the simulations, common to both the case study and the Monte Carlo study, is as follows. An acoustic scene that consists of a speaker and a rigid spherical microphone array in a shoe box room, was simulated using the image method (Allen et al. 1979). The speech signal is a 5 seconds sample, drawn randomly from the TSP Speech Database (P. Kabal, “Tsp speech database,” McGill University, D atabase Version, vol. 1, no. 0, pp. 09-02, 2002 (Kabal 2002)).

The array has 32 microphones, and a radius of 4.2 cm (similar to the Eigenmike (M. Acoustics, “Em32 Eigenmike Microphone Array Release Notes (v17.0),” 25 Summit Ave, Summit, N J 07901, USA, 2013)), facilitating plane wave decomposition with spherical harmonics order N=4. The microphone signals were sampled at 48 kHz. Sensor noise was added, such that the direct sound to noise ratio is 30 dB. The positions of the speaker and the array were chosen at random inside the room, with the restriction that the distance between each other, and to the room boundaries is no less than 1 m. Three different rooms sizes are considered. Their dimensions and several acoustic parameters, are presented in Table 3.

Methodology

The signals recorded by the microphones were used to compute a_nm(J). An STFT was applied to the PWAD coefficients signals using a Hanning window of 8192 samples, and an overlap of 75%. A frequency range of [500,5000] Hz was chosen for the analysis. The number of time bins used for averaging, J_t, was set to 6, while the number of frequency bins used for the phase alignment transform, J_F, was set such that J_fΔf=2000 Hz. The delay resolution Δτ was set to 83.33 microseconds (equivalent to setting M=2048 in Eq. (32)), while J_τ was chosen such that the maximal delay is 20 ms.

With these parameters, PHALCOR, was applied to the simulated data. The values of the different hyper-parameters of PHALCOR, including ρ_min, ∈_μ, S_max, γ_Ω and γ_τ, were set as exemplified herein above.

The MUSIC algorithm (Khaykin 2009) was applied as a reference method for DOA estimation, by selecting the peaks in the MUSIC spectrum ∥y(Ω)^HU∥, where U is a matrix whose columns are orthonormal eigenvectors that correspond to the L largest eigenvalues of the time and frequency smoothed SCM. The time and frequency smoothing parameters are the same as in PHALCOR. The dimension of the signal subspace L was determined using the SORTE method. (K. Han and A. Nehorai, “Improved source number detection and direction estimation with nested arrays and ulas using jackknifing,” IEEE Trans. Signal Process., vol. 61, no. 23, pp. 6118-6128, December 2013).

To reduce false positives, peaks for which the MUSIC magnitude height is lower than 0.75 were discarded. The local estimates were clustered using DBSCAN, to obtain global estimates. The delays of the detected reflections were estimated using the following method. First, each reflection signal was estimated by solving Eq. (6) for s in the least squares sense. Then, the delay of the k'th reflections was estimated by selecting the maximum of the cross correlation values between s_k and s₀. Note that in contrast to PHALCOR, the delays are estimated after the clustering.

For both PHALCOR and the reference method, a detected reflection as a true positive was considered if its delay and DOA matched simultaneously the delay and DOA of a true reflection. The matching tolerance was chosen to be 500 ρs for the delay, and 15 degrees for the DOA. The probability of detection (PD) at a given maximal delay t is defined as the fraction of true positive detections of reflections r with a delay smaller than or equal to t, out of the total number of reflections with a delay smaller than or equal to t:

$\begin{array}{cc}{\mathrm{PD}}_t\stackrel{\Delta }{=}\frac{❘\left\{r\in 𝒟:r\mathrm{is}\mathrm{true}\mathrm{positive}\mathrm{with}\mathrm{delay}\le t\right\}❘}{❘\left\{r\mathrm{ϵ𝒢𝒯}:r\mathrm{has}\mathrm{delay}\le t\right\}❘}& \left(33\right)\end{array}$

where is the set of all ground truth reflections, and is the set of all estimated reflections. The false positive rate (FPR) at a given maximal delay t is defined as the fraction of false positive detections with a delay smaller than or equal to t, out of the total number of detections with a delay smaller than or equal to t:

The false positive rate (FPR) at a given maximal delay t is defined as the fraction of false positive detections with a delay smaller than or equal to t, out of the total number of detections with a delay smaller than or equal to t:

$\begin{array}{cc}{\mathrm{FPR}}_t\stackrel{\Delta }{=}\frac{❘\left\{r\in 𝒟:r\mathrm{is}\mathrm{false}\mathrm{positive}\mathrm{with}\mathrm{delay}\le t\right\}❘}{❘\left\{r\mathrm{ϵ𝒟}:r\mathrm{has}\mathrm{delay}\le t\right\}❘}& \left(34\right)\end{array}$

Here, |⋅| denotes the cardinality of the set.

Results of a Case Study

The test case presented in this section is of a female speaker in the medium sized room. There are K=31 reflections with a delay less than 20 ms in this case. FIGS. 3a and 3b illustrate the delay detection routine. FIG. 3A shows ρ as a function of τ. Since ρ(τ) measures the similarity between v_τ, a first right singular vector of R(τ), and a steering vector (see Eq. (25)), it is high for values of τ that are close to a reflection's delay, indicated on the plot using red vertical lines. There are also values of τ that are not close to a reflection's delay, for which ρ(τ) is high.

These correspond to delays between two reflections (as opposed to delays between a reflection and the direct sound). For example, the peak near τ=2 ms, corresponds to the delay between the second and third reflections, whose delays are about 3 ms and 5 ms, respectively. Such cases may be identified by testing ∠({circumflex over (Ω)}′(τ), Ω₀), the angle between the DOA of the steering vector that is most similar to v_τ, and the DOA of the direct sound. As shown in FIG. 3B, ∠({circumflex over (Ω)}′(τ), Ω₀) is small for values of τ which are close to a reflection's delay, and high otherwise. Therefore, false detections such as τ=2 ms, will be discarded during the first step of the second part of the algorithm.

FIG. 4 schematically illustrates the process of DOA estimation, according to some embodiments. Each plot shows a different function on the sphere, which is projected onto the 2D page using the Hammer projection. The red markers correspond to ground truth DOAs Ω₀, . . . , Ω₉, with the numbers indicating the index.

In the top row, the function |y(Ω)u_τ| is shown, where each column corresponds to a different value of τ. When τ equals a reflection's delay, the direction that maximizes the response is excepted to be that of the direct sound. Indeed, as τ varies across columns, the location of the peak remains, and is equal to Ω₀, the DOA of the direct sound.

In the middle row, the function |y(Ω)u_τ| is shown. It is similar to the top row, except that a first left singular vector is used instead of a right one. When τ is a reflection's delay u_τ is approximately equal to a linear combination of the steering vectors of reflections with delays of approximately τ. When the DOAs are sufficiently separated, they can be identified as peaks in |y(Ω)u_τ|. For τ₁ and τ₂, only one such peak is apparent, and its location matches the DOA of the corresponding reflection. When τ=τ₄, it is apparent that there are two dominant peaks, at directions Ω₄ and Ω₅. This is due to the fact that the 4th and 5th reflections have similar delays. Similarly, since the 8th and 9th reflections have similar delays, when τ=τ₈ the two peaks correspond to Ω₈ and Ω₉.

FIG. 4 demonstrates the effectiveness of PHALCOR in separating reflections from the direct sound, as well as reflections with different delays. This is in contrast to the MUSIC spectrum (shown in bottom row of FIG. 4), which shows only a few peaks, which are much less separable; as a result, fewer reflections are potentially identified.

FIGS. 5 and 6 present the local estimates obtained with PHALCOR and the reference methods, respectively. DOA and delay estimates are colored by assigned cluster number.

In FIG. 5, the θ and ϕ axes are for elevation and azimuth, respectively. The red circles correspond to true DOAs and delays. In FIG. 6, The red circles correspond to Ω0, . . . , Ω31. The Hammer projection is used to project the sphere onto the plane.

It is apparent that, compared to the reference, in the example discussed, PHALCOR is able to detect significantly more reflections. PHALCOR detected successfully 29 reflections, while the MUSIC based method could only detect 8 (not including the direct sound). Furthermore, FIGS. 5 and 6 illustrate the advantage of simultaneously estimating DOA and delay for cluster analysis.

Finally, we note that while the performance of the proposed method is superior, the difference in computation time is quite significant: 303 seconds for the proposed method, and only 19 seconds for the reference method (as obtained using MATLAB 2020a, on a MacBook Pro 2019 with a 2.3 GHz 8-Core Intel Core i9 processor, 16 GB RAM).

Monte Carlo Analysis

TABLE 4
Average estimation errors for the entire Monte Carlo simulation
	DOA Error	Delay Error	Average Number of True
Method	RMS (deg)	RMS (μs)	Positive Detection
PHALCOR	4.3	77	24.8
Reference	6.5	43	3.8

The simulation described above is repeated 50 times for each of the 3 rooms, varying the speakers, their location, and the microphone array location.

FIG. 7 presents PD and FPR, for different values of t, the maximum delay of the identified reflections.

Compared with the reference method, the performance of PHALCOR is significantly better, both in terms of probability of detection and false positive rates, by a factor ranging from 3 to 20. As the delay of a reflection increases, the probability of detection decreases. This is since later reflections usually have lower amplitudes. Furthermore, the reflection density is higher as the delay increases, making it more difficult to separate the reflections spatially.

FIGS. 8 and 9 show how the distribution of reflection amplitudes and reflection temporal densities vary between detected and undetected reflections. Here, we define the temporal density of a reflection with a delay τ as the number of reflections with a delay within 1 ms of τ. The amplitudes are normalized such that the direct sound amplitude is 1. As is evident, detected reflections tend to be of higher amplitude, and lower temporal density, compared to undetected reflections. It is also apparent that this effect is more significant for temporal density.

The root mean square (RMS) for DOA and delay estimation errors for each method are computed and averaged for all the estimates in this Monte Carlo simulation, and are presented in Table 4. The RMS is calculated excluding the direct sound. Table 4 shows that the performance in terms of DOA estimation error is comparable between the two methods. In terms of delay estimation error, the reference method is superior, but note that the errors are calculated only on true positive detections, which are considerably more frequent in PHALCOR, as is evident from FIG. 7 and the last column of Table 4. Furthermore, in applications where high delay accuracy is important, it may be possible to apply cross correlation (as in the reference method) for delay estimation as a post processing step of PHALCOR.

Reference is now made to FIG. 10. In some embodiments, a system 1000 for estimating direction of arrival and delays of early room reflections may include a microphone array 1100 configured to capture sound waves 610 that are output by a speaker 600 located in a room 500. System 1000 may further include a processor 1200, a memory 1300 that is configured to store data 1310 and algorithm code 1320. Processor 1200 may be configured to execute algorithm code 1320 for the processing of data 1310 resulting in the implementation of a reflection and DOA analysis (RDOA) engine 1400. RDOA engine 1400 may be configured to implement the method outlined herein.

The term “processor”, as used herein, may additionally or alternatively refer to a controller. Processor 1200 may be implemented by various types of processor devices and/or processor architectures including, for example, embedded processors, communication processors, graphics processing unit (GPU)-accelerated computing, soft-core processors, quantum processors, and/or general purpose processors.

Memory 1300 may be implemented by various types of memories, including transactional memory and/or long-term storage memory facilities and may function as file storage, document storage, program storage, or as a working memory. The latter may for example be in the form of a static random access memory (SRAM), dynamic random access memory (DRAM), read-only memory (ROM), cache and/or flash memory. As working memory, memory 1300 may, for example, include, e.g., temporally-based and/or non-temporally based instructions. As long-term memory, memory 1300 may for example include a volatile or non-volatile computer storage medium, a hard disk drive, a solid state drive, a magnetic storage medium, a flash memory and/or other storage facility. A hardware memory facility may for example store a fixed information set (e.g., software code) including, but not limited to, a file, program, application, source code, object code, data, and/or the like.

System 1000 may further include an input/output device 1500 which may be configured to provide or receive any type of data or information. Input/output device 1500 may include, for example, visual presentation devices or systems such as, for example, computer screen(s), head mounted display (HMD) device(s), first person view (FPV) display device(s), device interfaces (e.g., a Universal Serial Bus interface), and/or audio output device(s) such as, for example, speaker(s) and/or earphones. Input/output device 1500 may be employed to access information generated by the system and/or to provide inputs including, for instance, control commands, operating parameters, queries and/or the like. For example, input/output device 1500 may allow a user of system 1000 to view and/or otherwise receive information of at least one (early) reflection and its associated DOA.

System 1000 may further comprise at least one communication module 1600 configured to enable wired and/or wireless communication between the various components and/or modules of the system and which may communicate with each other over one or more communication buses (not shown), signal lines (not shown) and/or a network infrastructure.

System 1000 may further include a power module 1700 for powering the various components and/or modules and/or subsystems of the system. Power module 1700 may comprise an internal power supply (e.g., a rechargeable battery) and/or an interface for allowing connection to an external power supply.

It will be appreciated that separate hardware components such as processors and/or memories may be allocated to each component and/or module of system 1000. However, for simplicity and without be construed in a limiting manner, the description and claims may refer to a single module and/or component. For example, although processor 1200 may be implemented by several processors, the following description will refer to processor 1200 as the component that conducts all the necessary processing functions of system 1000.

Functionalities of system 1000 may be implemented fully or partially by a multifunction mobile communication device also known as “smartphone”, a mobile or portable device, a non-mobile or non-portable device, a digital video camera, a personal computer, a laptop computer, a tablet computer, a server (which may relate to one or more servers or storage systems and/or services associated with a business or corporate entity, including for example, a file hosting service, cloud storage service, online file storage provider, peer-to-peer file storage or hosting service and/or a cyberlocker), personal digital assistant, a workstation, a wearable device, a handheld computer, a notebook computer, a vehicular device, a non-vehicular device, a robot, a stationary device and/or a home appliances control system. For example, microphone array 1100 may be part of a smartphone camera or of an autonomous or semi-autonomous vehicle, and some of RDOA engine 1400 functionalities may be implemented by the smartphone or the vehicle, and some by devices and/or system external to the smartphone or vehicle. Alternative configurations may also be conceived.

Reference is now made to FIG. 11. According to some embodiments, a method for estimating direction of arrival and delays of reflected and non-reflected sound waves in a room received by a plurality of sound capturing devices located in the room, may comprise obtaining sound data descriptive of non-reflected and reflected sound waves in a room, the sound waves containing a plurality of frequencies (block 11100). In some examples, the method sound data may be generated based on electronic signals generated by the plurality of sound capturing devices. In some examples, the plurality of sound capturing devices may be employed by a microphone array.

In some embodiments, the method may further include determining, based on the sound data, a correlation value between at least two audio channels of a plurality of audio channels to obtain, for a same frequency of the plurality of frequencies, a set of correlation values (e.g., a 3D correlation matrix) (block 11200).

In some embodiments, the method may include performing inverse weighted Fourier transform on each set of frequency-related correlation values to obtain a plurality of subsets of delay representations composing a set of values represented in a 3D delay matrix) (block 11300). Each subset of delay representation may relate to a different time delay between two received sound waves.

In some embodiments, the method may further include analyzing each subset of time delay representations to extract, for a selected delay, information about the non-reflected and reflected sources to obtain a set of linear combinations of the sources (block 11400). Each linear combination may pertain to a same or substantially same delay.

In some embodiments, the method may include analyzing the set to identify at least one reflection at a given delay (block 11500). In some examples, the set of values may represent linear combinations of a direct source and at least one reflection or of a plurality of reflections. In some examples, the analyzing may be performed to differentiate between at least two reflections of a selected linear combination having an about same delay.

In some embodiments, the method may include determining, for the at least one identified reflection at the given delay, a direction-of-arrival (DOA) estimate (block 11600). This may for instance be performed by determining or searching for a match for the at least one identified reflection with steering vectors to estimate the respective DOA. In some examples, the method may include determining which of the at least two obtained differentiated reflections and their associated delays match with steering vectors, to estimate their DOAs (e.g., respectively).

Additional Examples

Example 1 pertains to a method for estimating direction of arrival and delays of reflected and non-reflected sound waves in a room received by a plurality of sound capturing devices in the room, the method comprising:

• • obtaining sound data descriptive of received sound waves containing a plurality of frequencies;
  • determining, based on the sound data, a correlation value between each two audio channels of a plurality of audio channels to obtain, for a same frequency of the plurality of frequencies, a set of correlation values;
  • performing a transform on each set of frequency-related correlation values to obtain a plurality of subsets of delay representations of a 3D Matrix, each subset of delay representations relates to a different time delay between two received sound waves;
  • analyzing each set of time delay representations to:
  • extract information about the non-reflected sound wave; and
  • obtain a set of linear combinations of reflections having a same or substantially same delay relative to the non-reflected sound wave; and.
  • analyzing the set of linear combinations of reflections to identify at least one reflection at a given delay. Optionally, the method of example 1 may further include determining a match for the at least one reflection with a steering vector to estimate the direction-of-arrival (DOA) of the at least one reflection.

Example 2 includes the subject matter of example 1 and, optionally, performing clustering on the reflection data information for identifying outliers and clusters.

Example 3 includes the subject matter of example 1 and/or example 2 and, optionally, wherein analyzing each set of time delay representations is performed using SVD decomposition.

Example 4 includes the subject matter of any one or more of the examples 1 to 3 and, optionally, the method of any one or more of the preceding claims, wherein the step of analyzing the set of linear combinations of reflections is performed using a sparse recovery algorithm.

Example 5 includes the subject matter of example 4 and, optionally, wherein the sparse recovery algorithm is an orthogonal matching pursuit (OMP) algorithm.

Example 6 includes the subject matter of any one or more of the examples 1 to 5 and, optionally, wherein the sound data are created based on electronical signals generated by a microphone array located in the room.

Example 7 includes the subject matter of any one or more of the examples 1 to 6 and, optionally, wherein the sound waves are produced by at least one speaker located in the room.

Example 8 includes the subject matter of any one or more of the examples 1 to 7 and, optionally, comprising determining a match for a plurality of reflections with a plurality of steering vectors, respectively to estimate the direction-of-arrival (DOA) for the plurality of reflections.

Example 9 includes the subject matter of any one or more of the examples 1 to 8 and, optionally, wherein the reflections are early room reflections.

Example 10 includes the subject matter of example 9 and, optionally, wherein the early room reflections pertain to the 1st reflection (only), the first two room reflections (only), the first 3 room reflections (only), or up to the first 10 or 20 room reflections.

Example 11 pertains to a system for estimating direction of arrival and delays of reflected and non-reflected sound waves in a room received by a plurality of sound capturing devices in the room, the system comprising:

• • a processor; and
  • a memory storing instructions that are executable by the processor to cause the system to perform the following steps:
  • obtaining sound data descriptive of received sound waves containing a plurality of frequencies;
  • determining, based on the sound data, a correlation value between each two audio channels of a plurality of audio channels to obtain, for a same frequency of the plurality of frequencies, a set of correlation values;
  • performing a transform on each set of frequency-related correlation values to obtain a plurality of subsets of delay representations of a 3D Matrix, each subset of delay representations relates to a different time delay between two received sound waves;
  • analyzing each set of time delay representations to:
  • extract information about the non-reflected sound wave; and
  • obtain a set of linear combinations of reflections having a same or substantially same delay relative to the non-reflected sound wave; and
  • analyzing the set of linear combinations of reflections to identify at least one reflection at a given delay. Optionally, Example 11 may further include determining a match for the at least one reflection with a steering vector to estimate the direction-of-arrival (DOA) of the at least one reflection.

Example 12 includes the subject matter of example 11 and, optionally, performing clustering on the reflection data information for identifying outliers and clusters.

Example 13 includes the subject matter of examples 11 and/or 12 and, optionally, wherein the analyzing of each set of time delay representations is performed using SVD decomposition.

Example 14 includes the subject matter of any one or more of the examples 11 to 13 and, optionally, wherein the step of analyzing the set of linear combinations of reflections is performed using a sparse recovery algorithm.

Example 15 includes the subject matter of example 14 and, optionally, wherein the sparse recovery algorithm is an orthogonal matching pursuit (OMP) algorithm.

Example 16 includes the subject matter of any one or more of the examples 11 to 15 and, optionally, wherein the sound data are created based on electronical signals generated by a microphone array located in the room.

Example 17 includes the subject matter of any one or more of the examples 11 to 16 and, optionally, wherein the sound waves are produced by at least one speaker located in the room.

Example 18 includes the subject matter of any one or more of the examples 11 to 17 and, optionally, determining a match for a plurality of reflections with a plurality of steering vectors, respectively to estimate the direction-of-arrival (DOA) for the plurality of reflections.

Example 19 includes the subject matter of any one or more of the examples 11 to 18 and, optionally, wherein the reflections are early room reflections.

Example 20 includes the subject matter of example 19 and, optionally, wherein the early room reflections contain the 1st reflection, the first two room reflections, the first 3 room reflections, or up to the first 10 or 20 room reflections.

It is important to note that the methods described herein and illustrated in the accompanying diagrams shall not be construed in a limiting manner. For example, methods described herein may include additional or even fewer processes or operations in comparison to what is described herein and/or illustrated in the diagrams. In addition, method steps are not necessarily limited to the chronological order as illustrated and described herein.

Any digital computer system, unit, device, module and/or engine exemplified herein can be configured or otherwise programmed to implement a method disclosed herein, and to the extent that the system, module and/or engine is configured to implement such a method, it is within the scope and spirit of the disclosure. Once the system, module and/or engine are programmed to perform particular functions pursuant to computer readable and executable instructions from program software that implements a method disclosed herein, it in effect becomes a special purpose computer particular to embodiments of the method disclosed herein. The methods and/or processes disclosed herein may be implemented as a computer program product that may be tangibly embodied in an information carrier including, for example, in a non-transitory tangible computer-readable and/or non-transitory tangible machine-readable storage device. The computer program product may directly loadable into an internal memory of a digital computer, comprising software code portions for performing the methods and/or processes as disclosed herein.

The methods and/or processes disclosed herein may be implemented as a computer program that may be intangibly embodied by a computer readable signal medium. A computer readable signal medium may include a propagated data signal with computer readable program code embodied therein, for example, in baseband or as part of a carrier wave. Such a propagated signal may take any of a variety of forms, including, but not limited to, electro-magnetic, optical, or any suitable combination thereof. A computer readable signal medium may be any computer readable medium that is not a non-transitory computer or machine-readable storage device and that can communicate, propagate, or transport a program for use by or in connection with apparatuses, systems, platforms, methods, operations and/or processes discussed herein.

The terms “non-transitory computer-readable storage device” and “non-transitory machine-readable storage device” encompasses distribution media, intermediate storage media, execution memory of a computer, and any other medium or device capable of storing for later reading by a computer program implementing embodiments of a method disclosed herein. A computer program product can be deployed to be executed on one computer or on multiple computers at one site or distributed across multiple sites and interconnected by one or more communication networks.

These computer readable and executable instructions may be provided to a processor of a general purpose computer, special purpose computer, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions/acts specified in the flowchart and/or block diagram block or blocks. These computer readable and executable program instructions may also be stored in a computer readable storage medium that can direct a computer, a programmable data processing apparatus, and/or other devices to function in a particular manner, such that the computer readable storage medium having instructions stored therein comprises an article of manufacture including instructions which implement aspects of the function/act specified in the flowchart and/or block diagram block or blocks.

The computer readable and executable instructions may also be loaded onto a computer, other programmable data processing apparatus, or other device to cause a series of operational steps to be performed on the computer, other programmable apparatus or other device to produce a computer implemented process, such that the instructions which execute on the computer, other programmable apparatus, or other device implement the functions/acts specified in the flowchart and/or block diagram block or blocks.

The term “engine” may comprise one or more computer modules, wherein a module may be a self-contained hardware and/or software component that interfaces with a larger system. A module may comprise a machine or machines executable instructions. A module may be embodied by a circuit or a controller programmed to cause the system to implement the method, process and/or operation as disclosed herein. For example, a module may be implemented as a hardware circuit comprising, e.g., custom VLSI circuits or gate arrays, an Application-specific integrated circuit (ASIC), off-the-shelf semiconductors such as logic chips, transistors, and/or other discrete components. A module may also be implemented in programmable hardware devices such as field programmable gate arrays, programmable array logic, programmable logic devices and/or the like.

The term “random” also encompasses the meaning of the term “substantially randomly” or “pseudo-randomly”.

The expression “real-time” as used herein generally refers to the updating of information based on received data, at essentially the same rate as the data is received, for instance, without user-noticeable judder, latency or lag.

In the discussion, unless otherwise stated, adjectives such as “substantially” and “about” that modify a condition or relationship characteristic of a feature or features of an embodiment of the invention, are to be understood to mean that the condition or characteristic is defined to within tolerances that are acceptable for operation of the embodiment for an application for which it is intended.

Unless otherwise specified, the terms “substantially”, “about” and/or “close” with respect to a magnitude or a numerical value may imply to be within an inclusive range of −10% to +10% of the respective magnitude or value.

It is important to note that the method may include is not limited to those diagrams or to the corresponding descriptions. For example, the method may include additional or even fewer processes or operations in comparison to what is described in the figures. In addition, embodiments of the method are not necessarily limited to the chronological order as illustrated and described herein.

Discussions herein utilizing terms such as, for example, “processing”, “computing”, “calculating”, “determining”, “establishing”, “analyzing”, “checking”, “estimating”, “deriving”, “selecting”, “inferring” or the like, may refer to operation(s) and/or process(es) of a computer, a computing platform, a computing system, or other electronic computing device, that manipulate and/or transform data represented as physical (e.g., electronic) quantities within the computer's registers and/or memories into other data similarly represented as physical quantities within the computer's registers and/or memories or other information storage medium that may store instructions to perform operations and/or processes. The term determining may, where applicable, also refer to “heuristically determining”.

It should be noted that where an embodiment refers to a condition of “above a threshold”, this should not be construed as excluding an embodiment referring to a condition of “equal or above a threshold”. Analogously, where an embodiment refers to a condition “below a threshold”, this should not be construed as excluding an embodiment referring to a condition “equal or below a threshold”. It is clear that should a condition be interpreted as being fulfilled if the value of a given parameter is above a threshold, then the same condition is considered as not being fulfilled if the value of the given parameter is equal or below the given threshold. Conversely, should a condition be interpreted as being fulfilled if the value of a given parameter is equal or above a threshold, then the same condition is considered as not being fulfilled if the value of the given parameter is below (and only below) the given threshold.

It should be understood that where the claims or specification refer to “a” or “an” element and/or feature, such reference is not to be construed as there being only one of that element. Hence, reference to “an element” or “at least one element” for instance may also encompass “one or more elements”.

Terms used in the singular shall also include the plural, except where expressly otherwise stated or where the context otherwise requires.

In the description and claims of the present application, each of the verbs, “comprise” “include” and “have”, and conjugates thereof, are used to indicate that the data portion or data portions of the verb are not necessarily a complete listing of components, elements or parts of the subject or subjects of the verb.

Unless otherwise stated, the use of the expression “and/or” between the last two members of a list of options for selection indicates that a selection of one or more of the listed options is appropriate and may be made. Further, the use of the expression “and/or” may be used interchangeably with the expressions “at least one of the following”, “any one of the following” or “one or more of the following”, followed by a listing of the various options.

As used herein, the phrase “A,B,C, or any combination of the aforesaid” should be interpreted as meaning all of the following: (i) A or B or C or any combination of A, B, and C, (ii) at least one of A, B, and C; (iii) A, and/or B and/or C, and (iv) A, B and/or C. Where appropriate, the phrase A, B and/or C can be interpreted as meaning A, B or C. The phrase A, B or C should be interpreted as meaning “selected from the group consisting of A, B and C”. This concept is illustrated for three elements (i.e., A,B,C), but extends to fewer and greater numbers of elements (e.g., A, B, C, D, etc.).

It is appreciated that certain features of the invention, which are, for clarity, described in the context of separate embodiments or example, may also be provided in combination in a single embodiment. Conversely, various features of the invention, which are, for brevity, described in the context of a single embodiment, example and/or option, may also be provided separately or in any suitable sub-combination or as suitable in any other described embodiment, example or option of the invention. Certain features described in the context of various embodiments, examples and/or optional implementation are not to be considered essential features of those embodiments, unless the embodiment, example and/or optional implementation is inoperative without those elements.

It is noted that the terms “in some embodiments”, “according to some embodiments”, “for example”, “e.g.”, “for instance” and “optionally” may herein be used interchangeably.

The number of elements shown in the Figures should by no means be construed as limiting and is for illustrative purposes only.

“Real-time” as used herein generally refers to the updating of information at essentially the same rate as the data is received. More specifically, in the context of the present invention “real-time” is intended to mean that the image data is acquired, processed, and transmitted from a sensor at a high enough data rate and at a low enough time delay that when the data is displayed, data portions presented and/or displayed in the visualization move smoothly without user-noticeable judder, latency or lag.

It is noted that the terms “operable to” can encompass the meaning of the term “modified or configured to”. In other words, a machine “operable to” perform a task can in some embodiments, embrace a mere capability (e.g., “modified”) to perform the function and, in some other embodiments, a machine that is actually made (e.g., “configured”) to perform the function.

Throughout this application, various embodiments may be presented in and/or relate to a range format. It should be understood that the description in range format is merely for convenience and brevity and should not be construed as an inflexible limitation on the scope of the embodiments. Accordingly, the description of a range should be considered to have specifically disclosed all the possible subranges as well as individual numerical values within that range. For example, description of a range such as from 1 to 6 should be considered to have specifically disclosed subranges such as from 1 to 3, from 1 to 4, from 1 to 5, from 2 to 4, from 2 to 6, from 3 to 6 etc., as well as individual numbers within that range, for example, 1, 2, 3, 4, 5, and 6. This applies regardless of the breadth of the range.

The phrases “ranging/ranges between” a first indicate number and a second indicate number and “ranging/ranges from” a first indicate number “to” a second indicate number are used herein interchangeably and are meant to include the first and second indicated numbers and all the fractional and integral numerals there between.

While the invention has been described with respect to a limited number of embodiments, these should not be construed as limitations on the scope of the invention, but rather as exemplifications of some of the embodiments.

Claims

1. A method for estimating direction of arrival and delays of reflected and non-reflected sound waves in a room received by a plurality of sound capturing devices in the room, the method comprising: obtaining sound data descriptive of received sound waves containing a plurality of frequencies;determining, based on the sound data, a correlation value between each two audio channels of a plurality of audio channels to obtain, for a same frequency of the plurality of frequencies, a set of correlation values;performing a transform on each set of frequency-related correlation values to obtain a plurality of subsets of delay representations of a 3D Matrix, each subset of delay representations relates to a different time delay between two received sound waves;analyzing each set of time delay representations to: extract information about the non-reflected sound wave; andobtain a set of linear combinations of reflections having a same or substantially same delay relative to the non-reflected sound wave; andanalyzing the set of linear combinations of reflections to identify at least one reflection at a given delay.

2. The method of claim 1, further comprising determining a match for the at least one reflection with a steering vector to estimate the direction-of-arrival (DOA) of the at least one reflection.

3. The method of claim 1, further comprising performing clustering on the reflection data information for identifying outliers and clusters.

4. The method of claim 1, wherein analyzing each set of time delay representations is performed using SVD decomposition.

5. The method of claim 1, wherein the step of analyzing the set of linear combinations of reflections is performed using a sparse recovery algorithm.

6. The method of claim 5, wherein the sparse recovery algorithm is an orthogonal matching pursuit (OMP) algorithm.

7. The method of claim 1, wherein the sound data are created based on electronical signals generated by a microphone array located in the room.

8. The method of claim 1, wherein the sound waves are produced by at least one speaker located in the room.

9. The method of claim 1, comprising determining a match for a plurality of reflections with a plurality of steering vectors, respectively to estimate the direction-of-arrival (DOA) for the plurality of reflections.

10. The method of claim 1, wherein the reflections are early room reflections.

11. (canceled)

12. A system for estimating direction of arrival and delays of reflected and non-reflected sound waves in a room received by a plurality of sound capturing devices in the room, the system comprising: a processor, anda memory storing instructions that are executable by the processor to cause the system to perform the following steps:obtaining sound data descriptive of received sound waves containing a plurality of frequencies;determining, based on the sound data, a correlation value between each two audio channels of a plurality of audio channels to obtain, for a same frequency of the plurality of frequencies, a set of correlation values;performing a transform on each set of frequency-related correlation values to obtain a plurality of subsets of delay representations of a 3D Matrix, each subset of delay representations relates to a different time delay between two received sound waves;analyzing each set of time delay representations to: extract information about the non-reflected sound wave; andobtain a set of linear combinations of reflections having a same or substantially same delay relative to the non-reflected sound wave; andanalyzing the set of linear combinations of reflections to identify at least one reflection at a given delay.

13. The system of claim 12, wherein the steps further comprise determining a match for the at least one reflection with a steering vector to estimate the direction-of-arrival (DOA) of the at least one reflection.

14. The system of claim 12, further comprising performing clustering on the reflection data information for identifying outliers and clusters.

15. The system of claim 12, wherein analyzing each set of time delay representations is performed using SVD decomposition.

16. The system of claim 12, wherein the step of analyzing the set of linear combinations of reflections is performed using a sparse recovery algorithm.

17. The system of claim 16, wherein the sparse recovery algorithm is an orthogonal matching pursuit (OMP) algorithm.

18. The system of claim 12, wherein the sound data are created based on electronical signals generated by a microphone array located in the room.

19. The system of claim 12, wherein the sound waves are produced by at least one speaker located in the room.

20. The system of claim 12, comprising determining a match for a plurality of reflections with a plurality of steering vectors, respectively to estimate the direction-of-arrival (DOA) for the plurality of reflections.

21. The system of claim 12, wherein the reflections are early room reflections containing at least the 1st reflection, the at least first two room reflections, the at least first 3 room reflections, at least 10 room reflections, or at least 20 room reflections.

22. (canceled)