This repository contains the accompanying code and audio files for the Journal of The Acoustical Society of America paper: '3D Reflector Localisation and Room Geometry Estimation using a Spherical Microphone Array', Available at: INSERT THE URL HERE. These programs and audio files are distributed in the hopes that they will prove useful under the Creative Commons Attribution 4.0, with no warranty; or the implied warranty of merchantability or fitness for a particular problem. Please give appropriate credit for use of the material provided in this repository back to the author. 

All code has been tested using MATLAB R2018a.

Please read the corresponding paper for explanation of the whole process, and analysis of the results, Available at: https://doi.org/10.1121/1.5130569

---

The repository contents are as follows:

---

Folder Scenario One

MATLAB Scripts:

1.) HigherOrderAmbi_SmallerCuboid.m - This MATLAB script is used to run the test case for the Cuboid-Shaped Room. This scripts loads in the SRIRs 'Smaller_Cuboid_Normal_A0_01_3rd.WAV' and `Smaller_Cuboid_Normal_A1_01_3rd.WAV', and truncates them the 1500 samples in length. The SRIR are then temporally shifted to align the arrival of the direct sound with that expected given the source-to-receiver distance, and the Eigenbeam Detection and Evaluation of Simultaneously Arriving Reflections method is used to detect reflections within both SRIRs. From these detected reflections the Acoustic Reflection Cartographer method is used to infer the geometry of the room. The resulting plots are as follows: Figure 1 are the detected reflection for SRIR 1, Figure 2 are the detected reflections for SRIR 2, Figure 100 are all candidate boundaries inferred from the image-sources, Figure 3 shows the resulting inferred geometry, and Figure 4 show all of the candidate image-sources inferred. The key output variables are: echos - a structure containing all of the detected reflections for all SRIRs, inferredPlanes - A structure containing the boundaries that have been inferred by the proposed geometry method, and imageSources - a matrix containing all of the inferred candidate image-sources.

2.) HigherOrderAmbi_OctagonalPrism.m - This MATLAB script is used to run the test case for the Octagonal-Shaped Room. This scripts loads in the SRIRs 'Simple_Octagon_Normal_A0_01_3rd.WAV' and `Simple_Octagon_Normal_A2_01_3rd.WAV', and truncates them the 1500 samples in length. The SRIR are then temporally shifted to align the arrival of the direct sound with that expected given the source-to-receiver distance, and the Eigenbeam Detection and Evaluation of Simultaneously Arriving Reflections method is used to detect reflections within both SRIRs. From these detected reflections the Acoustic Reflection Cartographer method is used to infer the geometry of the room. The resulting plots are as follows: Figure 1 are the detected reflection for SRIR 1, Figure 2 are the detected reflections for SRIR 2, Figure 100 are all candidate boundaries inferred from the image-sources, Figure 3 shows the resulting inferred geometry, and Figure 4 show all of the candidate image-sources inferred. The key output variables are: echos - a structure containing all of the detected reflections for all SRIRs, inferredPlanes - A structure containing the boundaries that have been inferred by the proposed geometry method, and imageSources - a matrix containing all of the inferred candidate image-sources.


