The effect of hardware-computed travel-time on localization accuracy in the inversion of experimental (acoustic) waveform data

TL;DR

This paper investigates how hardware-based travel-time computation methods affect the accuracy and robustness of waveform data inversion in acoustic tomography, especially under noisy and incomplete measurements.

Contribution

It introduces hardware-level implementations for integrated and thresholded travel-time computation and compares their effects on inversion robustness using experimental data.

Findings

ITT improves robustness with noisy data

Hardware preprocessing affects inversion accuracy

Thresholded approach enhances stability

Abstract

This study aims to advance hardware-level computations for travel-time tomography applications in which the wavelength is close to the diameter of the information that has to be recovered. Such can be the case, for example, in the imaging applications of (1) biomedical physics, (2) astro-geophysics and (3) civil engineering. Our aim is to shed light on the effect of that preprocessing the digital waveform signal has on the inversion results and to find computational solutions that guarantee robust inversion when there are incomplete and/or noisy measurements. We describe a hardware-level implementation for integrated and thresholded travel-time computation (ITT and TTT). We compare the ITT and TTT approaches in inversion analysis with experimental acoustic travel-time data recorded using a ring geometry for the transmission and measurement points. The results obtained suggest that ITT is essential for maintaining the robustness of the inversion with imperfect signal digitization and sparsity. In order to ensure the relevance of the results, the specifications of the test setup were related to those of applications (1)-(3).

Tables

Table 1. TABLE I: Transmitters and receivers used in recording the experimental data.

Type	Item	Details
Transmitter	Two-way active speaker	20 Hz – 18 kHz frequency response.
Receiver	Dynamic microphone	Unidirectional 60 Hz – 13 kHz frequency response, -72 dB sensitivity.

Table 2. TABLE II: The diameter and center coordinates of the domain Ω Ω \Omega and the foam cylinders A – C .

Target	x (cm)	y (cm)	Diameter (cm)
A	-11.0	12.0	15.0
B	12.0	10.0	10.0
C	-2.5	-13.0	10.0
$\Omega$	0.0	0.0	59.0

Table 3. TABLE III: Synthesis results for the travel-time calculation modules implemented with Catapult C. The percentages given for area is the percentage of available logic elements used on the FPGA chip (EP2C35F672C6N) of the DE2 board.

Mode	Bits	Speed	Latency	Area
		(MHz)	(Cycles)	(Logic Elements)
Thresholding	8	106.0	13 594	821 (2.5 %)
Thresholding	16	99.4	13 558	1 075 (3.2 %)
Integrating	8	74.8	13 565	891 (2.7 %)
Integrating	16	82.6	13 581	1 253 (3.8 %)

Table 4. TABLE IV: The significance of the acoustic tomography results based on the parameters of other tomography applications. SSSB is a small solar system body, such as a comet or an asteroid. ( a superscript 𝑎 \mbox{}^{a} Signal frequency, b superscript 𝑏 \mbox{}^{b} Signal wavelength, c superscript 𝑐 \mbox{}^{c} Diameter of the target object, d superscript 𝑑 \mbox{}^{d} Diameter of the target of reconstruction)

Application	$\lambda$ (m)${}^b$	$D$ (m)${}^c$	$d$ (m)${}^d$	$D/\lambda$	$d/\lambda$	Significance
Test	$59\cdot {10}^{-3}$	0.58	0.10	10	1.0
MCT (T1)	$19\cdot {10}^{-3}$	0.15	0.02	7.9	1.1	Strong
MCT (T1a)	$19\cdot {10}^{-3}$	0.15	$5.0\cdot {10}^{-3}$	2.0	0.28	Medium
UCT (T1)	$1.5\cdot {10}^{-3}$	0.15	0.02	100	13	Medium
UCT (T1a)	$1.5\cdot {10}^{-3}$	0.15	$5.0\cdot {10}^{-3}$	25	3.3	Strong
SSSBs	$15$	$150$	$15$	10	1.0	Strong
Concrete	$35\cdot {10}^{-3}$	0.3	0.03	10	0.86	Strong
Ionosphere	1.5	$1.0\cdot {10}^6$	$0.1\cdot {10}^6$	$670$	100	Weak

Equations

$$\frac{1}{{\mathtt{c}}^2}\frac{{\partial }^{2}u}{\partial t^2}-{\Delta }_{\vec{x}}u$$

$$u(0,\vec{x})$$

$$\frac{1}{{\mathtt{c}}^2}\frac{\partial u}{\partial t}-\nabla \cdot \vec{g}$$

$$\frac{\partial \vec{g}}{\partial t}-\nabla u$$

$$u_m(t,\vec{x})=u(t,\vec{x})+\epsilon (t,\vec{x}).$$

$$T(\vec{x})=\frac{{\int }_{t_1}^{t_2}t{u}_m(t,\vec{x}{)}^{2}dt}{{\int }_{t_1}^{t_2}{u}_m(t,\vec{x}{)}^{2}dt}.$$

$$T={\int }_{\mathcal{C}}\mathtt{n}ds,$$

$$\mathbf{y}=\mathbf{L}\mathbf{x}+{\mathbf{y}}_0,$$

