ggiwphd

Gamma Gaussian Inverse Wishart (GGIW) PHD filter

Description

The ggiwphd object is a filter that implements the probability hypothesis density (PHD) using a mixture of Gamma Gaussian Inverse-Wishart components.

GGIW implementation of a PHD filter is typically used to track extended objects. An extended object can produce multiple detections per sensor, and the GGIW filter uses the random matrix model to account for the spatial distribution of these detections. The filter consists of three distributions to represent the state of an extended object.

1. Gaussian distribution — represents the kinematic state of the extended object.

2. Gamma distribution — represents the expected number of detections on a sensor from the extended object.

3. Inverse-Wishart (IW) distribution — represents the spatial extent of the target. In 2-D space, the extent is represented by a 2-by-2 random positive definite matrix, which corresponds to a 2-D ellipse description. In 3-D space, the extent is represented by a 3-by-3 random matrix, which corresponds to a 3-D ellipsoid description. The probability density of these random matrices is given as an Inverse-Wishart distribution.

For details about ggiwphd, see [1] and [2].

Note

ggiwphd object is not compatible with trackerGNN, trackerJPDA, and trackerTOMHT system objects.

Creation

Syntax

PHD = ggiwphd

PHD = ggiwphd(States,StateCovariances)

phd = ggiwphd(States,StateCovariances,Name,Value)

Description

PHD = ggiwphd creates a ggiwphd filter with default property values.

PHD = ggiwphd(States,StateCovariances) allows you to specify the States and StateCovariances of the Gaussian distribution for each component in the density. States and StateCovariances set the properties of the same names.

phd = ggiwphd(States,StateCovariances,Name,Value) also allows you to set properties for the filter using one or more name-value pairs. Enclose each property name in quotes.

Properties

States — State of each component in filter

P-by-N matrix

State of each component in the filter, specified as a P-by-N matrix, where P is the dimension of the state and N is the number of components. Each column of the matrix corresponds to the state of each component. The default value for States is a 6-by-2 matrix, in which the elements of the first column are all 0, and the elements of the second column are all 1.

If you want a filter with single-precision floating-point variables, specify States as single-precision vector variables. For example,

filter = ggiwphd(single(zeros(6,4)),single(ones(6,6,4)))

Data Types: single | double

StateCovariances — State estimate error covariance of each component in filter

P-by-P-by-N array

State estimate error covariance of each component in the filter, specified as a P-by-P-by-N array, where P is the dimension of the state and N is the number of components. Each page (P-by-P matrix) of the array corresponds to the covariance matrix of each component. The default value for StateCovariances is a 6-by-6-by-2 array, in which each page (6-by-6 matrix) is an identity matrix.

Data Types: single | double

PositionIndex — Indices of position coordinates in state

[1 3 5] | row vector of positive integers

Indices of position coordinates in the state, specified as a row vector of positive integers. For example, by default the state is arranged as [x;vx;y;vy;z;vz] and the corresponding position index is [1 3 5] representing x-, y- and z-position coordinates.

Example: [1 2 3]

Data Types: single | double

StateTransitionFcn — State transition function

@constvel (default) | function handle

State transition function, specified as a function handle. This function calculates the state vector at time step k from the state vector at time step k–1. The function can also include noise values.

• If HasAdditiveProcessNoise is true, specify the function using one of these syntaxes:

  x(k) = transitionfcn(x(k-1))
  

  x(k) = transitionfcn(x(k-1),dT)

  where x(k) is the state estimate at time k, and dT is the time step.

• If HasAdditiveProcessNoise is false, specify the function using one of these syntaxes:

  x(k) = transitionfcn(x(k-1),w(k-1))
  

  x(k) = transitionfcn(x(k-1),w(k-1),dT)

  where x(k) is the state estimate at time k, w(k) is the process noise at time k, and dT is the time step.

Example: @constacc

Data Types: function_handle

StateTransitionJacobianFcn — Jacobian of state transition function

@constveljac (default) | function handle

The Jacobian of the state transition function, specified as a function handle. This function has the same input arguments as the state transition function.

