Gaussian mixture (GM) PHD filter
The gmphd
object is a filter that implements the
probability hypothesis density (PHD) using a mixture of Gaussian components. The filter
assumes the target states are Gaussian and represents these states using a mixture of Gaussian
components. You can use a gmphd
filter to track extended objects or
point targets. In tracking, a point object returns at most one detention per sensor scan, and
an extended object can return multiple detections per sensor scan.
You can directly create a gmphd
filter. You can also initialize a
gmphd
filter used with trackerPHD
by
specifying the FilterInitializationFcn
property of trackingSensorConfiguration
. You can use the provided initcvgmphd
, initctgmphd
, initcagmphd
, and initctrectgmphd
as
initialization functions. Or, you can create your own initialization functions.
creates a phd
= gmphdgmphd
filter with default property values.
allows
you to specify the states and corresponding state covariances of the Gaussian distribution
for each component in the density. phd
= gmphd(states,stateCovariances)states
and
stateCovariances
set the States
and StateCovariances
properties of the filter.
also allows you to specify properties for the filter using one or more namevalue pairs.
Enclose each property name in quotes.phd
= gmphd(states,stateCovariances,Name,Value
)
States
— State of each component in filterState of each component in the filter, specified as a
PbyN 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 one component. The default value for
States
is a 6by2 matrix, in which the elements of the first
column are all 0, and the elements of the second column are all 1.
Data Types: single
 double
StateCovariances
— State estimate error covariance of each component in filterState estimate error covariance of each component in the filter, specified as a
PbyPbyN array, where
P is the dimension of the state and N is the
number of components. Each page (PbyP matrix) of
the array corresponds to the covariance matrix of each component. The default value for
StateCovariances
is a 6by6by2 array, in which each page
(6by6 matrix) is an identity matrix.
Data Types: single
 double
StateTransitionFcn
— State transition function@constvel
(default)  function handleState 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(k1))
x(k) = transitionfcn(x(k1),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(k1),w(k1))
x(k) = transitionfcn(x(k1),w(k1),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 handleJacobian 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)
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 a PbyP matrix at time
k
, where P 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)
w(k)
is a Qelement vector of the process
noise at time k
. Unlike the case of additive process noise,
the process noise vector in the nonadditive noise case doesn't need to have the
same dimensions as the state vector.Jw(k)
denotes the
PbyQ Jacobian of the predicted state
with respect to the process noise elements, where P 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 covarianceeye(3)
(default)  positive realvalued scalar  positive definite realvalued matrixProcess noise covariance:
When HasAdditiveProcessNoise
is true
,
specify the process noise covariance as a realvalued scalar or a positive
definite PbyP matrix.
P is the dimension of the state vector. When specified as a
scalar, the matrix is a multiple of the
PbyP identity matrix.
When HasAdditiveProcessNoise
is false
,
specify the process noise covariance as a
QbyQ 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 noisefalse
(default)  true
Option to model process 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
HasExtent
— Indicate if components have extentfalse
(default)  true
Indicate if components have extent, specified as true
or
false
. Set this property to true
if the filter
is intended to track extended objects. An extended object can generate more than one
measurement per sensor scan. Set this property to false
if the filter
is only intended to track point targets.
Example: true
MeasurementOrigin
— Origination of measurements from extended objects'center'
(default)  'extent'
Origination of measurements from extended objects, specified as:
'center'
— The filter assumes the measurements originate
from the mean state of a target. This approach is applicable when the state does
not model the extent of the target even though the target may generate more than
one measurement.
'extent'
— The filter assumes measurements are not centered
at the mean state of a target. For computational efficiency, the expected
measurement is often calculated as a function of the reported measurements
specified by the measurement model function.
Note that the function setups of MeasurementFcn
and
MeasurementJacobianFcn
are different for
'center'
and 'extent'
options. See the
descriptions of MeasurementFcn
and
MeasurementJacobianFcn
for more details.
To enable this property, set the HasExtent
property to
'true'
.
Data Types: double
Labels
— Label of each component in mixture[0 0]
(default)  1byN row vector of nonnegative integerLabel of each component in the mixture, specified as a 1byN row vector of nonnegative integers. N is the number of components in the mixture. 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)  1byN row vector of positive real valueWeight of each component in the mixture, specified as a 1byN
row vector of positive real values. N is the number of components in
the mixture. The weight of each component is given in the same order as the
Labels
property.
Example: [1.1 0.82 1.1]
Data Types: single
 double
