{'fraction_non_alphanumeric': 0.02966091640367321, 'fraction_numerical': 0.009871152979142444, 'mean_word_length': 4.965180467091295, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Nonlinear Beamforming for Enhanced Spatiotemporal Sensitivity in High Frame Rate Ultrasound Flow Imaging

Typically, an ultrasound flow imaging system employs the conventional delay and sum (DAS) beamformer due to its inherent low complexity. But the conventional DAS technique offers poor contrast, low imaging resolution, and limited spatiotemporal sensitivity. This article attempts to improve the spatiotemporal sensitivity of conventional flow imaging with a novel multiply and sum based nonlinear high-resolution (NLHR) beamforming approach. The major advantages of the proposed beamformer are the harmonic generation and the enhanced coherence in beamformed signals that improve the spatiotemporal sensitivity towards flow transients. We demonstrate the proposed beamformer for a directional cross-correlation as well as an autocorrelation based velocity estimator with simulated parabolic flow profiles of different velocities and flow directions, an in-vitro rotating disk dataset, and pulsatile flow experiments. The sensitivity of NLHR beamforming towards the flow transients is validated in-vitro with a sudden reversal of flow direction and air bubble tracking experiments. The comparison between the time-frequency plots of DAS and NLHR beamforming indicates that the impulsive spatiotemporal changes induced by the flow of air bubbles are clearly characterized by nonlinear beamforming than that of DAS beamforming. Furthermore, better spatiotemporal velocity tracking of a single air bubble and a clear distinguishability between the tracking of two proximal air bubbles are observed in-vitro. Preliminary studies on the in-vivo carotid data also show comparable, if not better, results than that of the DAS algorithm. Detailed results for each test case in simulation, phantom, and in-vivo studies are available as movies with the supplementary material and online.


Introduction
Ultrasound Doppler technique is the most widely used non-invasive technique for flow imaging and has a great diagnostic value in investigating vascular hemodynamics [1]. It uses high frequency non-ionizing sound waves to create multiple insonification of blood vessels which makes it the safest technique to quantify the complex flow dynamics. Traditional Doppler ultrasound-based flow imaging techniques are highly dependent on the beam-to-flow angle and least sensitive to the transverse flow component [1,2,3]. This makes it difficult to quantify the complex flows associated with the superficial vasculature as they are immediately beneath the skin surface and the flow is purely transverse in most cases.
The directional beamforming technique in [8] and [9] focuses the received signals along the flow direction and employs cross-correlation to estimate the velocity magnitude. But this approach requires the knowledge of the flow angle with respect to the emitted beam for beamforming. In this regard, correlation and numerical triangulationbased techniques have been developed to estimate the flow angle [20,21]. However, synthesizing the directional signals by rotating the grid would be difficult in the case of complex flows in which the flow is not in a welldefined direction. To address these limitations, vector triangulation-based flow imaging approaches are proposed in the literature [14,15,16,17,18]. However, the contrast, resolution and accuracy of the velocity estimation are mainly determined by the beamforming algorithm employed in the system.
Delay and sum (DAS) is the most commonly used beamforming technique in medical ultrasound imaging due to its inherent low complexity [22,23]. Conventional DAS beamformer in ultrasound systems involves delaying the received signals and then summing it up to obtain the dynamically focused receive beams. Usually, the beamformer uses an apodization window to weight the delay compensated signal before summing. The apodization enhances the energy of the signal components from the desired grid point while suppressing the signal components from the undesired locations [24]. This effectively improves the signal-to-noise ratio and provides a better beam at the output of the beamformer [23]. But the conventional DAS technique employed in most of the aforementioned flow imaging methods offers poor contrast, low imaging resolution, and limited spatiotemporal sensitivity. In this paper, the term sensitivity is intended as the capability to detect the flow velocity changes and the transient events.
Improving the spatial resolving capability of VFI helps in understanding more localized flow patterns and velocity distribution in various pathological conditions where VFI is used for diagnosis and monitoring. Enhanced sensitivity in VFI towards highly dynamic and impulse components in velocity aids in identifying new biomarkers or pathological features in medical diagnostics. Moreover, the importance of enhancing the sensitivity towards impulse velocity components and higher temporal resolution (reduced time window to detect velocity changes) is further emphasized in various applications such as functional ultrasound imaging of the brain [25], microbubble localization and tracking, remote drug delivery, neovascularization imaging [26,27].

Related Work
One of the early attempts towards improving the sensitivity of conventional Doppler has been the coherent plane wave compounding (CPWC) at high-frame rates as reported in [26,28]. The CPWC technique employs multiple steered plane wave transmit and coherently combines the beamformed data to obtain a high-resolution image and to estimate the velocity [28,29]. However, this technique was predominantly employed for power Doppler imaging to increase the sensitivity of axial velocity components. The influence of motion in CPWC has been studied by Denarie, Bastien et al., and a compensation scheme for radial motion is reported in [30], but the effect of multi-directional velocity components is not addressed. Moreover, with the increase in the number of steered plane wave transmit, the effective frame rate reduces and hence the CPWC based approach would be ideal for the investigation of low velocity flows as shown in [31] and [32].
Later, singular value decomposition (SVD) based filters were introduced to increase the spatiotemporal sensitivity in high frame rate Doppler imaging [33,34]. Subsequently, principal component analysis based approaches have been demonstrated in-vivo for vascularization imaging [35] and surgical management of open-brain tumors [36]. Recently, machine learning and deep learning based approaches are emerging as automated ultrasound image analysis tools and have shown promising results for various diagnostic imaging tasks like detection, segmentation, and classification of tumors and lesions [37,38,39]. The convolutional neural networks, despite being successful in optical flow estimations, have not been exploited much for ultrasound blood flow imaging and very limited efforts have been reported in the literature [40,41,42]. Moreover, most of the efforts in flow imaging have been coupled with ultrasound contrast agents and are focused on obtaining the vascular maps at a high spatial resolution but don't contribute much towards improving the temporal sensitivity or the temporal super-resolution [43,44,45]. However, it has been reported that the deep learning models could provide better results by incorporating temporal information of the microbubbles [45].
Towards improving the temporal resolution, a minimum variance based adaptive velocity estimation technique has been recently reported in [46], but with higher computational complexity. However, the efforts towards employing nonlinear beamformers of lower computational complexity (than minimum variance) like delay multiply and sum (DMAS) technique for flow imaging is limited. Delay multiply and sum, commonly known as DMAS algorithm is a nonlinear beamforming technique in B-mode ultrasound imaging which has a higher spatial and contrast resolution with a better clutter and noise rejection when compared to DAS [24,47]. In one of the recent studies reported in [47] the DMAS algorithm is modified to a filtered delay multiply and sum (F-DMAS) by employing an additional bandpass filter for Brightness-mode (B-mode) imaging. The simulation and the in-vivo studies with F-DMAS algorithm have demonstrated a better-quality image with a higher dynamic range and improved contrast resolution at the output [48]. This nonlinear technique is further investigated for plane wave and synthetic aperture high frame rate ultrasound imaging [48,49] as well as for photoacoustic imaging [50,51]. Similar to the formulation of DMAS, spatial coherence based techniques have been proposed for better clutter suppression in power Doppler imaging [52,53,54]. In this regard, the quality improvement in B-mode and power Doppler imaging with DMAS based beamforming [48] could be leveraged to improve the resolution and sensitivity of VFI techniques.

Objective and Contributions
Inspired by the F-DMAS beamforming in B-mode imaging, this paper attempts to address the spatiotemporal sensitivity of the conventional flow imaging techniques with a novel non-linear beamforming approach without the use of any contrast agents and deep learning based methods. The major contributions of this work are as follows.
1. A novel nonlinear beamforming technique is proposed for ultrasound flow imaging and to the best of the authors' knowledge, this is the first effort towards the application of nonlinear beamforming in flow imaging. 2. We validate the proposed nonlinear beamformer using typical parabolic flow simulations with a cross-correlation based velocity estimator and an autocorrelation based velocity estimator which are the common velocity estimation techniques in the literature. 3. Further, the proposed approach has been thoroughly investigated for velocity sensitivity with in-vitro datasets including a rotating disk, air bubble tracking, and flow direction reversal. Finally, we report the in-vivo performance evaluation of the proposed approach for a typical pulsatile flow in a carotid artery dataset. 4. The results, when compared to the state-of-the-art DAS based flow imaging approaches, suggest that the proposed method provides better spatiotemporal sensitivity towards the flow transients.
The rest of the article is organized as follows. The following section (Section 2) explains the proposed nonlinear high-resolution beamforming technique for flow imaging and describes how it can be used in conjunction with triangulation with auto-correlation estimator and directional cross correlation estimator. Section 3 discusses the implementation details along with the simulation and experimental setup. Section 4 presents the results obtained from the simulations and experimental investigations in comparison with DAS based VFI methods and the discussion of results is presented in Section 5. The article is concluded in Section 6 with a summary of the work.

Proposed Approach
A schematic of the proposed non-linear high resolution (NLHR) beamformer for flow imaging is shown in Fig.  1(a). A detailed description of each step in the approach is provided in the subsequent sections.

Nonlinear High-Resolution Beamformer
The first part of the NLHR beamforming deals with dynamic delay compensation of the received signals to enable channel-wise spatial focusing. The second part consists of the design of an appropriate apodization function to form the directional weighted focused beams at left and right receive sub-apertures. Further, these beams are multiplied pairwise and summed before it is filtered to obtain the higher harmonic beamformed signals in the final step. These beams are used to estimate the velocity magnitude and direction. It must be noted that these beams will be having higher lateral resolution and better contrast (clutter suppression) as discussed in [47]. In this work, non-steered plane waves are used to insonify the region of interest at very high frame rates to aid the improvement in velocity sensitivity. However, the proposed approach is also applicable to multi-angle plane waves, synthetic aperture and diverging wave transmit schemes. A detailed explanation of the approach is provided in the following subsections.

Channel Directive Beam Synthesis
The received radio frequency (RF) echo signals, ( ) (where represents the ℎ transducer), are delay compensated and mapped to a two-dimensional grid in the region of interest. This process is referred to as time to space mapping of the RF signals. The delay compensated signal, ( ) for a given pixel p is obtained as , , where, , is given by (1): where , , , are the transducer and the pixel coordinates (corresponding to pixel ) respectively. It is assumed that a linear array transducer is placed at the start depth (i.e., at = 0) and is the speed of sound. The size of the delay compensated signals for the full region of interest will be × × , where, is the number of channels, is the number of pixels and is the number of frames. The non-steered plane wave excitation is assumed. The channel directive beams for a pixel, ( ), are formed with each of the transducers as the beam centers according to (2).
where, ( ) is the apodized amplitude of the ℎ transducer for ℎ pixel, , is a Gaussian weight vector, centered at ℎ receive element. The width of the Gaussian function is dynamically sized according to the receive distance (i.e., ℎ = ∕ , where is the F-number), instead of depth as in traditional apodization. Thus, unlike traditional apodization, the channel directive beam synthesis as given by (2) creates beam observations, that are equivalent to the focused observations from different angles for each pixel. The size of the matrix would be × × as shown in Fig. 1 and obtained by repeating (2) for pixels and frames. This step converts the delay compensated channel data to the spatially focused channel beams and its significance in flow imaging has been evaluated in our previous work [55].

Sub-Apertures Formation
The next step is to apply suitable weight onto the channel directive beams to obtain multiple sub-apertures to the left and right of the selected transmit element following the principles of triangulation. The formation of sub-apertures is a traditional step involved in any the triangulation-based velocity estimation scheme [7,13,14,18]. These are referred to as left and right sub-aperture signals, , ( ) and , ( ), corresponding to different transmit-receive angles ( ) and are given by (3) and (4): where, , and , are the left and right sub-aperture apodization weight vector of length . For each pixel, a Gaussian weight vector centered at the receive element that corresponds to the transmit-receive angle ( ) under consideration is used. This provides maximum weightage to the beam received at an angle and thus ensures the directionality of the left and right receive sub-aperture signals corresponding to each pixel in the region of interest. Unlike in the traditional DAS beamformer, the weighted signals are not directly summed in our approach and hence the signals obtained with (3) and (4) are regarded as -weighted beam focused left and right sub-aperture signals. It should be noted that these signals are weighted in the direction defined by the and not the flow direction. This sub-aperture formation at receive enables the beamforming to be employed for most of the split aperture vector flow imaging techniques.

Multiply and Sum Operation
Once the -weighted beam focused left and right sub-aperture signals are formed, each beam is multiplied with every other beam and then summed up to obtain the equivalent DMAS beamformed signal in VFI,̂ ( ) and̂ ( ), as given by (5) and (6).
It can be observed that the multiply and sum operation is mathematically a cross-correlation function between all the receive channels in which the self-product terms ( = ) are excluded. Also, it is to be noted that either is considered in the summation of multiply and sum operation as both the coefficients are equal in magnitude and hence the information is redundant. This will halve the corresponding correlation coefficient and reduce the number of required multiplications by a factor of 2 [47,56]. The total number of multiplications required for an N channel receive aperture, , is given by (7).

Band Pass Filtering
The frequency spectrum of the delay compensated directional focused beams of each channel will be a band centered around the center frequency of the insonified plane wave, 0 . But the frequency spectrum of multiply and sum beamformed signals,̂ ( ) and̂ ( ), will have two peaks, one centered at 0 − 0 ≈ 0 and the other at 0 + 0 ≈ 2 0 because of the fact that (5) and (6) multiply two signals having almost similar frequency [47]. Now, a band pass filter with a passband centered around 2 0 is used in the final step to remove the DC component. Thus, the frequency spectrum of filtered left and right beamformed signals,̃ ( ) and̃ ( ), will be centered at 2 0 . In any Doppler technique, the higher the transmit frequency, the higher the Doppler shift and the better the velocity estimate. In this work, since the harmonics induced by multiply and sum is employed, it is expected to give a higher Doppler resolution.

Velocity Estimation
A truncated SVD based spatiotemporal clutter rejection filter [33] (as available in the Matlab UltraSound Toolbox [57]) is used to remove the wall motion. The details of the SVD filter are provided in the results section. As shown in Fig. 1(a), the flow velocity vectors can be derived either by employing an autocorrelation-based estimator or a cross-correlation based estimator. Both the methods are described briefly in the following subsections.

Triangulation with Autocorrelation (TAC)
The mean slow-time frequency of the clutter suppressed left and right beamformed signals, ( ) and ( ), is determined using the classical lag-one autocorrelation approach [58]. A window size of is used to estimate the autocorrelation. Once the slow-time frequency, and , are estimated, the orthogonal components of the velocity can be estimated using the following equations (8) and (9) as in conventional triangulation (9) [7] and that used in [17]: where ′ = 0 in the case of DAS based beamforming and ′ = 2 0 in the case of DMAS based beamforming because the center frequency of the DMAS based beamformed signal is around 2 0 . Unlike in [17] where only one pair of observations is considered for triangulation, the multiple sub-aperture formation in the proposed approach provides multiple observations enabled by multiple at receive. The velocity magnitude and flow direction are calculated as the mean of all such observations.

Directional Cross-Correlation (DCC)
Directional cross-correlation based velocity estimation as in [8] is performed for the observations enabled by multiple values at receive and the mean estimates among them are computed. However, to perform DCC, beamforming must be done along the flow direction. This is done by rotating the flow grid and the time to space mapping of the received RF signals is done for the rotated spatial grid. Apart from this, no changes are required in the beamformer architecture to perform DCC based velocity estimation. A correlation interval of −10 ∶ 10 (equivalently, a cross-correlation window size of = 20 ) and a sampling interval of 0.1 are chosen for DCC, as in [8].

Implementation and Validation
In this section, the details on the implementation of the proposed approach in MATLAB® as well as the simulation and experimental setup used for the evaluation and validation of the proposed approach are discussed.

Considerations to Avoid Aliasing
With the multiply and sum operation discussed in Section 2.1.3, the frequency spectrum of the beamformed signal would get shifted to 2 0 and the Doppler shift would also get doubled. To ensure no aliasing due to the frequency doubling, the received RF data is resampled along the axial direction as well in the temporal (frame) direction immediately after the acquisition to satisfy the Shannon-Nyquist sampling theorem. The size × × in Fig. 1 represents the size of the received RF signals after resampling. It is assumed that the acquired RF data is not aliased and has a sufficient frame rate to capture the flow velocities as in any Doppler system acquisition.

Simulation Setup
The performance of the proposed approach is initially evaluated and validated with extensive simulations using Field II [59,60] and MATLAB R2019a (The MathWorks Inc., Natick, MA, USA). Flow phantoms are generated using the Field II program (version 3.24) with the parameters shown in Table 1. A standard linear array probe with 128 active elements is modeled in the Field II simulator and no transmit apodization is applied for simulation studies. A 5-cycle sinusoid at 8 is used for transducer excitation and a Hanning weighted single cycle sinusoid is used as the transducer impulse response at the transmit and receive. Non-steered plane wave transmit is used for all the datasets in this paper. The simulation parameters closely follow the parameters chosen in [8].
The performance of the NLHR beamforming with TAC and DCC based velocity estimators is compared to the true/theoretical profiles and the corresponding velocity estimators with DAS beamforming. The DAS beamforming for the TAC velocity estimator adopts the single transmit dual-angle multi receive scheme (which is proposed in one of our recent works [18]). The DAS beamforming for the DCC velocity estimator adopts the directional beamforming scheme reported in [8]. To have a fair comparison, the transmit-receive angles, correlation interval, and other parameters are kept the same for DAS and NLHR beamforming. The resampled RF data corresponding to = 20 and = 0.8 is used for DAS as well as NLHR beamforming unless otherwise mentioned. It should be noted that the comparison of NLHR and DAS beamforming for DCC is only done for magnitude estimates as DCC estimates only provide the velocity magnitude. Three different test cases were investigated with simulations for the performance evaluation. Test case 1: To verify and validate the accuracy in velocity estimation, typical parabolic flow profiles, as given by (10) [8], having the peak velocities, = 5, 10, 25, and 50 ∕ is simulated for a transverse vessel of radius, = 5 with its axis at 25 depth from the probe surface.
where r is the radial position in the vessel. The velocity magnitude and angles are estimated with TAC and DCC for NLHR beamforming and are compared with the corresponding estimates with DAS beamforming and the true profiles. Test case 2: In the second test, datasets are simulated for a 5 radius vessel whose axis is inclined at −20 • , −10 • , 0 • , 10 • , and 20 • with the horizontal and having a peak velocity, = 50 ∕ . The performance comparison is done as in test case 1.
Test case 3: The third test aims to study the influence of correlation intervals in TAC and DCC estimators with DAS and NLHR beamforming. The correlation interval for DCC for comparison is chosen as −10 ∶ 10 and −5 ∶ 5 , denoted by = 20 and 10 , respectively whereas the averaging interval is chosen as = 0.8 and 1.6 . The correlation interval for TAC for comparison is chosen as = 0.8 and 1.6 whereas the averaging interval is chosen as = 20 and 10 . The results for the 4 combinations of the and values are compared for DAS and NLHR beamforming with TAC and DCC estimators.

Experimental Setup
In-vitro Rotating Disk: The feasibility of the proposed approach is verified in-vitro by using a rotating disk (diameter of 2 ) dataset [61,62]. The dataset was acquired using a Verasonics scanner (Verasonics Inc., Kirkland, WA, USA) and a 128-element linear array with non-steered plane wave insonification at a center frequency of 5 and a pulse repetition frequency of 10 . The RF signals for 32 frames with a sampling frequency of 20 were recovered from the bandpass sampled (at a sub-Nyquist sampling rate of 6.66 ) RF dataset before beamforming.

In-vitro Flow Phantom Experiments:
The proposed approach is thoroughly investigated for velocity sensitivity using a commercial Doppler ultrasound flow simulator from CIRS (CIRS Inc., Norfolk, VA, USA). The simulator includes a peristaltic flow Pump (Model 769), Doppler fluid (Model 769DF), and a tissue-mimicking Doppler ultrasound flow phantom (Model 069A) which has a peripheral vessel with a controlled flow velocity. Three test cases with pulsatile flow profiles were investigated in the flow phantom at a pulse repetition frequency of 8 . For the first test case, a typical pulsatile flow with a peak velocity of 50 ∕ is configured. In the second test case, the Doppler fluid is seeded with some air bubbles by gently shaking the fluid and the pulsatile flow with parameters same as above is configured using the flow pump. The velocity sensitivity and spatial localization capability are evaluated by tracking air bubbles over time with the proposed NLHR approach. This aims to test the sensitivity of the proposed approach towards the impulsive changes in the flow induced by the air bubbles. In the third test case, a pulsatile flow with a sudden reversal of direction is imposed during the data acquisition. The results are compared to that of DAS based method.
In-vivo Carotid Artery: Finally, the proposed beamforming technique is tested on an in-vivo carotid dataset captured from a healthy volunteer by following the principles outlined in the Helsinki Declaration of 1975, as revised in 2000. The instantaneous quiver plots and the velocity profiles along with the velocity estimates over time are compared to that of DAS based methods. For all the experimental studies, the TAC estimator is used with NLHR beamforming and DAS beamforming. The parameters used for the in-vitro flow phantom experiments and the in-vivo studies are shown in Table 1.

Performance Metrics
The bias and standard deviation in velocity magnitude ( % and % ) for the simulation and experimental studies, if not reported differently, are determined as: where , and are the measured, the true (or theoretical), and the peak velocity values, respectively; (•) denotes the standard deviation function. The accuracy of angle estimation is determined using the bias ( ) and standard deviation ( ) measures as given by (12): where and are the measured and the true (or theoretical) values, respectively. Neither spatial and temporal averaging nor any bias compensation is applied to the velocity estimates. It should be noted that the DAS beamforming is performed on the temporally and spatially resampled dataset to ensure a fair comparison.

Results
This section discusses the results obtained from the detailed simulation study and the experimental investigations.

Simulation Results
The estimated velocity for a transverse parabolic profile of peak velocity 50 ∕ with the proposed NLHR beamforming for DCC estimator compared with that of DAS beamforming is shown in Fig. 2(a) and the comparison plots for velocity magnitude and angle with TAC estimator are shown in Fig. 2(b) and (c) respectively. The quantitative performance comparison of NLHR beamforming based velocity estimations with that of DAS based methods in terms of bias and standard deviation in the magnitude and angle estimates for test cases 1, 2, and 3 are shown in the plots in Fig. 2 columns 2, 3, and 4 respectively. For the plots shown in Fig. 2 (d)-(l), the inner 90% region of the vessel radius is considered. In general, an underestimation in the velocity magnitude is observed with respect to the true profile for both the beamformers and estimators. The bias and standard deviation in the velocity magnitude and angle estimates with NLHR beamforming are observed to be similar to that of DAS beamforming for TAC based estimates. However, a reduction in the bias is observed for the velocity magnitude estimates with NLHR beamforming for DCC based estimates whereas the standard deviation is found to be similar, if not better.

Phantom Results
In-vitro Rotating Disk: Clutter filter was not required for this dataset [57,61]. The comparison of the color Doppler plots and vector flow images for DAS and NLHR beamforming is shown in Fig. 3(a)-(b) and (c)-(d), respectively. The Comparison plots for test case 3 with DCC (j), and TAC (k, l) estimators. The error bar is centered at median bias and the length of the error bar corresponds to one standard deviation. major observations from this experiment are: 1) the velocity estimates with NLHR beamforming matches with that of DAS beamforming 2) the color Doppler map of NLHR beamforming is slightly sharper than that of DAS and provide better velocity estimates at the disk edges as evident from the disk_1.gif (supplementary material). The improved spatial localization of the beamformed signals due to the coherence enhancement and the improvement in contrast resolution due to the suppression of incoherent echoes with NLHR beamforming could be seen in the disk_2.gif (supplementary material). The experiment, being an evaluation of different axial and lateral velocities spanning all the flow directions (i.e., 0 − 2 radians), confirms the feasibility of nonlinear beamforming in flow imaging.
In-vitro Flow Phantom Experiments: Unlike simulations and the rotating disk dataset, clutter filtering was required for the flow phantom experiments. The top 3 and top 5 singular value signal components were removed for DAS based estimation and NLHR based estimation respectively to suppress the clutter. The choice of these thresholds is in accordance with the power of the singular values observed for the beamformed signals as shown in Fig. 4(a). A correlation interval of 32 (corresponds to 256 frames of acquisition and 512 frames after resampling) is used for velocity estimation with DAS and NLHR based methods. Fig. 4(c)-(h) shows a comparison between DAS and NLHR beamforming for a typical pulsatile flow dataset in which the velocity estimates of DAS and the proposed nonlinear beamforming are observed to be very similar. It is interesting to note that the parabolic nature of the flow is more evident in the NLHR when compared to the slightly flat-topped result in DAS (referring to Fig. 4(e)). Fig. 5(a1)-(a3) and (b1)-(b3) shows the set of frames (captured at three instances during the flow of an air bubble) for DAS and NLHR beamforming, respectively. A spatially well resolved velocity tracking of the air bubble is observed for NLHR when compared to DAS beamforming. The velocity transients due to the air bubble is captured more concisely by the proposed NLHR beamforming than DAS beamforming as evident from the velocity-time plots in Fig. 5(c) and (d). This is further supported by the axial velocity plot that depicts the velocity variations in the presence of an air bubble as shown in Fig. 5(e) and (f). Thus Fig. 5 shows the improved spatial resolution (in terms of localized spatial velocity) and velocity sensitivity (as evident from the localized time variation of velocity).
To investigate further, we consider a scenario in which two proximal air bubbles are moving in the vessel as shown in Fig. 6(a). From the velocity-time trace during the flow of the air bubbles shown in Fig. 6(b), it is interesting to notice that the NLHR based estimates distinguish between the two air bubbles more distinctly. The beamformed slow time signals of the left and right sub-apertures for the correlation window (of 512 samples) during the flow of air bubbles  are shown in Fig. 6(c), (d) for DAS, and in Fig. 6(g), (h) for NLHR beamforming. A clearer separation is observed in the NLHR beamformed left sub-aperture signals as highlighted in Fig. 6(g) when compared to that of DAS in Fig. 6(c). In addition to this, the NLHR beamformed right sub-aperture signals corresponding to the air bubble are observed to be much more confined in time than that of DAS as indicated in Fig. 6(d) and (h). The time-frequency plots of the beamformed slow-time signals in Fig. 6(c), (d), (g) and (h) are shown in Fig. 6(e), (f), (i) and (j). These highly localized time-frequency plots were obtained with the reassigned pseudo Wigner-Ville distribution [63]. The velocity transients, induced by the flow of air bubbles, which are characterized by the excitation of all the frequencies in the spectrum are evident for the NLHR beamforming as shown in Fig. 6(i) and (j) than that of DAS in Fig. 6(e) and (f). Fig. 7(a) and (b) summarizes the result for the test case in which a sudden reversal of the flow direction is imposed. The velocity magnitude estimates were observed to be similar for DAS and NLHR beamforming. However, the angle estimates with the proposed NLHR beamforming appear to be converging with the possible flow angles according to  the vessel orientation in the phantom as shown in Fig. 7(b). The video results for all the test cases are available as supplementary files which represent the inferences quite comprehensively.

