%% 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:

% Sennal ventana
M = 3;
m=-M+1:1:M-1;
wt=(M-abs(m))';
wt=wt/M;

% Estimas de la funcion de autocorrelacion en  m = -M+1:1:M-1;
rx = [3, 6, 9, 6, 3];

% Estimador Blackman-Tukey
Nfft = 8;
Sx=abs(fft(rx.*wt',Nfft));
k=0:1:Nfft-1;
wk=2*pi*k/Nfft;
figure(1);plot(wk,Sx,'bo-')
grid on

%%
%
%% Ejercicio 2:

% Realizacion del proceso estocastico
x = [1,2,0,-1,1,2,3,1,2,3,-2,-3,-2,-2,-1];
N = length(x);

% Estimador de Bartlett
K = 5;  % numero de muestras de cada segmento(bloque)
Ns = N/K; % numero de bloques (segmentos)
Nfft = 128;

Xs = reshape(x,K,Ns);
Sx_s=zeros(Ns,Nfft);
for ns=1:Ns
    xs=Xs(:,ns);
    X=fft(xs,Nfft);
    Sx_s(ns,:) = (1/K)*abs(X).^2;    
end
Sx_b = mean(Sx_s);

k=0:1:Nfft-1;
wk=2*pi*k'/Nfft;
figure(3);plot(wk,Sx_b,'k');
grid on; axis tight; xlabel('\omega')

%
%
%% Ejercicio 3:
%
% Estimas de la funcion de autocorrelacion en  m =-2, -1, 0, 1, 2

rx = [1/4, 1/2, 2, 1/2, 1/4]; 

% Frecuencias donde se evalua la DEP

Nfft = 128; 
k=0:1:Nfft-1;
wk=2*pi*k/Nfft;

%
% DEP:

Sx = (8-4*cos(wk))./(5-4*cos(wk));

% Correlograma:

Sx_c=abs(fft(rx,Nfft));

% Estimador parametrico AR(2)

R = [rx(3) rx(4); rx(4) rx(3)];
r = [rx(4) rx(5)]';
a=-inv(R)*r; % coeficientes del filtro
v_r=rx(3)+a'*r; % varianza del ruido
v_r=abs(v_r);
H=freqz(1,[1;a],wk); 
Sx_AR2=v_r*abs(H).^2;

% Estimador parametrico AR(1)

R = rx(3);
r = rx(4);
a=-inv(R)*r; % coeficiente del filtro
v_r=rx(3)+a'*r; % varianza del ruido
v_r=abs(v_r);
H=freqz(1,[1;a],wk); 
Sx_AR1=v_r*abs(H).^2;

%
% Dibuja la DEP y las estimas

figure(4);plot(wk,Sx,'b-',wk,Sx_c,'r-',wk,Sx_AR1,'g-',wk,Sx_AR2,'m-')
legend('DEP','Correlograma','AR(1)','AR(2)')
grid on; axis tight; xlabel('\omega')

%
%
%% Ejercicio 4:
%
% a) En general, la función de autocorrelación de un proceso MA(2) es
% r_X(m) = b_0\ b_0 \ \sigma^2 \ \delta(m) + b_0 \ b_1 \ \sigma^2 \ \delta(m+1) + b_0 \ b_2 \ \sigma^2\ \delta(m+2)  + \\
% \quad \quad \quad \quad  b_1\ b_0 \ \sigma^2 \ \delta(m-1) + b_1 \ b_1 \ \sigma^2 \ \delta(m) + b_1 \ b_2 \ \sigma^2 \ \delta(m+1) + \\
% \quad \quad \quad \quad  b_2\ b_0 \ \sigma^2 \ \delta(m-2) + b_2 \ b_1 \ \sigma^2 \ \delta(m-1) + b_2 \ b_2 \ \sigma^2 \ \delta(m)
% donde \sigma^2 es la varianza del ruido a la entrada del filtro, y b_0, b_1 y b_2son los coeficientes del filtro.
% Por tanto, la funcion de autocorrelación en  m=0,1,2 es

rx = [3,2,1]';

% y vale cero para m > 2.
% DEP:

wk = -pi:pi/100:pi;
Sx = 3 + 4*cos(wk) + 2*cos(2*wk);

% Estimador parametrico AR(2)

R = [rx(1) rx(2); rx(2) rx(1)];
r = [rx(2) rx(3)]';
a = -inv(R)*r; % coeficientes del filtro
v_r=rx(1)+a'*r; % varianza del ruido
v_r=abs(v_r);
H=freqz(1,[1;a],wk); 
Sx_AR2=v_r*abs(H).^2;

%
% Estimador parametrico AR(p)

p = 8;
rx = [rx;zeros(p-2,1)];

R=toeplitz(rx(1:p));
r=rx(2:p+1);
a = -inv(R)*r; % coeficientes del filtro
v_r=rx(1)+a'*r; % varianza del ruido a la entrada del filtro
v_r=abs(v_r);
H=freqz(1,[1;a],wk); 
Sx_ARp=v_r*abs(H).^2;

%
% Dibuja la DEP y la estima AR(2)

figure(4);plot(wk,Sx,'b-',wk,Sx_AR2,'r-',wk,Sx_ARp,'g-')
legend('DEP','AR(2)','AR(p)')
grid on; axis tight; xlabel('\omega')

%
%%