%% Ejercicios de análisis espectral de procesos estocásticos ergódicos
% 2022 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
%
%
%% Ejercicio 1:
%
% Proceso estocástico coseno con fase aleatoria
%
% Datos:

A=1;
w=pi/8;
N=128; % longitud de la sennal
Nr=1e4;

% Realizaciones:

X=zeros(Nr,N);
Phi=2*pi*rand(Nr,1);
for nr=1:Nr
    x=A*cos(w*(0:N-1)+Phi(nr));    
    X(nr,:)=x;
end

% Dibuja una realizacion del proceso estocastico:

figure(1);stem((0:N-1),X(Nr,:))
figure(1);xlabel('n')
figure(1);ylabel('x[n]')
figure(1);title('Una realizacion del proceso estocastico tono con fase aleatoria')
figure(1);axis tight

% Dibuja la funcion de medias:

medias=mean(X);
figure(2);stem((0:N-1),medias)
figure(2);title('funcion de medias:  \mu_X[n]')
figure(2);axis tight

% Dibuja la funcion de autocorrelacion (suponiendo que es WSS):

n=N/2;
m=-(N/2-1):(N/2-1);
rx=zeros(length(m),1);
for im=1:length(m)
    rx(im)=mean(X(:,n+m(im)).*X(:,n));
end
figure(3);stem(m,rx)
figure(3);title('funcion de correlacion:   r_X[m]')
figure(3);axis tight

% Estima la densidad espectral de potencia como transformada de Fourier de la estima de la función de correlación:

Nfft=8*N;
Sx2=abs(fft(rx,Nfft));
wk2=2*pi*(0:Nfft-1)/Nfft;

% Estima la densidad espectral de potencia promediando los espectros de las realizaciones:

