Modeling the microdopler effect
The example introduces the concept of the microdopler effect and demonstrates the construction of a radar system capable of extracting a microdopler from a reflected signal. This feature can be used in the future to solve the problem of goal classification.
Connecting libraries and functions
let
installed_packages = collect(x.name for (_, x) in Pkg.dependencies() if x.is_direct_dep)
list_packages = ["DSP","FFTW","LinearAlgebra"]
for pack in list_packages
pack in installed_packages || Pkg.add(pack)
end
end
using DSP,FFTW,LinearAlgebra
gr()
default(titlefontsize=12,top_margin=5Plots.px,guidefont=10,
fontfamily = "Computer Modern",colorbar_titlefontsize=8,size=(800,400),
margin = 2Plots.mm
)
function helicopmotion(t,tgtmotion,BladeAng,ArmLength,BladeRate,prf = 2e4, radarpos = [0;0;0])
Nblades = size(BladeAng,2)
tgtpos,tgtvel = tgtmotion(1/prf)
RotAng = BladeRate*t;
scatterpos = [0 ArmLength*cos.(RotAng .+ BladeAng);0 ArmLength*sin.(RotAng .+ BladeAng);zeros(1,Nblades+1)].+tgtpos
scattervel = [0 -BladeRate*ArmLength*sin.(RotAng .+ BladeAng);0 BladeRate*ArmLength*cos.(RotAng .+ BladeAng);zeros(1,Nblades+1)] .+tgtvel
_,scatterang = rangeangle(scatterpos,radarpos)
return scatterpos,scattervel,scatterang
end
function spectrogramm_ped(x,t_step;down_lim = -35)
out,f,t = EngeePhased.Functions.spectrogram(sum(x;dims=1)[:];
window = DSP.kaiser(128,10),
noverlap=120,
nfft=256,fs=1/t_step,
centered=true,freqloc = "yaxis",out = :data
)
out_sig = map(x -> x < down_lim ? down_lim : x,DSP.pow2db.(abs.(out)))
fig = Plots.heatmap(
t[:],f,out_sig,color= :jet,
gridalpha=0.5,
colorbar_title = "Мощность, дБВт",
xticks=0:0.5:3,
yticks=-500:100:500,
)
Plots.xlabel!("Время, с")
Plots.ylabel!("Частота Доплера, Гц")
return fig
end
function calc_spectrogram_quad(Data,ts,T;limit = -150)
x = Data./maximum(abs.(Data))
dT = T/length(x)
F = 1/dT
spec_quad = DSP.spectrogram(x,128,120;window=DSP.kaiser(128,3.95),nfft=512,fs=F)
power_quad = map(x -> x < limit ? limit : x,pow2db.(spec_quad.power))
power_quad_quad_new = zeros(size(power_quad))
power_quad_quad_new[1:round(Int64,size(power_quad,1)/2),:] .= power_quad[round(Int64,size(power_quad,1)/2)+1:end,:]
power_quad_quad_new[round(Int64,size(power_quad,1)/2)+1:end,:] .= power_quad[1:round(Int64,size(power_quad,1)/2),:]
fig = Plots.heatmap(spec_quad.time*1e3, fftshift(spec_quad.freq)*1e-3, power_quad_quad_new,color= :jet,
gridalpha=0.5,ylims=(-10,10),xlims=(0,500),colorbar_title = "Мощность, дБВт",
xticks=0:50:500,yticks=-10:2:10)
Plots.xlabel!("Время, мс")
Plots.ylabel!("Частота Доплера, кГц")
Plots.title!("Спектрограмма отраженного сигнала от вертолета")
return fig
end
function helperAnnotateSpectrogram(fig)
Plots.plot!(fig,[20;50],[4;7],linewidth=2,c=:white,lab="",arrow=Plots.Arrow(:close, :tail, 0.5, 1))
Plots.annotate!(50, 7.5, ("Лопасть 1", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[70;120],[4;5.5],linewidth=2,c=:white,lab="",arrow=Plots.Arrow(:close, :tail, 0.5, 1))
Plots.annotate!(120, 6, ("Лопасть 2", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[20;50],[-4;-7],linewidth=2,c=:white,lab="",arrow=Plots.Arrow(:close, :tail, 0.5, 1))
Plots.annotate!(50, -7.5, ("Лопасть 3", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[70;120],[-4;-5.5],linewidth=2,c=:white,lab="",arrow=Plots.Arrow(:close, :tail, 0.5, 1))
Plots.annotate!(120, -6, ("Лопасть 4", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[0;250],[0;-0.1],linewidth=3,c=:white,lab="",arrow=Plots.Arrow(:close, :both, 0.5, 1))
Plots.annotate!(120, -1, ("Период", 10, :white, :center),annotationfontfamily="Roboto")
Plots.plot!(fig,[124;125],[0.2;4],linewidth=3,c=:white,lab="",arrow=Plots.Arrow(:close, :both, 0.5, 1))
Plots.annotate!(137, 2, ("Скорость вращения", 10, :white, :left),annotationfontfamily="Roboto")
end
function plotting_svd(x)
Plots.plot(10*log10.(x),lab="",ylabel="Сингулярные значения",xlabel="ранг")
Plots.plot!([55], seriestype="vline", label="",ls=:dashdot,color="red",lw=2)
Plots.plot!([101], seriestype="vline", label="",ls=:dashdot,color="red",lw=2)
Plots.plot!([109], seriestype="vline", label="",ls=:dashdot,color="red",lw=2)
Plots.annotate!(25, -10, ("A", 10, :black, :left),annotationfontfamily="Roboto")
Plots.annotate!(75, -10, ("B", 10, :black, :left),annotationfontfamily="Roboto")
Plots.annotate!(102, -10, ("C", 10, :black, :left),annotationfontfamily="Roboto")
Plots.annotate!(175, -10, ("D", 10, :black, :left),annotationfontfamily="Roboto")
end
1. The concept of the microdopler effect
The concept of micro-movements is inextricably linked with the microdopler effect*:
Micro–movement is a slight movement or deviation of a component part of an object, which is additional to the main volumetric movement of the whole object [1]. The source of micro-motion can be a vibrating surface, rotating rotor blades of a helicopter, a walking person making periodic movements with his arms and legs, flapping wings of a bird, or other reasons.
The microdopler effect is a phenomenon associated with the enrichment of the spectral composition of the reflected signal, which occurs during radar research as a result of micro-movements of objects.
Let's consider 2 scenarios for the movement of an object: a simple translational movement (top picture) and a complex translational-rotational movement (bottom picture).
In case 1, the object is moving towards us uniformly, so the Doppler frequency has a constant value, independent of time.
In the second case, the situation changes: rotational motion (for example, the rotation of helicopter blades) introduces periodically changing components into the spectrum of the reflected signal, which characterize the microdopler effect.
A spectrogram of a reflected signal containing Microdopler components in the spectrum is called Microdopler signatures (or abbreviated as MDS). An example of a bird's MDS is given below:

2. Simulation of pedestrian MDS
In the next section, we will consider the creation of an automotive radar based on a continuous signal from an LFM to identify a crossing on the road.
2.1 Descriptions of the model scenario
The model considers the following scenario:
- a car with a radar is moving along the inside lane at a speed of 30 km/h from the origin towards a parked car and a pedestrian;
- The parked car is located at a distance of 39.4 m
- The pedestrian is moving at a speed of 1 m/s
A visualization of the scenario is given below:
.png)
2.2 Model Parameters
High range and speed resolution is required to identify the transition. In automotive radars, a continuous linear frequency modulated (FM) signal is used to meet the characteristics.
bw = 250e6 # полоса сигнала, Гц
fs = bw # частота дискретизации, Гц
fc = 24e9 # несущая частота, Гц
tm = 1e-6 # время развертки, с
c = 3e8 # скорость света, м/с
egocar_pos = [0;0;0] # м, вектор положения автомобиля с радаром
egocar_vel = [30*1600/3600;0;0] # м/с, вектор скорости автомобиля с радаром
parkedcar_pos = [39;-4;0] # м, вектор положения впереди стоящего автомобиля
parkedcar_vel = [0;0;0] # м/с, вектор скорости впереди стоящего автомобиля
ped_pos = [40;-3;0] # м, вектор положения пешехода
ped_vel = [0;1;0] # м, вектор скорости пешехода
ped_heading = 90 # град, начальное направление движения пешехода
ped_height = 1.8; # рост пешехода
2.3 Formation of system objects
To formalize the reflected signal of a pedestrian object, use the system objects EngeeRadars.BackscatterPedestrian
wav_fmcw = EngeePhased.FMCWWaveform(
SampleRate=fs, # частота дискретизации
SweepTime=tm, # время развертки
SweepBandwidth=bw # полоса сигнала
)
# СО автомобиля с радаром
egocar = EngeePhased.Platform(
InitialPosition=egocar_pos, # начальное положение
Velocity=egocar_vel, # вектор скорости
OrientationAxesOutputPort=true # включение выходного порта ориентации в пространстве
)
# Формирование объекта впереди стоящего автомобиля
parkedcar = EngeePhased.Platform(
InitialPosition=parkedcar_pos, # начальное положение
Velocity=parkedcar_vel, # скорость
OrientationAxesOutputPort=true # включение выходного порта ориентации в пространстве
)
parkedcar_tgt = EngeePhased.RadarTarget(PropagationSpeed=c,OperatingFrequency=fc,MeanRCS=10)
# Формирование объекта - пешеход
ped = EngeeRadars.BackscatterPedestrian(
InitPosition=ped_pos, # начальное положение
InitHeading=ped_heading, # начальное направление движения
PropagationSpeed=c, # скорость распространения сигнала
OperatingFrequency=fc, # несущая частота
Height=ped_height, # рост
WalkingSpeed=1 # скорость движения
)
# каналы распространения: до пешехода и ближайшего автомобиля
chan_ped = EngeePhased.FreeSpace(
PropagationSpeed=c, # скорость распространения
OperatingFrequency=fc, # несущая частота
TwoWayPropagation=true, # учет двунаправленного распространения
SampleRate=fs # частота дискретизации
)
chan_pcar = EngeePhased.FreeSpace(
PropagationSpeed=c, # скорость распространения
OperatingFrequency=fc, # несущая частота
TwoWayPropagation=true, # учет двунаправленного распространения
SampleRate=fs # частота дискретизации
)
# приемник и передатчик
tx = EngeePhased.Transmitter(PeakPower=1,Gain=25)
rx = EngeePhased.ReceiverPreamp(Gain=25,NoiseFigure=10);
2.4 Simulation of the model
For the correct display of the pedestrian's MDS, it is necessary to accumulate the reflected signal during at least one cycle of human movement. The pedestrian's speed is 1 m/s, so it takes 2 seconds to complete the cycle. Let's simulate a pedestrian movement scenario for 2.5 seconds in 1 ms increments.:
Tsim = 2.5 # время симуляции
Tsamp = 0.001 # шаг расчета симуляции, с
npulse = round(Int64,Tsim/Tsamp) # количество импульсов
# Выделение памяти под:
xr = complex(zeros(round(Int64,fs*tm),npulse)) # суммарный отраженный сигнал
xr_ped = complex(zeros(round(Int64,fs*tm),npulse)) # отраженный сигнал от пешехода
@inbounds for m = 1:npulse
# расчет текущего положения, скорости, ориентации и угла визирования
# до пешехода и припаркованного автомобиля
pos_ego,vel_ego,ax_ego = egocar(Tsamp)
pos_pcar,vel_pcar,ax_pcar = parkedcar(Tsamp)
pos_ped,vel_ped,ax_ped = move(ped,Tsamp,ped_heading)
_,angrt_ped = rangeangle(pos_ego,pos_ped,ax_ped)
_,angrt_pcar = rangeangle(pos_ego,pos_pcar,ax_pcar)
# расчет отраженного сигнала
global x = tx(wav_fmcw()) # зондирующий сигнал
# прохождение сигналов в канале прямого и обратного распространения
xt_ped = chan_ped(repeat(x,1,size(pos_ped,2)),pos_ego,pos_ped,vel_ego,vel_ped)
xt_pcar = chan_pcar(x,pos_ego,pos_pcar,vel_ego,vel_pcar)
xt_ped = reflect(ped,xt_ped,angrt_ped) # отраженный сигнал от перехода
xt_pcar = parkedcar_tgt(xt_pcar) # отраженный сигнал от автомобиля
xr_ped[:,m] = rx(xt_ped) # принятый сигнал от пешехода
xr[:,m] = rx(xt_ped+xt_pcar) # принятый сигнал от пешехода
end
# перенос сигнала на нулевую частоту
xd_ped = conj(dechirp(xr_ped,x))
xd = conj(dechirp(xr,x));
2.5 Processing of the received signal
In the simulated signal xd_ped it contains only the pedestrian's response, and xd It contains the response of both a pedestrian and a parked car. If we create a spectrogram using only the reflected signal from a pedestrian, we get the following MDS:
spectrogramm_ped(xd_ped,Tsamp;down_lim = -35)
title!("Спектрограмма отраженного сигнала от пешехода")
The spectrogram of the total reflected signal from a pedestrian and a car looks like this:
spectrogramm_ped(xd,Tsamp;down_lim = -35)
title!("Спектрограмма отраженного сигнала от пешехода и автомобиля")
It can be noted that the vehicle's MDS dominates the pedestrian's MDS. Therefore, a pedestrian cannot be detected without additional processing.
2.6 SVD decomposition processing
To detect a weaker signal based on a strong one, you can use singular value decomposition to identify the area of interest.
out_svd = svd(xd,full=true) # расчет компонента разложения
uxd,sxd,vxd = out_svd.U,out_svd.S,out_svd.V; # считывание матриц S,U,V
Using the function plotting_svd let's display the main diagonal singular matrix S:
plotting_svd(sxd)
On the graph of the dependence of the signal values on the rank of the matrix, 4 areas can be conditionally distinguished:
- A - includes the main signal, which makes a significant contribution to the spectrogram;
- B - a combination of pedestrian and car signals, with the dominance of the second
- C - a combination of pedestrian and car signals, with the former dominating;
- D is the noise contribution of the receiving device.
Therefore, in order to effectively isolate the pedestrian's MDS, it is necessary to extract from the entire spectrogram the area C located within the range of the rank 103-110:
rk = 103:110 # диапазон рангов матрицы сигнулярных значений
# расчет обработанной спектрограммы
sxd_new = [sxd;zeros(size(vxd,1)-length(sxd))]
xdr = uxd[:,rk]*Diagonal(sxd_new)[rk,:]*vxd';
We will build the pedestrian and car MDS again after additional processing (for correct display, we will reduce the lower limit to -40 dB since the main power has been removed):
spectrogramm_ped(xdr,Tsamp;down_lim = -40)
title!("МДС пешехода и автомобиля после svd преобразования")
After filtering the signal from the car, a clear picture of the pedestrian's MDS is formed, while the level of the car's MDS has significantly decreased.
Thus, with the help of additional processing (svd decomposition), it was possible to isolate the spectral components of a pedestrian against the background of an interfering signal from a nearby car.
3. Helicopter MDS simulation
3.1 Description of the model scenario
Consider the case:
- The helicopter is flying at a distance of [-25,500,500] meters evenly in the x-z plane at a speed of 100 m/s along the x axis.
- The radar is stationary and located at the origin.
In this case, at the beginning of the movement, the helicopter approaches the radar, and after moving to the positive area of the x-axis, it will begin to move away. A visualization of this scenario is given below.:

3.2 Model parameters
Let's set the parameters of the radar and target (helicopter) movement scenario:
radarpos = [0;0;0] # вектор положения РЛС, м
radarvel = [0;0;0] # вектор скорости РЛС, м/с
tgtinitpos = [-25;500;500] # начальное положение цели, м
tgtvel = [100;0;0]; # скорость цели, м/с
Next, we will form the geometry of the helicopter, specifying its main parameters (number of blades, length of blades, speed of rotation of the propeller):
f_rot = 4 # частота вращения винта, c
Nblades = 4 # количество лопастей
bladeang = Matrix((0:Nblades-1)')*2*pi/Nblades # вектор углов лопастей
bladelen = 6.5 # длина лопасти, м
bladerate = deg2rad(360*f_rot); # скорость вращения винта рад/с
The next stage is the formation of a radar station. Let's set the fundamental parameters of the system:
c = 3e8 # скорость распространения сигнала, м/с
fc = 5e9 # несущая частота, Гц
fs = 1e6 # частота дискретизации, Гц
prf = 2e4 # частота следования импульсов, Гц
lambda = c/fc; # длина волны, м
3.3 Creating system objects
Now, using the parameters of the target scenario and the radar system, we will create system model objects using the library EngeePhased:
# Формирование цели
helicop = EngeePhased.RadarTarget(
MeanRCS=[10 .1 .1 .1 .1], # ЭПР элементов вертолета (фюзеляжа и 4 лопастей)
PropagationSpeed=c, # скорость распространения
OperatingFrequency=fc # несущая частота сигнала
)
tgtmotion = EngeePhased.Platform(InitialPosition=tgtinitpos,Velocity=tgtvel)
# Генератор зондирующего сигнала
wav = EngeePhased.RectangularWaveform(
SampleRate=fs, # частота дискретизации
PulseWidth=2/fs, # длительность импульса
PRF=prf # частота следования импульсов
)
# Формирование антенной решетки (прямоугольная 4х4)
ura = EngeePhased.URA(
Size=[4 4], # размер антенной решетки
ElementSpacing=[lambda/2 lambda/2] # расстояние между элементами
)
tx = EngeePhased.Transmitter() # передатчик
rx = EngeePhased.ReceiverPreamp() # предусилитель
# Cреда распространения
env = EngeePhased.FreeSpace(
PropagationSpeed=c, # скорость распространения, м/с
OperatingFrequency=fc, # несущая частота, Гц
TwoWayPropagation=true, # учет двунаправленного распространения
SampleRate=fs # частота дискретизации
)
# Передающая АР
txant = EngeePhased.Radiator(
Sensor=ura, # геометрия АР
PropagationSpeed=c, # скорость распространения, м/с
OperatingFrequency=fc # несущая частота, Гц
)
# Принимающая АР
rxant = EngeePhased.Collector(
Sensor=ura, # геометрия АР
PropagationSpeed=c, # скорость распространения, м/с
OperatingFrequency=fc # несущая частота, Гц
);
3.4 Running a simulation of the model
We implement the calculation of the radar operation scenario when the helicopter is moving by calling system objects in a loop.:
NSampPerPulse = round(Int64,fs/prf) # количество отсчетов в 1 импульсе
Niter = Int64(1e4) # количество итераций (число излученных импульсов)
y = complex(zeros(NSampPerPulse,Niter)) # выделение памяти для отраженного сигнала
@inbounds for m in 1:Niter
# обновление положение и скорости вертолета
t = (m-1)/(prf) # шаг симуляции по времени
global scatterpos,scattervel,scatterang = helicopmotion(t,tgtmotion,bladeang,bladelen,bladerate,prf,radarpos)
# расчет отраженного сигнала
x = txant(tx(wav()),scatterang) # зондирующий сигнал
xt = env(x,radarpos,scatterpos,radarvel,scattervel) # сигнал, прошедший среду
xt = helicop(xt) # отраженный сигнал от вертолета
global xr = rx(rxant(xt,scatterang)) # принятый сигнал
y[:,m] = sum(xr;dims=2) # формирование луча
end
3.5 Signal processing
Using a system object RangeDopplerResponse we implement range-Doppler processing:consistent range filtering and spectral processing based on speed FFT.
rdresp = EngeePhased.RangeDopplerResponse(
PropagationSpeed=c, # скорость распространения сигнала
SampleRate=fs, # частота дискретизации
DopplerFFTLengthSource="Property", # Cпособ задания БПФ (в параметрах СО)
DopplerFFTLength=128, # длина БПФ
DopplerOutput="Speed", # Доплероский выход (скорость)
OperatingFrequency=fc # несущая частота
)
mfcoeff = getMatchedFilter(wav)
fig1 = plotResponse(rdresp,y[:,1:128],mfcoeff)
Plots.ylims!(0,3000)
Plots.plot!(fig1,title="Дальностно-доплероский портрет",xlabel="Скорость, м/с",
ylabel="Дальность, м",titlefontsize=12,top_margin=5Plots.px,guidefont=10,
colorbar_title="Мощность, дБВт",fontfamily = "Computer Modern",colorbar_titlefontsize=8)
It can be noted that the range-Doppler portrait captures the presence of 3 components: the central component is the fuselage of the helicopter, and the lateral components are the rotation of the blades.
Let's estimate the radial velocity of the helicopter:
tgtpos = scatterpos[:,1] # положение вертолета на момент окончания симуляции
tgtvel = scattervel[:,1] # скорость вертолета на момент окончания симуляции
tgtvel_truth = radialspeed(tgtpos,tgtvel,radarpos,radarvel)
println("радиальная скорость фюзеляжа вертолета $(round(tgtvel_truth;sigdigits=3)) м/с")
We will also look at the maximum and minimum values of the radial speeds of the screws.:
maxbladetipvel = [bladelen*bladerate;0;0] # скорость лопастей винта
vtp = radialspeed(tgtpos,-maxbladetipvel+tgtvel,radarpos,radarvel)
vtn = radialspeed(tgtpos,maxbladetipvel+tgtvel,radarpos,radarvel)
println("максимальное значение радиальной скорости $(round(vtp;sigdigits=3)) м/с")
println("минимальное значение радиальной скорости $(round(vtn;sigdigits=3)) м/с")
Let's perform consistent filtering and construct a spectrogram of the reflected signal.:
# Cогласованная фильтрация
mf = EngeePhased.MatchedFilter(Coefficients=mfcoeff) # согласованный фильтр
ymf = mf(y) # обработка
ridx = findmax(sum(abs.(ymf);dims=2)[:])[2]; # обнаружение пика по дальности
# Построение МДС
T = Niter/prf # длительность спектрограммы
fig1 = calc_spectrogram_quad(ymf[ridx,:],1/prf,T;limit=-100)
The central component characterizes the reflected signal from the helicopter fuselage, the periodic components are the rotation of the blades. We visualize the components of the MDS using the function helperAnnotateSpectrogram:
helperAnnotateSpectrogram(fig1)
The oscillation period, which is approximately 250 ms, can be estimated from the graph of the spectrogram. In this case, the estimated number of blades is:
Tp = 250e-3 # период МДС
bladerate_est = 1/Tp # оценка количества лопастей
println("Оценка числа лопастей $(round(bladerate_est;sigdigits=2))")
The maximum deviation of the Doppler frequency is approximately 4 kHz, then the radial velocity of the helicopter blades is estimated as:
f_doppler = 4e3 # амплитуда доплероской частоты, Гц
Vt_detect = dop2speed(f_doppler,lambda)/2
println("Радиальная скорость лопастей $(Vt_detect) м/с")
To estimate the true speed, it is necessary to determine and solve the pelingation problem based on the angle of the location. Let's use the "Music" algorithm implemented in the system object. MUSICEstimator2D:
doa = EngeePhased.MUSICEstimator2D(
SensorArray=ura, #
OperatingFrequency=fc, # несущая частота
PropagationSpeed=c, # скорость распространения сигнала
DOAOutputPort=true, # включение выхода углового направления
ElevationScanAngles=-90:90 # диапазон сканирования по углу места
)
_,ang_est = doa(xr) # расчет углов визирования
Vt_est = Vt_detect/cosd(ang_est[2]) # оценка скорости движения лопастей
println("Оценка скорости лопастей $(round(Vt_est;sigdigits=3)) м/с")
Knowing the estimates of the speed and number of blades, we will determine the approximate length of the blade.:
bladelen_est = Vt_est/(bladerate_est*2*pi)
println("Оценка длины лопасти $(round(bladelen_est;sigdigits=3)) м")
The true length value is 6.5 m, hence the error was 4%.
Conclusion
In the example, the concept of the "microdopler effect" and Microdopler signatures is considered. Using the EngeePhased system objects, MDS reflected signals were modeled for 2 radar scenarios: a ground-based radar in the presence of an aerial target - a helicopter and an automobile radar in the presence of a pedestrian and a nearby car.