3.) HigherOrderAmbi_LShape.m - This MATLAB script is used to run the test case for the L-Shaped Room. This scripts loads in the SRIRs 'L_Shaped_Room_Normal_A0_01_3rd.WAV' and `L_Shaped_Room_Normal_A1_01_3rd.WAV', and truncates them the 2000 samples in length. The SRIR are then temporally shifted to align the arrival of the direct sound with that expected given the source-to-receiver distance, and the Eigenbeam Detection and Evaluation of Simultaneously Arriving Reflections method is used to detect reflections within both SRIRs. From these detected reflections the Acoustic Reflection Cartographer method is used to infer the geometry of the room. The resulting plots are as follows: Figure 1 are the detected reflection for SRIR 1, Figure 2 are the detected reflections for SRIR 2, Figure 100 are all candidate boundaries inferred from the image-sources, Figure 3 shows the resulting inferred geometry, and Figure 4 show all of the candidate image-sources inferred. The key output variables are: echos - a structure containing all of the detected reflections for all SRIRs, inferredPlanes - A structure containing the boundaries that have been inferred by the proposed geometry method, and imageSources - a matrix containing all of the inferred candidate image-sources.


4.) HigherOrderAmbi_TShaped.m - This MATLAB script is used to run the test case for the T-Shaped Room. This scripts loads in the SRIRs 'Non_Convex_Room_2_Normal_A0_3_01_3rd.WAV', `Non_Convex_Room_2_Normal_A1_3_01_3rd.WAV', and 'Non_Convex_Room_2_Normal_A2_3_01_3rd.WAV', and truncates them the 2000 samples in length. The SRIR are then temporally shifted to align the arrival of the direct sound with that expected given the source-to-receiver distance, and the Eigenbeam Detection and Evaluation of Simultaneously Arriving Reflections method is used to detect reflections within both SRIRs. From these detected reflections the Acoustic Reflection Cartographer method is used to infer the geometry of the room. The resulting plots are as follows: Figure 1 are the detected reflection for SRIR 1, Figure 2 are the detected reflections for SRIR 2, Figure 100 are all candidate boundaries inferred from the image-sources, Figure 3 shows the resulting inferred geometry, and Figure 4 show all of the candidate image-sources inferred. The key output variables are: echos - a structure containing all of the detected reflections for all SRIRs, inferredPlanes - A structure containing the boundaries that have been inferred by the proposed geometry method, and imageSources - a matrix containing all of the inferred candidate image-sources.


Folder Additional_Functions

1.) multiNorm.m - This MATLAB function computes the distance between a vector and a matrix, or two matrices. Input: Matrix/Vector 1 and Matrix/Vector 2. Output: Matrix of distances.

2.) plotLine3D.m - Plots a line in 3D (used to simplify plot calls) - Input: Matrix containing the line coordinates, Colour of the line (as defined by MATLABs plot commands), and the width of the line.

3.) plotPlane.m - Plots a plane in 3D (used to simplify plot calls) - Input: Matrix containing the coordinates for the corners of the boundary, Colour of the line (as defined by MATLABs plot commands), and the width of the line.

4.) plotPoint3D.m - Plots a matrix of points in 3D (used to simplify plot calls) - Input: Matrix containing the coordinates for each point to be plotted, Colour of the line (as defined by MATLABs plot commands), and the width of the line.

Folder EDESAR

1.) EDESAR.m - This is the main function for the Eigenbeam Detection and Evaluation of Simultaneously Arriving Reflections method for reflection detection. This function analyses short time-frames of the SRIR, and computes the directional spectrum for candidate time-frames that pass defined thresholds (as discussed in the paper) using either the Minimum-Variance Distortionless Resposne Beamformer (real-world data) or a plane wave decomposition beamformer, the steered response power map (simulated data). The directional Spectrum for each time-frame is then analysed to find candidate reflections arriving at the microphone array. Input: IR - Spherical harmonic domain SRIR measurement, Fs - The sampling frequency used when capturing the SRIR, c - The speed-of-sound, receiverRotation - Matrix containing the azimuth rotation (in degrees) needed to ensure that 0 degrees aligns with the positive going x-axis, receiverNumber - The number of the receiver being analysed (only relevant if it is desired that the reflections detected across multiple SRIR be kept seperate), isSimulated - Boolean used let the function know if the data is simulated or not 1 = simulated 0 = real-world, echos - Optional structure containing already detected reflections. Output: echos - Structure containing information about each detected reflection and detectedReflections - a vector of NaNs with the sample value of any detected reflections at the corresponding sample.

2.) sphMVDR.m - Implementation of the spherical harmonic domain Minimum-Variance Distortionless Response beamformer (based on: [1,2]). This function computes the directional spectrum for a real-world measured SRIR. Input: s - Signal being analysed, Fs - Sapling Frequency, and format - spherical harmonic format (only supports Furse-Malham at the moment 'FuMa'). Output: directionalSpectrum - The resulting directional spectrum for the time-frame of the SRIR.

