WO2018172882A1 - Doppler ultrasound using sparse emission pattern - Google Patents

Abstract

An apparatus for measuring target velocity includes a transducer (30) and a processor (40). The transducer is configured to transmit into a target a set of pulses, which is sparse but has a recoverable power spectrum, and to receive a signal reflected from the target in response to the set of pulses. The processor is configured to estimate a velocity of a selected region in the target by analyzing the reflected signal.

Description

DOPPLER ULTRASOUND USING SPARSE EMISSION PATTERN

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of U.S. Provisional Patent Application 62/474,164, filed March 21, 2017, and U.S. Provisional Patent Application 62/485,955, filed April 16, 2017, whose disclosures are incorporated herein by reference.

TECHNICAL FIELD

Embodiments described herein relate generally to processing sparse signals, and particularly to methods and systems for efficient ultrasound imaging.

BACKGROUND

Ultrasound imaging is a non-invasive imaging technique, commonly used for various medical purposes. Ultrasound techniques are based on transmitting ultrasound pulses toward the target and analyzing echoes of the transmitted pulses reflected by various objects in the target.

Ultrasound systems typically support both brightness mode (B-mode) imaging and spectral Doppler imaging. B-mode imaging mode enables the physician to visualize a relatively large area such as an entire organ of the patient, whereas the Doppler mode is used for evaluating tissue movement and blood flow velocity within blood vessels.

Doppler imaging is based on detecting frequency shifts caused to the reflected ultrasound signal by the movement of blood cells or other moving tissue. The Doppler power spectrum, representing the distribution of the scatterers velocities, is conventionally estimated by calculating peridodograms from the received ultrasound signal using the Welch method. A Doppler spectrogram is displayed to visualize the evolution of blood velocity over time.

For acceptable spatial resolution, B-mode imaging typically requires transmitting short-time pulses that are wideband with a high center frequency. Spectral Doppler imaging, on the other hand, requires transmitting narrowband pulses at a low center frequency for achieving high frequency-resolution and sufficient penetration depth. In practice, B-mode advantageously operates in parallel to the Doppler mode, e.g., for orientation.

Methods for combining between the B-mode and the Doppler mode are known in the art. In one approach, pulses used for the two modalities are transmitted interleaved, e.g., by repeatedly transmitting pairs of a B-mode pulse followed by a Doppler-mode pulse. Such a transmission pattern, however, reduces the effective Pulse Repetition Rate (PRF) by half, which according to the Nyquist sampling theorem also reduces the maximal measurable velocity by half. Interleaving between the Doppler-mode and B-mode pulses can also be done in a non-uniform manner, as described, for example, by Jensen in "Spectral velocity estimation in ultrasound using sparse data sets," Journal of the Acoustical Society of America (ASA), volume 120, number 1, 2006, pages 211-220.

Recovering the Doppler power spectrum typically requires taking the non-uniform transmission pattern into consideration. For example, iterative methods for estimating Doppler spectrogram using arbitrary transmission patterns are described by Gudmundson et al. in "Blood velocity estimation using ultrasound and spectral iterative adaptive approaches," Signal Processing, volume 91, number 5, 2011, pages 1275-1283 and by Gran at al. in "Adaptive spectral Doppler estimation," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, volume 56, number 4, 2009.

In "Blood velocity estimation using compressive sensing," IEEE Transactions on Medical Imaging, volume 32, number 11, 2013, pages 1979-1988, Richy et al. assume that the Doppler signal is sparse under the Fourier transform or in the wave atom domain, and propose to decompose the Doppler signal into several equal segments that each is recovered using Compressed Sensing (CS) recovery techniques.

Known methods in which the Doppler power spectrum is estimated from a nonuniform sequence of transmitted pulses, such as those given by example above, are typically computationally complex, recover the power spectrum with artifacts, or both. SUMMARY

An embodiment that is described herein provides an apparatus for measuring target velocity including a transducer and a processor. The transducer is configured to transmit into a target a set of pulses, which is sparse but has a recoverable power spectrum, and to receive a signal reflected from the target in response to the set of pulses. The processor is configured to estimate a velocity of a selected region in the target by analyzing the reflected signal.

In some embodiments, the set of pulses is based on a sequence containing, within an observation window, a smallest number of pulses for which the power spectrum is recoverable. In other embodiments, (i) respective start times of the pulses are selected multiples of a time unit, and are spaced non-uniformly, and (ii) for every possible multiple of the time unit up to a predefined maximal multiple, the set of pulses includes at least one pair of pulses having start times that differ by that multiple. In yet other embodiments, the set of pulses occupies a series of time intervals, such that in each time interval the set of pulses includes two or more pulse sequences. In an embodiment, the two or more pulse sequences have respective different Pulse- Repetition Frequencies (PRFs). In another embodiment, the two or more pulse sequences include a first pulse sequence having a first PRF followed by a second pulse sequence having a second PRF. In yet another embodiment, the two or more pulse sequences include at least first and second pulse sequences that at least partially overlap one another in time.

In some embodiments, the set of pulses includes at least one pulse that is common to at least two of the pulse sequences. In other embodiments, the set of pulses and the reflected signal include ultrasound signals.

In an embodiment, the processor is configured to estimate the power spectrum of the reflected signal, and to derive the velocity from the estimated power spectrum. In another embodiment, the processor is configured to estimate the power spectrum by calculating a covariance of the reflected signal, to produce a measurement vector by averaging elements of the covariance corresponding to a common inter-pulse interval in the set of pulses, and to estimate the power spectrum from the measurement vector. In yet another embodiment, the processor is configured to estimate the power spectrum over predefined Doppler frequencies by applying to the measurement vector a Fast Fourier Transform (FFT).

In some embodiments, the processor is configured to estimate the power spectrum by reordering the measurement vector as a Toeplitz matrix, and applying to the Toeplitz matrix eigenvalue decomposition. In other embodiments, the processor is configured to estimate an actual number of Doppler frequencies of the power spectrum based on the eigenvalue decomposition. In yet other embodiments, the processor is configured to estimate the power spectrum by applying a regularization functional to the measurement vector. In yet further other embodiments, the processor is configured to apply to the reflected signal a Finite Impulse Response (FIR) filter in a correlation domain.

