# Patent

## Title
Methods, systems, and computer readable media for utilizing adaptive rectangular decomposition (ARD) to generate head-related transfer functions

## Abstract

Methods, systems, and computer readable media for utilizing adaptive rectangular decomposition (ARD) to perform head-related transfer function (HRTF) simulations are disclosed herein. According to one method, the method includes obtaining a mesh model representative of head and ear geometry of a listener entity and segmenting a simulation domain of the mesh model into a plurality of partitions. The method further includes conducting an ARD simulation on the plurality of partitions to generate simulated sound pressure signals within each of the plurality of partitions and processing the simulated sound pressure signals to generate at least one HRTF that is customized for the listener entity.

## Background
PRIORITY CLAIM
This application claims the benefit of U.S. Provisional Patent Application Ser. No. 62/199,880, filed Jul. 31, 2015; the disclosure of which is incorporated herein by reference in its entirety.


GOVERNMENT INTEREST
This invention was made with government support under Grant Nos. IIS-0917040, IIS-1320644, IIS-1349074 awarded by the National Science Foundation and W911NF-10-1-0506, W911NF-12-1-0430, W911NF-13-C-0037 awarded by the U.S. Army Research Office. The government has certain rights in the invention.


TECHNICAL FIELD
The subject matter described herein relates to sound propagation. More specifically, the subject matter relates to methods, systems, and computer readable media for utilizing adaptive rectangular decomposition to generate head-related transfer functions.
BACKGROUND
Three dimensional (3D) Audio Systems often rely on Head-Related Transfer Functions (HRTFs) to add spatial characteristics to auditory images that the audio systems generate. Industrial implementations use “standard” datasets or use mathematical models to generate approximations of HRTFs, which might generate inaccurate spatialization since HRTFs vary from person to person. For this reason, researchers working on spatial sound or psychoacoustics often make physical measurements in an anechoic chamber to generate HRTFs specific to a person. While this produces better results, the process is expensive and time consuming.
Accordingly, there exists a need for systems, methods, and computer readable media for efficiently generating personalized HRTFs at low cost.
SUMMARY
Methods, systems, and computer readable media for utilizing adaptive rectangular decomposition to generate head-related transfer functions are disclosed herein. According to one method, the method includes obtaining a mesh model representative of head and ear geometry of a listener entity and segmenting a simulation domain of the mesh model into a plurality of partitions. The method further includes conducting an ARD simulation on the plurality of partitions to generate simulated sound pressure signals within each of the plurality of partitions and processing the simulated sound pressure signals to generate at least one HRTF that is customized for the listener entity.
A system for utilizing adaptive rectangular decomposition to generate head-related transfer functions is also disclosed. The system includes a preprocessing engine, an ARD simulation engine, and an HRTF engine, each of which are executable by a processor. In some embodiments, the preprocessing engine is configured to obtain a mesh model representative of head and ear geometry of a listener entity and segment a simulation domain of the mesh model into a plurality of partitions. Likewise, the ARD simulation engine is configured to conduct an ARD simulation on the plurality of partitions to generate simulated sound pressure signals within each of the plurality of partitions. Further, the HRTF engine is configured to process the simulated sound pressure signals to generate at least one HRTF that is customized for the listener entity.
The subject matter described herein can be implemented in software in combination with hardware and/or firmware. For example, the subject matter described herein can be implemented in software executed by one or more processors. In one exemplary implementation, the subject matter described herein may be implemented using a non-transitory computer readable medium having stored thereon computer executable instructions that when executed by the processor of a computer control the computer to perform steps. Exemplary computer readable media suitable for implementing the subject matter described herein include non-transitory devices, such as disk memory devices, chip memory devices, programmable logic devices, and application specific integrated circuits. In addition, a computer readable medium that implements the subject matter described herein may be located on a single device or computing platform or may be distributed across multiple devices or computing platforms.
As used herein, the terms “node” and “host” refer to a physical computing platform or device including one or more processors and memory.
As used herein, the terms “function” and “engine” refer to software in combination with hardware and/or firmware for implementing features described herein.



BRIEF DESCRIPTION OF THE DRAWINGS
Preferred embodiments of the subject matter described herein will now be explained with reference to the accompanying drawings, wherein like reference numerals represent like parts, of which:
FIG. 1 is a block diagram illustrating an exemplary system for utilizing adaptive rectangular decomposition to generate HRTFs according to an embodiment of the subject matter described herein;
FIG. 2 is a block diagram illustrating an exemplary HRTF computational pipeline according to an embodiment of the subject matter described herein;
FIG. 3 is a block diagram illustrating an exemplary mesh model acquisition pipeline according to an embodiment of the subject matter described herein; and
FIG. 4 is a diagram illustrating a flow chart of an exemplary method for utilizing adaptive rectangular decomposition to generate head-related transfer functions according to an embodiment of the subject matter described herein.