$${\mathbf{x}}^{(k)}={\left({\mathbf{L}}^T\mathbf{L}+\alpha \mathbf{D}{\mathbf{\Gamma }}^{(k)}\mathbf{D}\right)}^{-1}{\mathbf{L}}^T(\mathbf{y}-{\mathbf{y}}_0),$$

$$D_{i,j}=\beta {\delta }_{i,j}+\frac{{\int }_{P_i\cap P_j}\left(2{\delta }_{i,j}-1\right)ds}{{\max }_{i,j}\left({\int }_{P_i\cap P_j}ds\right)}$$

$$f(t)=\sin (36t)\exp [-35(t-0.60{)}^2]fort=[0,1.2](milliseconds),$$

$$=$$

$${ROA}_{min}$$

$${\hat{\mathbf{x}}}_{\ell +1}=({\hat{\mathbf{L}}}^T\hat{\mathbf{L}}+\alpha {\hat{\mathbf{\Gamma }}}_{\ell }{)}^{-1}{\hat{\mathbf{L}}}^T\mathbf{y},{\hat{\mathbf{\Gamma }}}_{\ell }=diag(|{\hat{\mathbf{x}}}_{\ell }|{)}^{-1},{\hat{\mathbf{\Gamma }}}_0=\mathbf{I}.$$

$$\frac{\partial H(\hat{\bf x},\hat{\bf z})}{\partial\hat{z}_{j}}\Big{|}_{\hat{\bf z}^{\prime}}=-{\alpha}\frac{\hat{x}^{2}_{j}}{(\hat{z}^{\prime}_{j})^{2}}+1=0,\quad\hbox{i.e.}\quad\hat{z}^{\prime}_{j}=|\hat{x}_{j}|\sqrt{\alpha}.$$

$$\hat{\Psi }(\hat{\mathbf{x}})$$

Full text

I Introduction

This paper concerns waveform tomography in which either an acoustic or electromagnetic wave travels through a target object and, based on the external measurements of the wave, the task is to estimate the distribution(s) of a given parameter within the target [1] such as the velocity of the wave or the absorbtion parameter. Tomographic imaging based on wave propagation requires computationally heavy mathematical inversion of the data in order to retrieve the relevant result set, such as an image or three-dimensional model of the test subject.

In this paper, we explore the problem of determining the travel-time of the signal [2, 3] and its relation to accuracy of the the inverse localization. In particular, we focus on the effect of hardware-level computations in applications where the signal wavelength may be expected to be close to the diameter of the details that are to be recovered. Such can be the case, for example, in biomedical microwave or ultrasonic tomography [4, 5, 6, 7, 8, 9, 10, 11], in astro-geophysical subsurface imaging [12, 13, 14, 15, 16] and in non-destructive testing of civil engineering materials and structures [17, 18, 19, 20]. Our aim is determine what effect the preprocessing of digital waveform signal has on the inversion results and in that context to find computational solutions that guarantee robust inversion even with incomplete and noisy measurements. We use travel-time data, as it is known to yield robust information of the unknown parameter [21], and also because it requires minimal data transfer between the hardware and the inversion routines. Data preprocessing performed on embedded hardware has its limitations, and therefore, we aim to find out how the hardware-level evaluation of the travel-time is related to the tomography results. To perform the preprocessing operations, a field programmable gate array (FPGA) [22] chip on a development board is used as a platform for the design. In order to make the implementation fast and flexible, the hardware on the FPGA was implemented by adapting Matlab scripts [23] to C code [24] and then utilizing Catapult C hardware synthesis [25] to generate the hardware from the adapted C code.

FPGA-based processing of tomographic travel-time data is utilized in all above-mentioned application fields. As specific examples we consider (1) the microwave and ultrasonic computed tomography (MCT and UCT) of the breast [4, 5, 6, 7, 8, 9], (2) tomography of small solar system bodies (SSSBs), e.g., the COmet Nucleus Sounding Experiment by Radiowave Transmission (CONSERT) [12, 13, 15, 16], and (3) the acoustic/electromagnetic imaging and testing of concrete structures [17, 18, 19, 20, 10, 11]. An M/UCT breast scan (Figure 1) can be performed by utilizing a 2D sensor ring [8, 4, 26] which records data slices sharing the direction of the plane normal with the ring. MCT and UCT have recently been shown to have the potential to detect and classify breast lesionsat least as reliably as other imaging methods, such as computed X-ray tomography (CT) and magnetic resonance imaging (MRI) [8, 27, 28]. Other recently studied methods for breast cancer diagnosis include, e.g., optical tomography [29, 30]. The CONSERT experiment took place as a part of the European Space Agency’s Rosetta mission. The objective of CONSERT was to recover the internal structure of the nucleus of the comet 67P/Churyumov-Gerasimenko based on sparse lander-to-orbiter signal transfer between the Rosetta spacecraft and a single comet lander Philae (Figure 1). Space technology also involves also the challenging space environment as a limitation [31, 32], leading, e.g., to sparse measurements. In concrete testing [17, 18, 19, 20], the task can be for example to detect interior defects within a concrete element (Figure 1) in a similar way to, e.g. the way voids are localized in geoscientific applications.