There is additionally provided, in accordance with an embodiment that is described herein, a method for measuring target velocity including transmitting into a target, by a transducer, a set of pulses, which is sparse but has a recoverable power spectrum. A signal reflected from the target in response to the set of pulses is received by the transducer. A velocity of a selected region in the target is estimated by analyzing the reflected signal.

These and other embodiments will be more fully understood from the following detailed description of the embodiments thereof, taken together with the drawings in which: BRIEF DESCRIPTION OF THE DRAWINGS

Fig. 1 is a block diagram that schematically illustrates an imaging system that transmits a sparse emission pattern in Doppler mode, in accordance with an embodiment that is described herein;

Fig. 2 is a diagram that schematically illustrates a uniform transmission pattern and two nested array transmission patterns, in accordance with an embodiment that is described herein;

Fig. 3 is a flow chart that schematically illustrates a method for estimating the power spectrum of an ultrasound signal in which the Doppler frequencies are confined to predefined discrete values, in accordance with an embodiment that is described herein; and

Fig. 4 is a flow chart that schematically illustrates a method for estimating power spectrum of ultrasound signal so that the Doppler frequencies are estimated over a continuous frequency range, in accordance with an embodiment that is described herein.

DETAILED DESCRIPTION OF EMBODIMENTS OVERVIEW

Doppler ultrasound is commonly used in various medical applications. In the description that follows we focus mainly on spectral Doppler techniques for estimating blood flow velocity. The disclosed techniques are similarly applicable, however, to imaging other types of moving tissue, as well as to other fields involving moving targets such as in radar and communications applications.

Embodiments that are described herein provide methods and systems for spectral Doppler imaging in which the transmitted sequence of ultrasound pulses is sparse, but allows recovery or estimating a good approximation of the Doppler power spectrum. Consider a uniform sequence of P transmissions, i.e., transmitting a stream of P pulses with uniform spacing. A sparse sequence with respect to the P-pulse uniform sequence contains a partial subset of N pulses among the P pulses, wherein N is strictly smaller than P.

A measurement cycle typically starts with the ultrasound probe transmitting toward the target a sequence of ultrasound pulses modulated using some center frequency. The measurement cycle is also referred to as an "observation window."

A naive approach would be to transmit the pulses at a constant rate referred to as a

Pulse Repetition Frequency (PRF). A scatterer in the target such as a red blood cell in a vessel would reflect such transmitted pulses at a frequency proportional to its axial velocity and to the center frequency. In practice, the ultrasound pulses are reflected from multiple blood cells whose velocities are distributed in accordance with a time-varying distribution that can be estimated by analyzing pulses reflected from a desired depth of interest. A signal containing samples of the received ultrasound signal sampled at multiples of the PRF interval is referred to as a "slow-time signal," and serves for evaluating the Power Spectral Distribution (PSD). The PSD is displayed as a spectrogram that visualizes the evolution of the velocity distribution over time on a display device.

In spectral Doppler imaging, the frequency resolution is inversely related to the number of pulses transmitted in an estimation cycle, which is in turn limited by (i) the measurement duration, which should be shorter than the Coherence Processing Interval (CPI) during which the velocities of the scatterers are assumed to be approximately constant, and (ii) the minimal PRF required for imaging a desired depth.

In the disclosed embodiments, an ultrasound apparatus comprises a transducer and a processor. The transducer transmits into a target tissue a set of pulses, which is sparse but has a recoverable power spectrum, and receives a signal reflected from the target tissue in response to the set of pulses. The processor derives the set of pulses from a sequence containing a small number of pulses (e.g., the smallest number of pulses) in an observation window for which the power spectrum is recoverable, and estimates the velocities of scatterers in a selected region in the target tissue by analyzing the reflected signal.

In some embodiments, the processor derives the set of pulses so that (i) respective start times of the pulses are selected multiples of a time unit (e.g., the PRF interval) and are spaced non-uniformly, and (ii) for every possible multiple of the time unit up to a predefined maximal multiple, the set of pulses comprises at least one pair of pulses having start times that differ by that multiple. The second property ensures that a covariance function of the reflected signal from which the power spectrum is derived, is fully recoverable for all required time lags. Alternatively, an approximated covariance function can be recovered at a high approximation quality.

The processor may define the set of pulses in various ways. In one embodiment, the set of pulses occupies a series of time intervals, such that in each time interval the set of pulses comprises two or more pulse sequences that may have respective different Pulse-Repetition Frequencies (PRFs). For example, the two or more pulse sequences comprise a first pulse sequence having a first PRF followed by a second pulse sequence having a second PRF. As another example, the two or more pulse sequences comprise at least first and second pulse sequences that at least partially overlap one another in time. In some embodiments, the set of pulses comprises at least one pulse that is common to at least two of the pulse sequences.

In some embodiments, the processor is configured to estimate the power spectrum of the reflected signal, and to derive the velocities of the scatterers in the target from the estimated power spectrum. The processor may estimate the power spectrum using any suitable method. In some embodiments, the processor estimates the power spectrum by calculating a covariance matrix of the reflected signal, and producing a measurement vector by averaging elements of the covariance matrix corresponding to a common inter-pulse interval in the set of pulses.

In one embodiment, the processor estimates the power spectrum over predefined

Doppler frequencies by applying to the measurement vector a Fast Fourier Transform (FFT). In another embodiment, the processor estimates an actual number of Doppler frequencies of the power spectrum by performing an eigenvalue decomposition to a Toeplitz matrix derived from the measurement vector. In yet other embodiments, the processor is configured to estimate a smoothed power spectrum by imposing a quadratic or any other suitable regularization functional to the measurement vector.

In some applications, the reflected signal contains strong clutter reflected from non- moving tissue. In some embodiments, the processor applies to the reflected signal a Finite Impulse Response (FIR) filter in the correlation domain, for removing the clutter from the reflected signal.

In the disclosed techniques, an ultrasound system uses an irregular sparse emission pattern. The sparse emission partem may be selected with a minimal number of pulses that allows enhanced-resolution recovery of the Doppler power spectrum in the correlation domain. Using the disclosed techniques, the ultrasound system may interrogate several blood regions simultaneously, while still updating B-mode imaging at a high frame rate. An efficient method for recovering the power spectrum, which is tailored to the sparse emission pattern is also presented.

