APRiL

Advanced Passive Radar Library

Detector Module

This notebook is intended to give an introduction on the application of the detector module implemented in pyAPRiL. This module encomapasses different algorithms to calculate the so called range-Doppler matrix using the reference and surveillance channels. Generally (but not exclusively), to obtain the range-Doppler matrix, it is required to have the time domain samples of these two channels and know the sampling frequency of the signals to prepare the proper scaling. In the next few sections we are going through the realized range-Doppler calculation functions in pyApril and demonstrate their operation on simulated data.

# Import third-party packages
import numpy as np
import plotly.graph_objects as go
from plotly.offline import iplot, init_notebook_mode
init_notebook_mode(connected = False)
import matplotlib as plt
%matplotlib inline

# Import pyAPRiL packages
from pyapril import detector

1. Test signal simulation

As a first step we are going to preapre simulated reference and surveillance channel signals. The simulated illuminator signal is drawn from normal distribution with unit variance $(\sigma_s = 1)$ and 0 mean ($\mathbb{E}\{s[n]\} = 0$) . \begin{equation} s[n] \sim \mathcal{N}(0,1) \end{equation} The simulated illuminator signal has $n=0 \dots N-1$ samples. After generating the transmitted signal the reference and surveillance channels are prepared accrodingly \begin{equation} x_r[n] = s[n] \\ x_sn] = s[n-\tau_0] e^{j 2 \pi f_0/f_s n}, \end{equation} where $\tau_0$ and $f_0$ are time delay and the Doppler shift of the simulated target reflection and $f_s$ is the sampling frequency of the signals.

This signal modell and simulation is of course a way too simple, but our aim here is not to create a sophisticated scenario but only to show the operation of the implemented algorithms.

# Simulation parameters
N     = 2**18 # Number of samples in a coherently processed block
fs    = 8*10**6 # Sampling frequency [Hz]
tau_0 = 12 # Target time delay [sample]
f_0   = 50 # Target Doppler frequency [Hz]

# Samples of the illuminator signal are drawn from normal distribution
s = np.random.normal(0,1,N+tau_0) +1j* np.random.normal(0,1,N+tau_0)

# Prepare reference and surveillance channels
x_r    = np.zeros(N, dtype=complex) # Reference channel allocation
x_s    = np.zeros(N, dtype=complex) # Surveillance channel allocation
x_s[:] = s[0:N] * np.exp(1j*2*np.pi*f_0/fs*np.arange(N)) # Target reflection is Doppler shifted has a time lag relative to the reference channel
x_r[:] = s[tau_0:tau_0+N]  # Prepare reference channel with cutting out the relevant signal portion

2. Detectors

Algorithms implemented herein realize the matched filtering. According to this procedure target reflection prototype signals are generated from the referece channel and each of these prototypes are cross-correlated with the surveillance channel. The detector module have three implementations for the same operation differing in their computational complexity. The detailed description of these algorithms can be found in [1] "Chapter 4. - Correlation Processing".

2.1 cc_detector_td:

rd_matrix = cc_detector_td(ref_ch, surv_ch, fs, fD_max, fD_s, r_max, verbose=0, Qt_obj=None)

This is the raw, time-domain implementation of the matched filter accroding to equation (4.9) in [1]. Using the notation introduced above: \begin{equation} \chi(\tau,f) = \sum_{n=0}^{N-1} x_r[n - \tau]^* e^{-j 2 \pi f/f_s n} x_s[n], \end{equation}

Parameters:

• ref_ch: 1D complex numpy array (1 x N)

  Reference channel samples (1..N)

• surv_ch: 1D complex numpy array (1 x N)

  Surveillance channel samples (1..N)

• fs: float ( fs >0 )

  Sampling frequency of the signal [Hz]

• fD_max: float (fD_max >0 )

  Maximum Doppler frequency for which the range-Doppler matrix will be calculated. This parameter should be set according to the maximum expected bistatic Doppler frequency of the target. The maximum Doppler frequency (assuming monostatic geometry) is $f_{max} = \frac{2v}{\lambda} $, where $\lambda$ is the wavelength of the center frequency. E.g.: If the maximum target velocity is $15 \frac{m}{s}$, and the center frequency of the transmitted illuminator signal is $f_c = 1 GHz$, then $\lambda = \frac{c}{f_c} \simeq 0.3 m $ and correspondingly $f_{max} = 100 Hz$.

• fD_s: float ( fD_max > fD_s > 0)

  Doppler frequency step size [Hz]. The Doppler frequency domain will be discretized using this step size. Increasing this number will result in finer sampling, more cells in the Doppler dimension, however note that the Doppler freuqency resolution will be the same as it is controlled by the length of the coherent processing interval, thus the time duration of the input reference and surveillance channel vectors. $T_{CPI} = \frac{1}{f_s} N $.

