Engee documentation
Notebook

Detection Curves (ROC) Part 2: Monte Carlo Simulation

This example shows how to plot the receiver detection curve (ROC curve) of a radar system using Monte Carlo simulation. The ROC determines how well the system can detect targets while rejecting large values of false signals when there is no target (false alarms). The detection system determines whether or not a target is present by comparing the received signal value to a predetermined threshold.

Probability of correct detection (Pd) is the probability that the instantaneous signal value is greater than the threshold value when the target is actually present.

The probability of false alarm (Pfa) is the probability that the signal value is greater than the threshold when the target is not present. In this case, the signal is due to noise and its properties depend on the noise statistics. Monte Carlo simulation generates a very large number of radar signals in the presence and absence of a target.

The ROC curve is a plot of the dependence of Pd on Pfa. The shape of the ROC curve depends on the signal-to-noise ratio (SNR) of the received signal. If SNR of the received signal is known, the ROC curve shows how reliable the system is by analysing the values of Pd and Pfa. By specifying Pd and Pfa, it is possible to determine how much power is needed to fulfil this requirement.

You can use the rocsnr function to calculate theoretical ROC curves. This example shows the ROC curve from a Monte Carlo simulation of a single-antenna radar system and compares it with the theoretical curve.

Tactical and technical characteristics of the radar

Set the desired Pd = 0.9 and Pfa = $10^{-6}$, set the maximum radar range to 4000 metres and range resolution to 50 metres. Set the actual target range to 3000 metres. Also set the effective scattering area (EDA) of the target to 1.5 $м^2$ and the operating frequency to 10 GHz.