Example experimental and simulated results are provided in U.S. Provisional Patent Applications 62/474,164 and 62/485,955, cited above. The results demonstrate accurate estimation of the blood velocity spectrum using only 12% of the available transmission slots.

SYSTEM DESCRIPTION

Fig. 1 is a block diagram that schematically illustrates an imaging system 20 that transmits a sparse emission pattern in Doppler mode, in accordance with an embodiment that is described herein. Imaging system 20 is typically used for imaging various types of moving tissues, e.g., blood flow velocity within a selected blood vessel 22 in an organ such as a kidney, brain or the heart. Spectral Doppler techniques may also be applied for estimating displacement of solid tissue such as cardiac-chamber wall, and for estimating the velocity of contrast agents injected into the circulatory system. The disclosed techniques are also applicable to other Doppler modes such as, for example, color Doppler and tissue Doppler.

Imaging system 20 comprises an ultrasound probe 30, which comprises a transducer array 32 of transducer elements 34. Transducer element 34 converts between electrical signals and ultrasonic signals in transmit and receive directions, respectively. During imaging, the physician typically couples transducer array 32 to the patient body. Transducer array 32 transmits an ultrasonic signal in the form of a sequence of pulses toward the target tissue, and receives reflected signals or "echoes" of the transmitted pulses from the tissue. The reflected signals received are processed so as to estimate the distribution of velocities of the red blood cells at the target, as will be described blow.

In the embodiment of Fig. 1, imaging system 20 comprises a processor 40, which is coupled to ultrasound probe 30 via an interface 44. The interface is coupled to the ultrasound probe using a suitable link 46, which may comprise any suitable cable, typically connected at both ends, electrically and mechanically, using suitable connectors. Processor 40 exchanges TX signals and RX signals with ultrasound probe 30 via interface 44.

In the transmit path, a TX beamforming module 48 generates TX signals for transmitting ultrasound waves via transducer elements 34 of the transducers array. In some embodiments, TX beamforming module 48 adjusts the amplitudes and phases of the TX signals so that transducer elements 34 together emit an ultrasound plane wave toward the target. In some embodiments, interface 44 comprises a Digital to Analog Converter (DAC) 52 for converting digital signals produced by the TX beamforming module into analog signals for controlling pulse generator 54, which generates a sequence of pulses toward ultrasound probe 30.

Pulse generator 54 modulates the transmitted pulses with a sinusoidal signal having a predefined center or carrier frequency denoted†Q . The carrier frequency is selected depending on the transducer elements used and on the desired imaging depth. In an example embodiment, the carrier frequency is configured to 3.5MHz, but other suitable carrier frequencies can also be used. In the receive path, ultrasound probe 30 receives ultrasound signals containing echoes of the ultrasound pulses reflected by red blood cells and the surrounding tissue in the target. Interface 44 processes the signals received from the ultrasound probe using an analog front- end module, which typically comprises elements such as a low-noise amplifier, a low pass filter and a sampler (not shown). An Analog to Digital Converter (ADC) 58 converts the signals output by the analog front-end into digital RX signals to be processed by processor 40.

Interface 44 provides the RX signals to a RX beamforming module 62, which focuses the RX signals at selected direction and depth using dynamic focusing methods as known in the art. In an embodiment, RX beamforming module 62 delays the RX signals of respective transducer elements of the transducer array and sums the delayed signals using suitable sum- weights. In an embodiment, RX beamforming module 62 applies to the summed signal a bandpass filter (not shown) that contains the carrier frequency, for removing noise outside the passband supported by transducer elements 34. In some embodiments, processor 40 demodulates the RX signals, prior to beamforming, using the carrier frequency of the transmitted pulses to produce In-phase and Quadrature (IQ) signals (not shown).

A Doppler analyzer 70 processes the beamformed Rx signals to estimate the velocity distribution of the imaged red blood cells. In the disclosed techniques, the Doppler analyzer evaluates the velocities by estimating the spectral Doppler content of the RX signals caused by the movement of the blood cells. Doppler processing techniques that are based on estimating an autocorrelation function of the Rx signals resulting from transmitting a sparse sequence of pulses will be described in detail below. Alternatively, any other suitable Doppler technique for estimating the blood velocity can also be used.

An image reconstruction module 74 receives from Doppler analyzer 70 multiple power spectra, estimated over multiple respective estimation cycles, and reconstructs from these spectra a Doppler spectrogram to be displayed on a display 78.

In some embodiments, imaging system 20 is a duplex imaging system configured to display both B-mode images and spectral Doppler data. In such embodiments, the imaging system is configured to transmit pulses suitable for producing B-mode images, in addition to the pulses suitable for generating the Doppler data. For example, the physician first images the target using B-mode, and in response to identifying a region of interest that contains a blood vessel, the physician switches to the Doppler mode for imaging the blood flow in that vessel. Alternatively, B-mode pulses are transmitted in available pulse locations of the sparse sequence defined for the Doppler-mode pulses, allowing simultaneous display of B-mode images and spectral Doppler data. The system configuration of Fig. 1 is an example configuration, which is chosen purely for the sake of conceptual clarity. In alternative embodiments, any other suitable system configuration can be used. The elements of imaging system 20 may be implemented using hardware. Digital elements can be implemented, for example, in one or more off-the-shelf devices, Application-Specific Integrated Circuits (ASICs) or FPGAs. Analog elements can be implemented, for example, using discrete components and/or one or more analog ICs. Some system elements may be implemented, additionally or alternatively, using software running on a suitable processor, e.g., a Digital Signal Processor (DSP). Some system elements may be implemented using a combination of hardware and software elements.

In some embodiments, some or all of the functions of imaging system 20 may be implemented using a general-purpose processor (e.g., processor 40,) which is programmed in software to carry out the functions described herein. The software may be downloaded to the processor in electronic form, over a network, for example, or it may, alternatively or additionally, be provided and/or stored on non-transitory tangible media, such as magnetic, optical, or electronic memory.

System elements that are not mandatory for understanding of the disclosed techniques, such as circuitry elements relating to interface 44, are omitted for the sake of clarity or only briefly discussed.