In the numerical experiments, we compared integrated and thresholded travel-time (ITT and TTT) approaches via inversion analysis utilizing experimental acoustic waveform data. A 16-bit and 8-bit analog-to-digital (A/D) conversions were tested together with two different threshold criteria and normalization levels. The results obtained suggest that ITT is essential for maintaining the robustness of the inversion if the A/D conversion is incomplete. In order to ensure the relevance of the results, the specifications of the test setup were related to those of applications (1)–(3).

This paper is organized as follows. Section 2 describes the materials and methods, including the mathematical forward modeling and inversion techniques, the test setup, the equipment and the FPGA implementation. Section 3 presents the inversion results. Section 4 sums up the findings of this work and discusses the results and the potential direction of future work.

II Materials and methods

II-A Forward modeling

II-A1 Signal wave

In this study, a waveform signal is modeled as a scalar field $u$ representing an electromagnetic or acoustic wave. The computational domain $\Omega$ is assumed to contain the target object of the tomography together with its immediate surroundings. During the measurements, that is, when $t\in[t_{1},t_{2}]$, the transmitters and receivers can be either fixed or moving or touching the surface or remote from it. In $(t,{\vec{x}})\,\in\,[t_{1},t_{2}]\times\Omega$, the scalar field $u$ is assumed to obey the following second-order wave equation system:

[TABLE]

where $\mathtt{c}$ denotes the signal velocity in $\Omega$, $f(t)$ is the time dependence of the signal transmission, $\vec{x}^{({0})}$ is the point (spatial location) of the transmission, and $\delta(\vec{x}-\vec{x}^{({0})})$ is a Dirichlet delta function, i.e., it is zero everywhere except at $\vec{x}=\vec{x}^{({0})}$ and satisfies the integral identity $\int_{\Omega}\delta(\vec{x}-\vec{x}^{({0})})\,d\Omega=1$. The left-hand side of Equation (1) is the standard hyperbolic operator of the wave equation [33], and the right-hand side represents a point source transmitting an isotropic signal pulse. The isotropic radiation pattern is altered, if the source is placed in front of a reflector (with $\mathtt{c}=0$). In this study, a loudspeaker profile is used as a reflector. Additionally, it is assumed that only the set of points belonging to $\Omega$ can transmit a signal and that there are no other signal sources present. Defining a new variable $\vec{g}(t,\vec{x})=\int_{0}^{t}\nabla u(\tau,\vec{x})\,d\tau$ and $h(t)=\int_{0}^{t}f(\tau)\,d\tau$, the resulting system is of the following first-order form:

[TABLE]

This system can be discretized spatially using the finite element method (FEM) [34] and temporally via the finite-difference time-domain (FDTD) method [35] leading to so-called leap-frog formulae which enable the simulation of the complete wave [15].

II-A2 Measurement model

In this study, we aim at localizing perturbations in the inverse of the velocity distribution $\mathtt{n}={1}/{\mathtt{c}}$ based on the travel-time $T$ of the signal, as this is known to provide robust information about the velocity [21, 2]. The wave $u_{m}(t,\vec{x})$ measured at $\vec{x}$ is a sum of the propagating wave $u(t,\vec{x})$ and a random noise term $\varepsilon(t,\vec{x})$, i.e.,

[TABLE]

The noise can include both modeling and measurement errors. In order to minimize the effect of the noise, a priori information should be used, for example, to determine the most relevant time interval for each measurement point. However, deriving a complete statistical model for $\varepsilon(t,\vec{x})$ can be difficult due to potential but unknown error sources, such as reflections, refractions and absorption. There is also no unique way to obtain $T$ based on the waveform measurement $u_{m}$ [2].

II-A3 Integrated and thresholded travel-time estimates

In this paper, we study integration and thresholding as two alternative techniques for estimating $T$. Given the measured wave $u_{m}$, the integrated travel-time (ITT) estimate for a signal received at $\vec{x}$ within the time interval $[t_{1},t_{2}]$ can be defined via the formula [36]

[TABLE]

To interface the mathematic model and real signals, some decisions have to be made. These include how to decide the time interval [$t_{1},t_{2}]$ and whether to use the ITT or an alternative strategy, e.g., the thresholded travel-time (TTT) estimate. ITT and TTT were calculated as follows:

Detect the time value $t_{0}$ (TTT) where the amplitude reaches a pre-selected threshold value. 2. 2.

Set $t_{1}=t_{0}-\tau_{1}$ and $t_{2}=t_{0}+\tau_{2}$, where $\tau_{1}$ and $\tau_{2}$ are auxiliary a priori parameters ensuring that the essential part of the signal pulse will be contained in $[t_{1},t_{2}]$. Then, obtain ITT through Equation (6).

Here, $\tau_{1}$ extends the inspected time interval in the reverse direction, that is, (before) the signal detection point $t_{0}$. Parameter $\tau_{2}$ cuts the signal based on the a priori information of the pulse length in order to prevent noise due to reflections. Figure 2 visualizes the calculation of ITT and TTT for measurement and simulated data. Figure shows that in the measurements, the travel-time can only be calculated if the starting point (control pulse) is also measured. The measured travel-time can be compared to the simulated case (Figure 2(b)) where the starting point (zero-point) is defined exactly.