In [ ]: 
Pkg.add(["DSP"])
using DSP

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`
In [ ]: 
using DSP

c = 299_792_458; # скорость распространения сигнала в пространстве, м/с
pd = 0.9; # вероятность правильного обнаружения
pfa = 1e-6; # вероятность ложной тревоги
max_range = 4000; # максимальная дальность, м
target_range = 3000.0; # дальность до цели, м
range_res = 50; # разрешение по дальности,м
tgt_rcs = 1.5; # ЭПР цели, м^2
fc = 10e9; # несущая частота, Гц
lambda = c/fc; # длина волны, м

Any modelling that calculates pfa and pd requires processing many signals. To avoid taking up a lot of memory, we will only process portions of the pulses. Let's set the number of pulses to be processed to 45_000 and the size of each part to 10_000.

In [ ]: 
Npulse = 45_000; # количество импульсов
Npulsebuffsize = 10_000; # количество отсчетов в импульсе

Selecting the type of probing signal (PS) and its parameters

Let's form the system object of the LFM signal generator. For this purpose it is necessary to calculate:

  • ** Pulse band**, using the range resolution;
  • ** Pulse repetition frequency (PRF)** based on the maximum radar range.
  • Sampling rate (fs) set equal to twice the bandwidth, since the signal is wideband;
  • ** Pulse duration** based on the bandwidth.
In [ ]: 
pulse_bw = c/(2*range_res); # полоса импульса
prf = c/(2*max_range); # частота следования импульсов
fs = 2*pulse_bw; # частота дискретизации
pulse_duration = 10/pulse_bw; # длительность импульса

# Формирование генератора ЛЧМ-сигнала
waveform = EngeePhased.LinearFMWaveform(
    PulseWidth=pulse_duration,
    SampleRate=fs,
    SweepBandwidth=pulse_bw,
    PRF=prf
);

Achieving the specified Pd and Pfa requires that sufficient signal power is delivered to the receiver after the target reflects the signal. Let us calculate the minimum SNR required to achieve a given false alarm probability and detection probability using the Albersham equation.

In [ ]: 
snr_min = albersheim(pd,pfa); # минимальное значение ОСШ
println("ОСШ = $(round(snr_min;sigdigits=4)) дБ")

ОСШ = 13.11 дБ

To achieve an OSR = 13 dB, sufficient power must be transmitted to the target. Use the radar equation to estimate the peak transmit power required to achieve a given ROS for a target at a range of 3000 metres. The received signal also depends on the EPR, which is assumed to follow a non-fluctuating model (Swerling 0). Let us set the radar to the same receive and transmit gains of 20 dB.

In [ ]: 
txrx_gain = 20; # усиление передатчика и приемника

# Расчет дальности с помощью основного уравнения радилокации
peak_power = ((4*pi)^3*noisepow(1/pulse_duration)*target_range^4*
    DSP.db2pow(snr_min))/(DSP.db2pow(2*txrx_gain)*tgt_rcs*lambda^2) # 
println("peak_power = $(round(peak_power;sigdigits=5)) Вт")

peak_power = 293.18 Вт

Setting the parameters of the transmitter system objects

Let's create system objects that define the transmitter part of the simulation: antenna motion model, antenna geometry, transmitter and radiator.

In [ ]: 
# модель движения радара
antennaplatform = EngeePhased.Platform(
    InitialPosition=[0; 0; 0], # начальное положение радара
    Velocity=[0; 0; 0] # компоненты вектора скорости радара
); 
# антенный элемент
antenna = EngeePhased.IsotropicAntennaElement(
    FrequencyRange=[5e9 15e9] # диапазон частот
);
# передатчик
transmitter = EngeePhased.Transmitter(
    Gain=txrx_gain, # усиление
    PeakPower=peak_power, # пиковая мощность
    InUseOutputPort=true # состояние радара
);
radiator = EngeePhased.Radiator(
    Sensor=antenna, # антенный элемент
    OperatingFrequency=fc # несущая частота
);

Setting up the target object

Create a Target object corresponding to a real reflective target with a given EPR. The reflections from this target will simulate real reflected signals. To calculate false alarms, create a second target with zero EPR. The reflections from this target are zero, except for noise.

In [ ]: 
target = [EngeePhased.RadarTarget() EngeePhased.RadarTarget()] # вектор целей
targetplatform = [ EngeePhased.Platform()  EngeePhased.Platform()] # вектор моделей движения целей
# задание параметров целей
target[1] = EngeePhased.RadarTarget(
    MeanRCS=tgt_rcs, # ЭПР 
    OperatingFrequency=fc # несущая частота
);
target[2] = EngeePhased.RadarTarget(
    MeanRCS=1e-16,
    OperatingFrequency=fc
)
# задание параметров движения целей
targetplatform[1] = EngeePhased.Platform(
    InitialPosition=[target_range; 0; 0] # начальное положение
)

targetplatform[2] = EngeePhased.Platform(
    InitialPosition=[target_range; 0; 0] # начальное положение
);

Setting the parameters of the propagation medium object

Let's set the signal propagation conditions - propagation channels from the radar to the targets and back.

In [ ]: 
channel = [EngeePhased.FreeSpace() EngeePhased.FreeSpace()]
channel[1] = EngeePhased.FreeSpace(
    SampleRate=fs, # частота дискретизации
    TwoWayPropagation=true, # модель двунаправленного распространения
    OperatingFrequency=fc # несущая частота
);
channel[2] = EngeePhased.FreeSpace(
    SampleRate=fs, # частота дискретизации
    TwoWayPropagation=true,  # модель двунаправленного распространения
    OperatingFrequency=fc # несущая частота
);

Setting the parameters of the transmitter system object

Let's generate noise by setting the NoiseMethod property to "Noise temperature " and the ReferenceTemperature property to 290 K.

In [ ]: 
# Геометрия приемной антенны
collector = EngeePhased.Collector(
    Sensor=antenna, # антенна
    OperatingFrequency=fc # Несущая частота
);
# формирование передатчика
receiver = EngeePhased.ReceiverPreamp(
    Gain=txrx_gain, # усиление
    NoiseMethod="Noise temperature", # источник шумов - тепловые шумы
    ReferenceTemperature=290.0, # шумовая температура
    NoiseFigure=0, # фактор шума
    SampleRate=fs, # частота дискретизации
    EnableInputPort=true # синхронизация с передатчиком
);
# Задание начального значения источника шума
receiver.SeedSource = "Property";
receiver.Seed = 2010;

Setting the fast time grid

A fast time grid is a set of time samples within the pulse duration. Each sample corresponds to an element of range resolution.

In [ ]: 
fast_time_grid = unigrid(0,1/fs,1/prf,"[)"); # сетка быстрого времени
rangebins = c*fast_time_grid/2; # элементы разрешения по дальности

Probing signal generation

Generate pulses of the probing signal you want to transmit.

In [ ]: 
wavfrm = waveform(); # формирования ЛЧМ-сигнала

Next, create a transmitted signal that includes the gain of the transmitted antenna.

In [ ]: 
sigtrans,tx_status = transmitter(wavfrm);

Create matched filter coefficients from the waveform object. Then create a system object of matched filter.

In [ ]: 
MFCoeff = getMatchedFilter(waveform); # коэффициенты фильтра
matchingdelay = size(MFCoeff,1) - 1; # заддержка на согласованную фильтрацию
mfilter = EngeePhased.MatchedFilter( # согласованный фильтр
    Coefficients=MFCoeff,
    GainOutputPort=false 
);

Calculating the range to the target

Let's calculate the range to the target, and then calculate the element number in the range resolution vector. Since the target and radar are stationary, the position and velocity values are constant throughout the simulation. We can assume that the number in the resolution vector is constant throughout the simulation.

In [ ]: 
# параметры движения радара и целей
ant_pos = antennaplatform.InitialPosition;
ant_vel = antennaplatform.Velocity;
tgt_pos = targetplatform[1].InitialPosition;
tgt_vel = targetplatform[1].Velocity;
# Расчет дальности и азимутальный угол от радара до целей
tgt_rng,tgt_ang = rangeangle(tgt_pos,ant_pos);
# поиск номера элемента разрешения по дальности
rangeidx = val2ind(tgt_rng,rangebins[2]-rangebins[1],rangebins[1])[1];

Pulse overlap

Let's create a signal processing cycle. Each step is executed by executing system objects. The processing cycle is called twice: once for the condition "target present" and once for the condition "target absent".

  1. Radiate the signal into space using EngeePhased.Radiator.
  2. Propagate the signal to the target and back to the antenna using EngeePhased.FreeSpace
  3. Reflecting the signal away from the target using EngeePhased.Target.
  4. Passing the received signal through the receive amplifier with EngeePhased.ReceiverPreamp.This step also adds random noise to the signal.
  5. Receiving the reflected signal at the antenna using EngeePhased.Collector.
  6. Matched filter the amplified signal using EngeePhased.MatchedFilter.
  7. Save the matched filter output in the target band bin index for further analysis.
In [ ]: 
rcv_pulses = zeros(ComplexF64,length(sigtrans),Npulsebuffsize);
h1 = zeros(ComplexF64,Npulse,1);
h0 = zeros(ComplexF64,Npulse,1);
Nbuff = floor(Int64,Npulse/Npulsebuffsize);
Nrem = Npulse - Nbuff*Npulsebuffsize;
for n = 1:2 # H1 and H0 Hypothesis
    trsig = radiator(sigtrans,tgt_ang);
    trsig = channel[n](trsig,
        ant_pos,tgt_pos,
        ant_vel,tgt_vel);
    global trsig  = trsig
    rcvsig = target[n](trsig);
    rcvsig = collector(rcvsig,tgt_ang);
    
    for k = 1:Nbuff
        for m = 1:Npulsebuffsize
            rcv_pulses[:,m] = receiver(rcvsig, .!(tx_status.>0));
        end
        rcv_pulses = mfilter(rcv_pulses);
        rcv_pulses = buffer(rcv_pulses[(matchingdelay+1):end][:,:],size(rcv_pulses,1));
        if n == 1
            h1[(1:Npulsebuffsize) .+ (k-1)*Npulsebuffsize] = transpose(rcv_pulses[rangeidx,:]);
        else
            h0[(1:Npulsebuffsize) .+ (k-1)*Npulsebuffsize] = transpose(rcv_pulses[rangeidx,:]);
        end
    end
    if (Nrem > 0)
        for m = 1:Nrem
            rcv_pulses[:,m] = receiver(rcvsig, .!(tx_status.>0));
        end
        rcv_pulses = mfilter(rcv_pulses);
        rcv_pulses = buffer(rcv_pulses[(matchingdelay+1):end][:,:],size(rcv_pulses,1));
         if n == 1
            h1[(1:Nrem) .+ Nbuff*Npulsebuffsize] = transpose(rcv_pulses[rangeidx,1:Nrem]);
        else
            h0[(1:Nrem) .+ Nbuff*Npulsebuffsize] = transpose(rcv_pulses[rangeidx,1:Nrem]);
         end
    end
end;

Building histograms of matched filter outputs

Let's calculate the histograms of the signals "target present" and "target absent". Use 100 samples to get a rough estimate of the scatter of the signals. Set the histogram limits from the smallest signal to the largest signal.

In [ ]: 
h1a = abs.(h1);
h0a = abs.(h0);
thresh_low = minimum([h1a;h0a]);
thresh_hi  = maximum([h1a;h0a]);
nbins = 100;
binedges = range(thresh_low,thresh_hi;length=nbins);
histogram(h0a,label="Target Absent",)
histogram!(h1a,label="Target Present")
title!("Target-Absent Vs Target-Present Histograms")
Out[0]:

Comparison of modelled and theoretical Pd and Pfa

To calculate Pd and Pfa, count the number of times the "target absent" return and the "target present" return exceed a given threshold. This set of thresholds has a finer granularity than the bins used to create the histogram in the previous simulation. We then normalise these counts by the number of pulses to obtain a probability estimate. The vector sim_pfa is the modelled false alarm probability as a function of threshold, thresh. The vector sim_pd is the modelled probability of detection, also a function of thresh. The receiver sets the threshold so that it can detect the presence or absence of a target. The histogram above shows that the best threshold is about 1.8.

In [ ]: 
nbins = 1000;
thresh_steps = range(thresh_low,thresh_hi;length=nbins)
sim_pd = zeros(1,nbins)
sim_pfa = zeros(1,nbins)
for k = 1:nbins
    thresh = thresh_steps[k]
    sim_pd[k] = sum(h1a .>= thresh)
    sim_pfa[k] = sum(h0a .>= thresh)
end
sim_pd = sim_pd./Npulse;
sim_pfa = sim_pfa./Npulse;

To plot the experimental ROC curve, you must invert the Pfa curve so that you can plot Pd versus Pfa. You can only invert the Pfa curve if you can express Pfa as a strictly monotonically decreasing function of Thresh. To express Pfa in this way, find all array indices for which Pfa is a constant over the neighbouring indices. Then remove these values from the arrays Pd and Pfa.

In [ ]: 
pfa_diff = diff(sim_pfa[:])
idx = (pfa_diff.!=0).*(1:length(pfa_diff))
sim_pfa = sim_pfa[idx[idx .!= 0]]
sim_pd = sim_pd[idx[idx .!= 0]];

Let's restrict the smallest value of Pfa to $10^{-6}$.

In [ ]: 
minpfa = 1e-6;
N = sum(sim_pfa .>= minpfa)
sim_pfa = reverse(sim_pfa[1:N])
sim_pd = reverse(sim_pd[1:N]);

Calculate the theoretical values of Pfa and Pd from the smallest Pfa to 1. Then plot the theoretical Pfa curve.

In [ ]: 
theor_pd,theor_pfa = rocsnr(snr_min,"Pd&Pfa",
    SignalType="NonfluctuatingNoncoherent",
    MinPfa=minpfa,NumPoints=N,NumPulses=1);
plot(theor_pfa,theor_pd,label="Theoretical",legend_position=:bottomright,
    legend_background_color =:transparent)
scatter!(sim_pfa,sim_pd,lab = "Simulated",xscale= :log10,ms=2,
    markerstrokecolor = :red)
title!("Simulated and Theoretical ROC Curves")
xlabel!("Pfa")
ylabel!("Pd")
Out[0]:

Improved modelling with one million pulses

In the previous simulation, the Pd values at low Pfa do not fall along a smooth curve or even reach the specified operating mode. The reason for this is that at very low Pfa levels, very few, if any, samples exceed the threshold value. To generate curves at low Pfa, it is necessary to use a number of samples of the order of the inverse of Pfa. This type of modelling is time consuming. The following curve uses a million pulses instead of 45,000. To run this simulation, repeat the example but set Npulse to 1000000.

image.png