DETAILED DESCRIPTION
The human auditory system's ability to localize the direction of incoming sound based on the sound signals received at a subject's ears is attributed to cues such as interaural time difference, interaural intensity difference and spectral modification due to the scattering of sound waves due to the body. Three dimensional sound systems often incorporate these cues into the audio rendering, which is usually accomplished through the use of head related transfer functions (HRTFs).
A significant challenge involving the use of HRTFs is the variation of head, pinna and torso geometries, and the corresponding variation in HRTFs across different individuals. The HRTF measurement techniques that have been traditionally used to obtain personalized HRTFs often require the use of specialized, expensive equipment as well as tedious processes where subjects must remain still for long periods of time. As a result, personalized HRTFs of individuals are very rarely available and virtual auditory displays usually resort to using generic HRTFs. The use of such non-personalized HRTFs can lead to problems, such as lack of externalization, front-back confusions and reversals, incorrect elevation perception, and overall unconvincing spatializations. These difficulties have motivated the need to develop efficient techniques to obtain personalized HRTFs for individuals.
One approach to solving this technical problem is based on the notion that HRTF measurement can be considered to be an acoustic scattering problem in free-field. Given the 3D mesh model of a human body and its acoustic properties, numerical sound simulation techniques can be used to compute HRTFs. Techniques such as the boundary element method and the finite-difference time-domain method may be used to compute HRTFs. The accuracy of these computed HRTFs has been demonstrated by comparing them with measurements. However, these techniques are computationally expensive and can take several hours or days to process.
In some embodiments, the disclosed subject matter presents an efficient technique for computing personalized HRTFs using a numerical simulation technique called adaptive rectangular decomposition (ARD). To reduce computation time, the disclosed system and technique may be configured to use of the acoustic reciprocity principle to reduce number of simulations required and the Kirchhoff surface integral representation (KSIR) to reduce the size of the simulation domain. In some instances, embodiments of the disclosed system and technique may only require approximately 20 minutes of simulation time to compute broadband HRTFs on an eight-core computing device machine compared to hours or days needed by other techniques. Further, the accuracy of the presented approach may be analyzed by computing the left-ear HRTF of the Fritz and KEMAR manikins. For example, the mean spectral mismatch between the HRTF computed by the pipeline disclosed in the subject matter and measurements was 3.88 dB for Fritz and 3.58 dB for KEMAR, within a linear frequency range from 700 Hz to 14 kHz.
The subject matter described herein discloses methods, systems, and computer readable media for utilizing adaptive rectangular decomposition (ARD) to generate head-related transfer functions (HRTFs). Reference will now be made in detail to exemplary embodiments of the subject matter described herein, examples of which are illustrated in the accompanying drawings. Wherever possible, the same reference numbers will be used throughout the drawings to refer to the same or like parts.
FIG. 1 is a block diagram illustrating an HRTF simulation system 101 for generating at least one HRTF that is customized for a listener entity according to an embodiment of the subject matter described herein. As used herein, a listener entity may include a human listener or a virtual listener entity. HRTF simulation system 101 may be any suitable entity (e.g., such as a computing device or platform) for generating a mesh model of head and ear geometry and/or generating an HRTF using a mesh model input and ARD simulation. In accordance with embodiments of the subject matter described herein, components, modules, engines, and/or portions of HRTF simulation system 101 may be implemented in a single computing device or platform, or alternatively distributed across multiple devices or computing platforms.
In some embodiments, HRTF simulation system 101 may comprise a special purpose computing platform that includes a plurality of processors 1021 . . . N that make up a central processing unit (CPU) cluster. In some embodiments, each of processors 102 may include a processor core, a physical processor, a field-programmable gateway array (FPGA), an application-specific integrated circuit (ASIC), and/or any other like processing unit. Each of processors 1021 . . . N may include or access memory 104 in HRTF simulation system 101, such as for storing executable instructions and/or software based constructs. Memory 104 may be any non-transitory computer readable medium and may be operative to communicate with processors 102. Memory 104 may include and/or store a mesh generation engine 106, preprocessing engine 108, an ARD simulation engine 110, an HRTF engine 112, and/or a surface integral formulation engine 114. The functions executed by engines 106-114 are described in greater detail below.
It will be appreciated that FIG. 1 is for illustrative purposes and that various components, their locations, and/or their functions may be changed, altered, added, or removed. For example, some engines and/or functions may be combined into a single entity.
FIG. 2 is a block diagram illustrating an exemplary HRTF computational pipeline process that may be supported and/or executed by system 101. In some embodiments, a 3D mesh model 202 of the head and/or torso of a listener entity/subject is available as input (e.g., a mesh model generated by mesh generation engine 106, or a mesh model generated by a 3D scanning device) for system 101. Embodiments and/or exemplary techniques in which a mesh model may be generated and/or acquired is described in detail below in FIG. 3. In general, system 101 utilizes ARD to perform a sound propagation simulation by solving the acoustic wave equation (see equation 3 below). In some embodiments, ARD may be utilized to divide a simulation domain into grid cells and compute sound wave pressure at each of those grid cells at each time step Compared to finite difference based methods, ARD has much less numerical dispersion error and exhibits the technical advantage of being up to two orders of magnitude faster for homogeneous media. The principle behind ARD's efficiency and accuracy is the use of the exact numerical solution of the acoustic wave equation within rectangular (e.g., cuboidal) domains comprising an isotropic, homogeneous, dissipation-free medium.
For example, in the domain preprocessing stage 206 shown in FIG. 2, system 101 may be configured to initiate the ARD simulation process by generating a rectangular (e.g., cuboidal in three dimensions) decomposition of the computation domain. In some embodiments, this decomposition is generated via preprocessing engine 108 in a series of steps. First, the domain is voxelized to generate a grid of voxels by preprocessing engine 108.
In some embodiments, two ARD simulations (e.g., a simulation for each of the left ear and the right ear) are then executed by ARD simulation engine 110 using this simulation domain. Notably, the principle of acoustic reciprocity is used to reverse the role (and/or position) of source and receivers. For example, the aforementioned receiver positions are designated and used by ARD simulation engine 110 as source positions for these simulations, while the original source positions are designated and used by ARD simulation engine 110 as receiver positions. To prevent reflections from domain boundaries, the simulation domain generated by preprocessing engine 108 is surrounded by perfectly absorbing layer.
The simulations generated by ARD simulation engine 110 produce pressure signals at each grid cell within the simulation domain, including the KSIR surface. The pressure signals at the KSIR surface are used as input by the Kirchhoff surface integral formulation (e.g., executed by surface integral formulation engine 114) to generate pressure signals at the reciprocal receiver positions. These signals are the pressure responses at the ear positions due to the original sources around the head. The signals are then used (e.g., processed and/or executed) by HRTF engine 112 to compute HRTFs using the following equations:




  
    
      
        
          
            H
            L
          
          ⁡
          
            (
            
              θ
              ,
              φ
              ,
              ω
            
            )
          
        
        =
        
          
            
              X
              L
            
            ⁡
            
              (
              
                θ
                ,
                φ
                ,
                ω
              
              )
            
          
          
            
              X
              C
            
            ⁡
            
              (
              
                θ
                ,
                φ
                ,
                ω
              
              )
            
          
        
      
    
    
      
        (
        1
        )
      
    
  
  
    
      
        
          
            
              H
              R
            
            ⁡
            
              (
              
                θ
                ,
                φ
                ,
                ω
              
              )
            
          
          =
          
            
              
                X
                R
              
              ⁡
              
                (
                
                  θ
                  ,
                  φ
                  ,
                  ω
                
                )
              
            
            
              
                X
                C
              
              ⁡
              
                (
                
                  θ
                  ,
                  φ
                  ,
                  ω
                
                )
              
            
          
        
        ,
      
    
    
      
        (
        2
        )
      
    
  




