Open access peer-reviewed chapter

Adaptive Clutter Cancellation Techniques for Passive Radars

Written By

Tamás Pető and Rudolf Seller

Submitted: 24 May 2017 Reviewed: 26 September 2017 Published: 20 December 2017

DOI: 10.5772/intechopen.71289

From the Edited Volume

Topics in Radar Signal Processing

Edited by Graham Weinberg

Chapter metrics overview

1,936 Chapter Downloads

View Full Metrics


In radar systems, the ambiguity function of the applied illuminator signal essentially determines the detection capabilities. Zero Doppler interference (ZDI) or close targets returns can mask weak target reflections from higher distances. This is particularly the case for passive radars where the illuminator signal is not under the control of the radar designer. In recent times, great efforts have been carried out to research and develop efficiently working filter algorithms. These adaptive algorithms aim to cancel the undesired interference components in order to enhance the useful dynamic range. A number of different algorithms are operating in the space and also in the time domain. Spatial algorithms apply beamforming techniques, while temporal algorithms utilize the available reference signal to suppress the interferences. The main goal of this chapter is to present and compare the available spatial and temporal adaptive interference cancellation techniques in terms of filtering efficiency and computation cost on real-life data.


  • passive radar
  • DPIS
  • ZDI suppression
  • clutter cancellation
  • adaptive algorithm

1. Introduction

Passive radar systems utilize the so-called illuminator of opportunity (IO) for target detection. These transmissions are inherently not designed for radar purposes. The direct consequence of the inappropriate illumination is the presence of ambiguities and high sidelobes in the range-Doppler map. In addition, the received surveillance signal is often dominated by the direct and multipath Zero Doppler interference (ZDI) signals from the illuminator. The sidelobes of these unwanted high power clutter components can mask the weak target echoes [1]. The main goal of the adaptive clutter cancellation algorithms in passive radar applications is to get rid of the masking effect.

In recent times a number of research and development works aim to find the most adequate filtering technique. The most important aspects are the filtering performance, the capability for real-time operation and the purity of the resultant range-Doppler map. Filtering is performed typically in two different domain. Space-domain suppression techniques apply beamforming methods to deal with the high power interferences. The adaptive optimum beamformer for passive radar application has been investigated in [2, 3]. Methods which utilize the eigen-structure of the spatial correlation matrix have been proposed in [4, 5]. Moscardini et al. [6] have made efforts to identify the most suitable multi-channel processing architecture. Their work focuses on the optimum beamformer for digital video broadcasting-terrestrial (DVB-T)-based systems. Villano et al. examined the non-adaptive beam steering and null-steering methods as well as the adaptive correlation matrix inversion-based methods [5].

Beside the application of spatial filtering by means of digital beamforming, the temporal algorithms also has great importance. These filtering algorithms utilize the available reference signal to cancel out the time delayed and Doppler shifted replicas of the illuminator signal on the surveillance channel. When it comes to time domain filtering, the Wiener filtering problem is solved without exception. Among the first, Cardinali et al. have investigated the application of the least mean square (LMS), normalized least mean square (NLMS), recursive least squares (RLS), extensive cancellation algorithm (ECA) and sequential cancellation algorithm (SCA) in passive radar scenario [7]. In their work they used simulated data with assuming FM transmitter. A similar study with more extensive examinations on real life data can be found in [8]. Garry et al. [8] examined the ECA, CLEAN, NLMS, fast block LMS (FBLMS) and RLS algorithms using DVB-T IO.

The main objective of this chapter is to evolve a common basis of discussion for the clutter cancellation algorithms. In this present study, the ECA, ECA-B, ECA-S, LMS, NLMS, RLS and block LMS methods are investigated. The supplemented autocorrelation matrix method for the block type algorithms and the weight inheritance technique for the iterative algorithm are also presented. Beside the evaluation of the clutter suppression capabilities, the survey of the required resources has fundamental importance. The assessment of the computational costs provides support for the selection of the adequate algorithm for every specific scenario.


2. Signal model

Before we start to discuss the operation of different algorithms, the proper signal model is constructed at the first place.

Figure 1 illustrates the considered signal model. The signal transmitted from the illuminators is received on a number of different signal paths. We can separate the individual signal contributions on the surveillance channel to three different types. These are the reference signal correlated interferences, the useful target reflections and the uncorrelated receiver noise. The signal processing steps are implemented digitally in all real cases, thus we use a discrete time domain description. Let us denote the illuminator signal with s[n], where n = 0…N is a discrete time index, and the signal is samples with fs sampling frequency. The number of samples in coherent processing interval (CPI) is denoted by N. Eq. (1) describes the signal received on the passive radar antenna array for one target reflection.

x a n = k = 0 K 1 α k s n k e j 2 π f k f s n a θ k + α t s n t e j 2 π f t f s n a θ t + ξ n E1

Figure 1.

Signal model for the passive radar scenario.

We are taking the following assumptions. The direct and multipath clutter components are received on the first K range cell. Each signal component has a propagation factor, which is described by the α parameter. The fk, k = 0…K − 1 parameters describe the Doppler frequency of the considered clutter component. For static reflections, fk = 0 and they are also called as the ZDI components. a(θ) denotes the spatial signature of the current signal component, that is arriving to the antenna array from the θ angle. Correspondingly, the target reflection has the αt, ft, θt parameters. Finally, the uncorrelated thermal noise vector collected by the receiver array is denoted with ξ[n]. The received signal array xa[n] is an M dimension array, where M is the number of antenna elements in the receiver array. For an M element equidistant linear antenna array the array response vector has the following form:

