# 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 an extended target with a collection of point scatterers, 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 rectangular array (URA) mounted in a fixed position above a sea surface. A set of point targets arranged on the surface of a cuboid will be used to emulate a medium-sized container ship.

Define the radar system parameters. Use an X-band radar at 10 GHz with a 1 degree azimuth beamwidth and 10 meter range resolution. Use a wide 60 deg elevation beamwidth which will result in a URA with two rows of elements.

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

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

```sinc3db = 0.8859; numelems(1) = round(2*sinc3db/deg2rad(elbw)); numelems(2) = round(2*sinc3db/deg2rad(azbw)); lambda = freq2wavelen(freq); array = phased.URA(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);```

Use the `MountingAngles` property on the radar to point the beam down so that it is centered between the radar nadir point and the edge of the sea surface.

```depang = atand(4*rdrHeight/seaLength); rdr.MountingAngles = [0 depang 0];```

Define a cuboid target then use the provided helper function to add a set of discrete scatterers to the scene for a basic representation of the target. Specify the target dimensions, total RCS, position, heading, and speed. The specified range resolution is used to determine the spacing of the scatterers along the surface of the cuboid.

```tgtdims = [120 18 22]; % length, width, height (m) tgtrcs = 10; % m^2 tgtpos = [seaLength/3 seaLength/16 0]; tgthdg = -30; % deg tgtspd = 8; % m/s helperMakeExtendedTarget(scenario,tgtdims,tgtrcs,tgtpos,tgthdg,tgtspd,rngres)```

Finally, enable clutter generation for the radar by calling the `clutterGenerator` method on the scenario. The clutter generator will include only the region of the surface that is within the 3 dB width of the beam. Surface shadowing is enabled by default, and can be disabled with the `UseShadowing` property.

`clut = clutterGenerator(scenario,rdr);`

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

`helperPlotScenarioGeometry(rdrHeight,seaLength,depang,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. One half is subtracted from the total number of pulses to ensure the simulation includes the exact number of pulses specified.

`scenario.StopTime = azcov/360*(npulses-1/2)/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; 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 a minimum ground range of 200 m and a maximum ground range equal to half the length of the sea surface, then find the range gate indices that correspond to those ranges.

```minRange = sqrt(200^2 + rdrHeight^2); maxRange = sqrt((seaLength/2)^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 360 degrees of azimuth.

```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

#### helperMakeExtendedTarget

```function helperMakeExtendedTarget(scene,dims,totrcs,pos,hdg,spd,rngres) % Add a set of point targets to a scenario to represent a cuboid target scatSpacing = rngres*.8; nScatsPerDim = ceil(dims/scatSpacing); platB = rotz(hdg); losrdr = -pos(:); losrdr = platB.'*losrdr; if losrdr(1) < -dims(1)/2 if losrdr(2) < -dims(2)/2 sideFaces = 'ws'; elseif losrdr(2) < dims(2)/2 sideFaces = 'w'; else sideFaces = 'wn'; end elseif losrdr(1) < dims(1)/2 if losrdr(2) < -dims(2)/2 sideFaces = 's'; elseif losrdr(2) < dims(2)/2 error('Radar is inside target') else sideFaces = 'n'; end else if losrdr(2) < -dims(2)/2 sideFaces = 'es'; elseif losrdr(2) < dims(2)/2 sideFaces = 'e'; else sideFaces = 'en'; end end xd = linspace(-dims(1)/2,dims(1)/2,nScatsPerDim(1)); yd = linspace(-dims(2)/2,dims(2)/2,nScatsPerDim(2)); zd = linspace(-dims(3)/2,dims(3)/2,nScatsPerDim(3)); x = []; y = []; z = []; for ind = 1:numel(sideFaces) switch sideFaces(ind) case 'e' [y0,z0] = meshgrid(yd,zd); x0 = dims(1)/2*ones(size(y0)); case 'n' [x0,z0] = meshgrid(xd,zd); y0 = dims(2)/2*ones(size(x0)); case 'w' [y0,z0] = meshgrid(yd,zd); x0 = -dims(1)/2*ones(size(y0)); case 's' [x0,z0] = meshgrid(xd,zd); y0 = -dims(2)/2*ones(size(x0)); end x = [x, x0(:).']; y = [y, y0(:).']; z = [z, z0(:).']; end propRcsSide = 0.5; % proportion of the RCS to put on the sides of the target rcs = totrcs*propRcsSide/numel(x)*ones(size(x)); % top face [x0,y0] = meshgrid(xd,yd); z0 = dims(3)/2*ones(size(x0)); x = [x, x0(:).']; y = [y, y0(:).']; z = [z, z0(:).']; rcs = [rcs, totrcs*(1-propRcsSide)/numel(x0)*ones(1,numel(x0))]; % Rotate [x,y,z] = deal(platB(1,1)*x+platB(1,2)*y+platB(1,3)*z,... platB(2,1)*x+platB(2,2)*y+platB(2,3)*z,... platB(3,1)*x+platB(3,2)*y+platB(3,3)*z); % Translate x = x + pos(1); y = y + pos(2); z = z + pos(3) + dims(3)/4; for ind = 1:numel(x) traj = kinematicTrajectory('Position',[x(ind) y(ind) z(ind)],'Velocity',spd*[cosd(hdg) sind(hdg) 0]); platform(scene,'Trajectory',traj,'Signatures',rcsSignature('Pattern',10*log10(rcs(ind)))); end end```

#### helperPlotScenarioGeometry

```function helperPlotScenarioGeometry(rdrHeight,seaLength,depang,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+depang); [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```