• If HasAdditiveProcessNoise is true, specify the Jacobian function using one of these syntaxes:

  Jx(k) = statejacobianfcn(x(k))
  

  Jx(k) = statejacobianfcn(x(k),dT)

  where x(k) is the state at time k, dT is the time step, and Jx(k) denotes the Jacobian of the state transition function with respect to the state. The Jacobian is an M-by-M matrix at time k, where M is the dimension of the state.

• If HasAdditiveProcessNoise is false, specify the Jacobian function using one of these syntaxes:

  [Jx(k),Jw(k)] = statejacobianfcn(x(k),w(k))
  

  [Jx(k),Jw(k)] = statejacobianfcn(x(k),w(k),dT)

  where w(k) is a Q-element vector of the process noise at time k. Q is the dimension of the process noise. Unlike the case of additive process noise, the process noise vector in the nonadditive noise case need not have the same dimensions as the state vector.

  Jw(k) denotes the M-by-Q Jacobian of the predicted state with respect to the process noise elements, where M is the dimension of the state.

If not specified, the Jacobians are computed by numerical differencing at each call of the predict function. This computation can increase the processing time and numerical inaccuracy.

Example: @constaccjac

Data Types: function_handle

ProcessNoise — Process noise covariance

eye(3) (default) | positive real-valued scalar | positive-definite real-valued matrix

Process noise covariance:

• When HasAdditiveProcessNoise is true, specify the process noise covariance as a scalar or a positive definite real-valued M-by-M matrix. M is the dimension of the state vector. When specified as a scalar, the matrix is a multiple of the M-by-M identity matrix.

• When HasAdditiveProcessNoise is false, specify the process noise covariance as a Q-by-Q matrix. Q is the size of the process noise vector. You must specify ProcessNoise before any call to the predict object function.

Example: [1.0 0.05; 0.05 2]

HasAdditiveProcessNoise — Model additive process noise

false (default)

Option to model processes noise as additive, specified as true or false. When this property is true, process noise is added to the state vector. Otherwise, noise is incorporated into the state transition function.

Example: true

Shapes — Shape parameter of Gamma distribution for each component

[1 1] (default) | 1-by-N row vector of positive real values

Shape parameter of Gamma distribution for each component, specified as a 1-by-N row vector of positive real values. N is the number of components in the density.

Example: [1.0 0.95 2]

Data Types: single | double

Rates — Rate parameter of Gamma distribution for each component

[1 1] (default) | 1-by-N row vector of positive real value

Rate parameter of Gamma distribution for each component, specified as a 1-by-N row vector of positive real values. N is the number of components in the density.

Example: [1.2 0.85 1.5]

Data Types: single | double

GammaForgettingFactors — Forgetting factor of Gamma distribution for each component

[1 1] (default) | 1-by-N row vector of positive real value

Forgetting factor of Gamma distribution for each component, specified as a 1-by-N row vector of positive real values. N is the number of components in the density. During prediction, for each component, the Gamma distribution parameters, shape (α) and rate (β), are both divided by forgetting factor n:

$$\begin{array}{l}{a}_{k+1\left|\right|k}=\frac{{\alpha }_k}{n_k}\\ {\beta }_{k+1|k}=\frac{{\beta }_k}{n_k}\end{array}$$

where k and k+1 represent two consecutive time steps. The mean (E) and variance (Var) of a Gamma distribution are:

$$\begin{array}{l}E=\frac{\alpha }{\beta }\\ Var=\frac{\alpha }{{\beta }^2}\end{array}$$

Therefore, the division action will keep the expected measurement rate as a constant, but increase the variance of the Gamma distribution exponentially with time if the forgetting factor n is larger than 1.

Example: [1.2 1.1 1.4]

Data Types: single | double

DegreesOfFreedom — Degrees of freedom parameter of Inverse-Wishart distribution for each component

[100 100] (default) | 1-by-N row vector of positive real value

Degrees of freedom parameter of Inverse-Wishart distribution for each component, specified as a 1-by-N row vector of positive real values. N is the number of components in the density.

Example: [55.2 31.1 20.4]

Data Types: single | double

ScaleMatrices — Scale matrix of Inverse-Wishart distribution for each component

