(PDF) Continuous Wavelet Transform-Based Frequency Dispersion Compensation Method for Electromagnetic Time-Reversal Imaging

The University of Manchester Research Continuous Wavelet Transform Based Frequency Dispersion Compensation Method for Electromagnetic Time-Reversal Imaging DOI: 10.1109/TAP.2016.2647594 Document Version Accepted author manuscript Link to publication record in Manchester Research Explorer Citation for published version (APA): Abduljabbar, A., Yavuz, M. E., Costen, F., Himeno, R., & Yokota, H. (2017). Continuous Wavelet Transform Based Frequency Dispersion Compensation Method for Electromagnetic Time-Reversal Imaging. IEEE Transactions on Antennas and Propagation, 1321-1329. https://doi.org/10.1109/TAP.2016.2647594 Published in: IEEE Transactions on Antennas and Propagation Citing this paper Please note that where the full-text provided on Manchester Research Explorer is the Author Accepted Manuscript or Proof version this may differ from the final Published version. If citing, it is advised that you check and use the publisher's definitive version. General rights Copyright and moral rights for the publications made accessible in the Research Explorer are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. Takedown policy If you believe that this document breaches copyright please refer to the University of Manchester’s Takedown Procedures [http://man.ac.uk/04Y6Bo] or contact

