function [xLocNow, weightsNow] = particleFilter(xLocPrev, deltaY, phi,...
         sigma_v, N,sideLength,weightsPrev,linkLocs,linkDist,sigma)
% PURPOSE : Perform particle filtering
% INPUTS  : - xLocPrev = particle estimations in previous run. 
%           - deltaY = RSS changes.
%           - phi = the attenuation paramter in pixel-free model.
%           - sigma_v = transition noise.
%           - N = number of particles per run.
%           - sideLength = side length of a square layout.
%           - weightsPrev = particle weights in previous run.
%           - linkLocs = locations of node pair for each link.
%           - linkDist = Distance of a link.
%           - sigma = variable controls rate of decay in pixel-free model
% OUTPUTS : - xLocNow = particle estimations in new run.
%           - weightsNow = particle weights in new run.
% AUTHORS  : Xi Chen and Yunpeng Li
%            @ECE Department, McGill University
%            http://www.tsp.ece.mcgill.ca/Networks/index.html
% DATE     : Jun 29, 2011


%Target motion dynamics 
xParticlesThisRound = xLocPrev + sigma_v *...
    randn(size(xLocPrev));
%Generate approximat measurement \y
yParticles = phi * simData(xParticlesThisRound,...
    linkLocs,linkDist,sigma); 
%Calculate weights
meanDeltaY = repmat(deltaY, [1 N]);
weightings = -(yParticles - meanDeltaY).^2;

logWeights= sum(weightings);
logWeights = logWeights - max(logWeights);
weightsNow = weightsPrev .* exp(logWeights);
%Normalized the weights
weightsNow = weightsNow/sum(weightsNow);
%Multinomial Resampling Process
outIndex = multinomialR(1:N,weightsNow');

%% Auxiliary part
xLocNow(:,:) = xParticlesThisRound(:,outIndex);
yNewAPFParticles(:,:) = yParticles(:,outIndex);

yAPFParticles = phi * simData(xLocNow,...
    linkLocs,linkDist,sigma);

% Calculate the final weights
weightings = -(yAPFParticles - meanDeltaY).^2;

logValueVec = sum(weightings);
maxLogValueVec = max(logValueVec);
logValueVec = logValueVec - maxLogValueVec;

weightings = -(yNewAPFParticles - meanDeltaY).^2;

logValueVecNew = sum(weightings);
maxLogValueVecNew = max(logValueVecNew);
logValueVecNew = logValueVecNew - maxLogValueVecNew;
likelihood = exp(logValueVec);
likelihoodNew = exp(logValueVecNew);
weightsNow = likelihood ./ likelihoodNew;
%Normalize the weights
weightsNow = weightsNow / sum(weightsNow);

%% Draw the dynamic tracking results
numPixelRow = 100;
numPixelCol = 100;
xShow = zeros(numPixelRow, numPixelCol);

for i = 1:N
    xPos = round(xLocNow(1,i)/sideLength * numPixelCol);
    yPos = round(xLocNow(2,i)/sideLength * numPixelRow);
    if xPos > numPixelCol
        xPos = numPixelCol;
    end
    if xPos < 1
        xPos = 1;
    end
    if yPos > numPixelRow
        yPos = numPixelCol;
    end
    if yPos < 1
        yPos = 1;
    end
    xShow(numPixelRow + 1 - yPos,xPos) = xShow(...
        numPixelRow + 1 - yPos,xPos) + 1;
end

figure(11)
movegui('center');
colormap summer;
set(gca,'units','points'); 
imagesc(xShow);
colorbar;