gmphd
Gaussian mixture (GM) PHD filter
Description
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.
Creation
Syntax
Description
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
)
Properties
States
— State of each component in filter
PbyN matrix
State 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.
If you want a filter with singleprecision floatingpoint variables, specify
States
as a singleprecision vector variable. For example,
filter = gmphd(single(zeros(6,4)),single(ones(6,6,4)))
Data Types: single
 double
StateCovariances
— State estimate error covariance of each component in filter
PbyPbyN array
State 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 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(k1))
wherex(k) = transitionfcn(x(k1),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(k1),w(k1))
wherex(k) = transitionfcn(x(k1),w(k1),dT)
x(k)
is the state estimate at timek
,w(k)
is the process noise at timek
, anddT
is the time step.For more details on the state transition model, see Algorithms for
trackingEKF
.
Example: @constacc
Data Types: function_handle
StateTransitionJacobianFcn
— Jacobian of state transition function
@constveljac
(default)  function handle
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 a PbyP matrix at timek
, where P 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 Qelement vector of the process noise at timek
. 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 covariance
eye(3)
(default)  positive realvalued scalar  positive definite realvalued matrix
Process noise covariance:
When
HasAdditiveProcessNoise
istrue
, 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
isfalse
, specify the process noise covariance as a QbyQ 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)  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 extent
false
(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.
Dependencies
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 integer
Label 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 value
Weight 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
— Detections
Delement cell array of objectDetection
objects
Detections, 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 handle
Measurement 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
isfalse
, orHasExtent
istrue
andMeasurementOrigin
is'center'
. In these two cases,If
HasAdditiveMeasurementNoise
istrue
, specify the function using one of these syntaxes:z = measurementfcn(x)
where the PbyN matrixz = measurementfcn(x,parameters)
x
is the estimated Gaussian states at timek
andx(:,i)
represents thei
th state component in the mixture. The MbyN matrixz
is the corresponding measurement, andz(:,i)
represents the measurement resulting from thei
th component.Parameters
areMeasurementParameters
provided in theobjectDetections
set in theDetections
property.If
HasAdditiveMeasurementNoise
isfalse
, specify the function using one of these syntaxes:z = measurementfcn(x,v)
wherez = measurementfcn(x,v,parameters)
v
is an Rdimensional measurement noise vector.
HasExtent
istrue
andMeasurementOrigin
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
istrue
, specify the function using:where the PbyN matrixz = measurementfcn(x,detections)
x
is the estimated Gaussian states at timek
andx(:,i)
represents thei
th state component in the mixture.detections
is a cell array ofobjectDetection
objects, andz
is the expected measurement. Note thatz(:,i,j)
must return the expected measurement based on thei
th state component and thej
thobjectDetection
indetections
.If
HasAdditiveMeasurementNoise
isfalse
, specify the function using:wherez = 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 handle
Jacobian 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
isfalse
, orHasExtent
istrue
andMeasurementOrigin
is'center'
. In these two cases:If
HasAdditiveMeasurmentNoise
istrue
, specify the Jacobian function using one of these syntaxes:Jmx = measjacobianfcn(x)
where the Pelement vectorJmx = measjacobianfcn(x,parameters)
x
is one state component at timek
andJmx
is the MbyP Jacobian of the measurement function with respect to the state. M is the dimension of the measurement.Parameters
areMeasurementParameters
provided in theobjectDetections
set in theDetections
property.If
HasAdditiveMeasurmentNoise
isfalse
, specify the Jacobian function using one of these syntaxes:[Jmx,Jmv] = measjacobianfcn(x,v)
where[Jmx,Jmv] = measjacobianfcn(x,v,parameters)
v
is an Rdimensional measurement noise vector, andJmv
is the MbyR Jacobian of the measurement function with respect to the measurement noise.
HasExtent
istrue
andMeasurementOrigin
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
istrue
, specify the Jacobian function using:whereJmx = measjacobianfcn(x,detections)
x
is one state estimate component at timek
.detections
is a set of detections defined as a cell array ofobjectDetection
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 ofobjectDetection
objects indetections
.If
HasAdditiveMeasurmentNoise
isfalse
, specify the Jacobian function using:where[Jmx,Jmv] = measjacobianfcn(x,v,detections)
v
is an Rdimensional measurement noise vector, andJmv
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 thej
thobjectDetection
indetections
.Jmv(:,:,j)
defines the measurement noise Jacobian corresponding to thej
thobjectDetection
indetections
.
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 noise
false
(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 detections
1000
(default)  positive integer
Maximum 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 components
1000
(default)  positive integer
Maximum number of components the gmphd
filter can maintain,
specified as a positive integer.
Data Types: single
 double
Object Functions
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 
Examples
Run gmphd Filter for Point Objects
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')
References
[1] Vo, B. T, 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, Karl, Christian Lundquist, and Omut Orguner. "Extended target tracking using a Gaussianmixture PHD filter." IEEE Transactions on Aerospace and Electronic Systems 48, no. 4 (2012): 32683286.
[3] Panta, Kusha, et al. “Data Association and Track Management for the Gaussian Mixture Probability Hypothesis Density Filter.” IEEE Transactions on Aerospace and Electronic Systems, vol. 45, no. 3, July 2009, pp. 1003–16.
Extended Capabilities
C/C++ Code Generation
Generate C and C++ code using MATLAB® Coder™.
The filter object supports nondynamic memory allocation code generation. For details on nondynamic memory allocation code generation, see Generate Code with Strict SinglePrecision and NonDynamic Memory Allocation.
The filter object supports strict singleprecision code generation with this restriction:
The motion model and measurement model used in the filter must support singleprecision.
For details, see Generate Code with Strict SinglePrecision and NonDynamic Memory Allocation.
Version History
Introduced in R2019b
See Also
trackerPHD
 trackingSensorConfiguration
 partitionDetections
 ggiwphd
 initctrectgmphd
 initctgmphd
 initcvgmphd
 initcagmphd
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)