Measurement uncertainty analysis of field-programmable gate-array-based, real-time signal processing for ultrasound flow imaging

Research in magnetohydrodynamics (MHD) aims to understand the complex interactions of electrically conductive fluids and magnetic fields. A promising approach for investigating complex instationary flow phenomena are lab-scale experiments with low-melting alloys. They require a noninvasive flow instrumentation for opaque liquids with a high spatiotemporal resolution, a low velocity uncertainty and a long measurement duration. Ultrasound Doppler velocimetry can achieve multiplane, multicomponential flow imaging with multiple linear ultrasound arrays. However the average raw data output amounts to 1.2 GBs⁻¹ at a frame rate of 33 Hz in a typical configuration for 200 transducers. This usually prevents long-duration measurements when offline signal processing is used.

In this paper, we propose an online signal-processing chain for pulsed-wave Doppler velocimetry that is tailored to the specific requirements of flow imaging for lab-scale experiments. The trade-off between measurement uncertainty and computational complexity is evaluated for different algorithmic variants in relation to the Cramér–Rao bound. By utilizing selected approximations and parameter choices, a prepossessing could be efficiently implemented on a field-programmable gate array (FPGA), enabling a typical reduction of the data bandwidth of 6.5:1 and online flow visualization. We validated the performance of the signal processing on a test rig, yielding a velocity standard deviation that is a factor of 3 above the theoretical limit despite a low computational complexity.

Potential applications for this signal processing include multihour flow measurements during a crystal-growth process and closed-loop velocity feedback for model experiments.

1 Introduction

Many important industrial processes, such as continuous steel casting and photovoltaic wafer production, involve metal or semiconductor melt flows. The quality of the product and the energy efficiency of the process strongly depends on the flow behavior of the liquid (Müller and Friedrich, 2010; Gardin et al., 1995; Yasuda et al., 2007). A noncontact way of influencing the flow of electrically conductive melts is the application of magnetic fields that introduce Lorentz forces to the fluid. Investigating the interaction of a magnetic field with the flow pattern and optimizing the spatiotemporal structure of the magnetic field for different applications are subjects of ongoing research in magnetohydrodynamics (MHD). Besides numerical simulations, low-temperature, model-scale experiments are important tools for MHD investigations (Eckert et al., 2007b). They often require advanced flow instrumentation for visualizing complex and instationary flows in opaque liquids. A typical set of requirements for MHD research are as follows:

• Noninvasiveness – the influence of the instrumentation to the flow should be negligible (Eckert et al., 2007a).

• Flow imaging capability – the fluid's velocity should be visualized in multiple planes (2D) with two or three velocity components (2c or 3c) in order to adequately represent complex flow patterns.

• Spatial resolution – the relevant flow structures have to be resolved, typically in the range of 10 mm (Timmel et al., 2011).

• Temporal resolution – fluctuations (typically at 1 … 5 Hz) have to be resolved in order to capture instationary flows (Timmel et al., 2011).

• Long measurement duration – flow phenomena on different timescales should be adequately captured; for instance, rapid spontaneous changes of the flow regime in a rotating flow (Galindo et al., 2017) or in multihour model experiments of the semiconductor crystallization process (Thieme et al., 2017).

• Capability of near-wall measurements – in typical MHD experiments, the metal melt is contained in a vessel. The vicinity of the wall is especially important because the Lorentz force is often concentrated in this region. Contrary to, for instance, medical applications, the walls can be seen as completely stationary in most cases.

• Online capability – conducting long-running MHD experiments requires the ability to examine the data during the duration of the measurement. Some model experiments in the semiconductor crystallization process even benefit from an active control of parameters, like magnetic field intensity and temperature gradient, based on the feedback from online velocity data to stabilize the flow (Thieme et al., 2017).

A measurement system for flow mapping of opaque liquids, namely the ultrasound array Doppler velocimeter (UADV; Nauber et al., 2013a, b), was presented in previous publications. It extends the pulsed-wave Doppler principle (Takeda, 1986; Baker, 1970) by employing multiple linear sensor arrays to achieve multiplane, two-componential flow imaging. The sensors are designed to achieve a lateral resolution of ≈3 mm in Galinstan (GaInSn). A combination of spatial- and time-division multiplexing allows one to parallelize the scanning process for a planar velocity map; hence increasing the temporal resolution compared to a strict sequential scan. However, online processing of the data for 200 transducer elements simultaneously on 32 channels at a temporal resolution typically of 33 Hz overburdens PC-based hardware with 1.2 GBs⁻¹. Therefore, only discontinuous offline measurements could be performed with a limited duration of a few seconds. This severely impedes the usability of the UADV in the context of MHD experiments and restricts the investigations into stationary or periodic flows.

Although several investigations on the measurement uncertainty of Doppler velocity estimation methods for laser-based instrumentation (Fischer et al., 2010), for flow-rate measurements in a pipe (Furuichi, 2013), and for blood-flow measurements in the human body (Lovstakken et al., 2007) have been performed, no comprehensive measurement uncertainty budget in the context of instrumenting an MHD experiment has been presented to the knowledge of the authors.

