Computer modeling of barrel-vaulted sanctuary exhibiting flutter echo with comparison to measurements

: Computer modeling in acoustics allows for the prediction of acoustical defects and the evaluation of potential remediations. In this article, computer modeling is applied to the case of a barrel-vaulted sanctuary whose architectural design and construction led to severe ﬂutter echoes along the main aisle, which was later mitigated through acoustical remediations. State-of-the-art geometrical acoustics and wave-based simulations are carried out to analyze the acoustics of this space, with a particular focus on the ﬂutter echoes along the main aisle, before and after remediations. Multi-resolution wavelet and spectrogram analyses are carried out to isolate and characterize ﬂutter echoes within measurements and computer-simulated room impulse responses. Comparisons of simulated responses to measurements are also made in terms of decay times and curves. Simulated room impulse responses from both geometrical acoustics and wave-based methods show evidence of ﬂutter echoes matching measurements, to varying degrees. Time-frequency analyses isolating ﬂutter echoes demonstrate better matches to measurements from wave-based simulated responses, at the cost of longer simulation times than geometrical acoustics simulations. This case study highlights the importance of computer modeling of acoustics in early design phases of architectural planning of worship spaces.


Introduction
Computer modeling has been used in the field of room acoustics to predict and analyze the acoustics of spaces for over fifty years [1,2].Traditionally, computer modeling in acoustics has been based on geometrical acoustics (GA) [3] with its associated ray-tracing and image-source algorithms [1,4] and its various extensions to include diffraction (for an overview of methods, see [5]).GA simulation methods have reached a point where various mature GA-based commercial softwares are available (e.g., [6,7]).However, the limitations of GA are relatively well-known, in that GA is generally only valid for high frequencies, and its extent of validity into low frequencies depends on the room being modeled and on the relative positions of sources and receivers with respect to any obstructions in the space (i.e., the significance of diffraction effects).
Wave-based methods-such as, e.g., finite-element [8], finite-difference [9][10][11][12], and boundaryelement methods [13]-offer another avenue to make use of computers in acoustical modeling [14].Wave-based methods are, in theory, valid for the entire range of audible frequencies and capture all instances of diffraction, but this comes at the cost of higher compute requirements.For example, computer memory storage requirements in wave-based methods generally scale with the cube of the frequency (in three spatial dimensions), and the number of floating-point calculations required can scale even more severely (e.g., to the fourth power of the frequency).As a consequence of these high compute requirements wave-based methods have seen less use over the years than geometrical approaches.On the other hand, wave-based methods tend to be highly parallelizable, and recent progress in parallel implementations [15][16][17][18][19][20], combined with steady increases in compute power from parallel hardware such as graphics processing units (GPUs), are bringing large-scale wave-based simulations to within reach of consumer hardware [21].
This study aims to make use of GA and wave-based acoustic computer modeling in the study of acoustics of the sanctuary of Christ the King (CTK) Charismatic Episcopal Church of New Paltz, NY.Built in 2009, this church has a barrel-vaulted ceiling over the nave portion of the sanctuary, with the floor of the nave carpeted and covered with rows of upholstered chairs, except for a marble-tiled center aisle along the length of the room.As a result of the design choices that left the barrel-vaulted ceiling centered above a rigid aisle, the room exhibited a strong flutter echo along the center aisle of the room.Acoustical remediations were thus called for to mitigate the strong flutter echo.Remediations were carried out in 2018 in the form of panel absorbers mounted on the ceiling above the aisle; see Figure 1.A previous study by the first author [43] details the measurements undertaken and the acoustical remediation carried out, and its effect on decay times in the room.The present study aims to model this space using state-of-the-art geometrical (Odeon v15) and wave-based (finite-difference time-domain [44]) computer simulations, and compare the obtained simulation results to measurements taken before and after remediations.The GA and wave-based frameworks are chosen because they are, arguably, the two room acoustic simulation methodologies receiving the most interest from the research community currently.This study builds upon a previous preliminary study [45] wherein GA simulations were carried out, but evidence of flutter echoes (as seen in measurements) were not found in simulation results.The purpose of this study is thus to (a) refine the GA computer model and (b) to see if a more general wave-based simulation is able to better reproduce the fine features of such non-diffuse effects.The simulation results from the GA and wave-based models are compared in terms of objective parameters (T(20),T(30)) and energy decay curves.Additionally, as the aforementioned flutter echoes are not well represented by traditional room acoustic objective metrics (which assume a diffuse field), a multi-resolution wavelet analysis is employed for the comparison of flutter echoes in measurements and computer simulation results, along with spectrogram time-frequency analyses.Evidence of flutter echoes are found in simulated room impulse responses, with the wave-based results showing significant similarities to measurements in multi-resolution analyses.This case study is significant because it comprises one of the few studies comparing both wave-based and geometrical acoustics simulations to measurements.Additionally, this case study highlights the need for computer modeling in early design phases of building construction.The contributions of this paper are thus: (a) comparison of measurements from a non-diffuse space with state-of-the-art GA and wave-based simulations of the space (including before and after acoustical remediations); (b) the use of multi-resolution wavelet analyses for the detection and analysis of flutter echoes; (c) discussion and analysis of strengths and weaknesses of computer modeling methodologies in regards to this unique space.
The outline of the rest of this paper is as follows.Section 2 details the methods of the study, including measurements, computer simulation methodologies, and signal processing and analyses.Section 3 provides an overview of the results and time-frequency signal analyses thereof, and a discussion of those results and their implications.Finally, Conclusions are found in Section 4.

