Real-Time Millimeter-Wave Imaging With Linear Frequency Modulation Radar and Scattered Power Mapping

We present a novel real-time image reconstruction method processing linear frequency modulated (LFM) signals. The method exploits the principles of the Fourier-space scattered power mapping (F-SPM). We show that F-SPM, originally developed for frequency-domain signals, can be easily modified to process time-domain data with the same reconstruction speed and image quality. To facilitate validation, we have developed an in-house time-domain radar simulator, which generates synthetic LFM data much faster than full-wave time-domain simulations, which are prohibitively slow. The new image-reconstruction method is validated through synthetic data generated by the radar simulator as well as experimental data acquired with off-the-shelf millimeter-wave (77 to 81 GHz) LFM radar. Comparisons in terms of reconstruction speed and accuracy are carried out with the method of microwave holography, which is deemed the fastest image-reconstruction method for LFM radar.

Imaging radars produce a 2-D or a 3-D image of a target (or scene).Qualitative images depict the target's reflectivity, i.e., the intensity of the scattering within its volume, whereas quantitative images depict the target's permittivity composition.The image-reconstruction algorithms depend on the type of data the radars provide.The frequency-modulated continuouswave [27], [28], [29] and the ultrawide-band (UWB) pulsed radars [30] provide time-domain data, and both are common in the microwave (low-GHz) frequency ranges.Their advantage is faster measurement compared to the wide-band frequencysweep (or stepped-frequency) systems.However, at mm-wave frequencies, pulsed radar is currently impractical due to the limitations of the direct time-sampling technology and its excessive cost.On the other hand, the linear frequency modulated (LFM) radar down-converts the received signal to the beat-frequency (kHz to MHz) range [31], where real-time sampling is performed by low-cost analog-to-digital converters.For this reason, LFM radars are currently the most common low-cost option in the mm-wave frequency range.
Most of the image-reconstruction algorithms developed for microwave imaging rely on coherent stepped-frequency continuous wave (SFCW) measurements, which can be time-consuming when taking many frequency samples across a wide frequency range.Also, at mm-wave frequencies, the equipment is costly.Since LFM radars offer faster and more cost-effective option [32], there is great interest in developing fast image-reconstruction methods to process the LFM data.
Back-projection is a classic synthetic aperture radar image-reconstruction approach, and back-projection algorithms (BPAs) have been developed for LFM mm-wave imaging [33], [34], [35], [36], [37], where they operate directly on the time-domain data.They compute the round-trip delays in the background medium between each imaged pixel and the receiving/transmitting antenna pairs in order to obtain a coherent sum of all measured signals specific to a pixel.The image depicts the energy of these pixel-specific sums, indicating the scattering intensity (or reflectivity) within the imaged scene.
Fourier-based imaging is a computationally efficient alternative to back-projection.The approach is often referred to as microwave holography [19], [38]. 1 The hallmark of microwave holography algorithms (MHAs) is the image reconstruction in the spatial-frequency domain (k-space).This requires the 2-D FT of the data.The MHAs, originally developed for frequencydomain data, are also applicable to LFM data.However, the latter application requires an approximation of the downconverted (de-chirped) signal that neglects the second-order time-delay term [41].
To achieve 3-D image reconstruction in k-space, most MHAs (see, e.g., [20], [41], [42]) rely on an analytical range-migration model, which provides the link between the frequency (ω) dependence of the data and the range (or depth) dependence of the image, along with Stolt's interpolation.The k-space result is then cast back to 3-D real space via the 3-D inverse Fourier transform (IFT).This approach is known as the range-migration algorithm (RMA).Stolt's interpolation is by far the most computationally intensive task but this drawback has been overcome by recent MHAs, which avoid this interpolation, e.g., the range-stacking algorithms [43], [44] and the near-field MHAs [21], [40], [45], [46].They perform the inversion in the mixed (k x , k y , z) space, where k x and k y are the Fourier variables corresponding to the lateral coordinates x and y, whereas z is the range.
Stolt's interpolation and the FTs introduce numerical errors, which may lead to image artifacts in MHA reconstructions unless filtering is applied [46].BPAs do not suffer from such artifacts.It is shown in [47] that the RMA yields 2-D images with better cross-range resolution compared to the BPA, but offset errors due to Stolt's interpolation may occur.In 3-D imaging, however, the BPA seems to offer better resolution.Overall, the BPAs are significantly slower than the MHAs [44], [47] but they are less prone to image artifacts [9], [48], [49], [50].
To improve the image accuracy and to enable quantitative reconstruction in near-field imaging, the measured system (or data) point-spread function (PSF) is used [40], [45], [51], [52], [53] in place of the analytical PSFs used in far-field imaging.The measured PSF provides the system-specific quantitatively accurate resolvent kernel of the linearized scattering model.Using measured PSFs, quantitative imaging of dielectric objects has been demonstrated by algorithms such as quantitative microwave holography (QMH) [21], [54], [55] and Fourier-space scattered power mapping (F-SPM) [56], [57].
Here, we propose a novel image-reconstruction method for processing LFM signals, which we refer to as Fourierspace scattered power mapping in the time domain responses (FSPM-TD).It is based on the F-SPM method, originally developed for SFCW data [56], [57], and it operates directly on time-domain data.The data spatial dependence is treated in k-space, leading to superior computational speed, shown to be better than that of the existing k-space algorithms.At the same time, unlike these algorithms, the FSPM-TD algorithm does not neglect the second-order time-delay term in the LFM signal.The algorithm is validated through simulated data (generated by an in-house radar simulator) as well as measured data obtained with an off-the-shelf LFM platform [58].Its speed and accuracy are compared with the fast QMH algorithm [21], [54], which does not employ Stolt's interpolation.
Next, Section II introduces the FSPM-TD method and its implementation with LFM data.Sections III and IV present validation examples with synthetic and measured data, respectively.Conclusion is drawn in Section V.