This paper provides a signal-processing chain that is tailored to the specific requirements of MHD model experiments and shows a real-time implementation using a field-programmable gate array (FPGA). It enables the UADV system to perform long-duration measurements with high frame rates and online flow visualization. Furthermore, we evaluate the measurement uncertainty of the whole UADV system in the context of MHD experiments and present an uncertainty budget according to the methodology proposed by the “Guide to the expression of uncertainty in measurement” (GUM; JCGM, 2008) for a typical configuration.

2 Pulsed-wave ultrasound Doppler velocimetry

2.1 Measurement principle

In pulsed-wave ultrasound Doppler velocimetry (PW–UDV), short bursts are emitted periodically with a pulse repetition frequency f_PR (Baker, 1970). The emission times $t_{\mathrm{s}}=n_{\mathrm{b}}/f_{\mathrm{PR}}$ span the so-called slow-time axis t_s, with n_b=0…N_EPP being the bursts number. The emitted bursts usually consist of N_periods periods of a sinusoidal wave, with the frequency f₀. As the bursts travel through the fluid, scattering particles reflect a fraction of the signal back to the ultrasound transceiver. The received echo signal z(t_f, t_s) is acquired, starting from the emission time along the fast-time axis t_f. Figure 1 depicts an example of the echo signal for a single moving scattering particle.

The movement of a scattering particle leads to a phase shift of the echo signal between multiple burst emissions (Kasai et al., 1985). The mean phase shift per time unit, expressed as mean frequency f_d, is related to the velocity v for a given speed of sound c by the following:

$$\begin{array}{}\text{(1)}& {\displaystyle }v=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{f_{\mathrm{d}}}{f_{\mathrm{tx}}}}c,\end{array}$$

with f_tx denoting the mean frequency of the received signal burst and c≫v. The mean phase shift per time unit f_d can be interpreted as a Doppler frequency shift f_d (Kasai et al., 1985); hence the name Doppler velocimetry.

The time since the burst emission t_f corresponds to the distance d between the scattering particle and the transducer according to the following equation:

$$\begin{array}{}\text{(2)}& d={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{t}_{\mathrm{f}}c.\end{array}$$

This allows a spatially resolved flow measurement along the axis of the transducer, given that the scattering particles follow the motion of the fluid with negligible slip. The axial resolution can be estimated with the following (Jensen, 1996):

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }d={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{N}_{\mathrm{periods}}{\displaystyle \frac{c}{f_{\mathrm{0}}}}.\end{array}$$

The lateral resolution Δx is given by the width of the ultrasound beam, which is a result of the transducer geometry, the frequency f₀, and the speed of sound c in the fluid. The temporal resolution Δt of the velocity measurement is determined through the following:

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }t={\displaystyle \frac{N_{\mathrm{EPP}}}{f_{\mathrm{PR}}}}.\end{array}$$

Figure 1 PW–Doppler principle: multiple bursts are emitted at $t_{\mathrm{s}}=n_{\mathrm{b}}/f_{\mathrm{PR}}$, with a repetition rate f_PR. This constitutes the slow-time axis t_s. After emission, the received echo signals are sampled with a frequency f_s along the fast-time axis t_f. The echo signal phase shift corresponds with the velocity of the scattering particles in the fluid. An example of a single particle moving away from the transducer is given.

2.2 Ultrasound array Doppler velocimeter

The ultrasound array Doppler velocimeter (UADV) is a modular research platform developed at the Laboratory of Measurement and Sensor System Technique (MST) for flow imaging in opaque liquids with PW–UDV. It is flexible and especially well suited for instrumenting a wide range of experiments in the field of MHD. The hardware of the UADV consists of individually configurable modules driving 25 ultrasound transducers each. It can be scaled to support up to 200 transducers in various configurations; for instance, in four linear arrays which can be individually parameterized regarding ultrasound frequency, pulse shape and length, and pulse-repetition frequency (Nauber et al., 2016; Büttner et al., 2013).

A module of the UADV consists of an arbitrary function generator and a power amplifier for generating parameterizable burst signals which are routed through a programmable switching matrix and a transmit/receive switch to the transducers. The received echo signals are amplified with a parametric gain and routed to the digitization unit. A single microcontroller-driven control unit provides the overall synchronization and the communication with the host PC. Using a combined spatial- and time-division multiplexing scheme, an ultrasound transducer array can scan a measurement plane at higher rates than a strict sequential scan. The UADV supports four independent digitization channels per module. The detailed description of the measurement system is given in Nauber et al. (2016).

3 Signal processing for velocity estimation

3.1 Overview

The signal processing for PW–UDV can generally be classified into wideband and narrowband techniques; a comprehensive comparison is given by Torp et al. (1993). While parts of the signal processing in the radio frequency (RF) band can be realized in analog circuitry (Shung, 2015), fully digital implementations have found widespread use in the last decades because of the availability of fast digitizers and the increased flexibility and robustness of such approaches. To simultaneously handle a high number (e.g., 32) of channels through a fully digital signal-processing chain, very large data bandwidths have to be processed. This can be achieved by utilizing the parallel-processing capability of a field-programmable gate array (FPGA). Especially narrowband algorithms are very suitable for FPGA-based implementations, due to their low computational complexity (Alam and Parker, 2003; Loupas et al., 1995a). Therefore, this paper focuses on investigating the most common narrowband method, the velocity estimator by Kasai et al. (1985) and the extensions proposed by Loupas et al. (1995b).