• r_max: float ( $N \geq r_{max} > 0$ )

  Maximum range for which the range-Doppler matrix will be calculated. This parameter should be set according to the maximum expected bistatic range of the target. This input value is interpreted in sample delays. Each of this range cells corresponds to $\frac{c}{fs}$ bistatic range from the radar. See section 2.2 [1].

• verbose: bool

  Verbose mode, 0 means disabled. If enabled the calcaultion progress is printed to the standard output

• Qt_obj: bool

  This object is used to transfer progress update to a Qt based GUI. In case a Qt progressbar object is specified, its emit(progress) function is progressively called during the calculation.

ref_ch and surv_ch input vectors must have the same size (N)

Returns:

• rd_matrix 2D complex numpy array (D x R)

  The output of the function is the calculated complex valued range-Doppler matrix, that has $D = 2 f_{Dmax}/f_{Ds} +1 $ Doppler cells and $R=r_{max}$ range cells.

The following code section shows the calling of this function for the simulated data channels.

# Detector parameters
fD_max  = 200 # The simulated target signal has 50 Hz Doppler frequency.
fD_s    = 1 
r_max   = 64 # The simulated sample delay of the target reflection is 12.
rd_matrix = detector.cc_detector_td(x_r, x_s, fs, fD_max, fD_s, r_max, verbose=0, Qt_obj=None)

After calculating the range-matrix we can display it. When we are considering detection, normally we are only interested in the absolute value of the matrix.

# Plot the range-Doppler matrix using plotly
fig = go.Figure()       
fig.add_trace(go.Heatmap(z=abs(rd_matrix)))
fig.update_layout(xaxis=dict(title="Time delay [sample] / Bistatic range cell"),
                  yaxis=dict(title="Doppler frequency [bin]"))

iplot(fig)

We can see that the calculated range-Doppler matrix has a maximum at the simulated target coordinates (bistatic range and Doppler frequency cells). We can of course use better scalling and display for inspecting the matrix. pyAPRiL has functions implemented to ease this part of the job. We are not dealing here with the operation of these functions, just only use them.

from pyapril.RDTools import plot_rd_matrix

dyn_range = 10*np.log10(N) # Dynamic range of the ambiguity function of the noise illuminator signal
plot_rd_matrix(rd_matrix,  
                   dyn_range=dyn_range,
                   interpolation='best', 
                   cmap='jet', 
                   scaling_mode="adaptive-floor-fix-range",
                   fs=fs,
                   max_Doppler=fD_max)

2.2 cc_detector_fd:

rd_matrix = cc_detector_fd(ref_ch, surv_ch, fs, fD_max, r_max, verbose=0, Qt_obj=None)

This function differs in a way from cc_detector_td, that it calculates the cross-correlation between the channels in frequncy domain. This modification greatly speeds up the overall calculation procedure. Beside this the Doppler frequency shift of the reference channel is also realized in the frequency domain. As direct consequence the Doppler frequency dimension of the range-Doppler matrix can only be sampled at the frequency bin size of the transformed channels. Due to zero-padding (required for the correlation processing), this bin size is equal to $\Delta_f = \frac{fs}{2N}$.

The implementation use the FFTW library (through pyFFTW) to improve the calculation speed of the forward and backward Fourier transforms.

Among others, this algorithm is described in [1], Section 4.3.2 as "Method 2".

Parameters:

• ref_ch: 1D complex numpy array (1 x N)

  Reference channel samples (1..N)

• surv_ch: 1D complex numpy array (1 x N)

  Surveillance channel samples (1..N)

• fs: float ( fs >0 )

  Sampling frequency of the signal [Hz]

• fD_max: float (fD_max > $\Delta_f$ )

  Maximum Doppler frequency for which the range-Doppler matrix will be calculated. This parameter should be set according to the maximum expected bistatic Doppler frequency of the target. The maximum Doppler frequency (assuming monostatic geometry) is $f_{max} = \frac{2v}{\lambda} $, where $\lambda$ is the wavelength of the center frequency. E.g.: If the maximum target velocity is $15 \frac{m}{s}$, and the center frequency of the transmitted illuminator signal is $f_c = 1 GHz$, then $\lambda = \frac{c}{f_c} \simeq 0.3 m $ and correspondingly $f_{max} = 100 Hz$.

• r_max: float ( $N \geq r_{max} > 0$ )

  Maximum range for which the range-Doppler matrix will be calculated. Maximum range for which the range-Doppler matrix will be calculated. This parameter should be set according to the maximum expected bistatic range of the target. This input value is interpreted in sample delays. Each of this range cells corresponds to $\frac{c}{fs}$ bistatic range from the radar. See section 2.2 [1].