Methods
The following section includes discussion of the measurement methodology, the modeling methodology, and the signal processing techniques used.

Measurements
Room impulse response (RIR) measurements were made before and after the installation of the acoustical panels using a logarithmic sine-sweep IR technique specified in ISO 3382-1 [46].The room under investigation falls into the category of performance space, therefore ISO 3382-1 was considered appropriate for developing the measurement technique.Of the two different room excitation methods supported by this standard, the integrated impulse response method was used, and the logarithmic sine-sweep was selected as the source sound signal.This was selected because of the improvement of the signal-to-noise ratio and reduction in harmonic distortion which it provides.A comprehensive set of measurements was taken using several different sound source and receiver devices and locations.The following describes only the measurements used for the comparison with modeling presented in this study.

Room Configuration
The room under investigation is a 1546 m 3 sanctuary which contains a barrel-vaulted ceiling over the nave portion of the room.The floor of the nave is completely carpeted and covered with angled rows of upholstered chairs except for a 1.83 m wide marble-tiled center aisle along the length of the room.The barrel-vaulted sanctuary ceiling was constructed with a 10.4 m radius arc with a 11.4 m span and 6.5 m height at the crown.All of the wall surfaces of the sanctuary are built with gypsum board and finished with a standard eggshell paint.
Due to complaints of unintelligibility along the center aisle of the nave, the room configuration underwent acoustical remediation in 2018.This treatment involved the installation of 12 Auralex Acoustics absorbing panels on the ceiling (ten on the barrel vault and two on the flat ceiling above the alter).Each of these panels is 0.05 m × 1.2 m × 1.2 m (2 in × 4 ft × 4 ft).The panels are suspended horizontally from the ceiling, 6 m from the floor, with the center of each panel located 1.5 m from the crown of the barrel-vaulted ceiling.

Measurement Equipment
The acoustical measurements used in this comparison were completed by broadcasting a logarithmic sine sweep into the room using a Brüel and Kjaer TYPE 4292-L dodecahedral loudspeaker, as shown in Figure 2. The input signal was a 4 s logarithmic sine sweep, with a bandwidth of 360 Hz chosen to excite the 64 Hz-8 kHz octave bands.Acoustical pressure measurements were collected during the sine sweep using a PCB model 378B20 omni-directional microphone, with a model 426E01 preamp.All data collection was performed at a sampling frequency of 51,200 Hz.

Source and Receiver Locations
Acoustical measurements were taken of both the original configuration and the acoustically treated room.Acoustical data was collected from a wide range of source and receiver locations, corresponding to possible sound source/listener locations during ordinary use of the space.
While a large number of data points were collected in the space, after reviewing the RIRs a small subset of those points which highlight the unique acoustical behavior of the room were strategically selected for comparison with the computer models.All measurements used were taken with the source located in the center of the room along the aisle.This source location was selected for use in the study because it most effectively excited the flutter echo.Additionally, six receiver locations were selected for comparison with the modeling because these locations highlight specific features of the room's acoustical behavior.This selection included two points located off the aisle, in the audience and four points along the aisle.The location and height of each of the receiver locations are listed on Table 1, and can be seen visually on Figure 3.

Computer Modeling
In this subsection, a brief overview of the computer modeling techniques (GA and wave-based) is given, along with specific parameters used for the simulations themselves.

Geometrical Acoustics (GA) Methods
Geometrical acoustics assumes that sound travels akin to light rays, i.e., sound travels along rays and these sound rays bounce specularly off objects and/or scatters non-specularly or diffusely [47].Sound rays can carry amplitude or energy (i.e., magnitude only), and interactions with walls in terms of sound absorption and/or scattering are typically described by frequency-band coefficients taking on values between zero and one [5,48].Diffraction is not strictly part of the GA framework (by some definitions, GA simply precludes diffraction [49]), but there has been a great deal of work on integrating early edge diffraction into GA simulations, following the idea that diffraction can be described by new sound rays, emanating from edges [50] and applying approximate solutions [49]-or analytic solutions [51] when available.GA is often coupled with statistical models for late decay in the impulse response-i.e., that sound decays evenly throughout the room and can be treated as an ensemble average (i.e., a diffuse field) [52].For computational reasons in GA simulations, it can be beneficial to take advantage of diffuse-field assumptions to reduce the number of specular rays that are needed to trace a full impulse response, and various popular algorithms have evolved to this effect, such as the "Ray Radiosity Method" employed in Odeon [53] and the "Diffuse Rain" algorithm [54] employed in RAVEN [52].It would follow, however, that the more reliance there is on scattering and diffuse late-field assumptions, the harder it can be to reproduce any fine features of non-diffuse reverberation (e.g., flutter echoes).