A typical narrowband signal-processing chain is shown in Fig. 2. In this fully digital realization, the slow time t_s is sampled for each burst n_b at $t_{\mathrm{s}}=n_{\mathrm{b}}\frac{\mathrm{1}}{f_{\mathrm{PR}}}$ and the fast time t_f is sampled with a frequency f_s as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{z}_{\mathrm{raw}}^{\prime }(k,n_{\mathrm{b}})& =z(t_{\mathrm{f}}=k/f_{\mathrm{s}},t_{\mathrm{s}}=n_{\mathrm{b}}/f_{\mathrm{pr}}),\\ k& =\mathrm{0},\mathrm{1},\mathrm{\dots }K,n_{\mathrm{b}}=\mathrm{0},\mathrm{1},\mathrm{\dots }{N}_{\mathrm{EPP}}.\end{array}\end{array}$$

The signals are then bandpass filtered to reduce noise contributions outside of the bandwidth of the transmitted ultrasound signal. A quadrature demodulation is performed, consisting of a Hilbert transform and a subsequent down sampling. Static echoes are removed through a clutter reduction filter (CRF) and the velocities are estimated by an autocorrelation.

Figure 2 Signal flow for a typical implementation of the Kasai algorithm: raw echo signals are bandpass filtered, quadrature demodulated and down sampled. The subsequent operations are performed on the complex demodulated signals, namely clutter-reduction filtering and velocity estimation.

3.2 Quadrature demodulation

In order to meet the assumptions of the narrowband signal processing and to reduce the influence of noise, a bandpass filtering is performed as follows:

$$\begin{array}{}\text{(6)}& z^{\prime }(k,n_{\mathrm{b}})=\sum _{n=\mathrm{0}}^{N_{\mathrm{periods}}}{c}_i\cdot z_{\mathrm{raw}}^{\prime }(k-n,n_{\mathrm{b}}),\end{array}$$

with the filter coefficients c_i. In order to maximize the SNR for signals with additive white Gaussian noise, a matched filter is used (Turin, 1960) as follows:

$$\begin{array}{}\text{(7)}& c_i=s_{\mathrm{tx}}(N_{\mathrm{tx}}-i),i\in [\mathrm{1},N_{\mathrm{tx}}],\end{array}$$

with the transmitted signal s_tx with N_tx samples.

The result of the quadrature demodulation is a complex signal ${h^{\prime }}_{\mathrm{unfilt}}(k,n_{\mathrm{b}})$ in the baseband, which can be sampled at a lower rate (reduction by a factor of n_sub) than the raw signal, as follows:

$$\begin{array}{}\text{(8)}& {\displaystyle }{h^{\prime }}_{\mathrm{unfilt}}(k/n_{\mathrm{sub}},n_{\mathrm{b}})=z^{\prime }(k,n_{\mathrm{b}})+j\cdot \widehat{z^{\prime }}(k,n_{\mathrm{b}}),\end{array}$$

with $N_{\mathrm{b}}=\mathrm{0},\mathrm{1},\mathrm{\dots }{N}_{\mathrm{epp}},k=\mathrm{0},\mathrm{1},\mathrm{\dots }K$ and the Hilbert transform signal $\widehat{z^{\prime }}(k,n_{\mathrm{b}})$ (with 90 ^∘ phase shift with respect to z^′).

3.3 Clutter-reduction filtering

A common problem of ultrasound Doppler flow measurements is distinguishing between static echoes originating from the walls (the so-called clutter) and echoes originating from scattering particles. Multiple reflections from the transmitted burst inside the wall superimpose the signal from scatter particles in the vicinity of the wall. For this problem, a multitude of signal-processing methods were proposed, most of them based on digital filters (finite impulse response (FIR) or infinite impulse response (IIR) filters) with various initialization techniques (Lee et al., 2009). With these methods, the clutter is distinguished from the particle echoes by a velocity close to zero, respectively, by a Doppler frequency shift close to zero. Because filtering will influence the spectrum of the signal, a bias may be introduced to the subsequent velocity estimation, depending on the frequency cutoff. For typical MHD experimental setups, the wall can be assumed to be completely stationary (in contrast to, e.g., medical applications where clutter is often constituted by slowly moving tissue; cf. Jensen, 1996); therefore, a steep cutoff at a frequency of zero is desirable. The simplest and computationally most efficient approach is to filter the constant component of the demodulated IQ signal by subtracting its mean value, which is the equivalent of applying a very narrow-band, high-pass filter (Thomas and Hall, 1994; Jensen, 1996; Torp, 1997; Bjaerum et al., 2002) as follows:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{h}^{\prime }(k/n_{\mathrm{sub}},n_{\mathrm{b}})=& {h^{\prime }}_{\mathrm{unfilt}}(k/n_{\mathrm{sub}},n_{\mathrm{b}})\\ & -{\displaystyle \frac{\mathrm{1}}{N_{\mathrm{EPP}}}}\sum _{n_{\mathrm{b}}^{\prime }=\mathrm{0}}^{N_{\mathrm{EPP}}-\mathrm{1}}{h^{\prime }}_{\mathrm{unfilt}}(k/n_{\mathrm{sub}},n_{\mathrm{b}}^{\prime }).\end{array}\end{array}$$