a θ = 1 e jβdcos θ e jβdcos θ 2 e jβdcos θ M 1 T C Mx 1 , E2

where θ is the incident angle of the impinging signal and β is the wavenumber. The detection stage is often preceded by the channel preparation algorithms, which aims to produce the clean reference and surveillance channels. On the reference channel, it is desired to receive only the direct path signal. Thus, in ideal case the reference signal has the following form:

x r n = s n + μ n , E3

where μ[n] is the thermal noise. This can be often guaranteed by the use of high-gain antennas or the application of reference signal reconstruction algorithms [9]. Signal reconstruction algorithms are mainly used for illuminators, that use digital modulation techniques such as the digital audio broadcasting (DAB), DVB-T, digital video broadcasting: satellite (DVB-S), universal mobile telecommunications system (UMTS), long-term evolution (LTE) signals.

On the other hand, the surveillance channel must contain only the signal component reflected from the moving target in order to prevent performance degradation at the detection stage. It can be proved by analytical investigations that the detection performance is affected by both the reference and surveillance signal purity [10]. Assuming white noise illuminator, the signal-to-interference plus noise ratio (SINR) of the target can be described as the function of the SINR of the reference and surveillance signals.

SINR target = 1 + SINR ref + SINR surv + N SINR ref SINR surv 1 + SINR ref + SINR surv + SINR ref SINR surv E4

The dependencies of the observed target SINR on the range-Doppler surface can be calculated according to Eq. (4) and illustrated in Figure 2. In Figure 2, the number of samples in a CPI is N = 220. In this present chapter, we focus our attention to the filtering and preparation of the surveillance channel.

Figure 2.

Target SINR (dB) in the function of the reference and surveillance signal purity.


3. Algorithm evaluation methodology

In this section, the way of comparing the performances of the inspected algorithms is described in detail. The performance of the different clutter cancellation techniques is evaluated on real-life measured data.

3.1. Reference receiver system

The receiver system used to obtain the raw data is a quad channel software defined receiver especially designed for passive radar application. The system receives the surveillance signals on a quad channel fractal patch antenna system. On this surveillance antenna array, the elements are placed equidistantly with having 0.628λ spacing between each other, where λ denotes the wavelength of the center frequency. Among the four antenna elements, only three are connected to the receiver and the remaining one receiver channel is dedicated to the reference channel. This channel is connected to a Yagi antenna that receives the direct path signal. The detailed description and the design considerations of the software defined receiver can be found in [11]. During the measurements, the receiver system was only able to record the received multichannel signal in data chunks instead of streaming. This limitation results in short signal gaps between the consecutive CPIs. The CPI used in these investigations is TCPI = 0.126s. The duration of the signal gaps is dependent on the sample download speed, but its average is approximately Tgap ≈ 0.6s.

3.2. Measurement scenario

The measured data were recorded near an airport using DVB-T signal as the illumination. During measurement, there were three illuminators operating in single frequency network mode at 634 MHz, and with a useful bandwidth of 7.61 MHz. In this analysis, the data acquired on Dec 21, 2015 have been processed and examined. The recorded data contains observations from several landing airplanes. After performing the initial processing steps, a long range target have been selected as the reference target for algorithm performance evaluations. For the analysis 25 CPIs have been processed.

3.3. Performance metrics

The clutter cancellation performance of the investigated algorithms can be measured on a wide variety of metrics. Choosing the proper metric is essential to obtain an objective comparison. Cardinali et al. [7] introduced the clutter attenuation (CA), which measures the averaged power ratio of the original and the filtered surveillance signal.

CA = n = 0 N 1 x f n 2 n = 0 N 1 x s n 2 E5

In some cases, the applied surveillance channel filter unintentionally suppresses the target reflection beside the clutter. As the CA metric does distinguish the target reflection and the clutter higher CA value may represent worst filter performance. Also, Garry et al. [8] pointed out that CA is not reliably for the estimation of the RD map noise floor reduction. Instead of the CA, they proposed to measure the ratio of the estimated noise floor reduction in the RD map. To obtain this metric, they designated certain area on the RD map where no target reflection and dominant clutter components are expected. The noise floor reduction can be estimated over the specified region using the following equation.

R NF = τ T , f F χ τ f x f 2 τ T , f F χ τ f x s 2 , E6

where T and F defines the sets of the considered time delays and Doppler frequencies of the region of interest. The range-Doppler cell having τ time delay and f Doppler shift is calculated as follows:

χ τ f x = n = 0 N 1 x n x r n e j 2 π f f s n E7

In contrast to the CA, the noise floor reduction metric isolates the clutter and the target echo contributions however, it still does not prove information about the impact of the filter on the useful target reflection. To overcome this limitation in [3, 8, 12], the SINR improvement of a detected target has been used as a metric. For the estimation of the target SINR we have

SINR target x = χ τ 0 , f 0 , x 2 1 Q τ T , f F χ τ f x 2 , E8

where τ0 and f0 are the time delay and the Doppler shift of the detected target correspondingly. T and F sets contains the time delay and Doppler shift values of the neighboring range-Doppler cells, and Q denotes the number of cells used to estimate the clutter power. The improvement of the target SINR is defined as the ratio of the estimated values between the filtered and the original surveillance signals.

R SINR = SINR target x f SINR target x s E9

The target SINR improvement can be used effectively to evaluate the operation of different algorithms, however, its limitation must also be taken into account. The SINR of a detected target is dependent both on the purity of the reference and the surveillance signals. In case the reference signal has low quality, the achievable target SINR is limited, regardless of the interference suppression performance on the surveillance channel. Beside this, the relation between the surveillance signal purity and the target SINR is not linear, thus even a dramatical increase in the clutter suppression may result in slight target SINR improvement. Both of the highlighted limitations can be observed in Figure 2.