Based on the measurements, exactly one travel-time value $T$ is calculated for each transmitter-receiver position pair within the measurement configuration. The resulting set of values is referred to as the measurement data. The simulated signal refers to the wave that can be obtained via the FDTD method [15] using an initial (constant) estimate for the unknown parameter.

II-A4 Path integrals and ray-tracing

In addition to FEM/FDTD computations, a forward simulation for the travel-time can be obtained via an integral of the form

[TABLE]

where $\mathcal{C}$ is the signal path [21]. In this approach, the prediction for $\mathcal{C}$ has a central role. This can be done, e.g., based on Snell’s law of reflection and reftraction [37]. We use the approximation (7) in order to find a reconstruction of $\mathtt{n}$. We assume that the signal paths $\mathcal{C}_{1},\mathcal{C}_{2},\ldots,\mathcal{C}_{M}$ are straight line segments and that $\mathtt{n}$ is of the form $\mathtt{n}=\mathtt{n}_{p}+\mathtt{n}_{0}$, where $\mathtt{n}_{0}$ is an a priori known constant background distribution and $\mathtt{n}_{p}$ is an unknown perturbation. A discretized version of (7) can be obtained by subdividing the computational domain into pixels $P_{1},P_{2},\ldots,P_{N}$ and estimating $\mathtt{n}_{p}$ with a pixelwise constant distribution. Using ${\bf x}$ to denote the vector of pixel values, the equation (7) can be written in the discretized form

[TABLE]

where $L_{i,j}=\int_{\mathcal{C}_{i}\cap P_{j}}\,ds$ and $y_{0}$ is a simulated vector estimated from the travel-time data corresponding to $\mathtt{n}_{0}$.

II-B Inversion procedure

To invert the data, we use a classical total variation based regularization technique, which aims at minimizing (Appendix) the regularized objective funcion $\Psi({\bf x})=\|{\bf y}-{\bf y}_{0}-{\bf Lx}\|_{2}+\alpha\|{\bf D}{\bf x}\|_{1}$ via the iterative recursion procedure

[TABLE]

where ${\bf\Gamma}^{(k)}$ is a diagonal weighting matrix defined as ${\bf\Gamma}^{(k)}=\hbox{diag}\left(\|{\bf D}{\bf x}^{(k)}\|\right)^{-1}$ and ${\bf D}$ is a regularization matrix whose entries satisfy

[TABLE]

with $\delta_{i,j}=1$ if $i=j$, and zero if otherwise. In ${\bf D}$, the first term penalizes the norm of ${\bf x}$ and the second one the total variation, i.e., the total sum $\sum_{i,j}\int_{P_{i}\cap P_{j}}|x_{i}-x_{j}|ds$ in which each non-zero term equals to the total jump of ${\bf x}$ between two adjacent pixels multiplied by the pixel’s side-length [1]. This simple inversion approach usually converges sufficiently in a relatively low number of iteration steps, e.g., five. The value of the parameter $\beta$ determines the balance between the regularization matrices. A small value for $\beta$ leads to inverse estimates with low total variation and larger values can be expected to result in well-localized estimates [37, 15, 36].

II-C Test setup and scenario

II-C1 Domain

An acoustic setup with one speaker and two microphones was utilized to gather the experimental data. The computational domain $\Omega$ was a 59 cm diameter disk drawn on a 0.5 cm thick soft foam covering (Figure 3). In the experiment, the locations of three foam cylinders A, B and C with diameters 15, 10 and 10 cm, respectively, were to be detected. These were placed in the disk $\Omega$ in the upright position and apart from each other. The perimeter of $\Omega$ contained 64 equally spaced control points 1–64 for localizing the transmitter and the receiver of the signal. A top-down view of the experiment setup is included in Figure 3. The diameters and positions of $\Omega$ and A–C can be found in Table II.

II-C2 Signal

The signal pulse transmitted (Figure 4) was of the form

[TABLE]

and $f(t)=0$, in all other times. The resulting signal with 5.8 kHz center frequency was transmitted from the points 1, 24 and 43 using an active two-way speaker. Corresponding to each point of transmission, the signal was recorded for the 57 centermost control points opposite the transmitter. Based on the measurements, two travel-time data sets (I) and (II) (Figure 3) were formed. In (I) (dense set), the 47 centermost points were utilized in the final data set consisting of $3\times 47=141$ individual travel-time values, while only every other point of the set (I) was included in (II) (sparse set) resulting in $3\times 24=72$ travel-time values. The measured waveform data can be found in Figure 5. The simulated data were obtained using the FDTD method using the constant $334$ m/s as the signal velocity (the speed of sound in air at 20 $\mbox{}^{\circ}$C). Both the measured and simulated data have been included as supplementary material.

II-C3 Equipment