• verbose: bool

  Verbose mode, 0 means disabled. If enabled the calcaultion progress is printed to the standard output

• Qt_obj: bool

  This object is used to transfer progress update to a Qt based GUI. In case a Qt progressbar object is specified, its emit(progress) function is progressively called during the calculation.

ref_ch and surv_ch input vectors must have the same size (N)

Returns:

• rd_matrix 2D complex numpy array (D x R)

  The output of the function is the calculated complex valued range-Doppler matrix, that has $D = 2 f_{Dmax}/\Delta_{f} +1 $ Doppler cells and $R=r_{max}$ range cells.

# Detector parameters
fD_max  = 200 # The simulated target signal has 50 Hz Doppler frequency.
r_max   = 64 # The simulated sample delay of the target reflection is 12.

# Run CC-Detector in frequency domain
rd_matrix = detector.cc_detector_fd(x_r, x_s, fs, fD_max, r_max, verbose=0, Qt_obj=None)

# Plot the obtained range-Doppler matrix
plot_rd_matrix(rd_matrix,  
                   dyn_range=dyn_range,
                   interpolation='best', 
                   cmap='jet', 
                   scaling_mode="adaptive-floor-fix-range",
                   fs=fs,
                   max_Doppler=fD_max)

2.3 cc_detector_ons:

rd_matrix = cc_detector_ons (ref_ch, surv_ch, fs, fD_max, r_max, verbose=0, Qt_obj=None)

This algorithm is very similar to what is usually done in FMCW and pulse radars. The overall processing is separated into two steps. In the first step the range processing is performed via the calculation of the cross-correlation functions on short time intervals ($N/r_{max} = k, k \in \mathbb{Z}^+ $ batches). In the second step Fourier transformation is performed across the range vectors, thus executing the Doppler processig. The main benefit of this algorithm is the low calculation complexity, however faster processing comes at the expense of loosing processing gain at higher Doppler frequencies. This loss is expressed in [2], in section 4.3, equation (25). The algorithm is described in more detail in [1], Section 4.3.3 as "Method 3" and in [2]

The Doppler frequency bin size is equal to $\Delta_f = \frac{fs}{N}$.

Parameters:

• ref_ch: 1D complex numpy array (1 x N)

  Reference channel samples (1..N)

• surv_ch: 1D complex numpy array (1 x N)

  Surveillance channel samples (1..N)

• fs: float ( fs >0 )

  Sampling frequency of the signal [Hz]

• fD_max: float ( $k \Delta_f$>fD_max > $\Delta_f$ )

  Maximum Doppler frequency for which the range-Doppler matrix will be calculated. This parameter should be set according to the maximum expected bistatic Doppler frequency of the target. The maximum Doppler frequency (assuming monostatic geometry) is $f_{max} = \frac{2v}{\lambda} $, where $\lambda$ is the wavelength of the center frequency. E.g.: If the maximum target velocity is $15 \frac{m}{s}$, and the center frequency of the transmitted illuminator signal is $f_c = 1 GHz$, then $\lambda = \frac{c}{f_c} \simeq 0.3 m $ and correspondingly $f_{max} = 100 Hz$. fDmax can not exceed the overall span of the unambigious Doppler frequency. With increase of $r{max}$, fD_max deacreses.

• r_max: float ( $N \geq r_{max} > 0, N/r_{max} = k, k \in \mathbb{Z}^+ $)

  Maximum range for which the range-Doppler matrix will be calculated. Maximum range for which the range-Doppler matrix will be calculated. This parameter should be set according to the maximum expected bistatic range of the target. This input value is interpreted in sample delays. Each of this range cells corresponds to $\frac{c}{fs}$ bistatic range from the radar. See section 2.2 [1]. This value should partionate the overall processing interval (with $N$ samples) to multiple batches ($k$)

• verbose: bool

  Verbose mode, 0 means disabled. If enabled the calcaultion progress is printed to the standard output

• Qt_obj: bool

  This object is used to transfer progress update to a Qt based GUI. In case a Qt progressbar object is specified, its emit(progress) function is progressively called during the calculation.

ref_ch and surv_ch input vectors must have the same size (N)

Returns:

• rd_matrix 2D complex numpy array (D x R)

  The output of the function is the calculated complex valued range-Doppler matrix, that has $D = 2 f_{Dmax}/\Delta_{f} +1 $ Doppler cells and $R=r_{max}$ range cells.

# Detector parameters
fD_max  = 200 # The simulated target signal has 50 Hz Doppler frequency.
r_max   = 64 # The simulated sample delay of the target reflection is 12.

