Modelling of test signals for a radar receiver

This example shows how to simulate the received signal of a monostatic pulse radar to estimate the range to a target. The monostatic radar has a transmitter located near the receiver. The transmitter generates a pulse that hits the target and causes an echo that is received by the receiver. By measuring the location of the echoes over time, it is possible to estimate the range to the target.

This case study focuses on the design of a pulsed radar system that can achieve a given specification. It describes the steps to translate design specifications, such as detection probability and range resolution, into radar system parameters such as transmitter power and pulse duration. The environment and targets are also modelled to obtain the reflected signal. Finally, signal processing techniques are applied to the received signal to determine the range to the targets.

Functions used

Pkg.add(["DSP"])

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

using DSP

# Функция для построения выходного отклика радара и порогового значения
function RadarPulsePlot(ReceivePulse,threshold,t,tp,num_plot::Int64)
    thresh = sqrt(threshold)*ones(length(t),num_plot)
    ReceivePulse =  Float64.(pow2db.((abs.(ReceivePulse)).^2))
    ReceivePulse[ReceivePulse.==0] .= eps(0.)
    thresh .=  pow2db.(thresh.^2)
    t = repeat(t,1,num_plot) .+ tp[1:num_plot]'
    text_title = reshape(["Импульс № $(i)" for i in 1:num_plot],1,num_plot)
    xlabel_text = [reshape(repeat([""],1,num_plot-1),1,num_plot-1) "Время, с"]
    plot(t,ReceivePulse[:,1:num_plot],layout=(num_plot,1),label="",title = text_title,
        fontfamily="Computer Modern", titlefontsize=10, guidefontsize=10, tickfontsize=10,gridalpha=0.3,lw=1.2)
    plot!(t,thresh[:,1:num_plot],layout=(num_plot,1),xlabel=xlabel_text,lw=1.2,
        ylabel = repeat(["Мощность, дБВт"],1,num_plot),label="",ls=:dash) # 
end;

# Функция для расчета мощности шума
function noisepow1(nbw::Number, nf::Number = 0, reftemp::Number = 290)
    function systemp(nf::Number,reftemp::Number = 290)
        nf >= 0 || error("Invalid input")
        nfLinear = db2pow(nf); 
        stemp = reftemp*(nfLinear); # Kelvin
        return stemp
    end
    B = 1.380649e-23
    npow =  B * systemp(nf,reftemp) * nbw
    return npow 
end;

Device specifications

The design objective of this pulsed radar system is to detect non-fluctuating targets with a radar cross section (RCS) of at least one square metre at distances up to 5000 metres from the radar with a range resolution of 50 metres. The desired performance metric is a probability of detection (Pd) of 0.9 and a probability of false alarm (Pfa) of less than 1e-6. Since coherent detection requires phase information and is therefore computationally expensive, we use a non-coherent detection scheme. In addition, this example assumes that the medium is in free space.

pd = 0.9; # Вероятность верного обнаружения
pfa = 1e-6; # Вероятность ложной тревоги
max_range = 5000; # Максимальная дальность
range_res = 50; # Требуемое разрешение по дальности
tgt_rcs = 1; # ЭПР цели

Designing a monostatic radar system

We need to define several characteristics of the radar system such as waveform, receiver, transmitter and antenna used to emit and receive the signal.

Generator

In this example, we have chosen a rectangular waveform. The desired range resolution determines the bandwidth of the waveform, which, in the case of a rectangular waveform, determines the pulse width.

Another important parameter of a pulsed waveform is the pulse repetition frequency (PRF). The PRF is determined by the maximum single digit range.

prop_speed = 299_792_458; # Скорость распространения
pulse_bw = prop_speed/(2*range_res);  # Полоса импульса в частотной области
pulse_width = 1/pulse_bw; # Длительность импульса
prf = prop_speed/(2*max_range); # Частота следования импульсов
fs = 2*pulse_bw; # Частота дискретизации системы
waveform = EngeePhased.RectangularWaveform(PulseWidth=1/pulse_bw,PRF=prf,SampleRate=fs); # системый объект

Note that we have set the sampling rate equal to twice the bandwidth.

Noise characteristics of the receiver

We assume that the only noise present in the receiver is thermal noise, so there is no interference in this simulation. The power of thermal noise depends on the receiver bandwidth. The receiver noise bandwidth is set equal to the bandwidth of the waveform. This is often the case in real systems. We also assume that the receiver has a gain of 20 dB and a noise figure of 0 dB.

