# Simulate a Maritime Radar PPI

This example shows how to simulate a plan position indicator (PPI) radar image for a rotating antenna array in a maritime environment. You will configure a radar scenario and spectral sea surface model, emulate a cargo ship with a cuboid extended target, generate return signals, and plot a PPI image of the result.

Set the RNG seed for repeatable results.

`rng default`

The scenario will consist of a large rotating uniform linear array (ULA) mounted in a fixed position above a sea surface. A large cuboid target will be placed on the sea surface.

Define the radar system parameters. Use an X-band radar at 10 GHz with a 1 degree azimuth beamwidth and 14 meter range resolution.

```freq = 10e9; % Hz azbw = 1; % deg rngres = 14; % m```

Calculate the number of uniformly-spaced elements required to achieve the desired beamwidth, and create a `phased.ULA`. Set the `BackBaffled` property of the array element to true so that only returns from the pointing direction are included.

```sinc3db = 0.8859; numelems = round(2*sinc3db/deg2rad(azbw)); lambda = freq2wavelen(freq); array = phased.ULA(numelems,lambda/2); array.Element.BackBaffled = true;```

Let the radar rotate at 50 rpm and express this in deg/sec. Calculate the rotation period from this rate.

```rotrate = 50/60*360; % deg/sec rotperiod = 360/rotrate; % sec```

Calculate the whole number of pulses required to minimally cover 360 degrees with the specified azimuth beamwidth. Use the rotation period to find the required pulse repetition frequency (PRF).

```npulses = ceil(360/azbw); prf = npulses/rotperiod; % Hz```

The fast-time sampling rate must be a whole multiple of the PRF. Use the desired range resolution to find the required sampling rate then adjust to match this constraint with the PRF. The radar will use a short pulse with one sample per resolution cell.

```sampleRate = 1/range2time(rngres); % Hz sampleRate = prf*round(sampleRate/prf); rngres = time2range(1/sampleRate); % m```

Now create the `radarTransceiver` and set the required parameters.

```rdr = radarTransceiver; rdr.TransmitAntenna.OperatingFrequency = freq; rdr.ReceiveAntenna.OperatingFrequency = freq; rdr.TransmitAntenna.Sensor = array; rdr.ReceiveAntenna.Sensor = array; rdr.Waveform.PRF = prf; rdr.Receiver.SampleRate = sampleRate; rdr.Waveform.SampleRate = sampleRate; rdr.Waveform.PulseWidth = 1/sampleRate;```

Create a `radarScenario` and set the `UpdateRate` to 0 to let the update rate be derived from the `radarTransceiver` configuration.

`scenario = radarScenario('UpdateRate',0);`

Now configure the sea surface. The surface will be defined on a square region of size 2 km-by-2 km.

`seaLength = 2e3; % m`

The `seaSpectrum` object is used to define the spectral model for the surface. The `Resolution` property defines the spacing of the underlying grid in the spatial domain, and the `WaveVectorSpacing` defines the spacing of the grid in the frequency domain. The resolution must be a factor of the length of the surface, so specify a desired surface resolution of 1/4th of the radar's range resolution (which is sufficient to capture the shape of the waves), then adjust to match this constraint. Set the `WaveVectorSpacing` to $2\pi \text{\hspace{0.17em}}$divided by the length of the surface.

```seaRes = rngres/4; % m seaRes = seaLength/round(seaLength/seaRes); spec = seaSpectrum('Resolution',seaRes,'WaveVectorSpacing',2*pi/seaLength);```

The `surfaceReflectivitySea` object is used to specify the reflectivity model that will be associated with the surface. Use the NRL reflectivity model with a sea state of 3 for horizontal polarization.

`refl = surfaceReflectivitySea('Model','NRL','SeaState',3,'Polarization','H');`

Create the 2 km-by-2 km sea surface using the reflectivity and spectral models. Specify the wind speed as 10 m/s, and the wind direction as 0 degrees so the wind is blowing in the +X direction.

```windSpeed = 10; % m/s windDir = 0; % deg bdry = [-1 1;-1 1]*seaLength/2; seaSurface(scenario,'Boundary',bdry,'RadarReflectivity',refl,'SpectralModel',spec,'WindSpeed',windSpeed,'WindDirection',windDir);```

Create the radar platform. Use the `kinematicTrajectory` object to mount the radar 24 meters above sea level, and use the `AngularVelocity` parameter to specify the rotation rate about the Z axis in radians per second.

```rdrHeight = 24; % m traj = kinematicTrajectory('Position',[0 0 rdrHeight],'AngularVelocitySource','Property','AngularVelocity',[0 0 deg2rad(rotrate)]); platform(scenario,'Sensors',rdr,'Trajectory',traj);```

Now define the cuboid target. Specify the target dimensions, total RCS, position, heading, and speed.

```tgtdims = struct('Length',120,'Width',18,'Height',22,'OriginOffset',[0 0 0]); % m tgtrcs = 40; % dBsm tgtpos = [seaLength/3 seaLength/16 5]; tgthdg = -30; % deg tgtspd = 8; % m/s```

Form the target velocity vector, create a kinematic trajectory, and add the target platform to the scene.

```tgtvel = tgtspd*[cosd(tgthdg) sind(tgthdg) 0]; tgttraj = kinematicTrajectory('Position',tgtpos,'Velocity',tgtvel,'Orientation',rotz(tgthdg).'); platform(scenario,'Trajectory',tgttraj,'Signatures',rcsSignature('Pattern',tgtrcs),'Dimensions',tgtdims);```