where XL(θ,φ,ω) and XR(θ,φ,ω) respectively represent the Fourier transforms of the left-ear and right-ear time-domain pressure signals for the original source at azimuth θ and elevation φ, and XC(θ,φ,ω) is the Fourier transform of the signal received at the point of origin due to the same source in the absence of the listener, all in free-field conditions.

In general, system 101 (and/or simulation engine 110) utilizes ARD, which is a numerical simulation technique that performs sound propagation simulation by solving the acoustic wave equation (see equation 3 below). Like finite difference based methods, system 101 utilizes ARD to divide the simulation domain into grid cells and computes sound pressure at each of those grid cells at each time step. However, compared to finite-difference-based methods, ARD processing conducted by system 101 has the technical advantage of having a much lower numerical dispersion error while being at least an order of magnitude faster. The principle behind ARD's efficiency and accuracy is system 101's use of the exact analytical solution of the wave equation within cuboidal domains comprising of a homogeneous, dissipation-free medium:




  
    
      
        
          
            p
            ⁡
            
              (
              
                x
                ,
                y
                ,
                z
                ,
                t
              
              )
            
          
          =
          
            
              ∑
              
                i
                =
                
                  (
                  
                    
                      i
                      x
                    
                    ,
                    
                      i
                      y
                    
                    ,
                    
                      i
                      z
                    
                  
                  )
                
              
            
            ⁢
            
              
                
                  m
                  i
                
                ⁡
                
                  (
                  t
                  )
                
              
              ⁢
              
                cos
                ⁡
                
                  (
                  
                    
                      
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        
                          i
                          x
                        
                      
                      
                        l
                        x
                      
                    
                    ⁢
                    x
                  
                  )
                
              
              ⁢
              
                cos
                ⁡
                
                  (
                  
                    
                      
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        
                          i
                          y
                        
                      
                      
                        l
                        y
                      
                    
                    ⁢
                    y
                  
                  )
                
              
              ⁢
              
                cos
                ⁡
                
                  (
                  
                    
                      
                        π
                        ⁢
                        
                          
                        
                        ⁢
                        
                          i
                          z
                        
                      
                      
                        l
                        z
                      
                    
                    ⁢
                    z
                  
                  )
                
              
            
          
        
        ,
      
    
    
      
        (
        3
        )
      
    
  