3.) sphSteeredResponsePowerMap.m - Implementation of the spherical harmonic domain steered response power map beamformer (based on: [2,3]). This function computes the directional spectrum for a real-world measured SRIR. Input: s - Signal being analysed, Fs - Sapling Frequency, and format - spherical harmonic format (only supports Furse-Malham at the moment 'FuMa'). Output: directionalSpectrum - The resulting directional spectrum for the time-frame of the SRIR.

Folder External_Code

1.) SPH_Toolbox This folder contains the required functions from [2] needed to run the code in the repository. getDiffuseness_CMD.m Computes the diffuseness profile of a signal using the Covariance Matrix Eigenvalue Diffuseness Estimation method [4], getSH.m - Computes the real- or complex-values spherical harmonics used to steer the response of the array towards a specific direction, grid2dirs.m - Creates a matrix of spherical grid points, and SPH_Toolbox_Licence.txt - The licence for the code. More information about these functions can be found at [2].

Folder Geometry Inference

1.) runGeometryInference.m - This is the main function for the geometry inference method, run this function to analyse the set of candidate detected reflections and infer the geometry of the room. Inputs: echos - Structure containing all the detected reflections for each SRIR, c- The speed-of-sound, receiverNumberList - A vector containing a number assigned to each receiver. Output: inferredPlanes - The inferred boundaries that ideally define the desired room, and echos - The structure containing all of the detected reflections (false-positives having been removed).

2.) computeDirectionalCosines.m - Functions that computes the directional cosines for a line. Input: line - Matrix containing the coordinates for the first and last point in the line. Output: dirCos - The directional cosines for the line. 

4.) constrainPlane.m - Constrain a boundary so that the corners of the boundary are at the intersection points between the two closest non-parallel boundaries. Input: intersections - Matrix containing all the intersection points between every candidate boundary, intersectionDirectionVector - The direction vector for each intersection, planeNumber - The number for the boundary being constrained (corresponds to the index in intersections that belongs to the boundary being tested), imageSource - The image-source location for the boundary being constrained, receiver - The receiver location that the reflection that produced the candidate boundary was detected for, maxZ - z-axis coordinate for the ceiling, and minZ - z-axis coordinate for the floor. Output: plane - The constrained boundary.

5.) constrainPlane2.m - Same as above but with the constraint that the reflection path ideally should pass through the constrained boundary.

6.) constrainRoomFromPlanes.m - This function is the main processing of the geometry inference method, and refines the candidate boundaries to ideally that of the desired room. This is achieved using the geometry validation process outlined in the paper. Input: echos - As before, unconstrainedPlane - A structure containing each candidate boundaries corners, planeInformation - Structure containing the echoNumber, receiverNumber, prevSource coordinates, prevSource number and image source list that created each of the candidate boundaries, planeNormal - Structure containing the boundary normal for each of the candidate boundaries, pointOnPlane - Structure containing a point on each candidate boundary, maxZ, minZ. Output: inferredPlanes - Structure containing the remaining boundary corners that ideal make up the desired room.

7.) findPreviousSource.m - This MATLAB function is used to find the most-likely previous-source for each image-source. Input: echos, receiverNumberList. Output: echos, confirmedImageSources - Matrix containing the image-sources for which a previous-source has been found, maxZ, minZ

8.) generateFloorCeiling.m - Function that generates the floor and ceiling from the corners of every candidate boundary (after constrainRoomFromPlanes.m has been run). Input: inferredPlanes, planeNo - Number of planes + 1. Output: inferredPlanes.

9.) generateInferredPlanes.m - This function is used to constrain every candidate inferred boundary based on the intersections between nearest non-parallel boundaries. Input: unconstrainedPlane - Structure containing all of the candidate inferred boundaries, pointOnPlane, planeInformation, echos, maxZ, minZ.