The data recording device was an ordinary laptop computer equipped with an external sound card interface. Two-channel audio data was recorded using the waveform audio file format (WAV), 24 bit resolution and a 48 kHz sample rate. To eliminate any possible hardware or software based delays, a control pulse was recorded by one microphone placed 3.5 cm directly above the speaker. The actual data pulse was recorded with the other microphone positioned on the perimeter of the target area in the radial direction. Figure 3 shows the test setup with both microphones near each other. The technical data of the hardware used in the setup can be found in Table I.

II-C4 Noise

Figure 4 includes a comparison between the transmitted signal pulse and the one received at the perimeter of the (empty) test domain $\Omega$ without the foam cylinders placed inside. The tail of the signal received was observed to contain noise due to echoes and inaccuracies in transmission (ringing). Based on a comparison with the original signal the noise peaks were estimated to be mainly 10 dB (greyed area) below the main peak.

II-C5 Relevance

The relevance of the test setup with respect to a real travel-time tomography application is the following. The transmitter sends a waveform signal pulse at a known position which is then recorded by the receiver in a known position, and the resulting signal recordings are then sent via a communication link. Further processing will consist of compression or some other processing technique, such as calculating travel-time values.

II-D Experimental FPGA hardware

The travel-time calculation was implemented using a high-level synthesis of hardware on an FPGA development board (Altera DE2) with the typical performance of an embedded signal acquisition system. FPGA enables fast processing of data which is essential in waveform imaging because of the massive data input needed to record a complete wave. A high-level data flow chart of the implementation has been included in Figure 6.

In this study, we investigated different bit resolution, threshold and normalization levels. These were set in a separate parameters file. A functioning Matlab script was first transformed directly into C code and the result was then further modified to accept parameters. The hardware was generated by using the Mentor Graphics Catapult High-Level Synthesis (HLS) tool [38]. HLS is a method where digital hardware is generated from a high level programming language, such as C [39]. The tool takes the modified, algorithmic C code and generates register transfer level (RTL) code, which is synthesized as digital logic on an Field Programmable Gate Array (FPGA) chip [40].

The basic hardware for signed 16-bit input (audio) data was implemented first for ITT and TTT, and then adapted separately for the signed 8-bit input data. The software on computer uses common datatypes [41], but efficient hardware requires accurately defined bit resolutions for all inputs, outputs and intermediate variables [42]. The HLS tool does not offer a way to calculate the maximum bit resolution for an arbitary integer number. An example of this is the divisor in Equation (6), which has the sum of the squares of the numbers for a part of the signal. Knowing the bit resolution of the element and thus the maximum values, the bit resolution for the square could be calculated. The 8-bit and 16-bit versions of the hardware required individual optimizations.

II-E Numerical experiments

In the numerical experiments, the initial guess for the signal velocity distribution was set to be $\mathtt{c}=334$ m/s (the speed of sound in air at 20 $\mbox{}^{\circ}$ C). The number of inverse iteration steps was set at three and the regularization parameters $\alpha$ and $\beta$ were both given the value 0.01 which, based on our preliminary tests, is a reasonable approximaitoion of the midpoint of the interval of the workable values. To evaluate the data processing artifacts, bit resolutions of the signed 16 and the signed 8 bits were used as in the A/D conversion units of practical applications [43].

The maximum amplitude in the measurement data set was normalized to a given level $\nu$ dB FS (full of the bit scale) and the other signals proportional thereto. The control signal was normalized to 0 dB FS for each individual measurement point. 16-bit preprocessing was evaluated at 100 % ($\nu=0$ dB FS) normalization and 8-bit at 100 % ($\nu=0$ dB FS) and 6 % ($\nu=-24$ dB FS) normalization. The 6 % level represents an extreme case where the signal noise is large, thus simulating either a weak reception of transmitted signals or an insensitive receiver. Two different threshold levels 90 % and 70 % for initial signal detection were tested corresponding to around -1 dB FS and -3 dB FS of the maximum value of a normalized signal.

The data vector ${\bf y}$ was obtained as the difference between the measurement and simulation based travel-time both of which were computed using either the ITT or TTT approach. The interval of the ITT was determined by $\tau_{1}=5$ and $\tau_{2}=250$ (samples at 48 kHz) resulting in a total length of 256 samples (5.3 ms). TTT was evaluated using two alternative strategies TTT 1 and TTT 2. In the former, ITT was applied to the control pulse and TTT to the data pulse. In the latter, the travel-time of both the control and data pulse was computed via TTT. The motivation for investigating TTT 1, was the potential situation in which the simulated and measured travel-time are obtained via different techniques, e.g., due to different suppliers of computer software and measurement equipment.

The reconstructions were analyzed by measuring the Relative Overlapping Area (ROA) and the minimum relative overlap $\hbox{ROA}_{\hbox{\scriptsize min}}$ between the foam cylinders $\mathcal{S}_{\hbox{\footnotesize\bf A}}$, $\mathcal{S}_{\hbox{\footnotesize\bf B}}$, $\mathcal{S}_{\hbox{\footnotesize\bf C}}$ the set $\mathcal{S}_{\hbox{\footnotesize rec}}$ in which the value of the reconstruction was less than a fixed limit such that $\hbox{Area}(\mathcal{S}_{\hbox{\footnotesize rec}})=\hbox{Area}(\mathcal{S}_{\hbox{\footnotesize\bf A}}\cup\mathcal{S}_{\hbox{\footnotesize\bf B}}\cup\mathcal{S}_{{\hbox{\footnotesize\bf C }}})$. ROA and $\hbox{ROA}_{\hbox{\scriptsize min}}$ were calculated as given by the equations

