In this paper we present a fast evaluation of the Rayleigh integral, which leads to fast and robust solutions in inverse acoustics. The method commonly used to reconstruct acoustic sources on a plane in space is Planar Nearfield Acoustic Holography (PNAH). Some of the most important recent improvements in PNAH address the alleviation of spatial windowing effects that arise due to the application of a Fast Fourier Transform to a finite spatial measurement grid. Although these improvements have led to an increase in the accuracy of the method, errors such as leakage and edge degradation can not be removed completely. Such errors do not occur when numerical models such as the Boundary Element Method (BEM) are used.Moreover, the forward models involved converge to the exact solution as the number of elements tends to infinity. However, the time and computer memory needed to solve these problems up to an acceptable accuracy is large. We present a fast (O(n log n) per iteration) and memory efficient (O(n)) solution to the planar acoustic problem by exploiting the fact that the transfer matrix associated with a numerical implementation of the Rayleigh integral is Toeplitz. In this paper we will address both the fundamentals of the method and its application in inverse acoustics. Special attention will be paid to comparison between experimental results from PNAH, IBEM and the proposed method

Boer, A. de

Wijnant, Y.H.

Wind, J.W.

English

Wind, Jelmer

Wijnant, Ysbrand H.

de Boer, Andries

University of Twente Research Information

