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
# A function for plotting the radar output response and threshold value
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(["Impulse # $(i)" for i in 1:num_plot],1,num_plot)
xlabel_text = [reshape(repeat([""],1,num_plot-1),1,num_plot-1) "Time, from"]
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(["Power, dBW"],1,num_plot),label="",ls=:dash) #
end;
# Function for calculating noise power
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; # The probability of correct detection
pfa = 1e-6; # The probability of a false alarm
max_range = 5000; # [m] Maximum range
range_res = 50; # [m] Required range resolution
tgt_rcs = 1; # [m^2] EPR targets
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"); # [m/s] The speed of propagation
pulse_bw = prop_speed/(2*range_res); # [Hz] Frequency domain pulse band
pulse_width = 1/pulse_bw; # [c] Pulse duration
prf = prop_speed/(2*max_range); # [Hz] Pulse repetition rate
fs = 2*pulse_bw; # [Hz] Sampling rate of the system
waveform = EngeePhased.RectangularWaveform(PulseWidth=1/pulse_bw,PRF=prf,SampleRate=fs); # system object
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; # [Hz] receiver noise band
receiver = EngeePhased.ReceiverPreamp(
Gain=20, # [dB] gain
NoiseFigure=0,# [dB] noise factor
SampleRate=fs, # [Hz] sampling rate
EnableInputPort=true # a port for synchronizing the transmitter and receiver
);
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 the receiver's 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]; # SNR values in dB
rocsnr(snr_db,SignalType="NonfluctuatingNoncoherent") # building detection characteristics
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 the pulse accumulation technique 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) # calculation of SNR by solving the Albrshem equation
print("Minimum SNR: $(round(snr_min;sigdigits=5)) dB")
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; # [dB] Transmitter gain
fc = 10e9; # [Hz] carrier frequency of the signal
lambda = prop_speed/fc; # [m] signal wavelength
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) # peak power
print("Peak power: $(peak_power*1e3) kW")
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, # [dB] gain
PeakPower=peak_power, # [W] peak power
InUseOutputPort=true # enabling the output port that defines
# The status of the transmitter (1 - on, 0 - off)
);
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.
# radar antenna element
antenna = EngeePhased.IsotropicAntennaElement(
FrequencyRange=[5e9 15e9]); # [Hz] frequency range of the radiation antenna
# radar motion model
sensormotion = EngeePhased.Platform(
InitialPosition=[0; 0; 0], # [m] initial position of the antenna
Velocity=[0; 0; 0]); # [m/s] antenna speed
Using the antenna and the operating frequency, we set the radiator (Radiator) and the receiving antenna (Collector).
radiator = EngeePhased.Radiator(
Sensor=antenna, # setting the antenna element
OperatingFrequency=fc); # [Hz] setting the carrier frequency
collector = EngeePhased.Collector(
Sensor=antenna, # setting the antenna element
OperatingFrequency=fc); # [Hz] setting the carrier frequency
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]]; # [m] components of the initial position of 3 targets [[px] [py] [pz]]
tgtvel = [[0;0;0] [0;0;0] [0;0;0]]; # [m/s] components of the velocity vector of 3 targets [[vx] [vy] [vz]]
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # goal movement model
tgtrcs = [1.6 2.2 1.05]; # m^2 EPR goals
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc); # The goal system object
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, # [Hz] sampling rate
TwoWayPropagation=true, # bidirectional propagation
OperatingFrequency=fc); # [Hz] carrier frequency
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); # [c] fast time grid
slow_time_grid = (0:num_pulse_int-1)./prf; # [c] slow time grid
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) # access to updating object parameters
receiver.SeedSource = "Property"; # changing the type of noise generation model
receiver.Seed = 2007; # setting the initial state of the thermal noise generator
# Allocating memory for the resulting data
rxpulses = zeros(ComplexF64,length(fast_time_grid),num_pulse_int);
for m = 1:num_pulse_int
# Updating the radar location and targets
sensorpos,sensorvel = sensormotion(1/prf);
tgtpos,tgtvel = tgtmotion(1/prf);
# Calculating angles based on the range from radar to targets
global tgtrng,tgtang = rangeangle(tgtpos,sensorpos);
# Simulation of signal propagation from radar to targets
pulse = waveform(); # rectangular pulse generator
txsig,txstatus = transmitter(pulse); # transmitting signal amplifier
txsig = radiator(txsig,tgtang); # signal emission into the propagation medium
txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel); # accounting for the distribution environment
# Reflection of the signal from targets, taking into account the EPR
tgtsig = target(txsig);
# Reception and pre-amplification of the reflected signal
rxsig = collector(tgtsig,tgtang); # receiving antenna
rxpulses[:,m] = receiver(rxsig,xor.(txstatus.>0,1)); # The preamp
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); # [W] noise power
threshold = npower * DSP.db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent")); # detection threshold
We plot the first two received pulses with a threshold
num_pulse_plot = 2; # number of graphs
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, which 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); # calculation of the coefficients of the agreed
# a filter based on a refined signal
matchedfilter = EngeePhased.MatchedFilter(
Coefficients=matchingcoeff, # filter coefficients
GainOutputPort=true); # gain output port
rxpulses, mfgain = matchedfilter(rxpulses); # response after consistent filtering
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 is accumulated continuously, there will be no end to the array.
matchingdelay = size(matchingcoeff,1)-1; # the length of the coefficients of the matched filter
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; # [m] range resolution
tvg = EngeePhased.TimeVaryingGain( # creating a VAR system object
RangeLoss=2*fspl.(range_gates,lambda), # [dB] range-adjusted losses
ReferenceLoss=2*fspl.(max_range,lambda)); # [dB] relative loss at maximum range
rxpulses = tvg(rxpulses); # the output signal after operation is WARU
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 integration of pulses that allows us to use a sequence of pulses of lower power.
Incoherent accumulation
We can further improve the SNR by incoherently integrating (video integrating) the received pulses.
rxpulses = pulsint(rxpulses,"noncoherent"); # incoherent accumulation
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)); # Peak detection
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 target range: $(true_range) m")
print("Target range estimate: $(range_estimates) m")
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.