View License

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video

Highlights from
Gaussian Mixture Probability Hypothesis Density Filter (GM-PHD)

4.2 | 5 ratings Rate this file 39 Downloads (last 30 days) File Size: 94.9 KB File ID: #42769 Version: 1.15
image thumbnail

Gaussian Mixture Probability Hypothesis Density Filter (GM-PHD)



22 Jul 2013 (Updated )

Implementation of the Gaussian mixture probability hypothesis density filter GM-PHD.

| Watch this File

File Information

This is an implementation of the Gaussian mixture probability hypothesis density filter (GM-PHD) described in:
B.-N. Vo, W.-K. Ma, "The Gaussian Mixture Probability Hypothesis Density Filter", IEEE Transactions on Signal Processing, Vol 54, No. 11, November 2006, pp4091-4104.

The submission includes a Matlab implementation of the linear Kalman filter GM-PHD filter and the extended Kalman filter GM-PHD filter algorithms described by Vo & Ma, as well as one of the simulated problems described in their paper. A few modifications were made from Vo & Ma's algorithm but they are for technical reasons and don't change the overall structure of the filter.

The GM-PHD filter is a means of estimating the number and positions of targets in measurement data. Its advantages include a representation of target position uncertainty (using a covariance matrix) as well as target existence uncertainty (using a weight) and the absence of data association in the update step. This implementation is quite heavily commented and will probably be helpful to people trying to learn about GM-PHD filtering, but the paper by Vo & Ma is essential to understand what is really going on.

The simulator creates noisy measurements of two moving targets in a cluttered environment, with a third target being created approximately halfway through the simulation. It's a fairly simple problem but it is effective in demonstrating filter performance. Vo & Ma use it to demonstrate the linear GM-PHD filter but I use it for the EKF version as well, rather than coding up the EKF simulator that they describe.

Ba-Ngu Vo has kindly allowed me to include his implementation of the Optimal Subpattern Assignment (OSPA) metric proposed by Schuhmacher et al in:
Schuhmacher, D.; Ba-Tuong Vo; Ba-Ngu Vo, "A Consistent Metric for Performance Evaluation of Multi-Object Filters," Signal Processing, IEEE Transactions on , vol.56, no.8, pp.3447,3457, Aug. 2008.
This isn't necessary for the GM-PHD filter to work but provides a nice way of visualising and quantifying performance. If people want to tweak this code, the metric is a nice way to see how well things are performing and it is fairly easy to modify this implementation to perform multiple test runs with the same data but different filter parameters and compare the metrics.

Read the README.txt, or just start with GM_PHD_Filter.m and work from there.


Error Ellipse inspired this file.

Required Products MATLAB
MATLAB release MATLAB 8.1 (R2013a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (16)
07 Aug 2016 Quang Lam

Hi Bryan,
Could you please email me ( the missing m file below:
Thank you,
Quang M Lam

Comment only
17 Jul 2014 Kai

Kai (view profile)

Hi Bryan,

To my understanding, new targets can only "birth" (assume no spawning) at predetermined positions defined in the birth intensity model.

As a result, how PHD can pick up the new targets "birth" at time step other than t=0. Assume measurements for these new targets are available.

Thanks in advance.

Comment only
25 Jun 2014 Fernando Iglesias

Just a very minor detail. As per MATLAB R2013b deg2rad (used in GM_EKF_PHD_Initialise_Jacobians) has been replaced by degtorad.

Comment only
12 Jun 2014 Kai

Kai (view profile)

Hi Bryan,

I am new to RFS and PDH filter. Could you please explain to me in the "GM_PHD_Initialisation.m" L147 and L147 about the target birth and spawn:
1) what does the weight "0.1 and 0.1" means in target birth.
2) what does the weight "0.05" means in target spawn.

Thanks a lot.

Comment only
13 Mar 2014 Kelin Lu

Hi, Bryan. I hava a question with the target modelling. In function `GM_PHD_Initialisation`, the standard deviation of process noise `sigma_v` is initialised. While in function `GM_PHD_Simulate_Measurements`, the functions for target 1, 2 and 3 are like `simTargetiState = F * simTarget1State`. It seems that the process noise has not been involved. And in function `GM_PHD_Predict_Birth`, there is also no process noise. Can you please tell me how did you treat with the process noise? Thank you.

Comment only
07 Jan 2014 Bryan Clarke

Hi Hannes,
Good find, it was used in the moving-target spawn of GM_PHD_Create_Birth at line 65 but was missed at the static-target spawn at line 116. It should be fixed now.

Comment only
17 Dec 2013 Hannes Pessentheiner

