II. RECEIVED SIGNAL MODEL OF A ROTATING OBJECT

A typical man-made object consists of a small number of dominant scatterers. For the narrow-band radar, echoes of the k-th scatterer in the slow time domain can be expressed as

(1) $$s({t_k}) = {\sigma _k}\exp \left( {j2\pi {f_0}\left( {t - \displaystyle{{2R({t_k})} \over c}} \right)} \right),$$

where t is the full time, t _k is the slow time, σ _k indicates the backscattering coefficient, f ₀ is the carrier frequency (λ = c/f ₀) and c expresses the speed of light. R(t _k ) denotes the instantaneous distance between the scatterer and the radar. In the imaging plane (x, y), y corresponds to the range coordinate and x to the cross-range coordinate. The radar is laid far from the object in negative direction of the y-axis and the distance R(t) is approximately expressed as R(t) ≈ R ₀ + lsinθ ₀ cosωt + lcosθ ₀ sinωt ≈ R ₀  + x ₀sinωt + y ₀cosωt = R ₀ + d, where l is the distance from the scatterer to the center of the imaging plane (x, y), ω is a rotation rate, θ ₀ is the initial angle, and x ₀, y ₀ are initial coordinates at the time t _k  = 0. The distance d is the new instantaneous y coordinate, because time passes. The angle is a function of time θ _k  = ωt _k .

After demodulation, down-conversion to the baseband and removing the constant phase term, the returned signal becomes

(2) $$s({t_k}) = \sigma \exp \left( { - j4\pi \displaystyle{{d({t_k})} \over \lambda}} \right).$$

In the presented paper some scaling is assumed for an easier derivation of the image algorithm from the mathematical point of view. It is assumed that the backscattering coefficient σ _k  = 1 and 4π/λ = 1 for the sake of mathematical clarity. The parameter sampling rate PRF ≈ 1 Hz, rotation at the rate RPM ≈ 0.0028 are next examples of such simplifications. It is a non-realistic situation but σ and λ are only constant scaling numbers in the equation. This simplification does not restrict the generality of the considerations. Taking simplification into account, the proposed imaging algorithm can be put in an advanced signal processing area. Practical application requires knowledge about real radar system parameters: finite PRF, the CPI used for calculation of a projection, as well as the rotation speed ω, the cross-range dimension of the object, and so on. Finally after simplifications the baseband signal takes the form

(3) $$s({t_k}) = \exp ( - jd({t_k})) = \exp ( - j({x_0}\sin {\theta _k} + {y_0}\cos {\theta _k})).$$

For totally K scatterers on the object the baseband echoes can be regarded as a summation of sub-echoes from independent scatterers.

(4) $$\eqalign{s({t_k}) = & \sum\limits_l {{\sigma _l}\exp \left( { - \displaystyle{{j4\pi} \over \lambda} ({x_{0l}}\sin {\theta _k} + {y_{0l}}\cos {\theta _k})} \right)} \cr = & \sum\limits_l {\exp \left( { - j({x_{0l}}\sin {\theta _k} + {y_{0l}}\cos {\theta _k})} \right)} \cr & \quad under\;conditions: l \in [1,K], {\sigma _l} = 1, 4\pi /\lambda = 1.} $$

III. DOPPLER PROFILES AS TIME-FREQUENCY PROJECTIONS

The synthetic signal derived in equations (3) and (4) is used to generate projections. The rotating motion forms the sine shape of time-frequency spectrum. As a rigid object rotates about an origin such as that depicted in Fig. 1, all points with the same cross-range distance (in the x dimension) from the origin will have the same Doppler shift. According to the above discussion the Doppler projection (the cross-range projection) can be represented as the magnitude of the time-varying spectrum. The Short Time Fourier Transform (STFT) S(t _k , f) is a good time-frequency representation for non-stationary spectrum of a rotating object but has undesired property related with trade-off between the time and frequency resolution [Reference Tran and Melino4, Reference Auger, Flandrin, Gonçalvès and Lemoine8].

(5) $$\eqalign{S({t_k},f;h) = & \int_{{t_k}}^{{t_k} + {T_{CPI}}} {h(t - {t_k}} )s(t,{t_k}){e^{ - j2\pi ft}}dt,\,\,\,f = \displaystyle{{2\omega \cos {\theta _k}} \over \lambda} r \cr {f_x} = &\,\, {f_r}\cos \theta \;\quad {f_y} = {f_r}\sin \theta \;\quad {f_r} = 2f/c.} $$