where p(x,y,z,t) represents the pressure field (or sound signal) at position (x,y,z) and at time t, (lx,ly,lz) are the extents of the cuboidal region, and mi(t) are time-varying mode coefficients. As this solution is composed of cosines, ARD (e.g., as executed by system 101) uses efficient Fast Fourier Transform (FFT) techniques to compute sound propagation within the cuboidal region. Below, each stage and/or engine of the HRTF computational pipeline executed by system 101 is described in detail.

In some embodiments, preprocessing engine 108 in system 101 may be configured to receive a mesh model (e.g., mesh model 202) generated by mesh generation engine 106 and subsequently establish a simulation domain of the mesh model. In other embodiments, the mesh model may be generated through other techniques, such as the use of a 3D scanner device. Description of mesh generation/acquisition techniques performed by mesh generation engine 106 is described below and illustrated in FIG. 3. Preprocessing engine 108 may be configured to execute a preprocessing stages (shown in FIG. 2) including a domain stage 206, voxelization stage 208, and rectangular decomposition stage 210. For example, preprocessing engine 108 may establish a simulation domain utilizing a mesh model 202 as input. In some embodiments, a 3D mesh (e.g., mesh model 202) of the head and torso is positioned by preprocessing engine 108 at the center of an empty cuboidal simulation domain. Preprocessing engine 108 may construct an offset surface around the head of mesh model 202 to serve as a KSIR surface. The size of this domain is selected by preprocessing engine 108 such that the domain closely fits the head and torso mesh model 202, as well as a cuboidal KSIR surface surrounding the head. A point close to each of the blocked ear canal entrances of the mesh model 202 may be designated by preprocessing engine 108 as the receiver positions for the HRTF computation. Further, the source positions are uniformly selected by preprocessing engine 108 at a fixed distance (e.g., one meter) away from the center of the head at different orientations. After the domain is established in stage 206, the domain is voxelized (e.g., decompose the simulation domain by dividing/apportioning the simulation domain into grid cells) in stage 208 by preprocessing engine 108.
For example, preprocessing engine 108 may subsequently be configured to generate a rectangular decomposition of the computation domain. This decomposition may be conducted via preprocessing engine 108 in a series of steps or stages. First, the domain is voxelized to generate a grid of voxels by preprocessing engine 108 (see stage 208). Preprocessing engine 108 may subsequently group the voxels (e.g., grid cells) into the plurality of partitions that include air partitions and perfectly matched layer (PML) partitions, which are separated and/or delineated by interfaces. More specifically, preprocessing engine 108 may subsequently group different voxels and/or grid cells together to form cuboidal regions called air partitions. Boundary conditions are established by preprocessing engine 108, which uses the PML partitions at the boundary to simulate both partially-absorbing and completely-absorbing surfaces. In other embodiments, the air partitions are formed by preprocessing engine 108, which is configured to group the voxels containing the isotropic, homogeneous, dissipation-free medium (e.g., air) together to form rectangular regions (i.e., air partitions). Finally, absorbing boundary conditions are applied by preprocessing engine 108, which uses PML partitions at the boundary to simulate free-field conditions (e.g., as indicated by the HRTF definition).
After the rectangular decomposition processing is conducted by preprocessing engine 108, ARD simulation engine 110 may be configured to conduct an ARD simulation on the plurality of partitions to generate simulated sound pressure signals within each of the plurality of partitions. Furthermore, ARD simulation engine 110 may execute a simulation process that includes using finite difference stencils to propagate sound across the interfaces of adjacent partitions. More specifically, ARD simulation engine 110 may be configured to initiate a number of simulation stages (e.g., current field stage 212, an interfacing handling and discrete cosine transform (DCT) stage 214, and an inverse DCT (IDCT) and modal update stage 216 shown in FIG. 2).
At current field stage 212, ARD simulation engine 110 processes different fields (or portions) of the rectangular decomposed simulation domain. For each of the different fields, ARD simulation engine 110 conducts an interface handling stage 214 in which finite-difference stencils are used to propagate sound across adjacent partitions. In some embodiments, interface handling stage 214 may involve ARD simulation engine 110 being used to propagate sound across two adjacent partitions, which can be either air-air partitions or air-PML partitions.
After conducting stage 214, ARD simulation engine 110 updates the time varying mode coefficients for each air partition based on the acoustic wave equation to propagate sound within partitions and subsequently updates pressure values for each PML partition based on the acoustic wave equation to propagate sound within the plurality of partitions. In some embodiments, ARD simulation engine 110 performs the modal update step by propagating sound within each air partition by updating FFT mode coefficients.
As previously mentioned, HRTFs are functions of source position and require multiple separate recordings of the signal at the ears due to different sound sources placed around the listener. Replicating this process through simulation typically requires multiple separate simulations, one for each source position (e.g., usually in the hundreds). In contrast, system 101 effectively avoids this cost by employing the acoustic reciprocity principle, which provides that the acoustic response remains the same if the sense (e.g., positioning) of source and receiver are reversed. Thus, sources are placed at the receiver positions (e.g., inside the ears) used in HRTF measurement. Similarly, receivers are placed at the various source positions used in HRTF measurement. Thus, system 101 effectively reduces the required number of simulations to only two, one for each ear.
In some embodiments, ARD simulation engine 110 may be further configured to modify the simulation domain of the mesh model to improve processing. For example, ARD simulation engine 110 may utilize surface integration formulation engine 114 to compute (e.g., using a surface integral representation, such as a Kirchhoff surface integral representation) a pressure value at a point outside of the simulation domain using pressure values on a cuboidal surface closely fitting the mesh model. Only pressure values at this surface need be computed by ARD simulation engine 110, thereby reducing the size of the simulation domain as well as computational costs. As a result, surface integration formulation engine 114 and/or ARD simulation engine 110 may output a set of responses that correspond to the mesh model's scattering of Gaussian impulse sound that can be provided to HRTF engine 112 for processing.
In some embodiments, HRTFs may be measured at a fixed distance from the center of the head of the subject. Therefore, in order to compute the full HRTF as described above, a simulation domain with a radius equal to this distance may be used. This distance is usually around 1.0 m (which is much greater than the typical size of the head), due to which the simulation domain is mostly empty as the size of the head and torso is relatively small. Since computation time required by ARD scales cubically with simulation domain dimension, this can lead to large computation times. To reduce the size of the simulation domain, surface integral formulation engine 114 may be configured to make use of the Kirchhoff surface integral representation (KSIR). By using KSIR, surface integral formulation engine 114 may be enabled to conduct the computation of pressure values outside the simulation domain by using pressure values at a tight-fitting surface that encloses the head and torso, resulting in a significantly smaller simulation domain and faster simulations. Notably, surface integral formulation engine 114 can be used to compute the pressure value at a point outside a simulation domain using pressure values on a cuboidal surface closely fitting the mesh. Thus, only pressure values at this surface need to be computed by system 101 and/or surface integral formulation engine 114, thereby significantly reducing the size of the domain as well as the computational cost.
After ARD simulation engine 110 generates simulated sound pressure signals within each of the plurality of partitions of the simulation domain, surface integration formulation engine 114 processes the simulated sound pressure signals. The sound pressure signals (e.g., represented as Fourier transforms of sound waves) are subsequently provided to HRTF engine 112, which may then perform digital signaling processing (DSP). In some embodiments, these HRTFs utilize Fourier transforms of sound pressure signals received at the entrance of the listening entity's left and right blocked ear canals as input variables. In such a scenario, the HRTFs are able to represent the sound signals from a signal as affected by the listener's body (particularly the head, torso, and pinnae of the ear(s) embodied in the mesh model) as measured at the entrance of the listener's ear canals. In addition, HRTF engine 112 may be further configured to determine head related impulse responses (HRIRs) respectively associated to the calculated HRTFs by performing and/or applying an inverse Fourier transform (IFT) on the HRTFs.
For example, ARD simulation engine 110 may be configured to utilize Gaussian impulse sources in the ARD simulations. As such, the output of the KSIR calculation conducted by surface integral formulation engine 114 include a set of responses that correspond to the mesh model's (e.g., head mesh) scattering of Gaussian impulse sound. In order to convert these Gaussian impulse responses to HRIRs, ARD simulation engine 110 utilizes a digital signal processing script that implements equations 1 and 2 presented above. For example, the frequency response of the Gaussian impulse signal at the center of the head in the absence of the head (e.g., XC(θ,φ,ω) in equation 1) is removed from the head responses by this script in the frequency domain, and the HRIR is obtained by ARD simulation engine 110 performing an inverse Fourier transform.
Lastly, in order to perform spatial sound rendering using HRTFs, three steps may be be performed by HRTF engine 112: (a) compute direction of incoming sound field at listener position, (b) model scattering of sound around the listening entity's head using HRTFs, and (c) incorporate listening entity's head orientation. To compute the direction of the incoming sound field at the listener position, system 101 and/or HRTF engine 112 may utilize a plane wave-decomposition approach that uses high-order derivatives of the pressure field at the listener position to compute the plane wave-decomposition of the sound field at interactive rates. Scattering of sound around the head is modeled using the personalized HRTFs computed by HRTF engine 112. Further, HRTF engine 112 may be configured to convert the HRTFs into spherical harmonic basis. By doing this, the listening entity's head rotation can be easily modeled by HRTF engine 112 using standard spherical harmonic rotation techniques. In some embodiments, the spatial sound for each ear can be computed by HRTF engine 112 as a simple dot product of the spherical harmonic coefficients of the plane-wave decomposition and the HRTF. This enables system 101 to generate spatial sound at interactive rates.
FIG. 3 is a block diagram illustrating an exemplary mesh model acquisition pipeline executed by mesh generation engine 106 according to an embodiment of the subject matter described herein. In some embodiments, mesh model 202 depicted in FIG. 2 is generated as output of mesh generation engine 106 depicted in FIG. 3. For example, engine 106 may execute a processing pipeline that includes capturing images (stage 304) of a subject (e.g., listening entity), determining a sparse point cloud (stage 306), generating a noisy mesh model (stage 308), and smoothening the mesh model (stage 310). The mesh model produced as a result from executing stages 304-310 is subsequently sent to an ARD solver 312 (e.g., preprocessing engine 108 in FIG. 1) by mesh generation engine 106.
In some embodiments, mesh generation engine 106 may be configured to generate a 3D mesh model of the head and ear geometry of a listener entity. For example, in stage 304, images of the listener entity's head and ears may be digitally captured (e.g., via a camera and/or video capture device) and subsequently provided to mesh generation engine 106 (e.g., as a set of digital files). In stage 306, mesh generation engine 106 may subsequently perform a Structure-from-Motion process that correlates the captured set of images using one or more distinctive features present in the images. Mesh generation engine 106 may be further configured to generate a sparse point cloud comprising of 3D locations of those distinctive features. In some embodiments, mesh generation engine 106 may be configured to process a set of captured images and compare any neighboring images to each other in order to identify a small set of distinctive “features” (e.g., freckle, mole, scar, etc.) that appear in at least two of the capture images. In some embodiments, multiple images that include a specific feature are taken at different angles (e.g., which are close to each other and can be used to identify the common feature). Notably, the specific feature that is common to the images may be used by mesh generation engine 106 to correlate the multiple images taken.
Next, in stage 308, mesh generation engine 106 may then perform dense modeling of the listener's head and ear geometry based on the sparse point cloud generated by stage 306 as well as the captured images in order to generate a mesh. For example, mesh generation engine 106 may be configured to utilize the sparse point cloud and the camera positions to initiate the generation of a denser mesh that combines all the rest of the parts of the images.
Mesh generation engine 106 may also be configured to apply various mesh cleanup steps (e.g., stage 310) on the mesh model prior to sending the mesh model to preprocessing engine 108 and/or ARD simulation engine 110 for further processing.
In other embodiments involving the generation of personalized HRTFs, mesh generation engine 106 may be configured to obtain accurate head and ear geometry of the user (e.g., stage 302). To facilitate easy acquisition and a highly accurate mesh model, system 101 may also be configured to use digital cameras for the acquisition of the head and ear geometry of the listener entity (e.g., stage 304). In some embodiments, images may be captured by a digital SLR camera (e.g., Canon 60D) with image resolution (3456×2304) and provided to mesh generation engine 106 as input. Such resolution allowed for observing details of the skin texture, which were leveraged by multi-view stereo estimation modules to determine reliable dense correspondences.
In some embodiments, in order to model the area of the head behind the ear (e.g., a critical area for computation of personalized HRTFs), the user may wear concealing headgear (e.g, a swim cap) to hide his or her hair during the data capture. For precise modeling, the user's (e.g., listening entity) head was densely captured all around with samples at approximately every 15 degrees. The selected angular separation between captures affords at least three samples within a 30 degree range, which enables both robust feature matching and precise geometric triangulation. Moreover, this sampling provides sufficient overlap between the views to enable high-accuracy multi-view stereo estimation. Empirically, it was found that sampling intervals larger than 15 degrees may introduce severe aberrations into the resulting 3D model. To increase the model resolution around the ear, 20 or more convergent close-up shots/images were captured for each ear. From the captured images SIFT features were calculated and matched for each image with its top K appearance nearest neighbors, as measured by the GIST descriptor. Using these matches, a structure from motion algorithm was leveraged to perform the incremental structure from motion and bundle adjustment using the cameras internal calibration as provided by the EXIF data of the images. This step provided for the camera registration needed for the dense modeling of the scene.
In some embodiments, dense modeling of the user's head may be performed by mesh generation engine 106 to obtain the desired mesh model required to compute personalized HRTFs. Using a two tier computation that first estimates two-view depths maps was opted. Besides limited accuracy from two view depth maps, highlights on the user's skin occur naturally, which can cause erroneous geometry. In some embodiments, mesh generation engine 106 may further perform smoothening processing on the mesh model (e.g., stage 308). For example, the two view depth maps may be combined by a depth map fusion performed by engine 106, which rejects the erroneous geometry resulting from highlights and produces a noisy mesh model In some embodiments, mesh generation engine 106 may be configured to apply a 3D Delaunay triangulation of dense point clouds and the construction of a graph based on the tetrahedrons from the Delaunay triangulation with weights set according to camera-vertex ray visibility. Mesh generation engine 106 may further refine the graph's t-edge weights and obtain a water-tight dense surface mesh by using a graph-cut based labeling optimization to label each tetrahedron as inside or outside.
Before the generated surface mesh is used as input for the processing pipeline 200 shown in FIG. 2, mesh generation engine 106 may perform some mesh cleanup steps in stage 310. First, since the generated mesh may not be to scale with the subject, mesh generation engine 106 may use the subject's measured head width and head depth (e.g., anthropometric measurements) to scale the generated mesh model. Next, mesh generation engine 106 may be configured to remove stray vertices and triangles from the main head mesh. Further, mesh generation engine 106 may also be configured to perform hole-filling using standard techniques to cover the holes existing in the mesh model. Finally, mesh generation engine 106 may align and orient the mesh model to match the alignment of the head during HRTF measurements and position the head mesh at the center of a cubical simulation domain. Notably, the cuboidal simulation domain (e.g., mesh model) is used as input for an ARD solver 312 (e.g., preprocessing engine 108 shown in FIG. 1). For example, the mesh model generated by mesh generation engine 106 may be embodied as mesh model 202, which is the input depicted in FIG. 2.
In some embodiments, in order to perform the disclosed methods and/or processes, system 101 may be configured to utilize scanned 3D mesh models of a KEMAR (e.g., with DB-60 pinnae) and/or Fritz mankin in order to generate HRTFs. Examples of pertinent simulation parameters that may be utilized by system 101 include the speed of sound within the homogeneous, dissipation-free medium of ARD simulation, which can be set to 343 ms−1 to match that of air. In some embodiments, second-order finite-difference stencils may be used in ARD for interface handling. The maximum simulation frequency for ARD can be set to 88.2 kHz, to have a small grid cell size of 1.94 mm. A Gaussian impulse source with a center frequency of 33.075 kHz can be used as source signal. The absorption coefficient of the mesh surface may be set to 0.02 to correspond to that of human skin. In some embodiments, simulations can be run to generate 5.0 ms pressure signals.
FIG. 4 is a diagram illustrating a flow chart of an exemplary method 400 for utilizing adaptive rectangular decomposition to generate head-related transfer functions according to an embodiment of the subject matter described herein. In block 402, a mesh model that is representative of head and ear geometry of a listener entity is obtained. For example, preprocessing engine 108 may be provided with a closed, accurate 3D mesh model of the head and torso of a listener entity/subject. In some embodiments, the mesh model may be created by mesh generation engine 106 in the manner described above.
In block 404, a simulation domain of the mesh model is segmented into a plurality of partitions. In some embodiments, preprocessing engine 108 uses the mesh model to generate a simulation domain that is subsequently voxelized into grid cells. Preprocessing engine 108 may subsequently group the grid cells into air partitions and/or PML partitions by performing a rectangular decomposition procedure.
In block 406, an ARD simulation is conducted on the plurality of partitions to generate simulated sound pressure signals within each of the plurality of partitions. In some embodiments, ARD simulation engine 110 utilizes the plurality of partitions as constituent rectangles subjected to a sound wave equation. Notably, ARD simulation engine 110 is able to determine the analytical solution of the sound wave equation in any rectangular domain. More specifically, since the spatial portion of the solution of the wave equation is composed of cosines, ARD simulation engine 110 may use a discrete cosine transform to obtain a simulation of the sound wave within a rectangular domain. ARD simulation engine 110 may also employ interfacing handling techniques to process (e.g., simulate) how a sound wave propagates across a boundary/interface between two partitions/rectangles. Using the above information, ARD simulation engine 110 is able to simulate sound pressure signals (e.g., Fourier Transforms of sound pressure waveforms) within each of the plurality of partitions.
In block 408, the simulated sound pressure signals are processed to generate at least one HRTF that is customized for the listener entity. In particular, HRTF engine 112 may receive the sound pressure signal as Fourier transform representations and calculate at least one HRTF. For example, HRTF engine 112 may receive i) Fourier transforms of the left-ear and right-ear time-domain sound pressure signals and ii) the Fourier transform of the signal received at the origin of the mesh model due to the same source in the absence of the listener and compute the HRTFs for the left and right ears using equations (1) and (2) listed above.
It should be noted that HRTF simulation system 101 and/or functionality described herein can constitute a special purpose computing system. Further, HRTF system 101, engines 106-112, and/or functionality described herein provides improvements toward the technological field of acoustic simulation. In particular, HRTF simulation system 101 presents a novel device and algorithm for performing efficient personalized HRTF computations that can be used to simulate high-fidelity spatial sound as perceived by a single listener entity. Notably, the present subject matter presents an advantageous alternative to (and/or obviates the need for) conducting physical measurements of subjects (e.g., in an anechoic chamber) to generate subject-specific HRTFs. Notably, these types of customized solutions can be both cost prohibitive and time consuming.
It will be understood that various details of the subject matter described herein may be changed without departing from the scope of the subject matter described herein. Furthermore, the foregoing description is for the purpose of illustration only, and not for the purpose of limitation, as the subject matter described herein is defined by the claims as set forth hereinafter.

