Configure Scenario for Clutter Generation

Configuration of a radar scenario to simulate surface clutter involves creating a radarScenario object, adding platforms with mounted radars, adding surface objects that define the physical properties of the scenario surface, and enabling clutter generation for a specific radar in the scene.

Select a Radar Model

Clutter generation is performed as part of the scenario detect and receive methods. These methods are used to generate simulated radar detections and IQ signals, respectively. For detections use the radarDataGenerator and for IQ signals use the radarTransceiver.

Define some typical medium-PRF pulse-Doppler parameters for a side-looking airborne radar. Use a 15 degree depression angle so the radar is pointed towards the surface. Use a -90 degree mounting yaw angle so the radar faces out to the right of the platform. The mounting roll angle can be 0 to indicate no rotation about the antenna boresight. Use a center frequency of 1 GHz, range resolution of 80 m, and an 8 kHz PRF with 128 pulses per coherent processing interval (CPI). The beam is a symmetric cone shape, spanning 10 degrees in azimuth and elevation.

mountingYPR = [-90 15 0];
fc = 1e9;
rangeRes = 80;
prf = 8e3;
numPulses = 128;
bw = 10;

Construct a radarDataGenerator from these parameters. The radar updates once per CPI. The mounting angles point the radar in a broadside direction. Let the field of view and angular resolution equal the beamwidth. Set DetectionCoordinates to 'Scenario' to output detections in scenario coordinates for easier inspection of the results. Calculate the unambiguous range and radial speed and enable range and range-rate ambiguities for the radar. The ambiguities can be calculated with the time2range and dop2speed functions.

c = physconst('lightspeed');
lambda = freq2wavelen(fc);
rangeRateRes = lambda/2*prf/numPulses;
unambRange = time2range(1/prf);
unambRadialSpd = dop2speed(prf/4,lambda);
cpiTime = numPulses/prf;
rdr = radarDataGenerator(1,'No scanning','UpdateRate',1/cpiTime,'MountingAngles',mountingYPR,...
    'DetectionCoordinates','Scenario','HasINS',true,'HasElevation',true,'HasFalseAlarms',false,'HasRangeRate',true,...
    'HasRangeAmbiguities',true,'HasRangeRateAmbiguities',true,'MaxUnambiguousRadialSpeed',unambRadialSpd,...
    'MaxUnambiguousRange',unambRange,'CenterFrequency',fc,'FieldOfView',[bw bw],...
    'AzimuthResolution',bw,'ElevationResolution',bw,'RangeResolution',rangeRes,'RangeRateResolution',rangeRateRes);

Create a Scenario

The radarScenario object is the top-level manager for the simulation. A radar scenario may be Earth-centered, where the WGS84 Earth model is used, or it may be flat. Use a flat-Earth scenario to enable use of simple kinematic trajectories for the platforms. Set UpdateRate to 0 to let the scenario derive an update rate from the objects in the scenario.

scene = radarScenario('UpdateRate',0,'IsEarthCentered',false);

Add a radar platform to the scenario. Use a straight-line kinematic trajectory starting 1 km up, and moving in the +Y direction at 200 m/s, at a dive angle of 10 degrees. Orient the platform so the platform-frame +X direction is the direction of motion. Use the Sensors property to mount the radar you created earlier.

