% amostragem.m
%
% este script demonstra a amostragem a diferentes taxas e a reconstrução do
% sinal por meio de um filtro passa-baixa ideal

% Universidade 
% André Noll Barreto, 18 abr 2011


% número de amostras a serem consideradas para a reconstrução do sinal
nsincs = 50;

close all

%%%%
% geração do sinal, composto de 4 senoides com frequências 0,25; 0,5; 0,75
% e 1 HZ (ou seja B=1Hz)

T0 = 1/1000;
t = -10:T0:10;
f0 = .25;
a1 = 1;
phi1 = -pi/2;
a2 = -.5;
phi2 = pi/6;
a3 = 2;
phi3 = pi/8;
a4 = -1;
phi4 = -3*pi/2;

x = a1*cos(2*pi*f0*t+phi1) +  a2*cos(2*pi*2*f0*t+phi2) + ...
    a3*cos(2*pi*3*f0*t+phi3) + a4*cos(2*pi*4*f0*t+phi4);

%%%%
% amostragem com taxa à frequência de Nyquist
figure(1)

plot(t,x)
hold on

B = 4*f0;
BNyq = 2*B;
Ba = BNyq;
Ta = 1/Ba; 

tsamp = -10:Ta:10;
os = Ta/T0;
xsamp = x(1:os:end);


plot(tsamp,xsamp,'*r')
hold on

%%%
% reconstrução do sinal
y = zeros(1,length(x));
for n = 1:min(nsincs, length(tsamp))
    aux = xsamp(n)*sinc((t-tsamp(n))/Ta);
    y = y + aux;
end


plot(t, y, 'r')
title('Amostragem à taxa de Nyquist')
hold off

%%%%
% amostragem com taxa inferior à frequência de Nyquist
figure(2)
plot(t,x)
hold on

B = 4*f0;
BNyq = 2*B;
Ba = BNyq*0.8;
Ta = 1/Ba; 

tsamp = -10:Ta:10;
os = Ta/T0;
xsamp = x(1:os:end);


plot(tsamp,xsamp,'*r')
hold on

%%%
% reconstrução do sinal
y = zeros(1,length(x));
for n = 1:min(nsincs, length(tsamp))
    aux = xsamp(n)*sinc((t-tsamp(n))/Ta);
    y = y + aux;
end


plot(t, y, 'r')
title('Amostragem a 80% da taxa de Nyquist')
hold off

%%%%
% amostragem com taxa superior à frequência de Nyquist
figure(3)
plot(t,x)
hold on

B = 4*f0;
BNyq = 2*B;
Ba = BNyq*1.25;
Ta = 1/Ba; 

tsamp = -10:Ta:10;
os = Ta/T0;
xsamp = x(1:os:end);


plot(tsamp,xsamp,'*r')
hold on

%%%
% reconstrução do sinal
y = zeros(1,length(x));
for n = 1:min(nsincs, length(tsamp))
    aux = xsamp(n)*sinc((t-tsamp(n))/Ta);
    y = y + aux;
end


plot(t, y, 'r')
title('Amostragem a 125% da taxa de Nyquist')
hold off




figure(4)
Xf = fft(x);
f = 0:length(Xf)-1;
f = f * 1/20;

%%%%
% filtro passa baixa na frequência
Bfil = 0.4;
filtro = (f<Bfil) | (max(f)- f < Bfil);
Zf = Xf .* filtro;

plot(f, abs(Xf),'--', f, filtro*length(Xf),'-', f, abs(Zf), '-');
xlabel('f(Hz)')
legend('sinal de entrada', 'filtro', 'sinal filtrado')

z = ifft(Zf);

%%%%
% filtro passa baixa no tempo
h_ideal = 2*Bfil*T0*sinc(2*Bfil*t);
zfil_t = conv(x, h_ideal);
zfil_t = zfil_t(10001:end-10000);

figure(5)
plot(t, x, t, real(z),t,zfil_t)
legend('sinal de entrada', 'sinal filtrado na frequência', 'sinal filtrado no tempo')






