% Proceso ARMA

% Datos
b=[1 1]; % coeficientes filtro MA
a=1; % coeficientes filtro MA
Lf=length(b); % Longitud del filtro
varianza=1; % varianza de ruido
N=100; % longitud de la señal
Nr=1e6;

% Realizaciones de 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

% Dibuja una realización del proceso estocastico x[n]
figure(1);subplot(2,1,1);stem((0:N-1),U(Nr,:))
figure(1);subplot(2,1,1);title('Realizacion del proceso a la entrada (WGN) : u[n]')
figure(1);subplot(2,1,2);stem((0:N-1),X(Nr,:))
figure(1);subplot(2,1,2);title('Realizacion del proceso a la salida (ARMA): x[n]')
figure(1);xlabel('n')

% Dibuja la función de medias
medias=mean(X);
figure(2);stem((0:N-1),medias)
figure(2);title('función de medias')

% Dibuja las funciones de autocorrelacióna la entrada y a la salida
n=50;
m=-20:20;
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
figure(3);subplot(2,1,1);stem(m,r_u)
figure(3);title('función de autocorrelación: r_u[m]')
figure(3);subplot(2,1,2);stem(m,r_x)
figure(3);title('función de autocorrelación: r_x[m]')

% Dibuja la función de correlación entre u[n] y x[n]
n=50;
m=-20:20;
r_ux=zeros(length(m),1);
for im=1:length(m)
    r_ux(im)=mean(U(:,n+m(im)).*X(:,n));
end
figure(4);stem(m,r_ux)
figure(4);title('función de correlación entrada-salida:  r_{ux}[m]')

% Densidad espectral de potencia
omega=-pi:pi/1000:pi;
S_x=freqz(r_x,1,omega);
S_u=freqz(r_u,1,omega);
figure(5);plot(omega,abs(S_u),'b',omega,abs(S_x),'r');
figure(5);xlabel('w')
figure(5);legend('S_u(w)','S_x(w)')
figure(5);title('Densidades espectrales de potencia a la entrada (WGN) y a la salida (ARMA)')