initialHeight = 1e3;
diveAngle = 10;
spd = 200;
vel = spd*[0 cosd(diveAngle) -sind(diveAngle)];
traj = kinematicTrajectory('Position',[0 0 initialHeight],'Velocity',vel,'Orientation',rotz(90).');
rdrplat = platform(scene,'Sensors',rdr,'Trajectory',traj);

Define the Scenario Surface

Physical properties of the scenario surface can be specified by using the scenario landSurface and seaSurface methods to define regions of land and sea surface types. Each surface added to the scene is a rectangular region with an associated radar reflectivity model and reference height. Land surfaces may additionally have associated static terrain data, and sea surfaces may have an associated spectral motion model. If no terrain data or spectral model is used, a surface may be unbounded, allowing for homogeneous clutter generation.

Create a simple unbounded land surface with a constant-gamma reflectivity model. Use the surfaceReflectivityLand function to create the reflectivity model. Use a gamma value of -10 dB.

refl = surfaceReflectivityLand('Model','ConstantGamma','Gamma',-10);
srf = landSurface(scene,'RadarReflectivity',refl)

srf = 
  LandSurface with properties:

    RadarReflectivity: [1×1 surfaceReflectivityLand]
      ReflectivityMap: 1
      ReferenceHeight: 0
             Boundary: [2×2 double]
              Terrain: []

In addition to the RadarReflectivity property, this surface object has a ReferenceHeight, which gives the constant height of the surface when no terrain is specified. The ReflectivityMap property is relevant only when a custom reflectivity model is used. The Boundary property gives the rectangular boundary of the surface in two-point form. Elements of Boundary can be +/-inf to indicate the surface is unbounded in one or more directions. Check the boundary of the surface created above to see that it is unbounded in all directions.

srf.Boundary

ans = 2×2

  -Inf   Inf
  -Inf   Inf

Access the SurfaceManager property of the scenario to see the surface objects that have been added, as well as any additional options related to the configuration of the scenario surface.

scene.SurfaceManager

ans = 
  SurfaceManager with properties:

    UseOcclusion: 1
        Surfaces: [1×1 radar.scenario.LandSurface]

The UseOcclusion property can be set false to disable line-of-sight occlusions by the surface, such as by terrain.

Enable Clutter Generation

Clutter Generator

Monostatic clutter generation is enabled for a specific radar by using the scenario clutterGenerator method. This method accepts parameter name-value pairs to set clutter generation properties. Clutter simulation configuration is performed on a radar-by-radar basis so that multiple radars can be simulated simultaneously, each with clutter generation settings appropriate for that radar.

The Resolution property specifies the nominal spacing of clutter patches for clutter return from scenario surfaces with no associated terrain or spectral model data (homogeneous surfaces). When clutter return is generated from surfaces with terrain data or spectral models, the spacing of the terrain data grid or spectral model resolution is used instead. The RangeLimit property is used to place an upper bound on the range of clutter generation. The UseBeam property is a logical scalar indicating whether or not the beam footprint clutter region should be used. Finally, the UseShadowing property is a logical scalar used to enable/disable shadowing (surface self-occlusion). Shadowing is only relevant for surfaces with terrain data or a spectral model.

Create a ClutterGenerator object, enabling clutter generation for the radar created above. Use a Resolution of 30 meters in order to get a couple clutter patches per range sample, and use a RangeLimit of 20 km to encompass the entire beam. The clutterGenerator method will return a handle to the ClutterGenerator object. Set UseBeam to true to enable automatic clutter generation within the field of view of the radar.

clut = clutterGenerator(scene,rdr,'Resolution',30,'RangeLimit',20e3,'UseBeam',true)

clut = 
  ClutterGenerator with properties:

      Resolution: 30
      RangeLimit: 20000
         UseBeam: 1
    UseShadowing: 1
         Regions: [1×0 radar.scenario.RingClutterRegion]
           Radar: [1×1 radarDataGenerator]

The clutter generator has two read-only properties. The Radar property stores a handle to the associated radar object, which was passed to the clutterGenerator method. The Regions property contains the set of custom clutter regions. The section below discusses clutter regions in detail.

Clutter Regions

Surface clutter is distributed over the entire range interval starting from the radar altitude and extending to the horizon range (or $$\infty$$ if using a flat-Earth scenario). It is distributed in elevation angle from -90 degrees up to the horizon elevation angle (or 0 if using a flat-Earth scenario), and over all 360 degrees of azimuth. Finally, surface clutter is distributed in Doppler as a result of platform motion.

There are two options available to specify regions on the surface from which to generate clutter return. The first is automatic mainlobe clutter, which generates clutter inside the footprint of the mainlobe of the radar antenna. This is only possible when the radar being used has a well-defined conical beam with azimuth and elevation two-sided beamwidths less than 180 degrees. For the radarDataGenerator, the width of the beam is defined by the FieldOfView property. For the radarTransceiver, a 3 dB beamwidth estimate is used.

The second option is to use the ringClutterRegion method of the clutter generator to specify a ring-shaped region of the scenario surface within which clutter is to be generated. This type of region is defined by a minimum and maximum ground range (relative to the radar nadir point) and an extent and center angle in azimuth. This region type is useful for capturing sidelobe and backlobe clutter return, mainlobe clutter return outside the 3 dB width, or to generate clutter from any other region of interest, such as at the location of a target platform.

The figure below illustrates these two region types. The beam footprint region is displayed as a magenta ellipse where the beam intersects the ground. Two ring regions are shown, one directly underneath the radar for capturing altitude return, and another for capturing some backlobe clutter return.

To demonstrate the utility of the customizable ring-shaped regions at capturing clutter return from arbitrary lobes of the antenna pattern, the above scenario is shown again below in a top-down view, with a gain pattern (such as from a uniform rectangular array) projected onto the ground. Note that a significant backlobe has been encompassed by a ring region. A circular region (which can be achieved by setting the minimum ground range to 0) can be used to capture altitude return. An additional region encompassing the mainlobe is also shown to demonstrate how this can be used to capture more of the mainlobe return, such as from null to null.

The radarDataGenerator is a statistics-based detectability simulator, and only simulates mainlobe detections within the field of view. As such, having UseBeam of the ClutterGenerator set to true is sufficient to capture the effect of clutter interference on the detectability of target platforms.

Visualize and Run Scenario

Theater Plotter

The theaterPlot object can be used along with a variety of theater plotters to create customizable visual representations of the scenario. Start by creating the theater plot.

tp = theaterPlot;

Now create plotter objects for the scenario surface, clutter regions, and resulting radar detections. The values specified for the DisplayName properties are used for the legend entries.

surfPlotter = surfacePlotter(tp,'DisplayName','Scenario Surface');
clutPlotter = clutterRegionPlotter(tp,'DisplayName','Clutter Region');
detPlotter = detectionPlotter(tp,'DisplayName','Radar Detections','Marker','.','MarkerEdgeColor','magenta');

Now that the scenario, clutter generator, and plotters are configured, use the detect method on the scenario to simulate a single frame and collect detections.

dets = detect(scene);

Plot the clutter region, which in this case is simply the beam footprint, along with the detection positions. Since the land surface used here is unbounded, the plotSurface call should come last so that the surface plot extends over the appropriate axis limits. The clutterRegionData method on the clutter generator is used to get plot data for the clutter region plotter. Similarly, for the surface plotter, the surfacePlotterData method on the scenario surface manager is used.

plotClutterRegion(clutPlotter,clutterRegionData(clut))
detpos = cell2mat(cellfun(@(t) t.Measurement(1:3).',dets,'UniformOutput',0));
plotDetection(detPlotter,detpos)
plotSurface(surfPlotter,surfacePlotterData(scene.SurfaceManager))

Note that the radarDataGenerator adds noise to the positions of detections, and may return detections with positions that fall outside the beam footprint.

Simulate Clutter IQ Signals

Now, create a radarTransceiver with similar radar system parameters. The helper function helperMakeTransceiver is provided to create a transceiver with the desired system parameters. The resulting radarTransceiver will have thermal noise disabled, and use a phased.CustomAntennaElement to approximate a uniform rectangular array with the specified beamwidth, which is recommended to speed up simulations when only a sum beam is needed.

rdr = helperMakeTransceiver(bw,fc,rangeRes,prf);
rdr.MountingAngles = mountingYPR;
rdr.NumRepetitions = numPulses;

Re-create the same scenario, using this new radar model.

scene = radarScenario('UpdateRate',0,'IsEarthCentered',false);
platform(scene,'Sensors',rdr,'Trajectory',kinematicTrajectory('Position',[0 0 initialHeight],'Velocity',vel,'Orientation',rotz(90).'));
landSurface(scene,'RadarReflectivity',refl);

Enable clutter generation for the radar. This time, disable the beam footprint clutter region in favor of a custom ring-shaped region.

clutterGenerator(scene,rdr,'Resolution',30,'RangeLimit',20e3,'UseBeam',false);

Note that if the clutterGenerator method was called without any output argument, the handle to the constructed ClutterGenerator may still be found with the scenario getClutterGenerator method by passing in a handle to the associated radar.

clut = getClutterGenerator(scene,rdr);

Use the ringClutterRegion method of the ClutterGenerator to create a custom region for clutter generation. Calculate the minimum and maximum radius for the region such that the mainlobe is fully encompassed in range and elevation angle. The null-to-null beamwidth in elevation can be estimated as about twice the 3 dB beamwidth, with an extra degree of margin on either side. Similarly, let the azimuth span of the region be 2.5 times the azimuth 3 dB beamwidth. The azimuth center is 0 degrees since the beam is pointing along the +X direction in scenario coordinates.

minr = initialHeight/tand(mountingYPR(2) + bw + 1);
maxr = initialHeight/tand(mountingYPR(2) - bw - 1);
azspan = bw*2.5;
azc = 0;
ringClutterRegion(clut,minr,maxr,azspan,azc)

ans = 
  RingClutterRegion with properties:

        MinRadius: 2.0503e+03
        MaxRadius: 1.4301e+04
      AzimuthSpan: 25
    AzimuthCenter: 0

Using the provided helper function, plot the ground-projected antenna pattern along with the ring clutter region you just created. The ring region created above nicely encompasses the entire mainlobe.

helperPlotGroundProjectedPattern(clut)

Run the simulation again for one frame, this time using the scenario receive method to simulate IQ signals.

iqsig = receive(scene);
PH = iqsig{1};

Since the radarTransceiver used a single custom element, the resulting signal will be formatted with fast-time samples along the first dimension, and pulse index (slow time) along the second dimension. This is the phase history matrix. Plot a DC-centered range-Doppler map (RDM) using the helperPlotRDM function.

helperPlotRDM(PH,rangeRes,prf,numPulses)

Re-create the theater plot visualization, this time just showing a representation of the surface and the ring-shaped clutter region.

tp = theaterPlot;
surfPlotter = surfacePlotter(tp,'DisplayName','Scenario Surface');
clutPlotter = clutterRegionPlotter(tp,'DisplayName','Clutter Region');
plotClutterRegion(clutPlotter,clutterRegionData(clut))
plotSurface(surfPlotter,surfacePlotterData(scene.SurfaceManager))

Clutter from Terrain Data

In the previous sections, you simulated homogeneous clutter from an unbounded flat surface. In this section, you will use a DTED file to simulate clutter return from real terrain data in an Earth-centered scenario. You will collect two frames of clutter return - one with shadowing enabled and one without shadowing, and compare the results.

Start by creating the scenario, this time setting the IsEarthCentered flag to true in order to use a DTED file, which consists of surface height samples over a latitude/longitude grid.

scene = radarScenario('UpdateRate',0,'IsEarthCentered',true);

Again use the landSurface method, passing in the name of the desired DTED file for the value of the Terrain parameter. The Boundary parameter can be used to constrain the domain of the added data. In general, as little of the terrain data should be loaded as needed for the specific application. In this case, use a 0.1-by-0.1 degree section of the example DTED data. For simplicity, continue using the same constant-gamma reflectivity model. Use an output argument to get a handle to the created surface object.

bdry = [39.5 39.6;-105.51 -105.41];
srf = landSurface(scene,'Terrain','n39_w106_3arc_v2.dt1','Boundary',bdry,'RadarReflectivity',refl);

When creating the radar platform, use the surface height method to place the platform 1 km above the surface. Earth-centered scenarios require trajectory information be specified with waypoints in latitude/longitude/altitude (LLA) format using the geoTrajectory object. The enu2lla function can be used to find the second waypoint in LLA format that results from the same velocity vector used in previous sections. Set the ReferenceFrame property of the trajectory to 'ENU' and use the same orientation as earlier so that our platform is flying North, with the radar pointing East and down at the specified depression angle.

waypoints = [39.53 -105.5];
srfHeight = height(srf,waypoints);
waypoints(1,3) = srfHeight + initialHeight;
toa = [0;1];
waypoints(2,:) = enu2lla(vel*diff(toa),waypoints(1,:),'ellipsoid');
traj = geoTrajectory('Waypoints',waypoints,'TimeOfArrival',toa,'Orientation',repmat(rotz(90).',1,1,2),'ReferenceFrame','ENU');
platform(scene,'Sensors',rdr,'Trajectory',traj);

Create the clutter generator and use the same ring-shaped region as before. The Resolution property is irrelevant here since clutter patches will be generated at each facet of the terrain data. Surface shadowing is enabled by default.

clut = clutterGenerator(scene,rdr,'RangeLimit',20e3,'UseBeam',false);
ringClutterRegion(clut,minr,maxr,azspan,azc);

Simulate clutter return for one frame and save the resulting phase history.

iqsig = receive(scene);
PH_withShadowing = iqsig{1};

The theater plot does not support visualizations of Earth-centered scenarios, so use the provided helper function to show a scenario overview. This time, the terrain is plotted along with the radar frame and clutter patches. Notice the shadowed region can be seen in the plot as a gap in the clutter patches (cyan). Those patches are shadowed by other parts of the terrain that are closer to the radar.

helperPlotClutterScenario(scene)
title('Clutter Patches - with terrain shadowing')

To see the difference when shadowing is not used, turn shadowing off, run the simulation for another frame, and save the result.

clut.UseShadowing = false;
iqsig = receive(scene);
PH_noShadowing = iqsig{1};

Plot the scenario overview and notice that some of the clutter region has been filled in with clutter patches. Notice there are still some gaps visible. This is because clutter patches that are facing away from the radar, such as if they are on the other side of a hill, will never be visible, regardless of the value of the UseShadowing property.

helperPlotClutterScenario(scene)
title('Clutter Patches - without terrain shadowing')

Now to see the effect of shadowing in the RDM, plot the clutter return signals recorded previously side-by-side.

subplot(1,2,1)
helperPlotRDM(PH_withShadowing,rangeRes,prf,numPulses)
title('RDM - with terrain shadowing')
subplot(1,2,2)
helperPlotRDM(PH_noShadowing,rangeRes,prf,numPulses)
title('RDM - without terrain shadowing')
set(gcf,'Position',get(gcf,'Position')+[0 0 560 0])

In the first plot there is a noticeable section around 3.4 km range with no clutter return. This is a result of surface shadowing. In the second plot, that region has been partially filled in.

Shadowing is an important phenomenon when simulating clutter return from real surfaces, but it may be disabled if not needed.

Simulate Clutter Range Profile

The automatic mainlobe clutter return supports scanning radars. In this section you will re-create the scenario to use a stationary scanning radarTransceiver that collects a single pulse per scan position, and view the resulting range profiles.

Start by re-creating the radar object and configuring for an electronic sector scan. Set number of pulses per scan position to 1. Set the scanning limits to cover 30 degrees in azimuth.

numPulses = 1;
rdr = helperMakeTransceiver(bw,fc,rangeRes,prf);
rdr.MountingAngles = mountingYPR;
rdr.NumRepetitions = numPulses;
rdr.ElectronicScanMode = 'Sector';
rdr.ElectronicScanLimits = [-15 15;0 0];
rdr.ElectronicScanRate = [prf; 0];

Re-create the scenario and platform. Set the scenario stop time to run about 1 full scan. Use a homogeneous unbounded land surface with the same reflectivity model as used earlier.

scene = radarScenario('UpdateRate',0,'IsEarthCentered',false,'StopTime',30/prf);
platform(scene,'Sensors',rdr,'Position',[0 0 initialHeight],'Orientation',rotz(90).');
landSurface(scene,'RadarReflectivity',refl);

Enable clutter generation, using only the 3 dB beam footprint for clutter generation.

clutterGenerator(scene,rdr,'Resolution',30,'RangeLimit',10e3,'UseBeam',true);

Run the simulation, collecting the range profile at each scan position. Use the info output from receive to record the look angle used by the radar on each frame.

rangeGates = 0:rangeRes:(unambRange-rangeRes);
frame = 0;
while advance(scene)
    frame = frame + 1;
    [iqsig,info] = receive(scene);
    lookAng(:,frame) = info.ElectronicAngle;
    rangeProfiles(:,frame) = 20*log10(abs(iqsig{1}));
end

The gif below shows the scenario overview with clutter patches and the range profile collected on each frame.

Plot the range profiles against range and azimuth scan angle.

figure
imagesc(lookAng(1,:),rangeGates/1e3,rangeProfiles);
set(gca,'ydir','normal')
xlabel('Azimuth Scan Angle (deg)')
ylabel('Range (km)')
title('Clutter Range Profiles (dBW)')
colorbar

The extent of clutter in range is greatest at +/-15 degree azimuth, and least at 0 degrees. The returned clutter power is greatest where the grazing angle is greatest, which occurs at 0 degrees azimuth.

Conclusion

In this example, you saw how to configure a radar scenario to include clutter return as part of the detect and receive methods, generating clutter detections and IQ signals with the radarDataGenerator and radarTransceiver, respectively. You saw how to define a region of the scenario surface with an associated reflectivity model, and how to specify regions of interest for clutter generation. Surface shadowing is simulated when generating clutter returns from surfaces with terrain.

Supporting Functions

helperMakeTransceiver

function rdr = helperMakeTransceiver( bw,fc,rangeRes,prf )
% This helper function creates a radarTransceiver from some basic system
% parameters.

c = physconst('lightspeed');

rdr = radarTransceiver;

rdr.TransmitAntenna.OperatingFrequency = fc;
rdr.ReceiveAntenna.OperatingFrequency = fc;

rdr.Waveform.PRF = prf;

sampleRate = c/(2*rangeRes);
sampleRate = prf*round(sampleRate/prf); % adjust to match constraint with PRF

rdr.Receiver.SampleRate = sampleRate;
rdr.Waveform.SampleRate = sampleRate;

rdr.Waveform.PulseWidth = 2*rangeRes/c;

if isempty(bw)
    % Use an isotropic element
    rdr.TransmitAntenna.Sensor = phased.IsotropicAntennaElement;
    rdr.ReceiveAntenna.Sensor = phased.IsotropicAntennaElement;
else
    % Create a custom element from the pattern of a URA
    N = round(0.8859*2./(bw(:).'*pi/180));
    N = flip(N);
    lambda = freq2wavelen(fc,c);
    array = phased.URA(N,lambda/2);
    az = -180:.4:180;
    el = -90:.4:90;
    G = pattern(array,fc,az,el,'Type','efield','normalize',false);
    M = 20*log10(abs(G));
    P = angle(G);
    E = phased.CustomAntennaElement('FrequencyVector',[fc-1 fc+1],...
        'AzimuthAngles',az,'ElevationAngles',el,'MagnitudePattern',M,'PhasePattern',P);
    rdr.TransmitAntenna.Sensor = E;
    rdr.ReceiveAntenna.Sensor = E;
end

% Turn off thermal noise
rdr.Receiver.NoiseMethod = 'Noise power';
rdr.Receiver.NoisePower = 0;

end

helperPlotGroundProjectedPattern

function helperPlotGroundProjectedPattern( clutter )
% Input a ClutterGenerator that has an associated radarTransceiver. Plots
% the clutter regions and ground-projected gain pattern. Assumes a flat
% infinite ground plane at Z=0.

Naz = 360*4;
Nel = 90*4;

% Force update the patch generator sensor data
clutter.PatchGenerator.updateSensorData(clutter.Radar,clutter.Platform,clutter.Scenario.SimulationTime);

pos = clutter.PatchGenerator.SensorData.Position;
fc = clutter.PatchGenerator.SensorData.CenterFrequency;
B = clutter.PatchGenerator.SensorData.SensorFrame;
maxGndRng = clutter.RangeLimit;

% Get azimuth/elevation grid in scenario coordinates
azScen = linspace(-180,180,Naz);
maxEl = -atand(pos(3)/maxGndRng);
elScen = linspace(-90,maxEl,Nel);

% Convert az/el to sensor frame
[losxScen,losyScen,loszScen] = sph2cart(azScen*pi/180,elScen.'*pi/180,1);
R = B.';
losx = R(1,1)*losxScen + R(1,2)*losyScen + R(1,3)*loszScen;
losy = R(2,1)*losxScen + R(2,2)*losyScen + R(2,3)*loszScen;
losz = R(3,1)*losxScen + R(3,2)*losyScen + R(3,3)*loszScen;
[az,el,~] = cart2sph(losx,losy,losz);
az = az*180/pi;
el = el*180/pi;

% Get gain pattern
sensor = clutter.Radar.TransmitAntenna.Sensor;
G = sensor(fc,[az(:) el(:)].');
G = reshape(G,size(az));
G = G/max(G(:));

rtg = -pos(3)./loszScen;
surf(losxScen.*rtg,losyScen.*rtg,zeros(size(az)),20*log10(abs(G)))
axis equal
shading flat
view(0,90)
clim([-60 0])

hold on
for reg = clutter.Regions
    helperPlotClutterRegion(reg,pos);
end
hold off

end

helperPlotRDM

function helperPlotRDM( PH,rangeRes,prf,numPulses )
% This helper function forms and plots an RDM from a phase history matrix

c = physconst('lightspeed');

% Form DC-centered RDM and convert to dBW
RDM = fftshift(fft(PH,[],2),2);
RDM = 20*log10(abs(RDM));

% Range and Doppler bins
rngBins = 0:rangeRes:c/(2*prf);
dopBins = -prf/2:prf/numPulses:prf/2-prf/numPulses;

% Plot
imagesc(dopBins,rngBins/1e3,RDM);
set(gca,'ydir','normal')
xlabel('Doppler (Hz)')
ylabel('Range (km)')
mx = max(RDM(:));
if ~isinf(mx)
    clim([mx-60 mx]);
end
colorbar

end