Patent Publication Number: US-8532951-B2

Title: Method for calibrating a transducer array

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     This application claims the benefit of U.S. Provisional Application No. 61/289,154, filed on 22 Dec. 2009, which is incorporated in its entirety by this reference. 
    
    
     TECHNICAL FIELD 
     This invention relates generally to the transducer field, and more specifically to an improved method for calibrating a transducer array. 
     BACKGROUND 
     In an array of transducers that emit and/or receive signals, each transducer may be characterized by a number of parameters that affect overall measurements of signals passed among the transducers, and the processing of such signals. For instance, in an array of N transducers having positions denoted as x i  (i=0, 1, . . . , N−1), each transducer may be assumed to have unknown but fixed emission and reception delays (denoted by e i  and r i , respectively). These and other calibration parameters affect overall time delay measurement m i,j , or the time elapsed between when a probing signal sent to an emitter and when a signal is recorded at a receiver. In various applications of transducer arrays, the accuracy of acquired data and the processing of such data rely on accurate estimates through calibration parameters characterizing the transducers. However, in practice when calibrating transducer arrays using time delay measurements, only a noisy version of these time delay measurements m i,j  (i.e., {circumflex over (m)} i,j ) is obtainable or accessible. The origin of the noise is manifold, but includes electronic noise, non-ideal sensor characteristics and in-homogeneities in the propagation medium. The measurement {circumflex over (m)} i,j  can also be considered as missing when, for example, the signal-to-noise ratio of the output signal is too weak to provide a relevant estimate, or if any transducers are faulty, thereby leading to a faulty time delay measurement. These and other complications typically inhibit accurate and robust calibration of the transducer array. Thus, there is a need in the transducer field to create an improved method for calibrating a transducer array using time delay measurements. This invention provides such an improved method for calibrating a transducer array. 
    
    
     
       BRIEF DESCRIPTION OF THE FIGURES 
         FIG. 1  is a schematic of the method for calibrating a transducer array of a preferred embodiment; 
         FIG. 2  is a detailed schematic of the step of forming a time delay measurement data set in the method of a preferred embodiment; 
         FIG. 3  is a detailed schematic of the step of modifying the time delay measurement data set in the method of a preferred embodiment; 
         FIG. 4  is a detailed schematic of the step of modifying estimated positions of transducers in the method of a preferred embodiment; 
         FIG. 5  is a detailed schematic of the step of determining a test value to compare to a first stopping criterion in the method of a preferred embodiment; 
         FIG. 6  is a detailed schematic of the step of denoising the distance data set in the method of a preferred embodiment; 
         FIG. 7  is a detailed schematic of the step of determining a test value to compare to a second stopping criterion in the method of a preferred embodiment; 
         FIG. 8  is a schematic of the system for calibrating a transducer array of a preferred embodiment; 
         FIG. 9  is a measurement matrix obtained with a ring transducer with a plurality of transducer elements in an example embodiment of the method; 
         FIGS. 10A and 10B  are plots illustrating the comparison between estimated transducer positions and assumed transducer positions, according to distance from the center of a transducer ring and angle between adjacent transducer elements, respectively, where the estimated transducer were obtained in an example embodiment of the method and the assumed transducer positions are based on known approximate geometry of the transducer ring; and 
         FIGS. 11A and 11B  are plots illustrating the comparison between a first set of estimated transducer positions and a second set of estimated transducer positions, according to distance from the center of a transducer ring and angle between adjacent transducer elements, respectively, where the first and second sets of estimated transducer positions were obtained in first and second instances, respectively, of an example embodiment of the method. 
     
    
    
     DESCRIPTION OF THE PREFERRED EMBODIMENTS 
     The following description of preferred embodiments of the invention is not intended to limit the invention to these preferred embodiments, but rather to enable any person skilled in the art to make and use this invention. 
     In a preferred embodiment, as shown in  FIG. 1 , the method  100  for calibrating a transducer array characterized by calibration parameters, based on time delay measurements between transducer pairs, includes the steps of: receiving an estimate of the positions of the transducers S 110 ; forming a time delay measurement data set S 120  that includes time delay measurements between transducer pairs, including non-faulty time delay measurements and any faulty time delay measurements; generating from the time delay measurement data set an estimate of at least one calibration parameter S 130 ; modifying the time delay measurement data set S 140  including modifying the faulty time delay measurements, based on the estimated calibration parameters; modifying the estimated positions of the transducers S 150  based on the estimated calibration parameters; determining a test value for comparing to a stopping criterion S 210 ; and step S 220 , which includes repeating steps S 130  through S 210  until a stopping criterion is met. The method  100  may further include the step of denoising a distance data set or testing matrix S 160  representing the modified estimated positions of the transducers, and modifying estimated transducer positions based on the denoised distance data set S 190 . The method  100  preferably further includes providing an output of at least one modified estimated calibration parameters S 230  and/or storing at least one of the modified estimated calibration parameters S 240 . In one preferred embodiment, the method stores or provides an output of estimates of one or more of the following calibration parameters: respective emission delays for each transducer in the transducer array, respective reception delays for each transducer in the transducer array, positions of the transducers in the transducer array (relative and/or absolute), and speed of signal propagation. The emission, reception, and/or transducer positions may be presented, for example, in the form of a vector or a table. The method  100  provides an accurate and robust array calibration using time delay measurements, using an iterative masking approach to recover missing or faulty data obtained for calibration. The time delay measurements are preferably obtained through a test operation of the transducer array, but may alternatively be retrieved from storage, received from a user, and/or received in any suitable manner. The method  100  is also robust to noise and model mismatch, including a unique iterative denoising algorithm to further improving accuracy of the calibration. 
     The method  100  is used to calibrate an array of transducers, preferably for an array of ultrasound transducers, but may be used to calibrate any suitable kind of transducers in any suitable arrangement. The transducer array is preferably used to propagate signals in a presumably homogeneous environment (with constant but unknown speed of signal propagation), but may alternatively be used in any suitable environment. In an array of N transducers with emit and receive capabilities, the position of transducer i may be denoted by x i  for (i=0, 1, . . . , N−1) and each transducer may be assumed to have unknown but fixed emission and reception delays (denoted by e i  and r i , respectively, with delay vectors e and r whose ith component is respectively given by e i  and r i ). These parameters affect overall time delay measurement m i,j , or the time elapsed between when a probing signal sent to a transmitter and when a signal is recorded at a receiver. This delay can be decomposed as
 
 m   i,j   =t   i,j   +e   i   +r   j   (1)
 
