Immersive Wave Propagation Experimentation: Physical Implementation and One-Dimensional Acoustic Results

We present a fundamentally new approach to laboratory acoustic and seismic wave experimentation that enables full immersion of a physical wave propagation experiment within a virtual numerical environment. Using a recent theory of immersive boundary conditions that relies on measurements made on an inner closed surface of sensors, the output of numerous closely spaced sources around the physical domain is continuously varied in time and space. This allows waves to seamlessly propagate back and forth between both domains, without being affected by reflections at the boundaries between both domains, which enables us to virtually expand the size of the physical laboratory and operate at much lower frequencies than previously possible (sonic frequencies as low as 1 kHz). While immersive boundary conditions have been rigorously tested numerically, here we present the first proof of concept for their physical implementation with experimental results from a one-dimensional sound wave tube. These experiments demonstrate the performance and capabilities of immersive boundary conditions in canceling boundary reflections and accounting for long-range interactions with a virtual domain outside the physical experiment. Moreover, we introduce a unique high-performance acquisition, computation, and control system that will enable the real-time implementation of immersive boundary conditions in three dimensions. The system is capable of extrapolating wave fields recorded on 800 simultaneous inputs to 800 simultaneous outputs, through arbitrarily complex virtual background media with an extremely low total system latency of 200 μs. The laboratory allows studying a variety of long-standing problems and poorly understood aspects of wave physics and imaging. Moreover, such real-time immersive experimentation opens up exciting possibilities for the future of laboratory acoustic and seismic experiments and for fields such as active acoustic cloaking and holography.

We present a fundamentally new approach to laboratory acoustic and seismic wave experimentation that enables full immersion of a physical wave propagation experiment within a virtual numerical environment.Using a recent theory of immersive boundary conditions that relies on measurements made on an inner closed surface of sensors, the output of numerous closely spaced sources around the physical domain is continuously varied in time and space.This allows waves to seamlessly propagate back and forth between both domains, without being affected by reflections at the boundaries between both domains, which enables us to virtually expand the size of the physical laboratory and operate at much lower frequencies than previously possible (sonic frequencies as low as 1 kHz).While immersive boundary conditions have been rigorously tested numerically, here we present the first proof of concept for their physical implementation with experimental results from a one-dimensional sound wave tube.These experiments demonstrate the performance and capabilities of immersive boundary conditions in canceling boundary reflections and accounting for long-range interactions with a virtual domain outside the physical experiment.Moreover, we introduce a unique high-performance acquisition, computation, and control system that will enable the real-time implementation of immersive boundary conditions in three dimensions.The system is capable of extrapolating wave fields recorded on 800 simultaneous inputs to 800 simultaneous outputs, through arbitrarily complex virtual background media with an extremely low total system latency of 200 μs.The laboratory allows studying a variety of long-standing problems and poorly understood aspects of wave physics and imaging.Moreover, such real-time immersive experimentation opens up exciting possibilities for the future of laboratory acoustic and seismic experiments and for fields such as active acoustic cloaking and holography.

I. INTRODUCTION
In physical wave propagation experiments, a truncation of the domain on which the propagation is studied is inevitable.In a wave propagation laboratory, this truncation is represented by the physical boundaries of the laboratory (e.g., water tank or solid object being examined).The boundaries cause undesired reflections and backscattering of waves propagating on the interior, which can mask or contaminate signals from the target features of interest.
To mitigate this problem, traditional small scale analogue modeling is done with wavelengths in the millimeter to centimeter range (frequencies of 100s of kilohertz to megahertz), which means that the size of the laboratory domain (e.g., a water tank with a side length of a few meters) is typically several hundred wavelengths and the studied targets are many wavelengths in size [1,2].In doing so, the unwanted reflections are separated in time from the signal of interest.However, these laboratory frequencies are several orders of magnitude higher than those actually used in many practical applications, such as acoustics (1-50 kHz) or seismology (sub-Hz to 1 kHz).Since velocity dispersion and frequency-dependent attenuation are common phenomena in many (Earth) materials, especially in porous media, it is clear that observations at the laboratory "micro" scale do not adequately represent observations made at the "macro" scale, for example, in seismic surveying (e.g., [3][4][5]).Some of these scaling effects could be reduced by acting on the ambient temperature [6][7][8].However, such a strategy seems impractical, particularly for large laboratories, and still requires the medium to be self-similar across multiple magnitudes of scales.Hence, to avoid serious upscaling issues, seismic laboratory experiments ought to be conducted at frequencies much closer to those used in actual field practice.One approach, which is not really practicable (especially in the depth dimension), is to use a far larger laboratory facility.A second approach is to virtually expand the laboratory size by making the enclosing walls effectively transparent.This can be achieved by passively attenuating the wave field at the boundaries using absorptive materials e.g., [9][10][11][12], as, for instance, in an anechoic chamber [13].While such materials can be very effective at suppressing reflections from the boundaries, the size of the laboratory affects the wave propagation in a much more fundamental way and imposes a limit on the maximum signal wavelength that can be usefully employed.When the experimentation domain becomes smaller than the Fresnel volume for finitefrequency transmitted waves, the waves start sensing the size and shape of the laboratory.Hence, the amplitudes of forward scattered waves will be affected, and the resulting (modal) wave propagation cannot be detached from specifics such as the geometry of the laboratory [14].As these limitations are not inherent to (active) radiation boundary conditions, an active control of the wave field on the edge of the domain e.g., [15][16][17][18][19][20][21] appears more appropriate to mitigate the influence of the boundary.Here, we present a fundamentally new laboratory for physical, real-time, acoustic immersive wave propagation experiments using active boundary control.It utilizes the recent idea of exact or immersive boundary conditions [22][23][24] to overcome the aforementioned size-related limitations and enables longwavelength, low-frequency wave propagation experimentation.At the same time, the laboratory permits interactions with scatterers placed in a virtual numerical world outside the laboratory.Moreover, the laboratory will facilitate experiments involving active broadband acoustic wave field cloaking and holography.In this work, we first summarize the concept of immersive boundary conditions (IBCs) in Sec.II and briefly review the theory of IBCs and immersive experimentation (Sec.III) to aid the understanding of the numerical implementation of IBCs on a real-time acquisition, computation, and control system (WaveLab system).Next, the practical requirements for a one-dimensional (1D) laboratory implementation of IBCs are analyzed and the attendant latency is discussed (Sec.IV).We then demonstrate the feasibility of a physical implementation of IBCs by presenting results of the first laboratory implementation of IBCs on one side of a 1D sound wave tube (Sec.V).Finally, we briefly discuss the experimental implementation of IBCs in three dimensions (3D) and introduce the low-latency, high-performance WaveLab system enabling real-time, 3D immersive wave experimentation (Sec.VI).