SIGNAL MODEL FOR SPARSE TRANSMISSION PATTERNS

Consider a transmission pattern in which ultrasound pulses are transmitted at time instances that are multiples of a time unit T defined as the reciprocal of the Pulse Repetition Frequency (PRF), i.e.,

Within an observation window of time-duration PT, ultrasound probe 30 transmits a sequence of ultrasound pulses comprising

pulses given by:

Equation 1:

wherein p_n gets values in the range 1 ... P. The pulse function is defined over

the interval with a center frequency given by:

Equation 2:

wherein g(t) is an envelope function (such as a Gaussian window or some other suitable envelope). In the description that follows we assume that

purely for the sake of clarity.

The signal received at the transducer corresponding to the

pulse reflected by a single scatterer in the target is given by:

Equation 3:

In Equation 3, d₀ denotes the initial depth of the scatterer, which is assumed to be moving at a velocity V along the direction of the ultrasound beam. In addition, C is the propagation speed of sound, and the amplitude a depends on the reflectivity of the scutterer in question.

In practice, each of the N transmitted pulses is reflected by multiple scutterers that each moves at one of M axial velocities corresponding to respective M unknown

Doppler frequencies given as

Equation 4 :

The signal acquired over the observation window, after beamformed using RX beamforming module 62 and demodulated by is organized in a two-dimensional array

given by:

Equation 5:

wherein denotes a complex-valued amplitude

associated with the scatterers moving at a respective velocity . In Equation 5, the index k

corresponds to samples of the received signal within a time interval T. The index k is associated with depth within the target. In some embodiments, The Doppler frequencies are assumed to belong to an unambiguous frequency range

For a given pulse index the respective column elements of y

form a "fast-time signal" of the received signal. For a given index k, the respective row elements of y[k, 7l], i.e., one sample per transmitted pulse, form a "slow-time signal" of the received signal corresponding to the depth associated with k. Equation 5 can be written in a matrix form as:

Equation 6:

wherein the slow-time signal contains N samples from N respective emissions,

a is a column vector of length M comprising the amplitudes ^m& the model matrix

A_N is given by:

Equation 7:

In the context of the present disclosure, the Doppler power spectrum is defined by a set of M Doppler frequencies f . M , and a set of M variances that are indicative of

the respective velocities of the scatterers. Recovering the power spectrum is thus equivalent to estimating the M Doppler frequencies and respective Af variances. In terms of the model of Equation 6, the goal of Doppler analyzer 70 is to recover from the measurements y the

set of frequencies f_m defining the model matrix A_N, and the variances °f the

respective elements of the vector a.

Note that in conventional Doppler ultrasound, the ultrasound probe typically transmits a uniform sequence of P pulses over the PT observation window. Assuming that the Doppler frequencies lie on the Nyquist grid, i.e., are multiples of P

, the matrix A_N for

wherein F is the

matrix. In this case, the power spectrum is conventionally estimated by averaging multiple periodograms, each calculated from a different measurement vector

resulting in a spectral resolution

Since the minimal is determined by the maximal imaging depth required, there is a

tradeoff between temporal resolution, which for a given PRF is proportional to P and spectral resolution, which is inversely proportional to P.

In the disclosed embodiments, the Doppler power spectrum is estimated based on an autocorrelation function of the slow-time signal. The term "autocorrelation function" is also referred to herein as an "autocorrelation matrix," "covariance function" or "covariance matrix." The autocorrelation matrices of the slow-time signal and of a are defined respectively as resulting in a correlation-based

model:

Equation 8 :

Assuming that the amplitudes of

are uncorrected,

is an M diagonal

matrix, whose diagonal elements comprise the vector of power spectrum. Equation 8 can thus be rewritten as:

Equation 9:

wherein denotes a measurement vector produced by stacking the columns of

and denotes the power spectrum vector to be estimated.

matrix given by wherein the symbol O denotes the Khatri-Rao product, defined as a

column-wise Kronecker product between its matrix arguments.

The elements of the

column of the matrix

have the form e

wherein are pulse locations in the transmission pattern.

It can be shown that for certain sparse transmission patterns for which the

matrix in the model of Equation 9 has a full column rank, in which case the power spectrum P is recoverable from the measurements

In some embodiments, processor 40 generates a sparse transmission pattern that contains a minimal number of pulses N < P that allows high-precision recovery of the power spectrum by Doppler analyzer 70. Moreover, the resulting frequency resolution improves by a factor of about two relative to the frequency resolution achievable in priodogram-averaging power spectrum estimation, at a comparable complexity. In some of the disclosed embodiments, Doppler analyzer 70 estimates the Doppler power spectrum based on the autocorrelation of the slow-time signal, as given in Equation 9, without explicitly recovering the slow-signal itself.

SPARSE TRANSMISSION PATTERNS BASED ON NESTED ARRAYS

In some disclosed embodiments, processor 40 generates a sparse pattern of pulses for transmission in an observation window based on nested arrays. Nested arrays are non-uniform sensor arrays used, for example in Multiple-In Multiple-Our (MIMO) radar, and in estimating Direction of Arrival (DOA). Nested arrays are produced by nesting (concatenating) two Uniform Linear Arrays (ULAs). Nested arrays are described, for example, by Pal and Vaidyanathan in "Nested arrays: A novel approach to array processing with enhanced degrees of freedom," IEEE Transactions on Signal Processing, volume 58, number 8, 2010, pages 4167-4181.

In the context of the present disclosure, the term "non-uniform" means that the pulses are not transmitted at a fixed PRF, i.e., not spaced at a fixed time delay from one another.

In adopting the concept of nested array to transmitting sparse pulse-sequences in Doppler ultrasound applications, consider a sequence S_N of N < P pulses defined as a union of two subsequences denoted and , given by:

Equation 10:

The subsequences and S contain N1 and N2 pulses, respectively, wherein the

integers Nl and N2 satisfy:

Equation 11:

In accordance with Equation 10, S_N comprises a number Nl of pulses uniformly inter-spaced at a unity interval, followed by N2 pulses uniformly inter-spaced at an interval of Nl + 1 units. Note that for a given N, different combinations ofNl and N2 correspond to different respective sequences S_N. Based on the nested-array sequence S_N, the sparse transmission pattern is created by transmitting iV pulses, wherein the transmission time of n pulse isp