A. Fourier-Space Scattered Power Mapping With Time-Domain Responses
Scattered power mapping (SPM) is a well-established method for fast (real-time) microwave imaging [56], [57], [59].It is a direct reconstruction method since it relies on a linearized model of scattering.With quantitatively accurate (measured) system PSFs, it can also reconstruct images of the real and imaginary parts of the object's complex permittivity (quantitative images).The method operates on frequencydomain signals.The most computationally efficient SPM algorithm is F-SPM [57], which performs the inversion in k-space.Since this algorithm serves as the basis for the current development, it is summarized in Appendix.
The SPM is a two-stage inversion procedure.To understand its new implementation with time-domain signals, we start with its formulation in real (x, y, z) space.
With frequency-domain responses, the first SPM stage constructs a complex-valued qualitative image M ω (r ′ ) (scatteredpower map, or simply, map) of an object as [51], [56] M where ω is frequency, ζ indicates an antenna pair associated with a response, N T is the number of responses acquired at each observation (receiver) position r on the aperture S a , r ′ is a position in the imaged domain, S sc ζ (r, ω) is the scattered portion of the response measured with the object in place, H sc ζ (r, ω, r ′ ) is the scattering response measured with a point scatterer at r ′ in the background medium (the system PSF), and * indicates conjugation.It is clear from (1) that the object's map M ω (r ′ ) is an inner product of the measured responses and the system PSFs in the data space spanned by r, ω, and ζ .
In the time domain, (1) can be written as where X ζ (r, τ, r ′ ) is the temporal cross correlation of S sc ζ (r, t) and H sc ζ (r, t, r ′ ) with the time shift τ , and F τ is the FT with respect to time.We next consider the integral over ω in conjunction with the FT of X ζ (r, τ, r ′ ).Assuming infinite frequency bandwidth, at any r and r ′ , we obtain The substitution of (4) into (2) results in the map of the object expressed in terms of time-domain responses Here, the scaling factor 2π has been omitted since it has no impact on the final image.In conclusion, the first stage of the SPM image reconstruction can employ time-domain instead of frequency-domain data to obtain the object under test (OUT) map.The comparison of ( 5) and (1) shows that the integration over ω is replaced by that over time t.Note that with UWB radar, the temporal sequences are real, but with LFM radar systems, they are complex, i.e., at each r, , where I and Q denote the in-phase and quadrature receiver outputs.With complex time-domain signals, the conjugation in (5) matters.
The direct computation of the OUT map M(r ′ ) with ( 5) is slow.It can be carried out much faster in the 2-D k-space under the assumption of a homogeneous background, where the dependence of the PSF on r and r ′ reduces to a subtraction, H sc ζ (r−r ′ , t).Note that the time variable t is also a function of (r − r ′ ) through its dependence on the distance R Rx = |r − r ′ | between the imaged point r ′ and the measurement (receiver) point r.Thus, the integration over r ∈ S a in (5) becomes a 2-D cross correlation of the OUT response and the system PSF in the lateral (cross-range) coordinates.In 2-D k-space, this cross correlation is a point-wise multiplication of the respective 2-D FTs.The computation is most efficient in the case of uniform sampling on canonical surfaces (planar, cylindrical) since this allows for the use of the 2-D fast Fourier transform (FFT).
Fig. 1 illustrates a single-sided multistatic measurement setup, where the scan is over a planar surface.This setup reflects all examples presented later.The measurement LFM-radar platform features 3 transmitting (Tx) and four receiving (Rx) antennas, all moving together over the acquisition plane S a at regular intervals along x and y.
Appendix describes the process of casting the OUT scattered-power map (1) for frequency-dependent data into 2-D k-space in the case of planar scanning, where the cross-range variables are x and y.The same process can be applied to the time-domain formulation of the OUT scattered-power map in (5), leading to its 2-D FT form where κ κ κ = (k x , k y ) is a point in k-space with k x and k y being the Fourier variables corresponding to x and y, Once the 2-D FT of the OUT map is computed with (6), it can be cast back in real space using 2-D inverse FT The absolute value of the so obtained OUT map |M(x ′ , y ′ , z ′ )| (usually normalized) provides a qualitative image of the object's reflectivity.However, a significant image improvement is achieved with the second SPM stage.As explained next, this stage operates directly on the k-space OUT map M(κ κ κ, z ′ ), thus bypassing the inverse FT operation in (7).
It is shown in [56], [57], [59] (for the case of frequencydomain responses) that the second SPM step provides an image with significantly improved spatial resolution compared to the OUT qualitative image (the map) obtained with (1) (real-space processing) or with the IFT of (28) (Fourier-space processing; see Appendix).It also enables the quantitative estimate of the complex permittivity of dielectric objects, provided the system PSFs are quantitatively accurate. 2Similar to the first stage, for best computational efficiency, the second SPM stage is performed in k-space.Since the processing is essentially the same as in the second stage of the F-SPM algorithm for frequency-domain data (see Appendix), only the computations relevant to time-domain responses are presented below, followed by a summary of the algorithm.
The second SPM stage operates on the OUT k-space map M(κ κ κ, z ′ n ), n = 1, . . ., N z , computed with (6).It also requires the computation of N z k-space maps, M(κ κ κ, z ′ n ; z ′′ m ), n, m = 1, . . ., N z , of the scattering probe (SP), when this probe resides at r ′′ m = (0, 0, z ′′ m ), m = 1, . . ., N z .This computation mirrors that of the OUT map; see (6).Specifically Note that the real-space maps of the point scatterers, M(x ′ , y ′ , z ′ n ; z ′′ m ), corresponding to the k-space maps computed with (8), are the image PSFs (IPSFs) resulting from the first SPM stage. 3The SP maps in (8) are independent of the imaged object and can be precomputed for faster execution of the image reconstruction.
As shown in Appendix, with the OUT and scattering-probe maps available in k-space, the second SPM stage solves a small N z × N z system of equations at each point κ κ κ to obtain the 2-D FT of the reflectivity function ρ provides a qualitative image of the object's reflectivity.A quantitative estimate of the OUT complex permittivity is possible, provided the system PSFs, H sc ζ (x, y, t; z ′ n ), scale properly with the probe's volume sp and relative-permittivity contrast ε r,sp .Then, the object's relative-permittivity contrast is computed with (36).
The proposed FSPM-TD algorithm is summarized in Fig. 2. It takes as inputs the measured OUT responses S sc ζ (x, y, t), ζ = 1, . . ., N T , and the system PSFs, H sc ζ (x, y, t, z ′ n ), n = 1, . . ., N z , the latter being obtained either through measurements, or simulations, or analytical models.Note that, in a multistatic system, each Tx/Rx antenna pair, indicated by ζ ≡ (i, j), has a dedicated system PSF H sc ζ .

