Doppler estimation
The example shows a monostatic pulse radar that detects the radial velocity of moving targets at certain ranges. The speed is determined by the Doppler shift caused by moving targets. First, we will determine the presence of a target at a given range, and then, using Doppler processing, we will estimate the radial velocity of the target at this range.
Used features and libraries
Pkg.add(["DSP"])
using DSP,FFTW
function BasicMonostaticRadarExampleData()
c = 299_792_458 # signal propagation speed
# Creating system objects:
# rectangular pulse generator
waveform = EngeePhased.RectangularWaveform(PulseWidth=100/c,
PRF=c/1e4,OutputFormat="Pulses",NumPulses=1,SampleRate= 2c/100);
# antenna element (type - isotropic)
antenna = EngeePhased.IsotropicAntennaElement(FrequencyRange=[5e9 1.5e10]);
# narrowband signal emitter
radiator = EngeePhased.Radiator(Sensor=antenna,
PropagationSpeed=c,OperatingFrequency=1e10);
# The amplifier in the transmission path
transmitter = EngeePhased.Transmitter(PeakPower=5226.4791642003665401716716587543,
Gain=20,LossFactor=0, InUseOutputPort=true,CoherentOnTransmit=true);
# receiving antenna
collector = EngeePhased.Collector(Sensor=antenna,
PropagationSpeed=c,Wavefront="Plane",
OperatingFrequency=1e10);
# The preamp is in the receiving path
receiver = EngeePhased.ReceiverPreamp(Gain=20,NoiseFigure=0,
ReferenceTemperature=290,SampleRate=2*c/100, PhaseNoiseInputPort=false,
EnableInputPort=true,SeedSource="Property",Seed=1000);
# radar motion model
sensormotion = EngeePhased.Platform(InitialPosition=[0;0;0],Velocity=[0;0;0]);
return waveform, antenna, radiator, collector, transmitter, receiver, sensormotion
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;
Configuring the radar system
First, we will define the radar system. Since this example focuses on Doppler processing, we use the radar system built in the example Simulation of test signals for radar приемника. Readers are invited to explore the design details of radar systems using this example.
waveform,antenna,radiator,collector,transmitter,receiver,sensormotion = BasicMonostaticRadarExampleData();
System modeling
Goals
Doppler processing uses the Doppler shift caused by a moving target. Now we will set three targets, specifying their position, effective scattering area (EPR) and velocities.
tgtpos = [[1200; 1600; 0] [3543.63; 0; 0] [1600; 0; 1200]]; # components of the initial position of goals
tgtvel = [[60; 80; 0] [0;0;0] [0; 100; 0]]; # speed components of goals
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # goal movement model
tgtrcs = [1.3 1.7 2.1]; # EPR goals
fc = radiator.OperatingFrequency;
target = EngeePhased.RadarTarget(MeanRCS=tgtrcs,OperatingFrequency=fc);
Note that the first target is located at a distance of 3,550, 1,600, and 1,200 m, and is moving at speeds of 60, 130, and 0 m/s. The difference is that the first target moves radially, while the second moves tangentially and radially. The third target is not moving.
Distribution environment
We also need to set up a distribution environment for each target. Since we are using a monostatic radar, we will define a two-way propagation model.:
fs = waveform.SampleRate; # sampling rate of the system
# distribution channel model
channel = EngeePhased.FreeSpace(
SampleRate=fs, # sampling rate of the system
TwoWayPropagation=true, # distribution model
OperatingFrequency=fc); # carrier frequency of the signal
Signal synthesis
Having determined the radar system, environment and targets, we will simulate the received echo signal reflected from the targets. The resulting data is a matrix of elements with fast time (time inside each pulse) in each column and slow time (time between pulses) in each row.
To reproduce the same results, we need to set the initial data for generating noise in the receiver.
release!(receiver) # resetting the internal state of an object (modification )
receiver.SeedSource = "Property"; # changing the type of noise generation model
receiver.Seed = 2024; # setting the initial state of the thermal noise generator
prf = waveform.PRF; # pulse repetition rate
num_pulse_int = 10; # number of accumulation pulses
fast_time_grid = (0:1/fs:1/prf-1/fs) # The fast time grid
slow_time_grid = (0:num_pulse_int-1)./prf; # slow time grid
# 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
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); # the emitter of the generated signal on Wednesday
txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel); # 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
Doppler estimation
Range detection
To estimate the Doppler shift of targets, you first need to determine their location using range detection. Since the Doppler shift distributes the signal power over both the I and Q channels, it is necessary to use the signal energy for detection. This means that incoherent detection schemes must be used.
The detection process is described in detail in the above example, so here we will simply perform the necessary actions to estimate the target range.
# Setting the probability of a false alarm (pfa)
pfa = 1e-6;
# In the system under consideration, the noise band is equal to half the sampling frequency.
noise_bw = receiver.SampleRate/2; # noise band
npower = noisepow1(noise_bw,
receiver.NoiseFigure,receiver.ReferenceTemperature); # noise power
threshold = npower * db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent")); # detection threshold
# Applying consistent filtering and updating the detection threshold
matchingcoeff = getMatchedFilter(waveform); # calculation of the coefficients of the agreed
# a filter based on a refined signal
matchedfilter = EngeePhased.MatchedFilter( # creating a consistent filter object
Coefficients=matchingcoeff, # filter coefficients
GainOutputPort=true); # gain output port
rxpulses, mfgain = matchedfilter(rxpulses); # response after consistent filtering
threshold = threshold * db2pow(mfgain); # calculation of the new detection threshold value
# Delay compensation in a matched filter
matchingdelay = size(matchingcoeff,1)-1; # the length of the matched filter
rxpulses = buffer(rxpulses[(matchingdelay+1):end],size(rxpulses,1));
# Application of temporary automatic gain control (VAR)
# to compensate for range-dependent losses.
prop_speed = radiator.PropagationSpeed; # signal propagation speed
range_gates = prop_speed*fast_time_grid/2; # range resolution
lambda = prop_speed/fc; # wavelength
tvg = EngeePhased.TimeVaryingGain( # creating a VAR system object
RangeLoss=2*fspl.(range_gates,lambda), # range-adjusted losses
ReferenceLoss=2*fspl.(prop_speed/(prf*2),lambda)); # range-adjusted losses
rxpulses = tvg(rxpulses); # WARU application
# Detection of peaks exceeding the detection threshold after integration
_,range_detect = findpeaks(pulsint(rxpulses,"noncoherent"),
"MinPeakHeight",sqrt(threshold)); # Finding peaks
range_estimates = round.(range_gates[range_detect]) # range estimates
These estimates suggest the presence of targets in the range of 2000 m and 3550 m.
The Doppler spectrum
After the target ranges have been successfully estimated, the Doppler information for each target can be estimated.
Doppler signal estimation is, in fact, the process of spectrum estimation. Therefore, the first step in Doppler processing is to calculate the Doppler spectrum of the received signal.
The received signal after the matched filter is a matrix, the columns of which correspond to the received pulses. Unlike range estimation, Doppler processing uses data for the pulse accumulation time (slow time), that is, along the rows of the resulting matrix. Since 10 pulses are used, 10 samples are available for Doppler processing. Since one sample is taken from each pulse, the sampling frequency of the Doppler samples is the pulse repetition rate (PRF).
As predicted by Fourier theory, the maximum single-digit Doppler shift that a pulse location system can detect is equal to half of its PRF. This also corresponds to the maximum single-digit speed that the radar system can detect. In addition, the number of pulses determines the resolution in the Doppler spectrum, which determines the resolution of the velocity estimate.
max_speed = dop2speed(prf/2,lambda)/2; # calculation
println("Maximum detection speed: $(round(max_speed;sigdigits=7)) m/s");
speed_res = 2*max_speed/num_pulse_int
println("Speed resolution: $(round(speed_res;sigdigits=7)) m/s");
As shown in the calculations above, in this example, the maximum detectable velocity is 225 m/s, either approaching (-225) or receding (+225). The resulting Doppler resolution is about 45 m/s, which means that for separation in the Doppler spectrum, the two speeds must differ from each other by at least 45 m/s. To improve the ability to distinguish between different target speeds, more pulses are needed. However, the number of available pulses is also limited by the radial velocity of the target. Since Doppler processing is limited to a set range, all pulses used in processing must be received before the target switches from one resolution element to another.
Since the number of Doppler samples is generally limited, it is common to reset the sequence to interpolate the resulting spectrum. This will not improve the resolution of the obtained spectrum, but it may improve the estimation of the location of peaks in the spectrum.
The Doppler spectrum can be obtained using a periodogram.
num_range_detected = length(range_estimates);
# Calculation of the periogram for a range of 2000 m
period_out_1 = DSP.periodogram(rxpulses[range_detect[1],:][:],nfft=256,fs=prf)
p1, f1 = period_out_1.power,period_out_1.freq # extracting power (p1) and frequency (f1) values
# Calculation of the periogram for a range of 3550 m
period_out_2 = DSP.periodogram(rxpulses[range_detect[2],:][:],nfft=256,fs=prf)
p2, f2 = period_out_2.power,period_out_2.freq; # extracting power (p2) and frequency (f2) values
Then you can calculate the speed corresponding to each frequency value in the spectrum. Please note that it is necessary to take into account the circular motion effect.
speed_vec = dop2speed.(fftshift(f1),lambda)/2; # recalculation of the frequency grid into the speed grid
Estimation of the Doppler shift
To estimate the Doppler shift associated with each target, we need to find the location of peaks in each Doppler spectrum. In this example, the targets are present on two different ranges, so the evaluation process must be repeated for each range.
First, we will construct the Doppler spectrum corresponding to a range of 2000 meters.
spec_sig1 = fft([rxpulses[range_detect[1],:][:];zeros(256-length(rxpulses[range_detect[1],:][:]))])
spec_power1 = fftshift(real(spec_sig1.*conj.(spec_sig1))./100)
plot(fftshift(fftfreq(256,prf)).*1e-3,pow2db.(spec_power1),lab="",gridalpha=0.15,margin=5*Plots.mm)
xlabel!("Frequency, kHz", fontfamily="Computer Modern")
ylabel!("SPM,dBW/Hz", fontfamily="Computer Modern")
title!("Spectral power density of the reflected signal at a range of 2000 m", fontfamily="Computer Modern")
Please note that we are only interested in detecting peaks, so the spectrum values themselves are irrelevant. It can be seen from the graph of the Doppler spectrum that 5 dB below the maximum peak is a good threshold. Therefore, we use the value of -5 as the threshold for the normalized Doppler spectrum.
spectrum_data = spec_power1./maximum(spec_power1)
_,dop_detect1 = findpeaks(pow2db.(spectrum_data),"MinPeakHeight",-5)
sp1 = speed_vec[dop_detect1]
println("Speed estimates: $sp1 m/s");
The results show that there are two targets at a range of 2000 m: one with a speed of 3.5 m/s and the other with a speed of 104 m/s. The value of -104 m/s can be easily associated with the first target, since the first target leaves at a radial velocity of 100 m/s, which, given the Doppler resolution of this example, is very close to the calculated value. The value of 3.5 m/s requires a more detailed explanation. Since the third target is moving in a tangential direction, there is no velocity component in the radial direction. Therefore, the radar cannot detect the Doppler shift of the third target. The true radial velocity of the third target is therefore 0 m/s, and the estimate of 3.5 m/s is very close to the true value.
Note that these two targets cannot be distinguished using range estimation alone, as their range values are the same.
Then the same operations are applied to the data corresponding to a range of 3550 meters.
spec_sig2 = fft([rxpulses[range_detect[2],:][:];zeros(256-length(rxpulses[range_detect[2],:][:]))])
spec_power2 = fftshift(real(spec_sig2.*conj.(spec_sig2))./100)
plot(fftshift(fftfreq(256,prf)).*1e-3,pow2db.(spec_power2),lab="",gridalpha=0.15,margin=5*Plots.mm)
xlabel!("Frequency, kHz", fontfamily="Computer Modern")
ylabel!("Power, dBW", fontfamily="Computer Modern")
title!("Spectral power density over a range of 3550 meters", fontfamily="Computer Modern")
spectrum_data2 = spec_power2./maximum(spec_power2)
_,dop_detect2 = findpeaks(pow2db.(spectrum_data2),"MinPeakHeight",-5)
sp2 = speed_vec[dop_detect2]
println("Speed rating: $sp2 m/s");
This result shows an estimated velocity of 0 m/s, which corresponds to the fact that the target is not moving at this distance.
Conclusion
This example shows a simple way to estimate the radial velocity of moving targets using a pulse radar system. The Doppler spectrum of the received signal was generated and the location of the peaks along the spectrum was estimated, which correspond to the radial velocity of the target. The limitations of Doppler processing caused by the characteristics of the emitted signal and the processing algorithm were also discussed.