Preliminary In-vivo Validation
For the in-vivo dataset, the top 6 and top 11 singular value signal components were removed for DAS based estimation and NLHR based estimation respectively. The choice of these thresholds is again in accordance with the power of the singular values observed for the beamformed signals as shown in Fig. 8(a). We also observe that the singular vectors corresponding to those singular values are associated with the tissue and wall motions that are  characterized by low frequency as shown in Fig. 8(b). The quiver plots obtained with DAS and NLHR beamforming during peak systole are shown in Fig. 8(c) and (d) respectively. The velocity-time trace for three different locations (one in the middle of the vessel and the other two closer to the vessel wall as indicated in Fig. 8(c) and (d)) in the carotid artery is shown in Fig. 8(e) and (f). The video result with more details of the in-vivo study is available as supplementary material (invivo_carotid.mp4).

Discussion
As discussed in Section 1, the limited sensitivity of the conventional flow imaging techniques were addressed either with coherent compounding techniques or an improved clutter filter design. In this paper, we attempt to introduce a novel multiply and sum based nonlinear beamformer in flow imaging that could provide better spatiotemporal sensitivity towards the flow transients. The multiply and sum operation in F-DMAS beamforming is proven to provide enhanced coherence in the beamformed signals with significantly better contrast and lateral resolution in B-mode imaging [47]. Owing to these advantages of multiply and sum operation, the proposed approach for flow imaging has shown better performance in terms of spatial resolution and velocity sensitivity as compared to state-of-the-art DAS based methods. The output of each step involved in the proposed nonlinear beamformer for an in-vitro rotating disk dataset is included as supplementary material.