II. CONCEPT OF IMMERSIVE BOUNDARY CONDITIONS
Inspired by the work of Ting and Miksis [25], van Manen et al. [22] introduced a novel set of radiation boundary conditions, which offer an exact solution for broadband, nonreflecting boundaries in numerical wave propagation simulations.The boundary conditions rely on the injection of a secondary wave field on the edge of the domain that destructively interferes with the primary wave field incident on the boundary, thus canceling any undesired reflections.The two-step approach requires (1) the prediction of the wave field incident on the boundary by extrapolating the wave field from a mathematically closed surface of sensors near the edge of the domain and (2) the injection of the extrapolated wave field on the boundary of a spatially and temporally staggered finite-difference grid [26].These socalled exact or immersive boundary conditions introduced by van Manen et al. [22] not only remove size-related limitations, but also allow the projection of a larger scale onto a smaller truncated domain, by using extrapolation Green's functions that include wave interactions with a larger background medium.In this case, IBCs accurately account for arbitrary-order scattering interactions between the truncated and the background medium.This also includes all wave phenomena emerging from these interactions, such as the correct modal behavior of the model under investigation.For these reasons, IBCs are useful for many numerical applications, including target-oriented modeling [24,27], target-oriented waveform inversion [28], reservoir transplantation [29], and elimination of free-surface-related multiples [30], as well as broadband cloaking and holography [31,32].Vasmel et al. [23] proposed to implement IBCs physically for immersive experimentation in a seismic wave propagation laboratory, where IBCs remove the influence of the laboratory boundary and link the physical wave propagation in the laboratory with a simulated propagation in a virtual background medium.We follow the idea of Vasmel et al. [23] and present here the first physical implementation of IBCs in 1D.We also introduce the required data acquisition, computation, and control system for a physical implementation of IBCs in 3D.

III. THEORY OF IMMERSIVE BOUNDARY CONDITIONS AND IMMERSIVE EXPERIMENTATION
Consider two different domains as shown in Fig. 1: (1) a truncated domain V phy , which coincides with the physical laboratory, and (2) a full domain V full , for which V phy can be viewed as being embedded in an arbitrarily complex, virtual background medium with a transparent boundary between V phy and the background.Both domains have a region in common, in which the material properties are identical.This region is bounded by a mathematically closed recording surface S rec on its inner boundary and by a mathematically closed emitting surface S emt on its outer boundary.The latter coincides with the physical boundary of the laboratory.The goal of IBCs is to find and apply boundary conditions on S emt that immerse V phy into V full by canceling reflections from the boundary of V phy , while accounting for all orders of scattering between both domains.This may be due to sources or reflectors within V phy or conversely reflectors or sources within the virtual background medium from which waves leave and/or enter the inner medium.The boundary conditions required on S emt can be found by differencing the representations for the pressure fields in the full and physical domain, denoted by p full ðx; tÞ and p phy ðx; tÞ, respectively.Following Vasmel [33] and Broggini et al. [24], we use Rayleigh's reciprocity theorem of the convolution type and assume a rigid boundary on S emt (i.e., the particle velocity normal to S emt vanishes) to formulate the difference as p full ðx 0 ; tÞ − p phy ðx 0 ; tÞ ½G p;q phy ðx 0 ; x; tÞ Ã v full;i ðx; tÞn i dS; where the asterisk denotes temporal convolution; G p;q phy ðx 0 ; x; tÞ represents the acoustic pressure impulse response, or Green's function, of the medium at x 0 , due to a point source of volume injection rate volume density at x; v full;i ðx; tÞ denotes the ith component of the particle velocity in V full ; and n i is the ith component of the outwardpointing normal on S emt .Hence, IBCs are attained on S emt by emitting monopole sources with the signature of the particle velocity normal to S emt measured in V full .This requires the prediction of v full;i ðx emt ; tÞn i ahead of time, which can be achieved due to causality by extrapolating the wave fields from S rec to S emt using a Kirchhoff-Helmholtz integral [22,34] v full;i ðx emt ; tÞn i Here, G v;q full;i ðx emt ; x; tÞ and G v;f full;i;m ðx emt ; x; tÞ represent the ith component of the particle velocity impulse response of V full due to a monopole source and an m-directed dipole source, respectively; and n m is the mth component of the outward-pointing normal on S rec .To facilitate the numerical implementation, we discretize Eq. (2) in time [24]: strength of the monopole sources at x emt .All future values vfull;i ðx emt ; k; l > kÞn i are summed to a buffer for use at later times of the computation.This is the reason why l also appears on the left-hand side of Eqs.(3) and ( 4).An example for the extrapolation is given in Table I.The core of the implementation of the IBCs on S emt is formed by Eqs. ( 1) and (4).Figure 2 shows a schematic of the numerical implementation of the integral extrapolation on the real-time data acquisition, computation, and control system introduced and further discussed in Sec.VI.