## Claims

1. A method for utilizing adaptive rectangular decomposition (ARD) to generate a head-related transfer function (HRTF), the method comprising:
obtaining a mesh model representative of head and ear geometry of a listener entity, wherein obtaining the mesh model includes capturing images of a head and ear geometry of the listener entity and conducting dense modeling processing on the captured images to generate a three-dimensional (3D) mesh model;
segmenting a simulation domain of the mesh model into a plurality of partitions;
conducting an ARD simulation on the plurality of partitions to generate simulated sound pressure signals within each of the plurality of partitions; and
processing the simulated sound pressure signals to generate at least one HRTF that is customized for the listener entity.

2. The method of claim 1 comprising utilizing the at least one generated HRTF to render spatial sound to the listener entity.

3. The method of claim 1 wherein the at least one HRTF includes a first HRTF and a second HRTF respectively associated with a right ear and a left ear of the listener entity.

4. The method of claim 1 comprising voxelizing the simulation domain of the mesh model into grid cells, wherein the grid cells are subsequently grouped into the plurality of partitions.

5. The method of claim 1 wherein the sound pressure signals are subjected to a surface integral representation after the ARD simulation.

6. The method of claim 1 wherein at least one head-related impulse response (HRIR) customized for the listener entity is determined by applying an inverse Fourier Transform (IFT) to the at least one generated HRTF.