Hey Bryan,
I'm currently implementing a GM-CPHD Filter, and - just to find out the similarities - I've read some old papers about the standard GM-PHD Filter, e.g., 'The Gaussian Mixture Probability Hypothesis Density Filter', Vo et al., 2006 (as mentioned in your description). Additionally, I've compared your GM-PHD implementation with mine. I've realized that you do not consider the old target weights when you calculate the spawning probabilities (see eq. 28 and tab. 1, step 1, 2nd loop in Vo's paper).

It should be

w_spawn = w_beta * w_old * N(·)

instead of

w_spawn = w_beta * N(·).

Please, let me know if I'm wrong. Thanks!

Comment only
10 Dec 2013 Bryan Clarke

Hi Kelin,
That was an untidy bug that I left in there because I kept forgetting to fix it. It didn't break anything, it just meant that the simTarget3Vel was half what it should have been. Well spotted, it should be gone now.

Comment only
05 Dec 2013 Kelin Lu

Hi Bryan. Thanks for your sharing code. I think it's great! At the same time, I have a question in function `GM_PHD_Simulate_Measurements`. Why does `simTarget3State = F * simTarget3State` been executed twice after the `simTarget3SpawnTime` during one iteration? At the first time, it was executed because `simTarget3State` will not be empty(since it will be constructed at the `simTarget3SpawnTime`) after the `simTarget3SpawnTime` at line 23-25, then it will be executed again since `k > simTarget3SpawnTime` at line 31-33. Can you please explain to me if I'm missing something? Thank you.

20 Nov 2013 Kevin Nickels

20 Nov 2013 Kevin Nickels

12 Sep 2013 Bryan Clarke

Hey Hannes,
Good finds! They are fixed now, thank you for pointing them out.

Comment only
11 Sep 2013 Hannes Pessentheiner

Hey Bryan,
I think I've discovered a small bug in your merging/pruning code. Let's check out the code first:

for i = 1:length(I)
thisI = I(i);
P_range = calculateDataRange4(i);
delta_m = m_k(:,thisI) - m_k(:,j);
mahal_dist = delta_m' * inv(P_k(:,P_range)) * delta_m;
if(mahal_dist <= mergeThresholdU)
L = [L, thisI];

The 'inv(P_k(:,P_range))' for distance computation causes some troubles: first, 'inv' isn't an accurate matrix-inversion operation, use 'mldivide' or '\' instead; second, and that's the main reason of this message, you use the correct 'm_k(:,thisI)' for your distance computation, but not the correct covariance-matrices 'P_k(:,P_range)'. You always (!) take the matrices with column indices '1:4','5:8',etc., instead of column indices corresponding to 'I(i)'. Please, let me know if I'm wrong.

16 Aug 2013 Ying

Ying (view profile)

Thanks again, I will try my best to digest the program. Best wishes.

Comment only
15 Aug 2013 Bryan Clarke

Hi Ying. Currently the sensor is simulated as a direct measurement of target position in cartesian space; the measurement vector Z = [X; Y] where X = [x1, x2, ... xN], Y = [y1, y2, ... yN] for detected targets/clutter 1:N. See the file GM_PHD_Simulate_Measurements.m for how this is done. This is in line with the simulation described in Vo & Ma, and allows a linear Kalman filter to be used.
The same paper describes an extended Kalman filter EKF-PHD algorithm which can be used in nonlinear problems, such as those using range-bearing sensors. I will upload an implementation of this as a separate submission in a couple of days.

Comment only
14 Aug 2013 Ying

Ying (view profile)

Hi, I want to know how did you simulate the radar sensor in the program.

26 Jul 2013 1.2

Fixed typos in the description.

31 Jul 2013 1.3

New targets can now be initialised with a 'spawned' weight rather than a 'birthed' weight. Whichever weight is higher is the one that is used. Previously all targets were birthed.
Extra output statements are now output when VERBOSE is set.

27 Aug 2013 1.4

Added a few comments to hopefully make it clearer how to modify this code for different applications. Also did a bit of variable/renaming deleting but nothing major.

04 Sep 2013 1.6

Enhanced spawning for better compliance with Vo & Ma paper. Fixed error in spawn weight calculation. Minor plotting enhancements.

12 Sep 2013 1.7

Fixed a bug in GM_PHD_Prune. Changed inv to \ in pruning and constructing update components.

21 Nov 2013 1.8

A few bug fixes. A few new comments. Added a couple of features to better fit the Vo & Ma implementation. Added an implementation of the Optimal Subpattern Assignment (OSPA) metric proposed by Schuhmacher et al to quantify filter performance.

02 Dec 2013 1.9

I have included the extended Kalman filter algorithm described by Vo & Ma (the filter can be run as either linear KF or EKF). I also included a way of drawing velocity arrows on the main plot to demonstrate estimated target velocities.

05 Dec 2013 1.11

Fixed a small bug in a print statement.

10 Dec 2013 1.12

Modified GM_PHD_Create_Birth and GM_EKF_PHD_Create_Birth so that static targets are spawned from iteration 1 rather than waiting until iteration 2.

12 Dec 2013 1.13

Removed an error that was passing the OSPA metric the wrong data to calculate with. Fixed a bug that meant costs were not being calculated in the OSPA for measurements in Y unmatched in X. The results are much better now but still a bit buggy.

13 Dec 2013 1.14

Removed CalculateOSPAMetric.m, added Ba-Ngu Vo's implementation ospa_dist.m which uses Alex Melin's Hungarian.m

07 Jan 2014 1.15

Fixed a bug where, for some spawning settings, the weight of the existing target was not used to scale the weight of a target it spawned.

Contact us