echo on
% generamos el ruido blanco
Fm=8000;
w=randn(4*Fm,1);
k=-1023:1023;
a=0.8;
%generamos proceso de salida filtro H(z)=1/(1-az^(-1))

x=filter(1,[1 -a],w);

% calculamos correlaciones

rxb0=xcorr(x(1:1024),'biased');
rxb1=xcorr(x(4560:4560+1023),'biased');
rxu0=xcorr(x(1:1024),'unbiased');
rxu1=xcorr(x(4560:4560+1023),'unbiased');

rxteo=1/(1-a^2).*a.^abs(-1023:1023);

figure(1);clf;
plot(k,rxteo); title('Correlacion teorica');grid on;
pause
hold on
%con sesgo segmento 0
plot(k,rxb0,'r'); title('Correlacion con sesgo segmento inicial');pause
%con sesgo segmento 1
plot(k,rxb1,'y'); title('Correlacion con sesgo segmento intermedio');pause
%sin sesgo segmento 0
plot(k,rxu0,'m'); title('Correlacion sin sesgo segmento inicial');pause
%sin sesgo segmento 1
plot(k,rxu1,'g');title('Correlacion sin sesgo segmento intermedio');pause
% para hacer la fft utilizamos la funcion ifftshift para desplazar la parte
% negativa de la autocorrelacion a la mitad superior de la seņal


Sxb0=fft(ifftshift(rxb0));
Sxb1=fft(ifftshift(rxb1));
Sxu0=fft(ifftshift(rxu0));
Sxu1=fft(ifftshift(rxu1));
f=[0:0.5/1024:0.5*(1-1/1024)]*Fm;
Sxteo=abs(1./fft([1 -a],2048)).^2;
figure(2);clf; %Teorica
plot(f,10*log10(real(Sxteo(1:1024))));grid on;title('Densidades Espectrales de Potencia Teorica');
pause; hold on
%con sesgo segmento 0
plot(f,10*log10(real(Sxb0(1:1024))),'r');title('Densidades Espectrales de Potencia con sesgo');
pause
%con sesgo segmento 1
plot(f,10*log10(real(Sxb1(1:1024))),'y');title('Densidades Espectrales de Potencia con sesgo');
pause
%sin sesgo segmento 0
plot(f,10*log10(real(Sxu0(1:1024))),'m');title('Densidades Espectrales de Potencia sin sesgo');
pause
%sin sesgo segmento 1
plot(f,10*log10(real(Sxu1(1:1024))),'g');title('Densidades Espectrales de Potencia sin sesgo');
echo off