Better Spatial Resolution and Velocity Contrast
The improvement in the spatial resolution is observed with the spatial distribution of the velocity vectors in Fig. 5, where the air bubble is tracked along with the typical blood flow velocity. This improved velocity contrast is due to the improved contrast resolution which is a characteristic property of DMAS beamformer. In line with the characteristic of F-DMAS in [47], the following characteristics can be observed for the proposed non-linear beamformer as well in flow imaging: 1. The widening of the "synthetic" receive aperture due to the correlation function in (5) and (6) that has 2 − 1 coefficients as in the case of the F-DMAS beamformer in [47] 2. The doubling of the center frequency of the beamformed signal provides twice the number of frequency components than in a DAS beamformed signal (as evident from the singular value plots in Fig. 3(e), (f), Fig.  4(a), (b) and Fig. 8(a), (b)) that could be considered as the improvement in frequency resolution or equivalently the velocity resolution.
Conversely, the signal components in DAS beamformed signal for the highlighted portion in Fig. 4(a) and Fig. 8(a) are arguably not relevant for velocity estimation for two reasons: 1) the frequency of all those components are observed to be the same and 2) the energy of those components are observed to be the same and very minimal. Moreover, this signifies that the closer frequencies (Doppler shifts) are better resolved with NLHR beamforming than DAS and a greater number of spatiotemporal components will be available for the SVD clutter filtering resulting in more efficient clutter suppression when compared to that of DAS. The improvement in the velocity and the spatial resolution may not be significant in typical parabolic flow evaluations as we have observed with the simulations (Fig. 2) and experiments ( Fig. 3 and Fig. 4). However, the improved velocity contrast and spatial resolution are highly appreciated in the presence of flow transients as demonstrated in-vitro in Fig. 5 and Fig. 6.