Detections
— DetectionsobjectDetection
objectsDetections, specified as a Delement cell array of objectDetection
objects. 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 handleMeasurement model function, specified as a function handle. This function specifies
the transition from state to measurement. Depending on the
HasExtent
and MeasurementOrigin
properties,
the measurement model function needs to be specified differently:
HasExtent
is false
, or
HasExtent
is true
and
MeasurementOrigin
is 'center'
. In these
two cases,
If HasAdditiveMeasurementNoise
is
true
, specify the function using one of these
syntaxes:
z = measurementfcn(x)
z = measurementfcn(x,parameters)
x
is
the estimated Gaussian states at time k
and
x(:,i)
represents the i
th state
component in the mixture. The MbyN
matrix z
is the corresponding measurement, and
z(:,i)
represents the measurement resulting from the
i
th component. Parameters
are
MeasurementParameters
provided in the
objectDetections
set in the Detections
property.If HasAdditiveMeasurementNoise
is
false
, specify the function using one of these
syntaxes:
z = measurementfcn(x,v)
z = measurementfcn(x,v,parameters)
v
is an Rdimensional measurement noise
vector.
HasExtent
is true
and
MeasurementOrigin
is 'extent'
. In this
case, the expected measurements originate from the extent of the target and rely on
the actual distribution of the detections:
If HasAdditiveMeasurementNoise
is
true
, specify the function
using:
z = measurementfcn(x,detections)
x
is the estimated Gaussian states at time k
and
x(:,i)
represents the i
th state
component in the mixture. detections
is a cell array of
objectDetection
objects, and z
is the
expected measurement. Note that z(:,i,j)
must return the
expected measurement based on the i
th state component and
the j
th objectDetection
in
detections
.If HasAdditiveMeasurementNoise
is
false
, specify the function
using:
z = measurementfcn(x,v,detections)
v
is an Rdimensional measurement
noise vector.
HasExtent  MeasurementOrigin  Measurement Function  Note  
false  NA 

 
true  'center'  
true  'extent' 


Data Types: function_handle
MeasurementJacobianFcn
— Jacobian of measurement function@cvmeasjac
(default)  function handleJacobian of the measurement function, specified as a function handle. Depending on
the HasExtent
and MeasurementOrigin
properties, the measurement Jacobian function needs to be specified differently:
HasExtent
is false
, or
HasExtent
is true
and
MeasurementOrigin
is 'center'
. In these
two cases:
If HasAdditiveMeasurmentNoise
is
true
, specify the Jacobian function using one of these
syntaxes:
Jmx = measjacobianfcn(x)
Jmx = measjacobianfcn(x,parameters)
x
is one state
component at time k
and Jmx
is the
MbyP Jacobian of the measurement
function with respect to the state. M is the dimension of the
measurement. Parameters
are
MeasurementParameters
provided in the
objectDetections
set in the Detections
property.If HasAdditiveMeasurmentNoise
is
false
, specify the Jacobian function using one of these
syntaxes:
[Jmx,Jmv] = measjacobianfcn(x,v)
[Jmx,Jmv] = measjacobianfcn(x,v,parameters)
v
is an Rdimensional measurement noise
vector, and Jmv
is the
MbyR Jacobian of the measurement
function with respect to the measurement noise.
HasExtent
is true
and
MeasurementOrigin
is 'extent'
. In this
case, the expected measurements originate from the extent of the target and rely on
the actual distribution of the detections. The measurement Jacobian function must
support one of these two syntaxes:
If HasAdditiveMeasurmentNoise
is
true
, specify the Jacobian function
using:
Jmx = measjacobianfcn(x,detections)
x
is one state estimate component at time
k
. detections
is a set of detections
defined as a cell array of objectDetection
objects.
Jmx
denotes the
MbyPbyD
Jacobian of the measurement function with respect to the state.
M is the dimension of the measurement,
P is the dimension of the state, and D
is the number of objectDetection
objects in
detections
.If HasAdditiveMeasurmentNoise
is
false
, specify the Jacobian function
using:
[Jmx,Jmv] = measjacobianfcn(x,v,detections)
v
is an Rdimensional measurement
noise vector, and Jmv
is the
MbyRbyD
Jacobian of the measurement function with respect to the measurement noise.
Note that Jmx(:,:,j)
must define the state
Jacobian corresponding to the j
th
objectDetection
in detections
.
Jmv(:,:,j)
defines the measurement noise Jacobian corresponding
to the j
th objectDetection
in
detections
.
HasExtent  MeasurementOrigin  Measurement Jacobian Function  Note  
false  NA 
 x is only one Gaussian component in the