Finally, enable clutter generation for the radar by calling the `clutterGenerator` method on the scenario. The 3 dB beamwidth of the beam is used for clutter generation by default. In order to capture more of the surface return, disable this by setting `UseBeam` to false, and instead use a `RingClutterRegion` with an azimuth coverage 16 times that of the beam to encompass some sidelobes. Set the clutter resolution to half the range resolution. Returns will only be considered from a minimum ground range of 200 meters out to the maximum ground range of the surface. The azimuth center can be 0 to indicate that the region is centered about the +X direction. The `ringClutterRegion` method returns a handle to this region so that the parameters can be modified during the simulation.

```clut = clutterGenerator(scenario,rdr,'UseBeam',false,'Resolution',rngres/2); minGndRng = 200; maxGndRng = seaLength/2; azCov = 16*azbw; azCen = 0; reg = ringClutterRegion(clut,minGndRng,maxGndRng,azCov,azCen);```

Use the provided helper function to create a visualization of the scenario.

`helperPlotScenarioGeometry(rdrHeight,seaLength,array,freq,tgtpos)`

### Run the Simulation and Collect Returns

Each frame of the simulation will generate a range profile for one azimuth pointing direction. After receiving the raw IQ data on each frame, a `phased.RangeResponse` object will be used to perform matched filtering. Create the response object now, specifying the sample rate used, and get the required matched filter coefficients from the radar's waveform object.

```resp = phased.RangeResponse('RangeMethod','Matched filter','SampleRate',sampleRate); mfc = getMatchedFilter(rdr.Waveform);```

Specify how much of the full 360 degrees you would like to simulate. The scan starts at 0 degrees azimuth, and 45 degrees of coverage is enough to see the target. The full 360 degrees takes about 10 minutes on a machine with 64 GB of RAM and a 3.6 GHz CPU.

`azcov = 45;`

Set the scenario stop time based on the desired azimuth coverage.

`scenario.StopTime = azcov/360*npulses/prf;`

Run the simulation. Keep track of the frame number of each loop, and get the vector of range bins from the range response object. The matrix `ppi` will contain the signal data formatted as range-by-azimuth.

```frame = 0; nrng = floor(time2range(1/prf)/rngres); ppi = zeros(nrng,npulses); while advance(scenario) frame = frame + 1; reg.AzimuthCenter = (frame-1)*rotrate/prf; iqsig = receive(scenario); [ppi(:,frame),rngbins] = resp(sum(iqsig{1},2),mfc); end```

### Create a PPI Image

A PPI image consists of a set of range profiles arranged in radial lines to form a circular image of the scene in Cartesian space.

The IQ data covers the full range ambiguity, so start by trimming the data to the ranges of interest. Use the min and max ground range defined earlier, then find the range gate indices that correspond to those ranges.

```minRange = sqrt(minGndRng^2 + rdrHeight^2); maxRange = sqrt(maxGndRng^2 + rdrHeight^2); minGate = find(rngbins >= minRange,1,'first'); maxGate = find(rngbins <= maxRange,1,'last');```

Convert the azimuth pointing angles and range bins to rectangular coordinates, then use the `surf` function to plot the image. When specifying the azimuth domain for the image, use one more point than the number of pulses so that the image fully encompasses all azimuth angles.

```az = rotrate*(0:frame)/prf; gndrng = sqrt(rngbins(minGate:maxGate).^2-rdrHeight^2); x = gndrng.*cosd(az); y = gndrng.*sind(az); surf(x,y,zeros(size(x)),20*log10(abs(ppi(minGate:maxGate,1:frame)))) shading flat view(0,90) colorbar axis equal axis tight colormap winter title('PPI Image')```

The gif below shows a recording of this scenario for 30 full rotations of the antenna over about 34 seconds.

### Conclusion

In this example, you saw how to generate clutter and target returns with a rotating radar in a maritime environment. You saw how to use a spectral model to get realistic sea heights and surface shadowing, and how to emulate an extended target with a set of point targets, which allows for partial occlusion of the target by the surface. The IQ data was converted from polar format to Cartesian and plotted with the color channel of the `surf` function to create a simple PPI image.

### Supporting Functions

#### helperPlotScenarioGeometry

```function helperPlotScenarioGeometry(rdrHeight,seaLength,array,freq,tgtpos) % Plot a visualization of the scenario fh = figure; % Surface t = 0:10:360; r = linspace(0,seaLength/2,10).'; x = r.*cosd(t); y = r.*sind(t); surf(x/1e3,y/1e3,zeros(size(x)),repmat(permute([0 1 1],[1 3 2]),[size(x) 1]),'EdgeColor','none') hold on % Antenna pattern az = linspace(-6,6,80); maxEl = -atand(2*rdrHeight/seaLength); el = linspace(-90,maxEl,80); G = pattern(array,freq,az,el); [az,el] = meshgrid(az,el); [x,y,~] = sph2cart(az*pi/180,el*pi/180,1); rtg = rdrHeight./sind(-el); surf(x.*rtg/1e3,y.*rtg/1e3,1e-2*ones(size(x)),G,'EdgeColor','none') % Radar plot3(0,0,rdrHeight,'o') line([0 0],[0 0],[0 rdrHeight],'color','black','linestyle','--') % Target plot3(tgtpos(1)/1e3,tgtpos(2)/1e3,tgtpos(3),'*r') hold off text(0,0,rdrHeight*1.1,'Radar') text(tgtpos(1)*.8/1e3,tgtpos(2)*1.1/1e3,rdrHeight/10,'Target') text(seaLength/4/1e3,-200/1e3,rdrHeight/10,'Beam Pattern') text(-seaLength/4/1e3,seaLength/4/1e3,rdrHeight/10,'Sea Surface') xlabel('X (km)') ylabel('Y (km)') zlabel('Z (m)') title('Scenario Geometry') cl = clim; clim([cl(2)-30 cl(2)]) fh.Position = fh.Position + [0 0 150 150]; view([-11 58]) end```