FAST EVALUATION OF THE RAYLEIGH INTEGRAL ANDAPPLICATIONS TO INVERSE ACOUSTICSJ. W. Wind∗1, Y H. Wijnant1 and A. de Boer11Department of Engeneering Technology, University of TwenteP.O. Box 217, 7500 AE Enschede, The Netherlandsj.w.wind@utwente.nlAbstractIn this paper we present a fast evaluation of the Rayleigh integral, which leads to fast androbust solutions in inverse acoustics. The method commonly used to reconstruct acousticsources on a plane in space is Planar Nearfield Acoustic Holography (PNAH). Some of themost important recent improvements in PNAH address the alleviation of spatial windo ingeffects that arise due to the application of a Fast Fourier Transform to a finite spatial mea-surement grid. Although these improvements have led to an increase in the accuracy of themethod, errors such as leakage and edge degradation can not be removd completely. Sucherrors do not occur when numerical models such as the Boundary Element Method (BEM)are used. Moreover, the forward models involved converge to the exact solution as the numberof elements tends to infinity. However, the time and computer memory needed to solve theseproblems up to an acceptable accuracy is large. We present a fast (O(n log n) per iteration)and memory efficient (O(n)) solution to the planar acoustic problem by exploiting the factthat the transfer matrix associated with a numerical implementation of the Rayleigh inte ralis Toeplitz. In this paper we will address both the fundamentals of the method and its applica-tion in inverse acoustics. Special attention will be paid to comparison between exp rimentalresults from PNAH, IBEM and the proposed method.INTRODUCTIONRayleigh’s second integral gives the pressure field for given normalvelocities on an infiniteplane [2]. It can be written as follows (see figure 1(a)):p(~y) = 2iωρ0∫SG(|~y − ~x|)vn(~x)dS with G(r) =e−ikr4πr(1)J. W. Wind, Y. H. Wijnant, A. de Boer~y~xSvn(a) Rayleigh integral (cross sec-tion)fieldsource(b) source and virtual planefieldvirtual plane(c) reconstructionFigure 1:Inverse acoustics using the Rayleigh integralWhere~x and~y denote the source and field points respectively,i, ω ρ0 andk are respectively,the imaginary unit, the frequency in radians per second, the density of the fluid at rest and thewave number= ω/c0, wherec0 is the speed of sound.p, vn andS denote respectively, thefield pressure, the source normal velocity and the set of all points on the infinite plane. IfS ishas finite area, then the equation gives the pressure field for a source inan infinite baffle – anacoustically hard surface that does not vibrate.The Rayleigh integral can be used for inverse acoustics as follows. Given a sourceradiating sound to afinitenumber of sensors, named the field (see figure 1(b)), one can definean infinite virtual plane between the source and the field. One can regard the virtual plane,rather than the physical source, as the radiator of sound. This changeof boundary conditionsdoes not change the measurement but reconstruction to the plane is simpler.We now have thesituation in figure 1(c). If it is assumed that the sensors are located close tothe physical source,then the normal velocities far away from the sensors are negligible. Neglectin these velocitiesis equivalent to placing a finite virtual plane in an infinite baffle. Hence, a finite virtual plane isused in the reconstruction. The accuracy of the numerical model can be increased by extendingthe area of the virtual plane and by decreasing the size of the elements. Likemost numericalmodels, this (forward) model converges to the exact solution if the number of degrees offreedom tends to infinity.Note that in this derivation, the number of sensors is not assumed to be infinite. In NAHtechniques, it is necessary to obtain the Fourier transform of the field pressure as a functionof the location. Naturally, this Fourier transform can not be found exactlybased on a finitenumber of measurements, and an approximation needs to be made by means of wi dowingtechniques. In the current methods, windowing techniques are not necessary and the forwardmodel converges to the exact solution of the Helmholtz equation irrespectiveof th numberof sensors. Finally, it is noted that calculating a rate of convergence forthe inversemodel as afunction of the source area remains difficult but numerical experiments demonstrate that thismodel converges as well [8].ICSV13, July 2-6, 2006, Vienna, AustriaFORWARD MODELIn the derivation below, the physical source is not significant. Therefore, the virtual planewill be referred to as the source for reasons of simplicity. In this section wedemonstratethe structure of a transfer matrix that follows from the a discretized Rayleighintegral. Thisproperty only occurs on regular grids. Hence, we assume that the source grid is equidistantand that a sensor is placed above every node. As demonstrated in a later section, the latterassumption can be weakened. For reasons of simplicity, we consider only constant elementswith a single Gauss point here, but the structure has been proven for a much wider range ofmodels. In this simple case, the model becomes:Hv = p with hij = 2iωρ0SeG(|~yi − ~xj |) (2)whereH, v andp are the transfer matrix and the vectors containing source velocities andfield pressures respectively.Se is the area of an element.~xi and~yi denote the locations of asource Gauss point and a field point respectively.Each element of the transfer matrix depends on the distance|~yi − ~xj | only: this isthe distance between the source point with numberj and the field point with numberi. TheToeplitz property (see equation 3) is revealed when it is noted that the same distances occurin a very specific structure. From figures 2(a) and 2(b), it can be deduced that the matrix isToeplitz. Furthermore, from figure 2(c), it can be seen that the matrix is symmetric.Toeplitz : hi−1,j = hi,j+1symmetric : hi,j = hj,i}H =h11 h12 h13 h14h12 h11 h12 h13h13 h12 h11 h12h14 h13 h12 h11(3)Due to this structure, the transfer matrix is fully determined by its first column. Hece, thetime and memory needed to calculate and store the transfer matrix can be reduced fromO(n2)toO(n). Although this property is not abundantly present in literature, it should benot d thatit is not a new discovery: Toeplitz matrices have been proven to exist for awide range ofconvolution operators on finite grids [1] [8].The matrix structure is similar for a two-dimensional source. A Toeplitz property canthen be identified in both directions. To be exact: when the nodes on the firstrow are numbered~x1 ~x2 ~x3 ~x4~y1 ~y2 ~y3 ~y4(a) occurrences of|~y1 − ~x1|~x1 ~x2 ~x3 ~x4~y1 ~y2 ~y3 ~y4(b) occurrences of|~y2 − ~x1|~x1 ~x2 ~x3 ~x4~y1 ~y2 ~y3 ~y4(c) occurrences of|~y1 − ~x2|Figure 2:Repeated distances for a Rayleigh integral model with constant elements and asingle Gauss pointJ. W. Wind, Y. H. Wijnant, A. de Boer1 to k and those in the second row are numberedk + 1 to 2k etcetera, such that the sourcevelocity distribution is not a matrix but a vector, then the structure of the transfer matrix isblock Toeplitz with Toeplitz blocks (BTTB) [1]. The memory usage and time needed fornumerical integration remainsO(n) (wheren is the total number of nodes).The model discussed here has been chosen for its conceptual simplicity but is notsuited to acoustic nearfield reconstruction. In practice, an adaptive quadrat re is used to calcu-late the numerical integrals accurately. Linear elements are used instead of constant elementsand it is possible to measure the acoustical particle velocity as well as the pressure. Impor-tantly, the source mesh is finer and larger than the field grid. Although these modifications donot change the concept of the method, it is stressed that they are crucialfor its accuracy.The method presented above will be referred to as TRIM (Toeplitz RayleighIntegralMethod).REGULARIZATIONThere is a wide range of efficient regularization methods available for Toeplitz matrices [1].However, if the sensors are placed above only some nodes – i.e. if the source mesh is largerand finer than the field grid – then part of the matrix structure is lost and not all efficientmethods can be applied. If the number of sensors is very large, a fast (O(n log n)) matrix-vector multiplication can be used to make the regularization more efficient. In order to exploitthe speed of this operation we apply the Krylov-subspace methods that arelso popular inIBEM, where the matrix-vector multiplications are performed with the fast algorithm.It is noted that the termfastrefers to the asymptotical speed of this operation. Numeri-cal experiments have indicated no increase in speed when the number of field points does notexceed 1000. In this case, it is equally efficient to apply an ordinary matrix-vector multipli-cation (keeping in mind that only the first column of the transfer matrix needs to be stored incomputer memory).The regularization methods we use are as follows. LSQR regularization is used formost problems. As an alternative regularization method, we solve the least-square form ofTikhonov regularization explicitly using a least-squares solver similar to LSQR.CASE STUDYIn this section, we compare the TRIM results with results from PNAH, IBEM and laservibrometer measurement. The noise source is a hard disk drive suspended in air with wires.The acoustic pressure is measured on an equidistant grid of 17×21 points with a mutualdistance of 10 mm. The field grid is 40 mm from the surface of the object. The reade isreferred to [3] for more detailed information on the measurement and for details concerningthe PNAH technique used.A reconstruction has been made with an implementation of the boundary elementmethod which has been especially developed for accurate solutions in inverse acoustics [6].ICSV13, July 2-6, 2006, Vienna, AustriaThe surface geometry is meshed with linear quadrilateral elements of 5×5 mm, which leadsto a total of 1872 elements. Boundary element calculations are performed for several frequen-cies of interest. The inverse problem is solved with LSQR. No systematic methodis usedfor selecting the regularization parameter. Instead, it is selected using the L-curve, where anattempt was made to select a solution close to the laser data. It is assumed that thissolu ioncan not be improved significantly without adding new information and may therefor serve asreference result.Linear quadrilateral elements of 5× mm are used for TRIM and the source is chosento be extended 0.1m beyond the projection of the measurement grid on all sides. This leadsto a total of 5760 elements. LSQR regularization is used and the regularizationp rameter ischosen in such a way that the differences between IBEM and TRIM are selected such thatthe results look similar. This way, the observed differences are caused by differences in themethod, not by differences in the selected regularization parameter. Finally, t is noted that thePNAH results are selected in the same manner.The results of the three inverse acoustical techniques are very similar on afrequencyof 1075Hz (see figure 3). The PNAH result is slightly closer to the IBEM result than TRIM.(a) Laser (b) IBEM (c) TRIM (d) PNAHFigure 3:Normal velocity (absolute value) at a frequency of 1075Hz(a) Laser (b) IBEM (c) TRIM (d) PNAHFigure 4:Normal velocity (absolute value) at a frequency of 9668HzJ. W. Wind, Y. H. Wijnant, A. de BoerThis effect may be caused by differences in the regularization methods use . Due to the smalldifference, it is concluded that all tested methods work well on this frequency.Significant differences between PNAH and TRIM can be seen on a frequencyof 9668Hz (see figure 4). The oscillatory behavior in PNAH is not caused by under-regularization: it also occurs in over-regularized solutions. Since a point-s urce is presentand the number of sensors per acoustic wavelength is low (i.e. 3) it is plausible that the in-accuracy is a sampling artifact. The results of IBEM and TRIM do not havethis problemand give similar, accurate, results. This is an indication that TRIM has practical value: thereconstruction of TRIM is accurate on this grid, where a finer grid is probably needed forPNAH.The fact that the results of IBEM and and TRIM are so similar raises a different ques-tion. The inverse solver for IBEM has made use of the information that all sound riginatesfrom an object with a known geometry. This information is not available to TRIM. It onlymakes use of the assumption that there the normal velocities are negligible outside the sourcemesh. However, the results are so similar that we may hypothesize that this lackof informa-tion is irrelevant for the current geometry. Theoretical results that quantify this effect lead tonew insight and possibly, more general efficient numerical techniques.Finally, we make a note on the computational cost of TRIM compared to IBEM. Thecomputational times involved in PNAH have not been researched but it is noted that computerresources needed are so small that computational power is not expectedto b a problem forPNAH, even for very large data sets. All numerical experiments have been p rformed on anIntel pentium 4 workstation with a clock speed of 3GHz and 1Gb RAM, running on WindowsXP. For IBEM, numerical integration requires approximately 6 minutes per frequency on thisvery fine mesh. In TRIM, the time needed is between 0.06 and 0.1 seconds, depending on theintegration method. The computational cost of regularization has also been studied. For bothmethods, LSQR regularization is performed for 80 iterations and the time used for IBEM is 2seconds, versus 3.3 seconds for TRIM. The fact that TRIM is slowerthan IBEM is presumablycaused by the larger number of degrees of freedom in the TRIM model. Itis recalled that theasymptotical order of TRIM is lower. The asymptotical order of memory usage is also lowerfor TRIM because only the first column of the transfer matrix is stored. In practice, the amountof memory used to store the transfer matrix is reduced so much that memory usage is oftendominated by the storage of regularization results, which is a problem that can be dealt with.It may be concluded that TRIM has proven to be a fast, robust and accurate method inthe current case study. Although all methods perform well for the frequency of 1050Hz, wehave also presented a case where TRIM and BEM give an accurate soluion and PNAH doesnot. Finally, it is noted that the comparison presented here is must not be seen a a summaryof an exhaustive survey: such a study has yet to take place.NUMERICAL EXAMPLEAn advantage of numerical regularization techniques versus FFT-based techniques is theirgreater flexibility. For example, the current method is not strictly bounded to equidistant grids.ICSV13, July 2-6, 2006, Vienna, Austria(a) Source distribution (b) Field grid (c) reconstructionFigure 5:Numerical example of source identification on a non-equidistant gridAlthough the symmetric Toeplitz structure only exists if a sensor is placed aboveevery sourcenode, it is possible to assign measurement data on only part of this ’virtual’ sensor grid andleave the sensor data on the other points undetermined. This means it is possible to have acoarser grid in one area of the field and a finer grid in another. This is illustrated with thefollowing example (see figure 5). A grid of Gaussian distributions is used assource (seefigure 5(a)) and the field pressures are calculated at a frequency of5500Hz, 20mm above thesurface, and 5% Gaussian noise is added to the calculated pressures. In the field grid, thedistance between the sensors varies spatially (see figure 5(b)). In figure 5(c), a reconstructionis given. It is clearly visible that the result is smoother towards the top left corner of thereconstruction, which is intuitive because there is less information available ther . We donot wish to claim that this variation of smoothness reflects the amount of informati n in any’optimal’ way. We merely state that it varies in a manner which is intuitive.CONCLUSIONSAn important advantage of PNAH is its conceptual simplicity. The main disadvantage is thefact that a Fourier transform of the spatial data must be obtained from a sall number ofsamples. Even with advanced PNAH techniques, the resulting errors are not always negligible.Also, the padding and windowing techniques needed may be considered less elegant from atheoretical point of view. From a practical point of view, PNAH is an established techniquethat has proven to give results engineers need.The forward model of TRIM converges to the exact solution if the element size tends tozero and the source size tends to infinity, which means no windowing techniques are needed.The fact that this method is based firmly in the field of numerical methods leads to severalother advantages. First of all, there is a strong link between these methods and the mathe-matical study of inverse problems. Second, there is a wide range of ’off the shelf’ solutionsavailable for common problems such as the automatic selection of a regularizationp rame-ter [4] [7]. The main disadvantage is the fact that numerical schemes are needed to obtain aphysical model. Although all of these algorithms are well documented [1] [5], they are notJ. W. Wind, Y. H. Wijnant, A. de Boeryet combined into a commercial package.Naturally, the main advantage of IBEM is the fact that it is suited to arbitrary geometriesand arbitrary distributions of sensors. Disadvantages lie in the computational cost nd memoryusage which are often a bottleneck for large structures and high frequencies. In the case study,we have seen that IBEM required approximately 6 minutes for integration andregularizationon a single frequency where TRIM required 4 seconds. This differenc will increase for largerproblems.Finally, it is noted that to the authors knowledge, the current approach can not be gen-eralized to arbitrary geometries, but problems with cylindrical and spherical geometries canpossibly result in transfer matrices that have an even more useful structure.ACKNOWLEDGEMENTSThis work has been performed within the framework of project TWO.6618 of the DutchTechnology Foundation (STW).Measurements were performed at the Dynamics and Control group of the faculty ofmechanical engineering of Eindhoven University of Technology. Rick Sholte is gratefullyacknowledged for his measurements and his supply of PNAH results.References[1] Hansen, P.C.Deconvolution and regularization with Toeplitz matrices, Numerical Algo-rithms, 29:323-378, 2002[2] Fahy, F.Foundations of Engineering Acoustics, Academic Press, 2001.[3] Scholte, R. Roozen, B. and Lopez Arteaga, I.Quantitative and qualitative verificationof pressure and velocity based Planar Near-field Acoustic Holography, In proceedingsof ICSV 13, Vienna, Austria, 2006.[4] Hansen, P.C.Regularization tools: A Matlab package for analys and solution of ill-posedproblems, Numerical Algorithms, 6:1-35,1994[5] Press, W.H and Flannary, B.P.Numerical Recipes in C++: The art of scientific comput-ing. Cambridge University Press, second edition, 2002.[6] Visser, R.A Boundary Element Approach to Acoustic Radiation and Source Identifica-tion, Ph.D. thesis, University of Twente, 2004[7] Engl, H.W. Hanke, M. and Neubauer, A.Regularization of Inverse Problems, Kluwer,1996.[8] Wind, J.W. iTrim software documentation, internal document, University of Twente,2006.

Fast evaluation of the Rayleigh integral and applications to inverse acoustics

NARCIS 

A Boundary Element Approach to Acoustic Radiation and Source Identiﬁcation,

Deconvolution and regularization with Toeplitz matrices,

Foundations of Engineering Acoustics,

iTrim software documentation, internal document,

Numerical Recipes in C++: The art of scientiﬁc computing.

Quantitative and qualitative veriﬁcation of pressure and velocity based Planar Near-ﬁeld Acoustic Holography,

https://ris.utwente.nl/ws/files/6164978/Wind06ICSV13.pdf

Fast evaluation of the Rayleigh integral and applications to inverse acoustics

Abstract

Similar works

Full text

Available Versions

University of Twente Research Information

NARCIS