In the following presented analysis, the target SINR improvement is displayed as the main performance metric. According to preliminary results, we can declare that the observed and estimated values have large variation along the evaluated CPIs. To exclude this effect the performance metrics are evaluated on several consecutive CPIs along the target trajectory and the mean value of the metrics are indicated and compared.

3.4. Computation cost analysis

The computation cost of the examined algorithms is determined by specifying the necessary floating point operations (FLOPs) to calculate the filtered output signal. The cost of a complex addition and a multiplication is taken into consideration as 2 and 6 FLOPs, respectively. Several algorithms perform matrix inversion to determine the coefficients used for the filtering. To estimate its requirements, 8J3 FLOPs are accounted for the inversion of a J dimensional matrix. The computational requirements are analyzed for each of the algorithms with considering their possible parameter settings.


4. Space domain filtering

During the process of space domain filtering, the M antenna channels are combined together. The coefficient vector of the beamformer (wsd ∈  C MxL ) is calculated in a way to add the interference component out of phase while summarize target energy in-phase from the separate antenna channels. The beam-space processed output signal is calculated using Eq. (10).

y sd n = w sd H x a n E10

The algorithms operating in the space-domain can be either data dependent adaptive algorithms or non-adaptive fixed beamformers. The main advantage of data independent methods is the fast output calculation. In this section, we only deal with beamformers that do not alter the coefficient vector during CPI.

4.1. Maximum signal-to-interference ratio

The maximum signal-to-interference ratio (SIR) is a data independent method, where the beampattern is manually calculated specifying the necessary constraints. These constraints are the incident angles and the desired responses of the antenna system in the specified directions. One can calculate the corresponding coefficients that fulfill the constraints using Eq. (11).

w sd fix H = u T A H A A H + σ n 2 I 1 , E11

where u is an L dimensional constraint vector and A is the array response matrix created from the specified directions, I is the identity matrix and σ n 2 is the power of added noise that is used to prevent instability [13].

A = a θ 0 a θ L 1 C MxL E12

Note that, the degrees of freedom of the beamformer is equal to the number of elements in the antenna system, thus only L constraints can be set (L ≤ M) For passive radar application, the constraint vector and the array response matrix can be set as follows:

u = 1 0 0 E13
A = a θ t a θ 1 a θ L 1 E14

The target is expected at θt angle and θl, l = 1…L − 1 are the incident angles of the clutter components.

4.1.1. Fixed max SIR

In this case the beampattern is not changing, thus the cofficient vector is set at the beginning of the processing and it used for all CPIs. In this case, the θl interference angles can be determined and set from preliminary information about the positions of the radar and the illuminators. This solution suffers from the inaccurate incident angle estimations and it is not able to deal with interference sources that change their positions over the surveillance time.

4.1.2. DOA estimation supported max SIR

Another way to configure the A array response matrix is to estimate the DOA of the dominant signal components in the received signal vector xa[n]. Thus the beamformer gets the opportunity to deal with all interference sources and to adapt to the changing environment. Note that in practical cases, the target reflection has low power relative to the ZDI and the other interferences, thus it does not affect the DOA estimation of these components. The choice of the DOA estimation technique has great importance as the performance of the interference suppression depends on its accuracy.

4.2. Adaptive optimum

Different optimization criteria such as the minimum variance distortionless response (MVDR) or the maximum signal to interference plus noise ratio (MSNR) leads to the well-known optimum weight vector calculation method described by Eq. (15).

w sd opt = R s 1 a θ t E15

In Eq. (15), Rs ∈  C MxM denotes the spatial correlation matrix of the interference components. It is often approximated by the correlation matrix of the full received signal, which is defined as

R s = 1 N n = 0 N 1 x a n x a n H E16

To calculate the corresponding beampattern, the expected incident angle of the target and the spatial autocorrelation matrix is required. In case the DOA of the target is not known precisely, the target reflection may be suppressed unintentionally. That is the so called pointing error. Another group of problem arise from the fact that the spatial correlation matrix contains other signal components than the interferences. Moscardini et al. [6] made progress on this problem with estimating the spatial autocorrelation matrix of the interferences from the corresponding range-Doppler cells after performing the cross-correlation. They report 2dB improvement over the conventional methods.

4.3. Subspace-based technique

Eigen-subspace-based beamforming techniques for passive radar application is proposed in [4, 5]. According to these works, the spatial correlation matrix can be decomposed to the signal subspace and the noise subspace. In practical cases, the target echo having relatively low power resides in the noise subspace. Also we can assume that the unwanted disturbances will determine the dominant eigenvalues of the spatial correlation matrix. Then it is possible to suppress the interferences by projecting the array response vector of the target to the noise subspace.

w sd eig = Q Q H a θ t , E17

where the columns of the Q matrix are the eigenvectors assigned to the smallest eigenvalues of the spatial correlation matrix, that are the vectors of the noise subspace. The computation intensive eigenvalue calculation can be avoided by using the Power of R (POR) method proposed in [4]. Villano et al. [5] suggested to select the noise subspace as the subspace orthogonal to the eigenvectors that is assigned to the principal eigenvalues. In the current investigation, it is found that the best result is achieved by constraining the principal eigen-subspace dimension to one.

4.4. Comparison of the space-domain filters

