Estimating the radii of air bubbles in water using passive acoustic monitoring
EEstimating the radii of air bubbles in waterusing passive acoustic monitoring
P. Hubert and L. Padovese Technology and Data Science Department, EAESP - FGVAcoustics and Environment Lab - EP - USP Acoustics and Environment Lab - EP - USP a The study of the acoustic emission of underwater gas bubbles is a subject of boththeoretical and applied interest, since it finds an important application in the devel-opment of acoustic monitoring tools for detection and quantification of underwatergas leakages. An underlying physical model is essential in the study of such emis-sions, but is not enough: also some statistical procedure must be applied in orderto deal with all uncertainties (including those caused by background noise). In thispaper we take a probabilistic (Bayesian) methodology which is well known in thestatistical signal analysis communitiy, and apply it to the problem of estimatingthe radii of air bubbles in water. We introduce the bubblegram , a feature extractiontechnique graphically similar to the traditional spectrogram but tailored to respondonly to pulse structures that correspond to a given physical model. We investi-gate the performance of the bubblegram and our model in general using laboratorygenerated data. © [http://dx.doi.org(DOI number)][XYZ] Pages: 1–11 I. INTRODUCTION
Underwater gas leakages are a matter of con-cern for several industries. Detecting and quantify-ing such leakages is thus a problem with theoreticalinterest but also with important engineering impli-cations.In this paper we are particularly interested inthe use of Passive Acoustic Monitoring (PAM) tech-nologies. This approach is cheaper compared toother, more complex options involving active mon-itoring, but presents the challenge of analyzing theacquired signal which will be contaminated by back-ground noise and other acoustic events at a widefrequency range.In this work we consider the single sensor caseand propose a probabilistic (Bayesian) approach tothe problem of modelling the acoustic behavior ofunderwater air bubbles and bubble plumes. Our ap-proach is independent of physical models; we adoptthe model of (Strasberg, 1956) as our first choice a The authors gratefully acknowledge University of S˜ao Pauloand support from SHELL Brazil (subsidiary company ofRoyal Dutch Shell) and FAPESP, through the Research Cen-tre for Gas Innovation (RCGI) hosted by the University of S˜aoPaulo (FAPESP Grant Proc. 2014/50279-4). We would alsolike to thank FAPESP and CNPq for their support, by grantsnumber FAPESP 2016/02175-0 and CNPq 303992/2017-4. in this paper, but generalization for other models isstraightforward.Based on this approach we propose a graphicalevaluation tool, the bubblegram , which bears resem-blance to the usual spectrogram graphs widely usedin signal processing. However, instead of represent-ing the signal’s energy at different frequency rangesas the spectrogram does, the bubblegram representsthe posterior density for the presence of bubbles ofany radius R at a given time t .The paper is organized as follows: section II de-scribes the Bayesian methodology for the analysisof exponentially decaying sinusoids in noisy signals.Secion III describe the physical model of (Strasberg,1956) and shows how to plug this model into theprobabilistic framework. Section IV presents thebubblegram and section V presents the results ofour model in laboratory data. Section VI concludesthe paper. II. PROBABILISTIC MODELLING
Consider a real-valued signal along a single timedimension, y ( t ), sampled (not necessarily uniformly)at points { t i } . We adopt the notation y t for thediscretely sampled signal. J. Acoust. Soc. Am. / 4 February 2021 Estimating radii of air bubbles in water 1 a r X i v : . [ s t a t . A P ] F e b efine a pulse as an exponentially decaying si-nusoidal function x t = A · e − λt · cos ( ωt + φ ) (1)where A is the (assumed constant) amplitude, λ is the decay factor, ω is the frequency and φ thepulse phase.This equation assumes that t = 0 is the startingtime of the pulse.Instead of writing the pulse equation dependingnon-linearly on the phase φ , we adopt the parame-terization x t = e − λ ( t − t ) ( A · cos ( ω ( t − t )) + B · sin ( ω ( t − t ))) · t>t (2)which is equivalent to the former but replacesthe nonlinear parameter φ by a new, linear parame-ter (the amplitude B ). Also we included the param-eter t to generalize the model for pulses starting attimes different from t = 0, and the function U isthe indicator function, taking the value 1 when U istrue, and 0 otherwise.The acquired signal y t is assumed to be formedby a single pulse starting at time t , corruptedwith Gaussian white noise r t , with E ( r t ) = 0 and var ( r t ) = σ . y t = x t + r t (3)In this general form the model can be used todescribe many different phenomena. Its main in-terest lies in the fact that it explicitly models thepulse in terms of decay constant and fundamentalfrequency, as opposed to traditional spectral analy-sis that focus on the frequency alone. This model,being probabilistic in nature, naturally incorporatesuncertainty into the analysis, and also allows the in-clusion of any prior information available about thesignal and the noise.In this work we are particularly interested inthe situation where the decay constant λ and thefundamental frequency ω are related by a known,deterministic function. This is the case when thepulse is taken to represent the acoustic emission ofa gas bubble in water, in which case both λ and ω depend on the bubble’s mean (equilibrium) radius.Directly estimating this radius can aid the designof signal detectors for underwater gas leakages, andalso serves as a tool for estimation of the leaked gasflow. Now taking the structural model of equation3, and assuming Gaussian white noise, the log-likelihood of the model can be written as (cid:96) ( y t | A, B, ω, σ, λ, t ) = − N log (cid:0) πσ (cid:1) − σ N (cid:88) i =1 ( y t − x t ) (4)Since our main interest lies in the estimationof λ , ω and t , we proceed to marginalize all otherfactors from the likelihood. By adopting uniformpriors for the amplitudes A and B , and a Jeffreys’prior for σ , this marginalization can be performedanalitically, resulting in the following expression forthe marginalized likelihood: P ( y t | t , ω, λ ) ∝ (cid:32) N (cid:88) i =1 d i (cid:33) − / (cid:32) N (cid:88) i =1 e i (cid:33) − / × N (cid:88) i =1 y i − (cid:16)(cid:80) Ni = t y i d i (cid:17) (cid:80) Ni = t d i − (cid:16)(cid:80) Ni = t y i e i (cid:17) (cid:80) Ni = t e i (5)with d i = e − λ ( t − t ) cos ( ω ( t i − t )) e i = e − λ ( t − t ) sin ( ω ( t i − t )) (6)This expression can further on be written in aconcise manner as P ( y t | t , ω, λ ) ∝ K ( ω, λ, t ) (cid:34) (cid:107) y (cid:107) − (cid:18) (cid:104) y, d (cid:105)(cid:107) d (cid:107) (cid:19) − (cid:18) (cid:104) y, e (cid:105)(cid:107) e (cid:107) (cid:19) (cid:35) (7)Notice that this expression is independent of theform of both d i and e i , given that these functions donot depend on the marginalized parameters.This methodology is not new, as it has been orig-inally proposed and studied by (Bretthorst, 1988 , , b , c) and applied to the analysis of NuclearMagnetic Ressonance data. Here we apply it to theanalysis of the acoustic emission of air bubbles inwater. III. ACOUSTIC EMISSION OF AIR BUBBLES
The study of the acoustic emission of air bub-bles in liquid media dates back to (Minnaert, 1933),is treated again by (Strasberg, 1956), and more re-cently has been extensively studied by (Ainslie andLeighton, 2011; Leighton et al. , 1996 , n (Strasberg, 1956) a formula relating the bub-ble’s fundamental frequency of oscilation to its ra-dius is given as: ω = (cid:115) γP ρ (cid:18) πR (cid:19) (8)here R is the bubble steady-state radius; ρ isthe liquid’s density; γ is the ratio of specific heats,and P the static pressure.The model for the pressure variations on the sur-face of a single bubble over time is of the same formas 2, given by: x t = ( ω R ) ρ r α e − ω λ ( t − t ) t>t cos ( ω ( t − t ))(9)with α the initial wall amplitude of the bub-ble, and r the distance between the bubble and thesensor.By assuming constant the liquid’s density, thedistance between the bubble and the sensor, and theinitial amplitude, we rewrite equation 9 as x t = C · ( ω R ) e − ω λ ( t − t ) t>t cos ( ω ( t − t ))(10)Now including a phase parameter and reparam-eterizing, we see that model 10 can be put in theform of 5 with d i = ( w R ) e − ω λ ( t − t ) cos ( ω ( t − t )) e i = ( w R ) e − ω λ ( t − t ) sin ( ω ( t − t )) (11)Using equation 8 that gives ω as a function of R , ω = Ω( R ), we can eliminate ω from the model: P ( y t | t , R , λ ) ∝ K ( t , R , λ ) (cid:34) (cid:107) y (cid:107) − (cid:18) (cid:104) y, d (cid:105)(cid:107) d (cid:107) (cid:19) − (cid:18) (cid:104) y, e (cid:105)(cid:107) e (cid:107) (cid:19) (cid:35) (12)with d i = (Ω( R ) R ) e − Ω( R ) λ ( t − t ) cos (Ω( R )( t − t )) e i = (Ω( R ) R ) e − Ω( R ) λ ( t − t ) sin (Ω( R )( t − t ))(13)Finally, (Strasberg, 1956) also provides an ex-pression for the decay constant as a function of thefundamental frequency. In our notation λ = 0 .
014 + 1 . · − Ω( R ) (14)By plugging this equation into our model we areable to eliminate one more parameter, writing λ = Λ( R ), and arriving at a final model depending onlyon the bubble’s radius R and its time of occurrence, t . This elimination depends crucially on the phys-ical models relating the frequency ω and decay fac-tor λ to the bubble’s radius. In this formulationwe accepted these models as groundtruth, and didnot include any uncertainty in these functions. It ispossible, however, to augmentate the model by in-cluding error terms in the equations for ω and λ ,and also to include prior distributions for all physi-cal constants involved in the model. For this work,however, we chose to adopt the simplest formulationto allow a direct comparison of our model with thetraditional tools of signal analysis. IV. ESTIMATING THE BUBBLE RADIUS
An immediate use for this model is the estima-tion of R and t from a given signal. This can beperformed by first adopting prior distributions foreach parameter in the model, and obtaining a jointposterior by application of Bayes’ theorem.To illustrate this process we generate randomdata for a single pulse representing a bubble start-ing at t = 0 . R = 1 mm . We addGaussian white noise with σ = 0 . ×
500 points in the time × radius plane appear as the color scale in figure 1,where brighter colors are associated with high poste-rior mass for a bubble with a given radius at a giventime. The true values of t and R are marked witha white cross.In figure 1 we also included a spectrogram of thesame data for comparison.The first and most important difference betweenthe visualization of the log-posterior and the spec-trogram is the noise. Since the log-posterior modelincorporates the noise term explicitly in the calcu-lations (via the marginalization of the noise power σ ) it is able to better recover the signal parameters.The spectrogram in itself is not a statistical proce-dure, and thus does not account for the presence ofnoise. The consequence is that signal and noise getblurred, as we see in the figure.It will be also interesting to analyze a signal con-taining a sequence of bubbles of various sizes. Thiswill occurr in most practical situations with highergas flows. The original model we proposed in sectionII assumes a single bubble; however, the presence ofthe parameter t allows the application of this samemodel for a longer signal. One can think of the usual J. Acoust. Soc. Am. / 4 February 2021 Estimating radii of air bubbles in water 3 a) Log-posterior (b) SpectrogramFIG. 1. Log-posterior and spectrogram of simulated, single pulse data.
Short-Term Fourier Transform approach, where a lo-cal model is applied to small sections of a long signal,possibly after applying a smoothing function in thetime domain. The smoothing, in our model, wouldbe represented by a prior distribution on t .What we show in the next figures is the plotof the posterior log-density for R , conditional on afixed value for t . In this sense we are looking forthe radii with larger evidence of appearing in a pulsestarting at time t , for t ∈ [ t a , t b ].When the gas flow is large enough and the bub-bles start to cluster together, the hydrodinamicalmodels for the acoustic emission of a single bubbleare no longer exactly valid (Ainslie and Leighton,2011). However, since the acoustic emission of a gasbubble in water depends only on the forming pro-cess of the bubble, spurious acoustic emissions willoccur only in the events of fragmentation of a bub-ble or the coalescence of two or more bubbles. Forthese reasons we believe that our model can be ap-plied succesfully even to the analysis of higher flowbubble plumes.To understand the behavior of the log-posteriorin this case we simulate a sequence of 10 pulsesoccurring along a one second interval. The bub-bles’ starting times and radii are simulated fromuniform distributions between [0 ,
1] and [0 . , . ×
500 points and R ∈ [0 . ,
2] and t ∈ [0 ,
1] appear in figure 2, where we includedempty circles centered at the true values of each bub- ble. Again the spectrogram estimated from the samesignal is included for comparison.Both figures 1 and 2 show that the model be-haves as expected. In the multiple bubbles case itappears that larger bubbles might be easier to de-tect (since the log-posterior is higher around thesebubbles), but local peaks in the log-likelihood can befound near the radius and time of occurrence of allpulses. By comparing the log-posterior plots withthe corresponding spectrograms we see that the log-posterior not only directly estimates the bubble’sradii but also is more precise, since the model ex-plicitly incorporates the background noise as alreadyobserved.These plots that evaluate the joint posteriorof R and t over a (not necessarily uniform) gridshow some visual resemblance with the traditionalspectrogram plots usually applied to the analysis ofacoustic signals. Inspired by this similarity we callthe plots of the log-posterior bubblegrams , since theyare designed to respond to the acoustic emission ofair bubbles in water, representing the posterior prob-ability associated with a given radius at a given time.One advantage of the bubblegram over the tradi-tional spectrogram is that it does not suffer fromthe usual resolution limitations of the spectrogram;to understand how that can be the case, one canimagine that by the use of the probabilistic modelone is effectively introducing an infinite zero-padding into the signal, thus allowing any arbitrary level ofresolution both in time and radius (frequency) do-mains. a) Log-posterior (b) SpectrogramFIG. 2. Log-posterior and spectrogram for simulated data. V. EXPERIMENTAL RESULTS
In this section we apply our model to the anal-ysis of experimental data. The experiments wereconducted by the
Acoustics and Environment Lab-oratory (LACMAM in the portuguese acronym) atthe Polytechnic School, University of S˜ao Paulo.We used sensors (hydrophones) developed bythe laboratory itself; these sensors have a frequencyband of 5 Hz to 60 kHz , with sensitivity of − ± dB rel µP a (preamplified) or − ± dB rel µP a (not preamplified). A picture of one sensor appearsin figure 3.Two experiments were conducted. The first wasdesigned to measure the size of single bubbles by us-ing a high speed CCD camera with a 25 mm lenswith F . cm (length) × cm (width) × cm (height) was used. Three nozzleswith diameters 2 .
5, 4 and 6 mm were positioned oneat a time under a 22 cm water column, and a flowof less than 1 l/min of air was induced by the useof a flow controller connected to a compressed aircylinder. A picture of the setup is seen on figure 4.The second experiment was designed to obtainacoustic signals emitted by bubble plumes in water.It was performed in a diving tank with a depth of 5 m . The bubbles are created by the injection of com-pressed air into the tank, with flow, pressure and exitdiameter orifice controlled. Figures 5 and 6 show theequipment used in this experiment. In figure 5 theblack lines show the position of the four sensors used, and the red circle indicates the approximate positionof the flow controller nozzle. A. Data acquisition and experimental results
1. Single bubble measurements
We recorded high speed (1000 frames per sec-ond) videos of the bubbles generated in the aquar-ium, with a ruler attached to the aquarium wall.Figure 7 below shows snapshots of bubbles gener-ated with the three different nozzle diameters. Thelabels in each picture indicate the diameter of thenozzle used in the air injection.The acoustic emissions were recorded with asampling frequency of 48 kHz . To create the bubble-grams we first passed the signals through a Butter-worth bandpass filter (with pass range 200 to 3000 Hz ). We then selected the sections of the signalcorresponding to the photographed bubbles. Thesepreprocessing steps were taken in order to acceleratefurther processing of the signal.The bubblegrams obtained from the prepro-cessed signals appear in figure 8. J. Acoust. Soc. Am. / 4 February 2021 Estimating radii of air bubbles in water 5
IG. 3. Hidrophone developed by LACMAM.FIG. 4. Setup for measuring bubble sizes.FIG. 5. Diving tank; see text for details.
IG. 6. Schematic drawing of the experimental setup.
J. Acoust. Soc. Am. / 4 February 2021 Estimating radii of air bubbles in water 7
IG. 7. Images of single bubbles (a)2 . mm nozzle (b) 4 mm nozzle(c) 6 mm nozzleFIG. 8. Bubblegrams for a single bubble with multiplenozzle diameters The Maximum Posterior (MAP) estimates forthe bubbles radii were 0 . . . mm for the 2 .
5, 4 and 6 mm nozzles respectively. Bycomparing this estimates to the pictures in figures 8we find that the estimates are accurate, consideringthe deformation of the bubbles in each snapshot.
2. Bubble plumes
To obtain the acoustic emission of bubbleplumes at higher flow rates the diving tank was used.The nozzle of the pneumatic system was positionedat the bottom of the tank, and several runs wereconducted with different nozzle diameters and gasflow rates. For each run, we open the flow controller nd let the gas escape for two minutes. Then weclose the controller and wait for one minute to openit again.In this section we are interested in studying boththe spectrogram and bubblegrams for different val-ues of the nozzle diameter and gas flow rates. We an-alyze signals obtained with diameter D ∈ { . , , } mm and flows of F ∈ { , , } l min, ie, we analyze9 different sections of the experimental signal.To conduct our analysis we extract small sec-tions (2 seconds) of the signals and start by obtain-ing the spectrograms of each section. These spectro-grams appear in figure 9. (a) d = 2 . mm ; f = 2 l/min (b) d = 6 mm ; f = 2 l/min (c) d = 11 mm ; f = 2 l/mind = 2 . mm ; f = 5 l/min d = 6 mm ; f = 5 l/min d = 11 mm ; f = 5 l/mind = 2 . mm ; f = 10 l/min d = 6 mm ; f = 10 l/min d = 11 mm ; f = 10 l/min FIG. 9. Full spectrograms for bubble plume data.
The spectrograms show that the signal’s energyis concentrated around smaller frequencies (between 250 and 1000 Hz ). However, for the reasons alreadyexposed, the graphs are noisy and do not allow easyor direct estimation of the probable sizes of the bub-bles.The bubblegrams of figure 10, on the other hand,are much clearer and easier to inspect. As we cansee, higher flow rates are associated with bigger bub-bles. (a) d = 2 . mm ; f = 2 l/min (b) d = 6 mm ; f = 2 l/min (c) d = 11 mm ; f = 2 l/mind = 2 . mm ; f = 5 l/min d = 6 mm ; f = 5 l/min d = 11 mm ; f = 5 l/mind = 2 . mm ; f = 10 l/min d = 6 mm ; f = 10 l/min d = 11 mm ; f = 10 l/min FIG. 10. Bubblegrams for bubble plume data.
The bubblegram is a useful representation of thesignal, directly obtaining a probability distributionover bubble sizes. From the figures we can observethat higher gas flow rates create bubble plumes withmore diverse sized bubbles (we see that from the re-gions with relevant evidence, colored blue, but far
J. Acoust. Soc. Am. / 4 February 2021 Estimating radii of air bubbles in water 9 rom the more concentrated regions that show thetypical bubble size in each case). Also the time be-tween the formation of each bubble changes accord-ing to the gas flow.The gas flow appears to have a more evidenteffect on the bubble size distribution than the noz-zle diameter, as we can see by comparing figures inthe same row (i.e., signals obtained with the sameflow) and then comparing figures in the same col-umn (signals obtained with the same nozzle diame-ter). However both gas flow and nozzle diameter arepositively correlated with bubble sizes, even thoughthe flow rate has a stronger effect.By aggregating and normalizing the log-posterior for the bubble sizes we arrive at the ap-proximate density for the distribution of bubble sizesin the entire signal. In figure 11 we plot these dis-tributions for a fixed flow rate, and varying nozzlediameter. d = 2 . mm d = 6 mmd = 11 mm FIG. 11. Log-posterior density for bubble sizes, fixedflow
With these figures we are able to observe moredirectly the effect of both nozzle diameter and flowrate in the estimated bubble size distribution.Our results indicate that this relationship iscomplex. For instance, from figure 11 we observefirst of all that the smallest nozzle induces more ho-mogeneous bubbles populations, since the curves areclose to being unimodal (specially for the lowest flow,blue line). For this smallest flow, the mode radius isalso the smallest, around 0 . mm ; the two biggerflows have very similar curves.Increasing the nozzle diameter we observe thatthe bubble populations start to break down into sub-populations with different sizes (that is, the curvesstart to show large deviations from unimodality).The curve for the smallest flow (blue lines in all threeplots) has always a prominent peak around 0 . mm ,but as the nozzle size increases bigger bubbles canalso be found according to our analysis. The in-termediate and higher flows that we tested behaveless monotonically; for a flow of 5 l/min , when thediameter is set at 6 mm (figure 11, red line) thereappears a high peak around a radius of 1 . mm ,indicating that increasing the nozzle allows the for-mation of bigger bubbles. However, when the diam-eter increases to 11 mm , a subpopulation of bubbleswith radii around 0 . mm appears, along with an-other subpopulation of bubbles with raddi around1 . mm . Therefore, for the bubble plumes withintermediate flow, the bubbles increase in size withnozzle diameter, but for the biggest diameter a pop-ulation of small bubbles is also formed.The curves obtained from the analysis of thesamples with biggest flow (green line in all threesubplots) all have an important peak around 1 mm .Increasing the nozzle diameter from 2 . mm causes the appearance of subpopulations of smallerbubbles (peaks around 0 .
25 and 0 . mm , green linein figure 11). When the nozzle increases to 11 mm (figure 11) the peak around smaller bubbles is muchattenuated, and is replaced by peaks around 1 . . mm . VI. CONCLUSION
We applied a known methodology of probabilis-tic (Bayesian) signal processing to the analysis ofthe acoustic emission of air bubbles in water. Ourmodel directly expresses the acoustic signal in thetime domain as a function of the bubble radii andtime of occurrence, thus providing a useful techniqueto analyze acoustic signals containing the emissionof bubble plumes.Testing our model on laboratory generated datawe found evidence that the model is accurate, pro-
10 J. Acoust. Soc. Am. / 4 February 2021 Estimating radii of air bubbles in water iding good estimates of the size of a single bubblefrom its acoustic emission.To analyze the acoustic emission of bubbleplumes we propose to evaluate the posterior of thesingle bubble model in small temporal segments ofthe signal in the time domain. This is done in a sim-ilar way as the usual Short Time Fourier Transform,and can be interpreted as imposing an informativeprior over the time of occurrence of a single bub-ble. This evaluation, performed over a grid of valuesfor both the time of occurrence and bubble radius,generates a visual representation we called bubble-gram , that is also similar to the usual spectrogram,but is built on a full probabilistic setup that embod-ies a physical model relating the bubble size to theparameters of an exponentially decaying sinusoidal.As such, the bubblegram provides a representationof the signal that is more accurate than the spectro-gram, and can be helpful in the study of the bubblesizes in bubble plumes by passive acoustic methods.
Ainslie, M. A., and Leighton, T. G. ( ). “Review of scat-tering and extinction cross-sections, damping factors, andresonance frequencies of a spherical gas bubble.,” The Jour-nal of the Acoustical Society of America , 3184–208.Bretthorst, L. ( ). Bayesian Spectrum Analysis and Pa-rameter Estimation , 1a ed. (Springer-Verlag). Bretthorst, L. ( a). “Bayesian analysis. i. parameter esti-mation using quadrature nmr models,” Journal of MagneticResonance , 533–551.Bretthorst, L. ( b). “Bayesian analysis. ii. signal detec-tion and model selection,” Journal of Magnetic Resonance , 552–570.Bretthorst, L. ( c). “Bayesian analysis. iii. applicationsto nmr signal detection, model selection and parameter es-timation,” Journal of Magnetic Resonance , 571–595.Leighton, T., Phelps, A., Ramble, D., and Sharpe, D. ( ).“Comparison of the abilities of eight acoustic techniques todetect and size a single bubble,” Ultrasonics (6), 661– 667, , doi: https://doi.org/10.1016/0041-624X(96)00053-4 .Leighton, T. G., Baik, K., and Jiang, J. ( ). “The use ofacoustic inversion to estimate the bubble size distribution inpipelines,” Proceedings of the Royal Society A: Mathemat-ical, Physical and Engineering Sciences (2145), 2461–2484, https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2012.0053 , doi: .Minnaert, M. ( ). “On musical air-bubbles and the soundsof running water,” The London, Edinburgh, and DublinPhilosophical Magazine and Journal of Science (104),235–248, https://doi.org/10.1080/14786443309462277 ,doi: .Strasberg, M. ( ). “Gas bubbles as sources of sound inliquids,” The Journal of the Acoustical Society of America (1), 20–26, https://doi.org/10.1121/1.1908212 , doi: ..