modeling head-related impulse responses

Abstract

A method (1900) for filtering a sound signal. The method comprises: generating (s1902) a filter pair for a particular location specified by an elevation angle θ and an azimuth angle φ, the filter pair consisting of a right filter and a left filter; filtering (s1904) a sound signal using the right filter; and filtering (s1906) the sound signal using the left filter. Generating the filter pair comprises: i) obtaining at least a first set of elevation basis function values at the elevation angle; ii) obtaining at least a first set of azimuth basis function values at the azimuth angle; iii) generating the right filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) right filter model parameters; and iv) generating the left filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) left filter model parameters.

Description

TECHNICAL FIELD

[0001] The present disclosure relates to rendering spatial sound. BACKGROUND

[0002] We have two ears that capture sound waves propagating towards us. Figure 1 The sound waves propagating towards a listener are shown as being specified by a direction of arrival (DOA) in terms of an elevation angle and an azimuth angle in a spherical coordinate system. On the path of propagation towards us, each sound wave interacts with our upper torso, head, outer ear, and the ambient medium before reaching our left and right eardrums. This interaction results in time and spectral changes to the waveforms reaching the left and right eardrums, some of which are DOA dependent. Our auditory system has learned to interpret these changes to infer various spatial properties of the sound waves themselves and the sound environment in which the listener finds him / herself. This ability, known as spatial hearing, involves how we evaluate spatial cues embedded in binaural signals (i.e., the sound signals in the right and left ear canals) to infer the location of the auditory events evoked by the sound events (physical sound sources) and the sound properties resulting from the physical environment in which we find ourselves (e.g., small room, tiled bathroom, concert hall, cave). This ability of humans (spatial hearing) can in turn be exploited to create spatial sound scenes by reintroducing spatial cues in the binaural signals that will result in a spatial perception of the sound.

[0003] The primary spatial cues include: 1) angle-dependent cues: binaural cues (i.e., interaural level difference (ILD) and interaural time difference (ITD)) and monaural (or spectral) cues; 2) distance-dependent cues: intensity and direct-to-reverberant (D / R) energy ratio. Figure 2 Examples of ITD and spectral cues for sound waves propagating towards a listener are shown. The two figures show the magnitude responses of a pair of HR filters obtained at 0 degree elevation and 40 degree azimuth (the data is from the CIPIC database: Object ID 28. The database is publicly available and can be accessed through the URL www.ece.ucdavis.edu / cipic / spatial-sound / hrtf-data / ). In the left figure, the ITD is shown as a function of frequency for the left and right ear filters. The right figure shows the magnitude responses of the left and right ear filters. Figure 1 and Figure 2In the middle, the convention of a right azimuthal direction is used and this is also the convention used in the remainder of this document. However, some HR filter sets do use another convention (where the right azimuthal direction is to the left). The mathematical representation of the time and spectral changes (1-5 milliseconds) of the DOA correlation of a short time waveform is a so-called head-related (HR) filter. The frequency domain (FD) representation of these filters is the so-called head-related transfer function (HRTF) and the time domain (TD) representation is the head-related impulse response (HRIR). A binaural rendering method based on HR filters has been established gradually, wherein a spatial sound scene is generated by directly filtering the sound source signals with the HR filters of the desired positions. This method is particularly attractive for many emerging applications, such as virtual reality (VR), augmented reality (AR), mixed reality (MR) or extended reality (XR), and mobile communication systems, where headphones are usually used.

[0004] Head-related (HR) filters are usually estimated from measurements of the impulse responses of a linear dynamic system that converts an original sound signal (input signal) into left and right ear signals (output signals) that can be measured in the ear canals of a listening subject at a predefined set of elevation and azimuth angles on the surface of a sphere of constant radius from the listening subject (e.g. an artificial head, a manikin or a human being). The estimated HR filters are usually set up as FIR filters and can be used directly in this format. To enable an efficient binaural rendering, the HRTF pairs can be converted into interaural transfer functions (ITF) or modified ITF to prevent steep spectral peaks. Alternatively, the HRTFs can be described by a parametric representation. These parametric HRTFs are easily integrated with parametric multi-channel sound coders, such as MPEG surround and spatial audio object coding (SAOC).

[0005] Rendering spatial sound signals to provide a convincing spatial perception of sound for arbitrary positions in space requires HR filter pairs for the corresponding positions, thus a set of HR filters for fine sampled positions on a 2D sphere. The minimum audible angle (MAA) characterizes the sensitivity of our auditory system as angular displacement of a sound event. With respect to azimuthal position, it is reported that for broadband noise bursts the MAA is smallest in forward and backward direction (about 1 degree) and much larger for lateral sound sources (about 10 degrees). The MAA in the median plane increases with elevation angle. It is reported that for broadband noise bursts the average MAA over elevation angles is as small as 4 degrees. Currently, there are some publicly available databases of spatially densely sampled HR filters, e.g. the SADIE database, the CIPIC database. However, they do not fully meet the MAA requirements, especially on the samples of the elevation angle. Even the SAIDE dataset of the artificial head Neumann KU100 and the KEMAR manikin contains more than 8000 measurements, its sampling resolution on the elevation angle between -15 and 15 degrees is 15 degrees, whereas 4 degrees are required according to the MAA study. It is unavoidable that an angular interpolation of HR filters is required, so that sound sources can be rendered at positions for which no actual measured filter is available. Figure 3 An example of sampling a grid on a 2D sphere is shown, where each point indicates a position for which an HR filter is measured.

[0006] Several different interpolation schemes have been developed for the angular interpolation of HR filters. Typically, M HR filter pairs are estimated from measurements at N points on the sphere where r denotes the right ear, l the left ear, denotes the elevation angle, denotes the azimuthal angle. The task is to find a function where which provides left and right filters at the unsampled angles that deliver a sound rendering with good perceptual accuracy. Once the left and right HR filters can be generated for any arbitrary position specified by It is to be noted that the superscripts l or r are sometimes omitted for simplicity without causing confusion.

[0007] Two main approaches for HRTF angular interpolation are as follows:

[0008] (1) Local-neighborhood approach: A commonly adopted approach is linear interpolation, where the missing HRTF is inferred by weighting the contributions of HRTFs measured at its closest locations around. The HRTFs can be pre-processed before interpolation, e.g., first converting the HRTFs measured at two or more closest locations to minimum phase, and then applying linear interpolation.

[0009] (2) Variational approach: A more sophisticated data-driven approach is to linearly transform the measured HRTFs to another space defined by a set of basis functions, where one set of basis functions covers the elevation angle dimension and the other set covers the frequency dimension. The basis functions can be obtained by eigen-decomposition of the covariance matrix of the measured HRTFs [1, 2]. In [3], spherical harmonics (SH) that are complete and orthogonal on a 2D sphere have been used to cover the elevation angle dimension and the azimuth angle dimension, and complex exponential functions have been used to cover the frequency dimension. The SH-based HRTF model has produced encouraging performance levels in terms of the mean square error (MSE) of the model and the perceptual loudness stability [4]. SUMMARY

[0010] The ability to accurately and efficiently render sound sources at spatial locations is one of the main features of HR-filter based spatial sound renderers. The spatial resolution of the set of HR filters used in the renderer determines the spatial resolution of the rendered sound sources. Using a set of HR filters that are coarsely sampled on a 2D sphere, VR / AR / MR / XR users often report spatial discontinuities of moving sounds. These spatial discontinuities lead to audio-video synchronization errors that significantly reduce the sense of immersion. Using a set of HR filters that are finely sampled on a sphere is one solution. However, estimating a set of HR filters from input-output measurements on a fine grid that meets the MAA requirements can be very time consuming and tedious for both the subject and the experimenter. Alternatively, it is more efficient to infer the spatially related information about missing HR filters given a sparsely sampled dataset of HR filters.