The signal s(t, t _k ) is obtained during one CPI with the corresponding angle θ _k . The period time t _k  < t < t _k  + T _CPI is the time for one CPI, which is also the time duration of the analysis window h(t). Frequencies f _x , f_y are spatial frequencies in directions x and y, respectively, and f _r is the spatial frequency along the line r at the angle θ. As the angle θ changes, the slope of the line r is also changed. The magnitude of the STFT of a sufficiently long segment of the signal s(t, t _k ) can be interpreted as a Doppler projection signed ${p_{{\theta _k}}}(r)$ and shown in equation (6).

(6) $${p_\theta} _{_k} (r) = \left \vert {S({t_k},f;h)} \right \vert. $$

A good frequency resolution requires a long analysis window h(t). On the other hand the length of a window h(t) is limited by the rotational motion with the rotation rate ω of a scatterer. The window length should be sufficiently short in order to keep the scatterer in the same cross-range bin during one CPI, otherwise the scatterer will move to a different bin causing a smeared image. The spectrogram SP(t _k , f, h) = |S(t _k , f, h)|² is the squared magnitude of the STFT. We can interpret the spectrogram as a measure of the energy of the signal contained in the time-frequency domain centerd on the point (t _k , f). The spectrogram being the squared magnitude of the STFT is limited exactly as it is for the STFT. The resolution of the final image strongly depends on the resolution of the projection, so the new time-frequency representation with better resolution is strongly required. The reassigned spectrogram signed SPR(t′, f′, h) is proposed for building tomographic projections. Each value of the spectrogram SP(t _k , f, h) of the signal s(t, t _k ) computed at any point (t, f) is moved to another point $(\hat t,\hat f)$ , which is the center of gravity (or the centroid) of the signal energy distribution around (t, f). This is precisely the essence of the reassignment technique [Reference Auger, Flandrin, Gonçalvès and Lemoine8].

(7) $$\eqalign{& SPR(t^{\prime},f^{\prime},h) \cr & \quad = \int_{ - \infty} ^{ + \infty} {\int_{ - \infty} ^{ + \infty} {SP}} (t,f;h)\delta (t^{\prime} - \hat t(s;t,f))\delta (\,f^{\prime} - \hat f(s;t,f))dtdf.} $$

The proposed new formula of projections ${p_{{\theta _k}}}(r)$ is expressed by the following equation

(8) $${p_\theta} _{_k} (r) = SPR(t^{\prime},f^{\prime},h).$$

The better resolution of projections, the clearer image of the object can be obtained and closely spaced scatterers can be resolved. The image in the Cartesian plane (x, y) is attained via the filtered backprojection algorithm, which is an inversion formula for the inverse Radon transform. The filtered backprojection algorithm in the polar formatting is described by the equation

(9) $$\eqalign{\sigma (x,y) & = \int_{-\pi}^{\pi} \int_{-\infty}^{+\infty}P_{\theta} (f_r)e^{j2{\pi}f_{r}(x\cos \theta + y\sin \theta )} f_{r}df_{r}d\theta \cr & = \int_{-\pi}^{\pi} \left[{\int_{-\infty}^{+\infty}P_{\theta} (f_{r})e^{j2{\pi}f_{r}(x\cos \theta + y\sin \theta)}\left\vert{f_{r}}\right\vert df_{r}}\right] d\theta \cr & = \int_{-\pi}^{\pi}F^{-1} \{ \left\vert{f_{r}}\right\vert P_{\theta} (f_{r})\}_{(x\cos \theta + y\sin \theta )}d\theta,}$$

where σ(x, y) is a recovered reflectivity function, |f _r | is a ramp filter, F ⁻¹ is the inverse Fourier transform. The term P _θ (f _r ) stands for the spatial Fourier transform of the projection ${p_{{\theta _k}}}(r)$ for the angle θ along the radial direction r.

(10) $${P_\theta} ({f_r}) = \int_{r = - \infty} ^{ + \infty} {{p_\theta} (r){e^{ - j{\kern 1pt} 2{\kern 1pt} \pi {\kern 1pt} {f_r}{\kern 1pt} r}}} dr.$$

The inverse Fourier transform F ⁻¹ is calculated exactly at r = x cosθ + y sinθ [Reference Tran and Melino4].

Fig. 1. The geometry of the rotating point target.

IV. LIMITATIONS OF IMAGE RESOLUTION INTRODUCED BY THE POINT SPREAD FUNCTION (PSF), THE REASSIGNMENT METHOD AND RADAR SYSTEM PARAMETERS