B. Forward Model of Scattering With LFM Signals
The LFM radar signal is a "chirp" waveform -a sine wave of frequency that increases or decreases linearly with time.A transmitted LFM chirp is expressed as [60], [61] s Tx (t) = A Tx P t/T p cos[2π( where A Tx is amplitude, f c is the center frequency, t is the fast time (the time within a single chirp), T p is the chirp duration (pulsewidth), γ = (B/T p ) is the frequency-modulation slope (chirp rate), B is the chirp's frequency bandwidth, and The spatial impulse response h sc (r, t, r ′ ) of the LFM radar describes the scattered signal from a differential scatterer (SP), dh sc (r, t, r ′ ) = ρd h sc (r, t, r ′ ), where ρ and d are the reflectivity and volume of the probe, respectively.For static objects in a homogeneous unbounded background, the LFM-radar impulse response (i.e., its analytical PSF) is modeled as a scaled and time-delayed version of s Tx (t) where r is the Rx position, r ′ is the probe's position, and is the time delay corresponding to the distance traveled by the signal.Here, c is the speed of light whereas R Tx = |r − r Tx | and R Rx = |r − r ′ | are the distances from the Tx antenna at r Tx to the probe and from the probe to the Rx antenna, respectively.The model in (11) accounts for the signal decay due to the spherical spread of the transmitted and scattered waves through the factor (R Tx R Rx ) −1 .On the other hand, it is a greatly simplified approximation of reality since it ignores the vector nature of the electromagnetic waves, the depolarization that may occur upon scattering, the gain and dispersion of the employed antennas, etc.Note that, at each scan position, R Rx and R Tx differ, depending on which antenna pair ζ ≡ (i.j) in the multistatic system the PSF describes.
Upon reception, the scattered signal is dechirped by quadrature down-conversion to produce the beat or baseband signal, which is used for the image reconstruction.The baseband output corresponding to h sc (r, t, r ′ ) in ( 11) is the analytical system PSF [41], [60] ) Note that the signal in ( 13) is complex, where its real and imaginary parts represent the I and Q Rx outputs, respectively.
The investigation of the LFM forward model in [41] points out that the spatial resolution of the images obtained with MHAs is negatively affected if the third term of the exponent in ( 13) is not negligible.This is due to the MHAs treating the LFM time-domain signals as frequency-domain signals of equivalent frequency f ′ = f c + γ t and wavenumber Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
k ′ = 2π f ′ /c [41].The proposed FSPM-TD algorithm does not suffer from this limitation since it takes the system PSFs in any form and it does not need an equivalent frequency.Nonetheless, if the phase contribution of the term π γ τ 2 d is kept below 2.5 • , as suggested in [41], then the LFM system PSF can be approximated as Note that the phase of H sc a is now proportional to τ d , leading to an exponential term of the form e −ik ′ (R Tx +R Rx ) , which matches that in a frequency-domain response.The PSF H sc a is employed by the QMH algorithm in the examples presented later, where the FSPM-TD algorithm is compared with QMH.
The linearized forward model of scattering views the signal from an object as a superposition of the scattering emanating from all differential scatterers that make up this object.Thus, using (13), the cumulative OUT signal is modeled as where ρ(r ′ ) is reflectivity, and The forward model in (15) is the basis of the LFM radar simulator used in the synthetic experiments presented next.
The FSPM-TD image reconstruction employs the analytical system PSF (13) with both synthetic and measured data.

III. VALIDATION WITH SYNTHETIC DATA
An LFM simulator is implemented in MATLAB [62] using the scalar scattering model (15) for the case of planar scanning with multistatic measurements.The PSF employs (13).Note that this PSF, along with (15), inherently assume scattering in an unbounded medium.The multistatic scenario allows for using any number of Rx and Tx antennas, which remain in a fixed configuration during the scan.Thus, at each scan position r, the number of acquired responses is N T = N Tx N Rx , where N Tx and N Rx are the number of Tx and Rx antennas, respectively.To match our experimental setup employing a single LFM board [58], the scans are single-sided.
In each synthetic experiment, the LFM imaging-system parameters are first set.For a list of these parameters, refer to Table I.Then the system-calibration simulations are performed using (13).These emulate the PSF measurements with an SP located at the center of each imaged slice (0, 0, z ′ n ), n = 1, . . ., N z .The SP volume sp = d is set equal to that of the imaged voxel.This process provides the system PSFs, H sc ζ (x, y, t; z ′ n ), ζ = 1, . . ., N T .This is followed by the computation of the OUT data S sc ζ (x, y, t) using (15).Table I summarizes the system parameters employed in all presented examples, except for the sampling step along x and y, which is 1 mm for the IPSF study in Section III-A.The first six system parameters describe the radar itself.These have been chosen to match the settings of the LFM radar [58] used in the experiments.The spatial sampling step x = y = ⊥ (2 mm in Table I) is always chosen to be somewhat smaller than the expected cross-range resolution limit δ ⊥ .This limit is given by δ ⊥ = (λ c /4 sin α) [51], [63] where λ c = c/ f c and α is the maximum viewing angle of the scan, α = min θ a , avg(0.5θh , 0.5θ e ) .Here, θ a = arctan(0.5A/R) is the viewing angle provided by the aperture, with A and R being the aperture width and the range distance to the object's center, respectively.θ h,e denote the antenna half-power beamwidths in two principal planes.For example, in the simulations, the antennas are isotropic, thus α is determined by the aperture.
Provided that A = 15 cm and the target is 22.5 cm away, α ≈ 18.4 • , leading to δ ⊥ ≈ 3 mm.It is worth commenting that obtaining synthetic LFM-radar data with full-wave simulators is prohibitive slow due to: 1) the extremely long chirp signals and 2) the need to simulate a large amount of illumination (Tx) positions associated with scanning a multistatic radar system over a large aperture. 4

