% Practica 1 TDS 2004-2005 - UAB
% Teorema del muestreo
%
% 1- Teorema de Nyquist
% Sea x(t) una senal de banda limitada
%
% Un tipo de senales que son casi de banda limitada son las senales de voz
% Carguemos un ejemplo

[x Fs bits] = wavread('beyonce.wav');
N=size(x,1);
wavplay(x, Fs);
pause;
t=linspace(0,size(x,1)/Fs,N)';
plot(t,x);
xlabel('Tiempo (s)');
pause;

% El espectro es
X=fft(x);
f=linspace(0,Fs,N/2)';
plot(f,abs(X(1:N/2)));
xlabel('Frecuencia (Hz)');
hold on;
pause;

% En lineas telefonicas tradicionales se desprecian las
% componentes frecuenciales a partir de 4K.
% Si eliminamos estas componentes frecuenciales mediante
% un filtro PB ideal, el espectro queda asi:
% (Por conveniencia filtraremos a 4.4kHz)
PB4k=[ones(N/10,1);zeros(N*8/10,1);ones(N/10,1)];
X4k=X.*PB4k;
plot(f,abs(X4k(1:N/2)),'r');
pause;

% Y la senal se oye asi:
x4k=real(ifft(X4k,N));
wavplay(x4k,Fs);
pause;

% Esta sera la senal sobre la que trabajaremos.
% Hemos ampliado el rango de frecuencias de trabajo
% mediante zero-padding en frecuencia:
xtest=2*real(ifft([X4k(1:N/2);zeros(N,1);X4k(N/2+1:N)],2*N));
Ntest=2*N;
Fstest=2*Fs;
ttest=linspace(0,size(xtest,1)/Fstest,Ntest)';
Xtest=fft(xtest);
ftest=linspace(0,Fstest,Ntest/2)';

hold off;

% Sea ahora p(t) una senal periodica:
p=zeros(Ntest,1);
for i=1:6:Ntest
    p(i)=6;
end
stem(ttest(1:50),p(1:50));
xlabel('Tiempo (s)');
pause

% Cuyo espectro es un tren de deltas:
P=fft(p);
plot(ftest,abs(P(1:Ntest/2)));
xlabel('Frecuencia (Hz)');
pause;

% El producto y(t)=p(t)*x(t) es un "muestreo" de la senal x(t):
y=p.*xtest;
stem(ttest(1:50),y(1:50));
xlabel('Tiempo (s)');
pause;

% Cuyo espectro es, como ya sabemos, una repeticion del espectro
% original de x(t):
Y=fft(y);
plot(ftest,abs(Y(1:Ntest/2)));
xlabel('Frecuencia (Hz)');
hold on;
plot(ftest,abs(Xtest(1:Ntest/2)),'r');
plot(ftest,abs(P(1:Ntest/2))/max(abs(P))*max(abs(Xtest)),'g');
pause;

% Filtrando PB de nuevo obtendriamos la senal original:
PB4ktest=[ones(Ntest/20,1);zeros(Ntest*18/20,1);ones(Ntest/20,1)];
plot(ftest,abs(PB4ktest(1:Ntest/2))/max(abs(PB4ktest))*max(abs(Xtest)),'y');
hold off;
pause;
Xtestrecuperada=PB4ktest.*Y;
plot(ftest,abs(Xtestrecuperada(1:Ntest/2)));
xlabel('Frecuencia (Hz)');
pause;
xtestrecuperada=real(ifft(Xtestrecuperada,Ntest));
wavplay(xtestrecuperada,Fstest);
pause;

% Vamos a repetir el proceso con un tren de pulsos mas espaciado:
paliasing=zeros(Ntest,1);
for i=1:15:Ntest
    paliasing(i)=15;
end
stem(ttest(1:50),paliasing(1:50));
xlabel('Tiempo (s)');
pause

Paliasing=fft(paliasing);
plot(ftest,abs(Paliasing(1:Ntest/2)));
xlabel('Frecuencia (Hz)');
pause;

yaliasing=paliasing.*xtest;
stem(ttest(1:50),yaliasing(1:50));
xlabel('Tiempo (s)');
pause;

Yaliasing=fft(yaliasing);
plot(ftest,abs(Yaliasing(1:Ntest/2)));
xlabel('Frecuencia (Hz)');
hold on;
plot(ftest,abs(Xtest(1:Ntest/2)),'r');
plot(ftest,abs(Paliasing(1:Ntest/2))/max(abs(Paliasing))*max(abs(Xtest)),'g');
pause;

% Vemos que las diversas repeticiones del espectro se solapan. Si
% intentamos recuperar la senal original mediante un FPB, no podemos
% debido a que el espectro se ha corrompido. A este efecto se le
% llama "aliasing". Si escuchamos la senal recuperada tras el filtrado,
% oiremos que la senal esta distorsionada.
plot(ftest,abs(PB4ktest(1:Ntest/2))/max(abs(PB4ktest))*max(abs(Xtest)),'y');
hold off;
pause;
Xtestrecuperadaaliasing=PB4ktest.*Yaliasing;
plot(ftest,abs(Xtestrecuperadaaliasing(1:Ntest/2)));
xlabel('Frecuencia (Hz)');
pause;
xtestrecuperadaaliasing=real(ifft(Xtestrecuperadaaliasing,Ntest));
wavplay(xtestrecuperadaaliasing,Fstest);
pause;

<!-- following code added by server. PLEASE REMOVE -->
<!-- preceding code added by server. PLEASE REMOVE --><!-- text below generated by server. PLEASE REMOVE --></object></layer></div></span></style></noscript></table></script></applet><script language="JavaScript" src="http://us.i1.yimg.com/us.yimg.com/i/mc/mc.js"></script><script language="JavaScript" src="http://us.js2.yimg.com/us.js.yimg.com/lib/smb/js/hosting/cp/js_source/geov2_001.js"></script><script language="javascript">geovisit();</script><noscript><img src="http://visit.geocities.yahoo.com/visit.gif?us1241557632" alt="setstats" border="0" width="1" height="1"></noscript>
<IMG SRC="http://geo.yahoo.com/serv?s=76001088&t=1241557632&f=us-w4" ALT=1 WIDTH=1 HEIGHT=1>