As the filter is noncausal, all N_EPP samples have to be acquired before the result can be computed.

3.4 One-dimensional autocorrelation algorithm

A widely used approach for velocity estimation is the autocorrelation method proposed by Kasai et al. (1985), which operates solely in the domain of IQ-demodulated echo signals and therefore can be implemented very efficiently (Alam and Parker, 2003). It uses the properties of the signals' discrete autocorrelation function as follows:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{R}^{\prime }(\mathrm{\Delta }k,\mathrm{\Delta }{n}_{\mathrm{b}})=& \sum _{m=\mathrm{0}}^{K/n_{\mathrm{sub}}-\mathrm{\Delta }k-\mathrm{1}}\sum _{n=\mathrm{0}}^{N_{\mathrm{EPP}}-\mathrm{\Delta }{n}_{\mathrm{b}}-\mathrm{1}}\\ & h^{\prime }(m,n)\cdot {h^{\prime }}^{*}(m+\mathrm{\Delta }k,n+\mathrm{\Delta }{n}_{\mathrm{b}}),\end{array}\end{array}$$

where its values at a lag of 1 relate to the center of mass of the signal's power density spectrum through the Wiener–Khinchin theorem. As shown by Kasai et al. (1985) and Jensen (1996), the mean Doppler shift f_d can be approximated through evaluating the autocorrelation function at a lag of Δn_b=1 slow-time samples as follows:

$$\begin{array}{}\text{(11)}& f_{\mathrm{d}}\approx \mathrm{1}/T_{\mathrm{PR}}\text{arg}\left(R\right(\mathrm{\Delta }k=\mathrm{0},\mathrm{\Delta }{n}_{\mathrm{b}}=\mathrm{1}\left)\right).\end{array}$$

This autocorrelation computation can be expressed solely by repeatedly multiplying accumulate operations and therefore can be implemented very efficiently. Kasai's method approximates the center frequency f_d of the received signal with the frequency of the emitted signal as follows:

$$\begin{array}{}\text{(12)}& f_{\mathrm{rx}}\approx f_{\mathrm{0}}.\end{array}$$

Being based on a phase estimation, the Kasai algorithm is inherently limited in the maximum measurable velocity. Given the 2π-phase ambiguity in Eq. (11), the measurable velocity range resulting from Eq. (1) is (Jensen, 1996) as follows:

$$\begin{array}{}\text{(13)}& v\in [\pm v_{\max }];v_{\max }={\displaystyle \frac{c{f}_{\mathrm{PR}}}{\mathrm{4}{f}_{\mathrm{0}}}}.\end{array}$$

3.5 Two-dimensional autocorrelation algorithm

An extension of Kasai's autocorrelation method is proposed by Loupas et al. to improve its performance in the following two regards (Loupas et al., 1995b):

1. The assumption of an unchanged center frequency of an ultrasound burst throughout emission, propagation inside the fluid and reception is discarded. This allows one to account for the effect of frequency-dependent attenuation, which is present in most relevant fluids. By explicitly estimating the center frequency of the received signal, a systematic velocity error stemming from the relationship in Eq. (1) $v\propto \mathrm{1}/{\stackrel{\mathrm{‾}}{f}}_{\mathrm{tx}}$ is avoided.

2. An information loss occurs if only a narrow-band part of a broadband echo signal is processed. Hence a better estimation of the velocity is achieved by including a larger part of the signal spectrum.

Both aspects are addressed by increasing the dimensionality of Kasai's autocorrelation; instead of just correlating along the slow-time axis, a 2D autocorrelation along the slow- and fast-time axis is performed. An autocorrelation with a lag of one fast-time sample yields the estimate of the center frequency as follows:

$$\begin{array}{}\text{(14)}& {\displaystyle }\begin{array}{rl}{f}_{\mathrm{rx}}& \approx f_{\mathrm{s}}\text{arg}\left(\sum _{m}R(\mathrm{\Delta }k=\mathrm{1},\mathrm{\Delta }{n}_{\mathrm{b}}=\mathrm{0})\right)\\ & \approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}{\displaystyle \frac{f_{\mathrm{s}}}{n_{\mathrm{sub}}}}\left(\mathrm{2}\mathit{\pi }\lfloor \mathrm{1}/\mathrm{2}+n_{\mathrm{sub}}{f}_{\mathrm{0}}/f_{\mathrm{s}}\rfloor +\text{arg}{R}^{\prime }(\mathrm{1},\mathrm{0})\right)\end{array}\text{(15)}& {\displaystyle }{f}_{\mathrm{rx}}\in ({\displaystyle \frac{f_{\mathrm{s}}}{n_{\mathrm{sub}}}}\lfloor \mathrm{1}/\mathrm{2}+n_{\mathrm{sub}}{f}_{\mathrm{0}}/f_{\mathrm{s}}\rfloor \pm {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{f_{\mathrm{s}}}{n_{\mathrm{sub}}}}].\end{array}$$