where t i,j  denotes the actual propagation time between transducers i and j, also referred to as “time of flight”. The signals are assumed to propagate in an homogeneous environment with unknown propagation speed c. Hence, the time of flight t i,j  satisfies
 
 d   i,j   =ct   i,j   (2)
 
where d i,j =∥x i −x j ∥ is the Euclidean distance between transducers i and j. A time of flight matrix T may be defined with (T) i,j =t i,j , a distance matrix D may be defined with (D) i,j =d i,j , and a time delay measurement matrix M may be defined with (M) i,j =m i,j . In this manner, the distance matrix only depends on the position of the sensors x i  and that D=cT. Estimated or noisy versions of these quantities are denoted with a hat symbol (e.g., {circumflex over (M)}). In one preferred embodiment in the context of these definitions, the method  100  aims to find one or more of the calibration parameters ĉ, ê, {circumflex over (r)} and {circumflex over (x)} i  that best fit the measurements in a mean square sense. More precisely, the method may solve the optimization problem of:
 
                         min     c   ,   e   ,   r   ,     x   i         ⁢            M   -     M   ^            2       =       min     c   ,   e   ,   r   ,     x   i         ⁢                c     -   1       ⁢   D     +     e   ⁢           ⁢     1   T       +     1   ⁢     r   T       -     M   ^            2         ,           (   3   )               