providing relevant details, so we can investigate your claim. Download date:01. Jun. 2020 1 Continuous Wavelet Transform Based Frequency Dispersion Compensation Method for Electromagnetic Time-Reversal Imaging Ammar M. Abduljabbar, Student Member, IEEE, Mehmet E. Yavuz, Fumie Costen, Senior Member, IEEE, Ryutaro Himeno, Hideo Yokota Abstract—The invariance of wave equations in lossless media of the wave equations is broken. While the additional phase allows the Time Reversal (TR) technique to spatio-temporally shift undergone due to dispersion can be compensated by the refocus back-propagated signals in a given ultrawideband imag- TR process itself [12], attenuation which affects the signals ing scenario. However, the existence of dispersion and loss in the propagation medium breaks this invariance and the resultant TR in both the forward and backward propagation stages are focusing exhibits frequency and propagation duration dependent not compensated. This attenuation ultimately degrades the degradation. We propose an algorithm based on the continuous focusing resolution of the localization. The earliest works on wavelet transform that tackles this degradation to improve focus- compensation of such attenuation involved a uniform absorp- ing resolution under such conditions. The developed algorithm tion model over the entire bandwidth for acoustic propagation has been successfully applied to the scenario for localization of lung cancer. in skull [14], multipaths in ocean communication [15] and at a single frequency with the so-called decomposition of Index Terms—Time reversal (TR), continuous wavelet trans- the TR operator method [16]. A more recent compensation form (CWT), dispersive media, attenuation compensation. technique using the Short-Time Fourier Transform (STFT) method takes frequency dependency into account and tries I. I NTRODUCTION to improve the resolution focusing of the TR technique in L OCALIZATION of targets behind structures or within obstructed media is of great interest in many microwave imaging applications such as detection of landmines [1], dispersive and lossy media [12]. Although this method proved to be a good candidate for compensation, the inverse filters utilized in the method also amplify unwanted oscillations explosives [2], subsurface imaging [3] or tumors in the body and perturbations in the received signals, which affect the [4] (e.g. breast [5], brain [6] and lung [7].) performance of TR technique. STFT-based method of [12] A recent method in the microwave imaging which has produced the inverse filters by comparing the solution of been introduced to detect and localize objects is the Time the wave equation in the applied dispersive medium against Reversal (TR) technique [8]. In a lossless medium, TR states a corresponding non-dispersive test medium by transmitting that for every wave component propagating away from a a pulse and observing it at various propagation distances source point following a certain path, there is a corresponding for both media. The process of creating the inverse filters time-reversed wave that can trace the same path back to the is empiric and time consuming. Later, the method in [12] original point of the source due to invariance of Maxwell’s is extended by adding a threshold approach to overcome equations to time. Conventional TR technique uses a set of the amplification of the noise in the inverse filter in [17]. antennas, called TR Array (TRA), to receive, record and time- Different window types and lengths are applied in [17] to reverse the reflected signal from the object(s). Time-reversed obtain the optimum settings for the STFT technique. However, signals at each TRA antenna are sent back to the medium the work only applied in dispersive and homogeneous media to focus around the original object location. One important which do not accurately represent the real-life applications. property of TR is that it achieves super-resolution by utilizing The inverse filters in [12] and [17] need the non-dispersive the multipath propagation in the medium as discussed in version of the propagation media to create the inverse filter [9], [10], [11]. On the other hand, as with other microwave which might be unknown in real-life scenarios. The STFT imaging techniques, TR suffers from dispersion and losses method uses a fixed time window length which provides a in the propagating medium [12], [13] where TR invariance lower resolution in determining attenuation when compared with other size-adjustable window methods [18], [19]. In A. M. Abduljabbar, and F. Costen are with the School of Electrical and Electronic Engineering, The University of Manchester, U.K. (email: this paper, we propose a compensation method that uses

) adaptive-window scheme based on the Continuous Wavelet M. E. Yavuz, Hillsboro, Oregon, USA. (email:

) Transform (CWT). Our compensation method uses an inverse R. Himeno is Advanced Center for Computing and Communication, RIKEN, Saitama, Japan filter in the wavelet domain to overcome the attenuation and H. Yokota and F. Costen are with the Image Processing Research Team, dispersion in the electromagnetic (EM) wave propagation due Center for Advanced Photonics, RIKEN, Saitama, Japan. to the dispersive medium. Our method uses adjustable-length Color versions of the figures in this paper are available online at http://ieeexplore.ieee.org. Additional research data supporting this publication windows by applying long time-windows at low frequencies are available from http://dx.doi.org/ repository at 10.15127/1.306001. and short time-windows at high frequencies [20]. Therefore 2 l 1 N ∆t m l m we do not need optimization of the time-window length as scale step, J is the largest scale J = log2 ( ) , was carried out in [17]. The proposed inverse filters need ∆j a0 represents the ceiling function, N is the sampling points in only the complex permittivity of one dominant medium at the time domain signal, k is frequency index and the angular the center frequency to create the inverse filter. creation of frequency ωk is [23] inverse filters in the STFT method requires more computation  compared to those in the CWT, hence another advantage of  0 k=1 2πk CWT over STFT. To the best knowledge of the authors, such ωk = 1 < k ≤ N2 + 1  N2πk∆t N an adaptive compensation algorithm has not been considered − N ∆t 2 + 1 < k ≤ N elsewhere before. Therefore introduction of such an algorithm for k = 1, 2, . . . , N . a0 is dimensionless and it is set to is considered as the main contribution of this paper. Numerical the value of ∆t and ∆j is set to 0.025 based on the value simulations have been chosen to validate the algorithm. The chosen in [23]. Both Ψ0 [n] in (1) and Ψ0 (aj ωk ) in (2) are outline of the paper is as follows: First, the details of the complex. In order for the analyzing function Ψ0 [n] to be called proposed compensation method is described in Section II. wavelet it needs to be admissible [24]. Admissibility is one of Section III and Section IV provide the radio environment the necessary conditions of the wavelet transform and it is settings and the analysis of the numerical simulations for the achieved when [24] canonical test and practical application respectively. Section V N presents the conclusions. X lim Ψ0 [n] = 0 (3) N →∞ n=−N II. C OMPENSATION OF DISPERSIVE ATTENUATION IN HUMAN TISSUES and (aj ωk −fΨ )2 ( A dispersive and lossy medium acts as a low-pass filter for √1 − 4 πe 2 ωk > 0 signals propagating through it. Thus, for compensation of such Ψ0 (aj ωk ) = 0 ωk ≤ 0. effects, inverse filters are needed. The real part of the medium dielectric permittivity causes a phase shift to the traveling We added white Gaussian noise to each observed signal x[n] waves during the forward propagation of the TR process. at each TRA antenna. The average power of the additive noise Thanks to the TR concept, this phase shift is inherently is 20% of that of x[n] which is chosen to mimic the system corrected in the backward propagation due to the fact that noise based on [25]. Our compensation method starts with the signals in the backward propagation are phase-conjugated conversion of the observed signal x[n] in time domain to the coherently along all bandwidth [21]. On the other hand, the wavelet domain. The CWT of x[n] is the convolution of the imaginary part of causes inevitable signal attenuation. The wave function Ψ∗ ( anj ) and x[n] which is equivalent to the attenuation is a function of time (in other words, duration the inverse Fourier transform of the product of Ψ∗ (aj ωk ) and signal has traveled in the dispersive medium) and frequency. x[k] in the frequency domain where ∗ represents the complex Note that the more the signal propagates in a dispersive conjugate, Ψ( anj ) is medium the more it is attenuated (duration dependency). Simi- s larly, for certain wave packet that has traveled in the dispersive n ∆t n Ψ = Ψ0 (4) medium, different frequency components undergo different aj aj aj attenuations (frequency dependency). Therefore the proposed method uses CWT technique which inherently incorporates and Ψ(aj ωk ) is both time and frequency variations as CWT suits well for r 2πaj this kind of problem. The wavelet function that we used for Ψ(aj ωk ) = Ψ0 (aj ωk ). (5) ∆t the CWT method is Morlet wavelet [22]. Morlet wavelet is considered in [19] as the proper choice for analyzing the Ψ(aj ωk ) represents the normalized Ψ0 (aj ωk ) to guarantee propagated wave in attenuating and dispersive media. The that the CWTs at each aj are comparable to each other and mother wavelet of Morlet in time domain [23] is to the CWTs of other time series [23]. Ψ(aj ωk ) satisfies 1 fΨ n∆t − (n∆t)2 N Ψ0 [n] = √ e e (1) X | Ψ(aj ωk ) |2 = N. 2 4 π (6) k=1 where n is the time index, ∆t is the time step in seconds and fΨ is the central frequency of the mother wavelet. Note The CWT for the observed signal x[n] at each TRA antenna that without loss of generality, we have chosen to work with is N N ! discrete signals instead of the continuous ones, hence the X X 2π 2π notation of Ψ0 [n] instead of Ψ0 (t). Ψ0 [n] is expressed in X[n, aj ] = x[i]e− N ik Ψ∗ (aj ωk ) e N nk . (7) scaled-frequency domain as k=1 i=1 1 (a ω −f )2 The CWT method applies a long window (large aj ) at the − j k2 Ψ Ψ0 (aj ωk ) = √ e (2) lower end of the spectrum (low frequency) and a short window 4 π (small aj ) at the higher end of the spectrum (high frequency) where aj is the dimensionless scaling factor of aj = a0 2j∆j to provide a perfect filter for time-dependent implementations for j = 0, 1, . . . , J − 1, a0 is the smallest scale, ∆j is the in terms of time and frequency localization. Furthermore 3 the CWT method provides us with an improved separation Ψ∗[aj ωk ] between the unwanted noise and our signal [19]. Next step is applying our inverse filter to X[n, aj ] in the wavelet domain to x[n] x[k] X[n, aj ] FFT IFFT compensate the attenuation caused by the dispersive medium. In order to derive the inverse filter, we need to calculate the Continous Wavelet Transform attenuation. The solution of Maxwell equations for a plane (CWT) Hs[n, aj ] wave propagating in a dispersive medium [26] is √ √1 aj E[n] = e2πf ·n∆t e−2πf µd[n] (8) Xc[n] J Y [n, aj ] Cδ−1 Σ j=1 ℜ[Y ] where f is the frequency of interest, d[n] is the distance between the excitation and the observation, µ = µ0 µr is the Inverse Continous Wavelet permeability of the medium, µ0 is the free space permeability Transform (ICWT) and µr is the relative permeability which is set to 1 as a non- magnetically dispersive medium is assumed, = 0 r is the Fig. 1. Compensation method for the dispersive attenuation in the wavelet domain where ⊗ means multiplication in frequency domain and corresponds permittivity of the medium, 0 is the free space permittivity, to convolution in time domain. S − ∞ σ r = ∞ + − is the complex relative 1 + 2πf τD 2πf 0 permittivity of medium, ∞ is the optical relative permittivity, stabilized compensated wave Y [n, aj ] in the wavelet domain S is the static relative permittivity, τD is the relaxation time, as σ is the static conductivity. Y [n, aj ] = X[n, aj ]Hs [n, aj ]. (15) (8) can be rewritten as √ d[n] To avoid the amplification of the high frequency noise, E[n] = e2πf ·n∆t e−2πf r C Hs [n, aj ] is set to 1 for 0 ≤ j ≤ J2 (short window) √ d[n] √ d[n] = e2πf ·n∆t e−2πf <[ r ] C e2πf =[ r ] C (9) which affects the higher part of the spectrum. Finally the compensated signal Xc [n] is obtained by applying the inverse , e2πf ·n∆t · Θ · Γ continuous wavelet transform (ICWT) to Y [n, aj ] as [23] where C is the speed of light C = √ 1 in the vacuum, J √ √ 0 µ0 1 X <[Y [n, aj ]] <[ ] and =[ r ] are the real and the imaginary parts of Xc [n] = √ (16) √ r Cδ j=1 aj r , Γ is the attenuation defined as √ d[n] where Cδ is a constant for each wavelet function. It is Γ[n, f ] = e2πf =[ r ] C (10) calculated as J and Θ is the phase shift defined as X <[Xδ [aj ]] Cδ = √ (17) √ d[n] aj Θ[n, f ] = e−2πf <[ r ] C . (11) j=1 for the Morlet wavelet where Xδ is the CWT of a delta Our filter is obtained from (10) and (11) as function (δ) as 1 H[n, aj ] = N Θ[n, aj ] · Γ[n, aj ] 1 X ∗ √ √ (12) Xδ [aj ] = Ψ (aj ωk ). (18) d[n] 2π fc <[ r ] C −2π a fc d[n] =[ r ] C N = e aj e j k=1 where aj is a scaling factor which controls the actual fre- The parameters for the CWT and the ICWT were selected quency f (, afcj ) of the filter in order to have a variable window based on the theory in [23]. Our proposed method is depicted in Fig. 1. size [19] and fc is the central frequency of x[n]. H[n, aj ] uses r of the dominant tissue in the propagating media at fc . H[n, aj ] is an exponential function which amplifies the noise III. C ANONICAL TEST in compensated signal and thus causes instability. Therefore We expect that the proposed method improves the localiza- H should be stabilized [27] as tion and the resolution of the conventional TR imaging because Γ[n, aj ] 1 the inverse filters should compensate the losses caused by the Hs [n, aj ] = · (13) dispersive media. We evaluate the accuracy and the resolution Γ[n, aj ]2 + T Θ[n, aj ] of our method in a random medium [28] to increase the where Hs is the stabilized filter and T is the stabilization complexity of the propagating medium. The random medium factor. The system is stable with a fixed value of T = 10−3 is of Gaussian distribution and the spatial random relative for all our simulations. Alternatively we can adaptively set T permittivity is defined as as [19] σg2 rnd [x, y, f ] = r [x, y, f ] + G[x, y] (19) T = 2 (14) σs where r is calculated using the one-pole Debye [29] pa- where σg and σs are the variances for the noise and the signal rameters of the muscle tissue [30][31], [x, y] is the Cartesian respectively. By applying Hs [n, aj ] to X[n, aj ], we obtain the coordinates in the finite difference time domain (FDTD) space, 4 150 ℜ[ǫr ] Conductivity PML Max Max Random Dispersive Distance(×∆y) Distance(×∆y) 120 Medium Distance(×∆y) Scatterer 90 PML PML (PEC) 60 Min Min 30 Distance(×∆x) Distance(×∆x) TRA aperture (a) (ls , η) = (4∆s, 0.14∞ ) (b) (ls , η) = (4∆s, 0.14∞ ) 0 PML 0 30 60 90 120 150 with minimum and maximum with minimum and maximum Distance(×∆x) values of 15.04 and 42.16 values of 0.38 [S/m] and 1.13 respectively. [S/m] respectively. Fig. 2. The simulations setups. ℜ[ǫr ] Conductivity Max Max Distance(×∆y) Distance(×∆y) G is a function of space defining the random fluctuation on rnd [x, y] and is calculated as v Rx XRy u Rx Ry 1 X uX X −2π( Rux vy +R ) G[x, y] = t L[x, y]e x y Rx Ry u=1 v=1 x=1 y=1 Min Min R R y ! Distance(×∆x) Distance(×∆x) x X ux vy ux vy −2π( R + ) 2π( R +R ) X · R·e x Ry e x y (c) (ls , η) = (10∆s, 0.14∞ ) (d) (ls , η) = (10∆s, 0.14∞ ) x=1 y=1 with minimum and maximum with minimum and maximum (20) values of 16.88 and 39.51 values of 0.44 [S/m] and 1.06 respectively. [S/m] respectively. where Rx × Ry is the 2D FDTD size, [u, v] is the Cartesian Fig. 3. The distribution of r for Fig. 2 at 3 GHz. coordinates in the spatial frequency domain, R is a Gaussian random number with zero mean and probability density func- ``` ``Resolution Rt Rs −ζ 2 √1 e (ls , η) ``` ` Rc Rc tion of 2πη , ζ is a random variable, η is the variance 2η (4∆s, 0.14∞ ) 3.13 3.76 and L is the correlation function between r at two different (10∆s, 0.14∞ ) 1.51 3.34 spatial points [x1 , y1 ] and [x2 , y2 ] given by a Gaussian function TABLE I T HE SPATIAL RESOLUTION RATIO . as |x1 −x2 |2 +|y1 −y2 |2 − l2 L[x, y] = ηe s (21) where ls is the transverse correlation length. Gaussian function B. Numerical Results is chosen for its generalization and mathematical proper- Fig. 3 shows the example of the spatial distribution of r ties [28]. The random media is characterized by ls and η. We for Fig. 2 at 3 GHz when ls is set to 4∆s and 10∆s. Fig. 4 chose muscle tissue because it is one of the most dispersive (a), (b) and (c) show |Ez | in the xy-plane with conventional tissues in the human body. TR, with the method in [12] and with our compensation method. Fig. 5 (a) and (b) show the cross-section of Fig. 4 A. Radio Environment Settings at y = 120∆y with (ls , η) = (4∆s, 0.14∞ ) and the one with (ls , η) = (10∆s, 0.14∞ ) respectively. Rc = (x4 − x3 )∆x, The simulation setup is depicted in Fig. 2. The simulation Rs = (x5 − x2 )∆x and Rt = (x6 − x1 )∆x represent uses the TMz polarized 2D FDTD method. The FDTD space the resolution of the spatial focusing with the proposed is 150 × 150 uniformly sampled with the spatial step of compensation method, the method in [12] and without any ∆s , ∆x = ∆y = 1mm. The time step ∆t is set to compensation method respectively where xi (1 ≤ i ≤ 6) 2.36ps. The CFS-PML absorbing boundary conditions [32] are on the cross-range and |Ez (xi )| = 0.5 for 1 ≤ i ≤ 6 are used with 9-cell layer to mimic open space settings. A where Ez is normalized and 0.5 is the cross-range at half TRA composed of 15 transceivers antennas are equally spaced maximum to calculate the resolution of the spatial focusing. from each other and located parallel to the x-axis. The TRA Rs has the aperture of 70∆s. The first TRA antenna is located The spatial resolution of the proposed method is = 3.76 Rc at (40∆x,25∆y) and the 15th is located at (110∆x,25∆y). A (Fig. 5(a)),3.34 (Fig.5(b)) times higher than the compensation round perfect electric conductor (PEC) scatterer with a 2mm Rt radius is centered at (75∆x,120∆y). The 8th antenna of the method in [12] and = 3.13 (Fig. 5(a)),1.51 (Fig. 5(b)) Rc TRA transmits a Gaussian pulse modulated at 3 GHz covering times higher than without using any compensation method from 1.72 GHz to 4.2 GHz. as summarized in Table I. Fig. 6 shows the cross-section 5 Normalised |Ez (x, y)| 1 140 1 Normalized |Ez | [V/m] Distance(×∆y) CWT 100 0.48 0.8 TR STFT 0.6 60 0.13 0.4 Rc 20 0 Rs 20 60 100 140 0.2 Distance(×∆x) Rt (a) Conventional TR. 0 x2 x1 x3 x4 x6x5 Normalised |Ez (x, y)| 20 40 60 80 100 120 140 1 140 Distance(×∆x) (a) (ls , η) = (4∆s, 0.14∞ ). Distance(×∆y) 0.61 100 1 Normalized |Ez | [V/m] CWT 0.29 TR 0.8 60 STFT 0.08 0.6 20 0 0.4 Rc 20 60 100 140 Distance(×∆x) Rs (b) TR with STFT. 0.2 Normalised |Ez (x, y)| Rt 1 140 0 x2 x1 x3 x4x6 x5 20 40 60 80 100 120 140 Distance(×∆y) 100 0.5 Distance(×∆x) (b) (ls , η) = (10∆s, 0.14∞ ). 60 Fig. 5. The cross-section of the normalized |Ez | distribution in Fig. 4 at 0.16 y = 120∆y in the FDTD space after applying CWT, STFT and without applying any compensation method. 20 0 20 60 100 140 Distance(×∆x) 1 (c) TR with CWT. (4∆s , 0.14ǫ∞ ) Normalized |Ez | [V/m] Fig. 4. The spatial distribution of |Ez | in the xy-plane for the canonical (10∆s , 0.14ǫ∞ ) test with conventional TR (without applying any compensation method), 0.8 (4∆s , 0.05ǫ∞ ) with the method in [12] (STFT) and with our compensation method (CWT) when (ls , η) = (4∆s, 0.14∞ ). The scatterer and the TRA antennas are represented by “◦” and “x” respectively. ∆x = ∆y = 1mm. 0.6 0.4 at y = 120∆y with our compensation method for the cases Scatterer Center (ls , η) = (4∆s, 0.14∞ ), (10∆s, 0.14∞ ) and (4∆s, 0.05∞ ). 0.2 TR techniques utilize the multipaths components generated by the random media to achieve more accurate focusing resolution 0 than in homogeneous media. Both the high correlation length 20 40 60 80 100 120 140 ls such as (10∆s, 0.14∞ ) and the low variance η such as Distance(×∆x) (4∆s, 0.05∞ ) reduce the randomness of the medium as both cases result in reduction of the overall fluctuations of the per- Fig. 6. The cross-section of the normalized |Ez | distribution at y = mittivity of the propagation medium, hence its scattererness. 120∆y in the FDTD space after applying CWT for the cases (ls , η) = (4∆s, 0.14∞ ), (10∆s, 0.14∞ ) and (4∆s, 0.05∞ ). Reducing the randomness decreases the multipath components which lowers the resolution of the TR focusing. 6 100 50 ℜ[ǫr ] lung tumor 10 ℜ[ǫr ] healthy lung Distance(×∆y) 150 1 250 conductivity [S/m] lung tumor conductivity [S/m] healthy lung 350 0.1 1 2 3 4 5 (a) Frequency [GHz] tumor 450 Fig. 8. Relative permittivity and the conductivity of the lung and the scatterer. 50 150 250 Distance(×∆x) (b) We shifted the frequency response of the healthy lung tissue Fig. 7. The human phantom with 5mm radius round tumor. The TRA antennas to match the complex permittivity of the lung cancer tissue are represented by “|”. ∆x = ∆y = ∆z = 1mm. presented in [42] and generated the Debye parameters for the lung tumor by data-fitting of the modified frequency response. Fig. 8 shows the relative permittivity and the conductivity of When the computational times between our proposed the lung and the tumor. The size and the Debye parameters of method and original STFT method of [12] are compared, the tumor are set to mimic an early stage cancer. The antenna it is not merely comparison between the STFT and CWT at (31∆x,281∆y,95∆z) transmits a Gaussian pulse modulated operations but the whole algorithm implementation starting at 3 GHz. The remaining simulation settings are the same as from the filter generations. The way [12] used to produce the in Section III-A. inverse filter took the major part of the simulation time because the creation of the inverse filters in [12] needs two FDTD calculations while our approach with CWT does not require B. Numerical Results any FDTD calculations. Applying different FDTD sizes with Fig. 9 shows |Ez | on z = 95∆z plane within the 3D DHP different simulation times shows that the proposed method is with conventional TR without applying any compensation, 1.8 times faster than the STFT method in [12] and [17] in one with the method given in [12] (STFT) and one with our average. method (CWT). Fig. 9 (b) shows that the method in [12] did not improve the conventional TR imaging or locate the tumor IV. P RACTICAL APPLICATION correctly. The inverse filters utilized in [12] not only amplify Lung cancer is identified as one of the deadliest cancers af- the signal of interest, but also the noise. This reduces the fecting people [33], [34]. Early-stage lung cancer is difficult to signal to noise ratio (SNR) and in return lowers the accuracy detect because the tumor is small and its electromagnetic prop- of the localization. The proposed CWT filters compensate the erties are close to the parameters of the healthy lung. We apply losses caused by the dispersive tissues. Therefore application the TR technique with our proposed compensation method to of the proposed CWT method to the TR technique gives more locate an early-stage tumor in a digital human phantom (DHP). accurate spatial focusing than both the conventional TR and The DHP was provided by RIKEN (Saitama,Japan) under the method in [12]. non-disclosure agreement between RIKEN and the University of Manchester. The usage was approved by RIKEN ethical V. C ONCLUSION committee. The DHP has 1 mm resolution and contains 52 The wave equation invariance of the time reversal technique segmented tissue. We fitted the one-pole Debye parameters is broken in human tissues due to their dispersive charac- of human tissues [35] using the measurement provided by teristics. Since the dispersive attenuation is both duration [36], [37]. The Debye media parameters for human tissues are and frequency dependent, we have introduced CWT based presented in [31]. The purpose of this numerical simulation is compensation method that uses inverse filters in the wavelet to evaluate the accuracy and the resolution of our method in domain to overcome this attenuation. The inherent time and the human phantom medium. frequency variations in a wavelet decompose the recorded signals into different time and frequency components to which A. Radio Environment Setting different filters would be applied for compensation. CWT does Fig. 7 shows the DHP with an early stage tumor repre- not need optimization of the time-window length as is the case sented by a round scatterer with a 5mm radius centered at with STFT. The proposed inverse filters need only the complex (115∆x,285∆y,95∆z). The FDTD space is 265 × 490 × 290 permittivity at the center frequency of one dominant medium including the CFS-PML cells. The 102 TRA z−directed 14 to create the inverse filter as opposed to the whole wide-band mm dipole antennas [38] are placed in the fat tissue assuming dispersion characteristic in the STFT method [12]. Therefore that we use the implantable antennas [39], [40], [41]. The hu- CWT technique can be applied in real-life scenarios while man tissues are represented using the one-pole Debye model. STFT method [12] [17] needs the non-dispersive version of the 7 Normalised |Ez (x, y)| propagation media to create the inverse filter (in real-life, such 1 an information might not be available readily). We applied our compensation method to the TR signals in different scenarios 50 to examine how our method improved the TR imaging. In Distance(×∆y) the canonical numerical simulation, our compensation method 150 improved the spatial focusing of the TR technique. The TR 0.5 focusing around the object location is more accurate with our proposed method than those with the method in [12] or without 250 applying any method (i.e. with conventional TR imaging). Furthermore, the resolution is 3.76 and 3.13 times higher with 0.16 the proposed method than with the method in [12] or without 350 any compensation method respectively. The tumor is located accurately after applying our proposed compensation method. 450 An additional advantage of the CWT-based method is that it 0 is 1.8 times faster than the STFT method in average. 50 150 250 Distance(×∆x) VI. ACKNOWLEDGMENT (a) TR with CWT. “Wavelet software was provided by C. Tor- Normalised |Ez (x, y)| rence and G. Compo, and is available at URL: 1 http://paos.colorado.edu/research/wavelets/”. The base code 50 to produce Fig. 7 (a) was provided by one of our research group members, Mr. K. Tekbas. Distance(×∆y) 0.6 150 R EFERENCES [1] P. D. Gader, M. Mystkowski, and Y. Zhao, “Landmine detection with ground penetrating radar using hidden markov models,” IEEE Trans. 250 0.3 Geosci. Remote Sens., vol. 39, No. 6, pp. 1231–1244, June 2001. [2] P. Pratihar and A. K. Yadav, “Detection techniques for human safety from concealed weapon and harmful EDS,” International Review of 350 Applied Eng. Research., vol. 4, No. 1, pp. 71–76, 2014. 0.1 [3] M. E. Yavuz, A. E. Fouda, and F. L. Teixeira, “GPR signal enhancement using sliding-window space-frequency matrices,” Progress In Electro- magnetics Research, vol. 145, pp. 1–10, 2014. 450 [4] W. Shao, B. Zhou, Z. Zheng, and G. Wang, “UWB microwave imaging 0 for breast tumor detection in inhomogeneous tissue,” Eng. in Med. and 50 150 250 Biol. 27th Annual Conference, pp. 1496–1499, 2005. Distance(×∆x) [5] M. Hossain, A. Mohan, and M. Abedin, “Beamspace time-reversal microwave imaging for breast cancer detection,” IEEE Antennas Wireless (b) TR with STFT. Propag. Lett., vol. 12, pp. 241–244, 2013. [6] N. Ray, B. Saha, and M. Graham Brown, “Locating brain tumors from Normalised |Ez (x, y)| MR imagery using symmetry,” in 2007 Conference Record of the Forty- 1 First Asilomar Conference on Signals, Systems and Computers (ACSSC), pp. 224–228, Nov 2007. 50 [7] A. Zamani, S. Rezaeieh, and A. Abbosh, “Lung cancer detection using frequency-domain microwave imaging,” Electronics Lett., vol. 51, Distance(×∆y) no. 10, pp. 740–741, 2015. 0.6 [8] M. Fink, D. Cassereau, A. Derode, C. Prada, P. Roux, M. Tanter, 150 J. Thomas, and F. Wu, “Time-reversed acoustics,” Rep. Prog. Phys., vol. 63, pp. 1933–1995, 2000. [9] C.-X. Li, W. Xu, J.-L. Li, and X.-Y. Gong, “Time-Reversal detection 250 0.3 of multidimensional signals in underwater acoustics,” IEEE J. Ocean. Eng., vol. 36, pp. 60–70, Jan 2011. [10] R. Solimene, G. Leone, and A. Dell’Aversano, “MUSIC algorithms for 350 rebar detection,” J. Geophysics and Eng., vol. 10, No. 6, 2013. 0.1 [11] A. E. Fouda, F. L. Teixeira, and M. E. Yavuz, “Time-reversal techniques for MISO and MIMO wireless communication systems,” Radio Science, vol. 47, 2012. 450 [12] M. E. Yavuz and F. L. Teixeira, “Frequency dispersion compensation in 0 time reversal techniques for UWB electromagnetic waves,” IEEE Geosci. 50 150 250 Remote Sens. Lett., vol. 2, no. 2, pp. 233–237, April 2005. [13] M. E. Yavuz and F. L. Teixeira, “Full time-domain DORT for ultra- Distance(×∆x) wideband electromagnetic fields in dispersive, random inhomogeneous (c) TR. media,” IEEE Trans. Antennas Propag., vol. 54, pp. 2305–2315, Aug 2006. Fig. 9. The spatial distribution of |Ez | on z = 95∆z plane for the human [14] M. Tanter, J. L. Thomas, and M. Fink, “Focusing and steering through phantom with our compensation method, the method in [12] and without any absorbing and aberrating layers: Application to ultrasonic propagation compensation method. The scatterer and the TRA antennas are represented through the skull,” J. Acoust. Soc. Am., vol. 103, pp. 2403–2410, May by “◦” and “x” respectively. ∆x = ∆y = ∆z = 1mm. In order to present 1998. |Ez | inside the lung clearly, |Ez | on the skin, muscle and fat are removed from the figures. 8 [15] T. Folegot, C. Prada, and M. Fink, “Resolution enhancement and [40] J. Kim and Y. Rahmat-Samii, “Implanted antennas inside a human separation of reverberation from target echo with the time reversal body: simulations, designs, and characterizations,” IEEE Trans. Microw. operator decomposition,” J. Acoust. Soc. Am., vol. 113, No. 6, pp. 3155– Theory Tech., vol. 52, pp. 1934–1943, Aug 2004. 3160, 2003. [41] S. Alamri, A. AlAmoudi, and R. Langley, “Gain enhancement of [16] M. Saillard, P. Vincent, and G. Micolau, “Reconstruction of buried ob- implanted antenna using lens and parasitic ring,” Electron. Lett., vol. 52, jects surrounded by small inhomogeneities,” Inverse Problems, vol. 16, no. 10, pp. 800–801, 2016. pp. 1195–1208, 2000. [42] A. Zamani, S. Rezaeieh, and A. Abbosh, “Lung cancer detection using [17] A. M. Abduljabbar, M. E. Yavuz, F. Costen, R. Himeno, and H. Yokota, frequency-domain microwave imaging,” Electron. Lett., vol. 51, no. 10, “Frequency dispersion compensation through variable window utiliza- pp. 740–741, 2015. tion in time reversal techniques for electromagnetic waves,” IEEE Trans. Antennas Propagat., vol. 64, pp. 3636–3639, Aug 2016. [18] C. Reine, M. Van Der Baan, and R. Clark, “The robustness of seis- mic attenuation measurements using fixed- and variable-window time- frequency transforms,” Geophysics, vol. 74, No. 2, pp. WA123–WA135, 2009. [19] I. Braga and F. Moraes, “High-resolution gathers by inverse Q filtering Ammar M. Abduljabbar was born in Baghdad, in the wavelet domain,” Geophysics, vol. 78, No. 2, pp. V53–V61, 2013. Iraq, in 1985. He received the B.Sc. degree in com- [20] I. Hostens, J. Seghers, A. Spaepen, and H. Ramon, “Validation of the puter engineering from the University of Baghdad, wavelet spectral estimation technique in biceps brachii and brachioradi- Baghdad, Iraq, in 2007, and the M.Sc. degree in alis fatigue assessment during prolonged low-level static and dynamic wireless communication systems from Brunel Uni- contractions,” J. Electromyography and Kinesiology, vol. 14, no. 2, versity, London, UK, in 2009. In 2013, he started his pp. 205 – 215, 2004. PhD study in Electrical and Electronics Engineering, [21] M. E. Yavuz and F. L. Teixeira, “Ultrawideband microwave sensing and University of Manchester, UK. His current research imaging using time-reversal techniques: A review,” Remote Sens., vol. 1, interests include microwave imaging techniques and no. 3, pp. 466–495, 2009. high performance computing. [22] J. Morlet, G. Arens, E. Fourgeau, and D. Giard, “Wave propagation and sampling theory-Part I: Complex signal and scattering in multilayered media,” Geophysics, vol. 47, No. 2, pp. 203–221, 1982. [23] C. Torrence and G. P. Compo, “A practical guide to wavelet analysis,” Bulletin of the American Meteorological Soc., vol. 79, No.1, pp. 61–78, 1998. [24] M. Farge, “Wavelet transforms and their applications to turbulence,” Mehmet E. Yavuz received the B.S. degree from the Annu. Rev. Fluid Mech, vol. 24, pp. 395–457, 1992. Middle East Technical University (METU), Ankara, [25] H. Q. Ngo, E. G. Larsson, and T. L. Marzetta, “The multicell multiuser Turkey in 2001, M.S. degree from the Ko University, MIMO uplink with very large antenna arrays and a finite-dimensional Istanbul, Turkey in 2003 and the Ph.D. degree from channel,” IEEE Trans. Commun., vol. 61, pp. 2350–2361, June 2013. the Ohio State University (OSU) in 2008, all in [26] A. Taflove and S. Hagness, Computational Electromagnetics. The finite- electrical engineering. He has worked as a Princi- difference Time-domain method. Artech House, 2005. ple Scientist at the Halliburton Energy Services in [27] Y. Wang, “Inverse Q-filter for seismic resolution enhancement,” Geo- Houston, TX, as a Visiting Assistant Professor at the physics, vol. 71, pp. V51–V60, 2006. Electrical Engineering Department of the American [28] C. Moss, F. Teixeira, Y. Yang, and J. A. Kong, “Finite-difference time- University of Sharjah in the U.A.E and as a Research domain simulation of scattering from objects in continuous random Scientist at Intel Corporation in Hillsboro, Oregon. media,” IEEE Trans. Geosci. Remote Sens., vol. 40, pp. 178–186, Jan His current research interests include computational electromagnetics for 2002. industrial applications, radar-based signal processing and microwave imaging [29] P. Debye, Polar Molecules. New York: Dover, 1929. and detection. [30] M. Abalenkovs, F. Costen, J.-P. Bérenger, R. Himeno, H. Yokota, and M. Fujii, “Huygens subgridding for 3-D frequency-dependent finite- difference time-domain method,” IEEE Trans. Antennas Propag., vol. 60, No. 9, pp. 4336–4344, September 2012. [31] RIKEN, Wako Saitama, Japan, “Media parameters for the Debye relax- ation model.” http://cfd-duo.riken.jp/cbms-mp/. [Online; accessed 26- April-2016]. Fumie Costen (M’07) received the B.Sc. de- [32] J.-P. Bérenger, “Numerical reflection from FDTD-PMLs: A comparison gree, the M.Sc. degree in electrical engineering and of the split PML with the unsplit and CFS PMLs,” IEEE Trans. Antennas the Ph.D. degree in Informatics, all from Kyoto Propag., vol. 50, No. 3, pp. 258–265, 2002. University, Japan. From 1993 to 1997 she was [33] W. Weber, Pharmacogenetics. Oxford University Press, Apr 2008. with Advanced Telecommunication Research Inter- [34] N. Emaminejad, W. Qian, Y. Guan, M. Tan, Y. Qiu, H. Liu, and national, Kyoto, where she was engaged in research B. Zheng, “Fusion of quantitative image and genomic biomarkers to on direction-of-arrival estimation based on Multiple improve prognosis assessment of early stage lung cancer patients,” IEEE SIgnal Classification algorithm for 3-D laser micro- Trans. Biomed. Eng., vol. 63, pp. 1034–1043, May 2016. vision. She filed three patents from the research in [35] T. Wuren, T. Takai, M. Fujii, and I. Sakagami, “Effective 2-debye-pole 1999 in Japan. She was invited to give 5 talks in fdtd model of electromagnetic interaction between whole human body Sweden and Japan during 1996-2014. From 1998 and uwb radiation,” IEEE Microw. Wireless Compon. Lett., vol. 17, to 2000, she was with Manchester Computing in the University of Manch- pp. 483–485, July 2007. ester, U.K., where she was engaged in research on metacomputing and has [36] C. Gabriel, S. Gabriel, and E. Corthout, “The dielectric properties of been a Lecturer since 2000. Her research interests include computational biological tissues: I. Literature survey,” Physics in Med. and Biol., electromagnetics in such topics as a variety of the finite difference time vol. 41, pp. 2231–2249, 1996. domain methods for microwave frequency range and high spatial resolution [37] R. W. L. S. Gabriel and C. Gabriel, “The dielectric properties of and FDTD subgridding and boundary conditions. She filed a patent from biological tissues: II. Measurements in the frequency range 10 Hz to the research on the boundary conditions in 2012 in the U.S.A. Her work 20 GHz,” Physics in Med. and Biol., vol. 41, pp. 2251–2269, 1996. extends to the hardware acceleration of the computation using general-purpose [38] Y. Chen and P. Kosmas, “Detection and localization of tissue malignancy computing on graphics processing units, Streaming Single Instruction Multiple using contrast-enhanced microwave imaging: Exploring information Data Extension and Advanced Vector eXtentions instructions. Dr. Costen theoretic criteria,” IEEE Trans. Biomed. Eng., vol. 59, pp. 766–776, received an ATR Excellence in Research Award in 1996 and a best paper March 2012. award from 8th International Conference on High Performance Computing [39] C. M. Furse and A. Chrysler, “A history & future of implantable and Networking Europe in 2000. antennas,” in 2014 IEEE Antennas and Propag. Soc. International Symposium (APSURSI), pp. 527–528, July 2014. 9 Ryutaro Himeno received his Doctor of Engineer- ing degree from the University of Tokyo in 1988. In 1979, he joined Central Research Laboratories, Nissan Motor Co., Ltd., Yokosuka, Japan, where he has been engaged in the research of apply- ing Computational Fluid Dynamics analysis to the car aerodynamic development. In 1998, he joined RIKEN and is the director of Advanced Center for Computing and Communication and had been the deputy program director of the Next Generation Computational Science Research Program at RIKEN till April, 2013. He is also a visiting professor at Hokkaido University, Kobe University and Tokyo Denki University. He currently studies Computational Bioengineering, High Performance Computing and blood flows of human bodies. He was a winner of Nikkei Science, Computer Visualization Contest in 2000 and Scientific Visualization Contest in 1996, and received JSME Computational Mechanics Division Award in 1997 and JSME Youth Engineer Award in 1988. He has also received the Paper Award by NICOGRAPH in 1993, Giga FLOPS Award by CRAY Research Inc. in 1990 and other awards. Hideo Yokota received his Doctor of Engineer- ing degree from the University of Tokyo in 1999. In 1993, he joined Higuchi Ultimate Mechatronics Project, Kanagawa Academy of Science and Tech- nology, Kawasaki, Japan. In 1999, he joined RIKEN and is the contract researcher of Computational biomechanics unit. 2004-2012 Bio-research Infras- tructure Construction Team Leader, VCAD system research program, RIKEN. 2007-2012 Cell-scale Research and Development Team Leader,Research Program for Computational Science, RIKEN. 2013- present Image Processing Research Team Leader, Center for Advanced Photonics,RIKEN. He is also a visiting Professor at Hokkaido University, Kobe University, Tokai University and the Tokyo University of Agriculture and Technology. He currently studies biomedical imaging and image processing to the biomedical simulation. Bioimaging Society, Best Image Award, 2005. The Commendation for Science and Technology by the Minister of Education, Culture, Sports, Science and Technology , Young Scientists Prize in 2008.