Table 1 summarizes the obtained SINR improvements for the different beamformers. As it can be seen, the fixed maximum SIR beamformer has the lowest performance. This beamformer has no ability to react to the changed environment. When the information from DOA estimation is utilized continuously this beamformer can clearly achieve gain over the non-adaptive case. The adaptive optimum method has the highest average gain. The adaptive principal eigenvalue beamformer in most cases performs worse than the optimum beamformer. This is due to the fact that only one dominant clutter component is canceled and most of the clutter power remains unsuppressed in the beam-space processed signal (Figure 3).

Beamformer SINR gain (dB)
Fixed max SIR 3.9
DOA supported max SIR 4.5
Adaptive optimum 5.4
Principal eigenvalue 5.2

Table 1.

Obtained SINR improvements of the different beamformer algorithms.

Figure 3.

Obtained radiation patterns of the different beamformers.

In [2] Di Lallo et al. have measured 13 − 25dB suppression using two antennas and FM illumination. A wide-ranging study of the beamforming algorithms application for the Metropolitan Beacon System (MBS) illuminator can be found in [14]. Navrátil et al. investigated the effectiveness of the direct signal separation utilizing the elevation angle difference of the target and the illuminator. Using four antennas they achieved 10 dB improvement measured on the SINR of a detected target.


5. Time domain filtering

The adaptive temporal algorithms rely on the fact that the ZDI signal components in the surveillance channel can be reproduced and subtracted by the use of the reference channel. According to the considered signal model, a FIR filter is able prepare the properly weighted and time delayed replicas of the reference signal. The n-th sample of the time domain filtered signal is calculated with

y td n = w td H x r n , y td n = x s n w td H x r n , E18

where wtd ∈  C Jx1 is the coefficient vector of the time domain filter and J is the tap size of the filter. J is set to the range of the farthest expected clutter component. We can obtain the temporal coefficient vector in the minimum mean squared error (MMSE) sense by using the well-known solution of the Wiener-Hopf equations [15]. In matrix form we have

w td = R t 1 r t E19

Rt ∈  C JxJ denotes the so called temporal autocorrelation matrix and rt ∈  C Jx1 is the temporal cross correlation vector. These are defined as follows

R t = E { x r n x r n H } r t = E x r n x s n , E20

where xr[n] contains the last J samples of the reference signal and xs[n] is the n th sample of the surveillance signal.

x r n = x r n x r n 1 x r n J + 1 T E21

In practical cases, the temporal autocorrelation matrix and cross correlation vector is not known preliminary and thus it is inevitable to estimate them from the measured data. The algorithms discussed thus far can be partitioned into two main groups. The block type algorithms estimate and apply the coefficient vector of the filter on larger data chunks, while the iterative type algorithms update the used coefficient vector from sample to sample. In this chapter, the investigated and presented block type algorithms are the Wiener-SMI, ECA, ECA-B and ECA-S methods and the iterative methods are the LMS, NLMS, RLS and BLMS. In the presented investigations, the dimension of the temporal filters is uniformly set to J = 128 tap. This choice ensures that most of the clutter contributions are removed and also grants fast calculations to the analysis. Note that, the tap size of the temporal filter must be fitted to the environment in every specific case.

5.1. Wiener: sample matrix inversion (SMI)

The most-simple algorithm is the direct application of Wiener filtering technique with the sample average estimation of the temporal autocorrelation matrix and cross-correlation vectors.

R t = 1 N n = 0 N 1 x r n x r n H E22
r t = 1 N n = 0 N 1 x r n x s n E23

The filter has only one parameter, which is the tap size. This parameter controls the number of considered clutter components. Figure 4 illustrates the mean value of the obtained target SINRs for different filter tap sizes. As we increase the filter tap size we can achieve higher filtering performance.

Figure 4.

Wiener-SMI filtering performance with different tap-sizes.

After reaching a certain filter tap-size the rate of the achievable further improvements decreases. In real systems, the optimal choice of the filter dimension can be set according to such an analysis. Note that the curve of the achieved improvements in the function of the filter dimension is fundamentally determined by the clutter profile of the environment. An exhaustive study on the clutter profile for FM-based passive radar can be found in [16].

The Wiener filter implemented with the SMI technique leads to solution of the LS filter. The LS filter minimizes the squared error of the filtered output signal array.

min w td x s X w td 2 E24

The filtered surveillance signal that satisfies Eq. (24) can be calculated using Eq. (25)

x f = x s X w td = I N X SMI X SMI H X SMI 1 X SMI H x s , E25

where XSMI is the reference signal subspace matrix, that is composed of the time delayed replicas of the reference signal array.

X SMI = x r Dx r D 2 x r D J 1 x r C NxJ E26

In Eq. (26) xr denotes the vector of the reference signal samples from the whole CPI.

x r = x r 0 x r 1 x r N 1 T C Nx 1 E27

The D matrix performs the time delay by shifting the elements of the reference signal vector xr. It has a size N × N and the elements are defined by the following expression:

d i , j = 1 0 , i = j + 1 , otherwise E28

In its form of the filter it can be considered as a projection. The Wiener-SMI algorithm projects the surveillance signal array to a subspace orthogonal to the reference signal and its time delayed replicas [17].

The processing demands for the Wiener-SMI algorithm can be written as

C Wiener SMI = 8 J 3 + J 2 + J 2 N + 2 JN + 2 N E29

Figure 5 illustrates the computational costs for different filter dimensions. The algorithm is relatively computational intensive, however using the method described in section the requirements can be greatly reduced.

Figure 5.

Computation costs of the wiener-SMI algorithm for different filter dimensions.

5.2. Extensive cancellation algorithm (ECA)

The ECA algorithm is first introduced by Colone et al. in [17, 18].They realized that beside the ZDI a portion of the disturbances resides in the low Doppler frequency region. The main reason is that vegetation and slowly varying environment spreads the clutter energy around the zero Doppler line. To cope with this effect, they proposed to extend the reference signal subspace matrix to the Doppler frequencies.