Furthermore, the estimation of the frequencies f_tx and f_d can be performed using M samples per gate, as follows:

$$\begin{array}{}\text{(16)}& f_{\mathrm{d}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}{f}_{\mathrm{PR}}\arg R^{\prime }(\mathrm{0},\mathrm{1}).\end{array}$$

The extension of the Kasai autocorrelation algorithm potentially improves the estimation performance while still preserving a low computational complexity.

4 Online-capable, FPGA-based signal-processing implementation

In order to provide online capability, the signal processing depicted in Fig. 2 has been realized on an FPGA (NI PXIe-7965R; National Instruments, Austin, Texas, USA). The FPGA communicates with a host PC through a peripheral component interconnect express (PCIe) bus and has the ability to stream data through direct memory access into the main memory of the PC.

The amplified echo signals ($U_{\mathrm{PP},\max }=\mathrm{1}\mathrm{V}$) are digitized through an A/D converter module (NI-5752; National Instruments, Austin, Texas, USA) for n_ch=32 channels at an externally provided sampling rate $\mathrm{32}\mathrm{MHz}>f_{\mathrm{s}}>\mathrm{50}\mathrm{MHz}$ with a quantization of 12 bit. The raw data rate r_ADC at this stage is as follows:

$$\begin{array}{}\text{(17)}& {\displaystyle }{r}_{\mathrm{ADC}}{\displaystyle }=n_{\mathrm{ch}}\cdot n_{\mathrm{sampbytes}}\cdot f_{\mathrm{frame}}\cdot N_{\mathrm{EPP}}\cdot N_{\mathrm{sw}}\cdot N_{\mathrm{gates}}\cdot K.\end{array}$$

Data are processed as signed 16 bit integer (n_sampbytes=2 B) and, for a typical configuration as listed in Table 1, the data rate amounts to 1.2 GBs⁻¹.

This data bandwidth is hardly suitable for continuous streaming to a storage device over a long duration (>1 h) with common PC hardware. Therefore, raw data are only briefly retrieved for debugging purposes or for low frame-rate measurements and are otherwise not transferred to the host.

Table 1 Overview of the parameters of the signal processing.

The signal-processing steps that perform an IQ demodulation (bandpass filtering, Hilbert transform and down sampling) are significantly reduced in their computational complexity by fixing the ratio of the sampling frequency f_s to the ultrasound center frequency f₀ at $f_{\mathrm{tx}}/f_{\mathrm{s}}=\mathrm{1}/\mathrm{4}$. The matched filter can be realized for a sinusoidal transmit signal at f_tx with N_periods periods, assuming f_tx≈f_rx with only trivial filter coefficients c_i, as follows:

$$\begin{array}{}\text{(18)}& c_i=\left\{\begin{array}{ll}\mathrm{1}& i=\mathrm{2}+\mathrm{4}n\\ \mathrm{0}& i=\mathrm{1}+\mathrm{4}n,i=\mathrm{3}+\mathrm{4}n\\ -\mathrm{1}& i=\mathrm{4}n\end{array}\right.;n\in [\mathrm{0},N_{\mathrm{periods}}].\end{array}$$

This allows one to implement the filtering without multiplication operations, only negations and additions are needed.

To provide a low computational complexity approximation of the Hilbert transform for a narrowband case, a fixed time delay can be employed (Kantz et al., 2012) as follows:

$$\begin{array}{}\text{(19)}& \widehat{z^{\prime }}(k,n_{\mathrm{b}})\approx z^{\prime }(k-\mathrm{1},n_{\mathrm{b}}),\end{array}$$

where $f_{\mathrm{tx}}/f_{\mathrm{s}}=\mathrm{1}/\mathrm{4}$. The signal processing up to this point contains just the summation, negation and storage primitives and therefore can be implemented on an FPGA with modest resources. The data rate r_IQ at this stage for a typical configuration is given by the following:

$$\begin{array}{}\text{(20)}& {\displaystyle }{r}_{\mathrm{IQ}}=r_{\mathrm{ADC}}\cdot \mathrm{2}\cdot \mathrm{1}/n_{\mathrm{sub}}.\end{array}$$

Through the data reduction of 6.5:1, the data rate at this stage is $r_{\mathrm{IQ}}=\mathrm{185}{\mathrm{MBs}}^{-\mathrm{1}}$ for a typical configuration, as listed in Table 1. A continuous data streaming to a storage device can be sustained for a long duration at this rate.

5 Performance evaluation of narrow-band signal-processing algorithms

5.1 Theoretical limit of measurement uncertainty

In order to characterize the performance of a signal-processing algorithm, it is not only helpful to have relative data compared to other algorithms but also to relate it to a fundamental limit of attainable precision. This absolute limit of uncertainty can be provided by means of the estimation theory using the Cramér–Rao bound (CRB; Radhakrishna Rao, 1945; Cramér, 1946). Given a suitable signal model, the CRB represents the lowest possible variance for estimating a parameter from the signal with an unbiased estimator. In the following, a simple signal model for a discrete time idealized ultrasound echo is described and a derivation of the CRB for velocity estimation is given.

A simple approximation of the ultrasound echo signal realizations $x(k,n_{\mathrm{b}},\mathit{\theta },{\mathit{\sigma }}_n^{\mathrm{2}})$ consists of a sinusoidal signal $s(k,n_{\mathrm{b}},\mathit{\theta })$ superimposed with additive white Gaussian noise $n(k,n_{\mathrm{b}},{\mathit{\sigma }}_n^{\mathrm{2}})$ sparsely and periodically sampled in the fast- (k) and slow-time (n_b) axis as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle }x(k,n_{\mathrm{b}},\mathit{\theta },{\mathit{\sigma }}_n^{\mathrm{2}})=s(k,n_{\mathrm{b}},\mathit{\theta })+n({\mathit{\sigma }}_n^{\mathrm{2}},k,n_{\mathrm{b}}),\text{(22)}& {\displaystyle }\begin{array}{rl}& \text{with}\\ & s(k,n_{\mathrm{b}},\mathit{\theta })=A\cos \left(\mathrm{2}\mathit{\pi }\right(f_{\mathrm{0}}+f_d)\left({\displaystyle \frac{k}{f_s}}+{\displaystyle \frac{n_{\mathrm{b}}}{f_{\mathrm{PR}}}}\right)+{\mathit{\phi }}_{\mathrm{0}}),\end{array}\end{array}$$