A. Image Point Spread Function and Spatial Resolution
In the first experiment, we image a single-voxel scatterer at 22.5 cm from the acquisition plane and obtain the IPSF of the FSPM-TD algorithm.From the IPSF, the cross-range and range resolution values are estimated and compared to the theoretical limits.Here, the scanning step is ⊥ = 1 mm, and sp = d = 1 mm 3 .The IPSF contains 29 range slices separated by z = 5 mm and centered on the z = 22.5 cm plane.Fig. 3(a) shows the IPSF slice at z = 22.5 cm.The IPSF width at −4 dB indicates the spatial resolution in the respective direction [51].Fig. 3(b) shows the line cuts of the IPSF along x and y at z = 22.5 cm, and along the line cut along z at x = y = 0.The results indicate cross-range resolution of 3 mm, which agrees with the theoretical limit δ ⊥ ≈ 3 mm.The range resolution is obtained as 22 mm whereas the theoretical limit is [63] δ z = (c/2B) ≈ 37 mm.    .This quantitatively accurate result is expected since both the PSFs and the OUT data are generated by the same "measurement system" emulated by the LFM radar simulator.Also, the radar simulator employs the simple superposition scattering model in (15), i.e., it does not model the mutual coupling and multiple scattering, which occurs in reality, and which is the main reason for image degradation in quantitative imaging.This example highlights the advantage of measuring the PSFs with the same system used to measure the OUT.Unfortunately, measuring the PSFs is not always possible, especially in far-zone measurements, where the SP signal may be too weak to detect with a sufficient signal-to-noise ratio.