[TABLE]

The reconstructions were computed using a laptop computer equipped with the 2.8 GHz Intel Core i7 processor 2640M and 8 GB of RAM. Computing a single reconstruction took 9 seconds of CPU time.

III Inversion results

Figures 7-10 show the test area recovery (inversion) results. These figures consist of image pairs. The outlines of the three foam cylinders are superimposed on the top of the images. For each image the left side image is the actual result image. From that image, an area with strongest values is collected so that it has the same area as that of the foam cylinders. The right side image shows only the solid black areas.

ITT was found to yield robust results with respect to the A/D conversion bit resolution, signal normalization, threshold, and sparsity of the measurements. The highest ROA ($\hbox{ROA}_{\hbox{\scriptsize min}}$) obtained with ITT was 64 % (56 %) and the lowest one 54 % (38 %). TTT was observed to be advantageous under optimized conditions with respect to the bit resolution, normalization, threshold, and density of the data. It was also significantly more sensitive than ITT to variations in any of these parameters. For TTT 1, the highest and lowest ROA ($\hbox{ROA}_{\hbox{\scriptsize min}}$) were 70 (57 %) and 17 % (0 %), respectively. For TTT 2, these values were 71 (68 %) and 25 % (0 %), respectively. The results obtained with TTT 2 were slightly superior to those achieved using TTT 1. The 16-bit signed integer travel-time calculation module produced essentially the same outcome as 64-bit floating point arithmetics.

The HLS method was found to work appropriately in developing the FPGA hardware. A complete prototyping cycle took a few hours, most of which time was spent implementing a new feature in C. It was found that since the algorithmic C-like code is much more maintainable than traditional designs using VHDL, further changes and reuse of the code can be done more easily via HLS. A summary of the synthesis results for the 8-bit and 16-bit designs for the ITT and TTT is presented in Table III. The hardware for ITT was almost as fast as that of TTT, and the difference in chip area between these two methods was not significant.

IV Discussion

This paper focused on developing the processing and inversion of waveform tomography data for applications in which the signal wavelength is close to the diameter of the details that are to be recovered. In particular, FPGA-based hardware was used. We compared the integrated and thresholded travel-time (ITT and TTT) in numerical experiments in which three foam cylinders were to be localized based on an experimental 5.8 kHz acoustic signal. We tested a 16-bit and 8-bit analog-to-digital (A/D) conversion together with two different threshold criterions and normalization levels. As reference applications of this study, we considered (1) microwave and ultrasonic computed tomography (MCT and UCT), (2) tomography of small solar system bodies (SSSBs), in particular, the CONSERT experiment and (3) ultrasonic/microwave detection of concrete defects.

IV-A Experiments

The recovery of the test object locations on the target area by using the described inversion methods was found to work appropriately. Our results show that ITT is more stable than TTT if the signal quality decreases. This can been seen from the stability of the relative overlapping area (ROA) percentages (53 - 64 %). These percentages are comparable to our previous research on waveform inversion within a 2D domain. In [15] the best ROA percentage 71 % was obtained by using the full wave data. Using TTT the best recovery result was exactly the same 71%. ITT was 7 % less accurate, but more stable when data preprocessing was performed using source data with lesser bit resolution.

The results show that ITT is invariant with respect to source data bit resolution reduction in the time-domain and the level of thresholding used to locate the signal pulse from recorded audio data. ITT was also found to be more reliable than TTT with respect to the normalization of the signal and the sparsity of the measurements. We normalized the 8-bit signal to two different levels, 100 % (0 dB FS) and 6 % (-24 dB FS) amplitude, to simulate weak receivers, i.e., to decrease signal-to-noise ratio (SNR). This is relevant in applications where the signal quality is reduced. For example, in astro-geoscientific applications, the signals can be weak in some directions. This was evident from the CONSERT experiment [13, 44] in which the power and quality of the received signal varied significantly depending on the direction of the measurement [12]. The result concerning the sparsity of the measurements is essential for CONSERT, and other applications in which not all the data can be gathered.

Based on the results, we suggest that ITT can be superior to TTT with respect to the robustness of the inversion. It also seems obvious that TTT might achieve a higher ROA than ITT for a high-quality signal and a well-chosen threshold parameter. Namely, under optimal conditions TTT filters out the noisy tail part of the signal that is present in the computation of the ITT. This advantage is utilized, e.g., in the first-arrival sound speed inversion [3]. The present results suggest that the ITT method does not significantly diminish the inversion quality, and indeed, it can increase the reliability of the results with respect to the uncertainty factors and incompleteness of the data.

IV-B Applications

Tomography applications differ from each other by the speed and length of the electromagnetic or sound waves. The feature size $d$ that can be detected is dependent on the wave length $\lambda$ (e.g. $\lambda/2$). In Table IV, the significance of the test setup with respect to the present reference applications (1)–(3) has been summarized based on $d$, $\lambda$, the diameter of the target domain $D$ and the ratios $d/\lambda$ and $D/\lambda$.