[0011] The nearest-neighbor HR filter interpolation approach assumes that at each sampling location, the HR filter only affects a region up to a certain finite distance. The HR filter at an unsampled location is then approximated as a weighted average of the HR filters at locations within a certain cutoff distance, or based on a specified number of nearest points on a straight 2D grid, e.g., where, is the HR filter vector estimated at the unsampled location and This approach is simple and has low computational complexity, which can lead to an efficient implementation. However, the interpolation accuracy can be insufficient to produce a convincing spatial sound scene. This is simply due to the fact that the variations between sample points can be more complex than what a weighted average of filters can produce.

[0012] The varying approach represents the HR filter as a linear combination of a set of basis functions (i.e. where ω p is the coefficient of the pth basis function The coefficients are typically least-squares estimates obtained by minimizing the sum of squared estimation errors over a set of measured points Given a set of basis functions, the coefficients are considered to be the "best" fit in the sense of solving a quadratic minimization problem. In principle, there is no restriction on the choice of basis functions. However, in practice, the actual choice is a set of basis functions that can effectively represent the set of HR filters in terms of estimation accuracy, in terms of the number of basis functions, and in terms of the complexity of the basis functions.

[0013] Early work on modeling HRTF magnitude responses used principal components (PCs) as basis functions, where the PCs were obtained by eigen-decomposition of the covariance matrix of HRTF magnitude responses measured for 10 listeners at 265 source positions. Using only five PCs, the resulting model achieved close to 90% of the variance in the original database. This model is efficient. It represents the original dataset well while not requiring a mechanism to interpolate HRTFs at missing positions. Recently, a hybrid approach combining principal component analysis (PCA) with a nearest-neighbor method was proposed, where the model coefficients are approximated by partial derivatives. However, this hybrid approach only achieved similar results to bilinear interpolation based on nearest neighbors.

[0014] SH has been used to model the angular dependence of a set of HRTFs. The resulting model produced encouraging levels of performance in terms of the mean squared error (MSE) of the model. However, unlike the basis functions in eigen-decomposition based models, which are fixed PC vectors, SH basis functions are complex and expensive to evaluate. An SH function of order p and degree q is written as where the associated Legendre polynomials are essentially triangular polynomials of degree P. For the entire model, the highest order (P+1) 2 SH of degree P needs to be evaluated.

[0015] ​To achieve high spatial resolution, the order of the SH representation should be as high as possible. The influence of the SH order on spatial aliasing has been investigated in the context of the perceived spatial loudness stability, which is defined as how stable the perceived loudness of a rendered sound scene is over different head orientations. Subjective results show that high order (P>10) SH HRTF representations are needed to promote high quality dynamic virtual sound scenes. This leads to L(P+1) 2 = 15488 coefficients, where L=128 corresponds to the number of frequency bins. Another study also modeled the HRTF frequency part with complex exponentials and the total number of coefficients is L(P+1) 2 where L is the truncation number of the frequency part representation. Results show that to represent the HRTF over the entire frequency range up to 20 kHz in terms of MSE, the order of the SH needs to be as large as P=30 and the truncated frequency part is L=40. The number of coefficients is 38440. Evaluating HRTFs using such high order SH HRTF models is basically impossible to implement in real-time VR / AR / MR / XR systems.

[0016] The present disclosure provides a process for generating HR filters that are sufficiently accurate and efficient for real-time VR / AR / MR / XR systems at any arbitrary position in space. In one embodiment, a varying approach is taken, where the spatial variation of a set of HR filters is modeled with B-spline basis functions and the filters are parameterized in terms of time-domain FIR filters or some mapping in the frequency domain, where the DFT is one such mapping. The resulting model is accurate in terms of MSE measures and perceptual evaluations. It is efficient in that the total number of basis functions and the amount of computation needed to evaluate the HR filters based on the model is much smaller than that of models using spherical harmonics or other such complex basis functions.

[0017] Thus, in one aspect, a method for sound signal filtering is provided. The method comprises generating a filter pair for a particular position specified by an elevation angle and an azimuth angle , the filter pair being given by a right filter and a left filter The method further includes filtering a sound signal using the right filter and filtering the sound signal using the left filter. Generating the filter pair includes: i) obtaining at least a first set of elevation basis function values at the elevation angle, ii) obtaining at least a first set of azimuth basis function values at the azimuth angle, iii) generating the right filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) right filter model parameters, and iv) generating the left filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) left filter model parameters.

[0018] In another aspect, a filtering device for filtering a sound signal is provided. The filtering device is adapted to perform a method comprising: generating a filter pair for a specific location specified by an elevation angle and an azimuth angle The filter pair is composed of a right filter and a left filter The method further includes filtering a sound signal using the right filter and filtering the sound signal using the left filter. Generating the filter pair includes: i) obtaining at least a first set of elevation basis function values at the elevation angle, ii) obtaining at least a first set of azimuth basis function values at the azimuth angle, iii) generating the right filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) right filter model parameters, and iv) generating the left filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) left filter model parameters.

[0019] The main advantages of the proposed procedure include: a) more accurate than bilinear PC-based solutions, b) more efficient than SH-based solutions, c) building a model without intensive sampling of the HR filter database, and d) the model occupies less space in memory relative to the original HR filter database. The above advantages make the proposed embodiments attractive for real-time VR / AR / MR / XR systems. BRIEF DESCRIPTION OF DRAWINGS

[0020] The accompanying drawings, which are incorporated herein and form part of the specification, illustrate various embodiments.

[0021] Figure 1 A sound wave propagating from a direction of arrival (DOA) specified by an elevation angle and an azimuth angle in a spherical coordinate system towards a listener is shown.

[0022] Figure 2An example of ITD and spectral cues of sound waves propagating towards a listener is shown.

[0023] Figure 3 An example of sampling a grid on a 2D sphere is shown.

[0024] Figure 4 An HR filter unit according to an embodiment is shown.

[0025] Figure 5 is a flowchart showing one embodiment of HR filter modeling.

[0026] Figure 6 is a flowchart of a process according to an embodiment describing pre-processing for obtaining zero-latency HR filters and ITD.

[0027] Figure 7A Delay estimates for right ear HRTF (solid line) and left ear HRTF (dashed line) on a horizontal plane with 0 degree elevation and from 0 to 360 degrees azimuth are shown.

[0028] Figure 7B Corresponding right ear HRTF (solid line) and left ear HRTF (dashed line) at 90 degrees azimuth are shown.

[0029] Figure 8 A block diagram depicting a modeling process according to an embodiment is depicted.

[0030] Figure 9 An example of B-spline basis functions is shown.

[0031] Figure 10 An example of periodic basis functions is shown.

[0032] Figure 11 A process according to an embodiment is shown.

[0033] Figure 12 An example of periodic B-spline basis functions is shown.

[0034] Figure 13A An example of B-spline basis functions is shown.

[0035] Figure 13B An example of standard B-spline basis functions is shown.

[0036] Figure 14A Another example of B-spline basis functions is shown.

[0037] Figure 14B Standard B-spline basis functions without smoothness conditions at knot-points 0 / 180 degrees are shown.

[0038] Figure 15A model representation of a HR filter dataset is shown, according to one embodiment.

[0039] Figure 16 is a block diagram of a system for generating zero-latency HR filter pairs and corresponding ITD, according to one embodiment.

[0040] Figure 17 A process for generating HR filter pairs at given HR filter model representations is shown, according to one embodiment.

[0041] Figure 18 A process for generating ITD at given ITD model representations is shown, according to one embodiment.

[0042] Figure 19 is a flowchart showing a process according to an embodiment.

[0043] Figure 20 is a flowchart showing a process according to an embodiment.

[0044] Figure 21 is a block diagram of a HR filter device 2100, according to one embodiment. DETAILED DESCRIPTION