GA Computer Modeling Methodology
The commercial software Odeon, version 15 [53], was used as a tool for GA computer modeling.A computer-aided design (CAD) model of the room was made using Trimble Sketchup, starting from architectural drawings and using refinements from laser measurements in the space (where architectural specifications were not followed exactly in construction).This CAD model of CTK church was built with a compromise between level of detail and suitability for GA modeling; that is, chairs and light fixtures were removed, as these are typically not included in CAD models used for GA simulation, and the curved ceiling was polygonized into twelve flat rectangular sections (see Figure 4).These rectangular sections were marked as "Fractional" surfaces within Odeon, which tells the program to treat each section as one part of a larger curved surface (as recommended [53]).Seats were replaced with two rectangular "audience" surfaces, to which absorption and scattering were applied in order to approximate the seating area (without persons) as best as possible within the chosen GA software framework.Absorption coefficients for each surface material are given in Table 2. Surface-roughness scattering coefficients for each material within Odeon were set to 0.05 (smooth surfaces) from which Odeon estimated octave-band scattering coefficients [53].These coefficients were chosen based on materials within Odeon's material database.It should be noted that a detailed calibration campaign was not carried out in this case, as it was understood that such calibrations (e.g., [23]) typically rely on diffuse-field objective parameter measurements, whereas this space was known to be highly non-diffuse, in particular along the aisle where the flutter echo was most pronounced.While there will be unavoidable input data uncertainty in the computer model [55], it was reasoned that a calibration to decay times would not have a significant effect on reproducing the non-diffuse effect of interest here.In order to obtain a good GA approximation from the Odeon software, 100k rays were used (per receiver) with the transition order of two, where the transition order specifies at which reflection order the software switches from its image-source algorithm to its ray-based algorithm.The number of rays was chosen in order to go beyond the "Precision" recommendations of Odeon, as it was hoped that the ray-based solution would converge to some best solution with a large number of rays.A detailed comparison of results with number of rays used is not within the scope or interest of this study, but it was found that 100k was a good compromise, and significantly more rays did not impact output results.Angle-dependent absorption (ADA) was used for all materials as it was found to have a significant effect on the authenticity of the room impulse response outputs, in particular where one expected to hear flutter echoes.While this choice was made based on perceptual observations, a perceptual similarity experiment was not within the initial scope of this study.

Wave-Based Methods
Wave-based methods aim to solve the wave equation in three spatial dimensions through numerical grids or meshes.The advantage of using a wave-based method is that the numerical solution naturally comprises all wave phenomena, which includes both geometrical behavior and diffraction.However, such numerical solutions can only be resolved in practice up to a given frequency, due to approximation errors and limitations in computational resources.Traditionally, wave-based methods have only been used to resolve modal behavior below the Schroeder frequency [56], but it is becoming increasingly possible to use wave-based methods beyond Schroeder frequencies of rooms thanks to increases in parallel compute power (see, e.g., [20]).In this study, a finite-difference time-domain (FDTD) method was used [10][11][12], which is a regular grid-based method that makes use of some time-marching scheme, in this case Störmer-Verlet integration [57], and finite-difference stencil operations to approximate spatial derivatives in the wave equation in three spatial dimensions [58].FDTD methods have similarities to simple finite-volume methods [10] and in some cases finite-volume methods can be used as generalizations of FDTD [59] in order to refine numerical boundary conditions to better adapt to non-Cartesian wall-boundaries [60].The FDTD implementation used in this study follows the finite-volume framework laid out in [60][61][62].This framework makes use of a frequency-dependent locally-reactive impedance boundary model, which allows for passive realizations of complex impedances [61].To devise such boundary models, wall absorption data for the wave-based model is also taken from Table 2.However, it should be pointed out that locally-reactive impedances are limited to approximately 0.95 Sabine absorption [47], so numerical clamping was employed for coefficients of acoustical panels whose random-incidence absorption coefficients went beyond the limit of the locally-reactive model.