The values $d$, $\lambda$, and $D$ utilized in Table IV can be reasoned as follows. (1) In MCT, signal frequencies 1–6 GHz are being used and the diameter of the sensor ring can be, e.g., $D=15$ cm [26, 4]. At 5 GHz, the wavelength is $\lambda\approx 1.9$ cm corresponding to the relative permittivity of the breast $\varepsilon_{r}\approx 10$ [45, 46]. Medical UCT utilizes frequencies in the range of 1 - 20 MHz [47]. Using the 1 MHz value and the speed of 1500 m/s for ultrasound propagation in human tissue [48], the wavelength is $\lambda=1.5$ mm. In UCT and MCT, the feature to be recovered can be, e.g., a small T1 or T1a tumor which can have a maximum diameter of $d=2$ cm and $d=0.5$ cm, respectively. (2) A signal frequency of 10 MHz has been suggested for the tomography of SSSBs [49] matching roughly with the wavelength $\lambda\approx 15$ m ($\varepsilon_{r}\approx 4$, e.g., for dunite and kaolinite [50]) which is the estimated resolution of the CONSERT experiment [13, 44]. The diameter of the SSSB can be for example $D=150$ m. (3) In ultrasonic material testing for concrete the velocity of the wave is 3500 m/s [51]. At a typical ultrasound frequency of 100 kHz the wavelength of the signal is $\lambda$ = $v$/$f$ = 3.5 m and the diameter of the concrete beam being tested could be, for example, $D=30$ cm. A fault or crack inside the beam could be $d=3$ cm.

Based on Table IV, the present experiment setup can be considered applicable to all application contexts (1)–(3). In MCT and UCT the best match to the test geometry is obtained with the T1 and T1a tumor size, respectively. Obviously, the difference between ITT and TTT can be less significant in applications, where the wave length is likely to be considerably smaller than the smallest detail to be detected. As an additional comparison, in the tomography of the ionosphere, the speed of the electromagnetic wave is around the speed of light and a typical frequency used in tomography is between 120 MHz and 400 MHz. Using a value of 200 MHz the wavelength is $\lambda=1.5$ m. The size of the ionosphere is up to $D=1000$ km from the surface of the Earth and a typical vertical feature to be recovered can be $d=100$ km, for example [52]. Consequently, $d/\lambda$ and $D/\lambda$ can be 100 and 670, respectively. This suggests that our test scenario might be too different from the tomography of the ionosphere to be able to draw inferences from the results.

The noise peaks in the experimental data were estimated to be mainly 10 dB below the maximum peak. This can be considered as appropriate for the applications (1)–(3). For a MCT imaging system, the relative reconstruction error has been shown to stay under 10 % for amplitude errors down to 10 dB SNR [45]. The total noise peak level of around 20 dB was observed in the CONSERT experiment [12]. In concrete testing, structural noise [53], e.g., echoes from walls, can be significant resulting in noise peaks that can be comparable to the main peak.

IV-C Hardware

The use of high-level synthesis (HLS) to develop the test hardware was found to be essential, as the specifications for the travel time calculation changed during the implementation phase and so changes to the hardware had to be made quickly. We aimed at a very direct workflow in implementing the hardware, and thus some optimization methods in the design partitioning and in the HLS tool were not utilized. In addition, having the different bit resolution implementations in separate files increased the work when changes had to be made. The resulting digital hardware was synthesized on an FPGA platform for demonstration and to facilitate further development of the data gathering system towards a laboratory instrument.

Many operations in the computation scripts, such as the normalization of values, require divisions. In the ITT calculation formula (6) there is a large-valued divisor that first sums the squares of the values and then divides the sum with that value. Catapult generates divisions as combinatorial logic if written directly as it is in typical C code [42]. Combinatorial logic [40] is not synchronized with the hardware clock, which makes it unreliable in use. This gives a false sense of flexibility in the HLS tool and the user has to know what is being generated. Catapult has a math library which has synthesizable, basic algorithmic division operators for integers. Changing all divisions to use division operators offered by this library improved the results, and resulted in a working design for most division operations whereas combinatorial dividers did not.

Travel-time calculation can be seen as an extreme form of compression in time domain. There are more compression methods such as filtering in the frequency domain. Our work used reduced the bit resolution for input audio data. Calculation of the travel times was also performed with integer arithmetics. These reduce the accuracy of the results. bit resolution reduction is relevant when using FPGAs in general. This is because vendor-provided hardware multipliers available for use on FPGA chips such as Altera’s Cyclone II have limited bit resolutions [54]. For example, a FIR filter can use these hardware multipliers on FPGAs as in [55]. Because the operator is a multiplication that can increase the required bit resolution to store the intermediate results, the original bit resolution of the data has to be reduced. One system [55] has implemented Fast Fourier Transform (FFT) filtering in the frequency domain and inverse FFT (IFFT) signal reconstruction back to the time domain. In currently available FPGA DSP chips, the available bit resolutions for multipliers are much larger [56], but reduction in bit resolution can still be required. This makes the invariance of the source data bit resolution of our ITT method an important point of interest. An example of decreased signal quality is found in Ground Penetrating Radar (GPR) applications. In [43] an instrumentation system with either 8-bit or 16-bit Analog-to-Digital conversion was employed. This shows that bit resolution limitations are also found in instrument hardware.