where 1 denotes the all-one vector and the dependency of the matrix D on the positions x i  is implicit. The equality in (3) follows directly from (1) and (2).
 
     Step S 110 , which includes receiving an estimate of the positions of the transducers, functions to establish an initial set of positions of the transducers for the iterative algorithm. The manner in which the positions of the transducers quantitatively represent the transducer arrangement may depend on the particular nature of the arrangement. In one variation, the positions of transducers in a transducer array may be quantified relative to one or more known reference points. In another variation, the positions of transducers in a transducer array may be quantified by the position between adjacent or consecutive transducer elements. For example, for a ring transducer having transducer elements arranged circumferentially around a ring, the positions of the transducer elements may be quantified by the distance of each transducer element to the center of the ring, and/or quantified by the angle between adjacent transducer elements. In another example, for a linear array of transducer elements, the positions of the transducer elements may be quantified by the distance of each transducer from one end of the linear array. However, the estimated positions of the transducers may be quantified in any suitable manner. 
     In an array of N transducers with emit and receive capabilities, the position of transducer i may be denoted by {circumflex over (x)} i  for (i=0, 1, . . . , N−1). Receiving an estimate of the positions of the transducers may include estimating the pairwise Euclidean distances S 112  between estimated positions of transducers (denoted by d i,j =∥x i −x j ∥ as the Euclidean distance between transducers i and j), and creating a distance matrix {circumflex over (D)} S 114  of estimated pairwise distances between estimated positions of transducers (denoted by (D) i,j =d i,j ). However, the estimated positions of transducers may be represented in any suitable distance data set or other manner. The method may further include storing the estimated positions of the transducers as a current set of estimated transducer positions {circumflex over (x)} i . The received estimated positions of transducers may be previously stored, received as input from a user, transferred from a storage unit, and/or be received in any suitable manner. Alternatively, the method may further include the step of estimating the positions of the transducers, such as by analyzing a photograph or other representation of the transducers or deriving from one or more previous estimates of the positions of the transducers (e.g., by averaging multiple sets of prior estimates). 
     Step S 120 , which includes forming a time delay measurement data set, functions to contain measured and/or estimated values of time delay measurements m i,j . The time delay measurement data set preferably includes non-faulty time delay measurements and faulty time delay measurements (e.g., due to high noise, faulty or missing transducers, etc.). As shown in  FIG. 2 , forming a time delay measurement data set S 120  preferably includes forming a time delay measurement matrix M of the time delay measurements S 122 , with (M) i,j =m i,j . Step S 120  preferably further includes flagging faulty time delay measurements S 124 , which may be suspicious due to their value relative to adjacent time delay measurements, known faults in particular transducers, and/or any suitable reason. Any other non-faulty time delay measurements may additionally and/or alternatively be flagged. In one variation, flagging faulty time delay measurements S 124  includes one or more of the following steps: removing selected rows of the time delay measurement matrix {circumflex over (M)} (such as those corresponding to faulty transducers) S 126 , removing selected columns of the time delay measurement matrix {circumflex over (M)} S 127 , and applying a mask to the time delay measurement matrix {circumflex over (M)} S 128  that sets the missing time delay measurements to zero. However, flagging faulty time delay measurements S 124  may additionally and/or alternatively include setting selected rows or columns to an alternative suitable value, or any suitable flagging step. 
     Step S 130 , which includes generating from the time delay measurement data set an estimate of at least one calibration parameter, functions to find values for at least one calibration parameter that best fit the time delay measurement data set or matrix. Generating an estimate of at least one calibration parameter S 130  preferably includes optimizing values for at least one calibration parameter using the non-faulty time delay measurements. The estimates of the optimized calibration parameters are preferably obtained using non-faulty time delay measurements, but may incorporate the faulty time delay measurements in some manner (e.g., weighing the importance of the non-faulty measurements greater than the faulty measurements). In a preferred embodiment, optimizing values for at least one calibration parameter includes computing mean square estimates of at least one calibration parameter. For example, optimizing values for at least one calibration parameter involves computing the mean square optimal estimates of the speed of signal propagation ĉ, as well as the emission and reception delays denoted ê and {circumflex over (r)}, using the non-zero elements of {circumflex over (M)}. Step S 130  preferably solves the optimization problem of (3) using the non-zero elements of {circumflex over (M)}. 
     The optimization of (3) may also be written as: 
                                      c     -   1       ⁢   D     +     e   ⁢           ⁢     1   T       +     1   ⁢     r   T       -     M   ^            2     =       ⁢                vecc     -   1       ⁢   D     +     vece   ⁢           ⁢     1   T       +     vec   ⁢           ⁢   1   ⁢     r   T       -     vec   ⁢           ⁢     M   ^              2                   =       ⁢                c     -   1       ⁢   d     +       K   e     ⁢   e     +       K   r     ⁢   r     -     m   ^            2       ,                 (   4   )               
where d=vecD and {circumflex over (m)}=vec{circumflex over (M)} are vectors obtained by stacking the columns of D and {circumflex over (M)}, respectively. The matrices K e  and K r  are given by
 
               K   e     =       1   ⊗   I     =     [         I           I           ⋮           I         ]                 and               K   r     =       I   ⊗   1     =     [         1       0       …       0           0       1       …       0           ⋮       ⋮       ⋱       ⋮           0       0       …       1         ]             
where   denotes the Kronecker product and I the identity matrix. Defining the matrix K and the vector p as
 
             K   =     [         I         K   e           K   r           ]               and               p   =     [           d   /   c             e           r         ]       ,         