Wave-Based Computer Modeling Methodology
For the FDTD wave-based model, the Sketchup CAD model of CTK church was modified to include a seating area overtop a carpet surface, since, rather than relying on scattering coefficients to emulate their scattering behavior, the FDTD model can naturally simulate the scattering from seats (see Figure 4).As mentioned above, the FDTD model used followed the framework laid out in [61,62], and a non-Cartesian cubic close-packed (CCP) grid was chosen for efficiency reasons [63], with finite-volume refinements at boundaries.The interior FDTD scheme thus followed the "CCP" scheme featured in [12].Boundary impedance conditions were set from absorption coefficients given in Table 2, where a passive complex impedance [60] was fit through numerical optimization [44].The grid spacing for the FDTD model was set such that all frequencies under 8 kHz would have less than 1% numerical dispersion errors, resulting in a grid spacing of approximately 5.5 mm.This FDTD model, with an overarching bounding box of 1950 m 3 and simulating 5 s of impulse response, was run on a cluster of NVIDIA GPUs hosted at the University of Edinburgh, with compute times of approximately one hour per second of impulse response.
The wave-based results, comprising frequencies below 8 kHz, were complemented with high-frequency ray-traced room impulse responses computed using a high-density ray tracing code developed at the University of Edinburgh.As such, the "wave-based" results are in fact "hybrid" wave/ray-based results, but for simplicity the term "wave-based" will be used henceforth, since the majority of the subsequent analyses focus on frequencies below 8 kHz, and in terms of octave bands, 90% of the hybrid RIRs are generated from the wave-based FDTD approach.Nevertheless, it is worth mentioning that for these GA results, 500k rays were traced throughout the scene with specular reflections only (i.e., a strict GA approximation), using angle-dependent absorption (based on a locally-reactive model), and using the same CAD model as for the FDTD simulations.Linear phase, partition-of-unity finite impulse response (FIR) filters were used to combine FDTD and ray-traced outputs.

Signal Processing of Room Impulse Response
The room impulse response (RIR) of the measured data was calculated for each sine sweep measurement through deconvolution of the recorded signal with the input sine sweep.These RIR waveforms were resampled to both 48 kHz and 45.254 kHz (following [64]), and used as a basis of all additional signal processing procedures as described.

•
Spectrogram Analysis: Spectrograms were generated from 48 kHz sampled RIR data using a 10 ms Hann window with 120 sample overlap.

•
Energy Decay Curve: The energy decay curves were generated from the 48 kHz sampled RIR data using the ITA toolbox.The curves were normalized and plotted in dB using P re f = 20 × 10 −6 over a 40 dB reduction in SPL.

Wavelet Multi-Resolution Analysis
Because of the highly non-diffuse nature of the sound field, it was necessary to look beyond the standard room acoustics metrics such reverberation times(T(20) or T( 30)) and energy decay curves in the evaluation of the sanctuary's acoustical characteristics.To this end, the use of the discrete wavelet transform (DWT) was employed to isolate and identify specific patterns relating to the (non-diffuse) flutter echo.In particular, multi-resolution analysis (MRA) performed using the DWT was utilized for the identification of the periodic patterns in the waveforms that constitute the flutter echo.

Wavelet Selection
Before completing the MRA, the Wavelet type (family) and length had to be selected.The type and length of the wavelet used can have a significant effect on the results of the MRA, therefore, these parameters needed to be selected based on the expected shape of the signal features associated with the flutter echo signal.It was demonstrated by Wiltschko et al. [66], that there is an improvement using wavelets over traditional filtering techniques for isolating spikes or sharp peaks in measured data.Specifically, wavelet transforms are able to isolate frequency bands in the signal with less distortion of the fast peaks in the signal.Their findings are significant for applications to non-diffuse room acoustics, because the identification of patterns in the flutter echo acoustical signals requires a similar process for the isolation of peaks.Zhang et al. [67] applied statistical ANOVA analysis to the evaluate the wavelet family type and length in evaluation of fMRI data.Three different types of orthogonal, compactly supported wavelet families evaluated were evaluated in their study: Daubechies Extremal Phase (dbN), Daubechies Least Asymmetric (symN), and Coiflet (coifletN).Based on their research, the following recommendations were made:

•
The differences between the wavelet types which they applied were insignificant, therefore, any of the tested wavelet families could be applied with similar results.

•
The length of the wavelet had a significant effect on the uniformity of the results, where wavelets with longer support gave rise to more consistent results.
The efficacy of applying these recommendations to non-diffuse room acoustics was confirmed for the flutter echo triplet patterns through a direct comparison of the amplitude and quality of the signal peaks for level 4 MRA of the original measured data at receiver location R04.A single spike in a triplet was the feature extracted from the using each of the selected wavelet support lengths.
After the assessment of the wavelets, it was determined that the MRA of the RIR be performed using the MATLAB Wavelet Toolbox modwt and modwtmra commands with a Daubechies 6 (db6) wavelet (containing 6 vanishing moments and a corresponding wavelet length of 11).The result of using this wavelet family can be seen in the evaluation of the MRA of an impulse response sampled at 45,254 Hz in Figure 5.
Performing a discrete wavelet transform alone, results in phase-shifted responses at each level.To avoid this time delay, the MRA uses the detail coefficients and the low resolution approximation developed by the discrete wavelet scaling function to project the wavelet transform of the signal onto different wavelet subspaces, or levels, relating to decreasing frequency ranges.For the wavelet analysis presented here, to aid the comparison with traditional octave-band filtering, the waveform signals were resampled to 45,254 Hz.Using this sample rate, the frequency ranges of the decomposition levels relate directly to octave band frequencies as seen in Table 3.The MRA levels of the dB6 wavelet can be seen in Figure 5 which shows the application of the wavelet to an impulse approximation (0.0442 ms triangular pulse) generated with a 45,254 Hz sample rate over −0.01 to 0.01 s.
-  The frequency bands of each of the MRA levels can be clearly identified in Figure 6, by comparison of the spectrum from the original signal with the spectrum of each of the MRA levels from Figure 5, demonstrating the MRA's behavior as a linear-phase filter bank with center frequencies corresponding to each of the octave bands.Using this wavelet MRA, provides an octave band filter which retains the specific fine features associated with the periodic peaks of the flutter echo signal.