[0045] Figure 4 A HR filter unit 400 according to an embodiment is shown. The HR filter unit 400 comprises a rendering unit 402. The unit 400 further comprises a HR filter generator 404 and an ITD generator 406 for generating HR filters and ITD respectively at arbitrary elevation and azimuth angles requested in real-time by the rendering unit 402. This requires efficient evaluation of left and right HR filter pairs based on HR filter models already loaded into the unit 400. This also requires efficient evaluation of ITD based on ITD models already loaded into the unit. This HR filter unit 400 will therefore have an interface 408 for loading HR filter models and ITD models from a database 410 of these models. The database of HR filter models is generated offline by evaluating HR filter models of different HR filter databases.

[0046] 1. HR filter set modeling

[0047] As mentioned above, a HR filter is a mathematical representation of angle-dependent spatial cues including ITD, ILD and spectral cues. ITD is defined as the difference in time of arrival of a sound signal at the two ears, as Figure 2 ​​The time delays that are not frequency dependent are removed from the HR filters and kept separately as a pure delay for each pair of HR filters. The remaining zero-time delay HR filters contain the interaural phase difference (IPD), the ILD and the spectral cues. The filters and the ITD are modeled according to the azimuth and elevation angles, respectively.

[0048] The method of variation is performed using existing HR filter databases that are mostly publicly available. The HR filters in these databases are estimated based on sound measurements done on different spatial sampling grids and are usually stored in different file formats that are naturally advantageous for the laboratory that provided the database. Recently, the spatially oriented format for acoustics (SOFA) format was developed for self-describing data with a consistent definition, which unifies the representation of different HR filter databases. Thus, in one embodiment, the SOFA format is employed, thus eliminating the need for additional effort for exchanging data formats prior to modeling. More information on the SOFA format can be found at www.sofaconventions.org / mediawiki / index.php.

[0049] Figure 5 A flowchart of one embodiment of HR filter modeling is described, in which a set of HR filters in SOFA format is loaded via the SOFA API. In a pre-processing unit, the time delays that are not frequency dependent are estimated for each HR filter, if not provided in the original database. The HR filters are then separated into zero-time delay HR filters and ITD. Finally, the zero-time delay HR filters and the ITD are modeled as a linear sum of continuous basis functions of the elevation and azimuth angles, respectively, when modeling the unit.

[0050] The steps of pre-processing, HR filter model estimation and ITD model estimation are described in more detail in the next three subsections. A description of the entire model representation is given thereafter.

[0051] 1.1 Pre-processing

[0052] The basic procedure for estimating a set of HR filters based on measurements includes the following steps:

[0053] (1) A known signal is emitted via loudspeakers placed at specified elevation azimuth and fixed distance from the subject's head;

[0054] (2) The subject's left and right ear signals are recorded using microphones placed in or at the entrance of the subject's ear canals;

[0055] (3) The recorded raw data is post-processed, which is mainly used to remove the response of the measurement system; and

[0056] (4) Using a known loudspeaker signal as an input signal and a preprocessed ear signal as an output signal, an HR filter, which is an impulse response of a linear dynamic system, is estimated based on the preprocessed data.

[0057] There is usually a frequency independent delay before the impulse response starts (onset). Some databases (e.g. CIPIC database) provide onset information. However, most databases do not provide such information. As mentioned above, the HR filter set can be modeled as a combination of a minimum phase system and a pure delay line. In this case, delay estimation is required. Given the delay information, the ITD is simply calculated by subtracting the delay of the left ear HR filter from the delay of the right ear HR filter. Secondly, the delay is removed by windowing the HR filter and a zero delay HR filter is obtained. Figure 6 A flowchart describing the preprocessing procedure for obtaining a zero-delay HR filter and ITD is shown in FIG.

[0058] In the temporal structure of the HRIR, it is easy to observe that the onset exhibits a sudden increase in amplitude after the onset occurs. Based on this temporal feature, one method for estimating the delay is to use an onset detection function that follows the energy envelope of the impulse response (IR). This onset detection function can be constructed as where {W(l): l = 1, ..., L} is a windowing function of L samples in length and R is the time step in samples between two windows. To avoid ambiguity, the angle parameter and the representation of the ear are omitted here for simplicity. The length of the window L can be chosen to cover 90% of the total energy of the HRIR. The above solution produces satisfactory results when there are strong impulse transients in the HRIR. However, this is not always the case, so by using the ratio of the cumulative energy to the total energy The solution is optimized for n=1,···N, where N is the length of the HRIR. The accumulated energy is defined as Where w(l) is a window of n points. The total energy is Another optimization takes the derivative of the ratio and finds the index of the first sample where the derivative exceeds a certain threshold as the onset index. TD It can be written as Where η is the threshold. Generally, the threshold of the ipsilateral HRTF is higher than that of the contralateral HRTF. Figure 7A and Figure 7BAn example of the delay of the HRTF estimated using the Princeton HRTF dataset - object ID 27 (URL of the database is www.princeton.edu / 3D3A / HRTFMeasurements.html) is shown. Figure 7A The plot in Fig. 1 shows the delay estimate of the right ear HRTF (solid line) and the left ear HRTF (dashed line) in the horizontal plane with 0 degree elevation and from 0 degree to 360 degree azimuth. The delay of the HRTF at 90 degree azimuth is shown in the data tip. The corresponding right ear HRTF (solid line) and left ear HRTF (dashed line) at 90 degree azimuth are shown in Fig. 2. The asterisk emphasizes the detected onset. Figure 7B The plot in Fig. 1 shows the delay estimate of the right ear HRTF (solid line) and the left ear HRTF (dashed line) in the horizontal plane with 0 degree elevation and from 0 degree to 360 degree azimuth. The delay of the HRTF at 90 degree azimuth is shown in the data tip. The corresponding right ear HRTF (solid line) and left ear HRTF (dashed line) at 90 degree azimuth are shown in Fig. 2. The asterisk emphasizes the detected onset.

[0059] When given the delay estimate, a zero time delay HR filter can be obtained by windowing the original HR filter. The most prominent position dependent influence on the spectral content of the HR filter can be traced to the outer ear or pinna, which lasts for about 0.3 ms. The "shoulder bounce" influence comes later. The total length of the position dependent IR usually does not exceed 1 ms. Therefore, a 1 ms rectangular window is long enough to preserve the main spectral related cues. A longer window can not be needed if no additional position dependent information is added.

[0060] 1.2 HR filter model estimation

[0061] The HR filters for the right ear and the left ear are modeled separately. In the following, a general time domain truncated (TD) FIR model of length N for the HR filter is given in two possible extensions (elevation extension and azimuth extension), where separate basis functions are used for elevation and azimuth.

[0062]

[0063] In the elevation extension, there is a single set of basis functions {Θ p : p = 1,..., P} for the elevation dimension and P sets of basis functions for the azimuth dimension, each set for one elevation index p, {φ p,q : q = 1,..., Q p}. K < N is the number of basis vectors of the N-dimensional vector space of filter parameter vectors, and e k is the standard orthonormal basis vector of length N:

[0064]

[0065] a = {a p,q,k : p = 1,..., P; q = 1,..., Q p; k = 1, …, K} is the set of model parameters to be estimated. Since the FIR model is truncated and K < N, the HR filter model values are equal to 0.

[0066] The azimuth spread form is a mirror image of the elevation spread form and has corresponding mirror terms. From now on we will show the properties of the elevation spread form. These properties also hold for the azimuth spread form in a mirror sense and those mirror properties can be derived by a person of ordinary skill in the art based on these properties of the elevation spread form.

[0067] The elevation spread form is very flexible in that it supports individual sets of azimuth basis functions for each elevation index p. Such comprehensive flexibility is not always needed, but it is indeed a good idea to use more than one set of azimuth basis functions. At the elevation angles of + / -90 degrees directly above and below the listener respectively, the HR filters at different azimuth angles are the same. This can be handled by using a single azimuth basis function with an elevation index p equal to 1, which has a basis function contributing to elevations + / -90 degrees. Other elevation indices can share a single different set of azimuth basis functions, where the number of basis functions Q > 1, or share several carefully chosen sets of azimuth basis functions to capture the elevation-azimuth variation of the modeled filter set.