d-by-d-by-N array of positive real value

Scale matrix of Inverse-Wishart distribution for each component, specified as a d-by-d-by-N array of positive real values. d is the dimension of the space (for example, d = 2 for 2-D space), and N is the number of components in the density. The default value for ScaleMatrices is a 3-by-3-by-2 array, where each page (3-by-3 matrix) of the array is 100*eye(3).

Example: 20*eye(3,3,4)

Data Types: single | double

ExtentRotationFcn — Rotation transition function of target's extent

@(x,varargin)eye(3) (default) | function handle

Rotation transition function of target's extent, specified as a function handle. The function allows predicting the rotation of the target's extent when the object's angular velocity is estimated in the state vector. To define your own extent rotation function, follow the syntax given by

R = myRotationFcn(x,dT)

where x is the component state, dT is the time step, and R is the corresponding rotation matrix. Note that R is returned as a 2-by-2 matrix if the extent is 2-D, and a 3-by-3 matrix if the extent is 3-D. The extent at the next step is given by

$$Ex(t+dT)=R\times Ex(t)\times R^T$$

where Ex(t) is the extent at time t.

Example: @myRotationFcn

Data Types: function_handle

TemporalDecay — Temporal decay factor of IW distribution

100 (default) | positive scalar

Temporal decay factor of IW distribution, specified as a positive scalar. You can use this property to control the extent uncertainty (variance of IW distribution) during prediction. The smaller the TemporalDecay value is, the faster the variance of IW distribution increases.

Example: 120

Data Types: single | double

Labels — Label of each component in mixture

[0 0] (default) | 1-by-N row vector of nonnegative integer

Label of each component in the mixture, specified as a 1-by-N row vector of nonnegative integers. N is the number of components in the density. Each component can only have one label, but multiple components can share the same label.

Example: [1 2 3]

Data Types: single | double

Weights — Weight of each component in mixture

[1 1] (default) | 1-by-N row vector of positive real value

Weight of each component in the density, specified as a 1-by-N row vector of positive real values. N is the number of components in the density. The weights are given in the sequence as shown in the labels property.

Example: [1.1 0.82 1.1]

Data Types: single | double

Detections — Detections

K-element cell array of objectDetection objects

Detections, specified as a K-element cell array of objectDetection objects, where K is the number of detections. You can create detections directly, or you can obtain detections from the outputs of sensor objects, such as radarSensor, monostaticRadarSensor, irSensor, and sonarSensor.

Data Types: single | double

MeasurementFcn — Measurement model function

@cvmeas (default) | function handle

Measurement model function, specified as a function handle. This function specifies the transition from state to measurement. Input to the function is the P-element state vector. The output is the M-element measurement vector. The function can take additional input arguments, such as sensor position and orientation.

• If HasAdditiveMeasurementNoise is true, specify the function using one of these syntaxes:

  z(k) = measurementfcn(x(k))
  

  z(k) = measurementfcn(x(k),parameters)

  where x(k) is the state at time k and z(k) is the corresponding measurement. The parameters argument stands for all additional arguments required by the measurement function.

• If HasAdditiveMeasurementNoise is false, specify the function using one of these syntaxes:

  z(k) = measurementfcn(x(k),v(k))
  

  z(k) = measurementfcn(x(k),v(k),parameters)

  where x(k) is the state at time k and v(k) is the measurement noise at time k. The parameters argument stands for all additional arguments required by the measurement function.

Example: @cameas

Data Types: function_handle

MeasurementJacobianFcn — Jacobian of measurement function

@cvmeasjac (default) | function handle

Jacobian of the measurement function, specified as a function handle. The function has the same input arguments as the measurement function. The function can take additional input parameters, such as sensor position and orientation.

• If HasAdditiveMeasurmentNoise is true, specify the Jacobian function using one of these syntaxes:

  Jmx(k) = measjacobianfcn(x(k))
  

  Jmx(k) = measjacobianfcn(x(k),parameters)

  where x(k) is the state at time k. Jmx(k) denotes the M-by-P Jacobian of the measurement function with respect to the state. M is the dimension of the measurement, and P is the dimension of the state. The parameters argument stands for all arguments required by the measurement function.