B. Three-dimensional Imaging With Synthetic Data
To quantify the image quality in Fig. 5, the structural similarity (SSIM) index is computed [64].The SSIM index ranges from 0 to 1, where 1 indicates perfect similarity and 0 indicates no similarity.Here, the SSIM is 0.9515 in the F-shape slice and 0.9104 in the bar-shape slice.Additionally, to evaluate the precision of the reconstructed permittivities, the root-mean-square error (RMSE) is calculated as [51] where ε r (r ′ n ) is the true distribution, ε r (r ′ n ) is the reconstructed distribution, and N v is the number of voxels.Here, the image RMSE is 0.0219.
The accuracy of the FSPM-TD reconstruction is compared with that of the QMH method.The QMH method is a fast MHA, which does not employ Stolt's interpolation.The QMH images are not shown here since there is no visible difference with those in Fig. 5.To compare better the two reconstructions, an RMSE is computed where the FSPM-TD result provides the reconstructed distribution ε r (r ′ n ) whereas the QMH result Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
The two algorithms are also compared in terms of running time.The FSPM-TD algorithm takes about 2.9 s whereas the QMH algorithms take about 5 s.Note that both algorithms are implemented in MATLAB without any code optimization and using the same direct and inverse fast FT (FFT) function calls.To understand the reason for the faster performance of FSPM-TD, we first point out that both methods share common initial steps, which involve the 2-D FFT of the PSFs and the OUT data.Also, both of them employ 2-D inverse FFT (slice by slice) on the reconstructed k-space reflectivity function.However, they differ in solving their respective linear systems of equations in k-space.FSPM-TD solves square N z × N z systems (31) whereas QMH solves tall (N ω N T ) × N z systems [46], [54].The number of slices N z rarely exceeds 10, and the number of response types N T is also on the order of 1 to 10.However, in the image reconstruction with QMH, the number of equivalent frequencies N ω equals that of the time samples N t , and that is on the order of ∼10 2 to ∼10 3 , depending on the length of the employed chirp sequence.In this example, which spans a single chirp, N t = 512.To solve the tall system of equations, QMH uses MATLAB's pseudo-inverse (pinv) function, which employs a singular value decomposition approach [65], the computational complexity of which is O((N t N T ) 2 N z + N 3 z ).FSPM-TD, on the other hand, can employ either LU decomposition or pseudo-inverse solvers.In either case, its computational complexity is about O(N 3 z ).It is now clear that the computational advantage of the FSPM-TD algorithm arises when N z < N t N T .

IV. VERIFICATION WITH MEASUREMENTS
The experiments are carried out in a planar raster-scanning chamber shown in Fig. 6(a).With the advent of system-on-chip mm-wave sensing technology, the market now offers various off-the-shelf radar modules.Here, we use the IWR1443Boost evaluation module [58] along with the real-time data-capture adapter board DCA1000EVM [66].The mm-wave sensor is equipped with three Tx and four Rx antennas as shown in Fig. 6(b) (from [58]).The LFM transceivers can accommodate up to 4 GHz bandwidth from 77 to 81 GHz.The configuration of the radar system is done via the TI mmWave Studio software suite.This includes activating/deactivating Tx and Rx channels, the choice of the chirp sequence, and the chirp settings.The chosen system parameters match those in Table I.All twelve radar channels, formed by the three Tx and four Rx on-board antennas, are used in the experiments.The OUT data S(r, t) are captured through measurements employing all available radar channels.
The radar module is mounted at the top of the chamber [see Fig. 6(a)] and it is stationary while the platform carrying the imaged object moves laterally along a raster-scan path with increments x = y = 2 mm.At each grid point, the radar takes measurements for about 2 s, during which time the platform does not move.
The relative positions of the Tx and Rx antennas are needed to calculate the time delay τ d associated with an imaged point Fig. 6.Photos of (a) acquisition chamber, and (b) on-board antenna array of the IWR1443 sensor (from [58]).and each Tx/Rx antenna pair.The center-to-center spacing between the Rx elements is 1.9 mm whereas between the Tx elements it is 3.8 mm.The center-to-center spacing from rx4 to tx1 is 4.75 mm [see Fig. 6(b)].The coordinate system is aligned so that Rx antenna #4 [rx4 in Fig. 6(b)] is at (0, 0, z) at the start of the scan, where z is the distance from the radar printed circuit board (PCB) to the center of the imaged object.
In all experiments presented next, background subtraction is not used to extract the scattered portion of a response from the total measured response.This subtraction is mandatory in near-field imaging and especially when forward-scattering signals are employed because the incident-field portion of the total object response is strong.Here, background de-embedding is unnecessary since the background signals are negligible compared to the backscattering from the objects.

