Xampling in Ultrasound Imaging

Xampling in Ultrasound Imaging Noam W agner a , Y onina C. Eldar a,b , Arie F euer a , Gilad Danin a and Zvi F riedman c a T ec hnion-Israel Institute of T ec hnology , T ec hnion Cit y , Haifa, Israel; b Visiting Professor, Electrical Engineering Departmen t, Stanford, CA; c GE Health-care, T ec hnion Cit y , Haifa, Israel ABSTRA CT Recen t dev elopments of new medical treatment techniques put challenging demands on ultrasound imaging systems in terms of b oth image quality and ra w data size. T raditional sampling metho ds result in very large amoun ts of data, th us, increasing demands on pro cessing hardware and limiting the flexibility in the p ost- pro cessing stages. In this paper, w e apply Compressed Sensing (CS) tec hniques to analog ultrasound signals, follo wing the recen tly dev elop ed Xampling framework. The result is a system with significan tly reduced sampling rates which, in turn, means significan tly reduced data size while main taining the quality of the resulting images. Keyw ords: Array Pro cessing, Beamforming, Compressed Sensing, Dynamic F o cus, Finite Rate of Innov ation (FRI), Matrix P encil, Ultrasound, Xampling 1. INTRODUCTION Mo dern ultrasound systems utilize an array of transducer elements in a pro cess known as b eamforming. 1, 2 An imaging cycle b egins when mo dulated acoustic pulses are transmitted from some or all array elements. Specific dela ys are applied to the transmitting elemen ts, such that the in terfering w a v es form a narro w beam, along whic h most energy propagates. As the energy propagates, the b eam gets narrow er, until reaching the fo cal p oint, after whic h it expands. Determining the p osition of the fo cal p oint is ac hieved by applying appropriate dela ys to the transmitting elemen ts. As the fo cused energy propagates along the tissue, ec ho es are scattered and reflected b y density and propagation v elo city perturbations. 3 The array elements detect the reflected energy . A second b eamforming pro cess is then p erformed, aimed at lo calizing the reflecting elements, while at the same time improving signal to noise ratio (SNR). This b eamforming pro cess is p erformed by summing delay ed samples of the received data. F or that purp ose, high rate A/D con version must b e first carried out at each of the receiving channels. Actually , optimizing the image resolution requires that each receiving channel will b e sampled at 3-5 times the cen ter frequency of the mo dulated pulse. 1 More explicitly , let us consider a B-Mo de scan, in which the central frequency of the transducer may v ary in b etw een 2 − 15 M H z , dep ending on the use. 2 Adv anced comp osite materials, often used in the transducer, can attain a relativ e bandwidth in excess of 100%. Therefore, if we assume a nominal center frequency of 5 M H z , we end up with a baseband bandwidth of approximately 10 M H z . Confined to the classic Nyquist-Shannon sampling theorem, 4 where the only prior on the signal is that it is bandlimited, this implies that standard ultrasound devices must sample the analog signal received in each active element at a t ypical rate of at least 20 M H z . Recen t works 5–9 sho w that, b y exploiting other priors regarding the signal structure, it is p ossible to design more efficien t sampling schemes, whic h actually break through the Nyquist barrier. Decreasing the sampling rate is of v ast interest, as it may obviously b e manifested in reduction of machinery size and p ow er consumption. F urther author information: (Send corresp ondence to Noam W agner) Noam W agner: E-mail: noamw a@tx.tec hnion.ac.il, T elephone: +972(0)73 725 2653 Y onina C. Eldar: E-mail: yonina@ee.tec hnion.ac.il, T elephone: +972(4) 829 3256 Arie F euer: E-mail: feuer@ee.technion.ac.il, T elephone: +972(4) 829 4648 Gilad Danin: E-mail: daning@tx.technion.ac.il, T elephone: +972(4) 829 4706 Zvi F riedman: E-mail: zvi.friedman@med.ge.com In their recent work, T ur, Eldar and F riedman 10 first prop osed to implement these ideas in ultrasound imaging. They suggested to regard the signal received in each of the array transducer elemen ts as having finite rate of inno v ation (FRI). 5 More specifically , they assume eac h such signal comprises at most L replicas of a kno wn-shape pulse, all received within the time interv al [0 , τ ). Dela y and gain parameters are asso ciated with eac h replica, suc h that the entire problem may b e characterized by 2 L degrees of freedom. They then dev elop ed a new sub- Nyquist sampling technique that exploits this structure, in order to reduce the sampling rate wa y b eyond that whic h is used in standard ultrasound devices, dictated b y the classic Nyquist-Shannon sampling theorem. T ur, Eldar and F riedman’s work follows the spirit of analog compressed sensing, also referred to as Xampling. 6 The latter is an emerging framew ork, which com bines classic metho ds from sampling theory together with recent dev elopments in compressed sensing, aimed at sampling analog signals far b elo w the Nyquist rate. Throughout this pap er, we will use the term Xampling whenever referring to the sub-Nyquist sampling schemes purp osed in Ref. 10, and in later work by Gedalyah u, T ur and Eldar. 11 It should b e noted, though, that b oth these schemes are sp ecial cases of Xampling, in which the signal’s FRI prop ert y is exploited in order to achiev e the goal of sub-Nyquist sampling. The Xampling scheme suggested by T ur, Eldar and F riedman 10 comprises the following steps: first, the received signal is filtered using the compactly supp orted Sum of Sincs kernel. The filtered signal is then sampled at nearly the rate of inno v ation, namely the num ber of unknown parameters 2 L , which is muc h smaller than the Nyquist rate of the pulse. The extracted samples are used for computing a finite set of F ourier co efficien ts, whic h corresp ond to the τ -p erio dic extension of the receiv ed signal. Ha ving obtained the finite set of F ourier co efficien ts, sp ectral analysis techniques, similar to these presen ted by V etterli et al., 5, 12 are utilized, in order to estimate the set of unkno wn dela ys and amplitudes characterizing the reflected pulses. Whereas T ur, Eldar and F riedman 10 adopt a filtering and sampling approach as a preliminary step for obtaining the required set of F ourier co efficients, Gedalyah u, T ur and Eldar’s 11 later w ork suggests a second approac h for obtaining the same co efficients: the single filtering and sampling c hannel is replaced by a bank of mo dulators and in tegrators. Both approaches are aimed at minimal rate sampling of a single received channel. Referring to an arra y of transducer elements utilized in the ultrasound imaging device, by pro cessing the signal receiv ed in eac h channel separately (using either approach), one may obtain a corresp onding set of delays and amplitudes from low-rate samples. The sets obtained from all receiving elements may then b e combined (via some geometric in terpretation), in order to estimate the tw o-dimensional co ordinates of the reflecting elemen ts. Note, how ever, that such schemes cannot achiev e the SNR impro vemen t which is an in tegral part of standard b eamforming tec hniques; this is b ecause the correlation b etw een signals received in differen t channels is not exploited throughout the pro cess of extracting the parameters (pulse amplitudes and dela ys) from each signal. More explicitly , the Xampling schemes prop osed in Refs. 10 and 11 b oth aim at accurately detecting strong, lo calized pulses, related with macroscopic p erturbations. Ho wev er, the actual signals received b y the array elemen ts also contain comp onents which, in the con text of our work, may b e regarded as noise. Actually , part of these noisy comp onents arises from constructive and destructiv e in terference of acoustic wa v es reflected from dense, sub wa velength scatterers in the tissue (these are t ypically manifested as gran ular texture in the ultrasound image, called speckle, after a similar effect in laser optics 1 ). Apparen tly , these noisy components induce erroneous results when Xampling the received signals. If we wish to obtain meaningful results by either Xampling sc heme, while main taining a rather low o versampling factor, the ov erall SNR improv emen t is indeed a crucial step. This pap er is aimed at generalizing the sc hemes of Refs. 11 and 10 to multiple antenna arrays. Our goal is to obtain a tw o-dimensional, fo cused ultrasound image, corresp onding to strong p erturbations in the scanned plane, while reducing the sampling rate in eac h active element, by a factor of 10-15 times relative to the rate used in standard ultrasound devices. F urthermore, w e aim at achieving this goal in the presence of noise in the received signals. In such case, straightforw ard implementation of either schemes of Refs. 10 or 11 on eac h c hannel indep endently , w ould require hard thresholding (whic h in turn cancels/atten uates desired pulses) and/or increasing the o versampling factor, such that the final sampling rate gro ws tow ards the Nyquist rate. Our Xampling scheme’s most exp ensiv e computational comp onen t regards the extraction of the pulses’ dela ys and amplitudes from the set of low rate samples. Referring to standard ultrasound devices, this component substitutes both the Hilb ert transform, which is applied throughout the pro cess of env elope detection, 1 and the exp ensive b eamforming computations, namely: summing samples obtained at Nyquist rate from all active elemen ts. In addition, referring to straigh tforward implementation of either systems of Refs. 10 or 11 on eac h c hannel separately , our scheme extracts the pulses’ parameters once p er image line, rather than once p er active elemen t participating in the image line generation (which equiv alently means tens of times p er image line, dep ending on the n umber of active elemen ts used for b eamforming). The Xampling system w e propose may b e summarized as follo ws: Let us assume that w e could someho w generate the b eamformed signal, corresp onding to a single image line, in the analog domain. Assuming that such a signal main tains the FRI prop erty of the received signals from which it was constructed, we may Xample it using the sc heme suggested by either Ref. 10 or 11, yielding the dela ys and amplitudes of pulses along the final beamformed image line. By mathematically formulating the pro cess of Xampling the b eamformed signal using the scheme suggested in Ref. 11, and then applying several algebraic manipulations to this form ulation, we end up with a new set of generalized modulation kernels. The latter ma y now b e applied directly to the analog signals whic h are received in each of the active elemen ts, thus bypassing the impractical step of actually generating the b eamformed signal in the analog domain. Due to space limitation, the following paper is not aimed at presenting a rigorous review of our results. Instead, w e presen t an outline of our approach, and preliminary results obtained using actual ultrasound data. The pap er is organized as follo ws: Section 2 briefly outlines principles of standard ultrasound imaging, namely the pro cess of b eamforming in b oth P olar and Linear scan metho ds, applying dynamic receive (Rx) fo cus. Section 3 reviews the one-dimensional Xampling scheme suggested by Refs. 10 and 11. In section 4 we presen t our system, whic h com bines the concepts of b eamforming with those of one-dimensional Xampling. Section 5 provides results obtained by applying our Xampling sc heme up on actual ultrasound data. In this section we also analyze the reduction in the necessary amoun t of samples, and the wa y this affects the o verall computational cost. Finally , conclusions are dra wn in Section 6. 2. DYNAMIC FOCUSING IN POLAR AND LINEAR SCAN This section is aimed at outlining the method by whic h an ultrasound image is generated using a linear array of transducer elements. Our discussion refers mainly to B-mo de scan, in which the array simultaneously scans a plane through the b o dy , resulting in a tw o-dimensional image, whic h ma y be view ed on screen. W e are sp ecifically in terested in the pro cess carried out b y the electronic b eamformer, in whic h multiple signals, received from a set of receivers, are fo cused into a single trace, known as the b eamformed signal. The latter is env elope detected, forming a single image line. The analysis reviewed throughout this section is based mainly on Ref. 3. An ultrasound image consists of roughly 100 lines. Each line is generated throughout a single transmit-receiv e cycle. In such a cycle, a set of active elemen ts first transmits acoustic pulses, mo dulated to central frequency of 2 to 15 M H z , dep ending on the use. An appropriate delay is applied to the pulse transmitted from each transducer elemen t, aimed at obtaining constructive interference at a sp ecific p oint in the plane, referred to as the fo cal p oin t. The in terfering acoustic w av es form an acoustic pulse, whic h propagates along a narrow b eam (containing most of the transmitted energy). The b eam gets narrow er, until reaching the fo cal p oint, after which it expands. The velocity at whic h the pulse propagates will b e denoted by c , and v aries b etw een 1446 m/sec (fat) to 1566 m/sec (spleen). An a verage v alue of 1540 m/sec is assumed by scanners for pro cessing purp oses. As the transmitted pulse propagates inside the tissue, it encounters densit y and propagation-v elocity perturbations. These cause scattered and reflected echoes, which are detected by the array elements. Applying the acoustic recipro city theorem, 13 a second b eamforming pro cess is now carried out, in which the received signals are combined into a single trace, whic h in a sense, visualizes structures in the tissue, along the transmitted b eam. Mo dern ultrasound devices t ypically handle 64 to 192 transducer elemen ts in the beamforming process. W e emphasize, that com bining the receiv ed signals is p erformed in the digital domain, implying that modern ultrasound devices must first sample the signal received in eac h of the active elements at the Nyquist rate (typically 20 M H z ). Refs. 10 and 11, in tro duced an approach, which allow ed reconstruction of the analog signal detected by an individual transducer elemen t from a very low num b er of samples. In contrast, our generalized scheme is aimed at obtaining the low rate samples from all active elements, in a manner which will enable to directly reconstruct the b eamformed signal. This is further discussed in Section 4. W e would no w like to form ulate the manner in which the b eamformer com bines the signals receiv ed in all active elemen ts into a single trace, kno wn as the beamformed signal, and b etter understand the significance of the latter. This formulation applies to conv en tional ultrasound imaging, and will b e necessary when we translate our theoretical scheme of Xampling the b eamformed signal, into an applicable sc heme where the signals received in the activ e transducer elements are sampled directly . Referring to Figure 1, we examine the tw o-dimensional plane XZ in which 2 M + 1 elements are aligned along the ˆ x axis (the center of the array coincides with the origin). W e analyze one cycle, in which a single image line is constructed. The image line corresp onds to a b eam lo cated within the XZ plane, emerging from the array cen ter. W e denote b y α , the angle b etw een the b eam and ˆ z axis (normal to the array). The cycle b egins, when eac h activ e elemen t transmits a single modulated pulse, suc h that the in terference pattern ma y b e observ ed as a concentrated pulse of energy , propagating along the b eam. W e regard the pulse as if it w as transmitted from the arra y center, at a known time, which we shall define as t = 0. Knowing the sp eed in whic h the pulse propagates (denoted by c , and assumed 1540 m/sec ), we may now estimate the distance which it tra veled along the b eam by the time instance t n , denoted by r ( t n ), and thereby its tw o-dimensional p osition, p n : p n =  cos α sin α  T r ( t n ) =  cos α sin α  T ct n . (1) Assuming that an echo was reflected due to some p erturbation lo cated at p n , w e may e asily estimate the time in whic h it will arrive bac k at the origin, τ 0 ( t n ): τ 0 ( t n ) = t n + 1 c || p n || = 2 t n . (2) Since the distance from p n to each of the active elements v aries, each elemen t will detect the reflected pulse at a different time instance. Namely , the time at which the pulse arrives at the m th receiver, p ositioned at x m , is giv en by: τ m ( t n ) = t n + 1 c || p n − x m || = t n + 1 c q ( ct n sin α − δ m ) 2 + ( ct n cos α ) 2 , (3) where δ m denotes the x co ordinate of the receiving element. Let us denote b y ϕ m ( t ), the analog signal detected by the activ e element indexed m . Beamforming is now achiev ed b y shifting each of the received signals ϕ m ( t ), in order to comp ensate for the time difference τ m ( t n ) − τ 0 ( t n ), and then summing the shifted versions. The acoustic recipro city theorem 13 implies that when we sum the shifted signals, constructive interference will o ccur at t = τ 0 ( t n ) = 2 t n , providing that an echo was indeed reflected from p n at time t n . Denoting by Φ( t ; α ) the sum of the delay ed signals, we are thus interested in the v alue whic h Φ( t ; α ) obtains at τ 0 ( t n ) = 2 t n . Shifting the signal ϕ m ( t ) so that the difference τ m ( t n ) − τ 0 ( t n ) is comp ensated, is obtained b y applying the (p ossibly negativ e) delay θ m ( t n ; α ) = τ 0 ( t n ) − τ m ( t n ) = t n − q t 2 n + ( δ m /c ) 2 − 2 t n ( δ m /c ) sin α (4) to the signal receiv ed in the m th elemen t. Summarizing the ab o ve, we hav e: Φ(2 t n ; α ) = Φ( t ; α ) | t =2 t n = M X m = − M ϕ m ( t − θ m ( t n ; α )) | t =2 t n = M X m = − M ϕ m  t n + q t 2 n + ( δ m /c ) 2 − 2 t n ( δ m /c ) sin α  , (5) 0 m  1 m  1 m  mM  ˆ z ˆ x 1  mM   n p  ' th f oc al zo ne n Figure 1: 2 M + 1 elements aligned along the ˆ x axis. The ra y along which the pulse propagates forms an angle α with ˆ z axis. W e analyze a reflection emerging from the p oint p n , which is asso ciated with the n th fo cal zone. whic h corresp onds to the in tensity of a reflection originating at time t n , from the co ordinate p n . F or purp oses of con venience, we will finally substitute 2 t n → t , obtaining an expression for the b eamformed signal: Φ( t ; α ) = M X m = − M ϕ m  1 2  t + p t 2 + 4 ( δ m /c ) (( δ m /c ) − t sin α )   . (6) Referring to (1), our last substitution implies that Φ( t ; α ) represents the intensit y of a reflection originating at time t/ 2 from a p oint distanced ct/ 2 from the origin, along the transmitted b eam. The assumption is that the transmitted pulse indeed intersected this p oint at t/ 2, and was p ossibly scattered. Φ( t ; α ) is obviously the result of v arying the receive focal p oint along time, and its construction is therefore referred to as “Dynamic F o cusing”. W e note that, instead of obtaining the b eamformed signal using dynamic fo cusing pro cess, b eamformers some- times divide the image line in to N segmen ts, called fo cal zones, such that an entire segment of Φ( t ; α ), corre- sp onding to the n th fo cal zone, is constructed using a single set of delays (one dela y p er element). The set is obtained using (4), whic h is calculated for a single, represen tative p oint within the fo cal zone. In this section, we hav e formulated the manner in which the b eamformer combines the signals received in the transducer elements, { ϕ m ( t ) } M m = − M in to a b eamformed signal, Φ( t ; α ), through the pro cess of dynamic fo cusing. The dynamically fo cused, b eamformed signal, may now b e used in order to generate a single image line. The b eamformer p erforms the computation formulated in (6) in the digital domain, using samples obtained from eac h of the transducer elemen ts at the Nyquist rate. Our goal is to retrieve a set of parameters from which the b eamformed signal Φ( t ; α ) formulated in (6) may b e reconstructed, by sampling the received signals far b elow the Nyquist rate. W e achiev e this by exploiting an approximately FRI structure c haracterizing Φ( t ; α ), within the Xampling metho ds of Refs. 10 and 11, which are outlined in the next section. Throughout the rest of this paper, we will limit ourselv es to the case of linear scan, in which all b eams are parallel to the ˆ z axis. This is achiev ed by setting α = 0 in (6). The parameter δ m no w represen ts the distance b et ween the m th receiver and the b eam pro cessed at the curren t cycle. Equation (6) then b ecomes: Φ( t ; α = 0) = M X m = − M ϕ m  1 2  t + q t 2 + 4 ( δ m /c ) 2  . (7) 3. XAMPLING THE SIGNAL OBT AINED IN A SINGLE TRANSDUCER ELEMENT Up until this p oint, we hav e outlined the pro cess in which a B-mo de ultrasound image line is generated. The cycle b egins by transmitting a mo dulated pulse along a narrow b eam. The device then captures the intensit y of ec ho es reflected along the b eam by applying dynamically fo cused b eamforming. Regarding the ultrasound signal detected by a single transducer elemen t indexed m (denoted b y ϕ m ( t )), T ur, Eldar and F riedman 10 assume that it consists of a set of kno wn-shap e pulses, which result from reflections of the transmitted pulse b y strong, macroscopic p erturbations in the tissue. The signal can hence b e appro ximated as an FRI of the form: ϕ m ( t ) = L X l =1 a l,m h ( t − t l,m ) , (8) where h ( t ) is a known-shape pulse, and there exists some τ > 0 such that h ( t − t l,m ) = 0, ∀ t / ∈ [0 , τ ), l = 1 ...L . The extremely short supp ort of h ( t ) implies that ϕ m ( t ) is of v ery wide band (the supp ort of h ( t ) is typically 600 nsec , with τ b eing typically 200 µsec ). Classic Nyquist-Shannon sampling theorem thus forces standard ul- trasound devices to sample ϕ m ( t ) at a high rate (t ypically 20 M H z ). Nevertheless, we may easily observe that ϕ m ( t ) actually has only 2 L degrees of freedom ( L b eing the n umber of macroscopic scatterers along the path of the transmitted pulse). The schemes in tro duced by Refs. 10 and 11 manage to exploit this prop ert y , enabling the reconstruction of ϕ m ( t ) from a m uch smaller n umber of samples (at least 2 L p er time interv al τ ). In this section we outline the system suggested by Ref. 11 for Xampling the signal ϕ m ( t ). This approach is the basis for our scheme, introduced in Section 4. W e emphasize that b oth Refs. 11 and 10 treat Xampling of a one dimensional signal which is received in a single transducer element. Moreov er, they do not treat integration of samples, obtained from multiple elements, into a tw o-dimensional ultrasound image. Our no velt y concerns a metho d for obtaining parameters of the b eamformed signal , Φ( t ; α = 0) (introduced in the previous section), whic h is directly related to the image line. W e obtain these parameters from low rate samples of the individual signals ϕ m ( t ). By extracting the parameters of the b eamformed signal w e also cop e with the noisy comp onents, whic h induce erroneous results when attempting to reconstruct the individual signal ϕ m ( t ), from its corresp onding lo w rate samples. Let us denote by H ( ω ) the CTFT of the kno wn-shap e pulse h ( t ), and by φ m [ k ] the k th F ourier co efficient of ϕ m ( t )’s τ -p erio dic extension. F urther denote by κ , a set of K consecutiv e indices for whic h H ( ω = 2 π t k ) 6 = 0, ∀ k ∈ κ . Ref. 10 sho ws that, as long as K ≥ 2 L and the unknown time delays are distinct, i.e. t i,m 6 = t j,m , ∀ i 6 = j , one may accurately estimate ϕ m ( t ), from the set { φ m [ k ] } k ∈ κ . Gedaly ah u, T ur and Eldar 11 suggest a practical approach for obtaining the set { φ m [ k ] } k ∈ κ , inv olving a bank of mo dulators and integrators. Referring to Figure 2, after ha ving chosen the set of indices κ , we construct p ≥ K branches, each comprising a modulating k ernel and an integrator. W e then set the mo dulating kernels to b e: s q ( t ) = P k ∈ κ  s q ,k e − j 2 π τ kt  , q = 1 , 2 , ..., p . (9) Ref. 11 pro ves that the follo wing relation holds: c = S φ, (10) where S is a p × K matrix with s q ,k as its ( q , k ) element, c denotes the length- p sample vector with the output of the q th branc h as its q th element, and φ denotes the length- K vector with the F ourier co efficient φ m [ k ] as its k th element. As long as S has full column rank, we can recov er φ from the samples by φ = S † c . Denote by H the K × K diagonal matrix with k th en try H ( ω = 2 π τ k ), k ∈ κ , and b y V ( t ) the K × L matrix with ( k , l ) element e − j 2 2 π τ kt l,m , where t = { t 1 ,m , ..., t L,m } is the vector of unknown pulse delays received in the individual transducer elemen t. In addition denote by a the length- L vector whose l th element is a l,m . Then: φ = HV ( t ) a . (11)   2 11 j kt k k s t s e          0 1 dt        m t    2 j kt p pk k s t s e        1 c   0 1 dt     p c  Figure 2: Gedaly ahu, T ur and Eldar’s 11 m ultichannel Xampling scheme. The resulting samples { c q } p q =1 constitute a mixture of F ourier co efficients corresp onding to ϕ m ( t )’s τ -p erio dic extension. The co efficients’ indices b elong to the set κ . Figure reprinted from Ref. 11. The matrix H is inv ertible b y construction. Generating the length- K vector y by left multiplying φ by H − 1 , we ha ve: y = V ( t ) a , (12) whic h is a standard problem of finding frequencies and amplitudes of a sum of L cisoids (complex sin usoids). The time-dela ys { t l,m } L l =1 ma y b e estimated using nonlinear tec hniques (e.g. annihilating filter, 14 or matrix p encil 15 metho ds). Having obtained the time delays, estimating the amplitudes { a l,m } L l =1 is a linear problem, whic h ma y b e easily solv ed using a least squares approac h. In the next section w e apply the sc heme of Figure 2 on the beamformed signal Φ( t ; α = 0), in order to reconstruct it from a small num b er of its samples. Recall that standard ultrasound devices digitally construct Φ( t ; α = 0), after sampling the individual signals received in the transducer elements at the Nyquist rate. Since our goal is to break the Nyquist barrier, w e b ypass the actual construction of Φ( t ; α = 0), by translating its Xampling to a sc heme which may b e applied directly on the analog signals { ϕ m ( t ) } M m = − M . 4. GENERA TING A 2D IMAGE BY XAMPLING THE BEAM FORMED SIGNAL A t the basis of our approac h is the assumption that the beamformed signal Φ( t ; α = 0) maintains the FRI prop ert y whic h characterizes the signals from whic h it is constructed. This prop ert y w as formulated in (8). More explicitly , Φ( t ; α = 0) may b e written in the following manner: Φ( t ; α = 0) = L X l =1 b l h ( t − t l ) . (13) This claim requires justification, which will not b e provided within the scop e of this pap er. The justification is concealed within the fact that the nonlinear scaling of ϕ m ( t ), formulated in (7), has little affect on the shap e of the pulses whic h construct it, due to their extremely short supp ort with resp ect to τ . This notion is well demonstrated in Figure 3. A single pulse is transmitted along the narro w b eam extending from the origin, along the ˆ z axis. Tw o scattering elements are illuminated b y the b eam (one is distanced 1 cm from the origin and the second is distanced 2 cm from the origin). As the pulse interacts with the elements, reflections are scattered and receiv ed in each of the 16 array elements. The traces well demonstrate the differen t 8 m  ˆ z ˆ x 1 m  16 m  9 m  1 cm z  2 cm z  Reflecting Elements 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 t [ µ sec] Receiver Index Received Traces 0 5 10 15 20 25 30 35 0 2 4 6 8 10 12 14 16 t [ µ sec] Receiver Index Received Traces after Applying Distortion Operator 0 5 10 15 20 25 30 35 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 t [ µ sec] Φ (t; α =0) Generated by Summing the Distorted Signals Figure 3: Generation of Φ( t ; α = 0) as defined in (7), by summing distorted versions of the signals received b y 16 elements. The setup is depicted in the top left: t w o scattering elements are p ositioned along the ˆ z axis, reflecting the pulse transmitted by an array of 16 elements. The 16 traces received by the elements are depicted in the top righ t. Notice that the reflected pulses are not aligned due to the differen t path trav eled to each elemen t. After applying the distortion suggested in (7), we obtain mo dified traces, in which the pulses app ear aligned, at the cost of slight distortion to the shap e of each replica (b ottom left). Finally , summing the distorted traces yields the signal Φ ( t ; α = 0) (b ottom right). The latter may b e approximated as a delay ed sum of tw o kno wn-shap e pulses. pulse dela ys obtained in eac h of the elements. Simply summing the 16 traces (namely - b eamforming with receiv e fo cus set to infinity) may yield certain SNR improv emen t, yet artifacts will b e formed due to the fact that the corresponding pulses are not aligned. Using standard imaging techniques, these artifacts are manifested in a non-focused image. Y et they may hav e ev en more profound implications when attempting to Xample the b eamformed signal. Let us now distort each of the received traces, parametrized b y m , as defined in (7): ˆ ϕ m ( t ) = ϕ m  1 2  t + q t 2 + 4 ( δ m /c ) 2  . (14) Observing Figure 3, one may notice that the pulses are no w aligned, although each may hav e undergone a slight distortion. Summing the distorted signals, we now obtain Φ ( t ; α = 0) which ma y evidently b e regarded as FRI. If Φ ( t ; α = 0) actually existed in the analog domain, then w e could Xample it using the scheme describ ed in Section 3, since it approximately satisfies (13). Namely , we could reconstruct it from a rather small subset of F ourier co efficients corresp onding to its τ -p erio dic extension. Let us no w assume that Φ ( t ; α = 0) existed in the analog domain, and feed it to the input of the scheme depicted in Figure 2, instead of the individual trace ϕ m ( t ). F orm ulating the op eration of a single branch, with Φ( t ; α = 0) as its input, w e hav e: c q = 1 τ Z τ 0 ( X k ∈ κ  s q ,k e − j 2 π τ kt  ) ( M X m = − M ϕ m  1 2  t + q t 2 + 4 ( δ m /c ) 2  ) dt. (15) After several algebraic manipulations, and exchange of v ariables inside the integral, (15) may b e brough t into a rather straigh tforward form: c q = M X m = − M ( 1 τ Z ˆ τ 0 ˆ s q ,m ( t ) ϕ m ( t ) dt ) , (16) where ˆ s q ,m ( t ) is defined as: ˆ s q ,m ( t ) ,  1 + ( δ m ct ) 2  ( X k ∈ κ s q ,k e − j 2 π τ k ( t − 1 t ( δ m c ) 2 ) ) u ( t − | δ m c | ) , (17) for 1 ≤ q ≤ p , − M ≤ m ≤ M , and u ( · ) is the unit step function: u( x ) =  1 x ≥ 0 0 else . (18) The form ulation of (16) may no w b e interpreted by the following algorithm, applied on the set { ϕ m ( t ) } M m = − M : 1. Having defined the full column-rank p × K matrix S , generate the extended set of p × (2 M + 1) mo dulating k ernels ˆ s q ,m ( t ), defined in (17). 2. Mo dulate each of the analog signals received by the active elements ϕ m ( t ), m ∈ {− M , ..., M } , using its corresp onding, size p , set of kernels, { ˆ s q ,m ( t ) } p q =1 , yielding a corresp onding size p set of co efficients: c q ,m = 1 τ R ˆ τ 0 ˆ s q ,m ( t ) ϕ m ( t ) dt q ∈ { 1 , ..., p } . (19) Applying the ab ov e step up on each of the signals ϕ m ( t ), m ∈ {− M , ..., M } yields a p × (2 M + 1) matrix of output samples. Note that the upp er integration b ound was mo dified to ˆ τ , where: ˆ τ = max m ∈{− M ,...,M }  1 2 ( τ + p τ 2 + 4( δ m /c ) 2 )  . (20) 3. Sum the samples along m , yielding the single, length- p vector, c , of whic h the q th element is c q , satisfying: c q = M X m = − M c q ,m . (21) 4. Obtain the vector y by: y = H − 1  S † c  . (22) Since we b egan the deriv ation b y injecting Φ( t ; α = 0) into the scheme depicted in Figure 2, then the v ector reco vered by S † c , holds, as its k th elemen t, the k th F ourier coefficient corresponding to the τ - p erio dic extension of the b eamformed signal Φ( t ; α = 0). 5. Solve the problem formulated in (12), extracting the 2 L unkno wns of (13). Here V ( t ) is the K × L matrix with ( k , l ) elemen t e − j 2 2 π τ kt l , t = { t 1 , ..., t L } is the vector of unknown pulse delays, and a is the length- L vector whose l th element is b l . The time dela ys are extracted by applying nonlinear techniques (e.g. annihilating filter, or matrix p encil metho ds). Estimating the amplitudes { b l } L l =1 is then a linear problem. Summarizing the ab ov e, Figure 4 sc hematically depicts the suggested generalized Xampling scheme. The active arra y elemen ts t ypically lie symmetrically with resp ect to the processed image line, resulting in symmetric k ernels, with resp ect to m (i.e. ˆ s q , − m ( t ) = ˆ s q ,m ( t )). Exploiting this symmetry enables reducing the total num ber of samples, by summing eac h pair ϕ m ( t ) and ϕ − m ( t ) prior to the mo dulation and integration branc h. This concept is implemen ted in the figure. Recall that the transmitted (and reflected) pulses are modulated to a high frequency carrier (t ypically f c = 5 M H z ), such that most of the pulse energy is concen trated far from the DC frequency . In order to a void singularit y of the matrix H it is necessary to pick the set κ , such that the frequencies  ω k = 2 π τ k  k ∈ κ are near 2 π f c . This implies, that when implementing the Xampling scheme depicted in Figure 2 using real filters, the set κ is doubled, so that if the index k is in the set κ , we m ust also add − k to κ . With the condition that p ≥ | κ | (so that S is full column rank), | · | denoting cardinality , we are forced to double the num ber of samples: for Φ( t = 0; α ) ha ving 2 L degrees of freedom ( { b l , t l } L l =1 ), reconstruction requires a minimum of 2 K samples within the interv al [0 , τ ), where K ≥ 2 L as b efore. Indeed, the actual ultrasound signal Xampled in Ref. 10 w as first demo dulated, yielding a “complex” analog signal. This practically means that the single, mo dulated real signal, w as split into tw o base-band signals, eac h Xam pled at least at the Rate of Innov ation - thereby the num b er of samples obtained within the interv al [0 , τ ) satisfied p ≥ 2 K ≥ 4 L . T o conclude, we Xample the signals { ϕ m ( t ) } M m = − M in a manner which enables reconstruction of the beamformed signal Φ( t ; α = 0), from which an image line ma y b e directly generated: the individual delay t l corresp onds to a co ordinate, lo cated along the line normal to the set of activ e elements, at a distance ct l 2 from its center. W e ma y simply set the pixel asso ciated with this co ordinate an intensit y prop ortional to b l . Alternatively , w e ma y generate a trace, by con volving the stream of pulses P L l =1 b l δ ( t − t l ) with the env elop e of h ( t ), and then set the pixel intensities along the corresp onding image line accordingly . The latter metho d was used for the sim ulations discussed in the next section, in order to obtain results comparable to those obtained by standard imaging tec hniques. F or our sim ulations, we assumed L = 30 macroscopic scatterers along a single b eam. Defining ρ ≥ 1 to be the o v ersampling factor, we then need at least K = 2 ρL consecutiv e indices, for which H ( ω = 2 π τ k ) 6 = 0. W e select the set of indices, κ , such that ω k = 2 π τ k are near 2 π f c , where f c is the frequency of the carrier wa v e (appro ximately 5 M H z ). In order to obtain real kernels, w e add the opp osite indices to the set κ , such that we finally hav e | κ | = 4 ρL . W e no w construct a scheme with p channels per element, − M ≤ m ≤ M , suc h that p = | κ | = 4 ρL . W e do this by choosing the matrix S in a very straightforw ard manner: S p × p =  1 2 I 1 2 I 1 2 j I − 1 2 j I  , (23) where I is a p 2 × p 2 iden tity matrix. Referring to (17), our selection of S forms the following set of p × (2 M + 1) k ernels: ˆ s q ,m ( t ) = (  1 + ( δ m ct ) 2  cos  2 π τ k q  t − 1 t ( δ m c ) 2  u  t − | δ m c |  1 ≤ q ≤ p 2 −  1 + ( δ m ct ) 2  sin  2 π τ k q − p 2  t − 1 t ( δ m c ) 2   u  t − | δ m c |  p 2 + 1 ≤ q ≤ p , (24) whic h were utilized in the sc heme depicted in Figure 4. W e c hose the matrix pencil 15 metho d in order to estimate the set of dela ys from the estimated F ourier co efficien ts obtained using our scheme. One significant adv an tage of the matrix p e ncil scheme is that it provides the ability to estimate the actual num ber of delays concealed within the noisy data, based on a SVD decomp osition pro cess. Recall that in real imaging we ha v e no prior knowledge regarding the n um b er of reflecting elements aligned along a single image line. Nevertheless, when we design the matrix p encil we must determine the pencil parameter, whic h w e shall denote by η . The matrix p encil metho d assumes η to b e greater than (or equal to) the num b er of complex exp onentials comprising the estimated signal (in our context, this is equiv alen t to the num ber of dela yed pulses, L ). As a result, when designing the matrix p encil blo ck we must first determine an upp er b ound on the n umber of reflected elements along a single image line and then set the p encil parameter η accordingly . This farther dictates a lo wer b ound on the num b er of samples used for estimating the delays (denoted here by K ). Summarizing the ab ov e constrain ts, the following must hold: L ≤ η ≤ K − L, (25) where L is an estimated upp er b ound on the num b er of reflecting elements. The need to guaran tee a non-empty in terv al [ L, K − L ], justifies the requirement that K ≥ 2 L .   1, ˆ M st   ˆ 0 1 dt        M t     1 M t     2 M t     0 t    2 M t     1 M t     M t       2, ˆ M st    , ˆ pM st    1, 1 ˆ M st     2, 1 ˆ M st     ,1 ˆ pM st     1, 2 ˆ M st     2, 2 ˆ M st     ,2 ˆ pM st     1,0 ˆ st    2,0 ˆ st    ,0 ˆ p st     0 c 1 c p c   ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt       ˆ 0 1 dt     Figure 4: Generalized Xampling sc heme yielding p co efficients, used for reconstructing Φ( t ; α = 0) from 2 M + 1 receiving elemen ts. 5. RESUL TS In the following section, we examine the result of applying our suggested Xampling scheme up on actual ra w RF ultrasound data. The data was acquired using a programmable imaging system (Mo del V-1-128, V erasonics, Inc., Redmond, W A), equipp ed with a 128-element 1-D linear transducer array (Mo del L7-4, Philips Healthcare, Bothell, W A). The imaging target w as a commercial multi-purpose gray-scale phantom (Mo del 403GS LE, Gammex, Inc., Middleton, WI) including 0.1-mm nylon wires embedded in tissue mimicking material. W e compare the follo wing images: 1. Standard Image - Ultrasound image generated using standard imaging tec hnique: the frequency of the carrier wa v e w as 5 . 142 M H z ; data was acquired at high rate (20 M H z ); the num b er of activ e transducer elemen ts used for fo cusing the b eam at transmit (Tx) and at receive (Rx) was 16. Setting the maximum imaging depth to 7 . 88 cm , the num b er of samples used in order to generate a single image line was th us 16 × 2048 = 32 , 768. W e used dynamic Rx fo cusing, comprising K = 100 fo cal zones. 2. Xampled Image, Dynamic Rx F o cusing - Ultrasound image generated using the exact scheme suggested in Section 4. 3. Xampled Image, Infinit y Rx F o cusing - Ultrasound image generated using the scheme suggested in Section 4, forcing δ m = 0, ∀ m ∈ {− M , ..., M } . Observing (7), one can easily see that this degenerate implementation forces Rx fo cusing to b e infinity: Φ ( t ; α ) is simply the sum of all signals obtained from the activ e elements, without applying an y delay . The resulting images are arranged in Figures 5-9, according to the following table: Figure Image T yp e Estimated No. of No. Elemen ts/Oversampling F actor 5 Standard Irrelev ant 6 Xampled (Dyn. F o cus and ∞ F o cus) 30/1 7 Xampled (Dyn. F o cus and ∞ F o cus) 30/2 8 Xampled (Dyn. F o cus and ∞ F o cus) 30/3 9 Xampled (Dyn. F o cus and ∞ F o cus) 30/4 T able 1: Prop erties of images display ed in Figures 5-9 Figure 5: Standard image comprising 113 Image Lines, generated with 100 fo cal zones. Figure 6: Xampled Images with L = 30 and ρ = 1. Dynamic F o cusing (left) and F o cus at Infinity (right). Figure 7: Xampled Images with L = 30 and ρ = 2. Dynamic F o cusing (left) and F o cus at Infinity (right). Figure 8: Xampled Images with L = 30 and ρ = 3. Dynamic F o cusing (left) and F o cus at Infinity (right). Figure 9: Xampled Images with L = 30 and ρ = 4. Dynamic F o cusing (left) and F o cus at Infinity (right). A t this p oint of our work we av oided deriving quan titative measures to the estimation quality; this is b ecause w e had no actual do cumentation regarding the true p ositions and in tensities of the reflecting elemen ts. Using SNR measurements with the “Standard Image” as reference may b e very deceiving, since the latter is obviously noisy b y itself - noise which was often eliminated by the Xampling sc heme, due to the fact that it extracts isolated delta excitations. The reader may receive qualitative impressions, observing the similarity b et ween the Xampled images and the “Standard” image, and the obvious difference b etw een a dynamically fo cused image and its equiv alent infinity fo cused one. In T able 2 we provide an estimate of: a) the num b er of samples required for generating a single image line using our suggested Xampling scheme (referring to L and ρ used in Figures 6-9) and using standard imaging tec hniques; b) the estimated computational cost of generating the same image line for b oth metho ds: Image T yp e Estimated Ov ersampling K Sampling Rate Cost No. of F actor [No. of Samples/ [MegaOps./Line] Reflectors Element/Line] Xampled 30 1 60 120 0.43 30 2 120 240 2.81 30 3 180 360 9.06 30 4 240 480 21.05 Standard Irrelev ant 2048 0.06 T able 2: Sampling rates and computational costs estimated for Xampled imaging (using v arious ov er-sampling factors) and for standard imaging. The computational cost of a Xampled image line roughly grows like O  | κ | 3  . The computational cost required for the standard imaging pro cess (assuming 2048 samples p er image line, and 16 activ e receiv ers used for b eamforming) comprises 2048 × 15 Add op erations, and additionally the cost of the Hilb ert transform which is utilized for env elope detection (typically implemented using tw o FFT op erations). T able 3 which app ears in App endix A details the blo c ks used for estimating the Xampling sc heme computational costs. 6. CONCLUSIONS AND FUTURE WORK This work fo cused on generalizing the Xampling metho d suggested in Ref. 11 to an array of m ultiple receiving elemen ts. A t the heart of our generalization is the observ ation, that the dynamic fo cusing and the filtering part of the Xampling can b e combined into a set of mo dulating kernels and p erformed directly on the analog signals. This, in turn, can b e sampled at a rate w ay b elow the Nyquist rate. Preliminary tests on actual ultrasound data yield results whic h are quite similar to an image obtained using standard techniques, while reducing the sampling rate by a factor of 5 to 15. Apparently , this is achiev ed at the cost of increased computational effort. By reducing the sampling rate, we hop e to simplify the front end hardw are (mainly in terms of size and p o wer consumption) while maintaining image quality . A CKNO WLEDGMENTS The authors would like to thank Dr. Omer Oralk an and Prof. Pierre Khuri-Y akub of the E. L. Ginzton Lab oratory at Stanford Univ ersity , for providing the RF ultrasound data and for many helpful discussions. APPENDIX A. ESTIMA TION OF XAMPLING SCHEME COMPUT A TIONAL COST Op eration Details Cost/ O ( · ) Sum Outputs S um  c L × (2 M +1) , 2  ( L − 1) × (2 M + 1) Matrix P encil Y K × 1 = A K × p c p × 1 K × P Obtain { t l } L l =1 S V D  Q 2 3 K × 1 3 K  2 3 K ×  1 3 K  2 Y 2 3 K × 1 3 K = U 2 3 K × L Σ 0 L × L V 0 T L × 1 3 K 2 3 K × L × L + L × L × 1 3 K Y 2 3 K × 1 3 K = U 2 3 K × L Σ 0 L × L V 0 T L × 1 3 K 2 3 K × L × L + L × L × 1 3 K pinv  Y 2 3 K × 1 3 K   2 3 K  3 +  1 3 K  3 Q 1 3 K × 1 3 K = A 1 3 K × 2 3 K B 2 3 K × 1 3 K 2 3 K ×  1 3 K  2 eig  Q 1 3 K × 1 3 K   1 3 K  3 Least Squares L K × L ? M K × L K × L Obtain { b l } L l =1 pinv ( V K × L ) K 3 + L 3 b L × 1 = A L × K c K × 1 K × L T able 3: Blo cks used for calculating Xampling scheme computational cost using the matrix p encil metho d. REFERENCES [1 ] Szab o, T. L., “Diagnostics ultrasound imaging: Inside out,” in [ A c ademic Pr ess Series in Biome dic al Engi- ne ering ], Bronzino, J., ed., Ch. 7, 10, Elsevier Academic Press, 200 Wheeler Road, 6th Flo or, Burlington, MA 01803, USA, first ed. (2004). [2 ] Jensen, J. A., “Ultrasound imaging and its mo deling,” T opics in Applie d Physics 84 , 135 – 165 (2002). [3 ] Jensen, J. A., “Linear description of ultrasound imaging systems.” Notes for the In ternational Summer Sc ho ol on Adv anced Ultrasound Imaging, T echnical Universit y of Denmark (1999). [4 ] Shannon, C. E., “Comm unication in the presence of noise,” Pr o c. IRE 37 , 10 – 21 (1949). [5 ] V etterli, M., Marziliano, P ., and Blu, T., “Sampling signals with finite rate of innov ation,” IEEE T r ansac- tions on Signal Pr o c essing 50, No. 6 , 1417 – 1428 (2002). [6 ] Mishali, M., Eldar, Y. C., Dounaevsky , O., and Shoshan, E., “Xampling: Analog to digital at sub-nyquist rates,” IET Journal of Cir cuits, Devic es and Systems 5, Issue 1 , 8 – 20 (2011). [7 ] Mishali, M. and Eldar, Y. C., “Sub-n yquist sampling: Bridging theory and practice,” submitte d (April 2010). [8 ] Mishali, M. and Eldar, Y. C., “F rom theory to practice: Sub-nyquist sampling of sparse wideband analog signals,” IEEE Journal of Sele cte d T opics on Signal Pr o c essing 4, No. 2 , 375 – 391 (2010). [9 ] Matusiak, E. and Eldar, Y. C., “Sub-nyquist sampling of short pulses: P art i,” submitte d to IEEE T r ans. Information The ory [Online] arXiv 1010.3132 (2010). [10 ] T ur, R., Eldar, Y. C., and F riedman, Z., “Innov ation rate sampling of pulse streams with application to ultrasound imaging,” to app e ar in IEEE T r ansactions on Signal Pr o c essing (2010). [11 ] Gedalyah u, K., T ur, R., and Eldar, Y. C., “Multic hannel sampling of pulse streams at the rate of innov ation,” ac c epte d to IEEE T r ans. on Signal Pr o c essing (2010). [12 ] Blu, T., Dragotti, P . L., V etterli, M., Marziliano, P ., and Coulot, L., “Sparse sampling of signal inno v ations,” IEEE Signal Pr o c ess. Mag. 25, No. 2 , 31 – 40 (2008). [13 ] Kinsler, L. E., F rey , A. R., Copp ens, A. B., and Sanders, J. V., [ F undamentals of A c oustics ], John Wiley and Sons, New Y ork, third ed. (1982). [14 ] Stoica, P . and Moses, R., [ Intr o duction to Sp e ctr al Analysis ], Prentice-Hall, Englewoo d Cliffs, NJ (2000). [15 ] Sark ar, T. K. and Pereira, O., “Using the matrix p encil metho d to estimate the parameters of a sum of complex exp onen tials,” IEEE Antennas and Pr op agation Magazine 37, No. 1 .