Doppler estimation

The example shows a monostatic pulse radar determining the radial velocity of moving targets at specific ranges. The velocity is determined from the Doppler shift caused by the moving targets. We first determine the presence of a target at a given range, and then use Doppler processing to estimate the radial velocity of the target at that range.

Functions and libraries used

Pkg.add(["DSP"])

   Resolving package versions...
  No Changes to `~/.project/Project.toml`
  No Changes to `~/.project/Manifest.toml`

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;

Setting up the radar system

First, we define a radar system. Since the focus of this example is on Doppler processing, we use the radar system built in the example Test Signal Simulation for Radar Receiver. Readers are encouraged to explore the details of radar system design using this example.

waveform,antenna,radiator,collector,transmitter,receiver,sensormotion = BasicMonostaticRadarExampleData();

System modelling

Objectives

Doppler processing utilises the Doppler shift caused by a moving target. Now define three targets by 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 are at distances of 3550, 1600 and 1200 m , and are travelling at 60, 130 and 0 m/s, . The difference is that the first target is moving in the radial direction and the second target is moving in the tangential and radial directions. The third target is not moving.

Propagation environment

We also need to configure the propagation environment for each target. Since we are using a monostatic radar, we will set up a two-way propagation model:

fs = waveform.SampleRate; # частота дискретизации системы
# модель канала распространения
channel = EngeePhased.FreeSpace(
    SampleRate=fs, # частота дискретизации системы
    TwoWayPropagation=true, # модель распространения
    OperatingFrequency=fc); # несущая частота сигнала

Signal Synthesis

Having defined the radar system, environment and targets, we model the received echo signal reflected from the targets. The resulting data is a matrix of elements with fast time (time within each pulse) in each column and slow time (time between pulses) in each row.

To reproduce the same results, we need to set the input data to generate 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 estimate

Range detection

To estimate the Doppler shift of targets, you must first locate them using range detection. Since the Doppler shift distributes the signal power over both the I and Q channels, the signal energy must be utilised for detection. This means that non-coherent detection schemes must be used.

The detection process is described in detail in the above example, so here we will simply perform the necessary steps to estimate the range of the target.

# Задание вероятности ложной тревоги (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]) # оценки дальности

2-element Vector{Float64}:
 2000.0
 3550.0

These estimates suggest targets in the range of 2000 metres and 3550 metres.

Doppler spectrum

Once the target ranges have been successfully estimated, the Doppler information for each target can be estimated.

Estimating the Doppler signal is actually a spectrum estimation process. Therefore, the first step in Doppler processing is to estimate the Doppler spectrum of the received signal.

The received signal after the matched filter is a matrix whose columns correspond to the received pulses. Unlike range estimation, Doppler processing uses data over the pulse accumulation time (slow time), that is, over 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 Doppler sampling rate is the pulse repetition frequency (PRF).

As predicted by Fourier theory, the maximum single-digit Doppler shift that a pulse-based locating system can detect is half its PRF. This also corresponds to the maximum single-digit velocity that a radiolocation system can detect. In addition, the number of pulses determines the resolution in the Doppler spectrum, on which the resolution of the velocity estimate depends.

max_speed = dop2speed(prf/2,lambda)/2; # расчет 
println("Максимальная скорость обнаружения: $(round(max_speed;sigdigits=7)) м/с");

Максимальная скорость обнаружения: 224.6888 м/с

speed_res = 2*max_speed/num_pulse_int
println("Разрешение по скорости: $(round(speed_res;sigdigits=7)) м/с");

Разрешение по скорости: 44.93776 м/с

As shown in the calculations above, in this example the maximum detectable velocity is 225 m/s, either approaching (-225) or moving away (+225). The resulting Doppler resolution is about 45 m/s, which means that the two velocities must differ from each other by at least 45 m/s to be distinguished in the Doppler spectrum. To improve the ability to distinguish between different target velocities, 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 given range, all pulses used in processing must be received before the target moves from one resolution element to another.

Since the number of Doppler samples is generally limited, it is common to zero the sequence to interpolate the resulting spectrum. This will not improve the resolution of the resulting spectrum, but 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)

The velocity corresponding to each frequency value in the spectrum can then be calculated. Note that the effect of circular motion must be taken into account.

speed_vec = dop2speed.(fftshift(f1),lambda)/2; # пересчет сетки частот в сетку скорости

Estimating the Doppler shift

To estimate the Doppler shift associated with each target, we need to find the location of the peaks in each Doppler spectrum. In this example, the targets are present on two different bands, so the estimation process must be repeated for each band.

First, we plot the Doppler spectrum corresponding to a range of 2000 metres.

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")

Note that we are only interested in peak detection, so the spectrum values themselves are not important. From the plot of the Doppler spectrum we can see that 5 dB below the maximum peak is a good threshold. Therefore, we use the value -5 as the threshold for the normalised Doppler spectrum.

spectrum_data = spec_power1./maximum(spec_power1)
_,dop_detect1 = findpeaks(pow2db.(spectrum_data),"MinPeakHeight",-5)
sp1 = speed_vec[dop_detect1]
println("Оценки скорости:  $sp1 м/с");

Оценки скорости:  [-103.56749129975046, 3.510762416940684] м/с

The results show that there are two targets at a range of 2000 m, one with a velocity of 3.5 m/s and the other with a velocity of -104 m/s. The -104 m/s value can be easily associated with the first target, as the first target is travelling away 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 travelling in the tangential direction, there is no velocity component in the radial direction. Therefore, the radar cannot detect the Doppler displacement 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 the two targets cannot be distinguished using only the range estimate because their range values are the same.

The same operations are then applied to the data corresponding to a range of 3550 metres.

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 м/с");

Оценка скорости: [0.0] м/с

This result shows an estimated velocity of 0 m/s, which corresponds to the fact that the target is not moving at this range.

Conclusion

This example shows a simple way to estimate the radial velocity of moving targets using a pulsed radar system. We generated a Doppler spectrum of the received signal and estimated the location of peaks on the spectrum that correspond to the radial velocity of the target. We also discussed the limitations of Doppler processing caused by the characteristics of the radiated signal and the processing algorithm.