IV. PHYSICAL IMPLEMENTATION OF IBCS IN 1D: PRACTICAL REQUIREMENTS AND LATENCY
The physical implementation of IBCs in a wave propagation laboratory requires a real-time evaluation of Eq. ( 4) and the emission of the predicted normal particle velocities as the signatures of closely spaced monopole sources on the rigid boundary at S emt .It is evident that the physical installation of IBCs entails a range of strict hardwareand latency-related requirements that need to be fulfilled for a successful experiment.In particular, t IBC , the latency up to the emission of vfull;i ðx emt ; k; lÞn i , has to be smaller than the least physical travel time between S rec and S emt .Many of the requirements prescribed by Eqs. ( 1) and ( 4) can be addressed with a range of filtering and processing steps, which we will discuss in detail in this section.However, as we will see, such corrections come at the expense of an increase in t IBC , which effectively dictates an increase in the distance between S rec and S emt .To experimentally validate the feasibility of a physical implementation of IBCs, they are installed on one side of a 1D acoustic waveguide using a portion of the full-scale WaveLab system consisting of two National Instruments (NI) PXIe-7976R and two  4).(a) The surface integral is replaced by a matrix-vector multiplication of the extrapolation Green's functions with the normal particle velocity and pressure recorded at S rec .Quantities N t , N AI , and N AO correspond to the number of time steps and the numbers of input and output channels (or sensors and sources), respectively.Note that the extrapolation Green's functions are scaled by the separation distance of the sources Δx to account for the integration step size.(b) The resulting column vector with dimensions N AO N t × 1 is reshaped to a matrix with dimensions N AO × N t (black arrow), which is used in the recursive evaluation of the time-convolution integral to yield the normal particle velocity at S emt .At each time step k, the leftmost column of the resulting matrix (green box) is applied as the boundary conditions on S emt .The values for l > k (purple box) are stored for later use.TABLE I. Example of the extrapolation using Eq. ( 4) with N t ¼ 3. Modified after Broggini et al. [24].The subscript full and the dependency of vfull;i ðx emt ; k; lÞ on x emt were omitted for improved readability.
PXIe-7965R field-programmable gate array (FPGA) modules supporting 16 input and 28 output channels (see Fig. 13).The waveguide and the experimental setup are discussed in detail in Appendix A and the full WaveLab system is introduced in Sec.VI.

A. The recording surface: Implicit wave field separation
The IBC method relies on a separation of the wave field in real time, which can be qualitatively explained as follows: Only outward propagating waves need to be extrapolated through the background medium and canceled on S emt .The (back)extrapolation of inward propagating waves would lead to an undesired feedback and increasing amounts of acoustic (and electric) energy in the system.Unfortunately, an explicit separation of the wave field in real time is not straightforward due to causality, because low frequencies are not captured in the data at early times of the experiment.Instead, Eq. ( 4) describes an implicit separation of the ingoing and outgoing wave fields on S emt , which constitutes a crucial aspect of the IBC implementation.However, this step requires measuring the pressure and normal particle velocity of V full on S rec .Alternatively, particle velocity can be derived from two or more parallel grids of pressure sensors, since particle velocity (vector quantity) is proportional to the gradient of the scalar pressure field.We achieve this by rearranging the equation of motion and replacing the spatial derivative of the pressure in the normal direction with an appropriate finite-difference approximation [35][36][37].For two pressure recording surfaces at x rec 1 and x rec 2 with a perpendicular separation distance from x rec of −d=2 and þd=2, respectively, we obtain where ρ is the medium mass density and the integral has been replaced by a cumulative sum and multiplication by the time step Δt s .Equations ( 5) and ( 6) hold as long as d is small compared to the minimum signal wavelength [35].We choose d ¼ 1.8 cm for the 1D experiments in air presented in this work.Considering a propagation velocity of c air ¼ 344 m s −1 and a sampling frequency of the acquisition system of f s ¼ 20 kHz (see Sec. VI), this corresponds to approximately two samples per minimum wavelength at the Nyquist frequency of the acquisition system.Experiments were conducted to further reduce d to improve the spatial sampling.However, this decreased the signal-to-noise ratio of the gradient estimate and reduced the travel time between x rec 1 and x rec 2 to below one sample, which led to inferior results.Additionally, pfull ðx rec ; kÞ can be approximated by linear interpolation of the pressure recordings [35,36]: pfull ðx rec ; kÞ ≈ 0.5½ pfull ðx rec 2 ; kÞ þ pfull ðx rec 1 ; kÞ: ð7Þ We follow an approach similar to the one described by Melbø et al. [38] and determine a frequency-dependent calibration filter αðkÞ that optimizes the wave field separation by implicitly accounting for differences in sensor coupling, inaccuracies in the estimated medium properties and sensor separation distance, and errors introduced due to the pressure interpolation and discretization of the integral in Eq. ( 6).To determine αðkÞ, the right-going pressure is calculated explicitly (i.e., offline) according to ⃗pðx rec ; kÞ ≈ 0.5½ pfull ðx rec ; kÞ þ αðkÞ Ã ρcv full;i ; ðx rec ; kÞ; where the seven-coefficient calibration filter αðkÞ is chosen such that the misfit between pfull ðx rec ; kÞ and ⃗pðx rec ; kÞ in a window around the first right-going arrival is minimized (which simultaneously minimizes the left-going energy in the same window).Figure 3 shows results of the real-time particle velocity estimation and wave field separation in the 1D sound wave tube using a portion of the WaveLab FIG. 3. Demonstration of real-time particle velocity estimation and wave field separation in 1D using a Ricker wavelet with a center frequency of 2 kHz as input.The top panel displays two pressure measurements recorded at x rec 1 and x rec 2 .Pressure and (estimated) particle velocity are extrapolated to x rec by evaluating Eq. ( 4) in real time on the WaveLab system.The scaling by the acoustic impedance of air and convolution with αðkÞ is included in the extrapolation Green's functions to obtain a real-time estimate of the right-going pressure ⃗pðx rec ; kÞ.The left-going pressure ⃖pðx rec ; kÞ is obtained by subtracting ⃗pðx rec ; kÞ from the total pressure field offline.The successful separation of the wave field is underlined by the low residual energy of 0.37% and 1.09% in the left-and right-going arrivals compared to the total field contained in the left and right shaded box, respectively.system and two pressure-field microphones.To avoid additional computational latency, the real-time computations according to Eqs. ( 6)-( 8) are carried out as part of the wave field extrapolation in the IBC computation (refer to Appendix B for details).