IV-D Outlook

Finally, the present results indicate various directions for future work. We will study hardware-level solutions regarding (i) processing and (ii) inverting waveform data. (i) Investigating harware constraints utilizing a more sophisticated statistical travel-time detection model, such as the akaike information criterion [57], would be an interesting goal. The calculation of travel-time values is only one data preprocessing method that can be used, and mathematical methods to accomodate other types of filtering, such as compressing in the frequency domain could be developed. The goal in preprocessing depends on the application: in the space environment minimal data transfer between the sensors and the computation unit is important for the limited communication capacity available, whereas in biomedical and civil engineering the main objective can be to optimize the speed of the procedure in order to allow recording as much data as possible. (ii) We aim to develop inversion approaches utilizing the FPGA environment, so extending the current study of hardware-level constraints to include inversion algorithms is essential.

Acknowledgements

M.T. and S.P. were supported by the Academy of Finland Key Project 305055 and the Academy of Finland Centre of Excellence in Inverse Problems Research.

IV-E Appendix

In the inversion procedure (9), the matrix ${\bf D}$ is symmetric and can thus be diagonalized. When $\beta>0$, ${\bf D}$ is also positive definite and also invertible. Hence, one can define $\hat{\bf x}={{\bf D}}{\bf x}$ and $\hat{\bf L}={\bf L}{{\bf D}^{-1}}$. Substituting $\hat{\bf x}$ and $\hat{\bf L}$ into (9) leads to the following form [15, 16, 1]

[TABLE]

which can be associated with alternating conditional minimization of the function $H(\hat{\bf x},\hat{\bf z})=\|\hat{\bf L}\hat{\bf x}-{\bf y}\|^{2}_{2}+{\alpha}\sum_{j=1}^{M}\frac{\hat{x}^{2}_{j}}{\hat{z}_{j}}+{\alpha}\sum_{j=1}^{M}\hat{z}_{j}$ in which $z_{j}>0$, for $i=1,2,\ldots,M$. As $H(\hat{\bf x},\hat{\bf z})$ is quadratic with respect to $\hat{\bf x}$, the conditional minimizer $\hat{\bf x}^{\prime}=\arg\min_{\hat{\bf x}}H(\hat{\bf x}\mid\hat{\bf z})$ is given by the least-squares solution of the form $\hat{\bf x}^{\prime}=(\hat{\bf L}^{T}\hat{\bf L}+\alpha\hat{\bf\Gamma}_{\hat{\bf z}})^{-1}\hat{\bf L}^{T}{\bf y}$, $\hat{\bf\Gamma}_{\hat{\bf z}}=\hbox{diag}(\hat{\bf z})^{-1}$ and $\hat{\bf\Gamma}_{0}={\bf I}$. At the point of the conditional minimizer $\hat{\bf z}^{\prime}=\arg\min_{\hat{\bf z}}H(\hat{\bf z}\mid\hat{\bf x})$, the gradient of $H(\hat{\bf z}\mid\hat{\bf x})$ vanishes with respect to $\hat{\bf z}$, i.e.,

[TABLE]

Hence, the global minimizer can be estimated via the following alternating iterative algorithm.

Set $\hat{\bf z}_{0}=(1,1,\ldots,1)$ and $\ell=1$. For a desired number of iterations repeat the following two iteration steps. 2. 2.

Find the conditional minimizer $\hat{\bf x}_{\ell}=\arg\min_{\hat{\bf x}}H(\hat{\bf x},\hat{\bf z}_{\ell-1})$. 3. 3.

Find the conditional minimizer $\hat{\bf z}_{\ell}=\arg\min_{\hat{\bf z}}H(\hat{\bf x}_{\ell},\hat{\bf z})$.

The sequence $\hat{\bf x}_{1},\hat{\bf x}_{2},\ldots$ produced by this algorithm is identical to that of (14) and ${\bf x}_{\ell}={\bf D}^{-1}\hat{\bf x}_{\ell}$ equals to the $\ell$-th iterate of (9). If for some $\ell<\infty$ the pair $({\bf\hat{x}}_{\ell},{\bf\hat{z}}_{\ell})$ is a global minimizer of $H({\bf x},{\bf z})$, then, since $(\hat{z}_{\ell})_{j}=|(\hat{x}_{\ell})_{j}|$, $j=1,2,\ldots,M$, then it is also the minimizer of

[TABLE]

Consequently, ${\bf x}_{\ell}={\bf D}^{-1}\hat{\bf x}_{\ell}$ is the minimizer of $\Psi({\bf x})$, as $\Psi({\bf x})=\hat{\Psi}(\hat{\bf x})$. Furthermore, the minimizer of $\hat{\Psi}(\hat{\bf x})$ is also the 1-norm regularized solution of the linearized inverse problem.