# 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.

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

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

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

### Description

creates a
`PHD`

= ggiwphd`ggiwphd`

filter with default property values.

allows you to specify the `PHD`

= ggiwphd(States,StateCovariances)`States`

and
`StateCovariances`

of the Gaussian distribution for each component in
the density. `States`

and `StateCovariances`

set the
properties of the same names.

also allows you to set properties for the filter using one or more name-value pairs.
Enclose each property name in quotes.`phd`

= ggiwphd(States,StateCovariances,`Name,Value`

)

## 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))

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

**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™.

Usage notes and limitations:

The code generation configuration must allow recursion to use

`merge`

method.

## See Also

`trackingSensorConfiguration`

| `trackerPHD`

| `partitionDetections`

| `gmphd`

**Introduced in R2019a**

## Open Example

You have a modified version of this example. Do you want to open this example with your edits?

## MATLAB Command

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

# Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)