% 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;