Λ ∈  C NxN is the transformation matrix of the Doppler shift. In order to perform Doppler shift on the reference signal vector with fd frequency, the Λ matrix should have a diagonal form with the following values:

Λ = diag 1 , e j 2 π f D f s e j 2 π f D f s 2 ,… e j 2 π f D f s N 1 E31

The P parameter controls the filter notch width in the Doppler domain. The wider filter notch means more clutter power to cancel and thus it ensures better target detection capabilities. At the same time, the overextended filter notch may prevent the detection of low Doppler frequency targets. The RD map of the ECA filtered surveillance signal can be seen in Figure 14. The filter notch for P = 5 can be inspected in Figure 15. The computational burden can be approximated by replacing J with J(2P + 1) in Eq. (29). For reasonable P values this can lead to extremely high requirements. For this reason, this algorithm is not applicable in practical cases.

5.3. Batched ECA algorithm (ECA-B)

The Wiener SMI and the ECA algorithm use all the samples in a CPI to estimate and apply the coefficient vector of the filter. In non-stationer environment, the optimal coefficients may change rapidly over the CPI and these algorithms are not able adapt fast enough. This case the calculated coefficient vector may be smoothed between the instantaneous optimal values. To overcome this limitation Colone et al. proposed the batched version of the ECA algorithm [17]. This method divides the CPI into shorter blocks and performs the filtering individually on these consecutive signal portions. This makes the filter more robust against the time-varying characteristic of the environment. Let us partition the N sample of the CPI to T number of batches. We can then describe the t-th signal fraction of surveillance channel as follows:

x s t = x s t N T x s t N T + 1 x s t + 1 N T 1 E32

The t-th batch of the reference signal x r t is defined in a similar way. Correspondingly the signal subspace for the t-th batch X ECA B t is constructed from the Doppler shifted and time delayed replicas of the x r t signal vector. Following the same filtering procedure, we have Eq. (33) to calculate the t-th batch of the filtered surveillance signal.

x f t = P t x s t = I N X ECA B t X t ECA B H X ECA B t 1 X t ECA B H x s t E33

A direct consequence of the fractioned processing is the widening of the filter notch around zero Doppler. Colone et al. have shown by analytical investigation in [19] that the notch width of the ECA-B filter is inversely proportional with the batch size and has the shape of sin(πfNB/fs)/ sin(πf/fs) where NB is the number of samples in a batch NB = N/T. According to this result the reduction of the batch size increase the filter notch width and enhance the domain of the canceled clutter components. However, the decreased batch duration may worsen the filter adaptivity. In this case the number of samples used for the coefficient estimation also starts to become insufficiently low that yield inaccurate estimation see Figure 6. Another remarkable nature of this filter is the introduction of unwanted sidelobes in the Doppler dimension. These undesired Doppler structures can highly affect the operation of the detection stage and thus resolved in the later development of the algorithm.

Figure 6.

Achieved SINR improvement of the ECA-B algorithm with different batch durations.

Figure 14 illustrates the RD map obtained when the CPI is partitioned to 32 batches. As it is apparent the ECA-B filter can significantly reduce the slowly moving clutter with its extended filter notch. We can also observe the appearing Doppler sidelobes that complicate the identification of true targets. The computation need of the algorithm has the following form.

C ECA B = 8 J 3 T + J 2 T + J 2 N + 2 NJ + 2 N E34

The calculated values for different batch intervals and filter dimensions are illustrated in Figure 7. It is clear from the figure that for long CPIs the required FLOPs are practically identical regardless of the batch duration. When T = 1, the ECA-B algorithm reduces to the Wiener-SMI. This also means that their calculation requirements for long CPIs are practically the same. However, it is worthwhile to mention that the memory requirement of the ECA-B algorithms is less in the direct implementations.

Figure 7.

Computational costs of the ECA-B algorithm with different filter dimensions and batch duration.

5.4. Sliding window ECA algorithm (ECA-S)

The previously introduced ECA-B algorithm suffers from the effect of parasitic Doppler sidelobes arising from the fractionated filtering. The theoretical investigations in [19] carried out that these Doppler ambiguities are separated by 1/TB, where TB = NBfs is the batch duration. One can say that by decreasing the batch size the Doppler ambiguities can be moved out of the interested region of range-Doppler map. However, as described earlier the direct consequence of decreasing the batch duration is the widening of the filter notch around the zero Doppler frequency. To resolve this contradiction Colone et al. [19] suggested to apply different window sizes for the coefficient estimation and the filtering. The coefficient estimation window is selected symmetrically around the filtering window. The batch size used for the filtering is denoted by Nf while we introduce the NE parameter that describes the number of samples used for the coefficient estimation.

x s , est t = x s t N f N e N f 2 x s t N f N e N f 2 + 1 x s t + 1 N f N e N f 2 1 E35
x s , filt t = x s t N f x s t N f + 1 x s t + 1 N f 1 E36

After all, the separation between the Doppler ambiguities can be controlled by the Nf parameter while the filter notch can be configured by choosing the proper Ne parameter. This filtering technique also improves the filtering performance of ECA-B as the coefficients are calculated on overlapped signals fractions thus we get a smoothed estimate. Figure 14 shows the range-Doppler map of the ECA-S algorithm. It is clearly visible that the unwanted Doppler structures are totally disappeared while the width of the filter notch is remained the same. Note that the algorithm reduces to the ECA-B method when Nf = Ne. This method utilize more processing power than the ECA-B algorithm as the windows used for the coefficient estimations are inherently wider. The exact requirements can be written as

