Edge diffraction toolbox Version EDB2toolbox Peter Svensson 24 April 2015 (svensson@iet.ntnu.no) OVERVIEW This toolbox contains a set of Matlab functions that calculate the impulse response (IR), and/or the transfer function (TF), for a point source in an environment of rigid, plane surfaces. Specular reflections (up to order 5-10, depending on number of planes and memory capability of your computer) and edge diffraction components (up to sixth order for IR, and up to 2nd order for TF) are included, also combinations of them to some extent. A single program, EDB2main, does all the calculations. It creates a number of intermediate result files, but also output files where IRs and/or TFs (total IR, direct sound IR, specular reflection IR, diffraction IR) are stored. Many sources and receivers can be run as a batch job, and one output file is stored for each source-receiver pair. The final IR could be transformed to a frequency response using FFT. Numerically, the TF result, and the FFT[ IR ], are very very similar for low frequencies. For increasing frequencies, the aliasing effect causes the IR to have a small error which grows with frequency. The aliasing error appears anytime we convert a continuous-time (analytical) expression to a sampled version: we must low-pass filter the analytical expression before we sample it. In the ED toolbox (all versions), the integration ±0.5 sample width before the sampling is a very crude LP filter. An ideal brickwall filter would remove the aliasing error entirely, but then the sampled IR would have infinite ripple before and after the extent of the IR which was finite in the continuous-time domain. Something more advanced than ±0.5 sample integration could be implemented for suppressed aliasing error. The IR can also be convolved with any source signal to get the sound pressure signal in a receiver position. Also, animations of sound propagation can be created quite easily by specifying a regular mesh of receivers and plot the instantaneous sound pressure time step by time step. Such instantaneous plots can then be put together into a movie using the supplied function EDB2makemovie. The calculation method use two distinct steps. Step 2 can give either time-domain calculations (2a) and/or frequency-domain (2b). 1. If the setupfile contains the line findpaths = 1; then all valid reflection/diffraction paths are found. Here, the classical image source method [Borish 1984] is used, with extensions for edge diffraction [Torres et al 2001]. 2a. If the setupfile contains the line calcirs = 1; each valid path generates an impulse response which is added to a total IR. Specular reflections give a pulse whereas diffraction IRs uses the expressions in [Svensson et al 1999] and [Svensson and Calamia 2006]. 2b. If the setupfile contains the line calctfs = 1; each valid path generates a transfer function, which is added to a total TF. Specular reflections give a contribution on the form ±exp(-jkR)/R for Neumann and Dirichlet surfaces, i.e., no absorption. Diffration TFs use the expressions in [Svensson et al 2009] but they are sensitive to shadow boundary singularities. As described in the section on Calculation steps below, the intermediate output files make it relative straightforward to use your own alternative method for any of the steps, as long as you generate a file which contains the right variables. In the setup file it is possible to specify, as an option, any input file with your own filename instead of the automatically generated ones: xxx_eddata for instance. CALCULATION TIME/MEMORY USE The code uses vectorizing to a very large extent in order to speed up the processing. This has the drawback that large matrices are created as intermediate variables inside some of the programs, especially EDB2SorRgeo, EDB2edgeo2x and EDB2findISEStree. There are some possibilities to limit the vectorizing in the EDB2SorRgeo and introduce for-loops but there are no such easy possibilities for the EDB2edgeo2x and EDB2findISEStree files. Consequently, higher than order 6 for specular reflections often leads to the "Out of memory" message even for small models. The file called xxx_ISEStree.mat contains all the potentially valid image sources and edge sources (that is, those that are on the right side of the last reflection plane) but visibility (that is, a check if the hit points are inside the finite planes) and obstruction (a check if there is no obstructing plane for any of the paths) is not tested until the program EDB2findpathsISESx is run. Therefore, the "ISEStree" is very large, and some huge matrices are not transferred via a file between the EDB2findISEStree and EDB2findpathsISESx programs but rather via global variables. The use of a beam tracing algorithm would make these steps much more efficient than the present implementation because the ISEStree could be reduced to a large degree. The calculation of the IRs is also somewhat time consuming but it should be noted that the calculation time depends on the sampling frequency chosen. Also, there is a compiled mex version available for Mac OSX for which the generation of first-order diffraction IRs is substantially faster than the Matlab version. CALCULATION STEPS The calculation flow of the main program, EDB2main, is described below. First, two input files have to be created. Both types are given in the attached example files. Step 1: You have to create a geometry file of the CAD format (specified by the CATT Acoustic software). This file is a simple corner-and-planes definition. Step 2: You have to create a setup file which assigns values to a number of calculation parameters, e.g., setting the name of the CADfile etc. Step 3: Run the EDB2main program, with the setup file name as input parameter. Then the following steps are run (note that a so-called Filestem, xxx, is used for all output files). NB! The brief description of each file is taken from the Contents file in the toolbox. EDsetupfile: zzz.m CADfile: yyyy.cad | | v v ----------------------------------- | | | EDB2readcad | | Reads a file of type .CAD and | saves all the geometry data in | a mat-file | | ----------------------------------- | v cadgeofile: xxx_cadgeo.mat EDsetupfile: zzz.m | v ----------------------------------- | | | EDB2readsetup | | Runs the EDsetup m-file and | saves all the settings in a | mat-file. | | ----------------------------------- | v xxx_setup.mat ________________________________________________ cadgeofile: xxx_cadgeo.mat | v ----------------------------------- | | | EDB2edgeox | | Calculates some plane- and | edge-related geometrical | parameters | | ----------------------------------- | v eddatafile: xxx_eddata.mat | v ----------------------------------- | | | EDB2SorRgeo | | Calculates some source- or | receiver-related geometrical | parameters | | ----------------------------------- | v sdatafile: xxx_sdata.mat rdatafile: xxx_rdata.mat ________________________________________________ (The step from here to the next horizontal line is done only for diffraction order 2 and higher) eddatafile: xxx_eddata.mat | v ----------------------------------- | | | EDB2ed2geox | | Calculates 2nd- and higher-order | edge-related geom. parameters | | ----------------------------------- | v eddatafile: xxx_eddata.mat ________________________________________________ if calcpaths == 1 eddatafile: xxx_eddata.mat | v ----------------------------------- | | | EDB2findISEStree | | Constructs a list ("an ISES | tree") of all possible specular- | diffraction combinations | | ----------------------------------- | v ISEStreefile: xxx_SS_ISEStree.mat (SS is source number) | v ----------------------------------- | | | EDB2findpathsISESx | | Finds all possible paths that | include direct sound, specular, | diffraction | | ----------------------------------- | v edpathsfile: xxx_SS_RR_edpaths.mat (RR is receiver number) ________________________________________________ if calcirs == 1 edpathsfile: xxx_SS_RR_edpaths.mat (RR is receiver number) | v ----------------------------------- | | | EDB2makeirs | | Constructs IRs from a list of | paths in an edpathsfile | | ----------------------------------- | v irfile: xxx_SS_RR_ir.mat ________________________________________________ if calctfs == 1 edpathsfile: xxx_SS_RR_edpaths.mat (RR is receiver number) | v ----------------------------------- | | | EDBmaketfs | | Constructs TFs from a list of | paths in an edpathsfile | | ----------------------------------- | v tffile: xxx_SS_RR_tf.mat #################################################### #################################################### Post-processing the results will include loading the results file and do whatever kind of processing. This can be done in a practical way by including the code below at the end of each setup file. Once calculations are finished, one can then type in the command "doplot = 1;" in the command window, and then execute the setup file. That way, the loading is done automatically, and one can insert plot commands or other processing etc. if exist('doplot') == 1 if doplot == 1 nreceivers = size(receivers,1); results = []; for ii = 1:nreceivers II = int2str(ii); eval(['load ',Filepath,Filestem,'_1_',II,'_tf.mat']) results = [results; tfdirect tfgeom tfdiff]; end end end #################################################### #################################################### LICENSE: All these programs are offered under the GNU GPL license. To my best knowledge, the programs implement the algorithms as intended, but I can not guarantee that they are bug free. I will do my best to correct bugs when they are found but I can not make any guarantees for that either. Anyone who finds bugs is very welcome to let me know and I will try to take care of them as best as I can. Suggestions for extensions/improvements are also very welcome. #################################################### #################################################### REFERENCES: A few benchmark examples are collected at http://www.iet.ntnu.no/~svensson/Benchmarks/index.html. You can compare your results with those presented there. Borish, J., "Extension of the image model to arbitrary polyhedra," J. Acoust. Soc. Am. 75, 1827-1836 [1984]. Svensson, U. P., Fred, R. I., Vanderkooy, J., "An analytic secondary source model of edge diffraction impulse responses," J. Acoust. Soc. Am. 106, pp. 2331-2344 [1999]. Torres, R. R., Svensson, U. P., Kleiner, M., "Computation of edge diffraction for more accurate room acoustics auralization," J. Acoust. Soc. Am. 109, pp. 600-610 [2001]. Svensson, U. P., Calamia, P. T., "Edge-diffraction impulse responses near specular- and shadow-zone boundaries," Acta Acustica/Acustica 92, pp. 501-512 [2006]. Svensson, U. P., Calamia, P. T., Nakanishi, S., "Frequency-domain edge diffraction for finite and infinite edges," Acta Acustica/Acustica 95, pp. 568-572 [2009].