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
= ggiwphdggiwphd
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
istrue
, 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 timek
, anddT
is the time step.If
HasAdditiveProcessNoise
isfalse
, 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 timek
,w(k)
is the process noise at timek
, anddT
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
istrue
, 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 timek
,dT
is the time step, andJx(k)
denotes the Jacobian of the state transition function with respect to the state. The Jacobian is an M-by-M matrix at timek
, where M is the dimension of the state.If
HasAdditiveProcessNoise
isfalse
, 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 timek
. 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
istrue
, 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
isfalse
, specify the process noise covariance as a Q-by-Q matrix. Q is the size of the process noise vector. You must specifyProcessNoise
before any call to thepredict
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:
where k and k+1 represent two consecutive time steps. The mean (E) and variance (Var) of a Gamma distribution are:
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
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
istrue
, 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 timek
andz(k)
is the corresponding measurement. Theparameters
argument stands for all additional arguments required by the measurement function.If
HasAdditiveMeasurementNoise
isfalse
, 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 timek
andv(k)
is the measurement noise at timek
. Theparameters
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
istrue
, 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 timek
.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. Theparameters
argument stands for all arguments required by the measurement function.If
HasAdditiveMeasurmentNoise
isfalse
, 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 timek
andv(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. Theparameters
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
See Also
trackingSensorConfiguration
| trackerPHD
| partitionDetections
| gmphd
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)