B. The emitting surface: Removal of the source transfer function
The physical implementation of IBCs requires matching an arbitrary waveform incident on S emt with an emitted waveform.An optimal match of waveforms requires removal of the impulse responses of all hardware components involved in the measurement, extrapolation, and emission process (e.g., sensors, sources, pulse amplifiers, sensor preamplifiers, etc.).For the 1D sound wave tube experiments, we only consider a frequency-dependent transfer function of the source, while the transfer functions of all other hardware components are assumed to be flat (i.e., their impulse responses are delta functions).Hence, we only determine ĥinv ðkÞ, the inverse impulse response of the source, and correct for it by an offline convolution with the extrapolation Green's functions.The signatures of the remaining hardware components are removed by a simple scaling of the extrapolation Green's functions (see Appendix B for more details).

C. The emitting surface: Deviations from rigid surface reflectivity
The derivation of Eq. ( 1) assumes a rigid boundary at S emt .For a free surface at S emt , a similar expression can be derived, which requires dipole (directed) instead of monopole (omnidirectional) sources.Of course, such rigid or freesurface boundary conditions are theoretical constructs impossible to achieve on the walls of a physical laboratory.Hence, in practice, S emt deviates from a perfectly rigid reflector.In principle, such deviations can be handled by using both monopole and dipole sources at S emt , driven with extrapolated particle velocity and pressure signatures, respectively [32].In a physical experiment, however, such a dual extrapolation implies either a doubling in channel count, also causing a significant increase in computational demand, or a reduction in the physical volume of the laboratory by a factor of 8 (in 3D) to comply with the Nyquist sampling requirements.Moreover, the combined use of monopole and dipole sources is rather impractical (e.g., because of the increased levels of acoustic wave scattering at the transducers and electromagnetic interactions of the additional sources).For these reasons, an IBC implementation with dual extrapolation is currently not feasible.Hence, it is indispensable to maximize the impedance contrast at the emitting boundary and to find a pragmatic approach that corrects for the imperfect boundary.As the emitted arrivals from the extrapolation through the background medium should not contain the signature of the boundary reflection, the correction for the imperfect boundary must be applied only to the extrapolation of the direct wave from S rec to S emt (and not to the waves scattered in the background medium).Hence, such a correction requires a decomposition of the extrapolation Green's functions into Ĝdir ðx emt ; x rec ; kÞ and Ĝscat ðx emt ; x rec ; kÞ, which facilitate wavefield continuation of the direct wave and the waves scattered in the background medium, respectively.We account for the imperfect boundary by convolving Ĝdir with an experimentally determined, frequency-dependent reflection coefficient remt ðkÞ, which is characterized by a filter with seven coefficients that optimally matches the incident waveform to the reflected waveform in a leastsquares sense.The filter is valid over the approximate frequency band from 0.5 to 7 kHz.

D. Latency
As we have seen in Sec.III, the normal component of particle velocity of the actual wave field arriving at the rigid boundary needs to be predicted ahead of time, since it is needed as the signature for the monopole sources emitting the IBCs on S emt .The prediction is enabled by evaluation of Eq. ( 4) with the help of the introduced auxiliary recording surface S rec .The total latency of the extrapolation, i.e., the time between physically recording the wave field on S rec and emitting the required boundary conditions on S emt , has to be smaller than the least physical travel time between the two closed surfaces.However, a multitude of delays and latencies contribute to the overall latency t IBC of the IBC implementation: (a) t A , the latency of the WaveLab system, including computational latency related to the wave field extrapolation; (b) t B , the latency related to all other electronic components; (c) t C , the delays from filters incorporated into the extrapolation Green's functions (see Appendix B for a summary of all included filters); and (d) t D , latency due to optional numerical operations (e.g., additional filtering) on the FPGAs of the WaveLab system.Despite a higher channel count and longer device path of the full WaveLab system compared to the subsystem used in this study, first tests and timing verifications have demonstrated that both systems reliably meet a  Hence, latency values will likely be significantly lower for hardware of higher quality (with flatter frequency responses), as is planned to be used for the 3D acoustic experiments.
V. RESULTS

