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

x=filter(1,conv([1 -a],[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=ifftshift(real(ifft(1./fft((conv(conv([1 -a],[1 -a']),conv([1 -1/a],[1 -1/a']))),2048))));
rxteo=ifftshift(real(ifft(1./fft([abs(a)^2, 0, (1+abs(a)^4), 0 ,abs(a)^2],2048))));
figure(1);clf;
plot([-1022:1025],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(conv([1 -a],[1 -a']),2048)).^2;
figure(2);clf; %Teorica
plot(f,real(Sxteo(1:1024)));grid on;title('Densidades Espectrales de Potencia Teorica');
pause; hold on
%con sesgo segmento 0
plot(f,real(Sxb0(1:1024)),'r');title('Densidades Espectrales de Potencia con sesgo');
pause
%con sesgo segmento 1
plot(f,real(Sxb1(1:1024)),'y');title('Densidades Espectrales de Potencia con sesgo');
pause
%sin sesgo segmento 0
plot(f,real(Sxu0(1:1024)),'m');title('Densidades Espectrales de Potencia sin sesgo');
pause
%sin sesgo segmento 1
plot(f,real(Sxu1(1:1024)),'g');title('Densidades Espectrales de Potencia sin sesgo');
echo off