, andp_n gets values in the range 1 ... P.

An important property of the nested array sequence of Equation 10 is that the set of differences between neighboring pulses contains the entire range of integers between

This property allows estimating the autocorrelation matrix R_yN of Equation 8 for all time-lags between - (P— 1) and P— 1, even though the transmitted pattern contains less than P pulses.

Fig. 2 is a diagram that schematically illustrates a uniform transmission pattern and two nested array transmission patterns, in accordance with an embodiment that is described herein. The observation window in Fig. 2 comprises up to P = 12 pulses at a minimal interspace interval of

time units.

In Fig. 2, the upper sequence, serving as a reference, is a uniform sequence containing all P = 12 pulses. In this case N

and

The middle sequence is a nested array sequence in which

and the bottom sequence is another nested array sequence in which

and N

. Note that in this example, in both the nested array sequences

, , only half of the possible number of pulses are transmitted for spectral Doppler estimation.

In some embodiments, for a given P, the sparse transmission pattern is configured so that N = Nl + N2 is minimal among all the possible nested array combinations satisfying Equation 11.

Table 1, depicts several combinations for constructing a nested array sequence, given an observation window of sizeP = 128.

Table 1 : various nested array patterns for P = 128

As seen in the table, the minimal solution is N = 23, which can be achieved by selecting

More generally, for a given P, define two sets:

and

], wherein the expression d\P means that d is a divisor of P. It can be shown that the minimal N is achievable by each of the two combinations:

Equation 12a:

Equation 12b:

In the example of P=128, Dl and D2 in equations 12a and 12b are given by D1={1, 2, 4, 8} and D2={16, 32, 64, 128}. Note that although both Equations 12a and 12b result in the same minimal N, the respective transmission patterns are different. The nested array pattern of Equation 10 contains N2— 1 gaps of size Nl + 1, which can be used, e.g., for B-mode imaging. For example, in coherent plane-wave compounding, the size of the gap determines the number of inclination angles, or equivalently the B-mode image quality, whereas the number of gaps relates to the imaging frame rate.

When the observation window is selected such that

has an integer value, the minimal number of pulses in the nested array transmission pattern is given by

1 wherein

In some embodiments, processor 40 generates a sparse transmission pattern based on Co-Prime arrays. In such embodiments, Nl and N2 are co-prime integers selected so that

Nl < N2. The co-prime sparse transmission sequence is defined as:

Equation 13:

In Equation 13, Nl and N2 should be selected so that

is satisfied.

Similarly to the nested array case, co-prime arrays also enable estimating the autocorrelation matrix of Equation 8 for all time-labs. In using co-prime arrays, S_N1 and S_N2 may share a common pulse. Alternatively or additionally, S_N1 and S_N2 may overlap in time, at least partially. Transmission patterns based ono-prime can be extended similarly to more than two subsequences of pulses.

The main drawback of co-prime arrays is that they require sending Doppler pulses at time instances beyond the observation window. Therefore, the reflected slow-time signal may not preserve its stationarity property, which is a key assumption in Doppler processing. In one embodiment, to overcome this, pulse transmissions beyond the observation window are omitted. In this embodiment, however, some of the time-lags for estimating the autocorrelation matrix may be missing, which may reduce the number of Doppler frequencies (or scatterer velocities) that can be recovered.

In some alternative embodiments, processor 40 constructs a sparse transmission pattern using a K-level nested array, which is an extension of the nested array describe above. A K- level nested array is defined by a number K of subsequences having respective lengths Ni so that the inter-spacing at the level is given by . It can be shown that if P can be factorized into multiple powers of distinct prime factors, a K-level nested sequence can be found such that the sum of Ni is minimal. When the size of the observation window P is a power of two, i.e.,

for some integer 71, the optimal K-level nested array contains

1 pulses with exponential inter-spacing pattern, as given by:

Equation 14:

Although the K-level nested array may result in a number of emissions N that is significantly smaller than in the nested array of Equation 10, higher-order statistics are involved, which requires larger amounts of measurement snapshots for estimating the autocorrelation matrix.

Another alternative to the nested array of Equation 10 is referred to as a "super-nested array." Such nested sequences are defined for Nl > 4 and N2 > 3 by concatenating six suitable ULAs. Super-nested arrays are described, for example, by Liu and Vaidyanathan in "Super nested arrays: Sparse arrays with less mutual coupling than nested arrays," the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2016, pages 2976-2980.

The super-nested arrays share the same properties as nested arrays in terms of the number of Doppler transmissions and their difference sets. In addition, super-nested arrays are advantageous over nested arrays in terms of a reduced undesirable effect of previous transmissions on a received signal corresponding to the current emission. This property may allow increasing the maximal depth examined. The super-nested arrays exhibit, however, complex geometry compared to the nested arrays. In particular, the inter-spacing gaps created may have different sizes, which is less suitable for B-mode imaging.

Additional sparse transmission patterns that processor 40 can use are based on the "minimum redundancy arrays" (MRAs) and the "minimum hole arrays" (MHAs). MRAs are described, for example, in "Minimum-Redundancy Linear Array," IEEE Transactions on antennas and propagation, volume AP-16, number 2, March 1968. MHAs are described, for example, in "Rulers, Part L" University of Southern California, Los Angeles, CA, USA, Technical Report 85-05-01, 1985.

In the embodiments described above, the sparse transmission patterns are designed so that the set of differences among the pulses in the observation window contain all the integers in a range that allows estimating the autocorrelation matrix over all time-lags, which allows recovery of the Doppler power spectrum. In alternative embodiments, other suitable sparse transmission patterns can also be used.

EFFICIENT POWER SPECTRUM RECOVERY

The model given in Equation 6 for solving the power spectrum is typically an ideal model. In practice, the received signal is noisy, which can be modeled as:

Equation 15:

wherein is the received signal containing a number Q of snapshots,

wherein Q is proportional

For example, may be configured as

Alternatively, set wherein or using any other suitable multiple of

^{ls a} vector of zero-mean white complex-valued Gaussian noise with an unknown covariance matrix The noise samples are assumed to be uncorrelated with the amplitudes

