Detection at a constant false alarm probability (FALT)¶
This example discusses Detection with Constant False Alarm Probability (CFAP) and shows how to use CFARDetector and CFARDetector2D from the EngeePhased library to implement the detection algorithm with CFAP
Functions used¶
Pkg.add(["DSP"])
include("$(@__DIR__)/functions.jl");
Introduction¶
One of the important tasks that a radar system performs is target detection. The detection procedure itself is quite simple. It compares the signal to a threshold value. Therefore, the real work of detection is to determine the appropriate threshold. In general, the threshold is a function of both the probability of correct detection and the probability of false alarm.
In many phased array systems, because of the costs associated with false detection, it is desirable to have a detection threshold that not only maximises the probability of correct detection, but also keeps the probability of false alarm below a given level.
There is an extensive literature on how to define a detection threshold. Readers may be interested in the examples of signal detection in white Gaussian noise and multi-sample signal detection, where some well-known results are given. However, all these classical results are based on theoretical probabilities and are limited to white Gaussian noise with known variance (power). In real devices, the noise is often coloured and its power is unknown.
The CFAR algorithm solves these problems. When a particular cell, often referred to as a Cell Under Test (CUT), needs to be detected, the noise power is estimated from neighbouring cells. In such a case, the detection threshold (T) is defined as follows:
$\begin{align} T = \alpha P_n \end{align}$
where $P_n$ is the estimation noise power, $\alpha$ is the scaling (threshold) factor
It can be seen from the equation that the threshold adapts to the data. It can be shown that with an appropriate threshold factor α, the resulting false alarm probability can be kept constant, hence the name PVLT.
Cell averaging. Detection with PVLT¶
The cell-averaged CFAR detector is probably the most widely used CFAR detector. It is also used as a baseline comparison for other CFAR methods. In the cell-averaged CFAR detector, noise samples are extracted from leading and lagging cells (called training cells) around the CUT. The noise estimate can be calculated as:
where $N$ is the number of training cells and $x_m$ is the sample in each training cell.
While $x_m$ is the output of the quadratic law detector, $P_n$ is the estimated noise power. In general, the number of leading and lagging training cells is the same. The guard cells are placed adjacent to the CUT, both ahead and behind it. The purpose of these guard cells is to avoid signal components leaking into the training cell, which can adversely affect the noise estimate.
The following figure shows the relationship between these cells for the 1-D case.
In the above cell-averaged CFAR detector, assuming that the data into the detector comes from a single pulse, i.e. without pulse integration, the threshold factor can be written as:
$\begin{align} \alpha = N(P_{fa}^{-1/N}-1) \end{align}$ where $P_{fa}$ is the required false alarm probability.
CFAR detection using automatic threshold factor¶
In the rest of this example, we show how to use EngeePhased to perform cell-averaged CFAR detection. For simplicity and without loss of generality, we still assume that the noise is white Gaussian. This allows us to compare CFAR with classical detection theory.
We can create a CFAR detector using the following command:
cfar = EngeePhased.CFARDetector(NumTrainingCells=20,NumGuardCells=2);
In this detector, we use 20 training cells and 2 guard cells. This means that there are 10 training cells and 1 guard cell on each side of the CUT. As mentioned above, if we assume that the signal comes from a quadratic law detector without pulse integration, the threshold can be calculated based on the number of training cells and the desired false alarm probability. Assuming that the desired false alarm probability is 0.001, we can configure the CFAR detector as follows to allow for this calculation.
exp_pfa = 1e-3 # вероятность ложной тревоги
cfar.ThresholdFactor = "Auto" # алгоритм определения порога (Автоматический)
cfar.ProbabilityFalseAlarm = exp_pfa;
The configured CFAR detector is shown below.
cfar
We will now simulate the input data. Since we want to show that the CFAR detector can keep the false alarm rate below a certain value, we will simply model the noise in these cells. Here are the settings:
- The data sequence is 23 samples long, and the CUT is cell 12. This leaves 10 training cells and 1 guard cell on each side of the CUT.
- The false alarm rate is calculated using 100,000 Monte Carlo tests.
rng = MersenneTwister(20)
npower = db2pow(-10); # Задание ОСШ - 10dB
Ntrials = Int(1e5);
Ncells = 23;
CUTIdx = 12;
# Модель шума после квадратичного детектора
rsamp = randn(rng,Ncells,Ntrials)+1im*randn(rng,Ncells,Ntrials);
x = abs.(sqrt.(npower/2)*rsamp).^2;
To perform detection, run the data through the detector. In this example, there is only one CUT, so the output is a logic vector containing the detection result for all trials. If the result is true, it means that the target is present in the corresponding trial. In our example, all detections are false alarms since we only pass in noise. As a result, we can calculate the false alarm rate based on the number of false alarms and the number of trials.
$\begin{align} \alpha = N(P_{fa}^{-1/N}-1) \end{align}$
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials;
println("act_pfa = $(act_pfa)")
The result shows that the resulting false alarm probability is 0.001, as we expected.
CFAR detection using custom threshold factor¶
As explained in the previous part of this example, there are only a few cases where the CFAR detector can automatically calculate the appropriate threshold factor. For example, in the previous scenario, if we use 10-pulse incoherent integration before the data arrives at the detector, the automatic threshold will no longer be able to provide the desired false alarm rate:
npower = db2pow(-10); # значение ОСШ - 10 дБ
xn = 0;
for m = 1:10
rsamp = randn(rng,Ncells,Ntrials) .+ 1im*randn(rng,Ncells,Ntrials);
xn = xn .+ abs.(sqrt.(npower/2)*rsamp).^2 # некогерентное накопление
end
x_detected = cfar(xn,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
It may be puzzling why we consider a false alarm rate of 0 to be worse than a false alarm rate of 0.001. After all, isn't a false alarm rate equal to 0 a great thing? The answer to this question lies in the fact that when the false alarm rate decreases, the probability of detection also decreases. In this case, since the true false alarm probability is well below the acceptable value, the detection threshold is set too high. The same probability of detection can be achieved with the desired probability of false alarm at a lower cost; for example, with less transmitter power.
In most cases, the threshold factor needs to be determined based on the specific environment and system configuration. We can configure the CFAR detector to use a custom threshold factor as shown below.
release!(cfar);
cfar.ThresholdFactor = "Custom";
Continuing with the pulse integration example and using empirical data, we find that a custom threshold factor of 2.35 can be used to achieve the desired false alarm rate. Using this threshold, we see that the resulting false alarm rate is in line with the expected value.
cfar.CustomThresholdFactor = 2.35;
x_detected = cfar(xn,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
CFAR detection threshold¶
CFAR detection occurs when the input signal level in a cell exceeds a threshold level. The threshold level for each cell depends on the threshold factor and the noise power in that cell derived from the training cells. To keep the false alarm rate constant, the detection threshold will increase or decrease in proportion to the noise power in the training cells. Configure the CFAR detector to output the threshold used for each detection using the ThresholdOutputPort property. Let's set the mode of automatic calculation of the threshold coefficient and 200 training cells.
release!(cfar);
cfar.ThresholdOutputPort = true;
cfar.ThresholdFactor = "Auto";
cfar.NumTrainingCells = 200;
Then generate a quadratic input signal with increasing noise power.
Npoints = 10_000;
rsamp = randn(rng,Npoints).+1im*randn(rng,Npoints);
ramp = Vector(range(1,10;length=Npoints));
xRamp = abs.(sqrt.(npower*ramp./2).*rsamp).^2;
Calculate the detections and thresholds for all cells in the signal.
x_detected,th = cfar(xRamp,Vector(1:length(xRamp)));
Then compare the CFAR threshold value with the input signal.
gr()
plot(1:length(xRamp),xRamp,label="Сигнал",line=(2,:blue,:sticks))
plot!(1:length(xRamp),th,label="Порог")
scatter!(findall(x_detected[:] .> 0),xRamp[findall(x_detected[:] .> 0)],label="Обнаружение",marker=(:circle,4,:yellow))
xlabel!("Временные отсчеты")
ylabel!("Амплитуда",legend=:topleft)
Here, the threshold increases as the noise power in the signal increases to keep the false alarm level constant. Detection occurs where the signal level exceeds the threshold value.
Comparison of CFAR with the classical Neyman-Pearson detector¶
In this section, we compare the performance of the CFAR detector with classical detection theory using the Neyman-Pearson principle. Returning to the first example and assuming that the true noise power is known, the theoretical threshold can be calculated as
T_ideal = npower*db2pow(npwgnthresh(exp_pfa));
The false alarm rate of this classical Neyman-Pearson detector can be calculated using this theoretical threshold.
act_Pfa_np = sum(x[CUTIdx,:] .> T_ideal)/Ntrials
println("act_Pfa_np = $act_Pfa_np")
Since we know the noise power, classical detection theory also allows us to obtain the desired false alarm rate. The false alarm rate achieved by the CFAR detector is similar.
release!(cfar);
cfar.ThresholdOutputPort = false;
cfar.NumTrainingCells = 20;
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
Next, suppose that both detectors are deployed in the field and the noise power is 1 dB higher than expected. In this case, if we use the theoretical threshold, the resulting false alarm rate will be four times what we want.
npower = db2pow(-9); # Assume 9dB SNR ratio
rsamp = randn(rng,Ncells,Ntrials).+1im*randn(rng,Ncells,Ntrials);
x = abs.(sqrt.(npower/2)*rsamp).^2;
act_Pfa_np = sum(x[CUTIdx,:] .> T_ideal)/Ntrials
println("act_Pfa_np = $act_Pfa_np")
In contrast, the performance of the CFAR detector is not affected.
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
Thus, the CFAR detector is robust to noise power uncertainty and is better suited for field applications.
Finally, we use CFAR detection in the presence of coloured noise. We first apply a classical detection threshold to the data.
npower = db2pow(-10);
fcoeff = [-0.00057693500095478268512028119374691, 0.0004446399942713034555974438433168, 0.023749844987587823141872434007382, 0.10475951999236174372320817838045, 0.22682709001336692766770397611253, 0.2895916800267339108465591834829, 0.22682709001336692766770397611253, 0.10475951999236174372320817838045, 0.023749844987587819672425482053768, 0.0004446399942713034555974438433168, -0.00057693500095478268512028119374691]
x = abs.(sqrt(npower/2)*filt(fcoeff,1,rsamp)).^2; # colored noise
act_Pfa_np = sum(x[CUTIdx,:] .> T_ideal)/Ntrials
println("act_Pfa_np = $act_Pfa_np")
Note that the resulting false alarm rate cannot meet the requirements. However, by using a CFAR detector with a custom threshold factor, we can obtain the desired false alarm rate.
release!(cfar);
cfar.ThresholdFactor = "Custom";
cfar.CustomThresholdFactor = 12.85;
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials
println("act_pfa = $act_pfa")
CFAR detection by range-Doppler portrait¶
In the previous sections, noise estimation was computed based on training cells leading and lagging the CUT in one dimension. We can also perform CFAR detection on images. The cells correspond to pixels in the images, and the guard and training cells are placed in bands around the CUT. The detection threshold is calculated using the cells in the rectangular training strip around the CUT.
In the above figure, the size of the guard band is [2 2] and the size of the training band is [4 3]. The size indices denote the number of cells on each side of the CUT in row and column respectively. The guard bar size can also be defined as 2, since its size is the same across rows and columns.
Next, create a two-dimensional CFAR detector. Use a false alarm probability of 1e-5 and specify a guard band size of 5 cells and a training band size of 10 cells.
cfar2D = EngeePhased.CFARDetector2D(GuardBandSize=5,TrainingBandSize=10,
ProbabilityFalseAlarm=1e-5);
Then load and construct a range Doppler image. The image includes returns from two stationary targets and one target moving away from the radar.
hrdresp = EngeePhased.RangeDopplerResponse(
DopplerFFTLengthSource="Property",
DopplerFFTLength=RangeDopplerEx_MF_NFFTDOP,
SampleRate=RangeDopplerEx_MF_Fs);
resp,rng_grid,dop_grid = step!(hrdresp,RangeDopplerEx_MF_X,RangeDopplerEx_MF_Coeff);
#Добавление белого шума
npow = db2amp(-95);
resp = resp .+ npow/sqrt(2)*(randn(rng,size(resp)) .+ 1im*randn(rng,size(resp)))
plot_1 = heatmap(dop_grid,rng_grid,mag2db.(abs.(resp)),colorbar=false,margin=5Plots.mm,
color = :rainbow1, titlefont = font(12,"Computer Modern"),
guidefont = font(10,"Computer Modern"))
resp = abs.(resp).^2;
xlabel!("Частота Доплера, Гц")
ylabel!("Дальность, м")
title!("Дальностно-доплеровский портрет")
rangeIndx = map(x -> x[1],(findmin(abs.(rng_grid[:] .- [1000 4000]);dims=1)[2]))
dopplerIndx = map(x -> x[1],(findmin(abs.(dop_grid[:] .- [-1e4 1e4]);dims=1)[2]))
columnInds,rowInds = meshgrid(dopplerIndx[1]:dopplerIndx[2],rangeIndx[1]:rangeIndx[2]);
CUTIdx = permutedims([rowInds[:] columnInds[:]])
Calculate the detection result for each cell under test. In this example, a cell is each pixel in the search area. Build a detection result map for the range-Doppler image.
detections = cfar2D(resp,CUTIdx)
detectionMap = zeros(size(resp))
detectionMap[rangeIndx[1]:rangeIndx[2],dopplerIndx[1]:dopplerIndx[2]] =
reshape(float(detections),rangeIndx[2]-rangeIndx[1]+1,dopplerIndx[2]-dopplerIndx[1]+1);
heatmap(dop_grid,rng_grid,mag2db.(abs.(detectionMap .+ 1eps(1.))),colorbar=false,c = :imola,
titlefont = font(12,"Computer Modern"),guidefont = font(10,"Computer Modern"),margin=5Plots.mm)
xlabel!("Частота Доплера, Гц")
ylabel!("Дальность, м")
title!("Дальностно-доплеровский портрет после обнаружения CFAR")
Three objects are detected. A cube of range Doppler image data over time can also be provided as input to cfar2D, and the detection will be calculated in one step.
Conclusion¶
In this example, we have presented the basic concepts behind CFAR detectors. In particular, we have examined how to use the EngeePhased toolkit to perform a cell-averaged CFAR detector on signal and range Doppler images. Comparing the performance of the cell-averaged CFAR detector and the detector with the theoretically calculated threshold clearly shows that the CFAR detector is more suitable for real field conditions.