Track Heterogeneous Aircraft using Heterogeneous Sensors
This example shows how to track heterogeneous aircraft using data collected by heterogeneous sensors. In the example, you track commercial passenger aircraft, general aviation aircraft, and helicopters flying in the airspace near the Logan International Airport in Boston, MA. You use active and passive radars and their outputs to track the aircraft using a JIPDATracker
.
Define Scenario
This example shows tracking of heterogeneous aircraft using heterogeneous sensors. The helper function helperCreateScenario
defines a scenario around Logan International Airport in Boston, MA metropolitan area. The scenario contains two ground-based sensors mounted on towers and seven aircraft.
scenario = helperCreateScenario;
The first sensor, which is mounted on the scenario platform number 1, is an airport primary surveillance radar, which is a monostatic, active, rotating radar. As it scans, the radar reports measurements that contain azimuth, elevation, range, and range rate. The second sensor is a passive secondary receiver, which scans the sky for signals emitted by the traveling aircraft. It measures the direction from which the signals arrive, in terms of azimuth and elevation angles.
The seven aircraft are listed in the table below.
ID | Type | Position | Speed |
1 | Primary surveillance radar tower | [0 0 -50] | Stationary |
2 | Secondary radar tower | [-5000 -3000 -50] | Stationary |
3 | Passenger aircraft | Traveling east towards Logan at 3000m altitude | 900 km/hr |
4 | Passenger aircraft | Traveling north away from Logan at 4000m altitude | 700 km/hr |
5 | Passenger aircraft | Traveling east tangential to Logan at 4000m altitude | 600 km/hr |
6 | General aviation aircraft | Traveling south tangential to Logan at 2000m altitude | 300 km/hr |
7 | General aviation aircraft | Traveling east tangential to Logan at 1000m altitude | 200 km/hr |
8 | Helicopter | Traveling north climbing from 700m to 1000m. | 50 km/hr |
9 | Helicopter | Traveling east descending from 600m to 320m | 10 km/hr |
The following code creates a trackingGlobeViewer
and visualizes the scenario.
mapOrigin = [42.39423231362 -70.95934958874 0]; mapViewer = trackingGlobeViewer('ReferenceLocation',mapOrigin,'Basemap','streets-dark'); campos(mapViewer, mapOrigin + [0 0 1e5]); drawnow; plotScenario(mapViewer,scenario); snapshot(mapViewer);
Define Tracker
Tracking heterogeneous aircraft using heterogeneous sensors is not trivial. You need to consider the varying data formats and sensor capabilities. Additionally, the system must be robust enough to handle real-time data integration and ensure accurate tracking across different aircraft types.
The tracker must be a multi-sensor data fusion algorithm capable of integrating inputs from all available sensors to provide a unified and accurate representation of each aircraft's state.
As new sensor data become available, the tracker takes into consideration existing tracks it maintains, associates the new sensor data to the existing tracks, updates the tracks using common filtering techniques, like Kalman filtering, and manages the track existence probability, taking into account the sensor instantaneous field of view. All these steps are enabled by defining the target and sensor specifications, which describe how the targets move and how the sensors observe the scenario. Therefore, the steps to define the tracker are shown in the following process:
Specify Aircraft Types
The first step in defining the tracker requires specifying the types of objects you want to track. In this example, these objects are passenger aircraft, general aviation aircraft, and helicopters. Use the trackerTargetSpec
function to create a target specification for a passenger aircraft and observe the properties of the specification.
passengerSpec = trackerTargetSpec("aerospace","aircraft","passenger"); disp(passengerSpec)
PassengerAircraft with properties: MaxHorizontalSpeed: 250 m/s MaxVerticalSpeed: 20 m/s MaxHorizontalAcceleration: 10 m/s² MaxVerticalAcceleration: 1 m/s²
The passenger aircraft target specification defines typical speeds and accelerations in the horizontal and vertical planes for a passenger aircraft. The speeds are used to define a prior distribution for a constant velocity model while the accelerations are used to define the process noise of the constant velocity model. This is all taken care of by the spec, when it integrates with the tracker shown below. In the passenger aircraft case, the speeds and accelerations represent a typical passenger jet flying at common cruise speeds. You can use this specification for tracking any passenger jet-propelled aircraft, including business jets. The speeds and accelerations are tunable and can be modified to fit specific aircraft type.
To add the general aviation and helicopter specifications, use the same process. Observe the differences between the relatively fast jet passenger aircraft, the slower general aviation fixed-wing aircraft, and the even slower helicopter aircraft.
generalAviationSpec = trackerTargetSpec("aerospace","aircraft","general-aviation")
generalAviationSpec = GeneralAviation with properties: MaxHorizontalSpeed: 150 m/s MaxVerticalSpeed: 15 m/s MaxHorizontalAcceleration: 1 m/s² MaxVerticalAcceleration: 5 m/s²
helicopterSpec = trackerTargetSpec("aerospace","aircraft","helicopter")
helicopterSpec = Helicopter with properties: MaxHorizontalSpeed: 80 m/s MaxVerticalSpeed: 10 m/s MaxHorizontalAcceleration: 5 m/s² MaxVerticalAcceleration: 1 m/s²
Specify Sensor Types
After specifying what you want to track, you must specify the sensors you use for tracking. In this example, you use the airport primary surveillance radar and a secondary radar. The airport surveillance radar is a monostatic rotating radar that scans the sky as it rotates. The secondary radar is a passive radar that looks for RF returns from the tracked objects. First, use the trackerSensorSpec
to create the monostatic radar specification for the airport surveillance radar and observe its properties.
activeRadarSpec = trackerSensorSpec("aerospace","radar","monostatic")
activeRadarSpec = AerospaceMonostaticRadar with properties: MaxNumLooksPerUpdate: 30 MaxNumMeasurementsPerUpdate: 10 IsPlatformStationary: 1 PlatformPosition: [0 0 0] m PlatformOrientation: [3⨯3 double] MountingLocation: [0 0 0] m MountingAngles: [0 0 0] deg HasElevation: 1 HasRangeRate: 1 FieldOfView: [60 20] deg RangeLimits: [0 1e+05] m RangeRateLimits: [-500 500] m/s AzimuthResolution: 1 deg RangeResolution: 100 m ElevationResolution: 5 deg RangeRateResolution: 10 m/s DetectionProbability: 0.9 FalseAlarmRate: 1e-06
You can divide the radar properties into the following broad categories:
Properties related to how the radar operates and how it reports to the tracker, for example
MaxNumLooksPerUpdate,
IsPlatformStationary,
PlatformPosition
, etc. These properties must be specified using information from the radar operator, for example Logan airport.Properties related to the radar measurement information, for example
HasElevation,
FieldOfView,
RangeLimits,
etc. These properties are specified by the radar manufacturer for a specific radar mode in use and can be found on a radar spec sheet.
Use the following code to define the Logan airport surveillance radar. The radar reports its measurements to the tracker once every two seconds.
towerPosition = [0 0 -50];
activeRadarSpec.MaxNumLooksPerUpdate = ceil(75*2/1.4); % Defines the number of looks in two seconds
activeRadarSpec.MaxNumMeasurementsPerUpdate = 20;
activeRadarSpec.PlatformPosition = towerPosition;
activeRadarSpec.PlatformOrientation = eye(3);
You follow a similar process to define the scanning passive radar.
passiveRadarSpec = trackerSensorSpec("aerospace","radar","direction-finder"); passiveRadarSpec.MaxNumMeasurementsPerUpdate = 10; passiveRadarSpec.MaxNumLooksPerUpdate = ceil(150/1.4); passiveRadarSpec.PlatformPosition = [-5000 -3000 -50]; passiveRadarSpec.PlatformOrientation = eye(3); passiveRadarSpec.RangeLimits = [0 1e5];
The following two lines of code ensure that the properties of the sensor specifications match the equivalent properties of the scenario sensors.
activeRadarSpec = helperSyncSensor2spec(scenario.Platforms{1}.Sensors{1},activeRadarSpec); passiveRadarSpec = helperSyncSensor2spec(scenario.Platforms{2}.Sensors{1},passiveRadarSpec);
Configure Tracker to Use Target and Sensor Specifications
The last step in specifying the tracker is to configure it to use the heterogeneous target and sensor specifications you have defined.
tracker = multiSensorTargetTracker({passengerSpec,generalAviationSpec,helicopterSpec},{activeRadarSpec,passiveRadarSpec},"jipda");
Using the target specifications, the tracker configures the right motion models. Using the sensor specifications, the tracker configures the right observability, measurement, clutter, birth, and state initialization models. This is all done within the tracker and does not require additional tuning. Note that all the target and sensor specification properties are tunable.
disp(tracker)
fusion.tracker.JIPDATracker with properties: TargetSpecifications: {[1×1 PassengerAircraft] [1×1 GeneralAviation] [1×1 Helicopter]} SensorSpecifications: {[1×1 AerospaceMonostaticRadar] [1×1 AerospaceESMRadar]} MaxMahalanobisDistance: 5 ConfirmationExistenceProbability: 0.9000 DeletionExistenceProbability: 0.1000
In addition to the target and sensor specifications, the tracker exposes three tuning parameters to help control data association, track confirmation, and track deletion.
The MaxMahalanobisDistance property controls the size of the association gate. Generally, this distance should be small to limit the association gate, but if you observe that the tracker fails to associate a measurement to the track, you can increase this value slightly.
Any multi-object tracker needs to handle the possibility of false alarms reported by the sensors. As a result, the tracker creates tracks that are initially tentative, meaning they are not confirmed to be real objects. The ConfirmationExistenceProbability
property controls how quickly the tracker confirms a track. A higher value of this property means that the track must have a higher probability of existence to be confirmed as a real object. A lower value means that tracks with lower probability can be confirmed and it may lead to a higher rate of false tracks.
Similarly, any multi-object tracker needs to handle the possibility that an object may not exist anymore. The DeletionExistenceProbability
property defines how quickly such tracks are removed from the tracker. A higher value means more aggressive deletion and can help improve the tracker performance at the cost of a higher likelihood of erroneously deleting a real track.
Understand Sensor Data Format
The sensor specifications allow you to pass sensor data into the tracker and act as an interpreter between the sensor data and the tracker. To understand this point better, display the data format for the active radar spec.
dataFormat(activeRadarSpec)
ans = struct with fields:
LookTime: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
LookAzimuth: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
LookElevation: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
DetectionTime: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
Azimuth: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
Elevation: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
Range: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
RangeRate: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
AzimuthAccuracy: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
ElevationAccuracy: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
RangeAccuracy: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
RangeRateAccuracy: [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
The active radar data format includes all the information needed to understand the sensor scanning and measurements. The LookTime
, LookAzimuth
, and LookElevation
fields provide information about the sensor beam azimuth and elevation angles as a function of time. Each of these fields contains as many zeros as the MaxNumLooksPerUpdate
property of the active radar specification. When populating the struct with your sensor data, these fields are variable sized and you can pass only the number of elements that you need, for example, [0 0.1 0.2] for time.
The rest of the fields provide information about the sensor detections: their time, measurement, and measurement accuracy. The measurement and measurement accuracy are further split into azimuth, elevation, range, and range rate. Furthermore, the number of zeros in each of these fields is equal to the MaxNumMeasurementsPerUpdate
property and you use them as variable sized.
A similar pattern is used for the passive radar data format.
The next line loads sensor data that was recorded in the active and passive radar data formats. The first cell in the sensorData
cell array contains the monostatic sensor data from the beginning to the end of the scenario. Similarly, the second cell contains the passive radar data. To recreate the data, uncomment the line below, but note that it may take a considerable amount of time to run.
load("sensorData.mat","sensorData"); % sensorData = helperRecordSensorData(scenario);
Run Tracker
The following controls allow tuning the tracker.
release(tracker); tracker.DeletionExistenceProbability = 0.3; tracker.ConfirmationExistenceProbability = 0.9; tracker.MaxMahalanobisDistance = 5;
You can define how much time passes between tracker updates using the control below. The helperSplitSensorData
and helperSavePlatformPoses
functions split the sensor and truth poses data based on this interval.
trackerUpdateInterval = 1; % seconds [activeRadarData,passiveRadarData] = helperSplitSensorData(sensorData,scenario.StopTime,trackerUpdateInterval); platformPoses = helperSavePlatformPoses(scenario,trackerUpdateInterval);
The following code prepares the map display.
clear(mapViewer);
plotPlatform(mapViewer,[scenario.Platforms{3:end}],TrajectoryMode="Full");
Update the tracker and display in a loop.
for ind = 1:numel(activeRadarData) % Update tracker tracks = tracker(activeRadarData(ind), passiveRadarData(ind)); % Visualize data and tracks (helper) plotActiveRadarData(mapViewer, activeRadarSpec,activeRadarData(ind)); plotESMRadarData(mapViewer, passiveRadarSpec,passiveRadarData(ind)); % Visualize tracks and platforms plotTrack(mapViewer, tracks, LabelStyle="ATC",LineWidth=6); plotPlatform(mapViewer,platformPoses(:,ind)); end
Take a final snapshot.
snapshot(mapViewer);
Summary
In this example, you learned how to track heterogeneous aircraft using a diverse set of radars. You saw how to specify the target types and the sensor types using the trackerTargetSpec
and trackerSensorSpec
functions, respectively. These functions allow you to define target and sensor specifications that are common in the aerospace and automotive domains. You also saw how to configure the tracker to use these specifications with the multiSensorTargetTracker
function.
This example was set up in a way that allows you to explore the tracker capabilities. You can use the controls in the "Run the Tracker" section to modify the tracker parameters and update rate. The tracker parameters, including those parameters of each specification, are tunable during tracker's operation.
To explore further, you can modify the scenario. Ensure that you use the correct target and sensor specifications to match the new scenario, and update the sensor specifications accordingly.
Utility Functions
helperCreateScenario - A helper to create the scenario
function scenario = helperCreateScenario() scenarioDuration = 60; % s scenario = trackingScenario(UpdateRate=0,StopTime=scenarioDuration); radarTower = platform(scenario, Position=[0 0 -50]); rpm = 12.5; fov = [1.4;10]; scanrate = rpm*360/60; % deg/s updaterate = scanrate/fov(1); % Hz activeRadar = fusionRadarSensor(1,"Rotator", ... UpdateRate=updaterate, ... % Hz FieldOfView=fov, ... % [az;el] deg MaxAzimuthScanRate=scanrate, ... % deg/sec AzimuthResolution=fov(1), ... % deg ReferenceRange=111e3, ... % m ReferenceRCS=0, ... % dBsm RangeResolution=135, ... % m HasElevation=true, ... HasRangeRate=true, ... RangeRateLimits=[-300 300], ... MountingLocation=[0 0 -15], ... MountingAngles=[0 0 0], ... DetectionCoordinates="Sensor spherical"); % Report spherical measurements % Mount radar at the top of the tower radarTower.Sensors = activeRadar; % Set mechanical elevation scan to begin at 2 degrees above the horizon elFov = fov(2); tilt = 2; % deg activeRadar.MechanicalElevationLimits = [-fov(2) 0]-tilt; % deg activeRadar.FieldOfView(2) = elFov+1e-3; esmTower = platform(scenario, Position=[-5000 -3000 -50]); esmfov = [1.4 15]; passiveRadar = fusionRadarSensor(2,"Rotator", ... DetectionMode="ESM", ... MechanicalElevationLimits=[-esmfov(2) 0]-tilt, ... % deg UpdateRate=updaterate, ... % Hz FieldOfView=esmfov, ... % [az;el] deg MaxAzimuthScanRate=scanrate, ... % deg/sec AzimuthResolution=esmfov(1), ... % deg HasElevation=true, ... MountingLocation=[0 0 -15], ... % m MountingAngles=[-90 0 0], ... % deg. Make sure it rotates after the active radar Sensitivity=-200 ... % dB ); % Reports spherical measurements % Set mechanical elevation scan to begin at 2 degrees above the horizon elFov = esmfov(2); passiveRadar.FieldOfView(2) = elFov+1e-3; % Mount radar at the top of the tower esmTower.Sensors = passiveRadar; emitter = radarEmitter(1,"No scanning",UpdateRate=updaterate,FieldOfView=[360 180]); % Inbound passenger aircraft ht = 3e3; spd = 900*1e3/3600; % m/s wp = [-5e3 -40e3 -ht;-5e3 -40e3+spd*scenarioDuration -ht]; traj = waypointTrajectory(wp,[0 scenarioDuration]); emitter1 = clone(emitter); emitter1.EmitterIndex = 3; platform(scenario,Trajectory=traj,Emitters=emitter1); % Outbound passenger aircraft ht = 4e3; spd = 700*1e3/3600; % m/s wp = [20e3 10e3 -ht;20e3+spd*scenarioDuration 10e3 -ht]; traj = waypointTrajectory(wp,[0 scenarioDuration]); emitter2 = clone(emitter); emitter2.EmitterIndex = 4; platform(scenario,Trajectory=traj,Emitters=emitter2); % Tangential passenger aircraft ht = 4e3; spd = 600*1e3/3600; % m/s wp = [-20e3 -spd*scenarioDuration/2 -ht;-20e3 spd*scenarioDuration/2 -ht]; traj = waypointTrajectory(wp,[0 scenarioDuration]); emitter3 = clone(emitter); emitter3.EmitterIndex = 5; platform(scenario,Trajectory=traj,Emitters=emitter3); % General aviation aircraft ht = 2e3; spd = 300*1e3/3600; % m/s wp = [spd*scenarioDuration/2 10e3 -ht;-spd*scenarioDuration/2 10e3 -ht]; traj = waypointTrajectory(wp,[0 scenarioDuration]); emitter4 = clone(emitter); emitter4.EmitterIndex = 6; platform(scenario,Trajectory=traj,Emitters=emitter4); % General aviation aircraft ht = 1e3; spd = 200*1e3/3600; % m/s wp = [-12e3 -4e3-spd*scenarioDuration/2 -ht;-12e3 -4e3+spd*scenarioDuration/2 -ht]; traj = waypointTrajectory(wp,[0 scenarioDuration]); emitter5 = clone(emitter); emitter5.EmitterIndex = 7; platform(scenario,Trajectory=traj,Emitters=emitter5); % % Helicopter taking off ht = 700; spd = 50*1e3/3600; % m/s wp = [-spd*scenarioDuration/2 -13e3 -ht;spd*scenarioDuration/2 -13e3 -ht-300]; traj = waypointTrajectory(wp,[0 scenarioDuration]); emitter6 = clone(emitter); emitter6.EmitterIndex = 8; platform(scenario,Trajectory=traj,Emitters=emitter6); % % Helicopter landing ht = 600; spd = 10*1e3/3600; % m/s wp = [-1e3 -5e3+spd*scenarioDuration/2 -ht;-1e3 -5e3+-spd*scenarioDuration/2 -ht+280]; traj = waypointTrajectory(wp,[0 scenarioDuration]); emitter7 = clone(emitter); emitter7.EmitterIndex = 9; platform(scenario,Trajectory=traj,Emitters=emitter7); end
helperSyncSensor2spec
- A helper to synchronize sensor properties to the sensor specification
function spec = helperSyncSensor2spec(sensor,spec) if isa(sensor,'fusionRadarSensor') && strcmpi(sensor.DetectionMode,'monostatic') spec = syncRadar2spec(sensor,spec); elseif isa(sensor,'fusionRadarSensor') && strcmpi(sensor.DetectionMode,'esm') spec = syncESM2spec(sensor,spec); end end function spec = syncRadar2spec(sensor,spec) spec.MountingLocation = sensor.MountingLocation; spec.MountingAngles = sensor.MountingAngles; spec.HasElevation = sensor.HasElevation; spec.HasRangeRate = sensor.HasRangeRate; spec.FieldOfView = sensor.FieldOfView; spec.RangeLimits = sensor.RangeLimits; spec.RangeRateLimits = sensor.RangeRateLimits; spec.AzimuthResolution = sensor.AzimuthResolution; spec.RangeResolution = sensor.RangeResolution; spec.ElevationResolution = sensor.ElevationResolution; spec.RangeRateResolution = sensor.RangeRateResolution; spec.DetectionProbability = sensor.DetectionProbability; spec.FalseAlarmRate = sensor.FalseAlarmRate; end function spec = syncESM2spec(sensor,spec) spec.MountingLocation = sensor.MountingLocation; spec.MountingAngles = sensor.MountingAngles; spec.HasElevation = sensor.HasElevation; spec.FieldOfView = sensor.FieldOfView; spec.AzimuthResolution = sensor.AzimuthResolution; spec.ElevationResolution = sensor.ElevationResolution; spec.DetectionProbability = sensor.DetectionProbability; spec.FalseAlarmRate = sensor.FalseAlarmRate; end
helperSavePlatformPoses - Record platform poses at tracker update timestamps
function platformPoses = helperSavePlatformPoses(scenario,trackerInterval) endind = ceil(scenario.StopTime/trackerInterval); samplePose = struct(PlatformID=1,ClassID=0,Position=zeros(1,3),Velocity=zeros(1,3),Acceleration=zeros(1,3),Orientation=eye(3),AngularVelocity=zeros(1,3)); platformPoses = repmat(samplePose,numel(scenario.Platforms),endind); sampleTimes = trackerInterval*(1:endind); for p = 1:numel(scenario.Platforms) thisPlatform = scenario.Platforms{p}; trajectory = thisPlatform.Trajectory; if isa(trajectory,"waypointTrajectory") [position,orientation,velocity,acceleration,angularVelocity] = lookupPose(trajectory,sampleTimes); orientation = rotmat(orientation,"frame"); for i = 1:endind platformPoses(p,i) = struct(PlatformID=thisPlatform.PlatformID,ClassID=thisPlatform.ClassID, ... Position=position(i,:),Velocity=velocity(i,:),Acceleration=acceleration(i,:), ... Orientation=orientation(:,:,i),AngularVelocity=angularVelocity(i,:)); end else % Assumes that non waypointTrajectory means stationary object platformPoses(p,:) = repmat(struct(PlatformID=thisPlatform.PlatformID,ClassID=thisPlatform.ClassID, ... Position=thisPlatform.Position,Velocity=zeros(1,3),Acceleration=zeros(1,3),Orientation=thisPlatform.Orientation, ... AngularVelocity=zeros(1,3)),1,endind); end end end
helperRecordSensorData - A helper to record sensor data from the scenario
function sensorData = helperRecordSensorData(scenario) %#ok<DEFNU> if isa(scenario,"trackingScenarioRecording") recording = scenario.RecordedData; % We only need the struct elseif isa(scenario,"struct") recording = scenario; elseif isa(scenario,"trackingScenario") || isa(scenario,"radarScenario") disp("Recording the scenario to get sensor data. This may take several minutes."); recording = record(scenario,"Rotmat",IncludeSensors=true,RecordingFormat="struct",InitialSeed=seed); else error("This example utility function can only be used with a trackingScenario a trackingScenarioRecording"); end if ~isfield(recording,"SensorConfigurations") || ~isfield(recording,"Detections") error("Scenario recording must contain sensor configurations and detections to record sensor data"); end % Assume all steps have the same sensors. sensors = recording(1).SensorConfigurations; detBuffer = vertcat(recording.Detections); for i = 1:numel(recording) for j = 1:numel(recording(i).SensorConfigurations) recording(i).SensorConfigurations(j).Time = recording(i).SimulationTime; end end configBuffer = reshape([recording.SensorConfigurations],[],1); sensorData = cell(1,numel(sensors)); for k = 1:numel(sensors) sensor = sensors(k); sensorIdx = sensor.SensorIndex; sensorMountingAngles = rotmat(recording(1).CoverageConfig(k).Orientation,"frame"); detSensorIndex = cellfun(@(x)x.SensorIndex,detBuffer); dets = detBuffer(detSensorIndex == sensorIdx); cfgSensorIndex = arrayfun(@(x)x.SensorIndex, configBuffer); cfgs = configBuffer(cfgSensorIndex == sensorIdx); if sensor.MeasurementParameters.HasRange sensorData{k} = recordMonostaticData(sensorMountingAngles, dets, cfgs); else sensorData{k} = recordESMData(sensorMountingAngles, dets, cfgs); end end end function activeRadarData = recordMonostaticData(sensorMountingAngles, detections, cfgs) % Reported config information lookMountRot = cell(numel(cfgs),1); for i = 1:numel(cfgs) if cfgs(i).MeasurementParameters(1).IsParentToChild lookMountRot{i} = cfgs(i).MeasurementParameters(1).Orientation; else lookMountRot{i} = cfgs(i).MeasurementParameters(1).Orientation'; end end beamRot = cellfun(@(x)x*sensorMountingAngles',lookMountRot,UniformOutput=false); ypr = cellfun(@(x)fusion.internal.frames.rotmat2ypr(x,"degrees"),beamRot,UniformOutput=false); lookTime = arrayfun(@(x)x.Time,cfgs); lookAzimuth = cellfun(@(x)x(1),ypr); lookElevation = cellfun(@(x)-x(2),ypr); % Reported measurement information detectionTime = cellfun(@(x)x.Time,detections); azimuth = cellfun(@(x)x.Measurement(1),detections); elevation = cellfun(@(x)x.Measurement(2),detections); range = cellfun(@(x)x.Measurement(3),detections); rangeRate = cellfun(@(x)x.Measurement(4),detections); azimuthAccuracy = cellfun(@(x)x.MeasurementNoise(1,1),detections); elevationAccuracy = cellfun(@(x)x.MeasurementNoise(2,2),detections); rangeAccuracy = cellfun(@(x)x.MeasurementNoise(3,3),detections); rangeRateAccuracy = cellfun(@(x)x.MeasurementNoise(4,4),detections); activeRadarData = struct(LookTime=lookTime(:)',LookAzimuth=lookAzimuth(:)', ... LookElevation=lookElevation(:)', DetectionTime=detectionTime(:)', ... Azimuth=azimuth(:)',Elevation=elevation(:)', Range=range(:)', ... RangeRate=rangeRate(:)',AzimuthAccuracy=sqrt(azimuthAccuracy(:)'), ... ElevationAccuracy=sqrt(elevationAccuracy(:)'),RangeAccuracy=sqrt(rangeAccuracy(:)'), ... RangeRateAccuracy=sqrt(rangeRateAccuracy(:)')); end function esmData = recordESMData(sensorMountingAngles, detections, cfgs) % Reported config information lookMountRot = cell(numel(cfgs),1); for i = 1:numel(cfgs) if cfgs(i).MeasurementParameters(1).IsParentToChild lookMountRot{i} = cfgs(i).MeasurementParameters(1).Orientation; else lookMountRot{i} = cfgs(i).MeasurementParameters(1).Orientation'; end end beamRot = cellfun(@(x)x*sensorMountingAngles',lookMountRot,UniformOutput=false); ypr = cellfun(@(x)fusion.internal.frames.rotmat2ypr(x,"degrees"),beamRot,UniformOutput=false); lookAzimuth = cellfun(@(x)x(1),ypr); lookElevation = cellfun(@(x)-x(2),ypr); lookTime = arrayfun(@(x)x.Time,cfgs); detectionTime = cellfun(@(x)x.Time,detections); azimuth = cellfun(@(x)x.Measurement(1),detections); elevation = cellfun(@(x)x.Measurement(2),detections); azimuthAccuracy = cellfun(@(x)x.MeasurementNoise(1,1),detections); elevationAccuracy = cellfun(@(x)x.MeasurementNoise(2,2),detections); esmData = struct(LookTime=lookTime(:)',LookAzimuth=lookAzimuth(:)', ... LookElevation=lookElevation(:)',DetectionTime=detectionTime(:)', ... Azimuth=azimuth(:)',Elevation=elevation(:)',AzimuthAccuracy=sqrt(azimuthAccuracy(:)'), ... ElevationAccuracy=sqrt(elevationAccuracy(:)')); end
helperSplitSensorData - A helper to split recorded sensor data to the tracker update interval
function varargout=helperSplitSensorData(sensorData,stopTime,trackerInterval) varargout = cell(1,nargout); if iscell(sensorData) varargout = cellfun(@(s) splitOneSensorData(s,stopTime,trackerInterval), sensorData, UniformOutput=false); else varargout{1} = splitOneSensorData(sensorData,stopTime,trackerInterval); end end function data=splitOneSensorData(sensorData,stopTime,trackerInterval) frameTime = 0; ind = 1; numFrames = ceil(stopTime/trackerInterval); data = repmat(sensorData,numFrames,1); while frameTime < stopTime withinInterval = sensorData.LookTime >= frameTime & sensorData.LookTime < frameTime+trackerInterval; data(ind).LookTime = sensorData.LookTime(withinInterval); data(ind).LookAzimuth = sensorData.LookAzimuth(withinInterval); data(ind).LookElevation = sensorData.LookElevation(withinInterval); withinInterval = sensorData.DetectionTime >=frameTime & sensorData.DetectionTime < frameTime+trackerInterval; data(ind).DetectionTime = sensorData.DetectionTime(withinInterval); if isfield(sensorData,"Azimuth") data(ind).Azimuth = sensorData.Azimuth(1,withinInterval); data(ind).AzimuthAccuracy = sensorData.AzimuthAccuracy(1,withinInterval); end if isfield(sensorData,"Elevation") data(ind).Elevation = sensorData.Elevation(1,withinInterval); data(ind).ElevationAccuracy = sensorData.ElevationAccuracy(1,withinInterval); end if isfield(sensorData,"Range") data(ind).Range = sensorData.Range(1,withinInterval); data(ind).RangeAccuracy = sensorData.RangeAccuracy(1,withinInterval); end if isfield(sensorData,"RangeRate") data(ind).RangeRate = sensorData.RangeRate(1,withinInterval); data(ind).RangeRateAccuracy = sensorData.RangeRateAccuracy(1,withinInterval); end if isfield(sensorData,"TargetIndex") data(ind).TargetIndex = sensorData.TargetIndex(1,withinInterval); end frameTime = frameTime + trackerInterval; ind = ind + 1; end end