[0068] In the following we will derive the properties of the general elevation spread form. However, it will be clear to a person of skill in the art how to modify these properties when the number of different sets of azimuth basis functions is less than P.

[0069] To estimate the model parameters {α p,q,k}, two things are needed.

[0070] (1) A minimization criterion needs to be specified, which typically takes the form of a measurement of the modeling error in the time domain, frequency domain, or a combination of both, and this criterion can even include: a regularization term for reducing the trend of overfitting the modeled data.

[0071] (2) An optimization method for estimating the parameters that minimize the minimization criterion.

[0072] Figure 8 depicts a block diagram of the modeling process when given a set of zero-delay HR filters associated with corresponding elevation angles and azimuth angles (i.e., ). Given a list of elevations and azimuths, basis functions on the elevation angles and azimuth angles are constructed respectively. Then the least squares method is used to estimate the model parameters. This model estimation process is described in more detail in subsection 1.2.1.

[0073] This model estimation process is described in more detail in subsection 1.2.1.

[0074] At this stage, the model specification is quite general, as no two sets of basis functions have been specified and The key to obtaining a model that can output accurate predictions and that evaluates HR filter models efficiently lies in the choice of these two sets of basis functions. After experimenting with different types of functions, we chose to use what we call periodic B-spline basis functions as the azimuthal basis functions and standard B-spline functions as the elevation basis functions. The chosen basis functions are explained in more detail in subsection 1.2.2.

[0075] 1.2.1 Model parameter estimation

[0076] Given a set of basis functions and and a set of zero-lag HR filters for the right or left ear sampled at M different angular positions The typical minimization criterion in the time domain is the sum of the norms of the modeling errors over the set of M HR filters (right or left ear):

[0077]

[0078] where,

[0079] and

[0080]

[0081] The number of estimated parameters is which should be much smaller than the number of available data samples M*N to avoid an ill-conditioned system.

[0082] Because the standard orthonormal basis vectors e k are orthogonal vectors, the parameters a k = {a p,q,k : p = 1,..., P; q = 1,..., Q p} can be solved independently for each sample k. For each sample k, the minimization criterion becomes:

[0083]

[0084] which can be expressed in matrix form as:

[0085]

[0086] where,

[0087]

[0088] J(a k) is the linear least squares criterion. The solution k = (B T B) -1 B T h k minimizes J(a k ). However, directly minimizing the above cost function leads to an exact solution for linear systems. This solution is sensitive to noise in the data and can lead to overfitting. Tikhonov regularization is then applied and the minimization criterion becomes:

[0089]

[0090] where I is the identity matrix of size and 0 is the 0 column vector with elements.

[0091] J(a k ) is also a linear least squares criterion. Similarly, the solution minimizes J(a ) is obtained by solving the normal equations where the value of λ can be determined to be such that the condition number of the matrix is less than 10 or some other value that leads to good model accuracy.

[0092] For better numerical accuracy, the actual solution for a k is obtained by means of singular value decomposition (SVD) of B:

[0093]

[0094] The columns and U span the same subspace. The projection onto the orthogonal matrix U is given by and is equal to This gives which leads to the solution This estimate is very efficient as it only requires one SVD of the matrix of smaller dimension, which is then used to evaluate the solution for k = 1,..., K, which can be done in parallel. Replacing by B and k by h this holds for J(a k ) as well.

[0095] Given the right ear HR filter measurements we obtain the set of model parameters represented by where each a r is a vector of dimension Similarly, given the left ear HR filter measurement result We obtain Represents the set of model parameters, where each α l The dimension is Column vector of .

[0096] The minimization criteria J(α) and They can be easily mapped to the frequency domain by converting the time domain vectors into the frequency domain using a DFT transform or similar processing (e.g., interaural transfer function (ITF)). and is mapped into a frequency domain vector, and the alternative criterion can easily use a combination of time domain components and frequency domain components.

[0097] The square norm of a vector v is defined as the inner product of the vector and itself ||v|| 2 =<v,v> The general form of the inner product is<v,v> =v T Γv, where Γ can be any positive definite matrix and in its simplest form Γ is the identity matrix. Using a Γ different from the identity matrix provides a mechanism for weighting different components in the time and frequency domains differently, which may be useful when some components are perceptually more important than others. It will be clear to those skilled in the art how to take advantage of these possible variations in the specification of the minimization criterion.

[0098] 1.2.2 Definition of elevation and azimuth basis functions

[0099] As explained before, after experimenting with different types of basis functions, we chose to use standard B-spline functions as elevation basis functions and what we call periodic B-spline basis functions as azimuth basis functions.