A. System Calibration
The measurements are susceptible to various types of uncertainties, of which the internal system delay t sys is the most detrimental to the ranging information carried by the scattered signal.Aside from ignoring signal dispersion due to the antennas, the PSF models in (13) or (14) assume that: 1) the signals at the Rx antenna terminals arrive at the input of the down-converting mixer without delays and 2) the signals transmitted by the Tx antennas are the same as those submitted to the mixer.The first assumption is not true due to the coplanar-waveguide (CPW) transmission lines connecting the Rx antennas to the radar chip [see Fig. 6(b)] along with signal pathways inside the chip.Similarly, the second assumption is not true due to the CPW and ON-chip interconnects to the Tx antennas.The cumulative effect of the delays along interconnects is represented by a constant t sys , which must be added to the signal-delay time variable τ d in the PSF model.The calibration method aims at extracting t sys .
The calibration measurement setup is illustrated in Fig. 7.It employs a 5 × 5 cm 2 copper plate serving as an ideal reflector, which lies parallel to the radar PCB and centered on the boresight of the Tx/Rx antenna set.The system delay does not depend on the distance between the radar and the plate but measurements at various distances should be carried out to verify the extracted τ sys .Here, distances anywhere between 20 and 40 cm have been employed, which are within the possible ranges in the chamber. The The response in (18) accounts for the phase reversal (the minus sign) upon reflection from the copper plate.The integration over the plate's surface is ignored since the plate's lateral size is much smaller than the range distance.The time delay τ d,ζ , needed to compute H sc ζ (t), is obtained from R Txζ and R Rxζ using (12), where R Txζ = |r Txζ − r 0 | and R Rxζ = |r ζ − r 0 | are the distances from the plate's center (at r 0 ) to the Tx and Rx antennas, respectively.
The alignment is done in the frequency ( f ) domain.The FT of S a ζ (t), Sa ζ ( f ), has a magnitude spectrum which peaks at the frequency f p = −γ τ d,ζ since [60], [67] Here, sinc(x) ≡ (sin(x)/x).The peak frequency f p is a crucial marker for the target's range.If the target motion is negligible (zero Doppler shift), f p is proportional to τ d,ζ [60], [67], and, therefore, to the distance to the target; see (12).We exploit this LFM signal feature to find t sys .First, we generate the time sequence of S a ζ (t) with the same sampling step and length as that of S m ζ (t).FFT is The internal system delays for the LFM radar employed here have been determined for all Tx/Rx channels first with the copper plate placed 355 mm away from the radar.The values are the same across all radar channels: τ sys = 0.26969 ns.The calibration has been repeated for various range positions of the copper plate and t sys has been confirmed to be the same.
To verify the calibration, the so obtained t sys is applied to the analytical response S a ζ (t) in ( 18) by replacing τ d with τ d +t sys .The magnitude spectrum of the calibrated analytical response showing the peak-frequency alignment.It is worth noting that although the peak-frequency misalignment between the uncalibrated analytical response and the measured response may appear small, it actually corresponds to about 4 cm difference in distance.Without calibration, this difference results in extremely unfocused images with the measured data.
Similarly, t sys is used to calibrate the analytical system PSFs in (13) and obtain H sc cal,ζ (r, t, r ′ ) by replacing τ d,ζ with τ d,ζ + t sys for each Tx/Rx antenna pair.These calibrated PSFs are used in the image reconstruction with the measured data.