and A being the amplitude of the scattering particles' echo, φ₀ a constant phase, and $n(k,n_{\mathrm{b}},{\mathit{\sigma }}_n^{\mathrm{2}})$ Gaussian white noise, with a variance ${\mathit{\sigma }}_n^{\mathrm{2}}$ and zero mean.

The unknown quantities are as follows:

$$\begin{array}{}\text{(23)}& \mathit{\theta }=\left(\begin{array}{l}A\\ f_d\\ {\mathit{\phi }}_{\mathrm{0}}\end{array}\right).\end{array}$$

The CRB provides the lower boundary for the variance of an estimator ${\widehat{\mathit{\theta }}}_i$ according to the inequality, as follows:

$$\begin{array}{}\text{(24)}& \mathrm{var}\left({\widehat{\mathit{\theta }}}_i\right)\ge \text{CRB}\left({\widehat{\mathit{\theta }}}_i\right)={\left[{\mathbf{I}}^{-\mathrm{1}}\left(\mathit{\theta }\right)\right]}_{ii},\end{array}$$

with I(θ) being the Fisher information matrix, as follows:

$$\begin{array}{}\text{(25)}& {\left[\mathbf{I}\left(\mathit{\theta }\right)\right]}_{ij}=-\mathrm{E}\left[{\displaystyle \frac{{\mathit{\delta }}^{\mathrm{2}}\ln \left(p\right(x,\mathit{\theta }\left)\right)}{\mathit{\delta }{\mathit{\theta }}_i\mathit{\delta }{\mathit{\theta }}_j}}\right].\end{array}$$

Kay (1993) provided a formula for the case when the probability density function p(x, θ) of the signal model x, Eq. (22), is a Gaussian joint probability function as follows:

$$\begin{array}{}\text{(26)}& {\left[\mathbf{I}\left(\mathit{\theta }\right)\right]}_{ij}={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_n^{\mathrm{2}}}}\sum _k\sum _{n_{\mathrm{b}}}{\displaystyle \frac{\mathit{\delta }s(k,n_{\mathrm{b}},\mathit{\theta })}{\mathit{\delta }{\mathit{\theta }}_i}}{\displaystyle \frac{\mathit{\delta }s(k,n_{\mathrm{b}},\mathit{\theta })}{\mathit{\delta }{\mathit{\theta }}_j}}.\end{array}$$

The differentiation of $s(k,n_{\mathrm{b}},\mathit{\theta })$ with respect to the unknown quantities is performed analytically, while the matrix inversion was performed numerically using MATLAB (The MathWorks, Inc., Natick, Massachusetts, USA). The resulting CRB for the velocity uncertainty as a function of the signal-to-noise ratio (SNR) is given in Fig. 4c–d. It has a slope of −20 dB/decade, which is consistent with the CRB of other Doppler-based signal-processing problems (Fischer et al., 2010; Chan et al., 2012; Demirli and Saniie, 2001).

5.2 UADV measurements on a reference experiment

For an experimental characterization of the measurement performance of the UADV system, a test rig based on the linear translation of a single scattering object is used (Fig. 3). It consists of a linear stage (41.121.102E; OWIS GmbH, Staufen, Germany) that is mounted over a glass tank with the dimensions of $\mathrm{212}\times \mathrm{81}\times \mathrm{135}{\mathrm{mm}}^{\mathrm{3}}$. It moves a scattering object (glass fiber with a spherical tip, and diameter of 0.6 mm, mounted in a hollow needle) with a constant velocity through water ($\mathit{\vartheta }=\mathrm{20}{}^{\circ }\mathrm{C}$; $c=\mathrm{1480}{\mathrm{ms}}^{-\mathrm{1}}$). The ultrasound sensor array is mounted on the front wall of the tank and therefore insonates through an 8 mm glass wall and a water-based ultrasound couplant.

