Simulation 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 the target. A monostatic radar has a transmitter located next to the receiver. The transmitter generates a pulse that hits the target and causes an echo to be received by the receiver. By measuring the location of the echo signals over time, it is possible to estimate the range to the target.
This example is devoted to the development of a pulse radar system that can achieve specified technical characteristics. 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 for receiving the reflected signal are also simulated. Finally, signal processing methods are applied to the received signal to determine the range to the targets.
Functions used
Pkg.add(["DSP"])
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;
Technical characteristics of the device
The purpose of developing this pulse radar system is to detect non-fluctuating targets with a radar cross-section (RCS) of at least 1 at a distance of up to 5,000 meters from the radar with a range resolution of 50 meters. The desired performance indicator is 0.9 detection probability (Pd) and less false alarm probability (Pfa). . Since coherent detection requires phase information and is therefore computationally expensive, we use an incoherent detection scheme. In addition, this example assumes that the environment is in free space.
pd = 0.9; # Вероятность верного обнаружения
pfa = 1e-6; # Вероятность ложной тревоги
max_range = 5000; # [м] Максимальная дальность
range_res = 50; # [м] Требуемое разрешение по дальности
tgt_rcs = 1; # [м^2] ЭПР цели
Designing a monostatic radar system
We need to identify several characteristics of the radar system, such as the waveform, receiver, transmitter, and antenna used to emit and receive the signal.
The 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 the pulsed waveform is the pulse repetition rate (PRF). The PRF is determined by the maximum single-digit range.
prop_speed = physconst("lightspeed"); # [м/с] Скорость распространения
pulse_bw = prop_speed/(2*range_res); # [Гц] Полоса импульса в частотной области
pulse_width = 1/pulse_bw; # [c] Длительность импульса
prf = prop_speed/(2*max_range); # [Гц] Частота следования импульсов
fs = 2*pulse_bw; # [Гц] Частота дискретизации системы
waveform = EngeePhased.RectangularWaveform(PulseWidth=1/pulse_bw,PRF=prf,SampleRate=fs); # системый объект
Please note that we have set the sampling rate to double the bandwidth.
Receiver noise characteristics
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 bandwidth of the receiver. The noise bandwidth of the receiver 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 factor of 0 dB.
noise_bw = pulse_bw; # [Гц] шумовая полоса приемника
receiver = EngeePhased.ReceiverPreamp(
Gain=20, # [дБ] усиление
NoiseFigure=0,# [дБ] фактор шума
SampleRate=fs, # [Гц] частота дискретизации
EnableInputPort=true # порт для синхронизации передатчика и приемника
);
Please note that since we are simulating a monostatic radar, the receiver cannot be turned on until the transmitter is turned off. Therefore, we set the flag EnableInputPort=true so that the synchronization signal can be transmitted from the transmitter to the receiver.
The transmitter
The most important parameter of the transmitter 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 purpose of the Pd and Pfa, as well as the detection scheme implemented in the receiver.
The relationship between Pd, Pfa, and SNR can best be represented by a receiver operational performance curve (ROC). We can construct a curve where Pd is a function of Pfa for various 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, Pfa = 1e-6 and Pd = 0.9; the SNR of the received signal must exceed 13 dB. This is a rather high requirement and not very practical. To make the radar system more feasible, we can use pulse accumulation techniques to reduce the required SNR. If we choose to integrate 10 pulses, then the curve can be constructed 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 reduction of SNR can be achieved by integrating more pulses, but the number of pulses available for integration is usually limited due to target movement or environmental heterogeneity.
In the above approach, the SNR value is read from the curve, but it is often desirable to calculate only the required value. For an incoherent detection scheme, calculating the required SNR is theoretically quite difficult. Fortunately, there are good approximations, 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)) дБ")
Once you have the desired SNR value at the receiver, you can calculate the peak power at the transmitter using the radar equation. Here we assume that the gain of the transmitter is 20 dB.
To calculate peak power using the 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) кВт")
Please note that the resulting power is about 5 kW, which is an acceptable value. For comparison, if we had not used the pulse integration method, the peak power would have been 33 kW, which is quite a large value.
After receiving all this information, we can set up the transmitter.
transmitter = EngeePhased.Transmitter(
Gain=tx_gain, # [дБ] усиление
PeakPower=peak_power, # [Вт] пиковая мощность
InUseOutputPort=true # включение выходного порта, который определяет
# состояние передатчика (1 - вкл, 0 - выкл)
);
Since this example simulates a monostatic radar system, the InUseOutputPort port is set to true to output the status of the transmitter. This status signal can be used to turn on the receiver.
Transmitting and receiving devices
In a radar system, the signal propagates as an electromagnetic wave. Therefore, the signal must be emitted and collected by the antenna used in the radar system. This is where the radiating and "collecting" antenna come into play.
In a monostatic radar system, the emitter and collector use the same antenna, so first we will define the antenna. To simplify the design, we chose an isotropic antenna. Please note that the antenna must operate at the operating frequency of the system (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]); # [м/с] скорость антенны
Using the antenna and the operating frequency, we set 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 necessary for modeling. Then we will simulate the return of the signal and perform a range determination based on the simulated signal.
System modeling
Goals
To test our radar's ability to detect targets, we must first identify the targets. Let's assume that there are 3 stationary, non-fluctuating targets in space. Their positions and radar cross-sections are shown 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]; # м^2 ЭПР целей
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc); # Системный объект целей
Distribution environment
To simulate the signal, we also need to determine 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 configured to simulate propagation delays in two directions.
Signal synthesis
Now we are ready to simulate the entire system.
The synthesized signal is a data matrix with fast time (time inside each pulse) in each column and slow time (time between pulses) in each row. To visualize a 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; # [c] сетка медленного времени
The next cycle simulates 10 pulses of the received signal.
Let's set the initial state of the noise generator in the receiver so that the same results can be reproduced.
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
Range detection
Threshold detection
The detector compares the signal strength with a preset threshold. In radar applications, the threshold is often chosen so that the 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 to integrate the pulses, the threshold signal power is determined 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 with a threshold
num_pulse_plot = 2; # количество графиков
RadarPulsePlot(rxpulses,threshold,
fast_time_grid,slow_time_grid,num_pulse_plot)
The threshold in these drawings is used for demonstration purposes only. Note that the signals from the second and third targets are much weaker than those from the first, as they are further away from the radar. Therefore, the power of the received signal depends on the range, and the threshold is unfair to targets located at different distances.
Consistent filtering
A matched filter provides enhanced processing that increases the detection threshold. It converts the received signal into a local, time-reversed and conjugated copy of the transmitted oscillator type. Therefore, when creating a matched filter, we must specify the type of probing signal. The received pulses are first passed through a 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 no longer matches the true location of the target. 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 data accumulates continuously, there will be no end to the array.
matchingdelay = size(matchingcoeff,1)-1; # длина коэффициентов согласованного фильтра
rxpulses = buffer(rxpulses[matchingdelay+1:end],200); # size(rxpulses,1) ->200
The threshold will then increase by the processing gain of the matched filter.
threshold = threshold * db2pow(mfgain)
The following graph shows the same two pulses after passing through a 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 power of the received signal depends on the range, the signal from a nearby target is still much stronger than from a target located at a greater distance. Therefore, as shown in the figure above, noise from a nearby bean has a significant chance of exceeding the threshold and overshadowing the target further away. To ensure that the threshold value is valid for all targets within the detectable range, we can use a time-varying gain factor to compensate for range-dependent losses in the received echo signal.
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. After receiving this information, we apply a time-varying gain to the received pulse so that the returns appear to be 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 let's plot the same two pulses after normalization of the range.
RadarPulsePlot(rxpulses,threshold,fast_time_grid,slow_time_grid,num_pulse_plot)
Changing the gain over time leads to an increase in the noise level. However, the return of the target is now independent of the 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. We need to integrate the pulses so that the power of the returned echoes from the targets exceeds the threshold, and the noise level remains below the bar. This is to be expected, since it is the 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 stage, the data is ready for the final detection stage. The figure shows that all three echo signals from targets exceed the threshold value, which means they can be detected.
Range detection
Finally, threshold detection is performed for integrated pulses. The detection circuit identifies 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) м")
Please 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 example, we have designed a radar system based on a set of set targets. Based on these goals, many design parameters of the radar system were calculated. The example also shows how to use the developed radar to perform the range detection task. In this example, the radar used a rectangular waveform.