The image resolution strongly depends on the accuracy of computing the variable expressed by equation (8). Integration limits and the frequency extent are finite in the reality. Variables f _r and θ belong to a limited set D.

(11) $$M({f_r},\theta ) = \left( {\matrix{ {1\;\quad {f_r},\theta \in \;D} \cr {0\quad \;{f_r},\theta \notin D} \cr}} \right).$$

The 2D inverse Fourier transform of the set M(f _r , θ) forms the PSF, which determines image resolution. The PSF should be like a delta function because this would imply that the recovered image would be identical to the reflectivity function σ(x, y). But this condition is never fulfilled because of limitation of the set D. A derivation of an analytic form of the PSF is not easy in general case, especially for an arbitrary-angle aperture. Only for the small-angle approximation i.e. for cosθ ≈ 1 and sinθ ≈ 0 and for the whole-angle aperture the compact relation for the PSF can be derived. The solution of the PSF takes the form of the sinc function in the range and the cross-range direction for a small aperture or the Bessel function of the first kind for the whole-angle aperture. In these cases the PSF is seen as the central symmetric function, where the resolution as null-to-null main lobe width is easily determined [Reference Sego, Griffiths and Wicks2]. The reassigned spectrogram with squeezing operation refocuses the time-frequency characteristics and is a strong nonlinear operation performed on the received echo. In consequence only simulations of one scatterer allow to examine the PSF. Limitation of imaging also results from limitation of the reassignment method. Due to the local action of the reassignment, the perfect localization property that holds for linear FM signals carries over to locally quasilinear situations. Radar echo from a rotating point is a signal with the sinusoidal FM modulation, so the localization is slightly degraded at those points where the instantaneous frequency trajectory cannot be treated as quasilinear (max and min points of the sinusoidal function). When the two components are sufficiently close, the spectrogram and the reassigned spectrogram are roughly equal to what we can obtain for a single component and the reassignment process may be imprecise. The readability of an image can be hampered by noise added to the signal. The reassignment treats components of a pure signal and noise equally and moves each value of the noisy spectrogram computed at any point (t, f) to another point (t ₁, f ₁), which is the center of gravity of the noisy signal energy distribution around (t, f), so that denoising procedure is required. Another limitation is related to radar system parameters, which should be properly established before applying the imaging process. The PRF needs to be sufficient to cover the entire Doppler frequency extent of the object expressed by the relation B _d  = 4rω/λ ≤ PRF in order to avoid aliasing. The expression of the full Doppler bandwidth B _d results directly from equation (5). For the established PRF, ω and λ, we can determine the maximal radial distance r of a scatterer without ambiguity (aliasing) of Doppler frequencies. Because the Doppler bandwidth B _d is a function of r, points with different radial distance r to the center of rotation will have different resolution of both projections and the PSF. The CPI should ensure no cross-range “walk-off” during the CPI. Otherwise degradation of the performance of the proposed methods is inevitable. In simulations the image is represented by the image matrix with dimension 360 × 360 and the rectangular coordinate system with the origin in the middle of the matrix (180, 180) is placed. Peaks (maximal value) after imaging correspond to the location of scatterers. Due to the PSF, peaks have non-zero vicinity characterized by a main lobe of the PSF. The PSF for the proposed imaging method possesses directional properties on the image plane because of strong nonlinearities. The resolution performance should be analyzed in two dimensions instead of the traditional 1D measure represented by the main lobe width.

A) Simulations

A helicopter is an example of a radar object with different rotating structures such as a rotor blade, a rotor hub or a tail rotor. After separation of a rotating part from the object body, only rotating structures on the target may induce additional periodical time-varying frequency modulations in the returned radar signal. An object of interest is rotor blades. Because of the shape of the rotor blade tips, returns from blade tips can be approximately treated as returns from rotating point scatterers. For two rotor blades only two rotating points symmetrically located, representing blade tips are interesting to build the image. In the well-focused image, the existence of two blade tips can be found, what enables to recognize some properties of the analyzed helicopter. Using an analogy to two rotor blades, an object simulated in experiments comprises two rotating points. In the first experiment two rotating points have been located at coordinates (30, 30) and at (−30, −30) in the imaging plane, symmetrically in respect to the center of rotation. The returned signal has been modeled according to equation (4) for two scatterers. In the proposed method of the imaging the reassigned spectrogram is converted to the binary matrix resulting in a gray-scale image. In spite of small inaccuracies in the vicinity of crossing of two frequency trajectories visible in Fig. 2, the proper locations have been obtained. The local maximum corresponds to the recovered point position.