[0100] variable The J-order univariate B-spline basis function set (where Located in the interval (in) is a set of J-1 degree piecewise polynomial functions defined on this interval. The range of these polynomial functions uses the so-called knot sequence θ=(θ1,…,θ U )(where θ1=θ A ,θ U =θ B ), and the subintervals of these polynomial functions are u=1,…,U-1. In each subinterval, each basis function is a J-1 degree polynomial function, which can be written as:

[0101] for

[0102] The smoothness of the functions (which are linear combinations of B-spline basis functions) at the knots is controlled using a so-called multiplicity sequence m = (m1,..., m U ) which is a sequence of integers greater than 0, where the value m u = i indicates that the (J-i)th derivative at the knot θ u is continuous. This means that i = 1 gives the maximum smoothness, while i = J gives only 0th derivative continuity. Given a sequence of knots and a multiplicity sequence, the polynomial model coefficients Details of this procedure can be found in the article "Bspline Basics" by Carl de Boor (ftp: / / ftp.cs.wisc.edu / Approx / bsplbasic.pdf).

[0103] An example of B-spline basis functions evaluated for an elevation angle of 3 degrees using the sequence of knots θ = (-90, -60, -30, 0, 30, 60, 90) and the multiplicity sequence m = (4, 1, 1, 1, 1, 1, 3) is shown in Figure 9 .

[0104] Azimuth angles of degrees are periodic (e.g. cyclic) in the sense that they are identical to points of azimuth angles of degrees (for any integer value of κ) in space (e.g. in the sense of the Euclidean distance), and it is important to use periodic basis functions (i.e. ) in the same way in order to obtain an efficient modeling in the azimuth dimension. An example of such a periodic basis function is shown in Figure 10 , where the part of the function in the range of angles from 0 to 360 is plotted using a solid line and the part of the function outside this range is plotted using a dashed line.

[0105] We have devised a method for generating a set of periodic B-spline basis functions on the azimuth range of 0 to 360 degrees. This method is shown in Figure 11 and comprises the following steps.

[0106] (Step 1) Specify a sequence of knots on the range of 0 to 360 degrees. Denote the length of this sequence of knots as L.

[0107] (Step 2) Extend this sequence of knots in a periodic way using J values less than 0 degrees and J-1 values greater than 360 degrees.

[0108] (Step 3) Use the extended node sequence and the extended multiplet sequence of the node to generate an extended set of B-spline basis functions, using standard methods to generate a set of B-spline functions.

[0109] (Step 4) Select L-1 consecutive extended basis functions from the extended basis functions, starting with index 2, and map them in a periodic fashion to the azimuth angle range from 0 to 360 degrees.

[0110] This method provides a set of L-1 periodic basis functions over the range from 0 to 360 degrees.

[0111] Each basis function over the azimuth angle is also a J-1th order polynomial function, and is written as:

[0112]

[0113] Figure 12 An example of a periodic B-spline basis function for an azimuth angle of 3 degrees, evaluated using a node sequence of length L = 11, φ = (0, 30, 70, 110, 150, 180, 210, 250, 290, 330, 360), is shown in FIG. 3.

[0114] 1.3 ITD model estimation

[0115] The general form of the ITD model is given by:

[0116]

[0117] and are B-spline basis functions for the elevation angle and the azimuth angle, respectively. is a set of model parameters.

[0118] 1.3.1 Model parameter estimation

[0119] The model parameters {c p′,q′},

[0120]

[0121] where,

[0122]

[0123] is the ITD for is the frequency-independent time delay provided by the original database or estimated using the method described in Section 1.1.

[0124] ​Applying Tikhonov regularization to avoid overfitting, the minimization criterion becomes

[0125]

[0126] where is the identity matrix of size and 0 is a 0-column vector with elements.

[0127] The value of can be determined such that the condition number of the matrix is less than 10 or some other value that leads to good model accuracy. With the help of the SVD of , the model parameters are obtained by c = V'S' -1 U' T τ, which is a column vector with elements.

[0128] 1.3.2. Specification of elevation and azimuth basis functions

[0129] When the elevation moves from -90 degrees upwards to 90 degrees, the ITD increases from zero to a maximum at elevation 0 degrees and decreases to zero afterwards. Based on this, it is natural to use basis functions that are zero at + / - 90 degrees elevation. This requirement is equivalent to at least one smoothness condition at + / - 90 degrees elevation. As explained in subsection 1.2.2, the smoothness of a function at a node is controlled by the multiplicity m. Write each basis function as

[0130]

[0131] An example of B-spline basis functions for an elevation angle of 3 degrees, evaluated using the node sequence Figure 13A and the multiplicity m = (3, 1, 1, 1, 2) is shown in

[0132] Considering that the ITD can not be exactly zero at + / - 90 degrees elevation due to asymmetries in the measurement setup and objects, it is still a good choice to use standard B-spline basis functions without the smoothness condition at the nodes + / - 90 degrees. An example of standard B-spline basis functions for an elevation angle of 3 degrees, evaluated using the node sequence Figure 13B and the multiplicity m = (4, 1, 1, 1, 3) is shown in

[0133] ​As the azimuth angle moves along the circumference, the change in ITD takes on a sine-like shape, where zero ITD occurs at azimuth angles of 0 / 180 / 360 degrees and maximum ITD occurs at azimuth angles of 90 / 270 degrees. Similarly, a smoothness condition can be satisfied at azimuth angles of 0 / 180 / 360 degrees. Furthermore, it can be considered that the ITD at azimuth angles between 180 and 360 degrees is the mirror image of the ITD at azimuth angles between 0 and 180 degrees. Therefore, we use a set of basis functions for azimuth angles in two intervals [0, 180] and [180, 360]. Each basis function is written as:

[0134]

[0135] exist Figure 14A The initial node sequence is shown in Example of a B-spline basis function for an azimuth angle of 3 degrees evaluated with the multiple sequence m = (3, 1, . . . , 1, 2).

[0136] Considering that the ITD may not be exactly zero at the azimuth angles of 0 / 180 / 360 degrees, the standard B-spline basis function can be used without the smoothing conditions at the knots 0 / 180 degrees. Figure 14B An example of such a basis function is shown in .

[0137] 1.4 Model Representation

[0138] Figure 15 A model representation of an HR filter dataset is shown. The representation includes a zero-delay HR filter model representation and an ITD model representation, each including basis functions and model parameters. Key to the modeling accuracy and computational efficiency of the modeling solution is a carefully constructed set of B-spline basis functions for modeling the angular variation of the HR filter set that is simple enough to give good computational efficiency and rich enough to give good modeling accuracy.

[0139] For the zero-delay HR filter model, there are: P elevation B-spline basis functions; all contain Q p P sets of azimuthal B-spline basis functions of functions; and Multiply the two sets of model parameters by the matrix K. For the ITD model, there are: elevation B-spline basis functions; all contain The azimuth B-spline basis function of the function a collection; and as a A set of model parameters of a vector of elements.

[0140] Each set of B-spline basis functions is represented by a sequence of knots and a polynomial model coefficient γ as a three-dimensional array. The first dimension corresponds to the order of the B-spline, the second dimension corresponds to the number of knot intervals, and the third dimension corresponds to the number of basis functions.

[0141] P or is much smaller than the number of elevation angles in the original HR filter dataset. is much smaller than the number of azimuth angles in the dataset. K is also smaller than the length or number of frequency bins of the original filter. Therefore, this model representation is efficient in representing HR filter datasets.

[0142] Furthermore, because the angular basis functions are continuous, the model representation can be used to generate HR filter pairs at any arbitrary location specified by elevation and azimuth angles.

[0143] 2. HR filter generation

[0144] Figure 16 This is a block diagram of a system for generating a zero-delay HR filter pair (i.e., right and left ear filters) and the corresponding ITD, given a model representation. The model representation can be written to a binary or text file. The file is loaded via the API to retrieve the model structure. The following describes how to use this model representation to obtain an HR filter pair and ITD at a specified location.

[0145] 2.1 Generating Zero-Delay HR Filter

[0146] Figure 17 It shows that when the HR filter model representation is given, the The process of generating zero-delay HR filter pairs. As described in Section 1.2.2, the model of the elevation B-spline basis function set {Θ p :p=1,···,P} includes: node sequence θ=(θ1,…,θ U ), which specifies the subinterval The function is a polynomial on the subinterval; and a 3-dimensional array indicating the model parameters Involves evaluating P elevation basis functions at elevation angles The value at The steps are as follows:

[0147] (1) Finding satisfaction The index u of ; and

[0148] (2) Evaluate the p-th elevation angle B-spline basis function at the elevation angle according to the following formula: The value at:

[0149]

[0150] The similar procedure is used to evaluate the values of the azimuth B-spline basis function set at a given azimuth angle

[0151] Once these sets of basis function values are obtained, the right ear zero-time delay HR filter at position is obtained according to:

[0152]

[0153] Based on this, the evaluation of is also clear according to:

[0154]

[0155] The left ear zero-time delay HR filter at position is obtained according to:

[0156]

[0157] The evaluation of is obtained according to:

[0158]

[0159] 2.2 Generating ITD

[0160] Figure 18 The procedure for generating ITD at position when given an ITD model representation is shown according to one embodiment.

[0161] Following the procedure described in subsection 2.1, the values of the elevation basis functions at elevation angle and the values of the azimuth B-spline basis function set at a given azimuth angle are evaluated. Once the values of the elevation basis functions and azimuth basis functions are evaluated, the ITD is obtained according to:

[0162]

[0163] As mentioned in subsection 5.1.1, we model the HR filter set as a combination of a minimum phase like system and a pure delay line. The delay of the right ear HR filter is:

[0164]

[0165] The delay of the left ear HR filter is:

[0166] ​

[0167] It is noted that, and the calculations should be consistent with the definition of the ITD used and the coordinate system.

[0168] Figure 19 is a flowchart illustrating a procedure 1900 according to an embodiment. The procedure 1900 can start from step s1902.

[0169] Step s1902 comprises generating a filter pair for a specific position specified by an elevation angle and an azimuth angle , the filter pair consisting of a right filter and a left filter .

[0170] Step s1904 comprises filtering a sound signal using the right filter.

[0171] Step s1906 comprises filtering a sound signal using the left filter.

[0172] As illustrated in Figure 20 , step s1902 comprises i) obtaining at least a first set of elevation basis function values at the elevation angle (step s2002), ii) obtaining at least a first set of azimuth basis function values at the azimuth angle (step s2004), iii) generating the right filter using a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) right filter model parameters (step s2006), and iv) generating the left filter using a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) left filter model parameters (step s2008).

