%clear
%RSB=0;              %défini comme le rapport vitesse (première position) sur écart type du bruit
%N=5;                % nb de points

% Indicateurs : 
MV=1;               % calcule les estimées du MV   
MAP=1;              % calcule les estimées du MAP
FIGS=0;             % Affiche les figures des différentes densités de proba
display_results=1;  % Affiche les résultats numériques


if display_results
disp('*****************************************************');
disp(['Nombre de points : ' num2str(N) ' RSB=' num2str(RSB) ' dB'])
disp('*****************************************************');
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Mise en place des données : trajectoire 1D à vitesse constante
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

vv=50*1000/3600;    % 50 km/h --> conversion en m/s
xo=0;               % position initiale

t=[1:N]';            
x=vv*t+xo;          %trajectoire non bruitée

%bruit gaussien sur la trajectoire
sigma=vv*10^(-RSB/10);
R=sigma^2*eye(N);   % bruit gaussien décorrélé
bgx=sigma*randn(N,1);

%trajectoire bruitée 
y=x+bgx;            % en mètres

if MV
%*************************************************************
% Critère du maximum de vraisemblance 
%*************************************************************

v=[0:.1:120]*1000/3600;         % Pour tracé des figures
for i=1:length(v)
%!!! A COMPLETER
end
% !!! NORMALISATION A COMPLETER

% Calcul de l'estimée (cf TD)
v_MV=%!!! A COMPLETER         % en m/s
sigma2_MV=%!!! A COMPLETER


% max de la courbe
[Val Ind]=max(p_y_si_v); %(en km/h)
v_MV2=v(Ind);

if display_results
disp(['Estimée de la vitesse au max de vraisemblance (calcul) : ' num2str(v_MV*3.6) '(km/h)'])
disp(['Estimée de la vitesse au max de vraisemblance (max de la courbe) :' num2str(v_MV2*3.6) '(km/h)'])
disp(['Estimée de la variance au max de vraisemblance (calcul) : ' num2str(sigma2_MV) ])
end; % du if display_results

end; % du if MV
%
if MAP
%*************************************************************
% Critère MAP
%*************************************************************
%Loi a priori: vitesse suit loi de Gauss centrée sur 50 km/h
sigmav=5/3.6;
m=50/3.6;
pv=exp(-(v-m).^2/(2*sigmav^2));     %~loi pour vitesse en ville
pv=pv/sum(pv);

%loi a posteriori
%!!! A COMPLETER
% !!! NORMALISATION A COMPLETER
%max critère 


% Calcul de l'estimée (cf TD)
v_MAP=    %!!A COMPLETER        % en m/s

% max de la courbe
[Val Ind]=max(p_v_si_y); %(en m/s)
v_MAP2=v(Ind);
if display_results
disp(['Estimée du max a posteriori (calcul) : ' num2str(v_MAP*3.6) '(km/h)'])
disp(['Estimée du max a posteriori (max de la courbe) :' num2str(v_MAP2*3.6) '(km/h)'])
end; % du if display_results
end; % du if MAP


%
if FIGS
% Graphiques
figure(1);
subplot(211);plot(t,x,t,y);legend('trajectoire non bruitée',['trajectoire bruitée avec RSB = ',num2str(RSB)])
subplot(212); plot(v*3.6,p_y_si_v); legend('Vraisemblance')
titre=['vitesse estimée MV = ',num2str(v_MV*3.6)];
title(titre);

figure(2)
subplot(211);plot(v*3.6,pv,v*3.6,p_y_si_v,':');
legend('loi vitesse','vraisemblance')
subplot(212);plot(v*3.6,p_v_si_y)
legend('loi a posteriori')
titre=['vitesse estimée MAP = ',num2str(v_MAP*3.6)];
title(titre);
end % du if figs