Results and Discussion
Using the methods described in Sections 2.3 and 2.4, the measured and modeled RIRs were used to develop a set of figures which collectively demonstrate the unique behavior of the flutter echo in the room under investigation.These figures are presented in such a manner that comparisons between measured and modeled data can be readily made.In the initial investigation all six selected receiver locations were evaluated, both with and without panels.However, to reduce the complexity for reporting purposes, the number of receiver locations retained for each of the curves generated was reduced to two or three receiver locations.To this end, while R02 and R06 were used within the study to evaluate the response at these locations, their waveforms were considered redundant for reporting purposes, and were not included in this results section.The receiver locations and room configuration (original or with panels) in each figure was chosen to highlight specific features of the signals which demonstrate characteristics of the room's acoustical behavior.
Rationale for the selection of receiver locations is as follows: • R01 and R04: As both of these figures are the same distance from the sound source, with one located on the aisle and the other located off the aisle, comparison between these two points demonstrates the non-diffuse nature of the room.

•
Original and Panel (with R01 and R04): In order to demonstrate the changes in the room acoustics as a result of the addition of acoustical panels, some figures feature side-by-side comparisons of the original room configuration and the room configured with panels as shown in Figure 4. • R03, R04 and R05: These three receiver locations represent three different heights at the same location in the room as shown in Figure 3.This comparison demonstrates the unique behavior of the flutter echo at different heights relative to the source height.

RIR Spectrogram
Spectrograms were generated for each of the selected receiver points.After evaluation of all points, receiver locations R01 and R4 were selected for comparison.Figures 7 and 8 highlight the difference of the RIR on the aisle and in the audience.Inspection of the spectrograms in each of these figures reveals key distinctions in the late decay (1.5-4.5 s) behavior between the R01 (audience) and R04 (aisle) RIR, especially in the 500-4000 Hz octave bands.These distinctions, which are most pronounced for the measured data and the wave-based model, will be further described in the following evaluation.This discussion of the spectrograms provides an overview of all of the different characteristics which will be evaluated throughout the entire results section.

•
Behavior of the room on and off the aisle: Comparison of Figures 7 and 8, demonstrates the non-diffuse nature of the room.As stated earlier, both receiver points R01 and R04 are the same distance from the sound source, however, because of their location with respect to the barrel-vaulted ceiling, their spectrograms display very different behaviors.
In high frequencies, above 4000 Hz, the measured original R01 and R04 spectrograms exhibit similar behavior, with the majority of the energy having died out by 0.5 s.Below 4000 Hz, the behavior of the two locations are quite different.In the audience, the signal dies out almost completely by 1.5 s, while in the aisle, it remains well beyond 4 s.Additionally, the signal is characterized by a periodic pattern representative of the flutter echo present in the aisle.Furthermore, along the aisle, it is possible to identify the comb filter behavior expected from a flutter echo.

•
Behavior of room before and after the addition of panels: In order to demonstrate the changes in the room acoustics as a result of the addition of acoustical panels, some figures feature side-by-side comparisons of the original room configuration and the room configured with panels as shown in Figure 1.
Similar to what was observed in the comparison of R01 and R04, at higher frequencies (above 4000 Hz), the addition of the panels had very little effect on the spectrogram.Below 4000 Hz, in the audience area, the addition of the panels reduces the energy, such that the signal is essentially eliminated after 1.5 s.On the aisle, the effect of the panels is more pronounced, with the primary effect being to eliminate the majority of the late decay energy in the 250-2000 Hz octave bands after 2 s.The panels also have a secondary effect of reducing the amplitude of the peaks from 0.5 to 1.5 s.
These results indicate that the panels had the desired effect of reducing the RT along the aisle without completely eliminating the spacious feel of the room promoted by the flutter echo.