• If HasAdditiveMeasurmentNoise is false, specify the Jacobian function using one of these syntaxes:

  [Jmx(k),Jmv(k)] = measjacobianfcn(x(k),v(k))
  

  [Jmx(k),Jmv(k)] = measjacobianfcn(x(k),v(k),parameters)

  where x(k) is the state at time k and v(k) is an R-dimensional sample noise vector. Jmx(k) denotes the M-by-P Jacobian matrix of the measurement function with respect to the state. Jmv(k) denotes the Jacobian of the M-by-R measurement function with respect to the measurement noise. The parameters argument stands for all arguments required by the measurement function.

If not specified, measurement Jacobians are computed using numerical differencing at each call to the correct function. This computation can increase processing time and numerical inaccuracy.

Example: @cameasjac

Data Types: function_handle

HasAdditiveMeasurementNoise — Model additive measurement noise

false (default)

Option to model measurement noise as additive, specified as true or false. When this property is true, measurement noise is added to the state vector. Otherwise, noise is incorporated into the measurement function.

Example: true

MaxNumDetections — Maximum number of detections

100 (default) | positive integer

Maximum number of detections the ggiwphd filter can take as input, specified as a positive integer.

Example: 50

Data Types: single | double

MaxNumComponents — Maximum number of components

1000 (default) | positive integer

Maximum number of components the ggiwphd filter can maintain, specified as a positive integer.

Note

When you use a ggiwphd object in a trackerPHD object by specifying the SensorConfigurations property of the tracker, the tracker determines the maximum number of components based on its own MaxNumComponents property instead of the MaxNumComponents property of the ggiwphd object.

Data Types: single | double

Object Functions

append	Append two phd filter objects
correct	Correct phd filter with detections
correctUndetected	Correct phd filter with no detection hypothesis
extractState	Extract target state estimates from the phd filter
labeledDensity	Keep components with a given label ID
likelihood	Log-likelihood of association between detection cells and components in the density
merge	Merge components in the density of phd filter
predict	Predict probability hypothesis density of phd filter
prune	Prune the filter by removing selected components
scale	Scale weights of components in the density
clone	Create duplicate phd filter object

Examples

Create ggiwphd Filter with Two 3-D Components

Creating a ggiwphd filter with two 3-D constant velocity components. The initial states of the two components are [0;0;0;0;0;0] and [1;0;1;0;1;0], respectively. Both these components have position covariance equal to 1 and velocity covariance equal to 100. By default, ggiwphd creates a 3-D extent matrix for each component.

states = [zeros(6,1),[1;0;1;0;1;0]];
cov1 = diag([1 100 1 100 1 100]);
covariances = cat(3,cov1,cov1);

phd = ggiwphd(states,covariances,'StateTransitionFcn',@constvel,...
             'StateTransitionJacobianFcn',@constveljac,...
             'MeasurementFcn',@cvmeas,'MeasurementJacobianFcn',@cvmeasjac,...
             'ProcessNoise',eye(3),'HasAdditiveProcessNoise',false,...
             'PositionIndex',[1;3;5]);

Specify information about extent.

dofs = [21 30];
scaleMatrix1 = 13*diag([4.7 1.8 1.4].^2);
scaleMatrix2 = 22*diag([1.8 4.7 1.4].^2); 
scaleMatrices = cat(3,scaleMatrix1,scaleMatrix2);
phd.DegreesOfFreedom = dofs;
phd.ScaleMatrices = scaleMatrices;
phd.ExtentRotationFcn = @(x,dT)eye(3); % No rotation during prediction

Predict the filter 0.1 second ahead.

predict(phd,0.1);

Specify detections at 0.1 second. The filter receives 10 detections at the current scan.

detections = cell(10,1);
rng(2018); % Reproducible results
for i = 1:10
   detections{i} = objectDetection(0.1,randi([0 1]) + randn(3,1));
end
phd.Detections = detections;

Select two detection cells and calculate their likelihoods.

detectionIDs = false(10,2);
detectionIDs([1 3 5 7 9],1) = true;
detectionIDs([2 4 6 8 10],2) = true;
lhood = likelihood(phd,detectionIDs)