10.) generateInferredPlanes2.m - As above but with the constraint that the reflection path ideally should pass through the constrained boundary. 

11.) generatePlane.m - Generate a boundary from a point on the boundary and the boundary's normal. Input: pointOnPlane - Vector containing the point of the boundary, planeNormal - Vector containing the boundary's normal vector. Output: plane: The corners of the resulting boundary.

12.) generateUnconstrainedPlanes.m - This functions generates all of the candidate boundaries from the image-source previous-source combinations. Input: echos. Output: echos, unconstrainedPlane, planeInformation, planeNormal, and pointOnPlane.

13.) inferGeometry.m - This function plots all of the unconstrained planes that have been inferred, calls the function that removes coincident boundaries, and calls the function `constrainRoomFromPlanes.m'. Input: echos, maxZ, and minZ. Output: inferredPlanes.

14.) inferImageSource.m - This function generates all candidate image-source coordinates based on the time- and direction-of-arrival for each detected reflection. Input: echos, c, and receiverNumberList. Output: echos, imageSources - Coordimates for each candidate image-source, receiverNumberList

15.) linePlaneIntersection.m - Compute the intersection between a boundary and line, when both are assumed to be infinite in length. Input: ray - Start and end point for a ray, plane - Corners of the boundary. Output: incidence - Coordinates for the point of incidence.

16.) linePlaneIntersection_Constrained.m - As above but constrained so the boundary and line are constrained to the length of the line and plane provided as input.

17.) linePlaneIntersection_Constrained_rayNotBound.m - As above with the plane constrained to be the length provided, but the line is infinite in length.

18.) newDirCos.m - Compute the directional cosines for a ray being reflected off a boundary specularly. Input: plane, ray, prevDirCos - The directional cosines of the ray going towards the boundary, planeNormal. Output: dirCos.

19.) planePlaneIntersection.m - Compute the intersection points between two non-parallel boundaries. Input: plane1 - Corners of the first boundary, pointOnPlane1 - Point on first boundary, plane2 - Corners of second boundary, pointOnPlane2 - Point on second boundary, maxZ, and minZ. Output: intersection,intersectionDirectionVector/

20.) receiverIntersection.m - Compute the intersection between a ray and a sphere around the receiver location. Input: receiver, ray. Output: DoA - The azimuth and elevation DoA of the ray arriving at the receiver.

21.) removePlanes_Coincident.m - This function is used to remove boundaries that are coincident - retaining only one. Input: echos, unconstrainedPlane, planeInformation, planeNormal, pointOnPlane. Output: echos, unconstrainedPlane, planeInformation, planeNormal, pointOnPlane.

22.) removePlanes_InteriorPathwayInvalidation.m - This function analyses all reflection paths that correspond to an inferred candidate boundary, removing any boundaries whose reflections paths are occluded by a boundary closer to the receiver. Input: unconstrainedPlane, planeNo, echos, inferredPlanes, planeInformation, pointOnPlane, planeNormal, planeList, and changeHappened - Boolean defining if any changes have happened when this function was run. Output: unconstrainedPlane, planeNo, echos, inferredPlanes, planeInformation, pointOnPlane, planeNormal, and changeHappened.

23.) removePlanes_lineOfSight.m - This function is used to remove any boundaries that are not in line-of-sight of at least one receiver location. Input: echos, inferredPlanes, unconstrainedPlane, planeList, planeInformation, pointOnPlane, and planeNormal. Output: echos, inferredPlanes, unconstrainedPlane, planeList, planeInformation, pointOnPlane, and planeNormal.

24.) removePlanes_notConnected.m - This function removes any boundaries that are not connected too at least two other boundaries. Input: echos, inferredPlanes, unconstrainedPlane, planeList, planeInformation, pointOnPlane, and planeNormal. Output: echos, inferredPlanes, unconstrainedPlane, planeList, planeInformation, pointOnPlane, and planeNormal.

Folder SRIR

This folder contains all of the SRIR used in this scenario, as described above.

---

Folder Scenario Two

1.) boundaryCombinationsLShaped.mat - MAT file containing the boundary combinations for each source-receiver pair for the first L-Shaped room. (These are used to compute the error metrics for each boundary).

2.) boundaryCombinationsLShaped2.mat - MAT file containing the boundary combinations for each source-receiver pair for the second L-Shaped room. (These are used to compute the error metrics for each boundary).

3.) L_Shaped_Multiple_SR_Combinations.m - This MATLAB scripts runs the geometry inference process for each source/receiver pair for the first L-Shaped room. To speed up the process the SRIR have been pre-analysed and the folder `Analysed_Reflection_Data' contains the `echos' structure for each measurement set-up. The structure outputResults contains the analysis for each measurement set-up.