C ECA S = 8 J 3 T + J 2 T + J 2 T N a + JT N a + JN + 2 N E37

Figure 8 depicts the emerging extra calculation compared to the ECA-B algorithm. The different curves belong to different window sizes used for the coefficient estimation. The tap-size is set to J = 128, and the CPI is partitioned to 32 filtering windows. Note that the curve corresponding to Ne = CPI/32 is identical with the computational need of the ECA-B algorithm when T = 32, see Figure 7.

Figure 8.

ECA-S algorithm computation cost with different estimation window lengths. (Nf = N/32).

5.5. Supplemented autocorrelation matrix technique

As it is well known the autocorrelation matrix of a wide sense stationary (WSS) process has Toeplitz and Hermitian properties. Ri, j = Ri + 1, j + 1 = ri − j and R = RH. This allows us to calculate only J elements from the J × J-sized R autocorrelation matrix and supplement the remaining elements. The obvious advantage of this technique is the fast coefficient calculation. At the same time, the limited information about the random process may lead to inaccuracies that worse the filtering performance. In the current investigations, the calculated J elements are selected as the first column of the autocorrelation matrix. Figure 9 shows the comparison of the achieved result on the Wiener-SMI and the ECA-B algorithm. As we can observe the difference is negligible when the sample count used for the coefficient estimation is relatively high. However, in case of the ECA-B algorithm, when short batch duration is applied the sample count is very low and the estimated correlation coefficients have still large variance. For this reason, the performance noticeably decreases.

Figure 9.

SINR loss arising from the usage of the supplemented autocorrelation matrix in the ECA-B algorithm.

As a conclusion, we can say that for algorithms, which use large signal batches for the coefficient estimation this technique has definitely relevance. The achievable speed up is dependent on the filter depth thus on the dimension of the temporal autocorrelation matrix.

The direct form of the modified computational costs for the Wiener-SMI, ECA-B and ECA-S algorithms are summarized as follows.

C S Wiener SMI = 8 J 3 + J 2 + 3 JN + 2 N C S ECA B = 8 J 3 T + J 2 T + 3 JN + 2 N C S ECA S = 8 J 3 T + J 2 T + 2 JT N e + JN + 2 N E38

The achieved speed-ups for different processing batch sizes are illustrated in Figure 10. For large CPIs, there is no significant difference in the speed up ratio for different batch durations. For the ECA-S algorithm the speed-up ratio is always greater than the ECA-B method has because, the simplified J2 term has larger coefficient.

Figure 10.

Achievable speed up with the supplemented autocorrelation matrix technique.

5.6. Least mean square (LMS)

In contrast to the block type algorithms the LMS algorithm updates the wLMS coefficient vector every time when a new signal sample is digitalized. The filter belongs to the family of the stochastic gradient-based algorithms. That means the filter updates the coefficient vector from the current estimation of the error gradient. The update equation for the coefficient vector can be written as

w n + 1 = w n + μ x r n e n , E39

where e[n] denotes the filter error at the n th time instant, which is defined as the difference between the n th sample of the filter output and the surveillance signal.

e n = x s n w LMS n H x r n E40

The coefficient vector is often initialized with the zero vector at the beginning of the processing interval. The filter uses the μ step size parameter to control the influence of the current modification. Higher step size values offer faster reactions to the changes in the environment, however the filter will suffer from misadjustment and will not able to suppress the ZDI properly. At the same time choosing the step size too small, result in sluggish filter response. In order to ensure stability Eq. (41) must also be satisfied.

0 < μ < 2 j = 0 J 1 E x s n j 2 E41

The optimal choice of the step-size is also dependent on the dimension of the filter and the power of the reference signal. Thus, the proper parametrization of the filter can be a complicated task. This greatly reduces the applicability of the filter. The obvious advantage of the LMS algorithm is that it calculates the filtered output very fast, thus it is capable of the real-time operation. The estimated target SINRs for different step-size values are shown in Figure 11.

Figure 11.

SINR improvements of the LMS algorithm variants for different step-sizes.

5.7. Normalized least mean square (NLMS)

The NLMS algorithm makes the LMS algorithm independent from the scale of the reference signal vector with normalizing the step-size parameter. For the coefficient vector update we have

w NLMS n + 1 = w NLMS n + μ n a + x r n H x r n x r n e n , E42

where a is small number used to prevent instability and μn is the step-size of the NLMS algorithm, which must be chosen in the range of 0 < μn < 2. With normalizing the energy of the reference signal vector the algorithm gains stability and better convergence properties. The price paid for this improvement is the increased computational cost. In every iteration, it is now necessary to determine the instantaneous energy of the reference signal and its time delayed replicas. The convergence behavior of the NLMS algorithm for a variety of step-size values are illustrated in Figure 12.The displayed curves represent the evaluation of the L2 norm of the coefficient vector.

Figure 12.

Convergence behavior of the NLMS algorithm for different step-sizes.

5.8. Recursive least squares (RLS)

Beside the instantaneous error the RLS algorithm also takes into consideration the weighted sum of the previous error values.

e RLS n = i = 0 n λ n i x s i w RLS n H x r n 2 E43

where λ is a small number, which is often referred to as the forgetting factor. With choosing the forgetting factor in the range of 0 < λ ≤ 1 the algorithm reduces the influence of the previous error values over time. The RLS algorithm recursively finds the optimal coefficient vector that minimizes the weighted sum of the error squares. To accomplish this, the algorithm utilize the Sherman-Morrison formula to the calculate the inverse of the updated autocorrelation matrix. The coefficient update equations of the RLS filter can be written as follows.