A. Nonreflecting boundary condition
In a first experiment, we demonstrate the capability of IBCs to obtain an active, nonreflecting boundary at S emt , which effectively renders the boundary transparent and replaces it by a homogeneous half-space for broadband acoustic signals.We choose extrapolation Green's functions representing a homogeneous background domain, i.e., we only extrapolate the direct acoustic wave from S rec to S emt , but taking the latency of the IBC implementation into account.The full processing flow to address the previously discussed practical requirements in the extrapolation Green's functions is shown in Appendix B. The IBCs are placed on the right boundary of the waveguide depicted in Fig. 5. Ricker wavelets with center frequencies between 1 and 5 kHz are emitted by the stimulus speaker on the left.The first five panels of Fig. 6 show the wave fields measured at x rec 1 with active and inactive IBCs, respectively.As can be observed, the reflections from the right boundary are substantially suppressed in all five traces, when the IBC source is active, whereas they are present when it is inactive.Panel 6 shows the ratio of the energy in the reflected wave for the active versus inactive IBC experiments.The panel illustrates a reduction in reflected energy of more than 95% over a frequency range from approximately 0.6-5.6 kHz, or approximately 3.2 octaves.In a second experiment, IBCs are positioned in a waveguide as depicted in the top part of Fig. 7.A Ricker wavelet of 2-kHz center frequency is used as the stimulus signal.To measure the resulting pressure field along the entire  waveguide, the experiment is repeated for various positions of the interrogation microphone.The first time-distance panel of Fig. 7 shows the pressure field with inactive IBCs.The two abrupt changes in tube diameter inside the waveguide cause reflections (reverberations), leading to a series of arrivals that are incident on the right boundary and reflect back into the tube.The center panel shows the pressure field with active IBCs on the right.All reflections from the right boundary are significantly reduced in amplitude as a consequence of the IBC implementation.Note that, for the experiment, the Green's functions were not changed compared to the previous experimentaltering the geometry (or medium parameters) of V phy outside the common region between S rec and S emt does not require an update of the extrapolation Green's functions.