[0173] In some embodiments, obtaining the first set of azimuth basis function values comprises obtaining P sets of azimuth basis function values, wherein the P sets of azimuth basis function values comprise the first set of azimuth basis function values and wherein, (p = 1 to P, q = 1 to Q p , and k = 1 to K) are left model parameter sets, (p = 1 to P, q = 1 to Q p , and k = 1 to K) are right model parameter sets, (p = 1 to P) defines a first set of elevation basis function values at an elevation angle , and (p = 1 to P, and q = 1 to Qp ) the P sets of azimuthal basis function values are defined at azimuthal angles e k (k = 1 to K) are a set of N-length orthonormal basis vectors.

[0174] In some embodiments, obtaining the first set of elevation basis function values includes obtaining Q sets of elevation basis function values, where the Q sets of elevation basis function values include the first set of elevation basis function values and where, (p = 1 to P q , q = 1 to Q, and k = 1 to K) are a set of left model parameters, (p = 1 to P q , q = 1 to Q, and k = 1 to K) are a set of right model parameters, (q = 1 to Q, and p = 1 to P q ) define the Q sets of elevation basis function values at elevation angles and (q = 1 to Q) define the first set of azimuthal basis function values at azimuthal angles e k (k = 1 to K) are a set of N-length orthonormal basis vectors.

[0175] In some embodiments, obtaining the first set of elevation basis function values includes evaluating, for each elevation basis function included in the first set of elevation basis functions, the elevation basis function at the elevation angles to produce an elevation basis function value corresponding to the elevation angles and the elevation basis function, and obtaining the first set of azimuthal basis function values includes evaluating, for each azimuthal basis function included in the first set of azimuthal basis functions, the azimuthal basis function at the azimuthal angles to produce an azimuthal basis function value corresponding to the azimuthal angles and the azimuthal basis function.

[0176] In some embodiments, each elevation basis function included in the first set of elevation basis functions is a B-spline basis function, and each azimuthal basis function included in the first set of azimuthal basis functions is a periodic B-spline basis function.

[0177] In some embodiments, the process further includes obtaining a model representing at least the first set of elevation basis functions, where the model includes a sequence (0), where 0 = (0i,..., 0 U ), which specifies subintervals the elevation basis functions are polynomials over the subintervals; and a three-dimensional array of model parameters

[0178] In some embodiments, the first set of elevation basis functions comprises a p-th elevation basis function, evaluating each of the elevation basis functions comprised in the first set of elevation basis functions at the elevation angle comprises evaluating a p-th elevation basis function at the elevation angle and evaluating the p-th elevation basis function at the elevation angle comprises the steps of finding an index u satisfying and evaluating the p-th elevation basis function at the elevation angle according to .

[0179] In some embodiments, the process further comprises obtaining a model representing at least the first set of azimuth basis functions, wherein the model comprises a sequence (f1) where which specifies sub-intervals the azimuth basis functions are polynomials over the sub-intervals; and a three-dimensional array of model parameters

[0180] In some embodiments, the first set of azimuth basis functions comprises a q-th azimuth basis function, evaluating each of the azimuth basis functions comprised in the first set of azimuth basis functions at the azimuth angle comprises evaluating the q-th azimuth basis function at the azimuth angle and evaluating the q-th azimuth basis function at the azimuth angle comprises the steps of finding an index 1 satisfying and evaluating the q-th azimuth basis function at the azimuth angle according to .

[0181] In some embodiments, the process further includes generating at least a first set of azimuthal basis functions, wherein generating the first set of azimuthal basis functions includes generating a set of periodic B-spline basis functions over a 0 to 360 degree azimuthal range. In some embodiments, generating the set of periodic B-spline basis functions over the 0 to 360 degree azimuthal range includes specifying a sequence of nodes of length L over a 0 to 360 degree range; generating an extended sequence of nodes based on the sequence of nodes of length L, wherein generating the extended sequence of nodes includes extending the sequence of nodes of length L in a periodic manner with J values less than 0 degrees and J-1 values greater than 360 degrees; obtaining an extended multiplicity sequence of nodes; generating an extended set of B-spline basis functions using the extended sequence of nodes and the extended multiplicity sequence of nodes; selecting L-1 consecutive extended basis functions from the extended basis functions starting from index 2; and mapping the selected extended basis functions to the 0 to 360 degree azimuthal range in a periodic manner.

[0182] In some embodiments, the process further includes determining an elevation-azimuth angle determining an interaural time difference In some embodiments, the process further includes determining a right delay based on determining a right delay and determining a left delay based on determining a left delay In some embodiments, filtering the sound signal using the right filter includes filtering the sound signal using the right filter and the right delay In some embodiments, filtering the sound signal using the left filter includes filtering the sound signal using the left filter and the left delay In some embodiments, filtering the sound signal using the right filter and the right delay includes computing In some embodiments, filtering the sound signal using the left filter and the left delay includes computing In some embodiments, filtering the sound signal using the left filter and the left delay includes computing In some embodiments, filtering the sound signal using the left filter and the left delay includes computing where u(n) is the sound signal.

[0183] In some embodiments,

[0184] and

[0185]

[0186] Figure 21 is a block diagram of an HR filtering device 2100 for implementing the HR filtering unit 400 according to some embodiments. That is, the device 2100 is operable to perform the processes disclosed herein. As shown in FIG. 21, the device 2100 includes a processor 2102, a memory 2104, and a communication interface 2106. The processor 2102 is configured to execute instructions stored in the memory 2104 to perform the processes disclosed herein. The communication interface 2106 is configured to communicate with other devices, such as the HR filtering unit 400. Figure 21As shown, the apparatus 2100 can include: processing circuitry (PC) 2102, which can include one or more processors (P) 2155, e.g., general-purpose microprocessors, and / or one or more other processing units such as an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or the like, which can collectively be comprised of a single processor or multiple processors in a single housing or in multiple data centers, or which can be geographically dispersed (i.e., the apparatus 2100 can be a distributed computing apparatus); a network interface 2148, including a transmitter (Tx) 2145 and a receiver (Rx) 2147 for enabling the apparatus 2100 to send and receive data to and from other nodes connected to a network 110 (e.g., an Internet Protocol (IP) network), where the network interface 2148 is connected (directly or indirectly) to that network 110 (e.g., the network interface 2148 can be connected wirelessly to the network 110, in which case the network interface 2148 is connected to an antenna arrangement); and a local storage unit (also referred to as a "data storage system") 2108, which can include one or more non-volatile storage devices and / or one or more volatile storage devices. In embodiments where the PC 2102 includes a programmable processor, a computer program product (CPP) 2141 can be provided. The CPP 2141 includes a computer readable medium (CRM) 2142 storing a computer program (CP) 2143 comprising computer readable instructions (CRI) 2144. The CRM 2142 can be a non-transitory computer readable medium (e.g., magnetic

[0187] The following is a summary of various embodiments described herein:

[0188] A1. A method for filtering of a sound signal, the method comprising: generating a filter pair for a specific position specified by an elevation angle and an azimuth angle , the filter pair comprising a right filter and a left filter constituting; filtering a sound signal using the right filter; and filtering the sound signal using the left filter, wherein generating the filter pair comprises: i) obtaining at least a first set of elevation basis function values at the elevation angle; ii) obtaining at least a first set of azimuth basis function values at the azimuth angle; iii) generating the right filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) right filter model parameters; and iv) generating the left filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) left filter model parameters.

[0189] A2. The method of claim Al, wherein obtaining the first set of azimuth basis function values comprises obtaining P sets of azimuth basis function values, wherein the P sets of azimuth basis function values include the first set of azimuth basis function values.

[0190] A3. The method of claim Al, wherein generating the right filter comprises computing and generating the left filter comprises computing wherein, (p = 1 to P, q = 1 to Q p , and k = 1 to K) are a set of right model parameters, (p = 1 to P, q = 1 to Q p , and k = 1 to K) are a set of left model parameters, (p = 1 to P) define the first set of elevation basis function values at the elevation angle and (p = 1 to P, and q = 1 to Q p ) define P sets of azimuth basis function values at the azimuth angle and e k (k = 1 to K) are a set of N-length orthonormal basis vectors.