noise_bw = pulse_bw; # шумовая полоса приемника
receiver = EngeePhased.ReceiverPreamp(
    Gain=20, # усиление 
    NoiseFigure=0,# фактор шума
    SampleRate=fs, # частота дискретизации
    EnableInputPort=true # порт для синхронизации передатчика и приемника
);

Note that since we are modelling a monostatic radar, the receiver cannot be switched on until the transmitter is switched off. Therefore, we set the EnableInputPort property to true so that the synchronisation signal can be sent from the transmitter to the receiver.

Transmitter

The most important transmitter parameter is the peak power. The required peak power is related to many factors, including the maximum single-digit range, the required SNR in the receiver, and the pulse width. Among these factors, the required SNR in the receiver is determined by the design objective Pd and Pfa and the detection scheme implemented in the receiver.

The relationship between Pd, Pfa and SNR can be best represented by the receiver operating characteristic (ROC) curve. We can plot the curve, where Pd is a function of Pfa for different SNRs, using the following command:

snr_db = [-Inf, 0.0, 3, 10, 13]; # значения ОСШ в дБ
rocsnr(snr_db,SignalType="NonfluctuatingNoncoherent") # построение характеристик обнаружения 

The ROC curves show that to meet the design goals of Pfa = 1e-6 and Pd = 0.9; the SNR of the received signal must exceed 13 dB. This is quite a high requirement and not very practical. To make the radiolocation system more feasible, we can use pulse accumulation technique to reduce the required SNR. If we choose to integrate 10 pulses, the curve can be plotted as follows:

num_pulse_int = 10;
rocsnr([0 3 5],SignalType="NonfluctuatingNoncoherent",NumPulses=num_pulse_int)

We can see that the required power has decreased to about 5 dB. Further SNR reduction can be achieved by integrating more pulses, but the number of pulses available for integration is usually limited due to target motion or environmental inhomogeneity.

In the above approach, the SNR value is read from the curve, but it is often desirable to calculate only the desired value. For an incoherent detection scheme, calculating the desired SNR is theoretically quite difficult. Fortunately, good approximations exist, such as the Albersham equation. Using the Albersham equation, the required SNR value can be derived as

snr_min = albersheim(pd, pfa,N= num_pulse_int) # расчет ОСШ путем решения уравнения Альбршема
print("Минимальное ОСШ: $(round(snr_min;sigdigits=5)) дБ")

Минимальное ОСШ: 4.9904 дБ

Having obtained the required SNR value at the receiver, we can calculate the peak power at the transmitter using the radar equation. Here we assume that the transmitter gain is 20 dB.

To calculate the peak power using radar equation, we also need to know the wavelength of the propagating signal which is related to the operating frequency of the system. Here we set the operating frequency to 10 GHz.

tx_gain = 20; # Усиление передатчика
fc = 10e9; # несущая частота сигнала
lambda = prop_speed/fc; # длина волны сигнала
peak_power = ((4*pi)^3*noisepow1(1/pulse_width)*max_range^4*
    DSP.db2pow(snr_min))/(DSP.db2pow(2*tx_gain)*tgt_rcs*lambda^2) # пиковая мощность
print("Пиковая мощность: $(peak_power*1e3) кВт")

Пиковая мощность: 5.22647386447291e6 кВт

Note that the resulting power is about 5 kW, which is an acceptable value. In comparison, if we did not use the pulse integration method, the peak power would be 33 kW, which is quite a large value.

Having obtained all this information, we can configure the transmitter.

transmitter = EngeePhased.Transmitter(
    Gain=tx_gain, # усиление
    PeakPower=peak_power, # пиковая мощность
    InUseOutputPort=true # включение выходного порта, который определяет 
                         # состояние передатчика (1 - вкл, 0 - выкл)
);

Since this example is modelling a monostatic radar system, the InUseOutputPort is set to true to output the transmitter status. This status signal can be used to switch on the receiver.

Transmitter and Receiver

In a radar system, the signal propagates as an electromagnetic wave. The signal must therefore be radiated and collected by the antenna used in the radar system. This is where the radiating antenna and the "collecting" antenna come into play.

In a monostatic radar system, the transmitter and collector use the same antenna, so we will first define the antenna. To simplify the design, we have chosen an isotropic antenna. Note that the antenna must operate at the system operating frequency (10 GHz), so we will set the antenna frequency range to 5-15 GHz.

We assume that the antenna is stationary.

# антенный элемент радара
antenna = EngeePhased.IsotropicAntennaElement(
    FrequencyRange=[5e9 15e9]); # диапазон частот антенны излучения