•
Comparison of the GA and wave-based modeling with the measured data: The spectrograms provide a means of comparison of the similarities and distinctions between the measured data and the two different modeling techniques, GA and wave-based.Comparison between measured and modeled response before 0.5 s highlights specific early behaviors (before 0.5 s) and later behaviors.Before 0.5 s, both models contain more energy above 8 kHz than the measured data, with Odeon predicting much more sustained high frequency energy.It is important to note that above 8 kHz, both Odeon and the wave-based model use ray-tracing techniques.Also, the measured data has a band of higher energy sustained past 0.5 s around 6500 Hz, which is not identified in the wave-based model.
Most significantly, the very distinct pattern of the flutter echo detected in measurements on the aisle in Figure 8, are clearly identified in the wave-based model, but are not present in the spectrogram of the Odeon model RIR.
In order to confirm and refine the results of the spectrograms, the RT, EDC, and fine features calculated from the RIR were observed.Results from each of these acoustic metrics are described in the following sections.

RIR Reverberation Times
Using the method described in Section 3.1, reverberation times (RT) were generated for each of the identified receiver points.Figure 9 shows the reverberation times (T (20) shown in dashed and T (30) in solid) calculated using the measured and modeled RIR data at specific receiver locations: R04, on the aisle; and R01, in the audience (both before and after the addition of the panels).Comparison of the measured RT (blue lines) for the original room locations on and off the aisle displays the non-diffuse nature of the room.Both receiver points, R01 and R04 are at the same height and the same distance from the sound source, however, because of their location with respect to the barrel-vaulted ceiling, the reverberation times differ substantially from one another, especially in the 250-4000 Hz octave bands.The addition of the acoustical panels greatly reduces all of the elevated RTs, and results in much lower discrepancies between the on and off aisle RTs across all of the octave bands.
There is little difference between the measured T(20) and T(30) RTs after the addition of the panels, both on and off the aisle, all of the RTs are below 2 s, and at the majority of frequencies, the T(20) and T(30) RTs are nearly identical.In the original room configuration, however, the T(20) and T (30) RTs show a number of striking dissimilarities.Comparison between the T(20) and T(30) RTs in the audience indicates that there is a strange behavior occurring in the 500 Hz octave band.The T(30) RTs are consistently slightly higher than the T(20) RTs, however, in the 500 Hz octave band, T(30) is almost two seconds greater than the T(20) RT.This discrepancy warrants a closer look at the EDCs from which the RTs are calculated.This substantial difference between T(20) and T (30) RTs is even more expansive on the aisle, where the difference between the two RTs are greater than one second for the 250-1000 Hz octave bands.For frequencies greater than 1000 Hz, however, the RTs become almost identical, showing less discrepancy than these frequencies in the audience measurement.
In comparing the RTs for the modeled RIR, the wave-based model exhibits very different behavior on and off the aisle.The wave-based model does not correlate with the measured data in two aspects: first, it does not anticipate the extremely high T(30)s around 500 Hz, and secondly, it overestimates the RTs above 1000 Hz in the audience, and in the aisle after the addition of the panels.In comparing the predicted effect of adding the panes, the audience receiver location suggests little difference in the RTs with and without the panels, which is not reflective of the measured behavior.On the aisle, however, there is a large difference, with a reduction in RT for almost all frequency bands, and especially in the 250-1000 Hz bands.Although the RTs are slightly higher than the measured results in the 1000-8000 Hz octave bands, the model still does a reasonable job of estimating the reduction in the RTs on the aisle due to the addition of the panels.
The Odeon model shows very little difference between the aisle and audience receiver positions, and greatly overestimates the low frequency RTs, however it provides a close approximation of the audience RTs above 1000 Hz.The Odeon model also indicates very little change between the RTs due to the addition of the panels.

Energy Decay Curves
Prompted by the long reverberation times, and the large discrepancies between T(20) and T(30) reverberation times, energy decay curves (EDC) were evaluated for each of the selected points.As with the reverberation time calculations, receiver locations R01 and R04 were used for comparison as shown in Figure 10.Based on the RT values, it was expected that the audience measurements would display a distinct difference between the slopes of the early and late portions of the EDC.As predicted by the RTs, the EDCs for the measured data in the audience almost always exhibit a linear decay, with the exception of the 500 Hz octave band, which demonstrates a distinct change in slope around 1 s, resulting in a much shallower curve during the late reflections, which is in keeping with the much higher T(30) prediction at that octave band.
As expected from the T(20) and T(30) results presented in Section 3.2, where the measurement T (30) values on the aisle are substantially higher than the T(20) in the 250-2000 Hz octave bands, the measured EDCs in these octave bands are also non-linear.Not only do they display an elbow around 1 s in the 500 and 1000 Hz octave bands, but they also exhibit two other non-linear characteristic behaviors.The first such characteristic is the broad ripples in the EDC, indicating variations in the late decay energy which cause the increase in the RTs in the 250 to 2000 Hz octave bands.The second non-linear characteristic feature of the flutter echo, a stepped decay in the EDC can be seen more clearly on the close-up plot of receiver R04's 2000 Hz octave band EDC in Figure 11b.A close-up of the EDCs shown in Figure 11 shows in greater detail how the stepped decay is overlaid on the broad ripples.This stepped decay behavior has been identified in rooms with strong flutter echo [41], and indicates that there are non-diffuse features of this sound field that need to be evaluated through direct observation of the waveform, rather than simply through the use of traditional room acoustic objective metrics which assume a diffuse field.To this end, a multi-resolution wavelet analysis will be used in the following section for further comparison of flutter echoes in measurements and computer simulation results.The highlighted area of the EDC in Figure 11 indicates the time which will be highlighted in the comparisons of on and off the aisle receiver locations as well as before and after the addition of the panels.