7. A system for utilizing adaptive rectangular decomposition (ARD) to perform head-related transfer function (HRTF) simulations, the system comprising:
a processor;
a preprocessing engine executable by the processor, wherein the preprocessing engine is configured to obtaining a mesh model representative of head and ear geometry of a listener entity and segmenting a simulation domain of the mesh model into a plurality of partitions;
a mesh generation engine configured to capture images of a head and ear geometry of the listener entity and conduct dense modeling processing on the captured images to generate a three-dimensional (3D) mesh model;
an ARD simulation engine executable by the processor, wherein the ARD simulation engine is configured to conduct an ARD simulation on the plurality of partitions to generate simulated sound pressure signals within each of the plurality of partitions; and
an HRTF engine executable by the processor, wherein the HRTF engine is configured to process the simulated sound pressure signals to generate at least one HRTF that is customized for the listener entity.

8. The system of claim 7 wherein the preprocessing engine is further configured to utilize the at least one generated HRTF to render spatial sound to the listener entity.

9. The system of claim 7 wherein the at least one HRTF includes a first HRTF and a second HRTF respectively associated with a right ear and a left ear of the listener entity.

10. The system of claim 7 wherein the ARD simulation engine is further configured to voxelize the simulation domain of the mesh model into grid cells, wherein the grid cells are subsequently grouped into the plurality of partitions.