# модель движения радара
sensormotion = EngeePhased.Platform(
    InitialPosition=[0; 0; 0], # начальное положение антенны
    Velocity=[0; 0; 0]); # скорость антенны

With antenna and operating frequency we define the radiator (Radiator) and the receiving antenna (Collector).

radiator = EngeePhased.Radiator(
    Sensor=antenna,             # задание антенного элемента 
    OperatingFrequency=fc);   # задание несущей частоты

collector = EngeePhased.Collector(
    Sensor=antenna,             # задание  антенного элемента 
    OperatingFrequency=fc);     # задание несущей частоты

This completes the configuration of the radar system. In the following sections we will define other objects such as the target and the environment that are needed for the simulation. Then we will simulate the signal return and perform range determination from the simulated signal.

System modelling

Objectives

To test our radar's ability to detect targets, we must first define targets. Assume that there are 3 stationary, non-fluctuating targets in space. Their positions and radar cross sections are given below.

tgtpos = [[2024.66;0;0] [3518.63;0;0] [3845.04;0;0]]; # компоненты начального положения 3 целей [[px] [py] [pz]]
tgtvel = [[0;0;0] [0;0;0] [0;0;0]]; # компоненты вектора скорости 3 целей [[vx] [vy] [vz]]
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # модель движения целей

tgtrcs = [1.6 2.2 1.05]; # ЭПР целей
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc); # Системный объект целей

Propagation environment

To model the signal, we also need to define the propagation channel between the radar system and each target.

channel = EngeePhased.FreeSpace(
    SampleRate=fs, # частота дискретизации
    TwoWayPropagation=true, # двунаправленное распространение
    OperatingFrequency=fc); # несущая частота

Since this example uses a monostatic radar system, the channels are set up to simulate propagation delays in two directions.

Signal synthesis

Now we are ready to model the whole system.

The synthesised signal is a data matrix with fast time (time within each pulse) in each column and slow time (time between pulses) in each row. To visualise the signal, it is useful to define both a fast time grid and a slow time grid.

fast_time_grid = (0:1/fs:1/prf-1/fs); # сетка быстрого времени
slow_time_grid = (0:num_pulse_int-1)./prf; # сетка медленного времени

The next cycle simulates 10 pulses of the received signal.

We set the primer to generate noise in the receiver so that we can reproduce the same results.

release!(receiver) # доступ к обновлению параметров объекта
receiver.SeedSource = "Property";  # изменение тип модели генерации шума
receiver.Seed = 2007; # задание начального состояния генератора тепловых шумов

# Выделение памяти для результирующих данных
rxpulses = zeros(ComplexF64,length(fast_time_grid),num_pulse_int);

for m = 1:num_pulse_int
    
    # Обновление пложения радара и целей
    sensorpos,sensorvel = sensormotion(1/prf);
    tgtpos,tgtvel = tgtmotion(1/prf);

    # Расчет углов по дальности от радара до целей  
    global tgtrng,tgtang = rangeangle(tgtpos,sensorpos);
    
    # Моделирование распространения сигнала от радара до целей
    pulse = waveform(); # генератор прямоугольных импульсво
    txsig,txstatus = transmitter(pulse); # усилитель передающего сигнала
    txsig = radiator(txsig,tgtang); # излучатель сгенерированного сигнала в среду
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel); # среда распространения
    
    # Отражение сигнала от целей с учетом ЭПР
    tgtsig = target(txsig);
    
    # Прием и предварительного усиление отраженного сигнала
    rxsig = collector(tgtsig,tgtang); # приемная антенна
    rxpulses[:,m] = receiver(rxsig,xor.(txstatus.>0,1)); # предусилитель
end

Determining the range

Detection threshold

The detector compares the signal strength to a predetermined threshold. In radar applications, the threshold is often chosen so that Pfa is below a certain level. In this case, we assume that the noise is white Gaussian and the detection is incoherent. Since we also use 10 pulses for pulse integration, the threshold signal power is defined as follows

npower = noisepow1(noise_bw,receiver.NoiseFigure,
    receiver.ReferenceTemperature);  # мощность шума
threshold = npower * DSP.db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent")); # порог обнаружения

We plot the first two received pulses against the threshold

num_pulse_plot = 2; # количество графиков
RadarPulsePlot(rxpulses,threshold,
    fast_time_grid,slow_time_grid,num_pulse_plot)

The threshold in these figures is used for demonstration purposes only. Note that the signals from the second and third targets are much weaker than the first target because they are further away from the radar. Therefore, the received signal strength is range dependent and the threshold is unfair to targets at different distances.

Matched filtering

