TP Signaux aléatoires - SRT4

TP6 Fonction de corrélation et densité spectrale de puissance

Nom - Prénom :
Soit la séquence aléatoire , où , ϕ est une phase aléatoire (variable aléatoire) uniformément répartie sur (fonction rand), et est une séquence aléatoire à valeurs indépendantes gaussiennes (fonction randn).
Calculez l'espérance mathématique ainsi que la séquence de corrélation de cette séquence aléatoire . La séquence aléatoire est-elle stationnaire ? Précisez votre réponse.
Calculez numériquement et représentez la séquence de corrélation de (fonction xcorr). Le résultat obtenu est-il conforme à la théorie ? En particulier, la valeur de obtenue est-elle en accord avec la valeur théorique. Que représente cette dernière valeur ?
Calculez numériquement et représentez la densité spectrale de puissance moyenne du processus en passant par la relation de Wiener-Kintchine (en utilisant la fonction fft). Que met en évidence cette densité spectrale de puissance ?
 
 
close all;
clear all;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fonction de corrélation et densité spectrale de puissance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Compte tenu du premier exercice il vient que
En effet
D'autre part
avec
et
La statistique du bruit nous indique que si , et si .
Si ,
Si ,
La séquence aléatoire est par conséquent faiblement stationnaire.
% Paramètres
nbpoints = 1024; % Nombre de points
freq = 1/128; % Fréquence fondamentale
Points = 0:nbpoints-1; % Indices de temps
FPoints = freq * Points; % Points de fréquence
 
% Nombre de réalisations
nbrealisations = 4;
 
% Préallocation pour stocker les résultats des réalisations
Y_total = zeros(nbrealisations, nbpoints); % Stocke toutes les réalisations
Y = 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
Y(i, :) = cos(2 * pi * FPoints + phi); % génération du signal
end
 
Y, Y_c = zeros(nbrealisations, nbpoints); % Stocke toutes les réalisations
Y = 4×1024
-0.6984 -0.7327 -0.7652 -0.7959 -0.8246 -0.8514 -0.8761 -0.8987 -0.9191 -0.9374 -0.9533 -0.9670 -0.9783 -0.9873 -0.9939 -0.9981 -0.9999 -0.9993 -0.9963 -0.9909 -0.9831 -0.9729 -0.9604 -0.9456 -0.9285 -0.9091 -0.8876 -0.8639 -0.8382 -0.8104 -0.7807 -0.7491 -0.7157 -0.6805 -0.6438 -0.6054 -0.5657 -0.5245 -0.4821 -0.4385 -0.3939 -0.3483 -0.3019 -0.2548 -0.2070 -0.1588 -0.1101 -0.0612 -0.0122 0.0369 -0.7233 -0.7563 -0.7875 -0.8168 -0.8441 -0.8694 -0.8926 -0.9136 -0.9325 -0.9491 -0.9634 -0.9754 -0.9850 -0.9923 -0.9972 -0.9997 -0.9997 -0.9974 -0.9927 -0.9855 -0.9760 -0.9642 -0.9500 -0.9335 -0.9148 -0.8939 -0.8708 -0.8456 -0.8184 -0.7892 -0.7581 -0.7252 -0.6906 -0.6543 -0.6164 -0.5770 -0.5362 -0.4941 -0.4509 -0.4065 -0.3612 -0.3150 -0.2681 -0.2205 -0.1724 -0.1238 -0.0750 -0.0260 0.0231 0.0721 -0.5229 -0.4805 -0.4369 -0.3922 -0.3466 -0.3001 -0.2530 -0.2052 -0.1569 -0.1083 -0.0594 -0.0103 0.0388 0.0877 0.1365 0.1850 0.2330 0.2804 0.3272 0.3731 0.4182 0.4623 0.5052 0.5470 0.5874 0.6264 0.6639 0.6998 0.7340 0.7664 0.7970 0.8257 0.8524 0.8770 0.8995 0.9199 0.9380 0.9539 0.9675 0.9787 0.9876 0.9941 0.9982 0.9999 0.9992 0.9961 0.9906 0.9827 0.9725 0.9599 -0.8247 -0.7960 -0.7653 -0.7328 -0.6985 -0.6626 -0.6250 -0.5860 -0.5455 -0.5037 -0.4607 -0.4166 -0.3715 -0.3255 -0.2787 -0.2313 -0.1833 -0.1348 -0.0860 -0.0370 0.0120 0.0611 0.1100 0.1586 0.2069 0.2546 0.3018 0.3482 0.3938 0.4384 0.4820 0.5244 0.5655 0.6053 0.6436 0.6804 0.7156 0.7490 0.7806 0.8103 0.8381 0.8639 0.8875 0.9091 0.9284 0.9455 0.9604 0.9729 0.9831 0.9909
 