lhood = 2×2

    1.5575   -0.3183
    0.1513   -0.7616

Correct the filter with the two detection cells and associated likelihoods.

correct(phd,detectionIDs, exp(lhood)./sum(exp(lhood),1));
phd

phd = 
  ggiwphd with properties:

                         States: [6x4 double]
               StateCovariances: [6x6x4 double]
                  PositionIndex: [3x1 double]
             StateTransitionFcn: @constvel
     StateTransitionJacobianFcn: @constveljac
                   ProcessNoise: [3x3 double]
        HasAdditiveProcessNoise: 0

                         Shapes: [6 6 6 6]
                          Rates: [2 2 2 2]
         GammaForgettingFactors: [1 1 1 1]

               DegreesOfFreedom: [25.9870 34.9780 25.9870 34.9780]
                  ScaleMatrices: [3x3x4 double]
              ExtentRotationFcn: @(x,dT)eye(3)
                  TemporalDecay: 100

                        Weights: [0.8032 0.1968 0.6090 0.3910]
                         Labels: [0 0 0 0]

                     Detections: {1x10 cell}
                 MeasurementFcn: @cvmeas
         MeasurementJacobianFcn: @cvmeasjac
    HasAdditiveMeasurementNoise: 1

Merge components in the filter.

merge(phd,5);
phd

phd = 
  ggiwphd with properties:

                         States: [6x2 double]
               StateCovariances: [6x6x2 double]
                  PositionIndex: [3x1 double]
             StateTransitionFcn: @constvel
     StateTransitionJacobianFcn: @constveljac
                   ProcessNoise: [3x3 double]
        HasAdditiveProcessNoise: 0

                         Shapes: [6 6.0000]
                          Rates: [2 2]
         GammaForgettingFactors: [1 1]

               DegreesOfFreedom: [25.9870 34.9780]
                  ScaleMatrices: [3x3x2 double]
              ExtentRotationFcn: @(x,dT)eye(3)
                  TemporalDecay: 100

                        Weights: [1.4122 0.5878]
                         Labels: [0 0]

                     Detections: {1x10 cell}
                 MeasurementFcn: @cvmeas
         MeasurementJacobianFcn: @cvmeasjac
    HasAdditiveMeasurementNoise: 1

Extract state estimates and detections.

targetStates = extractState(phd,0.5);
tStates = targetStates.State

tStates = 6×1

    0.1947
    0.9733
    0.8319
    4.1599
   -0.0124
   -0.0621

    
d = [detections{:}];
measurements = [d.Measurement];

Visualize the results.

figure()
plot3(measurements(1,:),measurements(2,:),measurements(3,:),'x','MarkerSize',10,'MarkerEdgeColor','b');
hold on;
plot3( tStates(1,:),tStates(3,:),tStates(5,:),'ro');
xlabel('x')
ylabel('y')
zlabel('z')
legend('Detections','Components')

References

[1] Granstorm, K., and O. Orguner." A PHD filter for tracking multiple extended targets using random matrices." IEEE Transactions on Signal Processing. Vol. 60, Number 11, 2012, pp. 5657-5671.

[2] Granstorm, K., and A. Natale, P. Braca, G. Ludeno, and F. Serafino."Gamma Gaussian inverse Wishart probability hypothesis density for extended target tracking using X-band marine radar data." IEEE Transactions on Geoscience and Remote Sensing. Vol. 53, Number 12, 2015, pp. 6617-6631.

Extended Capabilities

C/C++ Code Generation

Generate C and C++ code using MATLAB® Coder™.

• The filter object supports non-dynamic memory allocation code generation. For details on non-dynamic memory allocation code generation, see Generate Code with Strict Single-Precision and Non-Dynamic Memory Allocation.

• The filter object supports single-precision code generation with these limitations:

  • The motion model and measurement model used in the filter must support single-precision.

  • In the filter, some intermediate variables can possibly use double-precision.

  For details, see Generate Code with Strict Single-Precision and Non-Dynamic Memory Allocation.

Version History

Introduced in R2019a