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 # cкорость распространения сигнала
# Создание системных объектов:
# генератор прямоугольных импульсов
waveform = EngeePhased.RectangularWaveform(PulseWidth=100/c,
PRF=c/1e4,OutputFormat="Pulses",NumPulses=1,SampleRate= 2c/100);
# антенный элемент (тип - изотропный)
antenna = EngeePhased.IsotropicAntennaElement(FrequencyRange=[5e9 1.5e10]);
# излучатель узкополосного сигнала
radiator = EngeePhased.Radiator(Sensor=antenna,
PropagationSpeed=c,OperatingFrequency=1e10);
# Усилитель в передающем тракте
transmitter = EngeePhased.Transmitter(PeakPower=5226.4791642003665401716716587543,
Gain=20,LossFactor=0, InUseOutputPort=true,CoherentOnTransmit=true);
# принимающая антенна
collector = EngeePhased.Collector(Sensor=antenna,
PropagationSpeed=c,Wavefront="Plane",
OperatingFrequency=1e10);
# Предусилитель в приемной тракте
receiver = EngeePhased.ReceiverPreamp(Gain=20,NoiseFigure=0,
ReferenceTemperature=290,SampleRate=2*c/100, PhaseNoiseInputPort=false,
EnableInputPort=true,SeedSource="Property",Seed=1000);
# модель движения радара
sensormotion = EngeePhased.Platform(InitialPosition=[0;0;0],Velocity=[0;0;0]);
return waveform, antenna, radiator, collector, transmitter, receiver, sensormotion
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;
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]]; # компоненты начального положения целей
tgtvel = [[60; 80; 0] [0;0;0] [0; 100; 0]]; # компоненты скорости целей
tgtmotion = EngeePhased.Platform(InitialPosition=tgtpos,Velocity=tgtvel); # модель движения целей
tgtrcs = [1.3 1.7 2.1]; # ЭПР целей
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; # частота дискретизации системы
# модель канала распространения
channel = EngeePhased.FreeSpace(
SampleRate=fs, # частота дискретизации системы
TwoWayPropagation=true, # модель распространения
OperatingFrequency=fc); # несущая частота сигнала
Signal synthesis
Having identified 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) # сброс внутреннего состояния объекта (измненеие )
receiver.SeedSource = "Property"; # изменение тип модели генерации шума
receiver.Seed = 2024; # задание начального состояния генератора тепловых шумов
prf = waveform.PRF; # частота следования импульсов
num_pulse_int = 10; # количество импульсов накопления
fast_time_grid = (0:1/fs:1/prf-1/fs) # сетка быстрого времени
slow_time_grid = (0:num_pulse_int-1)./prf; # сетка медленного времени
# Выделение памяти для результирующих данных
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);
# Расчет углов по дальности от радара до целей
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
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.
# Задание вероятности ложной тревоги (pfa)
pfa = 1e-6;
# В рассматриваемой системе, шумовая полоса равна половине частоты дискретизации
noise_bw = receiver.SampleRate/2; # шумовая полоса
npower = noisepow1(noise_bw,
receiver.NoiseFigure,receiver.ReferenceTemperature); # мощность шума
threshold = npower * db2pow(npwgnthresh(pfa,num_pulse_int,"noncoherent")); # порог обнаружения
# Применение согласованной фильтрации и обновление порога обнаружения
matchingcoeff = getMatchedFilter(waveform); # расчет коэффициентов согласованного
# фильтра на основе изулченного сигнала
matchedfilter = EngeePhased.MatchedFilter( # создание объекта согласованного фильтра
Coefficients=matchingcoeff, # коэффициенты фильтра
GainOutputPort=true); # выходной порт усиления
rxpulses, mfgain = matchedfilter(rxpulses); # отклик после согласованной фильтрации
threshold = threshold * db2pow(mfgain); # расчет нового значения порога обнаружения
# Компенсация задержки в согласованном фильтре
matchingdelay = size(matchingcoeff,1)-1; # длина согласованного фильтра
rxpulses = buffer(rxpulses[(matchingdelay+1):end],size(rxpulses,1));
# Применение временного автоматического регулирования усиления (ВАРУ)
# для компенсации потерь, зависящие от дальности.
prop_speed = radiator.PropagationSpeed; # скорость распространения сигнала
range_gates = prop_speed*fast_time_grid/2; # разрешение по дальности
lambda = prop_speed/fc; # длина волны
tvg = EngeePhased.TimeVaryingGain( # создание системного объекта ВАРУ
RangeLoss=2*fspl.(range_gates,lambda), # потери с учетом дальности
ReferenceLoss=2*fspl.(prop_speed/(prf*2),lambda)); # потери с учетом дальности
rxpulses = tvg(rxpulses); # применение ВАРУ
# Обнаружения пиков, привышающих порог обнаружения, после интегрирования
_,range_detect = findpeaks(pulsint(rxpulses,"noncoherent"),
"MinPeakHeight",sqrt(threshold)); # нахождение пиков
range_estimates = round.(range_gates[range_detect]) # оценки дальности
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; # расчет
println("Максимальная скорость обнаружения: $(round(max_speed;sigdigits=7)) м/с");
speed_res = 2*max_speed/num_pulse_int
println("Разрешение по скорости: $(round(speed_res;sigdigits=7)) м/с");
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);
# Расчет периограммы для дальности 2000 м
period_out_1 = DSP.periodogram(rxpulses[range_detect[1],:][:],nfft=256,fs=prf)
p1, f1 = period_out_1.power,period_out_1.freq # излечение значений мощности (p1) и частот (f1)
# Расчет периограммы для дальности 3550 м
period_out_2 = DSP.periodogram(rxpulses[range_detect[2],:][:],nfft=256,fs=prf)
p2, f2 = period_out_2.power,period_out_2.freq; # излечение значений мощности (p2) и частот (f2)
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; # пересчет сетки частот в сетку скорости
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!("Частота,кГц", fontfamily="Computer Modern")
ylabel!("СПМ,дБВт/Гц", fontfamily="Computer Modern")
title!("Спектральная плотность мощности отраженного сигнала по дальности 2000 м", 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("Оценки скорости: $sp1 м/с");
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!("Частота,кГц", fontfamily="Computer Modern")
ylabel!("Мощность,дБВт", fontfamily="Computer Modern")
title!("Спектральная плотность мощности по дальности 3550 метров", fontfamily="Computer Modern")
spectrum_data2 = spec_power2./maximum(spec_power2)
_,dop_detect2 = findpeaks(pow2db.(spectrum_data2),"MinPeakHeight",-5)
sp2 = speed_vec[dop_detect2]
println("Оценка скорости: $sp2 м/с");
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.