Wavelet Multi-Resolution Analysis
As described in Section 2.4, when applied to the measured RIR using an appropriate wavelet, the MRA of the signal can be used to isolate features that occur at specific frequency bands, making the DWT useful in peak detection and pattern recognition in acoustical signals.Based on the spectrograms shown in Section 3.1, it was identified that the flutter echo was comprised of repeating patterns of peaks which persisted most strongly in the 1 and 2 kHz frequency bands.
To support the observation from the spectrogram that the flutter echo was most distinct in these frequency bands, wavelet MRA scale levels 1 through 7 of R04 measured data were evaluated during mid and (t > 1 s) and late reflections (t > 2 s) for patterns related to the flutter echo in Figure 12.Observation of the MRA confirms that the features of the flutter echo are most pronounced in levels 4 through 6, relating to the 500-2000 Hz octave bands (for data sampled at 45.254 kHz).In particular, based on the distinct peaks of the triplets identified in the level 4 MRA scale, representing the 2000 Hz frequency band, this level was selected to represent the triplet pattern generated by the flutter echo in each of the figures generated using the wavelet MRA.Throughout the rest of this analysis, Level 4 MRA will be used to demonstrate the fine features of the flutter echo and provide comparisons between measured and modeled results for signals on and off the aisle, before and after the addition of acoustical panels, and at various receiver heights.

Comparison of Flutter Echo On and Off the Aisle
In Figure 13, the RIR is plotted at 1 s both on and off the aisle, for the measured and modeled responses.The presence of the flutter echo is much more distinctly identifiable on the aisle (R04) than it is off of the aisle (R01), with the fine details of the flutter echo apparent as a triplet pattern with the middle peak characteristically being the largest amplitude.The dotted lines on the plots for R04 identify the middle peak of the triplets in the measured signal.The locations of these peaks were determined automatically using the MATLAB findpeaks command with a minimum peak prominence of 10 ×10 −7 , and a minimum peak distance of 1600.
Comparison with the modeled results show that the peaks of the triplets in the wave-based model line up closely with the dotted lines that were identified from the measured data.The Odeon model, on the other hand, does not exhibit the clear triplet peaks with the correct timing in the same way that the wave-based model does.As observed in Figure 13, there is a distinct flutter echo pattern in the measured data from the original room configuration.This triplet pattern is also clearly present in the measured data after the addition of the panels.Based on this comparison, it can be seen that the primary effect of the treatment is a substantial reduction of the amplitude of the flutter echo in the RIR, without eliminating the triplet structure.
Comparison of the behavior of both of the models demonstrates a similar result, with each of the models displaying a marked decrease in the amplitude of the signal, while still retaining the predicted structure of the echoes.The fine features of the triplet pattern are more clearly identified in the wave-based model in the original configuration, and therefore, are more clearly identifiable after the addition of the panels as well.

Structure of Triplets Based on Receiver Height
Although the structure of the flutter echo is identifiable throughout the RIR, it can be seen most clearly in the late reflections t > 2 s.Observation of the RIR waveform in Figure 15 demonstrate the changes in the fine structure of the flutter echo at different receiver points.This structure is directly related to the location of the receiver relative to the source in relation to the reflective surfaces.Specifically in this case, because the reflective surfaces are the tiled aisle and the curved ceiling, when the receiver is at nearly the same height as the sound source (at R04 the receiver is 0.2 m below the source), the flutter echo pattern takes on the very distinct triplet pattern as seen in the R04 plots of Figure 15, with the center peak being higher in amplitude than the outer peaks, due to constructive interference from the ceiling and floor reflections.
When the receiver is below the source, the peaks become spread out, as the two wave fronts which caused the center peak in the R04 triplet pattern, are now separated, to form a wider structure.This structure can be seen most clearly in the wave-based model for R03 which is located 0.7 m from the source.In R05, the distinct center peak of the triplet pattern is still visible, but it is slightly wider because its height, 0.3 m above the source is slightly further from the source than R04.Once again, the fine features of the structure of the flutter echo at this height is clearly replicated by the wave-based model, whereas the Odeon model exhibits some of these features, but they are greatly diminished in amplitude, and are not as accurately reproduced.
Results of spectrograms, T(20) and T(30) calculations, energy decay curves, and wavelet multi-resolution analysis have been presented in this section, providing a means to discuss both the unique features of the flutter echo itself, and the efficacy of GA and wave-based modeling to reproduce these features.