# Run CC-Detector in frequency domain
rd_matrix = detector.cc_detector_ons(x_r, x_s, fs, fD_max, r_max, verbose=0, Qt_obj=None)

# Plot the obtained range-Doppler matrix
plot_rd_matrix(rd_matrix,  
                   dyn_range=dyn_range,
                   interpolation='best', 
                   cmap='jet', 
                   scaling_mode="adaptive-floor-fix-range",
                   fs=fs,
                   max_Doppler=fD_max)

2.4 Comparison

Runtime

import timeit
fD_max  = 200
fD_s    = fs / (2*N) # Resolutions will be the same
r_max   = 64 

runtime_td  = timeit.timeit("detector.cc_detector_td(x_r, x_s, fs, fD_max, fD_s, r_max, verbose=0, Qt_obj=None)", globals=globals(), number=20)/20
runtime_fd = timeit.timeit("detector.cc_detector_fd(x_r, x_s, fs, fD_max, r_max, verbose=0, Qt_obj=None)", globals=globals(), number=20)/20
runtime_ons = timeit.timeit("detector.cc_detector_ons(x_r, x_s, fs, fD_max, r_max, verbose=0, Qt_obj=None)", globals=globals(), number=20)/20

print("CC-Detector Time Domain: {:.3f} ms".format(runtime_td*1000))
print("CC-Detector Frequency Domain: {:.3f} ms".format(runtime_fd*1000))
print("CC-Detector ONS (Batched): {:.3f} ms".format(runtime_ons*1000))

CC-Detector Time Domain: 20774.744 ms
CC-Detector Frequency Domain: 435.173 ms
CC-Detector ONS (Batched): 69.832 ms

References

[1] Mateusz Malanowski: Signal Processing for Passive Bistatic Radar, 2019 ARTECH HOUSE 685 Canton Street Norwood, MA 02062

[2] Michał Meller: Efficient Signal Processing Algorithms for Passive Radars, NATO STO-MP-SET-187

Further Work

Range and Doppler walk compensations: Using such techniques one can compensate the energy spread of the target in the range-Doppler matrix when the Doppler frequency or the range of the target change rapidly relative to the duration of coherent processing interval.

• [F-RWC1] Karl Erik Olsen, Jonas Myhre Christiansen: Multiband Passive Bistatic Radar Coherent Range and Doppler-Walk Compensation, 2015 IEEE Radar Conference (RadarCon), Arlington, VA, USA, 10-15 May 2015
• [F-RWC2] Graeme E. Smith, Kevin Chetty, Christopher John Baker, Karl Woodbridge3: Extended time processing for passive bistatic radar, IET Radar, Sonar and Navigation, ISSN 1751-8784
• [F-RWC3] Mateusz Malanowski: Signal Processing for Passive Bistatic Radar, 2019 ARTECH HOUSE 685 Canton Street Norwood, MA 02062, (Section 4.6)

GLRT detector: The Generalized Likelihood Ratio Test detectors can outperform the classic cross-correlation detector when the reference signal is considered to be noisy.

• [F-GLRT1] Guolong Cui, Jun Liu and Hongbin Li, Braham Himed:Target Detection for Passive Radar with Noisy Reference Channel, IEEE Transactions on Aerospace and Electronic Systems ( Volume: 56, Issue: 6, December 2020)
• [F-GLRT2] Xin Zhang, Hongbin Li, Braham Himed A Direct-Path Interference Resistant Passive Detector, IEEE Signal Processing Letters ( Volume: 24, Issue: 6, June 2017)

Reciprocal filter: "providing for high precision delay-time estimation with intensely reduced sidelobes. Using a FFT of several consecutive correlation functions, even the echo Doppler shift can be estimated in spite of the Doppler-tolerance of the autocorrelation maximum. Using the reciprocal filter, especially signals with a merely small spectral amplitude variation turn out to be favourable concerning the SNR. Varying the threshold β , however, even a higher SNR can be achieved at the cost of a reduced MSR and vice versa."[FR-1]

• [F-R1] Dr. Martin Glende: PCL-Signal-Processing for Sidelobe Reduction in Case of Periodical Illuminator Signals, 2006 International Radar Symposium, Krakow, Poland , 24-26 May 2006

• [F-R2] Philipp Wojaczek, Fabiola Colone, Diego Cristallini, Pierfrancesco Lombardo, Heiner Kuschel: The Application of the Reciprocal Filter and DPCA for GMTI in DVB-T – PCL, International Conference on Radar Systems (Radar 2017), Belfast, 23-26 October 2017

Authors

dr. Tamas Peto, PhD
Initial version: 2023 03 05