[0191] A4. The method of claim Al, wherein obtaining the first set of elevation basis function values comprises obtaining Q sets of elevation basis function values, wherein the Q sets of elevation basis function values include the first set of elevation basis function values.

[0192] A5. The method of claim Al, wherein generating the right filter comprises computing and generating the left filter comprises computing wherein, (p = 1 to P q , q = 1 to Q, and k = 1 to K) are a set of right model parameters, (p = 1 to P q , q = 1 to Q, and k = 1 to K) are a left model parameter set, (q = 1 to Q, and p = 1 to P q ) define Q sets of elevation basis function values at elevation angles , and (q = 1 to Q) define the first set of azimuth basis function values at azimuth angles , and e k (k = 1 to K) is a set of N-length orthonormal basis vectors.

[0193] A6. The method of any one of claims Al to A5, wherein each of the elevation basis function values depends on the azimuth angle, and / or each of the azimuth basis function values depends on the elevation angle.

[0194] A7. The method of any one of claims Al to A5, wherein obtaining the first set of elevation basis function values comprises, for each elevation basis function included in the first set of elevation basis functions, evaluating the elevation basis function at the elevation angle to produce an elevation basis function value corresponding to the elevation angle and the elevation basis function, and obtaining the first set of azimuth basis function values comprises, for each azimuth basis function included in the first set of azimuth basis functions, evaluating the azimuth basis function at the azimuth angle to produce an azimuth basis function value corresponding to the azimuth angle and the azimuth basis function.

[0195] A8. The method of claim A7, wherein each elevation basis function included in the first set of elevation basis functions is a B-spline basis function, and each azimuth basis function included in the first set of azimuth basis functions is a periodic B-spline basis function.

[0196] A9. The method of claim A7 or A8, further comprising obtaining a model representing at least the first set of elevation basis functions, wherein the model comprises a sequence (0) where 0 = (0i,..., 0 U ), which specifies subintervals the elevation basis functions are polynomials over the subintervals; and a three-dimensional array of model parameters

[0197] A10. The method of claim A9, wherein the first set of elevation basis functions includes a Pth elevation basis function, evaluating each elevation basis function included in the first set of elevation basis functions at the elevation angles comprises evaluating the Pth elevation basis function at the elevation angles evaluating the p-th elevation basis function at the elevation angle evaluating the p-th elevation basis function at the elevation angle includes finding an index u that satisfies and evaluating the p-th elevation basis function at the elevation angle .

[0198] A11. The method of claim A7 or A8, further comprising obtaining a model representing at least the first set of azimuthal basis functions, wherein the model comprises a sequence (φ1) where designates a sub-interval the azimuthal basis functions are polynomials over the sub-intervals; and a three-dimensional array of model parameters

[0199] A12. The method of claim A11, wherein the first set of azimuthal basis functions comprises a q-th azimuthal basis function, and evaluating each azimuthal basis function included in the first set of azimuthal basis functions at the azimuthal angle includes evaluating the q-th azimuthal basis function at the azimuthal angle and evaluating the q-th azimuthal basis function at the azimuthal angle includes finding an index 1 that satisfies and evaluating the q-th azimuthal basis function at the azimuthal angle .

[0200] A13. The method of any one of claims A7 to A12, wherein the step of obtaining the first set of azimuthal basis function values further comprises generating the first set of azimuthal basis functions.

[0201] A14. The method of claim A13, wherein generating the first set of azimuthal basis functions comprises generating a set of periodic B-spline basis functions over a 0 to 360 degree azimuthal range.

[0202] ​A15. The method of claim A14, wherein generating the set of periodic B-spline basis functions over the 0 to 360 degree azimuth angle range comprises: specifying a sequence of nodes of length L over the 0 to 360 degree range; generating an extended sequence of nodes based on the sequence of nodes of length L, wherein generating the extended sequence of nodes comprises extending the sequence of nodes of length L in a periodic manner using J values less than 0 degrees and J-1 values greater than 360 degrees; obtaining an extended multiplicity sequence of nodes; generating a set of extended B-spline basis functions using the extended sequence of nodes and the extended multiplicity sequence of nodes; selecting L-1 consecutive extended basis functions of the extended basis functions starting from index 2; and mapping the selected extended basis functions to the 0 to 360 degree azimuth angle range in a periodic manner.

[0203] A16. The method of any one of claims A1-A15, further comprising: determining an elevation-azimuth angle determining an interaural time difference

[0204] A17. The method of claim A16, further comprising: determining a right delay based on determining a left delay based on

[0205] A18. The method of claim A17, wherein filtering the sound signal using the right filter comprises filtering the sound signal using the right filter and the right delay ; and filtering the sound signal using the left filter comprises filtering the sound signal using the left filter and the left delay .

[0206] A19. The method of claim A18, wherein filtering the sound signal using the right filter and the right delay comprises calculating filtering the sound signal using the left filter and the left delay comprises calculating where u(n) is the sound signal.

[0207] A20. The method of any one of claims A17-A19, wherein

[0208] ; and

[0209] ​​​​A21. The method of any one of claims A7 to A15, wherein the azimuthal basis functions have a periodicity with a period of 360 degrees.

[0210] While various embodiments have been described herein, it should be understood that they have been described by way of example only and not by way of limitation. Thus, the breadth and scope of the present disclosure should not be limited by any of the above described exemplary embodiments, but should be defined in accordance with the following claims and their equivalents.

[0211] Furthermore, while the processes described above and illustrated in the drawings are performed by a series of steps, it is to be understood that such is only for illustration purposes and therefore should not be viewed as a limitation. It is possible that some steps can be added, some steps can be omitted, the order of some steps can be changed, and some steps can be performed in parallel.

[0212] Abbreviations:

[0213] AR augmented reality

[0214] DOA direction of arrival

[0215] FIR finite impulse response

[0216] HR head-related

[0217] HRIR head-related impulse response

[0218] HRTF head-related transfer function

[0219] ILD interaural intensity difference

[0220] IPD interaural phase difference

[0221] ITD interaural time difference

[0222] ITF interaural transfer function

[0223] MAA minimum audible angle

[0224] MPEG moving picture experts group

[0225] MR mixed reality

[0226] MSE mean squared error

[0227] PCA principal component analysis

[0228] SAOC spatial audio object coding

[0229] SH spherical harmonics

[0230] SOFA spatially oriented acoustic format

[0231] SVD singular value decomposition

[0232] VR virtual reality

[0233] XR extended reality

[0234] References:

[0235] [1] Doris J. Kistler, Frederic L. Wightman, "A model of head-related transfer functions based on principal components analysis and minimum-phase reconstruction," Journal of Acoustical Society of America, 91(3): 1637-1647, March 1992.

[0236] [2] Fabio P. Freeland, Luiz W. P. Biscainho, and Paulo S. R. Diniz, "Interpolation of Head-Related Transfer Functions (HRTFS): A multi-source approach," in 12th European Signal Processing Conference, pp. 1761-1764, Vienna, September 2004.

[0237] [3] Mengqiu Zhang, Rodney A. Kennedy, and Thushara D. Abhayapala, "Empirical determination of frequency representation in spherical harmonics-based HRTF functional modeling," IEEE / ACM Transactions on Audio, Speech, and Language Processing, vol. 23(2), pp. 351-360, Feb. 2015.

[0238] [4] Zamir Ben-Hur, David Lou Alon, Boaz Rafaely, and Ravish Mehra, "Loudness stability of binaural sound with spherical harmonic representation of sparse head-related transfer functions," EURASIP Journal on Audio, Speech, and Music Processing, May 2019, 2019.

Claims

