% Practica 1 TDS 2004-2005 - UAB
% Teorema del muestreo
%
% 2- Variacion de la frecuencia de muestreo
% Sea x(n) una senal muestreada
%
% Carguemos un ejemplo
[x Fs bits] = wavread('beyonce.wav');
N=size(x,1);
wavplay(x, Fs);
pause;
plot(x);
pause;
% Su TFSD es:
X=fft(x);
f=linspace(0,2*pi,N)';
plot(f,abs(X));
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
pause;
% Realizar un diezmado en factor D consiste en reducir la frecuencia de
% muestreo en un factor entero D. Para ello, en principio es suficiente
% con tomar 1 de cada D muestras y descartar el resto.
D=4;
xdownsampled=zeros(N/D,1);
for i=1:N/D
xdownsampled(i)=x((i-1)*D+1);
end
Fsdownsampled=Fs/D;
stem(x(10001:10050));
ylim([-0.2 0.2]);
hold on;
pause;
stem(xdownsampled(10000/D+1:10000/D+50),'g');
pause;
hold off;
% Como vemos el espectro es el doble de ancho. Esto significa que la
% frecuencia de muestreo es la mitad.
Xdownsampled=fft(xdownsampled);
fdownsampled=linspace(0,2*pi,N/D)';
subplot(2,1,1);
plot(f,abs(X));
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal original');
subplot(2,1,2);
plot(fdownsampled,abs(Xdownsampled),'g');
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal diezmada');
pause;
clf;
% Si hacemos zoom en la zona de pi radianes, observamos que ha habido
% aliasing. Por eso, cuando se realiza un diezmado primero se pasa la
% senal por un filtro anti-aliasing (paso-bajos),
PB=[ones(N/2/D,1);zeros(N*(1-1/D),1);ones(N/2/D,1)];
Xaa=PB.*X;
xaa=real(ifft(Xaa));
plot(f,abs(X));
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
hold on;
pause;
plot(f,PB/max(PB)*max(abs(X)),'y')
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
pause;
plot(f,abs(Xaa),'r');
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
pause;
hold off;
% y luego se reduce el numero de muestras. Podemos comparar los
% espectros con el resultado anterior (downsample sin filtrado
% anti-aliasing) y oir las diferencias (si tienes buen oido).
% Nota: fijemonos que la amplitud del espectro ha disminuido en un
% factor D. Esto es logico, ya que para conservar la energia de la
% senal, si el ancho de banda aumenta en factor D, la amplitud tiene
% que disminuir por el mismo factor.
xdiezmado=zeros(N/D,1);
for i=1:N/D
xdiezmado(i)=xaa((i-1)*D+1);
end
Xdiezmado = fft(xdiezmado);
subplot(2,1,1);
plot(fdownsampled,abs(Xdownsampled));
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal diezmada sin filtro anti-aliasing');
subplot(2,1,2);
plot(fdownsampled,abs(Xdiezmado),'g');
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal diezmada con filtro anti-aliasing');
pause;
wavplay(xdownsampled, Fsdownsampled);
pause;
wavplay(xdiezmado, Fsdownsampled);
pause;
clf;
% Realizar un interpolado de factor I consiste en aumentar la frecuencia de
% muestreo en un factor entero I. Para ello, lo primero que haremos sera
% intercalar ceros. Usaremos la senal anterior diezmada y haremos una
% interpolacion por el mismo factor que hemos diezmado previamente.
I=D;
xintercalado=zeros(N,1);
for i=1:N/D
xintercalado((i-1)*I+1)=xdiezmado(i);
end
stem(xdiezmado(10001:10050));
hold on;
pause;
stem(xintercalado(10000*I+1:10000*I+50),'g');
pause;
hold off;
% El espectro de esta senal es el espectro anterior comprimido y repetido
% periodicamente:
Xintercalado=fft(xintercalado);
plot(f,abs(Xintercalado));
subplot(2,1,1);
plot(fdownsampled,abs(Xdiezmado),'g');
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal antes de intercalar ceros');
subplot(2,1,2);
plot(f,abs(Xintercalado));
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal despues de intercalar ceros');
pause;
clf;
% Para obtener la senal interpolada, a continuacion aplicaremos sobre la
% senal en la que hemos intercalado los ceros, un FPB de frecuencia de corte
% igual al ancho de banda de la senal interpolada dividido por I. En
% nuestro caso, resulta ser el mismo filtro utilizado en el proceso de
% diezmado, pero con ganancia I para mantener la energia de la senal.
Xinterpolado=I*PB.*Xintercalado;
xinterpolado=real(ifft(Xinterpolado));
plot(f,abs(Xintercalado));
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
hold on;
pause;
plot(f,PB/max(PB)*max(abs(Xintercalado)),'y')
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
pause;
plot(f,abs(Xinterpolado),'r');
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
pause;
hold off;
% Este procesado frecuencial tiene el efecto temporal de una interpolacion
% de la senal mediante una funcion sinc. Es un proceso similar al de
% conversion D/A.
stem(xdiezmado(10001:10050));
hold on;
pause;
stem(xintercalado(10000*I+1:10000*I+50),'g');
pause;
stem(xinterpolado(10000*I+1:10000*I+50),'r');
pause;
hold off;
% Vemos que el proceso de diezmado y el de interpolacion son simetricos,
% aunque solo son completamente reversibles en el caso de que la senal
% procesada sea de banda limitada y de frecuencia suficientemente baja para
% no ser afectada por el filtro anti-aliasing del proceso de diezmado.
% Este no es el caso de nuestro ejemplo, de ahi la diferencia que vemos
% entre la senal original y la final.
subplot(3,1,1);
plot(f,abs(X));
ylim([0 1500]);
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal original');
subplot(3,1,2);
plot(fdownsampled,abs(Xdiezmado),'g');
ylim([0 750]);
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal diezmada');
subplot(3,1,3);
plot(f,abs(Xinterpolado));
ylim([0 1500]);
set(gca,'XTick',0:pi/2:2*pi)
set(gca,'XTickLabel',{'0','pi/2','pi','3*pi/2','2*pi'})
title('Senal diezmada e interpolada');
pause;
wavplay(x, Fs);
pause;
wavplay(xinterpolado, Fs);
pause;
clf;
