Help with particle filter for measurement denoising

I'm attempting to implement a simple particle filter for measurement denoising and compare the result with least mean squares.

The particle filter is tracking, more or less, the measurement but the performance is not good enough. What am I doing wrong?

clear all;
close all;
clc;

% 1000 measurements of a variable, presenting REAL VALUE 7 and corrupted
% by gaussian noise with standard deviation 3
M = 1000;
stddev = 3;
Y = randn(M,1)*stddev + 7;

% LEAST MEAN SQUARES
fi = ones(M,1);     % constant
teta=(fi'*fi)\fi'*Y % RESULT

% PARTICLE FILTER
N = 10;               % number of particles
P = randn(N,1)*1 + 1; % initial particle set with variance and value

Yhist = []; % History of the average value of Y
Phist = []; % History of the average value of P

for i = 1:M-N
    % Particle generation
    Ph = randn(N,1)*stddev + P;

      % Weight calculated based on the measured error between the particle
      % set and the current mean value of the measurement
      Pw = normpdf(mean(Y(1:i))-Ph, 0, stddev);

      % Weight normalized to form a probability distribution (i.e. sum to 1).
      Pw = Pw./sum(Pw);

      % Multinomial resampling
      idx = randsample(N,N,true,Pw);      
      P(idx) = Ph(idx);

      % Outputs
      Yhist = [Yhist; mean(Y(1:i))];
      Phist = [Phist; mean(P)];
  end

% RESULT
mean(P)

figure;
plot(1:size(Phist,1),Phist,'k');
hold all;
plot(1:size(Phist,1),Yhist,'r');
legend('Particle Filter', 'Measurement');

return

[JP Norair]

This comment is mostly for posterity given the time lag, but you may want to generate more than 10 particles. Try 100 or 1000.

Answer 1

[Hitesh]

Hi Daniel,

There are some conceptual and implementation issues that are causing the poor performance. The performance is enhanced by addressing the following points.

Particle Propagation

Ph = randn(N,1)*stddev + P;

• This was adding measurement noise (stddev) to each particle every iteration.
• In a typical particle filter, you need to add processed noise to simulate the evolution of the hidden state.

Weighting

Pw = normpdf(mean(Y(1:i))-Ph, 0, stddev);

• You were comparing the mean of all measurements up to time i to each particle. But the particle filter is a sequential (online) estimator.

Pw = normpdf(Y(i) - Ph, 0, stddev);

• You need to use the current measurement (Y(i)) only.

Resampling

The resampling logic was incorrect as it was not updating P correctly

idx = randsample(N,N,true,Pw);      
P(idx) = Ph(idx);

Corrtected :

P = Ph(idx);

Averaging

• You were tracking the mean of all measurements up to i (mean(Y(1:i))), but for fair comparison, you need to compare the filtered estimate at each step to the current measurement and the true value.

Particle Count & Initialization:

• Particle Count: N=10 is very low. Increase it to N=100 or more for better results.
• Initialization: You initialized particles around 1, but the true value is 7. This makes convergence slower.

For more information regarding the "Particle Filter Workflow", kindly refer to the following MATLAB documentation:

https://www.mathworks.com/help/nav/ug/particle-filter-workflow.html