4.) L_Shaped_Multiple_SR_Combinations2.m - This MATLAB scripts runs the geometry inference process for each source/receiver pair for the second L-Shaped room. To speed up the process the SRIR have been pre-analysed and the folder `Analysed_Reflection_Data' contains the `echos' structure for each measurement set-up.  The structure outputResults contains the analysis for each measurement set-up.

Folders EDESAR, Additional_Functions, External_Code, and Geometry_Inference are as described above.

---

Folder Scenario Three

1.) RealWorld_GeometryInference_Set1.m - This MATLAB script is used to run the test case for the T-Shaped Room. This scripts loads in the SRIRs 'Speaker_Position_One_EigenMike_TakeOne_3rd_Order_Impulse_Response_Set_1_Averaged.wav' and `Speaker_Position_Two_EigenMike_TakeOne_3rd_Order_Impulse_Response_Set_1_Averaged.wav', and 'Non_Convex_Room_2_Normal_A2_3_01_3rd.WAV', and truncates them the 2000 samples in length. The SRIR are then temporally shifted to align the arrival of the direct sound with that expected given the source-to-receiver distance, and the Eigenbeam Detection and Evaluation of Simultaneously Arriving Reflections method is used to detect reflections within both SRIRs. From these detected reflections the Acoustic Reflection Cartographer method is used to infer the geometry of the room. The resulting plots are as follows: Figure 1 are the detected reflection for SRIR 1, Figure 2 are the detected reflections for SRIR 2, Figure 100 are all candidate boundaries inferred from the image-sources, Figure 3 shows the resulting inferred geometry, and Figure 4 show all of the candidate image-sources inferred. The key output variables are: echos - a structure containing all of the detected reflections for all SRIRs, inferredPlanes - A structure containing the boundaries that have been inferred by the proposed geometry method, and imageSources - a matrix containing all of the inferred candidate image-sources.

All sub-folders are as described for Scenario One.

---

References:

[1] D. P. Jarrett, E. A. Habets, and P. A. Naylor, "Spherical harmonic domain noise reduction using an MVDR beamformer and DOA-based second-order statistics estimation," ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings 654-658 (2013) doi: 10.1109/ICASSP.2013.6637729

[2] A. Politis, "Spherical Array Porcessing Toolbox," (2016), https://github.com/polarch/Spherical-Array-Processing.

[3] D. P. Jarrett, E. A. P. Habets, and P. A. Naylor, "3D Source localization in the spherical harmonic domain using a pseudointensity vector," in Signal Processing Conference, 2010 18th European, 5, IEEE, Aalborg, Denmark (2010), pp. 442-446, https://ieeexplore.ieee.org/document/7096575.

[4] N. Epain and C. T. Jin, "Spherical Harmonic Signal Covariance and Sound Field Diffuseness," 24(10), 1796-1807 (2016) doi: 10.1109/TASLP.2016.2585862

----

All code and audio (unless otherwise stated) produced by: Michael Lovedee-Turner, PhD candidate in Music Technology at the Audio Lab, Department of Electronic Engineering, University of York. Corresponding distribution licences are provided for code that is not produced by the author.

Contact: mjlt500@york.ac.uk