A matched filter provides processing gain that increases the detection threshold. It converts the received signal into a local, time reversed and contiguous copy of the transmitted waveform. Therefore, while creating a matched filter, we must specify the transmitted waveform. The received pulses are first passed through the matched filter to improve the SNR before performing pulse integration, threshold detection, etc.

matchingcoeff = getMatchedFilter(waveform); # расчет коэффициентов согласованного 
                                            # фильтра на основе изулченного сигнала          
matchedfilter = EngeePhased.MatchedFilter(
    Coefficients=matchingcoeff, # коэффициенты фильтра
    GainOutputPort=true); # выходной порт усиления
rxpulses, mfgain = matchedfilter(rxpulses); # отклик после согласованной фильтрации

The matched filter introduces an internal delay, so the location of the peak (the output sample with the maximum SNR) no longer matches the true target location. To compensate for this delay, we move the output of the matched filter forward and add zeros at the end. Note that in real systems, since the data is accumulated continuously, there will be no end of the array.

matchingdelay = size(matchingcoeff,1)-1; # длина коэффициентов согласованного фильтра
rxpulses = buffer(rxpulses[matchingdelay+1:end],200); # size(rxpulses,1) ->200

The threshold will then be increased by the processing gain of the matched filter.

threshold = threshold * db2pow(mfgain);

The next graph shows the same two pulses after passing through the matched filter.

RadarPulsePlot(rxpulses,threshold,
    fast_time_grid,slow_time_grid,num_pulse_plot)

After the matched filter stage, the SNR improves. However, since the received signal power is range dependent, the signal from a close target is still much stronger than that from a target that is farther away. Therefore, as shown in the figure above, the noise from a close bin has a significant chance of exceeding the threshold and overshadowing a target further away. To ensure that the threshold is valid for all targets within the detectable range, we can use time-varying gain to compensate for range-dependent losses in the received echo.

To compensate for range-dependent losses, we first calculate the strobe pulse corresponding to each signal sample and then calculate the free-space loss corresponding to each strobe pulse . Once this information is obtained, we apply a time-varying gain to the received pulse so that the returns are as if from the same reference range (maximum detectable range).

range_gates = prop_speed*fast_time_grid/2; # разрешение по дальности
tvg = EngeePhased.TimeVaryingGain( # создание системного объекта ВАРУ
    RangeLoss=2*fspl.(range_gates,lambda), # потери с учетом дальности
    ReferenceLoss=2*fspl.(max_range,lambda)); # относительные потери при максимальной дальности
rxpulses = tvg(rxpulses); # выходной сигнал после работы ВАРУ

Now plot the same two pulses after range normalisation

RadarPulsePlot(rxpulses,threshold,fast_time_grid,slow_time_grid,num_pulse_plot)

The change in gain over time results in an increase in noise. However, the target return is now independent of range. A constant threshold can now be used for detection over the entire detectable range.

Note that at this stage the threshold is above the maximum power level contained in each pulse. Therefore, nothing can be detected at this stage yet. We need to perform pulse integration so that the power of the returned echoes from the targets exceeds the threshold and the noise level remains below the bar. This is expected because it is pulse integration that allows us to use a sequence of lower power pulses.

Incoherent accumulation

We can further improve the SNR by incoherently integrating (video-integrating) the received pulses.

rxpulses = pulsint(rxpulses,"noncoherent"); # некогерентное накопление
RadarPulsePlot(rxpulses,threshold,fast_time_grid,slow_time_grid,1)

After the video integration step, the data is ready for the final detection step. In the figure, we can see that all three target echoes exceed the threshold value, and thus can be detected.

Range Detection

Finally, threshold detection is performed for integrated pulses. The detection circuitry identifies the peaks and then translates their position into target ranges.

_,range_detect = findpeaks(rxpulses,"MinPeakHeight",sqrt(threshold)); # обнаружение пиков

The true and detected target ranges are shown below:

true_range = round.(tgtrng;sigdigits=4)
range_estimates = round.(range_gates[range_detect];sigdigits=4)
println("Истинная дальность целей: $(true_range) м")
print("Оценка дальности целей: $(range_estimates) м")

Истинная дальность целей: [2025.0 3519.0 3845.0] м
Оценка дальности целей: [2025.0, 3525.0, 3850.0] м

Note that these range estimates are accurate only up to the range resolution (50 m) that can be achieved by the radar system.

Conclusion

In this case study, we designed a radar system based on a set of specified targets. Based on these targets, many design parameters of the radar system were calculated. The example also shows how to use the designed radar to perform a range detection task. In this example, the radar used a rectangular waveform.