respectively, the optimization problem (3) can be expressed as
 
                     min     c   ,   e   ,   r   ,     x   i         ⁢              Kp   -     m   ^            2     .             (   5   )               
The above optimization depends on the positions of the transducers only through the distance vector d. The vector d is invariant with respect to translation, rotation and symmetry of the transducers; therefore the positions derived from d can be relative, but absolute positioning can easily be obtained by specifying a few references (e.g., the center of the array, the position of a reference transducer). The minimization (5) can be thus be carried in terms of this vector instead of the individual positions, that is,
 
                     min     c   ,   e   ,   r   ,   d       ⁢              Kp   -     m   ^            2     .             (   6   )               
Optimal values ĉ and {circumflex over (d)} can only be determined up to a positive scaling factor, as apparent from the definition of the vector p. The solution to the minimization (6) is obtained by orthogonal projection as
 
 {circumflex over (p)}=K   \   {circumflex over (m)},   (7)
 
where K \  is the Moore-Penrose pseudoinverse given by K \ =K T (KK T ) (−1)  in the full rank. Given an estimate of the propagation speed ĉ, the estimated delay vectors ê and {circumflex over (r)} can be retrieved from {circumflex over (p)}.
 
     Step S 140 , which includes modifying the time delay measurement data set, functions to modify faulty time delay measurements using the estimated calibration parameters. As shown in  FIG. 3 , modifying the time delay measurement data set preferably includes interpolating updated values for the faulty time delay measurements S 142  using estimated calibration parameters and modifying the faulty time delay measurements to the interpolated values S 144 . In a preferred embodiment, interpolating updated values includes interpolating updated values based on the current estimated positions of the transducers (e.g., relative positions in the distance matrix {circumflex over (D)}) and an estimated speed of propagation obtained in step S 130 . Any suitable method of interpolation may be used; several methods of interpolation are familiar to one ordinarily skilled in the art. 
     Step S 150 , which includes modifying the estimated positions of the transducers based on the estimated calibration parameters, functions to update the estimated transducer positions as part of a successive iteration as necessary in the method  100 . Modifying the estimated positions of the transducers S 150  may include modifying the estimated relative positions based on estimated values of at least emission and reception delay parameters e i  and r i . In particular, as shown in  FIG. 4 , modifying the estimated positions of the transducers S 150  preferably includes modifying the distance matrix {circumflex over (D)} S 152 . In one variation, modifying the estimated positions of the transducers S 150  further includes modifying the estimated positions of the transducers {circumflex over (x)} i  based on the modified distance matrix S 154 . In another variation, the method may further include denoising the modified distance matrix S 160  and modifying the estimated positions of the transducers based on the denoised modified distance set S 210  (e.g., modifying the estimated positions of the transducers based on the modified distance matrix S 154  without denoising the modified distance matrix). In yet another variation, the method may include modifying the estimated positions of the transducers based on the modified distance matrix S 154 , denoising the modified distance matrix S 160 , and modifying the estimated positions of the transducers based on the denoised modified distance data set S 210 , such as to compare the effect of denoising the modified distance data set). 
     Modifying the distance matrix {circumflex over (D)} S 152  may include subtracting estimated emission and reception delay values e i  and r i  from the time delay measurement matrix and multiplying the result of the subtraction by the estimated speed of propagation ĉ. However, the distance matrix {circumflex over (D)} may additionally and/or alternatively include any suitable steps for modification. In particular, based on the initial estimate d 0  of the distance vector, the modifying step S 152  may solve the optimization problem 
                       min     c   ,   e   ,   r       ⁢            Kp   -     m   ^            2       ,           (   8   )               
where the matrix K and the vector p are now defined as
 
             K   =     [           d   0           K   e           K   r           ]               and               p   =     [           1   /   c             e           r         ]       ,         