Fig. 2. (a) The spectrogram of scatterers located at (30, 30) and (−30, −30), (b) the reassigned spectrogram of scatterers located at (30, 30) and (−30, −30), and (c) the tomographic image (zoomed) based on the reassigned spectrogram of scatterers located at (30, 30) and (−30, −30).

In other experiments (shown in Figs 3 and 4) closely and far located two scatterers, not necessarily symmetrically with respect to the center of rotating, have been analyzed. Accurate positions of another two scattering points relatively closely located each other at (10, 10) and at (30, 30) in the image plain and at (190, 190) and at (210, 210) in the image matrix, respectively, are easily recovered.

Fig. 3. (a) The spectrogram of scatterers at (10, 10) and (30, 30), (b) the reassigned spectrogram of scatterers at (10, 10) and (30, 30), and (c) the tomographic image (zoomed) from the reassigned spectrogram of scatterers at (10, 10) and (30, 30) with enhanced directional properties of the PSF.

Fig. 4. (a) The spectrogram of scatterers located at (10, 10) and (70, 70), (b) the reassigned spectrogram of scatterers located at (10, 10) and (70, 70), and (c) the tomographic image (zoomed) from the reassigned spectrogram of scatterers located at (10, 10) and (70, 70).

Similarly, for two scatterers located at coordinates (10, 10) and (70, 70) with relatively long distance between points, imaging results also turned out to be correct.

The experiment presented in Fig. 4 shows the large spatially variant frequency resolution of the spectrogram of two scatterers located in different distance from the rotating center, whereas the reassigned spectrogram mitigates effects of this phenomenon. The full Doppler bandwidth B _d  = 4rω/λ is strongly dependent on the rotation radius r. This causes the spatially variant frequency resolution of the STFT as well as the spectrogram for two rotating points with the different rotation radius. Despite squeezing properties of the reassigned spectrogram, an image of the scatterer with the larger rotation radius can create slightly larger spread than the point with a small radius (shown Fig. 4). The improvement given by the reassignment method is obvious but we never get a sharp curve of the ideal time-frequency distribution. Different resolution of an image of two points at (10, 10) as well at (70, 70) has been evaluated by the PSF. These locations correspond to points (190, 190) and (250, 250) in the imaging matrix, respectively.

Both the mainlobe and sidelobes turned out to be asymmetric, irregular with irregular tapering. This irregularity is a result of squeezing of the frequency spread by the reassignment method. As expected, the PSF is slightly dependent on the scatterer location and has the sharper peak for the closer located scatterer. The slice of the PSF through the column with the maximal peak shows the achievable resolution in this direction. In the experiment shown in Fig. 5 the first peak of the sidelobe has been below −22 dB relative to the maximal peak for the scatterer at (70, 70) and below −30 dB for the scatterer at (10, 10). Similar properties are attained for the slice through the row with the maximal peak. The reassignment method cannot resolve two too closely located components and cannot be reliable when Doppler ambiguity appear, what is depicted in Fig. 6.

Fig. 5. Slice through the column with the maximal peak in the imaging matrix (a) of the normalized PSF for the point at (10, 10), (b) of the normalized magnitude of the PSF (in dB) for the point at (10, 10), and (c) of the normalized magnitude of the PSF (in dB) for the point at (70, 70).

Fig. 6. The reassigned spectrogram (a) of two points located too closely: at (10, 10) and (15, 15), (b) of two points located too far: at (10, 10) and (150, 150), in experiments with ambiguity of Doppler frequencies, and (c) of two points located properly for the effective imaging process.

Experiments with simulated data confirm the effectiveness of the described idea of imaging of rotating points for assumed parameter of simulations.

V. TOMOGRAPHIC IMAGING IN PRESENCE OF NOISE

The model of the noisy signal can be written as x(n) = s(n) + w(n), where w(n) is circular complex Gaussian noise (the real and imaginary part of noise are independent and have the same probability distribution) with the variance σ _w ². A deterministic, non-stationary signal s(n) and noise w(n) are independent. The reassignment procedure used for noisy signals can create the false energy concentration about the center of gravity, because it uses both desired signal samples and undesired samples of noise. Sophisticated methods can be performed for extracting a signal from noise. One of many examples is using the time-frequency map to select the time-frequency area that is believed to belong to the signal, and thus to reject all other areas. In this paper a denoising procedure is performed by the rejection of values of the STFT below a selected level called a threshold and signed Th. Due to the linearity of the STFT, the transformation of noise and the transformation of signals can be considered separately. Random noise tends to spread its energy over the entire time-frequency domain, while a signal concentrates its energy within limited time intervals and frequency bands. To reduce the redundancy of the continuous STFT, we can sample it in the time-frequency plane on the rectangular grid {nt ₀, mv ₀}.

