% TP Signaux aléatoires - SRT4
%% TP5 Propriété ergodique d'un signal aléatoire
%
% Nom - Prénom : Avenel - Théophile
%
% On considère le signal aléatoire en temps discret X(n)=\cos(2\pi f n + \phi) où f=1/128 et \phi est une phase aléatoire (variable aléatoire) uniformément répartie sur [0,2\pi[(fonction rand).\\
% Déterminez si X(n) est stationnaire au premier ordre. Calculer la moyenne statistique d'ordre 1 de X(n). La séquence aléatoire X(n) est-elle ergodique relativement à sa moyenne ?
% Générez 4 réalisations de ce signal sur 1024 points. A partir de ces 4 réalisations, estimez la moyenne statistique à chaque instant n (en utilisant la fonction mean), puis déterminez les moyennes temporelles (toujours avec la fonction mean) pour chacune des 4 réalisations. Que constatez-vous ?
% Générez maintenant 40 puis 400 réalisations de 1024 points de ce signal X(n). Dans chaque cas, estimez à nouveau les moyennes statistiques à chaque instant, et calculez les moyennes temporelles. Que constatez-vous ? Est-ce en accord avec la théorie ?
%


close all;
clear all;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Ergodisme
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%

% E(X(n))=\int_{-\infty}^\infty \cos(2\pi f n + \phi)f_\phi(\phi)d\phi=\frac{1}{2\pi}\int_{0}^{2\pi} \cos(2\pi f n + \phi)d\phi=0.
% E(X(n)) ne dépend pas de n, la séquence aléatoire est donc stationnaire au premier ordre.
% L'ergodisme relativement à la moyenne pour un processus stationnaire se traduit par l'égalité de la moyenne temporelle et de la moyenne statistique. Il est facile de vérifier que la moyenne temporelle d'une réalisation est bien égale à zéro dans notre cas :
% \left\langle X(n)\right\rangle=\lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_{-N}^N\cos(2\pi f n + \phi)=0
%
%


% Paramètres du signal
nbpoints = 1024;  % nombre de points (longueur du signal)
freq = 1/128;     % fréquence du signal
Points = [0:nbpoints-1];  % instants de temps
FPoints = freq * Points;  % pour calculer l'argument de la cosinus

% Génération de la première phase aléatoire uniforme
nbrealisations = 4;  % on commence avec 4 réalisations

% Initialisation des réalisations
realizations = zeros(nbrealisations, nbpoints);

% Génération des 4 réalisations
for i = 1:nbrealisations
    phi = 2 * pi * rand();  % phase aléatoire uniformément répartie entre 0 et 2*pi
    realizations(i, :) = cos(2 * pi * FPoints + phi);  % génération du signal
end

% Tracé des 4 réalisations
figure;
for i = 1:nbrealisations
    subplot(2, 2, i);
    plot(Points, realizations(i, :));
    title(['Réalisation ', num2str(i)]);
    xlabel('Temps (n)');
    ylabel('Amplitude');
end

% Calcul de la moyenne statistique à chaque instant
mean_statistique = mean(realizations, 1);  % moyenne sur les réalisations pour chaque n

% Tracé de la moyenne statistique
figure;
plot(Points, mean_statistique);
title('Moyenne statistique à chaque instant');
xlabel('Temps (n)');
ylabel('Moyenne statistique');

% Calcul de la moyenne temporelle pour chaque réalisation
mean_temporelle = zeros(1, nbrealisations);  % pour stocker les moyennes temporelles

for i = 1:nbrealisations
    mean_temporelle(i) = mean(realizations(i, :));  % moyenne temporelle de chaque réalisation
end

% Affichage des moyennes temporelles
disp('Moyenne temporelle pour chaque réalisation :');
disp(mean_temporelle);

% Génération de 40 réalisations
nbrealisations = 40;
realizations_40 = zeros(nbrealisations, nbpoints);

for i = 1:nbrealisations
    phi = 2 * pi * rand();
    realizations_40(i, :) = cos(2 * pi * FPoints + phi);
end

% Calcul de la moyenne statistique pour 40 réalisations
mean_statistique_40 = mean(realizations_40, 1);

% Calcul des moyennes temporelles pour 40 réalisations
mean_temporelle_40 = zeros(1, nbrealisations);
for i = 1:nbrealisations
    mean_temporelle_40(i) = mean(realizations_40(i, :));
end

% Génération de 400 réalisations
nbrealisations = 400;
realizations_400 = zeros(nbrealisations, nbpoints);

for i = 1:nbrealisations
    phi = 2 * pi * rand();
    realizations_400(i, :) = cos(2 * pi * FPoints + phi);
end

% Calcul de la moyenne statistique pour 400 réalisations
mean_statistique_400 = mean(realizations_400, 1);

% Calcul des moyennes temporelles pour 400 réalisations
mean_temporelle_400 = zeros(1, nbrealisations);
for i = 1:nbrealisations
    mean_temporelle_400(i) = mean(realizations_400(i, :));
end

% Affichage des résultats pour 40 et 400 réalisations
disp('Moyenne statistique pour 40 réalisations :');
disp(mean_statistique_40);

disp('Moyenne temporelle pour 40 réalisations :');
disp(mean_temporelle_40);

disp('Moyenne statistique pour 400 réalisations :');
disp(mean_statistique_400);

disp('Moyenne temporelle pour 400 réalisations :');
disp(mean_temporelle_400);

% Tracé de la moyenne statistique pour 400 réalisations
figure;
plot(Points, mean_statistique_400);
title('Moyenne statistique à chaque instant (400 réalisations)');
xlabel('Temps (n)');
ylabel('Moyenne statistique');


% On remarque que la moyenne temporelle est proche de la valeur théorique zero tandis que l'estimation de la moyenne d'ensemble est loin d'etre égale à zéros pour tout n lorsque nbrealisations=4.
% L'estimation de cette moyenne d'ensemble n'est pas fiable puisqu'elle ne repose pas sur un nombre suffisant de réalisations.
% Plus ne nombre de réalisations est grand, plus la moyenne est proche de zero quel que soit n.