respectively. The estimated sound speed ĉ and delay vectors ê and ŝ can be retrieved from {circumflex over (p)} computed using (7), and the distance vector {circumflex over (d)} can further be estimated as
 
 {circumflex over (d)}=ĉ ( {circumflex over (m)}−K   e   ê−K   r   {circumflex over (r)} ).  (9)
 
     Denoising the modified distance matrix S 160  may include storing a testing matrix that represents elements in the modified distance matrix and iteratively reinforcing at least one desired property of the testing matrix. In the method  100 , denoising the modified distance matrix S 160  cleans the distance matrix obtained in (9) prior to modifying the estimated positions of the transducers based on the denoised modified distance matrix. The denoising step S 160  reduces the effect of the noise on calibration accuracy and hence permits the use of relatively cheap, possibly off-the-shelf electronic components in or in conjunction with the transducer array. In one variation, the testing matrix preferably includes the square of every element in the modified distance matrix {circumflex over (D)}. In another variation, the testing matrix includes the square of every element in a time-of-flight matrix defined as T=c −1  D. The reinforced properties of the testing matrix are preferably at least one of matrix symmetry, zero diagonal elements, non-negative elements, and rank equal to or less than four, although any other suitable properties may additionally and/or alternatively be reinforced. 
     As shown in  FIG. 6 , denoising the modified distance matrix {circumflex over (D)} S 160  preferably includes storing an initial testing matrix Q S 164 , manipulating the initial testing matrix Q to possess at least one desired property S 166 , determining a test value to compare to a second stopping criterion S 168  (separate from the stopping criterion in the higher-level calibrating method  100 ), and repeating the steps of storing, manipulating, and determining a test value S 174  until the second stopping criterion is met. In each iteration the repeated step of storing an initial matrix includes storing an initial testing matrix that is equivalent to the manipulated testing matrix of the previous iteration. 
     In a more detailed preferred embodiment, the step of denoising S 160  manipulates a testing matrix Q to reinforce the following key properties: 
     1. Q is symmetric, i.e., Q=Q T . 
     2. Q has zero diagonal elements, i.e., (Q) i,i =0 for i=0, 1, . . . , N−1. 
     3. Q has non-negative elements, i.e., (Q) i,j ≧0 for i,j=0, 1, . . . , N−1. 
     4. Q is of rank at most 4, i.e., rank Q≦4. 
     The property of rank of most four can easily be shown by one ordinarily skilled in the art. Due to the presence of noise in the measurements, the matrix Q obtained by squaring every elements of the distance matrix {circumflex over (D)} estimated using (9) will fail to have some of these properties. The key idea is to successively enforce these properties as a means to denoise {circumflex over (D)}. This is achieved by means of the following mappings: 
     1. φ 1 (Q)=(½)(Q+Q T ). 
     2. φ 2 (Q) sets the diagonal elements of Q to zero. 
     3. φ 3 (Q) sets the negative elements of Q to zero. 
     4. φ 4 (Q) sets the N−4 smallest singular values of Q to zero. 
     It can be shown that the sequence of matrices obtained by iteratively applying the above mappings converges, in the limit of a large number of iterations, to a matrix that possesses the four required properties (see, e.g., J. A. Cadzow, “Signal enhancement—A composite property mapping algorithm,”  IEEE Transactions on Acoustics, Speech and Signal Processing , vol. 36, no. 1, pp. 49-62, January 1988, which is incorporated in its entirety by this reference). Note that symmetry is not violated by any of the mappings. Hence, it just needs to be enforced at the beginning of the iterations. 
     As shown in  FIG. 7 , determining a test value to compare to a second stopping criterion S 168  may be one or more of several variations. In a first variation, determining a test value S 168  includes determining the difference between the initial and manipulated testing matrices Q and Q′ S 169   a . In a second variation, determining a test value S 168  includes determining the rate of change in the difference S 169   b  between the initial and manipulated testing matrices between subsequent iterations. In a third variation, determining a test value S 168  includes determining the number of iterations S 169   c  that have been performed in the current denoising step of method  100 . The repeating step S 174  then checks whether this test value satisfies a second stopping criterion (e.g., threshold in difference, threshold in rate in change of difference, predetermined number of iterations). However, determining a test value S 168  may include any suitable steps to determine any suitable test value for comparing any suitable stopping criterion. 
     This procedure, as exemplified in one preferred embodiment, is summarized below: 
     Distance Matrix Denoising 
     Input: A noisy distance matrix {circumflex over (D)}. 
     Output: A cleaned distance matrix D. 
     Procedure: 
     1. Compute {circumflex over (Q)} 1 =φ 1 ({circumflex over (Q)}), where {circumflex over (Q)} is obtained by squaring every element of {circumflex over (D)}. 
     2. Compute {circumflex over (Q)} 2 =φ 2 ({circumflex over (Q)} 1 ). 
     3. Compute {circumflex over (Q)} 3 =φ 3 ({circumflex over (Q)} 2 ). 
     4. Compute {circumflex over (Q)} 4 =φ 4 ({circumflex over (Q)} 3 ). 
     5. If the norm of the difference ∥{circumflex over (Q)} 4 −{circumflex over (Q)} 1 ∥&gt;δ for some prescribed threshold δ, set {circumflex over (Q)} 1 ={circumflex over (Q)} 4  and go to step 2. 
     6. Output D obtained by taking the square root of every element of {circumflex over (Q)} 4 . 
     In other words, the step of denoising S 160  may include: defining a second stopping criterion; storing a current testing matrix {circumflex over (Q)} representing the distance matrix S 164 ; making {circumflex over (Q)} symmetric S 182  by summing {circumflex over (Q)} with its transpose and dividing the summed result by two; setting the diagonal elements of {circumflex over (Q)} to zero S 184 ; setting all the negative elements of {circumflex over (Q)} to zero S 186 ; computing a singular value decomposition of {circumflex over (Q)} and setting all the singular values of {circumflex over (Q)} to zero except the four largest ones S 188 ; determining a test value to compare to a second stopping criterion S 168 ; comparing the test value to the second stopping criterion S 172 ; repeating the storing and various modifying steps S 174  until the norm of the difference satisfies the second stopping criterion; deriving a denoised distance matrix from the final modified {circumflex over (Q)} matrix after the second stopping criterion is satisfied S 176 ; and providing the denoised distance matrix as an output S 178 . The modifying steps may be performed in any suitable order. Or instance, the steps of computing a singular value decomposition of {circumflex over (Q)} and setting all the singular values of {circumflex over (Q)} to zero except the four largest values may be performed before the steps of setting the diagonal and negative elements of {circumflex over (Q)} to zero. 
     Step S 154  and S 190 , which includes modifying the estimated positions of the transducers {circumflex over (x)} i  based on the modified distance matrix or denoised modified distance matrix, respectively, function to provide a modified estimate of the positions of the transducers, which may be useful for providing a comparative point to determine when to stop the iteration and/or for providing a calibration parameter. Modifying the estimated positions of the transducers preferably include using a multi-dimensional scaling (MDS) algorithm, but may include any optimization or other step for obtaining updated estimated positions of the transducers from the modified distance matrix or other data set. In a preferred embodiment, the transducers positions estimates {circumflex over (X)}=[{circumflex over (x)} 0 , {circumflex over (x)} 1 , . . . , {circumflex over (x)} N-1 ] T  are computed using a MDS algorithm described below: (e.g., the MDS algorithm may be similar to that described in P. Drineas et al., “Distance matrix reconstruction from incomplete distance information for sensor network localization,” 3 rd Annual IEEE Communications Society on Sensor and Ad Hoc Communications and Networks , September 2006, pp. 536-544, which is incorporated in its entirety by this reference): 
     MDS Localization 
     Input: A distance matrix D. 
     Output: Modified position estimates {circumflex over (X)}. 
     Procedure: 
     1. Compute τ(Q)=−(½)LQL, where L=I−(1/N)11 T  and Q is obtained by squaring every element of D. 
     2. Compute τ 2 (Q), the best rank two approximation to τ(Q) using its singular value decomposition, τ 2 (Q)=U 2 Λ 2   2 U 2   T . 
     3. Compute {circumflex over (x)}=U 2 Λ 2 . 
     Step S 210 , which includes determining a test value to compare to a stopping criterion, functions to provide a measure for evaluating whether the estimated calibration parameters have a sufficient level of apparent accuracy. As shown in  FIG. 5 , step S 210  may be one or more of several variations. In a first variation, step S 210  includes determining the difference between the estimated and modified estimated positions of transducers S 212   a . Determining the difference between the estimated and modified estimated positions of transducers S 212   a  preferably includes computing the norm of the difference between a vector of estimated positions of transducers and a vector of modified estimated positions of transducers. However, any other suitable indicator may be used to represent the difference between the estimated and modified estimated positions of transducers, such the mean of the difference between respective elements in the two vectors, or any suitable indicator. In a second variation, step S 210  includes determining the rate of change in the difference S 212   b  between the estimated and modified estimated positions of transducers, where the difference is calculated similarly to that in the first variation of step S 210 . In a third variation, step S 210  includes determining the number of iterations of method  100  performed. However, step S 210  may include any suitable steps for computing a test value for any suitable stopping criterion. 
     Step S 220 , which includes repeating at least steps S 130  through S 210  until the test value satisfies a stopping criterion, functions to iterate the modifications of the time delay measurement matrix and the distance matrix until particular estimated calibration parameters have a sufficient level of apparent accuracy. As shown in  FIG. 1 , step S 220  includes comparing the test value to the stopping criterion S 222 . The stopping criterion may be a fixed numerical value or a dynamic numerical value (such as one that changes depending on the rate of change in the difference between multiple iterations). In a first case, if the test value is greater than the stopping criterion, at least a portion of the method, more preferably steps S 120  through S 210 , is repeated using the modified estimated positions of transducers and modified time delay measurement data set as the initial values for the next iteration. In a second case, if the test value is less than the stopping criterion, then the iterations preferably stop and the modified estimated positions of transducers and modified time delay parameters may be considered sufficiently accurate. In a third case, if the test value is equal to the stopping criterion, the method may perform the action of the first or the second case. 
     The method may further include step S 230 , which includes providing an output of at least one of the estimated calibration parameters. Providing an output S 190  may include displaying on a screen or printing one or more of the estimates. 
     The method may further include step S 240 , which includes saving at least one of the estimated calibration parameters. The estimates may be saved to a hard drive, on a removable storage device, a server, or any suitable storage system. For instance, one or more estimates obtained through the method may be stored for future use during an experimental run with the transducer array. As another example, the estimates may be stored for comparison to previous or future estimates obtained through the method or other means, such as for tracking changes in the estimated transducer positions or time delay parameters over time or comparing parameters corresponding to other propagation mediums. 
     In one preferred embodiment, the method for estimating speed of propagation ĉ of an homogeneous medium and positions {circumflex over (x)} i , emission delays ê i  and reception delays {circumflex over (r)} i  of transducers based on time delay measurements {circumflex over (m)} i,j  between pairs of said transducers (i,j=0, 1, . . . , N−1) includes the following steps: 
     1. building a measurement matrix comprising time delay measurements {circumflex over (m)} i,j  between any pair of transducers, 
     2. removing rows and columns of the measurement matrix that correspond to faulty transducers, 
     3. defining a stopping criterion, 
     4. building a mask that sets missing entries in the measurement matrix to zero, 
     5. choosing initial positions {circumflex over (x)} i  of the transducers, 
     6. storing the current positions, 
     7. building a matrix that contains the pairwise distances between the transducers, 
     8. applying the mask to the measurement matrix, 
     9. computing mean square optimal estimates of the propagation speed ĉ, the emission delays ê i  and the reception delays {circumflex over (r)} i  using the non-zero measurements {circumflex over (m)} i,j , 
     10. computing interpolated values for the zero entries of the measurement matrix using the current speed of propagation estimate ĉ and the distance matrix, 
     11. computing a new distance matrix by subtracting the estimated delays ê i  and {circumflex over (r)} i  to the measurement matrix and multiplying by the estimated speed of sound ĉ, 
     12. computing a denoised version of the distance matrix using an iterative algorithm, 
     13. updating the positions {circumflex over (x)} i  of the transducers using a multi dimensional scaling algorithm on the denoised distance matrix, 
     14. computing the norm of the difference between the current positions and those stored in step 6, 
     15. repeating the method from step 6 if the norm of the difference is larger that the stopping criterion, and 
     16. choosing the current estimates of at least one of the calibration parameters of ĉ, {circumflex over (x)} i , ê i  and {circumflex over (r)} i  as the output. 
     In a preferred embodiment, computing a denoised version of the distance matrix may include the steps of: 
     1. defining a second stopping criterion, 
     2. storing the current distance matrix {circumflex over (D)}, 
     3. making the distance matrix {circumflex over (D)} symmetric by adding its transpose and dividing the result by two, 
     4. computing a singular value decomposition of the distance matrix {circumflex over (D)}, 
     5. setting all the singular values of the distance matrix {circumflex over (D)} to zero except the four largest ones, 
     6. setting the diagonal elements of the distance matrix {circumflex over (D)} to zero, 
     7. setting all the negative elements of the distance matrix {circumflex over (D)} to zero, 
     8. computing the norm of the difference between the current distance matrix {circumflex over (D)} and that stored in step 2 of computing a denoised version of the distance matrix, 
     9. repeating the method of computing a denoised version of the distance matrix from step 2 if the norm of the difference is larger that the second stopping criterion, and 
     10. choosing the current distance matrix {circumflex over (D)} as the output. 
     In a preferred embodiment, a system  300  for calibration of a transducer array characterized by calibration parameters, based on time delay measurements between transducer pairs, includes: a plurality of transducers  310  that perform at least one of emitting and receiving a signal; a data controller  320  that forms a distance data set including pairwise distances between estimated positions of transducers and a time delay measurement data set including time delay measurements between transducer pairs, including non-faulty time delay measurements and any faulty time delay measurements; and estimation engine  330  that generates from the time delay measurement data set an estimate of at least one calibration parameter  312 . The data controller  320  modifies the faulty time delay measurements in the time delay measurement data set and the estimated positions of transducers in the distance in the distance data set based on the estimated calibration parameters. The data controller  320  and estimation engine  330  preferably iteratively cooperate until a test value satisfies a stopping criterion. The test value may represent the difference between the estimated positions of transducers and final modified estimated positions of transducers, the rate of change of the difference between multiple iterations, the number of iterations performed up to that point, and/or any suitable test value. The system  300  preferably generates estimates of at least one of calibration parameters including: positions of transducers, emission delay, reception delay, and speed of signal propagation, but may additionally and/or alternatively generate estimates of any other suitable calibration parameter. The system may further include a signal processor  340  that denoises the modified distance data set with an iterative algorithm, preferably performing steps similar to step S 160  in method  100 , but alternatively any suitable denoising algorithm. In a preferred embodiment, the data controller and/or the estimation engine outputs (e.g., displaying or printing) or stores at least one of the estimated calibration parameters  312 . 
     EXAMPLE 
     The calibration method may be applied to the calibration of a set of N=256 ultrasound transducers used in a clinical prototype for breast cancer detection. The transducers are initially assumed to be uniformly spaced with an angular separation a of 0.024 radians&lt;a&lt;0.025 radians on a ring of radius r=0.0999 meters. Furthermore, the overall initial emission/reception delay of each transducer pair is assumed to be a constant τ. This amounts to set the matrices K e  and K r  in (4) as the all-one vector and choose e and r as the scalar τ/2. Measurements are computed from a water shot of a patient using a time delay estimation method based on the arrival of the first peak of the transmitted signal. The obtained measurement matrix is depicted in  FIG. 9 . 
     Due to the directional response of the transducers, the signals transmitted between adjacent pairs of emitters and receivers have a low signal-to-noise ratio and the corresponding time delays cannot be computed. In this case, as shown in  FIG. 9 , the 20 first off-diagonal elements of the measurement matrix are missing. Moreover, the transducers  30 ,  96 ,  144  to  151 , and  243  to  246  are considered as too noisy and are discarded. 
     The algorithm was run on the time delay data including the measurement matrix. As shown in  FIGS. 10A and 10B , the estimated positions differ from the assumed uniform ring geometry. In  FIG. 10A , the distance of every element to the center of the array can be plotted, and the values can be clearly seen to depart from the assumed constant radius r=0.0999 meters. A block structure can also be observed that reflects the way the array was built, that is, by concatenating small line arrays. In  FIG. 10B , the angle differences between two consecutive transducers can be plotted and compared to the uniform angular spacing 0.024 radians&lt;a&lt;0.025 radians of the assumed ring geometry. Here, significant differences can be observed every 32 elements, which probably correspond to consecutive transducers lying on two different line arrays. Note that the positions of the bad transducers have been linearly interpolated (in radius and angle) between the closest good transducers. For this water shot, the estimated sound speed is ĉ=1509.5 meters per second, and the delay of the transducers, including the systematic delay introduced by the time delay estimation method (i.e., the time to the first peak), is about 2.68 microseconds. The temperature of the water during the experiment is T=30.3740 degrees Celsius, which translates into a sound speed of c=1510.1 meters per second according to the formula c=1405.3+4.624 T−0.0383 T 2 . Note that this formula for c is accurate for temperatures between 10 and 40 degrees Celsius. The water sound speed is thus estimated with high accuracy. 
     To check consistency of the calibration parameters, the calibration was performed again using a second water shot of the same prototype with another patient. In this case, the estimated sound speed is ĉ=1500.5 meters per second, and the delay is about 2.63 microseconds. The temperature of the water during experiment is T=26.6 degrees Celsius which corresponds to a sound speed of c=1500.9 meters per second. As shown in  FIGS. 11A and 11B  (distance from the center and angle between adjacent transducer elements, respectively), the estimates of the positions of the transducers using the first and second water shots show a good consistency between the two different calibrations, which suggest that the calibration parameters are well estimated and consistent between water shots. 
     As a person skilled in the art will recognize from the previous detailed description and from the figures and claims, modifications and changes can be made to the preferred embodiments of the invention without departing from the scope of this invention defined in the following claims.