B. Imaging Experiments
The initial validation is conducted using the same 3-D F-shape/bar-shape object described in Section III-B.However, in this experiment, copper tapes of thickness 1 oz (34.8 µm) and width 4 mm are used in crafting the shapes, as shown Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. in Fig. 9(a).Each arm of the F and bar shapes comprises 4 stacked layers of copper tape to ensure large reflectivity.The F-shape contains arms of lengths 5 and 2.5 cm.The bar shape is 5 cm long.The shapes are affixed to paper sheets [see Fig. 9(a)].The imaging setup is shown in Fig. 9(b).The F-shape is placed at the reference plane z = 0, which is 22.5 cm away from the radar.The bar shape is at z = −12 cm (34.5 cm away from the radar).The paper sheets holding the shapes are placed on Styrofoam (ε r ≈ 1.175) slabs of thickness 12.7 mm.The system and sampling parameters are listed in Table I.
In the image reconstruction, we employ the calibrated analytical PSFs computed for a cubical SP of volume d = 1 mm 3 positioned at the three slices of interest: z = −120, 0, 120 mm.Therefore, in this and all subsequent experiments, the images are qualitative.
We briefly mention that we have attempted the measurement of the PSFs using probes of size (λ c /4) ≈ 1 mm.Unfortunately, at the employed ranges, the scattering from such probes is too weak to rise above the uncertainty of our measurement system.
The FSPM-TD reconstructed image of the normalized reflectivity is presented in Fig. 10(a).The F and bar shapes are reconstructed with good structural accuracy.Even a small air-gap in the top horizontal arm of the F-shape [see Fig. 9(a)] is reconstructed well in the middle layer (z = 0) of Fig. 10(a).The top layer at z = 120 mm [air in the setup shown in Fig. 9(b)] accurately shows the absence of objects.
Regarding the image quality, the image slices are well resolved along range, i.e., there is no "range bleeding".This is expected since the slice separation is well beyond the range resolution limit of δ z ≈ 37 mm.However, artifacts are present in the bottom layer where the bar shape is.These are due to reflections coming from the plastic rods as well as the scanning platform itself.Note that, as shown in Fig. 9(b), the F-shape at the top is far from any structural components of the scanner and its image in the slice z = 0 shows practically no artifacts.To compare the accuracy of FSPM-TD images with an MHA algorithm, the QMH image is provided in Fig. 10(b).The two images are almost identical, validating the accuracy of the FSPM-TD algorithm.We note that the QMH algorithm, although using the approximate PSF in ( 14), provides focused images since ( 14) is accurate in this example.With a frequency-modulation slope of γ = 72.42× 10 12 Hz/s and distance to target of about 35 cm, the phase contribution of the π γ τ 2 d term does not exceed 0.09 • , which is well below the limit of 2.5 • recommended in [41].As in the example with the synthetic data, since all sampling rates are similar, the FSPM-TD algorithm is faster than QMH (2.9 s versus 5.0 s).
The second imaging experiment addresses a scenario featuring realistic items.The object includes a metallic key, a penny, and a liquid lipstick, see Fig. 11(a).In an initial experiment, all three objects are lying on a Styrofoam sheet, which is 22.5 cm away from the radar.In a second experiment, the same objects are enclosed within a toy bag shown in Fig. 11(b) and the bag  is placed on the same Styrofoam sheet.The radar and sampling parameters remain the same as those in Table I.
Fig. 12(a) shows a 2-D image of the unobstructed key, penny, and lipstick experiment.The image is obtained by a maximum value projection of six slices within a volume of 2.5 cm range thickness, i.e., confined between the planes z = 19 cm and z = 21.5 cm.The reason for presenting the 2-D projection images is that the objects have different thicknesses and their reflectivity is best represented in a projection.The image in Fig. 12(a) shows all items with excellent resolution and no visible artifacts.The 2-D projection image of the same objects concealed in the bag is shown in Fig. 12(b), following the same procedure.It is clear that the bag has a negative impact on the structural accuracy of the reconstructed objects, likely due to the fact that the materials from which the bag is made are not entirely transparent to the mm-wave radiation.In fact, the outline of the bag is visible in Fig. 12(b).Moreover, the hello-kitty plush toy attached to the bag is relatively large and thick [see Fig. 11(b)].
The last experiment reported here addresses the realized cross-range resolution of the experimental setup.To this end, two benchmark targets are fabricated in PCB technology, each consisting of five copper strips of thickness 2 oz (69.6 µm) and length 2.5 cm; see Fig. 13.The PCBs employ FR-4 substrates (ε r ≈ 4.3) of size 8 × 8 cm 2 .The strip width in Benchmark #1 is 3 mm whereas in Benchmark #2 it is 2 mm.In both benchmark targets, the strip edge-to-edge spacing varies from 2 to 5 mm at 1 mm increment.The reconstructed 2-D image of Benchmark #1 is shown in Fig. 14  The proposed inversion algorithm employs a linearized integral scattering model whose kernel (the radar's spatial impulse response) is the PSF.For fast k-space inversion, the assumption of a uniform unbounded background medium is made, which renders the scattering model a 2-D convolution in the lateral coordinates.Unlike conventional direct-inversion methods, which rely on analytical PSFs, the FSPM-TD algorithm can operate with analytical, simulated, or measured PSFs without any modifications and with no impact on its speed.Since measured PSFs enable near-field and quantitative imaging, this capability is an important advantage.
In its first inversion stage, the FSPM-TD algorithm is a projection algorithm, since it employs the inner product of the measured responses with the system PSFs to produce a reflectivity image.In its second inversion stage, it performs image enhancement (and quantitative imaging, if the PSFs are measured) by deconvolving the object's reflectivity image with that of the SP This approach of projection forming an inner Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.from ( 21) and ( 22) as [57] S sc ζ (r, ω) ≈ where ρ(r ′ ) = ε r (r ′ )/( ε r,sp sp ) is termed the object's reflectivity function.H sc ζ (r, ω) is the response to a probe at the center of the imaged volume, r ′ = 0.With the data, S sc ζ (r, ω), and the PSFs, H sc ζ (r, ω), available, the F-SPM method solves (23) for ρ(r ′ ) through a computationally efficient two-stage procedure.The first SPM stage constructs the 3-D scattered-power map M(r ′ ) of the OUT as the inner product of the data and the system PSFs.In the case of a planar scan at z = z, the measurement position is given by r = (x, y, z), and the explicit map expression is M(r ′ ) in ( 25) is a 2-D cross correlation in x and y.The most efficient way of computing it is in 2-D Fourier space (k x , k y ), where k x and k y are the Fourier variables corresponding to x and y, respectively.For brevity, a point in Fourier (or k) space is denoted as κ κ κ = (k x , k y ).The k-space processing requires the FTs of the data at all frequencies, The second SPM stage is also performed in k-space [56], [57].In addition to the OUT map in (28), it requires the 2-D FTs of the SP maps.These are obtained analogously to ( 28 With the OUT and SP maps available, the 2D FT of the reflectivity function, ρ(κ κ κ, z ′ n ) = F 2D {ρ(x ′ , y ′ , z ′ n )}, is extracted using the linear map relation [57]: The N z equations in (30) form a small N z × N z system of equations at each k-space point written as: Since M(κ κ κ) is a small square matrix, (31) can be efficiently solved using LU decomposition.The real-space reflectivity function ρ(x ′ , y ′ , z ′ n ) is recovered via the inverse 2D FT of ρ(κ, z n ): The plot of |ρ(x ′ , y ′ , z ′ n )| provides a qualitative image of the object's reflectivity.Quantitative image is also possible, provided the system PSFs scale properly with the probe's volume sp and relative-permittivity contrast ε r,sp .As per (24), the quantitative estimate of the object's relative-permittivity contrast is obtained as

Fig. 1 .
Fig. 1.Illustration of the single-sided multistatic measurement setup with a planar aperture denoted as S a .The red triangles and the blue points represent Tx and Rx positions, respectively.The response acquired with the jth Tx antenna ( j = 1, 2, 3) and the ith Rx antenna (i = 1, . . ., 4) is denoted by ζ ≡ (i, j).The array of 3 Tx and 4 Rx antennas moves along a raster-scan path indicated by the gray dashed line.Thus, the positions of the Tx antennas, r Tx, j , and the Rx antennas, r i , are all incremented with a common sampling step along x and y during the scan.The imaged position is denoted as r ′ .respectively, z ′ is the range position of an imaged slice, t k = k t is the kth time sample as determined by the time-sampling step t, and N t is the number of time samples.Ssc ζ (κ κ κ, t) and H sc ζ (κ κ κ, t, z ′ ) are the 2-D FTs of S sc ζ (x, y, t) and H sc ζ (x, y, t, z ′ ), respectively, where H sc ζ is the system response acquired with a scattering probe (point scatterer) at position (x ′ , y ′ , z ′ ) = (0, 0, z ′ ).It should also be pointed out that the range variable z ′ belongs to a discrete set of range slices, z ′ n = n z, n = 1, 2, . . ., N z , where z is the range step size.The implementation of (6) in the case of LFM radar is detailed later in Section II-B.Once the 2-D FT of the OUT map is computed with (6), it can be cast back in real space using 2-D inverse FT

A 3 -
D object is implemented in the LFM radar simulator as shown in Fig.4.All structural components are built of

Fig. 3 .
Fig. 3. FSPM-TD reconstructed image of a cubical probe 1 mm on a side, at the range distance 22.5 cm.(a) Two-dimensional IPSF at the z = 22.5 cm in terms of normalized reflectivity ρ.(b) Range and cross-range profiles of the IPSF.

Fig. 5 .
Fig. 5. Reconstructed images of the real and imaginary parts of the relative permittivity of the object in Fig. 4 using synthetic data.(a) Bar shape in the z = −120 mm slice (true permittivity ε r = 3).(b) F-shape in the z = 0 slice (true permittivity ε r = 1.8).(c) Slice at z = 120 mm where there are no embedded targets.
echo signal S m ζ (t) is captured by every (ζ ) Tx/Rx antenna pair.The goal is to align the measured signals S m ζ (t) with an analytical model S a ζ (t) based on (13), namely

Fig. 8 .
Fig. 8.Comparison between the magnitude spectra of a measured calibration response Sm ζ ( f ) and the respective analytical response Sa ζ ( f ) for a copper plate placed 355 mm away from the radar.Sa ζ ( f ) is the analytical result before calibration whereas Scal ζ ( f ) is the result after calibration.

Fig. 9 .
Fig.9.Imaging setup for the 3-D reconstruction of two copper-tape objects depicted.(a).Both shapes consist of four layers of copper tape carefully applied to a paper surface.The F-shape is positioned at the uppermost layer, while the tilted bar shape resides at the bottom layer, as illustrated.(b) F-shape plane is situated 22.5 cm away from the radar, whereas the plane of the bar shape extends an additional 12 cm.All dimensions are in cm.

Fig. 10 .
Fig. 10.Reconstructed 3-D images in terms of normalized reflectivity ρ enabling.(a) FSPM-TD method.(b) QMH method using measured data with the F-shape/bar-shape object.The bar-shape and the F-shape are correctly located at z = −120 mm and z = 0 mm, respectively.The top layer at z = 120 mm correctly shows a slice without any targets.

Fig. 11 .
Fig. 11.Photos of (a) key, penny, and liquid lipstick lying on a Styrofoam sheet.(b) Small toy bag used to "conceal" the key, the penny, and the lipstick.Dimensions are in mm.

Fig. 12 .
Fig. 12. FSPM-TD image of the normalized 2-D projection of the reflectivity ρ.(a) Key, penny, and lipstick on the scanning platform.(b) Same objects inside the bag.
(a).All strips are resolved well, even the two strips with a 2 mm spacing.On the other hand, the image of Benchmark #2, shown in Fig.14(b), fails

Fig. 13 .
Fig. 13.Photos of benchmark PCB targets composed of copper strips of 2 oz (69.6 µm) thickness and length of 2.5 cm.(a) Benchmark #1 with strip width of 3 mm.(b) Benchmark #2 with strip width of 2 mm.In both benchmark targets, the strip edge-to-edge spacing varies from 2 to 5 mm at 1 mm increment.The PCBs employ FR-4 (ε r ≈ 4.3) substrates of size 8 × 8 cm 2 .

TABLE I SYSTEM
PARAMETERS IN THE IMAGING EXPERIMENTS WITH SYNTHETIC mm-WAVE LFM DATA