This example introduces the basic concept of micro-Doppler effect present in the radar return of a target due to the target's rotation. Such micro-Doppler signature can be used in the processing stage to help identify the target.
It is well known that a moving target introduces a frequency shift in the radar return due to Doppler effect. However, most targets are not rigid bodies. Therefore, there are often other vibrations and rotations in different parts of the target in addition to the platform movement. For example, when a helicopter is flying, its blades rotate, or when a person is walking, his or her arms swing naturally. All these micro-scale movements produces additional Doppler shifts, referred to as micro-Doppler effects, and prove to be very useful to identify target features. This example shows how micro-Doppler signatures can be simulated and processed to determine a helicopter's blade speed.
Consider a helicopter with four rotor blades, Assume the radar is located at the origin and the helicopter is at [500;0;500] meters away with a velocity of [60;0;0] m/s.
radarpos = [0;0;0]; radarvel = [0;0;0]; tgtinitpos = [500;0;500]; tgtvel = [60;0;0]; tgtmotion = phased.Platform('InitialPosition',tgtinitpos,'Velocity',tgtvel);
In the following simulation, the helicopter is modeled by five scatterers: the rotation center and the tips of four blades. The rotation center moves with the helicopter body. The four blades are 90 degrees apart and are rotating at a constant speed of 4 revolutions per second. The arm length of each blade is 6.5 meters.
Nblades = 4; bladeang = (0:Nblades-1)*2*pi/Nblades; bladelen = 6.5; bladerate = deg2rad(4*360); % rps -> rad/sec
All four tips are assumed to have identical reflectivities while the reflectivity for the rotation center is stronger.
c = 3e8; fc = 5e9; helicop = phased.RadarTarget('MeanRCS',[10 .1 .1 .1 .1],'PropagationSpeed',c,... 'OperatingFrequency',fc);
Assume the radar operates at 5 GHz with a simple pulse. The pulse repetition frequency is 20 kHz. For simplicity, assume the signal propagates in free space.
fs = 1e6; prf = 2e4; lambda = c/fc; wav = phased.RectangularWaveform('SampleRate',fs,'PulseWidth',2e-6,'PRF',prf); ura = phased.URA('Size',4,'ElementSpacing',lambda/2); tx = phased.Transmitter; rx = phased.ReceiverPreamp; env = phased.FreeSpace('PropagationSpeed',c,'OperatingFrequency',fc,... 'TwoWayPropagation',true,'SampleRate',fs); txant = phased.Radiator('Sensor',ura,'PropagationSpeed',c,'OperatingFrequency',fc); rxant = phased.Collector('Sensor',ura,'PropagationSpeed',c,'OperatingFrequency',fc);
At each pulse, the helicopter moves along its trajectory. Meanwhile, the blades keeps rotating and the tips of the blades introduce additional displacement and angular speed.
NSampPerPulse = round(fs/prf); Niter = 1e4; y = complex(zeros(NSampPerPulse,Niter)); rng(2018); for m = 1:Niter % update helicopter motion t = (m-1)/prf; [scatterpos,scattervel,scatterang] = helicopmotion(t,tgtmotion,bladeang,bladelen,bladerate); % simulate echo x = txant(tx(wav()),scatterang); % transmit xt = env(x,radarpos,scatterpos,radarvel,scattervel); % propagates to/from scatterers xt = helicop(xt); % reflect xr = rx(rxant(xt,scatterang)); % receive y(:,m) = sum(xr,2); % beamform end
The following figure shows the range-Doppler response using first 128 pulses of the received signal. You can see that there are three returns at the target range of approximately 700 meters.
rdresp = phased.RangeDopplerResponse('PropagationSpeed',c,'SampleRate',fs,... 'DopplerFFTLengthSource','Property','DopplerFFTLength',128,'DopplerOutput','Speed',... 'OperatingFrequency',fc); mfcoeff = getMatchedFilter(wav); plotResponse(rdresp,y(:,1:128),mfcoeff); ylim([0 3000])
At first sight, it is easy to mistake the returns as if they were from three different targets, but in fact all returns are from the same target. The center blob is from the rotation center. It is much stronger compared to the other two blobs because the reflection is stronger from the helicopter body. The plot shows a speed of -40 m/s for the rotation center. This matches the truth of the target radial speed given by
tgtpos = scatterpos(:,1); tgtvel = scattervel(:,1); tgtvel_truth = radialspeed(tgtpos,tgtvel,radarpos,radarvel)
tgtvel_truth = -43.6435
The other two blobs are from the tips of the blades when they approach or depart from the target at maximum speed. From the plot, the speeds corresponding to these two detections are about 75 m/s and -160 m/s, respectively. These can be verified as
maxbladetipvel = [bladelen*bladerate;0;0]; vtp = radialspeed(tgtpos,-maxbladetipvel+tgtvel,radarpos,radarvel) vtn = radialspeed(tgtpos,maxbladetipvel+tgtvel,radarpos,radarvel)
vtp = 75.1853 vtn = -162.4723
It is beyond the scope of this example, but it is possible to associate all three detections to the same target via further processing.
The time-frequency representation of micro-Doppler can reveal more information. The following code constructs a time-frequency representation in the detected target range bin.
mf = phased.MatchedFilter('Coefficients',mfcoeff); ymf = mf(y); [~,ridx] = max(sum(abs(ymf),2)); % detection via peak finding along range pspectrum(ymf(ridx,:),prf,'spectrogram')
The figure shows the micro-Doppler modulation caused by blade tips around a constant Doppler shift. The image suggests that each blade tip introduces a sinusoid-like Doppler modulation. As noted in the figure, within each period of the sinusoid, there are three extra sinusoids appearing at equal distance. This suggests that the helicopter is equipped with four equally spaced blades.
hanno = helperAnnotateMicroDopplerSpectrogram(gcf);
In addition to the number of blades, it can also be read from the image that the period of each sinusoid, Tr, is about 250 ms. This means that it takes 250 ms for a blade to return to its original position, which in turn suggests that the angular speed of the helicopter is about 4 revolutions per second, matching the simulation parameter.
Tp = 250e-3; bladerate_est = 1/Tp
bladerate_est = 4
In addition, it is shown that the tip velocity, Vt, can be derived from the maximum Doppler. The maximum Doppler is about 4 kHz away from the constant Doppler introduced by the bulk movement. Thus, the detected maximum tip velocity is given by
Vt_detect = dop2speed(4e3,c/fc)/2
Vt_detect = 120
It is worth noting that this is the maximum tip velocity along the radial direction. To obtain the correct maximum tip velocity, the relative orientation must be taken into consideration. Looking at the scene, because blades are spinning as a circle, the detection is not affected by the azimuth angle. Therefore, the result only needs to be corrected for the elevation angle.
doa = phased.MUSICEstimator2D('SensorArray',ura,'OperatingFrequency',fc,... 'PropagationSpeed',c,'DOAOutputPort',true,'ElevationScanAngles',-90:90); [~,ang_est] = doa(xr); Vt_est = Vt_detect/cosd(ang_est(2))
Vt_est = 164.0793
Based on the corrected tip velocity and the blade spinning rate, the blade length can be computed as
bladelen_est = Vt_est/(bladerate_est*2*pi)
bladelen_est = 6.5285
Note that the result matches the simulation parameter of 6.5 meters. Information such as number of blades, blade length, and blade rotation rate can help identify the model of the helicopter.
This example introduces the basic concept of the micro-Doppler and shows its effect on the target return. It also shows how to extract micro-Doppler signature from the received I/Q signal and then derive relevant target parameters from the micro-Doppler information.
 Victor Chen, The Micro-Doppler Effect in Radar, Artech House, 2011
 Chen, Victor, et al. Micro-Doppler Effect in Radar: Phenomenon, Model, and Simulation Study, IEEE Trans. on Aerospace and Electronic Systems, Vol 42, No. 1, 2006
helicopmotion models the motion of multiple scatterers of the helicopter.
function [scatterpos,scattervel,scatterang] = 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); end