B. Immersive boundary conditions including scatterers in the virtual background medium
IBCs are particularly powerful to realize immersive wave experimentation, where the physical domain is immersed in a larger, virtual domain that contains scatterers.In this case, the extrapolation Green's functions of V full contain, and account for, interactions between both domains.This also includes higher-order interactions due to the recursive character of the IBCs.To demonstrate immersive wave experimentation, we modify the extrapolation Green's functions from the previous example by virtually extending the waveguide by 30 cm (i.e., including a perfect rigid reflector in the extrapolation Green's functions as shown in the top part of Fig. 8).As discussed in Sec.IV C, this requires a decomposition of the extrapolation Green's function into a direct wave and a scattered wave part.Only the direct part is convolved with r emt ðkÞ, the reflectivity of the emitting boundary.We emit a Ricker wavelet of 2-kHz center frequency into the waveguide and record the pressure field along the tube.A reference experiment and the IBC results are shown in the upper and lower panels of Fig. 8.The experiment demonstrates the interactions between the physical and background domain: After approximately 2.7 ms, the Ricker pulse reaches the physical boundary, where its physical reflection is largely suppressed.Instead of reflecting at the physical boundary, the pulse is extrapolated through the background medium, reflects from the virtual reflector, and is emitted back into the physical domain.Because of the implicit wave field separation of IBCs, the left-going arrival at approximately 6 ms is not extrapolated back to the boundary.The second  right-going arrival is again extrapolated, its physical reflection suppressed, and its virtual reflection reemitted at approximately 10.2 ms.As anticipated, the IBCs virtually extend the tube by 30 cm.For an observer (a listener) inside the 60-cm-long physical waveguide, the received wave field is an accurate approximation of the wave field in a 90-cm-long tube.The example can be extended for arbitrarily complex virtual media by altering Ĝscat ðx emt ; x rec ; kÞ without changing the experimental setup.These first 1D results illustrate the ability of IBCs for immersive wave experimentation including higher-order interactions between the physical and the virtual domains.Moreover, the results underline the capability of a physical implementation of IBCs to overcome the discussed sizerelated limitations and scaling issues of conventional wave propagation laboratories, thus enabling long-wavelength, low-frequency immersive wave experimentation.

VI. PHYSICAL IMPLEMENTATION OF IBCS IN 3D
A. Latency and processing power requirements The primary goal of the WaveLab system is to physically implement IBCs on the boundaries of a water tank with side lengths of 2 m.To maximize the usable volume for immersive experimentation (i.e., the volume enclosed by S rec ), the perpendicular separation distance between S rec and S emt needs to be minimized.This imposes very strict requirements on the latency of the extrapolation, as t IBC has to be smaller than the least physical travel time between the recording and emitting surface.Besides these latency requirements, IBCs are computationally demanding, as they are nonlocal in time and space, representing an all-to-all problem: The computation of the required output of each individual source on S emt requires the input data from all sensors on S rec .Hence, the physical real-time implementation of IBCs poses severe latency and processing power requirements on the system that implements the recording, extrapolation, and emission, particularly in the 3D case, where the Nyquist sampling requirements on the closed recording and emitting surfaces necessitate a great number of input and output channels.For that reason, a tailor-made data acquisition, computation, and control system was constructed for the physical implementation of IBC in 3D.

B. The WaveLab system: A low-latency acquisition, computation, and control system
The planning, design, and construction of the lowlatency, high-performance WaveLab system was carried out in close collaboration with NI and Viviota.The latency requirements and the memory-bandwidth-bound nature of the extrapolation dictated a departure from conventional compute engines based on central or graphical processing units and a move towards the use of FPGAs, supported by technologies such as peer-to-peer (P2P) streaming to handle the computational load.During the study of the extrapolation kernels of an IBC implementation with 800 input and 800 output channels at a sample rate of 20 kHz, it was found (not shown in this work) that one compute node, consisting of one NI PXIe-7976R FPGA, is capable of computing the outputs for two sources.Hence, a total of 400 compute nodes is required.Subsequently, a massive parallelization of the system without a central computation unit and a placement of the compute nodes towards the end of the device path was the objective.A high-level overview of the resulting system architecture and a photograph of the WaveLab system are provided in Figs. 9 and 10.The lowlevel architecture as well as the device path are shown in Appendix C. To reduce the amount of data transferred between devices at runtime, the precomputed extrapolation Green's functions are transferred to the compute nodes prior to execution.During experiments, efficient lowlatency data acquisition, distribution, and IBC computation are enabled by more than 500 NI FlexRIO FPGA modules.Data acquisition of 800 analog input (AI) channels is distributed over four AI chassis and enabled by fifty 14-bit, FPGA-enabled digitizer modules that each sample 16 channels at 50 MHz.The input data are interpolated to 16-bit, resampled to 20 kHz, and the samples of all 800 channels are aggregated on the AI data master via  9. High-level view of the WaveLab system architecture.From top to bottom: The input data of 800 channels are collected over four analog input (AI) chassis and aggregated on the AI data master.The full input data set is distributed to six analog output (AO) data masters, which further pass the data set to 29 analog output chassis.Each AO chassis consists of up to 14 compute nodes, each capable of computing the output for two sources.The start of an experiment is triggered and synchronized by a central timing and stimulus chassis (blue dashed lines).MXI-Express (MXIe)-a bus-extension technology based on the cabled PCI-Express specification.Six identical copies of the full input data are distributed to six analog output (AO) data masters using novel 10 Gbit s −1 highspeed serial instruments with fiber-optic interconnections.Each AO data master further distributes the sensor data to five (four for AO master 6) AO chassis via MXIe.In each of the 29 AO chassis, up to 14 PXIe-7976R FPGA modules perform the demanding matrix-vector multiplication at each time step of the experiment to compute the required boundary conditions for a total of 800 sources.Finally, each AO chassis consists of a FPGA-enabled signal generator module that resamples the output signals from 20 kHz to 1 MHz and drives out the 16-bit digital signal for up to 28 channels simultaneously.During resampling, a reconstruction filter, or anti-imaging, is applied, which delays the output by 50 μs in order to linearly interpolate between samples.This avoids the emission of stair-cased signals, reducing the generation of spurious high frequencies during the IBC emission.To minimize computational overheads, increase throughput, and hence effectively minimize system latency, P2P streaming and pipelining are employed.P2P streaming allows the FPGA devices in the system to transmit data directly over the MXIe bus to one another, completely bypassing the host processor and its memory controller(s).Prior to the start of an experiment, the required P2P connections between FPGA devices are established.At runtime, the targets can directly read or write data without further communication with, or intervention from, the host.Pipelining makes use of the parallelization capabilities of the WaveLab system to optimize its performance.Instead of sequentially executing (1) data acquisition, (2) computation, and (3) data output, all three steps run in parallel, enabling an increase in the system clock rate and data throughput.This leads to a pipeline depth of three, requiring that the initial two samples of each experiment have zero wave amplitude by default.The full WaveLab system with 800 input and 800 output channels is operating at 12.8 TFLOPs in real time and has successfully passed extensive system verification tests.Moreover, a timing verification proved that the system reliably meets a latency of 200 μs from analog input to analog output, including 50 μs for anti-imaging.Each experiment is triggered from the central timing and stimulus chassis and is 12.5 ms in duration (i.e., 250 samples).Before considering additional latencies (such as the filter delays discussed in Sec.IV), the system thus enables immersive wave propagation experimentation in an approximately 2.7-m 3 physical underwater laboratory inside a virtual acoustic medium with a minimum volume of approximately 3400 m 3 at frequencies as low as 1 kHz and up to 10 kHz, or equivalently a virtual medium with a diameter of 12.5 to 125 wavelengths.However, the use of the WaveLab system is, in principle, not restricted to a specific shape of the experimentation domain or experiments in water.As we have shown in this study, it can also be used for immersive experiments in air-filled waveguides (or any acoustic medium that allows the implementation of the recording and receiving surfaces a sufficient distance apart and instrumented according to the Nyquist sampling requirements).Moreover, the system opens up possibilities for a large range of acoustic wave field cloaking and holography experiments by placing the recording and emitting surfaces in the interior of the experimentation domain instead of at the boundary.

VII. DISCUSSION AND CONCLUSIONS
Immersive boundary conditions have been used in a wide range of applications involving numerical simulations of wave propagation.In this study, we demonstrated the feasibility of a physical implementation of IBCs to actively control the wave field on the boundary of a 1D wave propagation laboratory in real time.This allowed us to virtually remove the imprint of the laboratory boundary on the wave field and to embed, or immerse, the physical laboratory in a larger, virtual background medium, while maintaining all wave interactions between the physical and the virtual medium.We physically implemented IBCs using a portion of a recently developed, massively parallelized, FPGA-enabled, high-performance data acquisition, computation, and control system.The unique system architecture with more than 500 FPGA modules realizes a real-time computing performance of 12.8 TFLOPs, enabling the system to gather the input data from 800 simultaneous input channels and compute the boundary conditions for 800 simultaneous output channels with a verified system latency of less than 200 μs.This system will also be the key enabler for a physical implementation of IBCs and immersive wave experimentation in 3D.In a first experiment, we demonstrated the use of IBCs to obtain nonreflecting boundaries, with an observed reduction in reflected energy above 95% over more than three octaves.Residual energy in the suppressed reflections from the emitting boundary is attributed to imperfect removal of the impulse response of the IBC source, the frequency-independent correction for attenuation between S rec and S emt , which overestimates losses at low frequencies and underestimates them at high frequencies, and to errors resulting from limited spatial and temporal sampling, which degrade the results of pressure interpolation and particle velocity estimation.The sampling limitations also explain the increase of residual energy with increasing center frequency, for which the effective number of samples per wavelength decreases.Additionally, the increased residual relative energy below 0.6 kHz and above 5.6 kHz is likely related to lower signal-to-noise ratios in these frequency regions, because the used input signals contain negligible energy for those frequencies.Despite these minor shortcomings, the IBC implementation proved robust, very repeatable, and effective in suppressing broadband boundary reflections.The obtained results are comparable to reported absorption coefficients for passive absorptive materials e.g., [10][11][12]39,40] and conventional active boundary control [17][18][19][20].Such approaches are, however, still affected by the Fresnel zone volume of the transmitted wave (passive materials) and do not address immersive experimentation (passive materials and conventional active boundary control).In a second experiment, we illustrated the immersive abilities of IBCs by virtually extending the physical waveguide.In this case, the IBCs provide the boundary conditions to virtually remove the physical boundary, while accounting for long-range interactions of the wave field between the physical and the virtual background domain.For a listener inside the IBC-enabled waveguide, the received pressure field is an accurate approximation of the pressure field in a physically extended waveguide.These first 1D experimental results pave the way for novel seismic laboratory experimentation at frequencies much lower than previously possible (as low as 1 kHz), for virtual acoustic immersion, and for applications involving active acoustic cloaking and holography.Practical findings include solutions addressing an imperfect emitting boundary, the removal of the impulse responses associated with the recording and emission hardware, as well as the real-time estimation of particle velocity from pressure recordings.Furthermore, we demonstrated that all proposed solutions can be implemented by an offline manipulation of the extrapolation Green's functions without causing additional computational latency.Most of the proposed steps are adoptable also for higher dimensions and will be of great benefit for the 2D and 3D implementation of IBCs, where frequency-and angle-dependent corrections need to be considered.Such experiments will require the introduced full data acquisition, computation, and control system and are currently under way.