Enhanced Sensitivity Towards Flow Transients
Flow velocity transients (ideally, impulsive changes in the velocity) are characterized by the excitation of almost all the frequencies in the spectrum of the beamformed signals. The introduction of the air bubbles to the Doppler fluid would effectively induce such velocity transients and it is interesting to note that the proposed beamformer was sensitive enough to capture the impulsive transients as evident from the time-frequency plots in Fig. 6. Moreover, the comparison of beamformed signals (Fig. 6(c) and (g)) and its time-frequency plot suggests that the transients created by each of the air bubbles are better resolved (temporally) with the proposed nonlinear beamforming. This is due to the coherence enhancement and the harmonic generation with themultiply and sum operation. Also, better flow angle estimates with the NLHR beamforming observed during the flow direction reversal (Fig. 7) further signifies the sensitivity of the proposed approach towards the flow transients (both in terms of the magnitude and direction of the flow).

Applications
The proposed technique with improved spatial resolution and velocity contrast would be of great relevance in microbubble localization and tracking, remote drug delivery, neovascularization, and microvasculature imaging applications. Furthermore, the enhanced sensitivity of the proposed approach could probably help in identifying new biomarkers or pathological features in medical diagnostics.

Limitations and Future Scope
The advantages of the proposed approach are at the cost of increased computational complexity when compared to DAS beamforming. So the proposed approach may not be suitable for portable imaging systems in its current form and would require further optimization. To address the computational complexity of the proposed approach, neural network accelerators which have found success in B-mode image reconstruction will be employed to speed up the proposed workflow as an immediate future work. Further, the proposed approach must be investigated with detailed in-vivo animal studies.
Since the major objective of this work is to improve the spatiotemporal sensitivity in high frame rate flow imaging, the proposed technique is demonstrated and validated only for non-steered plane transmit schemes and other transmit schemes like focused, synthetic aperture, and diverging wave transmit were not considered during the formulations of the proposed approach. Hence the proposed nonlinear approach must be re-visited for the transmit schemes other than plane waves. In this regard, further developments of this work also include the analysis of this approach for synthetic aperture and diverging wave vector flow imaging and power Doppler.
With reference to the improvement in the velocity contrast offered by the proposed nonlinear beamformer, as discussed in Section 5.1, the proposed approach can provide a better separation between closer velocities. This would have a significant impact on the evaluation of slow flows and the velocity estimates towards the periphery of blood vessels because such low-velocity magnitudes overlap with that of the tissue motion. In this regard, the scope of this work could be further expanded to a systematic study on the impact of the proposed nonlinear beamforming in the design of clutter filter.