(12) $$\eqalign{& {F_x}[n,m,h] = {F_x}(n{t_0},m{v_0},h) \cr & \quad = \int_{ - \infty} ^{ + \infty} {s(u)h*(u - n{t_0})\exp ( - 2j\pi {\kern 1pt} m{\kern 1pt} {v_0}u)du},} $$

where s(u) is a signal, h(u) is an analysis window of the STFT transformation. The coefficients F _x [n, m, h] are coefficients of the STFT assigned to the signal s(u) (shortly signal coefficients). For the noisy signal x(t) two kinds of coefficients are obtained after the STFT: noise coefficients and signal coefficients. We should distinguish those coefficients that belong to the signal from the ones that belong to noise. It is natural to assume that coefficients assigned to noise are smaller than signal coefficients. By rejecting points smaller than an appropriate threshold, noise coefficients can be attenuated and signal coefficients can easily emerge. Thresholding process is very sensitive to the choice of the threshold. If the threshold Th is too large, the signal will be highly distorted, if the threshold Th is too small, then the result will still be noisy. After a denoising procedure the signal can be extracted by taking the inverse of the STFT using only extracted coefficients above the threshold. Theoretically STFT noise-only coefficients are distributed as a circular complex Gaussian variable, so the threshold above the variance value σ _w ² is a good choice [Reference Swiercz7].

A) Simulations in noisy environment

An iterative approach in the scheme STFT ≧ thresholding process ≧ inverse STFT is proposed to improve denoising results. The iterative denoising procedure results from the specific properties of 2D filtering in the time-frequency domain. The Fourier transform has a one-to-one mapping between the frequency domain and the time domain, so any spectrum in the frequency domain corresponds to a unique signal in the time domain. For time-frequency transforms, there is no guarantee that such a one-to-one mapping exists between the time-frequency domain and time domain. Setting up the threshold for STFT coefficients and taking ones above the threshold we make nonlinear 2D filtering by the filter localized in the region of occurrence of STFT signal coefficients, which are higher than the threshold. We should find a time-domain signal that has corresponding time-frequency characteristics in the localized region. There are two approaches for solving this problem. The first solution uses a minimization error between the time-frequency transform of the signal and the desired distribution after filtering. The second solution uses iterations, which allows to get easier the desirable signal without noise. The threshold is dependent on the variance of noise. After analysis of noise the threshold is set up at Th = 1.7*σ _w ² [Reference Swiercz7]. The accuracy of the denoising procedure is evaluated by the error between a pure signal and a synthesized signal after the iteration procedure. Noisy simulations characterized by a signal to noise ratio (SNR) are performed on the example of two scattering points located at (10, 10) and (70, 70) in the imaging plane. The SNR is calculated as the ratio of a sum of squared signal samples to a sum of squared noise samples. For the SNR = 4.1536 dB the best recovering of the noiseless signal is attained for 17 iterations. Denoising is the necessary preprocessing step in a tomographic imaging based on the reassigned spectrogram (Fig. 7).

Fig. 7. (a) The modulus of the STFT for the signal with SNR = 4.1536 dB, (b) the reassigned spectrogram after 17 recursions, and (c) the denoising result after 17 recursion.

Final imaging results of the tomographic imaging for the SNR = 4.1536 dB, the SNR = 1.6548 dB, and the SNR =−0.2834 dB are presented in Fig. 8.

Fig. 8. Imaging results of the tomographic imaging (a) for the SNR = 4.1536 dB, (b) for the SNR = 1.6548 dB, and (c) for the SNR = −0.2834 dB.

In Fig. 8 two peaks from two scattering points are visible separately. The tomographic image of the point located at a longer distance from the rotating center is seen more smeared and has the increasingly smaller value for the lower SNR ratio. The second peak is located at the proper cross-range bin for the SNR = 4.1536 dB and the SNR = 1.6548 dB. For the SNR = −0.2838 dB the image of the second point is moved to the next cross-range bin. The degradation quality of the image, especially of the second point, is the result of both the slightly worse reassignment process for the point with longer radius of rotation and of worse efficiency of the denoising procedure.