ACKNOWLEDGMENTS
The first author is grateful to Patrick Elison for many fruitful discussions.Furthermore, we would like to thank National Instruments and Viviota for delivery and maintenance of a cutting-edge piece of engineering.This work was funded by the European Union's Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Grant No. 641943.

APPENDIX A: 1D EXPERIMENTAL SETUP IN A SOUND WAVE TUBE
A 1D acoustic waveguide (Fig. 11) was used for the immersive wave propagation experiments demonstrated in this study.The waveguide consists of an air-filled, circular tube of variable diameter.For signals below the cutoff frequency of the planar fundamental mode, the wave propagation in the tube can be regarded as one-dimensional due to the exponential decay of nonplanar modes [41].Assuming that mass density and propagation velocity are constant along the tube, impedance contrasts only depend on the cross-sectional area A. When passing from a tube segment of area A 1 to a segment of area A 2 , the reflection coefficient is given by R To create a transient pressure field inside the tube, a loudspeaker is mounted on one end.A second loudspeaker on the opposite end (at x emt ) applies the IBCs.The pressure field is measured using three pressure field microphones that are mounted flush with the tube wall.Two microphones measure at x rec 1 and x rec 2 at distances of −0.9 and þ0.9 cm from x rec .These microphones provide the necessary pressure recordings to predict the wave field at x emt .A third microphone is separated from the IBC system and is used to interrogate the wave field in V phy .The 1D character of the wave propagation in the waveguide was validated, and closely matching experimentally derived and theoretical reflection coefficients, as well as high signal-to-noise ratios, were observed in previous studies [43].Data acquisition, IBC computation, and IBC emission for the 1D experiments were enabled using the portion of the full WaveLab system shown in Fig. 13.

APPENDIX B: MANIPULATION OF THE EXTRAPOLATION GREEN'S FUNCTIONS
In Sec.IV we saw that a number of processing steps and corrections are required to physically implement IBCs.Additional computational latency can be avoided by including all required filters and operations in the extrapolation Green's functions, instead of explicitly processing the wave field recordings on the FPGAs at runtime.This is possible due to the associative property of all required mathematical operations, including the convolution with the extrapolation Green's functions in Eq. ( 4). Figure 12 illustrates the processing flow for the manipulation of the extrapolation Green's functions, which can be expressed by the following formula: ĜðkÞ ¼ a ĥinv ðkÞ Ã bðkÞ Ã δðk þ t IBC Þ Ã αðkÞ Ã ½r emt ðkÞ Ã Gv;q dir;i ðkÞ þ Gv;q scat;i ðkÞ þ a ĥinv ðkÞ Ã bðkÞ Ã δðk þ t IBC Þ Ã ½r emt ðkÞ Ã Gv;f dir;i;m ðkÞ þ Gv;f scat;i;m ðkÞ; ðB1Þ where particle velocity estimation and pressure interpolation have been included in the extrapolation Green's functions (denoted by a tilde) and b represents a 20-to 6000-Hz bandpass filter to stabilize the numerical integration in the presence of low-frequency noise and to remove high-frequency noise.Scalar a accounts for frequencyindependent hardware signatures (amplifier gains, microphone sensitivities) and for amplitude attenuation between S rec and S emt .In principle, viscothermal losses in waveguides are frequency dependent e.g., [44][45][46].However, we found that including a frequency-independent attenuation coefficient of 3.7 dB m −1 in the extrapolation Green's functions leads to satisfactory results for the IBC implementation across the examined frequency range.