Conclusions
This article investigated the acoustics of a barrel-vaulted sanctuary with a significant flutter echo along the aisle, and its later acoustical remediation.Through the use of spectrograms, T(20) and T(30) calculations, energy decay curves, and wavelet multi-resolution analysis, the broad behaviors and the fine details of the acoustical behavior of a non-diffuse room exhibiting a distinct flutter echo have been evaluated.This evaluation included comparison of the behavior on and off the center aisle of the room, before and after the addition of acoustical panels, and identification of similarities between measured data and two different acoustical modeling techniques.Within these analyses, GA and wave-based simulated RIRs were compared to measurements.There were various uncertainties in computer model absorption and scattering coefficients, but this non-diffuse room did not warrant a detailed calibration through diffuse-field objective metrics.Despite this, the simulated results showed agreement with fine features of flutter echoes along aisle, with wave-based results modeling the fine details of the echo more accurately than GA results.However, GA results showed better agreement in the differences between off-aisle decay times and curves before and after remediation, as compared to measurements.The discrepancy in off-aisle decay time modeling could be attributed to the limitations of a locally-reactive boundary condition (angle-dependent absorption) in the modeling of panel absorbers (which themselves deviate from the locally-reactive model) by the wave-based model.The applied absorption in the GA model, however, is a combination of angle-dependent and random-incidence absorption [53], which is perhaps better suited to modeling of panel absorbers.The particular behaviors of both models in regards to off-aisle decay times in this space could be the topic of future work.Future work could also include a perceptual evaluation of measured and simulated responses, to better refine quantitative similarities in reproduction of flutter echoes, along with the use of alternative commercial GA softwares.

Figure 1 .
Figure 1.(left) photograph of Christ the King (CTK) church taken after remediation.The main aisle is in view, along with absorbing panels mounted under the barrel-vaulted ceiling.(right) Simplified computer-aided design (CAD) model of CTK church after remediation with axis units in meters.

Figure 2 .
Figure 2. Brüel and Kjaer TYPE 4292-L dodecahedral loudspeaker in the source location used for measurements.

Figure 3 .
Figure 3. Layout of source (S) and six receivers (R01-R06), from x-y (top) and x-z (bottom) cross-sectional views.Note: R03-R05 are coincident in the top view, and axes units are in meters.

Figure 7 .
Figure 7. Spectrograms for receiver location R01 (in audience area, 3 m from source); before and after the addition of the acoustical panels.

Figure 8 .
Figure 8. Spectrograms for receiver location R04 (on aisle 3 m from source); before and after the addition of the acoustical panels.Note the presence of acoustical signals from 1.5-4.5 s not observed in Figure 7 above.

Figure 9 .
Figure 9. Reverberation times for receiver locations R01 (audience, 3 m from source) and R04 (aisle, 3 m from source); before and after the addition of the acoustical panels.

Figure 10 .
Figure 10.Energy decay curves for receiver location R01 (audience, 3 m from source) and R04 (aisle, 3 m from source); before and after the addition of the acoustical panels.For conciseness, only the 250 Hz through 4000 Hz octave bands are shown.

Figure 11 .
Energy decay curves for measured and modeled data before and after acoustical remediation.(a) Receiver located off the aisle 3 m from source, (b) receiver located on the aisle 3 m from source.Highlighted area reflects the data used for wavelet MRA in Sections 3.4.1 and 3.4.2.

Figure 13 .
Figure 13.Level 4 MRA of original measured and modeled room impulse response (RIR), R01-in the audience, 3 m from source, R04-on the aisle 3 m from source.Dotted lines identify the middle peak of the triplets in the measured signal.

3. 4 . 2 .
Comparison of Flutter Echo Before and After the Addition of PanelsComparison of the change in the flutter echo behavior as a result of the addition of the acoustical panels can be observed in Figure14.This figure compares the response after 1s on the aisle at R04, both in the original configuration and after the addition of the panels.The red lines included in each figure are calculated using an exponential curve fit of the peaks of each of the triplets, and are provided simply to aid in the visualization of how the peak amplitudes of the signal decay over the time window shown.

Figure 14 .
Figure 14.MRA of measured and modeled RIR at R04 located on the aisle 3 m from source.Comparison of the original room and the room with the addition of the acoustical panels.Red lines represent exponential decay of peaks and are included to aid in visualization of the signal decay.

Figure 15 .
Figure 15.Structure of triplet pattern measured for the original room configuration on the aisle for three different receiver heights.R03 is located at a height of 1m, R04 at 1.5m and R05 at 2m. Highlighted area demonstrates the compression of the triplet for receiver R04 (height near source height).

Table 1 .
Receiver locations included in this study (locations specified in m in terms of the coordinates shown in Figure1).