k n = λ 1 P n 1 x r n 1 + λ 1 x r n H P n 1 x r n α n = x s n w RLS H n 1 x r n w RLS n = w RLS n 1 + k n α H n P n = λ 1 P n 1 λ 1 k n x r n H P n 1 E44

The starting values of the inverse of the autocorrelation matrix must be initialized for the RLS algorithm. It is often estimated on a small number of samples.

P 0 = i = 0 i 1 x r i x r i H 1 E45

The filter is far more complex than the LMS or the NLMS algorithm, thus the real-time operation is difficult to implement. However, the RLS algorithm has faster convergence than the LMS or NLMS algorithm.

5.9. Weight inheritance technique

The iterative algorithms always have to initialize the coefficient vector to start the iteration. In most applications, the zero vector is an adequate choice. Observing Figure 12 we can notice that the re-initialization of the coefficient vector at the beginning of the CPI greatly degrades the convergence and lengthens the transient sections. Bárcena-Humanes realized that the large transient sections are responsible for the smoothing effects in range-Doppler map. This effect can be observed in Figure 14 on the result of the LMS algorithm. To shorten or even eliminate the transient sections he proposed to initialize the coefficient vector of the current CPI with last value of the coefficient vector of the previous CPI.

w 0 p = w N 1 p 1 , E46

where p denotes the index of the currently processed CPI. The success of this modification depends on the stationarity of the environment and the elapsed time between the consecutively processed CPIs. In case the time gap is reasonably large the parameters of the environment may change to a great extent and the filter has to compensate the accumulated error.

Figure 13 shows the result on the LMS filter convergence with and without the weight inheritance technique. On the results obtained from the measured data we can see some improvements, however the time gap between the CPIs was too large. Thus, the necessarily evolving transient sections still greatly affects the range-Doppler map. As a result, no improvement is realized.

Figure 13.

Effect of the weight inheritance technique on the LMS algorithm.

5.10. Block least mean square (BLMS)

In the LMS algorithm the fast variation of the instantaneous error signal can lead to the poor estimation of the gradient vector. These disadvantages can be avoided with averaging the gradient vector over a block samples. The block LMS filter only updates the coefficient vector after averaging L gradient estimations. The coefficient vector is then updated at the t + 1-th block according to Eq. (47).

w BLMS t + 1 = w BLMS t + μ B 1 N B i = 0 N B 1 e t N B + i x r t N B + i E47

The coefficient vector is only updated N/NB times instead of N, where NB denotes the number of samples in a block. The block size of the filter is often set to be equal to the filter length NB = J. The μB parameter gives the step-size of the algorithm. The block LMS filter has a fast implementation in the frequency domain. The fast block LMS (FBLMS) algorithm use the overlap-save method to calculate the linear convolution. The fast block algorithm derived from the NLMS algorithm is first proposed for passive radar application by Zhao et al. [20]. They have shown that the fast block NLMS (FBNLMS) algorithm can realize improvement over the traditional NLMS algorithm both in terms of cancellation and computational cost.


6. Summary

The main objective of this chapter is to provide support on the proper selection of the clutter cancellation algorithms. In the previous sections many of the so far proposed algorithms were investigated, their advantages and weaknesses were analyzed. Figure 14 shows the calculated range-Doppler maps.

Figure 14.

Range-Doppler maps of the differently filtered surveillance channels.

In Figure 15 filter notch with can be inspected for the different time domain algorithms. The dimension of the applied filters is uniformly set to J = 128 expecting the ECA algorithm, where it has been set to J = 2.

Figure 15.

Zero range slice of the range-Doppler map with the different time domain algorithms.

It can be seen that the ECA-S algorithm has a fairly deep null in contrast to the other algorithms. With applying the supplemented autocorrelation matrix technique (S-ECA-S) this wide and deep null is filled up. Also observe that the RLS algorithm has evolved a relatively wide and shallow notch around the zero Doppler. This can be also seen in Figure 14.

Figure 16 shows the location of the investigated algorithms on the map of the filtering performance versus computational cost. The picture provides relevant information on the proper selection of the clutter cancellation algorithm. The variation of the achieved SINR improvements over the different algorithms is relatively small, however the essential features are readable. Recall that even high clutter cancellation can result in slight improvement, see in Figure 2. It must be also strictly emphasized that resultant ambiguities on the calculated range-Doppler map should also be taken into account to the algorithm selection. The following section aims to classify these algorithms based on the application criteria.

Figure 16.

Map of the characterized algorithms, generated according to their filter performance and computation costs (J = 128).

Slowly moving targets

For detecting slowly moving targets the information preservation in the low Doppler frequency domain is essential. The Wiener-SMI filter cancels only the zero Doppler contribution while provides great filtering performance. Using the supplemented autocorrelation matrix technique, the rapid coefficient calculation can also be guaranteed. In case the ECA-S algorithm is properly configured it can also be a reasonable choice, however its computation is more complex. The investigated iterative algorithms can only be applied in case the range of interest lies outside the dimension of the filter or the weight inheritance technique can be properly applied.

High filter performance

Among the examined filtering algorithms, the ECA-S method achieved the best filtering performance. The price paid for the superior filtering performance is the high computational burden.

Low computation cost

