Estimation of range and Doppler using pulse compression
This example demonstrates the pulse compression effect when the transmitted pulse correlates with the received signal. Radar and sonar systems use pulse compression to improve signal-to-noise ratio (SNR) and range resolution by reducing the duration of reflected signals. This example also demonstrates Doppler processing, in which the radial velocity of the target is determined by the Doppler shift that occurs when the target moves.
Determining the range and speed of targets
The radar system presented below transmits a pulsed sequence of linear frequency-modulated (LFM) waves to a moving target and receives echo signals. By applying consistent filtering and Doppler processing, the radar system can effectively determine the range and speed of a target.
Let's set the requirements for the radar system. In this example, the carrier frequency fc = 3 GHz and the sampling frequency fs = 1 MHz are used.
using Pkg; Pkg.add("DSP")
using FFTW,DSP
fc = 3e9; # carrier frequency, Hz
fs = 1e6; # sampling rate, Hz
c = 299_792_458; # m/s
Create system objects to simulate the radar system. The system is monostatic. The transmitter is located at (0,0,0) m and is stationary, and the target is located at (5000,500,0) m at a speed of (25,25,0) m/s and has an effective scattering area (ESR) of 1 sq. m.
antenna = EngeePhased.IsotropicAntennaElement(FrequencyRange=[1e8 10e9]); # antenna element
transmitter = EngeePhased.Transmitter(Gain=20,InUseOutputPort=true); # The transmitter
txloc = [0;0;0]; # position of the transmitter
tgtloc = [5000;5000;0]; # target coordinates, range ~= 7071 m
tgtvel = [25;25;0]; # Target speed, relative velocity ~= 35.4 m/s
target = EngeePhased.RadarTarget(Model="Nonfluctuating",MeanRCS=1,OperatingFrequency=fc);
antennaplatform = EngeePhased.Platform(InitialPosition=txloc); # the transmitter object in space
targetplatform = EngeePhased.Platform(InitialPosition=tgtloc,Velocity=tgtvel); # the target object in space
radiator = EngeePhased.Radiator( # The radiator
PropagationSpeed=c, # The speed of propagation
OperatingFrequency=fc, # carrier frequency of the signal
Sensor=antenna # antenna element
);
channel = EngeePhased.FreeSpace( # distribution channel
PropagationSpeed=c, # The speed of propagation
OperatingFrequency=fc, # carrier frequency
TwoWayPropagation=false # bi-directional direction
);
collector = EngeePhased.Collector(
PropagationSpeed=c, # The speed of propagation
OperatingFrequency=fc, # carrier frequency of the signal
Sensor=antenna # antenna element
);
receiver = EngeePhased.ReceiverPreamp( # The transmitter
NoiseFigure=0, # The noise factor
EnableInputPort=true, # input port for synchronization with the transmitter
SeedSource="Property", # the source of the initial value of the PSP generator for phase noise
Seed=2024 # the value of the PSP generator for phase noise
);
Create an FM signal with a pulse duration of 10 microseconds, a pulse repetition rate of 10 kHz, and a scan bandwidth of 100 kHz. Matched filter coefficients are generated based on this signal.
waveform = EngeePhased.LinearFMWaveform( # LFM signal generator
PulseWidth=10e-6, # pulse duration
PRF=10e3, # pulse repetition rate
OutputFormat="Pulses", # output type
NumPulses=1, # number of pulses
SweepBandwidth=1e5 # The growth time of the saw
);
wav = waveform(); # LFM signal generation
maxrange = c/(2*waveform.PRF); # Maximum one-digit target detection range
SNR = npwgnthresh(1e-6,1,"noncoherent"); # calculation of SNR based on the detection threshold
lambda = c/fc; # wavelength
tau = waveform.PulseWidth; # pulse duration
Ts = 290; # Receiver noise temperature
# Consistent filter
Mfilter = EngeePhased.MatchedFilter(
Coefficients=getMatchedFilter(waveform), # filter coefficients
GainOutputPort=true # transmission ratio output port
);
To improve the Doppler resolution, the system emits 64 pulses. The signal data will be stored in the variable rxsig, and the Data Matrix stores fast time samples (time within each pulse) in each column and slow time samples (time between pulses) in each row.
numPulses = 64; # number of pulses
rxsig = zeros(ComplexF64,length(wav),numPulses); # memory allocation for the output signal
for n = 1:numPulses
tgtloc,tgtvel = targetplatform(1/waveform.PRF); # goal movement parameters
tgtrng,tgtang = rangeangle(tgtloc,txloc); # radar motion parameters
txsig, txstatus = transmitter(wav); # formation of the transmitting signal
txsig = radiator(txsig,tgtang); # emission of an EM wave into space
txsig = channel(txsig,txloc,tgtloc,[0;0;0],tgtvel); # Channel: transmitter-target
txsig = target(txsig); # Reflection of the signal based on the target's EPR
txsig = channel(txsig,tgtloc,txloc,tgtvel,[0;0;0]); # Channel: target-receiver
txsig = collector(txsig,tgtang); # signal generation in the receiving antenna
rxsig[:,n] = receiver(txsig,xor.(txstatus,1)); # pre-amplification of the received signal
end
prf = waveform.PRF; # pulse repetition rate
fs = waveform.SampleRate; # sampling rate
response = EngeePhased.RangeDopplerResponse(
DopplerFFTLengthSource="Property", # setting the FFT length source
DopplerFFTLength=2048, # Doppler FFT length
SampleRate=fs, # sampling rate
DopplerOutput="Speed", # output type
OperatingFrequency=fc, # carrier frequency
PRFSource="Property", # the source of the HSI task
PRF=prf # pulse repetition rate
);
filt = getMatchedFilter(waveform); # coefficients of the matched filter
resp,rng_grid,dop_grid = response(rxsig,filt); # range and speed response
rng_grid = range(extrema(rng_grid)...,length = length(rng_grid)) # range grid
dop_grid = range(extrema(dop_grid)...,length = length(dop_grid)) # The speed grid
# building a dolero-range portrait
heatmap(dop_grid,rng_grid,DSP.pow2db.((abs.(resp)).^2),ylims=(0,12e3),
fontfamily=font(14,"Computer Modern"), xticks=(-200:50:200),
colorbar_title = "Power, dBW",
title=("Range-Doppler portrait"),
margin=5Plots.mm
)
xlabel!("Speed, m/s")
ylabel!("Range, m")
The response indicates that the target is moving at a speed of about -40 m/s, and that the target is moving away from the transmitter because the speed is negative.
We will calculate the range estimate corresponding to the propagation velocity of the signal. The graph of the first slow time sample shows the highest peak in the region of 7000 m, which corresponds to the graph of the response to the range speed.
fasttime = unigrid(0,1/fs,1/prf,"[)"); # The fast time grid
rangebins = (c*fasttime/2); # range resolution
# Plotting the dependence of amplitude on range
plot(rangebins.*1e-3,
abs.(rxsig[:,1].*1e6),
lab="",xlabel="Range, km",
ylabel="Amplitude, MV"
)
Range detection
Let's set a detection threshold at which the probability of a false alarm will be lower. . Let's perform an incoherent accumulation of 64 pulses, assuming that the signal is in white Gaussian noise. Based on the largest peak exceeding the threshold value, we will determine the target range estimate.
pfa = 1e-6; # the probability of a false alarm
NoiseBandwidth = 5e6/2; # noise band
# Calculation of noise power
npower = noisepow(
NoiseBandwidth, # Noise band
receiver.NoiseFigure, # The noise factor
receiver.ReferenceTemperature # Noise temperature
);
thresh = npwgnthresh(pfa,numPulses,"noncoherent"); # Detection threshold
thresh = npower*DSP.db2pow(thresh); # receiver sensitivity
pks,range_detect = findpeaks(pulsint(rxsig,"noncoherent"),"MinPeakHeight",thresh); # Peak detection
range_estimate = rangebins[range_detect[1]]; # determining the range estimate
print("Range: $(round(range_estimate;sigdigits=6)), m")
Speed detection
Let's select a range sample that includes the detected target, and plot the power spectral density estimate using the periodagram function.
out = DSP.periodogram(rxsig[range_detect[1],:][:]; onesided=false,
nfft=256, fs=prf, window=nothing)
spec_power = fftshift(DSP.pow2db.(out.power))
freq = fftshift(fftfreq(256,prf))
plot(freq.*1e-3,spec_power,lab="",gridalpha=0.15)
xlabel!("Frequency, kHz", fontfamily="Computer Modern")
ylabel!("Power, dBW", fontfamily="Computer Modern")
title!("Spectral power density", fontfamily="Computer Modern")
The peak frequency corresponds to the Doppler shift divided by 2, which can be converted to the target speed. A positive velocity means that the target is approaching the radiator, and a negative velocity means that the target is moving away from the radiator.
Y,I = maximum(spec_power),argmax(spec_power);
lambda = c/fc;
tgtspeed = dop2speed(freq[I]/2,lambda);
print("Doppler shift estimate: $(freq[I]/2) Hz")
print("Speed estimate: $(round(tgtspeed;sigdigits=4)) m/s")