for realisation = 1:nbrealisations
Y_c(i,:) = xcorr(Y(i,:),'unbiased')
end
Unable to perform assignment because the size of the left side is 1-by-1024 and the size of the right side is 1-by-2047.
 
gammay = mean(Y_total,1);
xaxe = -1024:1:1013
figure;
plot(xaxe,Y_total)
 
% Boucle sur le nombre de réalisations
for realisation = 1:nbrealisations
phi = 2*pi*rand(1); % Phase aléatoire uniformément répartie sur [0, 2*pi]
% Générer un bruit gaussien N(1,1)
B = 1 + randn(1, nbpoints); % Bruit gaussien N(1,1)
% Générer la séquence aléatoire Y(n) pour chaque réalisation
Y = cos(2*pi*freq*Points + phi) + B; % Signal cosinus + bruit gaussien
Y_total(realisation, :) = Y; % Stocker la réalisation
% Calculer et afficher l'espérance mathématique E[Y(n)] pour cette réalisation
E_Y = mean(Y);
disp(['Espérance mathématique pour la réalisation ', num2str(realisation), ' : ', num2str(E_Y)]);
% Calculer et tracer la fonction d'autocorrélation pour cette réalisation
% Ici, on utilise la fonction de corrélation unbiased pour obtenir l'autocorrélation
[autocorr_Y, lags] = xcorr(Y, 'unbiased');
figure;
plot(lags, autocorr_Y);
title(['Fonction d''autocorrélation pour la réalisation ', num2str(realisation)]);
xlabel('Décalages');
ylabel('\Gamma_Y(n,n'')');
% Vérifier la valeur de l'autocorrélation au décalage 0 pour cette réalisation
Gamma_0 = autocorr_Y(lags == 0);
disp(['Autocorrélation au décalage 0 (Gamma_Y(0)) pour la réalisation ', num2str(realisation), ' : ', num2str(Gamma_0)]);
end
 
 
Le résultat obtenu est conforme à la théorie. La valeur en zero est la puissance totale moyenne du signal aléatoire, c'est à dire la puissance moyenne du cosinus qui vaut à laquelle s'ajoute celle du bruit gaussien qui vaut ; soit au total . On remarque que cette valeur est obtenue avec une imprécision qui se réduit si on augmente le nombre de points sur lequel le calcul est réalisé.
 
%densité spectrale de puissance
dsp = fftshift(fft(Gamma_0,2048))
freqs = linspace(-0.5, 0.5, length(dsp)); % Axe des fréquences allant de -0.5 à 0.5
 
% Tracer la DSP moyennée
figure;
plot(freqs, abs(dsp));
title('Densité spectrale de puissance (DSP) moyennée sur les réalisations');
xlabel('Fréquence normalisée');
ylabel('DSP');
grid on;
On voit que la puissance du signal aléatoire se répartit sur exactement deux fréquences : la fréquence nulle (puissance du bruit) et la fréquence . Sur la figure, l'axe des fréquences est donné relativement à la valeur (Une raie spectrale à l'abscisse serait à la fréquence ).