11. The system of claim 7 wherein the sound pressure signals are subjected to a surface integral representation after the ARD simulation.

12. The system of claim 7 wherein the HRTF engine is further configured to generate at least one head-related impulse response (HRIR) customized for the listener entity is determined by applying an inverse Fourier Transform (IFT) to the at least one generated HRTF.

13. A non-transitory computer readable medium having stored thereon executable instructions that when executed by a processor of a computer cause the computer to perform steps comprising:
obtaining a mesh model representative of head and ear geometry of a listener entity, wherein obtaining the mesh model includes capturing images of a head and ear geometry of the listener entity and conducting dense modeling processing on the captured images to generate a three-dimensional (3D) mesh model;
segmenting a simulation domain of the mesh model into a plurality of partitions;
conducting an ARD simulation on the plurality of partitions to generate simulated sound pressure signals within each of the plurality of partitions; and
processing the simulated sound pressure signals to generate at least one HRTF that is customized for the listener entity.

14. The computer readable medium of claim 13 comprising utilizing the at least one generated HRTF to render spatial sound to the listener entity.

15. The computer readable medium of claim 13 wherein the at least one HRTF includes a first HRTF and a second HRTF respectively associated with a right ear and a left ear of the listener entity.

16. The computer readable medium of claim 13 comprising voxelizing the simulation domain of the mesh model into grid cells, wherein the grid cells are subsequently grouped into the plurality of partitions.

17. The computer readable medium of claim 13 wherein the sound pressure signals are subjected to a surface integral representation after the ARD simulation.