Figure 3 A measurement setup in which the tip of a glass fiber is mounted on a needle (ND) is insonified by an ultrasound transducer (US) and moved by a linear translation stage (LS). A laser vibrometer (LV) measures its velocity (v_ref) and position (s_ref).

In order to trace back the measurement results of the UADV to the definitions of the respective units in the SI system, a simultaneous measurement of the relative position and velocity was done with a vibrometer (OFV-503; Polytech, Waldbronn, Germany; displacement decoder DD-900 and velocity decoder VD-09). A retroreflective tape (3M Scotchlite) was attached to the shaft of the scattering object's mount. For a velocity set point of 10 mms⁻¹, a standard deviation of the velocity ${\mathit{\sigma }}_{\mathrm{v},\mathrm{ref},\mathrm{rel}}=\mathrm{0.178}\mathrm{\%}$ was determined for the linear stage–vibrometer combination (for the same averaging time as the UADV system).

Table 2 Overview of the ultrasound parameters and the signal-processing algorithms.

A total of 130 measurement cycles were conducted, consisting of a constant translation away from the front wall of the tank with a velocity set-point $v_{\mathrm{ref}}=\mathrm{10}{\mathrm{mms}}^{-\mathrm{1}}$ and the respective backward motion. Of the continuously obtained UADV measurements, only those that originate from two defined positions near to and far from the wall during the movement away from the ultrasound transducer (Richter Sensor and Transducer Technology, Germany) are selected in the postprocessing. The clutter-to-signal ratio (CSR) is ${\text{CSR}}_{\mathrm{1}}=-\mathrm{7.3}\mathrm{dB}$ and ${\text{CSR}}_{\mathrm{2}}=-\mathrm{19.0}\mathrm{dB}$, respectively. To ensure a common time base for vibrometer and UADV measurements, the trigger signal of the UADV is acquired simultaneously with the velocity and position signals. To test the performance under different SNR conditions, white Gaussian noise was added to the raw digitized signals to achieve $\text{SNR}=-\mathrm{6},-\mathrm{3},\mathrm{\dots },\mathrm{12}\mathrm{dB}$. Four algorithmic variants were compared as follows:

• (DEF) – the 1D Kasai velocity estimator without clutter filtering, as described in Sect. 3.4

• (CRF) – the 1D Kasai velocity estimator with a clutter filtering according to Sect. 3.3

• (CRF 2D) – the 2D velocity estimator as described in Sect. 3.5 with clutter filtering but without an estimation of f_tx

• (CRF 2D RF) – the 2D velocity estimator as with clutter filtering including the estimation of f_tx.

The parameterization of the experiment and of the algorithms is listed in Table 2.

Figure 4 Relative systematic deviation (a, b) and relative standard deviation (c, d) of the velocity versus SNR for reference measurements far from the wall (a, c) and near to the wall (b, d); the relative systematic deviation of (DEF) (b) is outside of the axis, with $\mathrm{\Delta }v/v_{\mathrm{ref}}<-\mathrm{42}\mathrm{\%}$; the error bars denote the 95 % confidence interval from 130 measurement cycles.

Figure 5 Example of a flow image of magnetically stirred GaInSn in the central horizontal plane of a cubic vessel. Panel (a) depicts the experimental setup, (b) the mean flow velocity along the d axis, and (c) the standard deviation.

Morley et al., 2008Proffit and Carome (1962)(Proffit and Carome, 1962; Popel et al., 2005)(Galindo et al., 2017)

Table 3 Measurement uncertainty budget for typical MHD experiments in liquid GaInSn.

Figure 4 shows the relative systematic deviation from the reference velocity and the relative velocity standard deviation of the tested algorithms. For the low-CSR case (far from the wall) at SNR=12 dB, it can be seen that a slight negative bias of (DEF) is turned into a positive bias through clutter filtering (CRF) and (CRF 2D). This is compensated by the frequency estimation of (CRF 2D RF), which shows the lowest deviations of all variants for SNR≥3 dB. For the high-CSR case (near to the wall), the variant (DEF) without clutter filter has a relative deviation of $\mathrm{\Delta }v/v_{\mathrm{ref}}<-\mathrm{42}\mathit{\%}$. Through clutter filtering this strong negative bias is turned into a positive bias, which increases with lower SNR for (CRF) and (CRF 2D). The RF estimation of (CRF 2D RF) gives the lowest systematic bias for SNR≥0 dB. The relative standard deviations of all variants of the Kasai's algorithm do not reach the CRB for the given signal model, which is consistent with the findings of Chan et al. (2012). The lowest standard deviations are consistently provided by the variants (CRF 2D) and (CRF 2D RF), which come as close as a factor of 3 to the CRB by using more samples per gate than (DEF) and (CRF). For the given experimental data, the algorithm variant (CRF 2D RF) provides a suitable trade-off between systematic and standard deviation and computational complexity.

