Sum-of-harmonics method for improved narrowband and broadband signal quantification during passive monitoring of ultrasound therapies

Passive Acoustic Mapping (PAM) enables real-time monitoring of ultrasound 2 therapies by beamforming acoustic emissions emanating from the ultrasound focus. Reconstruction of the narrowband or broadband acoustic emissions component 4 enables mapping of different physical phenomena, with narrowband emissions arising 5 from non-linear propagation and scattering, non-inertial cavitation or tissue boiling, 6 and broadband (generally, of significantly lower amplitude) indicating inertial 7 cavitation. Currently, accurate classification of the received signals based on pre- 8 defined frequency-domain comb filters cannot be guaranteed because varying levels 9 of leakage occur as a function of signal amplitude and the choice of windowing 10 function. This work presents a time-domain parametric model aimed at enabling 11 accurate estimation of the amplitude of time-varying narrowband components in the 12 presence of broadband signals. Conversely, the method makes it possible to recover a 13 weak broadband signal in the presence of a dominant harmonic or other narrowband 14 component. Compared to conventional comb filtering, the proposed sum-of- 15 harmonics method enables PAM of cavitation sources that better reflect their physical 16 location and extent. 17


I. INTRODUCTION
in a setup may have overlapping bandwidths, allowing the simultaneous recording of 70 broadband, as well as narrowband emissions that are harmonics of the driving frequency of the therapeutic pulse. The relative magnitude of narrowband to 72 broadband emissions varies significantly depending on conditions such as the 73 amplitude of the transmitted signal, the degree of tissue nonlinearity, the distance 74 travelled by the transmitted acoustic wave until it is recorded, and the bandwidth of 75 the receiver relative to the frequency of the transmitter. Therefore, an adaptive 76 filtering method is required to accurately distinguish the frequency band of interest 77 regardless of the relative magnitude of the simultaneously recorded interfering signal. 78 The most commonly experimentally used method is filtering in the frequency 79 domain using a pre-defined comb filter (Choi and Coussios, 2012). However, the 80 accuracy of this technique cannot be guaranteed due to the fixed characteristics (i.e.

170
As already mentioned, the channel data, i.e. the signals recorded passively by 171 one or multiple sensors on an ultrasound array, are filtered as an initial signal 172 processing step in order to distinguish the broadband and narrowband components.

173
Assuming that the signal recorded by one sensor is y, it could be expressed as the sum 174 of a broadband signal b and a harmonic (narrowband) signal h at each time sample n: where ω 0 =2πf 0 T is the driving frequency of the system and therefore the fundamental these matrices has dimension N x q, with (n, k) elements cos(kω 0 n) and sin(kω 0 n), 201 respectively. The notation emphasizes that these matrices depend on the fundamental 202 frequency ω 0 and the number of harmonics q.

203
The harmonic signal parameters that must be estimated from measured data are The negative log-likelihood function ( 0 , , 2 ) = − ( ( | 0 , , 2 )) 215 has the following form after dropping terms that are independent of the parameters, 216 ( 0 , , 2 ) = ( 2 ) ( 2 ) + 1 2 2 ( 0 , ) The maximum likelihood (ML) parameter estimates are the values that  The MLE of the fundamental frequency for a fixed number of harmonics q is then the 232 value of ω0 that minimizes the nonlinear least-squares cost function ⊥ ( 0 ) .

233
However, since q is not known, the number of harmonics and fundamental frequency    2q(m+1)+1 from (13), and using the negative log-likelihood function in (7), the MDL 284 criterion is (dropping terms that are independent of q and 0) As soon as ω 0 and q are defined, the estimate of x can be obtained by ̂= increases. Therefore, it is proposed that the signal is segmented into several cycle 295 parts and that the proposed method is applied independently on each part.

298
An initial computational investigation was carried out in order to enable 299 understanding of the relationship between the parameters of the model including the 300 sampling frequency of the signal, the narrowband (harmonic)-to-broadband signal 301 ratio, the number of cycles per segment, and the order of polynomial used to 302 approximate the signal amplitude changes.

303
The input signal to these simulations comprises a fixed sum-of-harmonics 304 signal and white Gaussian noise with variance 2 that is chosen to set the harmonic to  The range of values for each of the parameters investigated in these 333 simulations is summarized in Table 1:   By selecting the order of the three comb filters to be equal to 11, the only 492 parameter allowed to vary amongst them is the bandwidth. In Fig.4  Since no inertial cavitation takes place in the pipette, the generated maps 578 should display no activity. In Fig.7 the maps corresponding to all filtering techniques 579 are displayed, using the same linear color scale, at three different pressures (1MPa,   place. Due to its wide bandwidth, it removes a significant part, but not the whole, of 682 the harmonic energy, as was already demonstrated, and simultaneously a part of the 683 broadband energy (Fig.6). Therefore, it could lead to underestimation of the intensity 684 and extent of inertial cavitation and potentially cause the onset of inertial cavitation 685 activity to be missed, as was observed in the PAM maps in Fig.8