Conclusion
Inspired by the advantages of F-DMAS beamforming in B-mode imaging, this article investigates the performance of a nonlinear beamformer in high frame rate ultrasound flow imaging. Towards this, a novel multiply and sum based nonlinear beamformer is introduced for high frame rate flow imaging and is demonstrated using DCC and TAC based velocity estimators. Further, the proposed approach has been investigated in detail for velocity sensitivity with invitro datasets including a rotating disk, air bubble tracking, and flow direction reversal. The proposed approach has performed similar, if not better, for parabolic flow simulations and typical pulsatile flows in-vitro and in preliminary invivo studies. However, owing to the enhanced coherence and the harmonic generation by the multiply and sum operation in the beamformer, an enhanced velocity contrast, spatial resolution, and sensitivity towards the flow transients were observed in-vitro. The results suggest that the method has the potential to greatly impact the studies on microbubble tracking, cavitation, and other applications in which velocity sensitivity and resolution are of prime interest. However, the improved spatiotemporal sensitivity offered by the proposed beamformer is at the cost of increased computational complexity and hence the proposed approach in its current form may not be suited for portable flow imaging systems that are constrained by power and computational resources. Moreover, in its current form, we have not considered any transmit schemes other than plane waves, and hence more studies are required to investigate the feasibility of the proposed approach for other transmit schemes.