a.

The covariance matrix R_yN can be approximated based on a finite set of measurements as given by:

Equation 16:

Incorporating these factors to the model of Equation 8 gives:

Equation 17:

and in a vector form similarly to Equation 9:

Equation 18:

Let D denote the set of differences between the pulses in the nested array transmission pattern S_N of Equation 10. The set D is given by D ] AS noted

above, the elements of A_N have the form wherein is the element of

D. Since the set D has duplicate elements, the system of equations in Equation 18 is redundant, i.e., some of the rows in are identical.

To reduce the equation system we define a set D_u of the unique elements of D, and for each element d G D_u the set

contains the indices in which the element d appears in D. Next we average elements of the measurement vector sharing the

same d to produce a compact measurement vector Z of size I

Equation 19:

Wherein denotes the index of d in D and | is the number of times d occurs in D_u. Using Equation 19, the model of Equation 19 becomes:

Equation 20:

In Equation 20,

is a matrix of

rows and M columns, and the elements of

A are given by:

Equation 21:

wherein being the element of D_u. The vector e has a size

and all its elements are zero expect the element in position P, which equals

In some embodiments, Doppler analyzer 70 solves Equation 21 under the constraint according to which the Doppler frequencies confined to the Nyquist grid, i.e,

for every 1 < m≤ M, and i_m is an integer in the range 0

. Note that the spectral resolution of this grid is compared to in conventional

priodogram-based spectral Doppler estimation, i.e., an improvement by a factor of almost two.

When the Doppler frequencies are confined to the Nyquist grid,

wherein F is the

matrix. In this case the power spectrum P is related to the Fourier transform F

Equation 22:

wherein 1 is a ones-vector of size P. Equation 22 implies that the power spectrum approximately equals z.

Fig. 3 is a flow chart that schematically illustrates a method for estimating the power spectrum of an ultrasound signal in which the Doppler frequencies are confined to predefined discrete values, in accordance with an embodiment that is described herein. The method will be described as being executed by various elements of processor 40 and interface 44 of Fig. 1.

The method begins with processor 40 transmitting a sparse pulse sequence S_N of N pulses over an observation window P > N using Tx beamforming module 48, DAC 52 and pulse generator 54, at a transmission step 100. In the present example we assume that the sequence of pulses is transmitted in accordance with a nested array pattern such as S_N defined in Equation 10 above. The number of pulses N may be selected to be the minimal number of pulses that allows power spectrum reconstruction, as given, for example, in Equations 12a and 12b, above. Alternatively, any other number of pulses smaller than the observation window can also be used.