Nfft=8*N;
Xf=fft(X',Nfft);
Xf=conj(Xf');
Sx1=mean((1/N)*abs(Xf).^2);
wk1=2*pi*(0:Nfft-1)/Nfft;

% Dibuja las estimas de la DEP:

figure(4);plot(wk1,Sx1,'r',wk2,Sx2,'b')
figure(4);title('Estima de la DEP')
figure(4);legend('Promediando espectros','TF de la estima de la funcion de correlacion')

%%
%
%% Ejercicio 2:
%
% Suma de procesos estocasticos: Z[n] = X[n] + Y[n], donde X[n] es un tono con fase aleatoria e Y[n] es WGN.
%
% Datos:

A=1; % amplitud coseno
w=pi/16; % frecuencia coseno
varianza=1; % varianza ruido
N=100; % longitud de la se?al
Nr=1e6;

% Realizaciones deX[n]:

X=zeros(Nr,N);
Phi=2*pi*rand(Nr,1);
for nr=1:Nr
    x=A*cos(w*(0:N-1)+Phi(nr));    
    X(nr,:)=x;
end

% Realizaciones de Y[n]:

Y=sqrt(varianza)*randn(Nr,N);

% Realizaciones de Z[n]:

Z=X+Y;

% Dibuja una realizacion del proceso estocastico Z[n]:

figure(1);stem((0:N-1),Z(Nr,:))
figure(1);xlabel('n')
figure(1);ylabel('z[n]')
figure(1);title('Una realizacion del proceso estocastico Z[n]')

% Dibuja la funcion de medias:

medias=mean(Z);
figure(2);stem((0:N-1),medias)
figure(2);title('funcion de medias:   \mu_Z[n]')

% Dibuja la funcion de autocorrelacion de Z[n] (suponiendo que es WSS)

n=50;
m=-20:20;
rz=zeros(length(m),1);
for im=1:length(m)
    rz(im)=mean(Z(:,n+m(im)).*Z(:,n));
end
figure(3);stem(m,rz)
figure(3);title('funcion de correlacion:  r_Z[m]')

% Dibuja la funcion de correlacion entre Z[n] y X[n]:

n=50;
m=-20:20;
rzx=zeros(length(m),1);
for im=1:length(m)
    rzx(im)=mean(Z(:,n+m(im)).*X(:,n));
end
figure(4);stem(m,rzx)
figure(4);title('funcion de correlacion entre z[n] y x[n]:  r_{Z,X}[m]')

% Dibuja la funcion de correlacion entre Z[n] e Y[n]:

n=50;
m=-20:20;
rzy=zeros(length(m),1);
for im=1:length(m)
    rzy(im)=mean(Z(:,n+m(im)).*Y(:,n));
end
figure(5);stem(m,rzy)
figure(5);title('funcion de correlacion entre z[n] y y[n]:  r_{Z,Y}[m]')

% Estima la densidad espectral de potencia de Z[n] como transformada de Fourier de la estima de la funcion de correlacion

w=-pi:pi/128:pi;
Sz=freqz(rz',1,w);
figure(6);plot(w,abs(Sz),'b')
figure(6);xlabel('\omega')
figure(6);ylabel(' S_Z(\omega) ')
figure(6);title('Estima de la DEP de z[n]:   S_Z(\omega)')
figure(6);axis tight

%
%
%% Ejercicio 3:
%
% Proceso ARMA
%
% Datos:

b=[1/3 1/3 1/3]; % coeficientes filtro 
a=[1]; % coeficientes filtro 
varianza=1; % varianza de ruido a la entrada
N=512; % longitud de las realuzaciones
Nr=1e4; % no. de realizaciones

% Realizaciones de U[n] y X[n]:

U=zeros(Nr,N);
X=zeros(Nr,N);
for nr=1:Nr
    u=sqrt(varianza)*randn(1,N);
    x=filter(b,a,u);    
    X(nr,:)=x;
    U(nr,:)=u;
end

% Estima la funcion de medias:

medias=mean(X);

% Estima las funciones de autocorrelacion a la entrada y a la salida:

n=N/2;
M=20;
m=-M+1:M-1;
r_x=zeros(length(m),1);
r_u=zeros(length(m),1);
for im=1:length(m)
    r_x(im)=mean(X(:,n+m(im)).*X(:,n));
    r_u(im)=mean(U(:,n+m(im)).*U(:,n));
end

% Estima la funcion de correlacion cruzada entre U[n] y X[n]:

r_ux=zeros(length(m),1);
for im=1:length(m)
    r_ux(im)=mean(U(:,n+m(im)).*X(:,n));
end

% Estima las densidades espectrales de potencia haciendo la TF de la estima de la funcion de autocorrelacion:

w2=-pi:pi/128:pi;
Sx2=freqz(r_x,1,w2);
Su2=freqz(r_u,1,w2);

% Estimo la densidad espectral de potencia promediando los espectros de las realizaciones:

w1=-pi:pi/128:pi;
Sx1=zeros(1,length(w1));
Su1=zeros(1,length(w1));
for nr=1:Nr
    Xf = freqz(X(nr,:),1,w1);
    Sx1 = Sx1 + (1/N)*abs(Xf).^2;
    Uf = freqz(U(nr,:),1,w1);
    Su1 = Su1 + (1/N)*abs(Uf).^2;
end
Sx1=Sx1/Nr;
Su1=Su1/Nr;

% Respuesta en frecuencia del filtro:

w=-pi:pi/1000:pi;
H=freqz(b,a,w);

% Dibuja una realizacion de los procesos estocasticos:

nfig=1;
figure(nfig);subplot(2,1,1);stem((0:N-1),U(Nr,:))
figure(nfig);subplot(2,1,1);title('Realizacion del proceso a la entrada (WGN) : u[n]')
figure(nfig);subplot(2,1,2);stem((0:N-1),X(Nr,:))
figure(nfig);subplot(2,1,2);title('Realizacion del proceso a la salida (ARMA): x[n]')
figure(nfig);subplot(2,1,1);xlabel('n')
figure(nfig);subplot(2,1,1);ylabel('u[n]')
figure(nfig);subplot(2,1,2);xlabel('n')
figure(nfig);subplot(2,1,2);ylabel('x[n]')
figure(nfig);subplot(2,1,1);axis tight
figure(nfig);subplot(2,1,2);axis tight

% Dibuja la funcion de medias:

nfig=nfig+1;
figure(nfig);stem((0:N-1),medias)
figure(nfig);title('funcion de medias')
figure(nfig);axis tight

% Dibuja las funciones de autocorrelacion a la entrada y a la salida:

nfig=nfig+1;
figure(nfig);subplot(2,1,1);stem(m,r_u)
figure(nfig);subplot(2,1,1);title('funcion de autocorrelacion: r_u[m]')
figure(nfig);subplot(2,1,1);axis tight
figure(nfig);subplot(2,1,2);stem(m,r_x)
figure(nfig);subplot(2,1,2);title('funcion de autocorrelacion: r_x[m]')
figure(nfig);subplot(2,1,2);axis tight

% Dibuja la funcion de correlacion entre U[n] y X[n]:

nfig=nfig+1;
figure(nfig);stem(m,r_ux)
figure(nfig);title('funcion de correlacion entrada-salida:  r_{ux}[m]')
figure(nfig); axis tight

% Dibuja la densidad espectral de potencia obtenida a partir de las TF de las estimas de la autocorrelacion :

nfig=nfig+1;
figure(nfig);plot(w2,abs(Su2),'b',w2,abs(Sx2),'r');
figure(nfig);xlabel('w')
figure(nfig);legend('S_u(w)','S_x(w)')
figure(nfig);title('Estima de las DEP a partir de las TF de las estimas de la autocorrelacion')
figure(nfig);axis tight

% Dibuja la densidad espectral de potencia promediando espectros de las realizaciones del proceso:

nfig=nfig+1;
figure(nfig);plot(w1,abs(Su1),'b',w1,abs(Sx1),'r');
figure(nfig);xlabel('w')
figure(nfig);legend('S_u(w)','S_x(w)')
figure(nfig);title('Estima de las DEP promediando espectros')
figure(nfig);axis tight

% Dibuja la respuesta en frecuencia del filtro:

nfig=nfig+1;
figure(nfig);plot(w,abs(H).^2,'b');
figure(nfig);xlabel('w')
figure(nfig);title('Respuesta en frecuencia del filtro:    |H(\omega)|^2')
figure(nfig);axis tight

%
%
%% Ejercicio 3 explotando la ergodicidad:
%
% Proceso ARMA. Calcula sus parametros considerando ergodicidad
%
% Datos:


b=[1/3 1/3 1/3];
a=[1]; % coeficientes filtro AR
varianza=1; % varianza de ruido
N=1e5; % longitud de la senal

% Realizacion:

u=sqrt(varianza)*randn(1,N);
x=filter(b,a,u);    

% Estima de las medias de los procesos a la entrada y a la salida del filtro considerando que son procesos ergodicos:

mx=mean(x);

% Estima de las funciones de autocorrelacion a la entrada y a la salida del filtro considerando que son procesoso ergodicos:

M=20;
rx=xcorr(x,M-1,'biased');
ru=xcorr(u,M-1,'biased');
rux=xcorr(u,x,M-1,'biased');
m=1-M:M-1;

% Estima las densidades espectrales de potencia haciendo la TF de la estima de la funcion de autocorrelacion:

w=-pi:pi/1000:pi;
Sx=freqz(rx,1,w);
Su=freqz(ru,1,w);

% Respuesta en frecuencia del filtro:

H=freqz(b,a,w);

% Dibuja una realizacion de los procesos:

nfig=1;
n=0:N-1;
ng=n(100:200);
figure(nfig);subplot(2,1,1);plot(ng,u(ng),'o-','MarkerFaceColor','r','MarkerEdgeColor','none')
figure(nfig);subplot(2,1,1);plot(ng,u(ng),'o-','MarkerFaceColor','r','MarkerEdgeColor','none')
figure(nfig);subplot(2,1,1);xlabel('n')
figure(nfig);subplot(2,1,1);ylabel('u[n]')
figure(nfig);subplot(2,1,1);axis tight
figure(nfig);subplot(2,1,2);plot(ng,x(ng),'o-','MarkerFaceColor','r','MarkerEdgeColor','none')
figure(nfig);subplot(2,1,2);plot(ng,x(ng),'o-','MarkerFaceColor','r','MarkerEdgeColor','none')
figure(nfig);subplot(2,1,2);xlabel('n')
figure(nfig);subplot(2,1,2);ylabel('x[n]')
figure(nfig);subplot(2,1,2);axis tight

% Dibuja las estimas de las funciones de correlacion:

nfig=nfig+1;
figure(nfig);subplot(3,1,1);stem(m,ru)
figure(nfig);subplot(3,1,1);title('Estima de la funcion de autocorrelacion a la entrada:  r_U[m]')
figure(nfig);subplot(3,1,1);xlabel('m')
figure(nfig);subplot(3,1,2);stem(m,rx)
figure(nfig);subplot(3,1,2);title('Estima de la funcion de autocorrelacion a la salida:  r_X[m]')
figure(nfig);subplot(3,1,2);xlabel('m')
figure(nfig);subplot(3,1,3);stem(m,rux)
figure(nfig);subplot(3,1,3);title('Estima de la funcion de correlacion entrada-salida:  r_{U,X}[m]')
figure(nfig);subplot(3,1,3);xlabel('m')

% Dibuja la estima de las DEP:

nfig=nfig+1;
figure(nfig);plot(w,abs(Su),'b',w,abs(Sx),'r');
figure(nfig);xlabel('w')
figure(nfig);legend('S_U(w)','S_X(w)')
figure(nfig);title('Estima de las DEP promediando espectros')
figure(nfig);axis tight

% Dibuja la respuesta en frecuencia del filtro:

nfig=nfig+1;
figure(nfig);plot(w,abs(H).^2,'b');
figure(nfig);xlabel('w')
figure(nfig);title('Respuesta en frecuencia del filtro:    |H(\omega)|^2')
figure(nfig);axis tight

%