6 Measurement uncertainty budget of the UADV in liquid metal

A measurement uncertainty budget according to the GUM (JCGM, 2008) is used to assess the contributions of measurement uncertainty for the UADV system. Based on Eq. (1), the measurand v is derived from the quantities f_d, c and f_rx. Furthermore, the direct influence of the spatial averaging over the flow within the ultrasound beam width is considered. In Table 3, the uncertainties contribution of these quantities are given for a typical MHD experiment.

For the uncertainties of f_d and f_rx, the results of Sect. 5.2 are transferred from the reference experiments in water to typical measurement conditions in low-melting liquid metals. The maximum relative systematic deviation and standard deviation for both investigated CSR and a typical SNR of SNR=5 dB are used to calculate the equivalent uncertainty of the velocity. The influence of an uncertainty in the fluid's speed of sound, c, is estimated by the uncertainty of the measurement of this quantity, in the literature and the temperature dependence, assuming a temperature gradient of $\mathrm{\Delta }\mathit{\vartheta }=\pm \mathrm{10}\mathrm{K}$. The uncertainty arising from spatial averaging through the ultrasound beam characteristics is calculated by assuming a lateral averaging of Δx=3 mm and velocity gradients of numerical simulations of typical MHD experiments (Galindo et al., 2017).

It can be seen that the biggest contribution to the velocity uncertainty of the UADV measurement system for typical MHD settings with ${\mathit{\sigma }}_{\mathrm{v},\mathrm{rel}}=\mathrm{13.9}\mathit{\%}$ stems from the spatial averaging over lateral resolution given by the ultrasound beam width of the unfocused transducers. This provides the most promising starting point for further improvements regarding the measurement uncertainty of the UADV system. Furthermore, it justifies the approximations taken for computationally efficiently implementing the signal processing, even though lower uncertainty algorithms exist that approach the CRB (Chan et al., 2012) because signal processing is not the limiting factor in the measurement uncertainty budget.

7 Example of liquid metal flow imaging

To demonstrate the capabilities of the ultrasound array Doppler velocimeter (UADV) with the proposed signal processing, it is applied to a simple MHD experiment. A cubic vessel with the dimensions of $\mathrm{67}\times \mathrm{67}\times \mathrm{67}{\mathrm{mm}}^{\mathrm{3}}$ is filled with GaInSn and a 25-element linear transducer array (Richter Sensor and Transducer Technology, Germany) is attached to insonify the central horizontal plane (cf. Fig. 5a). With the application of a horizontally counterclockwise rotating magnet field, a counterclockwise central vortex forms. The UADV measures the velocity component along the axis of the transducers (d axis) with the parameterization given in Table 2 and with f_PR=200 Hz. The resulting planar flow image, using signal-processing variant (CRF 2D RF), is shown in Fig. 5b and c.

8 Conclusions

Experimental research in the field of MHD can benefit from online, noninvasive flow imaging for investigating fundamental phenomena, such as flow instabilities and optimizing industrial processes. We describe an online-capable signal processing for pulsed-wave Doppler velocimetry that is tailored to the specific requirements of lab-scale model experiments. It is based on a 2D autocorrelator, which allows for a reduction of systematic and stochastic errors through explicitly estimating the RF and utilizing multiple samples per gate. We optimized the signal processing for low computational complexity and implemented substantial parts on an FPGA. A typical reduction of the data bandwidth of 6.5:1 enables continuous data streaming to PC hardware.

We evaluated the performance of the implemented signal processing in a water test rig with a single scattering object and a reference velocity obtained through a laser vibrometer. Two different clutter signal levels emulate a measurement close to and far from a wall. A velocity standard deviation of ${\mathit{\sigma }}_{\mathrm{v},\mathrm{rel}}=\mathrm{0.4}\mathrm{\%}$ was found, which is about 3 times the fundamental limit of the uncertainty, the CRB, for velocity estimation. The systematic deviation is $\mathrm{\Delta }{v}_{\mathrm{rel}}=\pm \mathrm{0.08}\mathrm{\%}$.

We investigated the measurement uncertainty budget for flow velocity measurements in a typical MHD experimental setup for the low-melting alloy GaInSn. The total measurement uncertainty of ${\mathit{\sigma }}_{\mathrm{v},\mathrm{rel}}=\mathrm{13.9}\mathrm{\%}$ almost solely stems from the effect of spatial averaging over the lateral resolution for flows with high-velocity gradients. This justifies the approximations taken for lowering the computational complexity of the signal processing.

A measurement uncertainty budget of a typical MHD experiment at laboratory scale suggests improvements towards a better lateral resolution. In the context of flow imaging, this can be provided by the focusing and steering of the ultrasound beam using the phased-array principle.

The presented signal processing enables online, multiplane flow visualization with the UADV research platform. A long measurement duration (>1 h), combined with a high frame rate (>10 Hz), allows one to investigate complex, instationary flows such as instability phenomena in cubes.