At a reception step 104, the processor receives via analog front end 56, ADC 58 and Rx beamforming module 62 a signal comprising reflections of the transmitted pulses from scatterers in the target. The slow-time signal of the received signal (after beamformed using

RX beamforming module 62) is { , which is assumed to satisfy the model of

Equation 15 above.

At a pre-processing step 108, Doppler analyzer 70 generates a measurement vector z by applying steps 112, 116 and 120 as described herein. At a covariance matrix calculation step 112, the Doppler analyzer calculates the covariance matrix of the received signal as

given by Equation 16. The Doppler analyzer stacks the columns of to produce the vector

of Equation 18, at a vectorization step 116. At a measurement vector generation step 120, the Doppler analyzer generates the vector z by averaging relevant elements of r_N as given in Equation 19. At a transformation step 124, the Doppler analyzer transforms Z into a transformed vector Z by applying a Fourier transform as given in Equation 22, wherein F is a

P-by-P Fourier matrix, and wherein

At a power spectrum estimation step 128, the Doppler analyzer calculates the power spectrum P of the received signal by applying to z a soft thresholding operator . The soft thresholding operator

reduces the noise variance and the effect of spurious frequencies resulting by estimating the covariance matrix from a finite set of snapshots.

Following step 128 the method loops back to step 100 for transmitting a sequence of N pulses during a subsequent observation window. In an embodiment, the processor transmits the same sequence over successive observation windows. In alternative embodiments, the processor may transmit different pulse sequences over different observation windows.

The complexity of the method of Fig. 3 depends on the size P of the observation window, the number of pulses N transmitted during the window and the number of snapshots Q acquired during the observation window. For given P, N and Q,the complexity is

. For a pulse sequence having a minimal N value, N² is proportional to the observation window size, and the complexity is which is suitable for

real-time spectral Doppler applications.

When assuming that the power spectrum to be estimated is smooth, instead of smoothing over multiple estimated power spectra as a post processing phase, the problem to solving a smooth power spectrum can be formalized, for example, as a quadratic regularization functional that imposes smoothness on the power spectrum, as given by:

Equation 23:

Wherein μ is a predefined smoothing parameter, and D_c is a circulant matrix in which the first row is defined as (1, -1, 0,...,0), and each of the other rows is a cyclically shifted to right version of its preceding row. Alternatively, any other regularization functional can also be used. A closed form solution to the problem of Equation 23 is given by

Equation 24:

Wherein G is a circulant matrix whose elements are given by:

Equation 25:

In some embodiments, the processor executes a variant of the method of Fig. 3, by

Doppler analyzer 70 calculating Z of step 124 using Equation 23. Note that setting μ = 0 in Equation 23 results in the same solution as in Equation 22. The complexity of this variant method is

In the method of Fig. 3, the power spectrum is estimated under the constraint that the Doppler frequencies are confined to predefined discrete values, i.e., on the Nyquist grid. In practice, however, the actual Doppler frequencies typically get values on a continuous frequency range, and are not limited to predefined discrete values.

Next we describe a method for estimating the power spectrum based on a model given by:

Equation 26:

wherein R is a Toeplitz matrix given by:

Equation 27:

and A is a whose entries are given by A

It can be shown that at noise-free conditions R and A share the same range space, which can be exploited to recover the Doppler frequencies using suitable subspace methods such as those described by Eldar in "Sampling Theory: Beyond Bandlimited Systems," Cambridge University Press, 2015. In Fig. 4 below we describe a method for spectral Doppler estimation based on a subspace method referred to as the "Esprit" method, which is described, for example, by Roy and Kailath in "Esprit-estimation of signal parameters via rotational invariance techniques," IEEE Transactions on Acoustics, Speech, and Signal Processing, volume 37, number 7, 1989, pages 984-995.

Fig. 4 is a flow chart that schematically illustrates a method for estimating power spectrum of ultrasound signal so that the Doppler frequencies are estimated over a continuous frequency range, in accordance with an embodiment that is described herein. The method will be described as being executed by various elements of processor 40 and interface 44 of Fig. 1.

The method begins with processor 40 transmitting a sequence of ultrasound pulses at a transmission step 200, receiving an ultrasound signal at a reception step 204, and producing, using Doppler analyzer 70, a measurement vector Z at pre-processing step 208. Steps 200, 204 and 208 are essentially the same as respective steps 100, 104 of Fig. 3, including sub-steps 112, 116 and 120 of step 108.

At a Toeplitz matrix generating step 212, the Doppler analyzer generates the matrix R of Equation 27, and decomposes this matrix using eigenvalues decomposition as given by:

Equation 28:

wherein J? is a matrix whose columns comprise the eigenvectors and b is a vector comprising the respective eigenvalues, in a non-increasing order.

The Esprit method typically assumes that the number M of Doppler frequencies is known. In Doppler ultrasound, however, the number of underlying Doppler frequencies is usually unknown, and may change over multiple observation windows. In some embodiments, the processor estimates M using any suitable method such as using the Minimum Description

Length (MDL) method. In the present example, the processor estimates Λί, at a number of frequencies estimation step 216, using the following expression:

Equation 29:

wherein is a thresholding operator whose threshold parameter λ is configured

empirically, or using any suitable other method. The f | - 11 o norm is the IQ semi -norm that counts the number of nonzero elements in its argument vector.

At a partial eigenvector matrix calculation step 220, the Doppler analyzer constructs a matrix E_M comprising the M eigenvectors in E corresponding to the M largest eigenvalues in b, wherein E and b were derived at step 212 above. Further at step 220, the Doppler analyzer constructs a matrix denoted E₁ comprising the first P— 1 rows of E_M , and a matrix E₂ comprising the last P— 1 rows of E_M. It can be shown that if the condition M < P— 1 is satisfied, the following expression holds:

Equation 30:

Wherein Λ is an

diagonal matrix whose diagonal elements are given by

and W is a suitable invertible matrix. Therefore, the Doppler frequencies can be estimated from the eigenvalues of as described herein. The

superscript operator applied to means taking the pseudo-inverse matrix of E<±.

At a Doppler frequencies estimation step 224, the processor calculates the eigenvalues and estimates the Doppler frequencies by calculating:

Equation 31:

At a model matrix construction step 228, the Doppler analyzer constructs the model matrix A whose elements are as given in Equation 21 above.

At a power spectrum estimation step 232, the Doppler analyzer calculates the power spectrum at the estimated Doppler frequencies as:

Equation 32:

Following step 232 the method loops back to step 200 for transmitting another sparse sequence of pulses in a subsequent observation window.

The method of Fig. 4 identifies the Doppler frequencies at a resolution much higher than the method of Fig. 3, especially at high Signal to Noise Ratio (SNR). The complexity of the method of Fig. 4 is dominated by the eigenvalue decomposition operation applied to the P- by-P Hermitian matrix E at step 224, which complexity is In alternative

embodiments, other suitable methods whose complexity is lower than the complexity of the Esprit method can also be used, such as methods that are described, for example, by Xu and Kailath in "Fast subspace decomposition," IEEE Transactions on Signal Processing, volume 42, number 3, 1994, pages 539-551.

CLUTTER FILTERING AND WINDOWING IN THE CORRELATION DOMAIN

In some embodiments, the received ultrasound signal is degraded by a high level of clutter reflected by nonmoving tissue or by tissue that moves very slowly relative to the moving scatterers. The clutter level may be on the order of 40dB-60dB higher than the desired signal and should be eliminated or reduced significantly by using a suitable clutter filter.

In some embodiments, Doppler analyzer 70 may apply a clutter filter (not shown) to the received signal for separating between echoes reflected by blood cells flowing in the blood stream, and the static tissue - generally referred to as clutter. In some embodiments, the clutter filter is designed with a low cutoff frequency (e.g., 0.03-PRF) in order to include slow flowing red blood cells while removing clutter artifacts. Efficient clutter filtering in the correlation domain will be described in detail below.

Note that a Finite Impulse Response (FIR) clutter filter is typically designed for a uniform sampling rate, and is not suitable for filtering the slow-time signal resulting from a non-uniform sequence of pulses transmitted during the observation window. In some embodiments, since the non-uniform pattern of the transmitted pulses is designed so that the autocorrelation function is recoverable over all lags, the clutter filtering operation is advantageously applied in the correlation domain rather than in the time domain, as explained herein.

Consider a signal that is filtered by a FIR filter having an impulse response

function to produce an output signal X Let and denote the

autocorrelation functions of the input and output signals, respectively. Using these notations, the following expression holds: Equation 33:

wherein the operator '*' denotes a convolution operation. Since the measurement signal z of Equation 19 is derived from an autocorrelation matrix of the received signal, a filtered version of Z is given by:

Equation 34:

wherein h[n] can be any suitable FIR filter, e.g., designed as a high pass clutter filter. Specifically, in the method of Fig. 4 that involves eigenvalue decomposition, any suitable eigen-based clutter filter can be used as h[n] in Equation 34, e.g., designed using methods as described by Alfred and Lovstakken, in "Eigen-based clutter filter design for ultrasound color flow imaging: a review," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, volume 57, number 5, 2010.

Adopization, also referred to as "windowing" is typically required for reducing spectral lobes resulting due to processing finite length intervals of the received signal. In some embodiments, instead of applying windowing in the time domain, the Doppler analyzer applies windowing to the signal Z in the correlation domain as:

Equation 35:

wherein denotes the windowed signal, is the autocorrelation function of

the underlying window and the operator '·' denotes an element-wise multiplication.

The embodiments described above are given by way of example, and other suitable embodiments can also be used.

Although the embodiments described herein mainly address efficient Doppler imaging, the methods and systems described herein can also be used in other applications, such as in Direction of Arrival (DOA) estimation in applications such as radar and communications.

It will be appreciated that the embodiments described above are cited by way of example, and that the following claims are not limited to what has been particularly shown and described hereinabove. Rather, the scope includes both combinations and sub-combinations of the various features described hereinabove, as well as variations and modifications thereof which would occur to persons skilled in the art upon reading the foregoing description and which are not disclosed in the prior art. Documents incorporated by reference in the present patent application are to be considered an integral part of the application except that to the extent any terms are defined in these incorporated documents in a manner that conflicts with the definitions made explicitly or implicitly in the present specification, only the definitions in the present specification should be considered.

Claims

1. An apparatus for measuring target velocity, the apparatus comprising:

a transducer, configured to:

transmit into a target a set of pulses, which is sparse but has a recoverable power spectrum; and

receive a signal reflected from the target in response to the set of pulses; and a processor, configured to estimate a velocity of a selected region in the target by analyzing the reflected signal.

2. The apparatus according to claim 1 , wherein the set of pulses is based on a sequence containing, within an observation window, a smallest number of pulses for which the power spectrum is recoverable.

3. The apparatus according to claim 1 or 2, wherein:

(i) respective start times of the pulses are selected multiples of a time unit, and are spaced non-uniformly; and

(ii) for every possible multiple of the time unit up to a predefined maximal multiple, the set of pulses comprises at least one pair of pulses having start times that differ by that multiple.

4. The apparatus according to claim 1 or 2, wherein the set of pulses occupies a series of time intervals, such that in each time interval the set of pulses comprises two or more pulse sequences.

5. The apparatus according to claim 4, wherein the two or more pulse sequences have respective different Pulse-Repetition Frequencies (PRFs).

6. The apparatus according to claim 4, wherein the two or more pulse sequences comprise a first pulse sequence having a first PRF followed by a second pulse sequence having a second PRF.

7. The apparatus according to claim 4, wherein the two or more pulse sequences comprise at least first and second pulse sequences that at least partially overlap one another in time.

8. The apparatus according to claim 4, wherein the set of pulses comprises at least one pulse that is common to at least two of the pulse sequences.

9. The apparatus according to claim 1 or 2, wherein the set of pulses and the reflected signal comprise ultrasound signals.

10. The apparatus according to claim 1 or 2, wherein the processor is configured to estimate the power spectrum of the reflected signal, and to derive the velocity from the estimated power spectrum.

11. The apparatus according to claim 1 or 2, wherein the processor is configured to estimate the power spectrum by calculating a covariance of the reflected signal, to produce a measurement vector by averaging elements of the covariance corresponding to a common inter-pulse interval in the set of pulses, and to estimate the power spectrum from the measurement vector.

12. The apparatus according to claim 11, wherein the processor is configured to estimate the power spectrum over predefined Doppler frequencies by applying to the measurement vector a Fast Fourier Transform (FFT).

13. The apparatus according to claim 11, wherein the processor is configured to estimate the power spectrum by reordering the measurement vector as a Toeplitz matrix, and applying to the Toeplitz matrix eigenvalue decomposition.

14. The apparatus according to claim 13, wherein the processor is configured to estimate an actual number of Doppler frequencies of the power spectrum based on the eigenvalue decomposition.

15. The apparatus according to claim 11, wherein the processor is configured to estimate the power spectrum by applying a regularization functional to the measurement vector.

16. The apparatus according to claim 1 or 2, wherein the processor is configured to apply to the reflected signal a Finite Impulse Response (FIR) filter in a correlation domain.

17. A method for measuring target velocity, the method comprising:

transmitting into a target, by a transducer, a set of pulses, which is sparse but has a recoverable power spectrum;

receiving, by the transducer, a signal reflected from the target in response to the set of pulses; and

estimating a velocity of a selected region in the target by analyzing the reflected signal.

18. The method according to claim 17, wherein the set of pulses is based on a sequence containing, within an observation window, a smallest number of pulses for which the power spectrum is recoverable.

19. The method according to claim 17 or 18, wherein: (i) respective start times of the pulses are selected multiples of a time unit, and are spaced non-uniformly; and

(ii) for every possible multiple of the time unit up to a predefined maximal multiple, the set of pulses comprises at least one pair of pulses having start times that differ by that multiple.

20. The method according to claim 17 or 18, wherein the set of pulses occupies a series of time intervals, such that in each time interval the set of pulses comprises two or more pulse sequences.

21. The method according to claim 20, wherein the two or more pulse sequences have respective different Pulse-Repetition Frequencies (PRFs).

22. The method according to claim 20, wherein the two or more pulse sequences comprise a first pulse sequence having a first PRF followed by a second pulse sequence having a second PRF.

23. The method according to claim 20, wherein the two or more pulse sequences comprise at least first and second pulse sequences that at least partially overlap one another in time.

24. The method according to claim 20, wherein the set of pulses comprises at least one pulse that is common to at least two of the pulse sequences.

25. The method according to claim 17 or 18, wherein the set of pulses and the reflected signal comprise ultrasound signals.

26. The method according to claim 17 or 18, wherein estimating the velocity comprises estimating the power spectrum of the reflected ultrasound signal, and deriving the velocity from the estimated power spectrum.

27. The method according to claim 17 or 18, and comprising estimating the power spectrum by calculating a covariance of the reflected signal, producing a measurement vector by averaging elements of the covariance corresponding to a common inter-pulse interval in the set of pulses, and estimating the power spectrum from the measurement vector.

28. The method according to claim 27, wherein estimating the power spectrum comprises estimating the power spectrum over predefined Doppler frequencies by applying to the measurement vector a Fast Fourier Transform (FFT).

29. The method according to claim 27, wherein estimating the power spectrum comprises reordering the measurement vector as a Toeplitz matrix, and applying to the Toeplitz matrix eigenvalue decomposition.

30. The method according to claim 29, wherein estimating the power spectrum comprises estimating an actual number of Doppler frequencies of the power spectrum based on the eigenvalue decomposition.

31. The method according to claim 27, and comprising estimating the power spectrum by applying a regularization functional to the measurement vector.

32. The method according to claim 17 or 18, and comprising applying to the reflected signal a Finite Impulse Response (FIR) filter in a correlation domain.