% este script mostra a estimativa da densidade
% espectral de potência por meio do cálculo da 
% autocorrelação.
% é gerado um sinal binário com pulsos retangulares
% cuja autocorrelação deveria ser uma função triangular
% e cuja densidade espectral de potência um sinc^2
%
% Universidade de Brasília
% André Noll Barreto, 05 maio 2011

close all

%%%%%%%%%%%%%
% parâmetros
rs = 8; % amostras por pulso
ncorr = 512; % número de amostras da autocorrelação


%%%%%%
% é inicialmente gerado um sinal com apenas 1000 pulsos
% e são calculadas a autocorrelação e a PSD

N=1000;

bits = 1.0*(rand(1,N) > 0.5);

p_nrz = ones(rs,1);

t = 0:1/rs:(N-1/rs);
t_s = 1/rs;


s_nrz_pol = p_nrz * (2*bits-1);
s_nrz_pol = reshape(s_nrz_pol, 1 , rs * N);



figure(1)

plot(t,s_nrz_pol)
axis([0,40,-1.1,1.1]);
title('sinal')
xlabel('t(s)')


figure(2)
autoc = xcorr(s_nrz_pol, s_nrz_pol, ncorr);
plot((-ncorr:ncorr)*t_s, autoc);
title('autocorrelação')
xlabel('t(s)')

f0 = 1 / (ncorr * 2 * t_s);
f = (0:(ncorr*2)) * f0; % frequências amostradas na DFT



figure(3);

psd_nrz_polar = fft( autoc );

plot(f, abs(psd_nrz_polar))
title('Densidade Espectral de Potência')
xlabel('frequência (Hz)');


%%%%%%%
% é gerado agora um sinal com 10000 pulsos
% e são calculadas a autocorrelação e a PSD

N=10000;


bits = 1.0*(rand(1,N) > 0.5);

p_nrz = ones(rs,1);

t = 0:1/rs:(N-1/rs);
t_s = 1/rs;


s_nrz_pol = p_nrz * (2*bits-1);
s_nrz_pol = reshape(s_nrz_pol, 1 , rs * N);



figure(4)

plot(t,s_nrz_pol)
axis([0 40 -1.1 1.1])
title('sinal')
xlabel('t(s)')



figure(5)
autoc = xcorr(s_nrz_pol, s_nrz_pol, ncorr);
plot((-ncorr:ncorr)*t_s, autoc);
title('autocorrelação')
xlabel('t(s)')

f0 = 1 / (ncorr * 2 * t_s);
f = (0:(ncorr*2)) * f0; % frequências amostradas na DFT



figure(6);

psd_nrz_polar = fft( autoc );

plot(f, abs(psd_nrz_polar))
title('Densidade Espectral de Potência')
xlabel('frequência (Hz)');



%%%%
% é agora gerado um sinal periódico
% e calculadaa a autocorrelação e a PSD

N=1000;

g = cos(2*pi*t) + .5*sin(4*pi*t);


figure(7)

plot(t,g)
axis([ 1 40 -1.1 1.1]);
title('sinal')
xlabel('t(s)')


figure(8)
autoc = xcorr(g, g, ncorr);
plot((-ncorr:ncorr)*t_s, autoc);
title('autocorrelação')
xlabel('t(s)')

f0 = 1 / (ncorr * 2 * t_s);
f = (0:(ncorr*2)) * f0; % frequências amostradas na DFT


figure(9);

psd_g = fft(autoc);

plot(f, abs(psd_g))
title('Densidade Espectral de Potência')
xlabel('frequência (Hz)');