mixture.  
true  'center'  
true  'extent' 


Data Types: function_handle
HasAdditiveMeasurementNoise
— Model additive measurement noisefalse
(default)  true
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 detections1000
(default)  positive integerMaximum number of detections the gmphd
filter can take as
input, specified as a positive integer.
Example: 50
Data Types: single
 double
MaxNumComponents
— Maximum number of components1000
(default)  positive integerMaximum number of components the gmphd
filter can maintain,
specified as a positive integer.
Data Types: single
 double
predict  Predict probability hypothesis density of phd filter 
correctUndetected  Correct phd filter with no detection hypothesis 
correct  Correct phd filter with detections 
likelihood  Loglikelihood of association between detection cells and components in the density 
append  Append two phd filter objects 
merge  Merge components in the density of phd filter 
scale  Scale weights of components in the density 
prune  Prune the filter by removing selected components 
labeledDensity  Keep components with a given label ID 
extractState  Extract target state estimates from the phd filter 
clone  Create duplicate phd filter object 
Create a filter with two 3D constant velocity components. The initial state of one component is [0;0;0;0;0;0]. The initial state of the other component is [1;0;1;0;1;0]. Each component is initialized with position covariance equal to 1 and velocity covariance equal to 100.
states = [zeros(6,1) [1;0;1;0;1;0]]; cov1 = diag([1 100 1 100 1 100]); covariances = cat(3,cov1,cov1); phd = gmphd(states, covariances, 'StateTransitionFcn', @constvel,... 'StateTransitionJacobianFcn',@constveljac,... 'MeasurementFcn',@cvmeas,... 'MeasurementJacobianFcn',@cvmeasjac,... 'ProcessNoise', eye(3),... 'HasAdditiveProcessNoise',false);
Predict the filter 0.1 time step ahead.
predict(phd,0.1);
Define three detections using ojbectDetection
.
rng(2019); detections = cell(3,1); detections{1} = objectDetection(0,[1;1;1] + randn(3,1)); detections{2} = objectDetection(0,[0;0;0] + randn(3,1)); detections{3} = objectDetection(0,[4;5;5] + randn(3,1)); phd.Detections = detections;
Calculate the likelihood of each detection. For a pointtarget filter, the partition of detections is unnecessary, and each detection occupies a cell. Therefore, detectionIndices
is an identity matrix. The resulting likelihood of detection 1 and 2 is higher than that of detection 3 because they are closer to the components.
detectionIndices = logical(eye(3)); logLikelihood = likelihood(phd,detectionIndices)
logLikelihood = 2×3
5.2485 4.7774 22.8899
4.5171 5.0008 17.3973
Correct the filter with the scaled likelihood.
lhood = exp(logLikelihood); lhood = lhood./sum(lhood,2); correct(phd,detectionIndices,lhood);
Merge the components with a merging threshold equal to 1.
merge(phd,1);
Extract state estimates with an extract threshold equal to 0.5.
minWeight = 0.5; targetStates = extractState(phd,minWeight); [ts1,ts2]= targetStates.State;
Visualize the results.
% Extract the measurements. d = [detections{:}]; measurements = [d.Measurement]; % Plot the measurements and estimates. figure() plot3(measurements(1,:),measurements(2,:),measurements(3,:),'x','MarkerSize',10,'MarkerEdgeColor','b'); hold on; plot3(ts1(1),ts1(3),ts1(5),'ro'); hold on; plot3(ts2(1),ts2(3),ts2(5),'ro'); xlabel('x') ylabel('y') zlabel('z') hold on; legend('Detections','Components')
[1] Vo, B. N., and W. K. Ma. "The Gaussian Mixture Probability Hypothesis Density Filter."IEEE Transactions on Signal Processing, Vol. 54, No. 11, pp. 4091–4104, 2006.
[2] Granstrom, K., C. Lundquist, and O. Orguner. "Extended Target Tracking Using a Gaussianmixture PHD filter." IEEE Transactions on Aerospace and Electronic Systems, Vol. 48, No. 4, pp. 3268–3286, 2012.
Usage notes and limitations:
The code generation configuration must allow recursion to use the merge
method.
ggiwphd
 initcagmphd
 initctgmphd
 initctrectgmphd
 initcvgmphd
 partitionDetections
 trackerPHD
 trackingSensorConfiguration
A modified version of this example exists on your system. Do you want to open this version instead?
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.
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: .
Select web siteYou can also select a web site from the following list:
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.