APPENDIX C: LOW-LEVEL ARCHITECTURE OF THE WAVELAB SYSTEM
The low-level architecture of the data acquisition, computation, and control system is shown in Fig. 13.FIG.12. Demonstration of the manipulation of the Green's functions for the extrapolation from two pressure recordings at x rec 1 and x rec 2 to one source at x emt .During the extrapolation, the extrapolated values for the recordings on channels AI 1 and AI 2 are added, which concludes the particle velocity estimation and pressure interpolation.The chosen extrapolation Green's functions only act as an example and are not related to the experiments presented in this study.

FIG. 2 .
FIG.2.Schematic of the numerical implementation of Eq. (4).(a) The surface integral is replaced by a matrix-vector multiplication of the extrapolation Green's functions with the normal particle velocity and pressure recorded at S rec .Quantities N t , N AI , and N AO correspond to the number of time steps and the numbers of input and output channels (or sensors and sources), respectively.Note that the extrapolation Green's functions are scaled by the separation distance of the sources Δx to account for the integration step size.(b) The resulting column vector with dimensions N AO N t × 1 is reshaped to a matrix with dimensions N AO × N t (black arrow), which is used in the recursive evaluation of the time-convolution integral to yield the normal particle velocity at S emt .At each time step k, the leftmost column of the resulting matrix (green box) is applied as the boundary conditions on S emt .The values for l > k (purple box) are stored for later use.

FIG. 4 .
FIG. 4. Estimated time delay of the direct wave at five locations along the waveguide (dots).The line of best fit (solid line) intersects with the ordinate at approximately 1.3 ms, which represents t IBC .The standard error of regression is 0.01 ms or 0.75%.

6 FIG. 6 .
FIG. 6. Panels 1-5 show pressure recordings at x rec 1 for emitted Ricker wavelets with center frequencies from 1 to 5 kHz.Results with an active IBC source (red dashed trace) are compared to the reference with an inactive IBC source (blue solid trace).The first reflections from the right tube boundary are contained in the grey boxes.Numbers represent the energy (sum of the squared pressure values) in the window relative to the reference trace.Panel 6 depicts the energy in the window relative to the reference as a function of frequency, summed over all five experiments.

FIG. 7 .
FIG. 7. Top: Geometry of the three-layer waveguide.Timedistance panel (1): Reference pressure field resulting from the emission of a Ricker wavelet with 2-kHz center frequency with inactive IBCs.Panel (2): As above, but with an active IBC source on the right.Panel (3): Difference of both pressure fields.Positions immediately next to the diameter changes and recording locations do not contain data due to practical limitations.

FIG. 8 .
FIG. 8. Top: Geometry of the homogeneous waveguide.A perfectly rigid virtual reflector 30 cm away from the physical termination of the waveguide (orange dashed line) is included in the extrapolation Green's functions.Time-distance panel (1): Reference pressure field obtained by emitting a Ricker wavelet with 2-kHz center frequency with inactive IBCs.Panel (2): Wave field with an active IBC source.Positions immediately next to the recording locations do not contain data due to practical limitations.

FIG. 10 .
FIG. 10.Photograph of the full-scale acquisition, computation, and control system.

FIG. 13 .
FIG. 13.Device path from analog input to analog output (panel 1) and Low-level view of the WaveLab system architecture (panel 2).Green and blue arrows represent data transfer and trigger signals, respectively.The numbers in parentheses and in the device path correspond to the National Instruments module names.The hardware components and device path used for the 1D experiments demonstrated in this study are denoted by red dashed lines.
At each time step k, the pressure and normal particle velocity fields for which l ≥ k are extrapolated from S rec to S emt and vfull;i ðx emt ; k; l ¼ kÞn i is applied as the source emt ; x; l − kÞ pfull ðx; kÞn m dS: ð4Þ of t A < 0.2 ms from analog input to analog output (not shown here).This is facilitated by the unique architecture and high degree of parallelization of the full WaveLab system described in Sec.VI.Since the isolation and determination of t B , t C , and t D are not straightforward, we estimate the sum of the three latencies by measuring t IBC and subtracting t A .We use a zero-delay extrapolation Green's function (i.e., a delta pulse in the first sample) and convolve it with all previously introduced filters.A 1-kHz Ricker wavelet is driven to a single analog input channel of the WaveLab subsystem.The corresponding output is fed to the source amplifier and further to the stimulus loudspeaker mounted on the side of the waveguide.The time difference of the first arrival measured at five locations along the tube with respect to the input to the WaveLab system is shown in Fig.4.The intersection of the best-fit line with the ordinate at approximately 1.3 ms represents t IBC .Considering that t ).The line of best fit (solid line) intersects with the ordinate at approximately 1.3 ms, which represents t IBC .The standard error of regression is 0.01 ms or 0.75%.latencyD¼ 0, since no optional operations were carried out on the FPGAs, this yields t B þ t C ≈ 1.1 ms.A summary of the different latencies is given in TableII.The overall latency of 1.3 ms corresponds to a distance of 0.45 m in air, requiring a minimum separation between S rec and S emt of the same distance.However, it must be noted that the values of t B and t C are a consequence of the specific hardware used in the 1D experiments.In fact, a similar experiment without removal of the source impulse response led to values of t IBC ≈ 0.35 ms, corresponding to a distance of 12 cm (not shown here).

TABLE II .
Latencies of the 1D IBC implementation.