1. A method (1900) for filtering of a sound signal, the method comprising: A filter pair is generated (s1902) for a particular location specified by an elevation angle and an azimuth angle The filter pair consists of a right filter and a left filter ​ filtering (s1904) a sound signal using the right filter; and filtering (s1906) the sound signal using the left filter, wherein generating the filter pair comprises: i) obtaining at least a first set of elevation basis function values at the elevation angle; ii) obtaining at least a first set of azimuth basis function values at the azimuth angle; iii) generating the right filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) right filter model parameters; and iv) generating the left filter using: a) at least the first set of elevation basis function values, b) at least the first set of azimuth basis function values, and c) left filter model parameters, wherein obtaining the first set of elevation basis function values comprises, for each elevation basis function included in a first set of elevation basis functions, evaluating the elevation basis function at the elevation angle to produce an elevation basis function value corresponding to the elevation angle and the elevation basis function, and obtaining the first set of azimuth basis function values comprises, for each azimuth basis function included in a first set of azimuth basis functions, evaluating the azimuth basis function at the azimuth angle to produce an azimuth basis function value corresponding to the azimuth angle and the azimuth basis function, each elevation basis function included in the first set of elevation basis functions is a B-spline basis function, and each azimuth basis function included in the first set of azimuth basis functions is a periodic B-spline basis function.

2. The method of claim 1, wherein 3. The method of claim 1, wherein Obtaining the first set of azimuthal basis function values includes obtaining a set of azimuthal basis function values wherein the set of azimuthal basis function values includes the first set of azimuthal basis function values. generating the right filter comprises the following calculation: generating the left filter comprises the following calculation: and 4. The method of claim 1, wherein wherein , = 1 to , = 1 to , and = 1 to is a right model parameter set, where p is an index of an elevation basis function in a set of elevation basis functions and an index of a set of azimuth basis functions corresponding to the elevation basis function, q is an index of an azimuth basis function in a set of azimuth basis functions, k is an index of a basis vector, P is a number of elevation basis functions in a set of elevation basis functions and a number of sets of azimuth basis functions corresponding to the elevation basis functions, is a number of azimuth basis functions corresponding to the index p, K is a number of basis vectors of an N-dimensional vector space of a filter parameter vector, , = 1 to , = 1 to , and = 1 to is the left model parameter set, , =1 to , define the first set of elevation basis function values at the elevation angle , and , = 1 to , and = 1 to , defines a set of azimuthal basis function values at the azimuthal angle , and , = 1 to is a set of orthonormal basis vectors of length .

5. The method of claim 1, wherein Obtaining the first set of elevation basis function values includes obtaining a set of elevation basis function values wherein the set of elevation basis function values includes the first set of elevation basis function values. generating the right filter comprises the following calculation: generating the left filter comprises the following calculation: ,as well as obtaining a model representing at least the first set of elevation basis functions, wherein the model comprises: wherein , = 1 to , = 1 to , and = 1 to is a right model parameter set, where q is an index of an azimuthal basis function in a set of azimuthal basis functions and an index of a set of elevation basis functions corresponding to the azimuthal basis function, p is an index of an elevation basis function in a set of elevation basis functions, k is an index of a basis vector, Q is a number of azimuthal basis functions in the set of azimuthal basis functions and a number of the set of elevation basis functions corresponding to the azimuthal basis function, is a number of elevation basis functions corresponding to the index q, K is a number of basis vectors of an N-dimensional vector space of a filter parameter vector, , = 1 to , = 1 to , and = 1 to is the left model parameter set, , =1 to , and =1 to , defines a set of elevation basis function values at the elevation angle , and , =1 to , defines the first azimuthal basis function value set at the azimuthal angle and , = 1 to is a set of orthonormal basis vectors of length .

6. The method of claim 1, further comprising:

7. The method of claim 6, wherein Sequence wherein, Designated sub-interval }, the elevation basis function is a polynomial over the sub-interval, and Three-dimensional array of model parameters , Where J is the highest degree of the polynomial coefficients, U is the number of elevation subintervals, P is the number of elevation basis functions, j is the polynomial term index, u is the subinterval index, and p is the elevation basis function index. are the polynomial coefficients of the elevation basis function set model. obtaining a model representing at least the first set of azimuth basis functions, wherein the model comprises: The first set of elevation basis functions includes a first elevation basis function elevation basis functions, at the elevation angle includes evaluating each elevation basis function included in the first set of elevation basis functions at the elevation angle includes evaluating the first elevation basis function at the elevation angle at the elevation angle the first The elevation basis functions comprise the steps of: finding an index that satisfies the index ; and According to The first The value of the elevation basis function at the elevation angle The value of the elevation basis function at the elevation angle 8. The method of claim 1, further comprising:

9. The method of claim 8, wherein Sequence wherein, Designated sub-interval }, the azimuthal basis function being a polynomial over the sub-interval, and Three-dimensional array of model parameters , where J is the highest degree of polynomial coefficients, L1 is the number of azimuth subintervals of the first azimuth basis function set, Q1 is the number of the first azimuth basis function set, j is a polynomial term index, is a subinterval index, and q is a basis function index.

10. The method of any one of claims 1-9, The first set of azimuthal basis functions comprises a first azimuthal basis function azimuthal basis function, evaluating each azimuthal basis function included in the first set of azimuthal basis functions at the azimuthal angle includes evaluating each azimuthal basis function included in the first set of azimuthal basis functions at the azimuthal angle includes evaluating each azimuthal basis function included in the first set of azimuthal basis functions at the azimuthal angle and at the azimuth angle the first The azimuthal basis function comprises the steps of: finding an index that satisfies the index ; and According to evaluating the first azimuth angle function at the azimuth angle of the value. the step of obtaining the first set of azimuth basis function values further comprises generating the first set of azimuth basis functions. wherein generating the first set of azimuth basis functions comprises generating a set of periodic B-spline basis functions over a 0 to 360 degree azimuth range.

11. The method of claim 10, wherein, generating a set of periodic B-spline basis functions over a 0 to 360 degree azimuth range comprises:

12. The method of claim 11, wherein, obtaining an extended sequence of nodes; designate to a sequence of nodes of length in a degree range generating an extended node sequence based on the length node sequence, wherein generating the extended node sequence includes expanding the length node sequence in a periodic manner using values less than values greater than ​​​ generating an extended set of B-spline basis functions using the extended sequence of nodes and the extended sequence of multipliers; 14. The method of claim 13, further comprising: selecting a number of consecutive extended basis functions starting from index 2 in the extended basis functions and The selected extended basis functions are mapped to to degree azimuthal range.

13. The method of any one of claims 1 to 9, further comprising: For elevation-azimuth angles Determining interaural time difference .

15. The method of claim 14, wherein based on determining right delay ; and based on determining the left delay .

16. The method of claim 15, wherein filtering the sound signal using the right filter includes using the right filter and the right delay filtering the sound signal; and filtering the sound signal using the left filter includes using the left filter and the left delay filtering the sound signal. ​ using the right filter and filtering the sound signal comprises calculating ), using the left filter and filtering the sound signal comprises calculating ), wherein, is the sound signal, where u is a signal identifier and n is a discrete time index.

17. The method of any one of claims 1 to 9, wherein, The azimuthal basis functions have a periodicity with a period of 360 degrees.

18. A computer program product comprising instructions (2144) which, when executed by processing circuitry (2102) of a filtering device (2100), causes the filtering device (2100) to perform the method according to any one of claims 1 to 17.

19. A computer-readable storage medium (2142) containing instructions which, when executed by processing circuitry (2102) of a filtering device (2100), causes the filtering device (2100) to perform the method according to any one of claims 1 to 17.

20. A filtering device (2100) for filtering of a sound signal, the filtering device (2100) comprising: processing circuitry (2102); and a memory (2142) containing instructions (2144) executable by the processing circuitry, whereby the filtering device is operatively enabled to perform the method according to any one of claims 1 to 17.

Images