The iterative algorithms are essentially developed to ensure real-time operation. Several research paper reports about their successful implementations in FPGAs or GPUs [20, 21]. In passive radar application, the block LMS algorithm not only has low computation costs but also grants good filtering performance. The main drawback of these filters is the strong distortion of the range-Doppler map.


  1. 1. Kulpa K. Signal Processing in Noise Waveform Radar. 1st ed. Norwood, MA: Artech House; 2013
  2. 2. Di Lallo A, Fulcoli R, Timmoneri L. Adaptive spatial processing applied to a prototype passive. 2-5 Sept 2008; Adelaide, SA, Australia; IEE; 2008. pp. 1-5. DOI: 10.1109/RADAR.2008.4653906
  3. 3. Peto T, Seller R. Quad channel DVB-T based passive radar. 17th International Radar Symposium (IRS); 10-12 May 2016; Krakow, Poland; IEEE; 2016. pp. 1-4. DOI: 10.1109/IRS.2016.7497383
  4. 4. Tao R, Wu HZ, Shan T. Direct-path suppression by spatial filtering in digital television terrestrial broadcasting-based passive radar. IET Radar Sonar and Navigation. 2010;4(6):791-805. DOI: 10.1049/iet-rsn.2009.0138
  5. 5. Villano M, Colone F, Lombardo P. Antenna Array for passive radar: Configuration design and adaptive approaches to disturbance cancellation. International Journal of Antennas and Propagation. 2013;2013:16
  6. 6. Moscardini C, Conti M, Berizzi F, Martorella M, Capria A. Spatial adaptive processing for passive Bistatic radar. IEEE Radar Conference 2014; 19-23 May 2014; Cincinnati, OH, USA; IEEE; 2014. p. 6. DOI: 10.1109/RADAR.2014.6875751
  7. 7. Cardinali R, Colone F, Ferretti C, Lombardo P. Comparison of clutter and multipath cancellation techniques for passive radar. IEEE Radar Conference 2007; 17-20 April 2007; Boston, MA, USA; IEEE; 2007. p. 6. DOI: 10.1109/RADAR.2007.374262
  8. 8. Garry JL, Baker CJ, Smith GE. Evaluation of direct signal suppression for passive radar. IEEE Transactions on Geoscience and Remote Sensing. 2017;5(7):14. DOI: 10.1109/TGRS.2017.2680321
  9. 9. Baczyk MK, Malanowski M. Reconstruction of the reference signal in DVB-T-based passive radar. International Journal of Electronics and Telecommunication. 2011;57(1):43-48. DOI: 10.2478/v10177-011-0006-y
  10. 10. Liu J, Li H, Himed B. Analysis of cross-correlation detector for passive radar application. IEEE Radar Conference, 2015; 10-15 May 2015; Arlington, VA, USA; IEEE; 2015. p. 1-5. DOI: 10.1109/RADAR.2015.713110
  11. 11. Pető T, Seller R. Quad channel software defined receiver for passive radar application. Archives of Electrical Engineering. 2017;66(1):5-16. DOI: 10.1515/aeee-2017-0001
  12. 12. Pető T, Seller R. Time domain filter comparison in passive radar systems. 18th International Radar Symposium (IRS) 2017; 28-30 June 2017; Prague, Czech Republic; IEEE; 2017. p. 1-10. DOI: 10.23919/IRS.2017.8008108
  13. 13. Gross FB. Smart Antennas with MATLAB. 2nd ed. New York, Chicago, San Francisco, Athens, London, Madrid, Mexico City, Milan, New Delhi, Singapore, Sydney, Toronto: McGraw-Hill Education; 2015 433 p
  14. 14. Navrátil V, O’Brien A, Garry JL, Smith GE. Demonstration of space-time adaptive processing. International Radar Symposium 2018; 28-30 June 2017; Prague, Czech Republic; IEEE; 2017. p. 10. DOI: 10.23919/IRS.2017.8008146
  15. 15. Simon H. Adaptive Filter Theory. 3rd ed. Upper Saddle River, New Jersey: Prentice Hall Inc.; 1996
  16. 16. Malanowski M, Haugen R, Greco MS. Land and sea clutter from FM-based passive bistatic radars. IET Radar, Sonar and Navigation. 2013;8(2):160-166. DOI: 10.1049/iet-rsn.2013.0186
  17. 17. Colone F, O'Hangan DW, Lombardo P, Baker CJ. A multistage processing algorithm for disturbance removal and target detection in passive Bistatic radar. IEEE Transaction on Aerospace and Electronic Systems. 2009;45(2):698-722. DOI: 10.1109/TAES.2009.5089551
  18. 18. Colone F, Cardinali R, Lombardo P. Cancellation of clutter and multipath in passive radar using a sequential approach. IEEE Conference on Radar, 2006; 24-27 April, 2006; Verona, NY, USA; IEEE; 2006. p. 393-399. DOI: 10.1109/RADAR.2006.1631830
  19. 19. Colone F, Palmarini C, Martelli T, Tilli E. Sliding extensive cancellation algorithm for disturbance removal in passive radar. IEEE Transactions on Aerospace and Electronic Systems. 2016;52(3):1309-1326. DOI: 10.1109/TAES.2016.150477
  20. 20. Zhao YD, Zhao YK, Lu XD, Xiang MS. Block NLMS cancellation algorithm and its real-time implementation for passive radar. 14-16 April 2013; Xi'an, China; IET; 2013. pp. 1-5. DOI: 10.1049/cp.2013.0341
  21. 21. Oba H, Kim M, Arai H. FPGA implementation of LMS and N-LMS processor for adaptive array applications. International Symposium on Intelligent Signal Processing and Communications ISPACS '06; 12-15 Dec 2006; Tottori, Japan; IEEE; 2007. p. 485-488. DOI: 10.1109/ISPACS.2006.364703

Written By

Tamás Pető and Rudolf Seller

Submitted: 24 May 2017 Reviewed: 26 September 2017 Published: 20 December 2017