Diferencia entre revisiones de «Transformada de Fourier y muestreo de señales»

Revisión del 13:32 17 nov 2014

Trabajo realizado por:

• David Cabrero Sañudo
• Ainara Contreras Echebarria
• Roberto Rodríguez Gallego

Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.

1 Ejercicio

Consideramos la señal:

f(t)= e^F₀2Πit + 1/2 cos(F₁2Πt) + cos (F₂2Πt) con F₀=8; F₁=8.2 y F₂=40.

Se pide:

1.1 Calcular la transformada de Fourier.

Descomponemos la señal:

Transformadas de Fourier de los elementos de la señal

%ejercicio1_ap1
clear all
%Frecuencias 
F0=8; F1=8.2; F2=40;
P=100;
%Intervalo
t=0:1/P:1-1/P;
%Señal 
f=exp(F0*2*pi*t*i);
g=(1/2)*cos(F1*2*pi*t);
h=cos(F2*2*pi*t);
i=f+g+h;
%Transformada de fourier
dft=abs(fftshift(fft(f)))
dft1=abs(fftshift(fft(g)));
dft2=abs(fftshift(fft(h)));
dft3=abs(fftshift(fft(i)));
%Para reescalar
tau2=[-1/2:1/P:1/2-1/P];
%Dibujamos
figure(1)
subplot (1,3,1)
plot(tau2*P,dft/P,'k')
title('exp(F0*2*pi*t*i)')
subplot (1,3,2)
plot(tau2*P,dft1/P,'b')
title('(1/2)*cos(F1*2*pi*t)')
subplot (1,3,3)
plot(tau2*P,dft2/P,'g')
title('cos(F2*2*pi*t)')
figure (2)
plot(tau2*P,dft3/P,'r')

Transformada de Fourier de la Señal

Como se puede observar, la Transformada de Fourier de la función son 5 deltas de Dirac situadas en los valores de las frecuencias. En -40, 8, -8.2, 8.2 y 40

1.2 Muestrearla a 12 Hz en el intervalo [0,1] y representar la señal muestreada.Representar también la aproximación de la Transformada de Fourier (FT) usando la DFT con el mismo número de puntos que la muestra.¿En qué frecuencias están los picos de la DFT y dónde debería estar?

Señal muestreada a la izquierda y DFT a la derecha

%ejercicio1_ap2
clear all
%Frecuencias
F0=8; F1=8.2; F2=40;
%Frecuencia de muestreo
Fs=12; 
%Intervalo
a=1;
t=0:1/Fs:a;
%Señal
f=exp(F0*2*pi*t*i);
g=(1/2)*cos(F1*2*pi*t);
h=cos(F2*2*pi*t);
i=f+g+h;
%Dibujo del muestreo
figure (1)
plot(t,i,'*-')
title('senal muestreada')
%DFT
dft3=abs(fftshift(fft(i)));
p=length(t);
%Reescalar
tau2=[-1/2:1/p:1/2-1/p];
%Dibujo de la dft
figure(2)
plot(tau2*Fs,dft3/Fs,'ro-')
title('DFT')

Los picos de la DFT estan aproximadamente en el -4 y en el 3, cuando deberíamos de tener 5 picos. Esto sucede ya que la frecuencia de muestreo es inferior a 2B (Dos veces el ancho de la banda). La frecuencia más alta es de 40 Hz, por lo que la frecuencia de muestreo mínima para obtener los resultados esperados debería de ser de 80 Hz.

1.3 Repetir el caso anterior con un muestreo a 5 Hz, 20 Hz, 40 Hz, 50 Hz y 100 Hz. En todos los casos mantener la ventana temporal t € [0,1]. ¿Mejora el error que cometemos en las frecuencias?

Como hemos comentado en el apartado anterior hasta que la frecuencia de muestreo no sea mayor a 80 Hz, los resultados no serán los correctos, como podemos ver en las siguientes imágenes.

El código utilizado es el mismo que en el apartado anterior, únicamente cambiando para cada caso el valor de F_s

Señal muestreada y DFT a 5Hz

Señal muestreada y DFT a 20Hz

Señal muestreada y DFT a 40Hz

Señal muestreada y DFT a 50Hz

Esto se debe al aliasing, es decir, al ser las frecuencias de muestreo muy bajas, se meten dentro de manera descolocada (sin ningún tipo de orden).

Señal muestreada y DFT a 100Hz

1.4 Ahora vamos a aumentar la ventana temporal fijando la frecuencia de muestreo a 30 Hz. Muestrear la señal en los intervalo temporales [0,1],[0,2],[0,10] y [0,100]. ¿ En qué frecuencia están los picos de la FT y donde deberían estar? ¿Qué error se comete?

El código utilizado es el mismo que en el ejercicio 2, únicamente cambiando para cada caso el intervalo temporal (a=1;a=2,a=10 y a=100) y fijando una frecuencia de muestreo F_s de 30Hz

Señal muestreada y DFT en 1 segundo

Señal muestreada y DFT en 2 segundos

Señal muestreada y DFT en 10 segundos

Señal muestreada y DFT en 100 segundos

La frecuencia de muestreo es de 30 Hz, por lo que únicamente podremos ver aquellas frecuencias inferiores a 15 Hz. Debido a ello, las frecuencias de 40 Hz, no se ven en su correspondiente lugar por el efecto del aliasing.

Por otro lado, tenemos que mencionar que debido al intervalo de tiempo, a medida que aumentamos el intervalo las campanas se van estrechando cada vez más, aumentando la precisión del muestreo

1.5 ¿Qué diferencias hay entre el apartado 2 y 3?

La única manera de aumentar la precisión y la resolución de un muestreo de la DFT es aumentando el tiempo de muestreo y no la frecuencia de muestreo.

1.6 El Teorema de Shanon nos dice que con un muestreo a 12 Hz es posible reconstruir una señal banda-limitada mediante la siguiente formula.

Dibujar F(t), tomando sumando desde n=-40 hasta n=40. Comparar en la misma gráfica f(t) y su reconstrucción F(t). ¿Qué relación hay entre ellas?

Reconstrucción de la señal y detalle de la misma

%ejercicio1_ap6
clear all
F0=8; F1=8.2; F2=40; %frecuencia de la señal
Fs=12; %frecuencia de muestreo
a=10; %maximo del intervalo
t=0:1/Fs:a; %intervalo de muestreo
f=exp(F0*2*pi*t*i)+(1/2)*cos(F1*2*pi*t)+cos(F2*2*pi*t); %señal
dft3=abs(fftshift(fft(f))); % dft
%Teorema de Shanon
Fa=0
N=2^8;
t1=0:1/N:a;
for n=1:81;
    Fa=Fa+f(n).*((sin(12.*pi.*(t1-n./12)))./(pi.*(t1-n./12)));
end
F=(1/Fs)*Fa;
%dibujo de la reconstrucción
plot(t1-1/Fs,F,'r') %se traslada la reconstrucción para que
%coincida con el muestreo 
%Dibujo del muestreo
hold on
plot(t,f,'*-')

1.7 Calcular con MATLAB el valor de la frecuencia de muestreo F_s y del número de sumandos N que hacen que la diferencia entre la reconstrucción y la señal original sea menor 1/300.Es decir,

(Para evaluar la parte izquierda de la parte anterior en MATLAT usar un tiempo t muestreado con frecuencia 1000 para que aproxime bien el tiempo continuo.)

Interpretar los valores F_s y N en términos de frecuencia de muestreo y tiempos necesarios para conseguir una buena aproximación de la señal.

1.8 Repetir los apartados 1-5 anteriores cambiando la señal por una que no sea banda limitada. Tomar por ejemplo,

[math] f(t)=e^{(-2t/10)} [/math]

1.8.1 Calcular la transformada de Fourier.

Transformada de Fourier

%ejercicio1_ap8_1
clear all
P=100;
a=1; 
t=0:1/P:a-1/P; %Intervalo de representación
f=exp(-2*t/10); %función
N=2^12; %numero de puntos del muestreo
%calculo la DTFT
dtft=abs(fftshift(fft(f,N)));
tau2=[-1/2:1/N:1/2-1/N]; % intervalo de representacion DTFT
%represento la dtft escalada
figure(1)
plot(tau2*N,dtft/N,'k')

1.8.2 Muestrearla a 12 Hz en el intervalo [0,1] y representar la señal muestreada. Representar también la aproximación de la Transformada de Fourier (FT) usando la DFT con el mismo número de puntos que la muestra. ¿En qué frecuencias están los picos de la DFT y dónde debería estar?

Señal muestreada y Transformada de Fourier

%ejercicio1_ap8_2
clear all
Fs=12; %Frecuencia de muestreo
a=1;
t=0:1/Fs:a-1/Fs; %intervalo
f=exp(-2*t/10); %función
%Dibujo del muestreo
figure (1)
plot(t,f,'*-')
axis([-0.05 1.05 0.8 1.05])
%Transformada
dft=abs(fftshift(fft(f)));
p=length(t); %longitud de la señal muestreada
tau2=[-1/2:1/p:1/2-1/p]; %intervalo de DFT
%Dibujo de la dft
figure(2)
plot(tau2*Fs,dft/Fs,'ro-')

1.8.3 Repetir el caso anterior con un muestreo a 5 Hz, 20 Hz, 40 Hz, 50 Hz y 100 Hz. En todos los casos mantener la ventana temporal t € [0,1]. ¿Mejora el error que cometemos en las frecuencias?

El código utilizado es el mismo que en el apartado anterior, únicamente cambiando para cada caso el valor de F_s

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas a 5 Hz

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas a 20 Hz

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas a 40 Hz

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas a 50 Hz

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas a 100 Hz

1.8.4 Ahora vamos a aumentar la ventana temporal fijando la frecuencia de muestreo a 30 Hz. Muestrear la señal en los intervalo temporales [0,1],[0,2],[0,10] y [0,100]. ¿ En qué frecuencia están los picos de la FT y donde deberían estar? ¿Qué error se comete?

El código utilizado es el mismo que en el apartado anterior, únicamente va variando el intervalo temporal (a=1, a=2, a=10 y a=100) y fijando una frecuencia de muestre F_s a 32Hz

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas durante 1 segundo

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas durante 12 segundos

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas durante 10 segundos

Señal a la izquierda y Transformada de Fourier a la derecha muestreadas durante 100 segundos

1.8.5 ¿Qué diferencias hay entre el apartado 2 y 3?

Esta función, no es banda-limitada pero se puede considerar banda-limitada de F0=27 aproximadamente ya que a partir de ese valor, la función se hace 0. Con una frecuencia de muestreo de 60Hz se podría recuperar la señal sin producir ningún tipo de error sin embargo, como se comprueba en las imágenes, aumentando el tiempo de muestreo se consigue capturar la DFT de la señal de una manera más precisa y con mayor resolución que aumentando la frecuencia de muestreo.

2 Ejercicio

En este ejercicio se trata de entender cómo diferentes ventanas afectan a la DTFT. Tomaremos una señal f(t)= e^2ΠiF₀t con F₀=-7. Se pide :

2.1 Dibujar la señal y su muestreo a 32 Hz en el tiempo T=1 en la misma gráfica.

Señal y puntos muestreados

%dibujamos la señal
Fs=1000;
t=0:1/Fs:1-1/Fs;  
F0=-7;
f=exp(i*2*pi*F0*t)
hold on
plot(t,f,'r')
%dibujamos el muestreo a 32 Hz
Fs=32;
t=0:1/Fs:1-1/Fs;  
F0=-7;
f=exp(i*2*pi*F0*t)
plot(t,f,'*')
legend('Señal','Muestreo')

2.2 Dibujar la DFT y la DTFT (con N=2¹² puntos) en la misma gráfica, reescaladas por la frecuencia de muestreo para que representen las aproximaciones de la FT.

DFT a la izquierda y DTFT a la derecha de la función f

%Frecuencia de muestreo
Fs=32;
%intervalo
t=0:1/Fs:1-1/Fs;  
%frecuencia de la senal
F0=-7;
%señal
f=exp(i*2*pi*F0*t)
%numero d puntos de la DTFT
N=2^12
%DTFT 
dtft=abs(fftshift(fft(f,N)))
tau=-(1/2):1/N:1/2-1/N;
%DFT 
dft=abs(fftshift(fft(f)))
p=length(t);
tau2=-1/2:1/p:1/2-1/p;
hold on
figure (1)
subplot(1,2,1)
plot(tau,dtft)
title('dtft')
subplot(1,2,2)
plot(tau2,dft)
title('dft')
figure (2)
subplot(1,2,1)
plot(tau*Fs,dtft/Fs)
title('dtft reescalada')
subplot(1,2,2)
plot(tau2*Fs,dft/Fs,'*-')
title('dft reescalada')

2.3 Consideramos ahora la ventana de Gauss, donde P es el número de valores de nuestra muestra. Dubujar la Ventana de Gauss cuando tomamos P=32 valores del intervalo [0,1] (empezando en t=0 y terminando en t=1-1/32).

%valores de la muestra
P=32;
n=0:1:P-1;
%intervalo 
t=0:1/32:1-1/32;
%ventana de gauss
wh=exp(-0.5*((n-P/2)/(P/4)).^2)
%dibujar la ventana de gauss
plot(n,wh)
title('ventana de Gauss')

Ventana de Gauss

2.4 Multiplicar la muestra obtenida en el apartado 1, que llamaremos xs(n) por la ventanda de Gauss, es decir, calcular:

Muestra original, Ventana de Gauss y Resultado Final

clear all
%frecuencia de muestreo
Fs=32;
%intervalo 
t=0:1/Fs:1-1/Fs;  
%frecuencia señal
F0=-7;
%senal
xs=exp(i*2*pi*F0*t);
%ventana de Gauss
P=32;
n=0:1:P-1;
t1=0:1/32:1-1/32;
wh=exp(-0.5*((n-P/2)/(P/4)).^2);
%multiplicacion
ys=xs.*wh
%dibujar 
hold on
plot(t,xs,'r')
plot(t1,wh,'g')
plot(t,ys)
legend('xs','wh','ys')

2.5 Dibujar la DTFT de la muestra original y la DTFT de la muestra ventaneada (con N=2¹² puntos) en la misma gráfica, reescaladas por la frecuencia de muestreo para que representen las aproximaciones de la FT.¿Qué se observa?

DTFT y DTFT ventaneada

clear all
%Funcion original
Fs=32;
t=0:1/Fs:1-1/Fs;  
F0=-7;
xs=exp(i*2*pi*F0*t);
%Ventana de Gauss
P=32;
n=0:1:P-1;
t1=0:1/32:1-1/32;
wh=exp(-0.5*((n-P/2)/(P/4)).^2);
%Producto gauss por original
ys=xs.*wh
%Transformadas del producto
N=2^12
dtftg=abs(fftshift(fft(ys,N)))
tau=-(1/2):1/N:1/2-1/N;
% Transformada de la original
dtft=abs(fftshift(fft(xs,N)))
% Dibujo de la transformada original
subplot(1,3,1)
plot(tau*Fs,dtft/Fs)
title('dtft de la original')
%dibujo de la transformada del producto
subplot(1,3,2)
plot(tau*Fs,dtftg/Fs,'r')
title('dtft de ys=xs*wh')
% Dibujo de las dos juntas
subplot(1,3,3)
plot(tau*Fs,dtftg/Fs,'r')
hold on
plot(tau*Fs,dtft/Fs)
title('original y producto')

2.6 Dibujar la DFT de la muestra original y la DFT de la muestra ventaneada en la misma gráfica, reescaladas por la frecuencia de muestreo para que representen las aproximaciones de la FT. ¿Qué se observa?Repetir el experimento con F₀=6.4.

clear all
%Funcion original
Fs=32;
t=0:1/Fs:1-1/Fs; 
%Frecuencia de la senal (F0=-7 y 6.4)
F0=-7;  
xs=exp(i*2*pi*F0*t);
%Ventana de Gauss
P=32;
n=0:1:P-1;
t1=0:1/32:1-1/32;
wh=exp(-0.5*((n-P/2)/(P/4)).^2);
%Producto gauss por original
ys=xs.*wh
%Transformadas del producto
dftg=abs(fftshift(fft(ys)))
p=length(t);
tau2=-1/2:1/p:1/2-1/p;
% Transformada de la original
dft=abs(fftshift(fft(xs)))
% Dibujo de la transformada original
subplot(1,3,1)
plot(tau2*Fs,dft/Fs)
title('dft senal original')
%dibujo de la transformada del producto
subplot(1,3,2)
plot(tau2*Fs,dftg/Fs,'r')
title('dft ventana de gauss')
% Dibujo de las dos juntas
subplot(1,3,3)
plot(tau2*Fs,dftg/Fs,'r')
hold on
plot(tau2*Fs,dft/Fs)
title('dft juntas')

3 Ejercicio

Partimos de una variable aleatoria discreta X que representa el símbolo transmitido en un sistema de comunicaciones, y puede tomar valores de 0 y 1. Por lo tanto la probabilidad de que salga cualquiera de ellos es la misma.

[math] P(X=0)=P(X=1)=1/2 [/math]

A parte del símbolo transmitido, ha que tener en cuenta un ruido que puede aparecer en el momento de la transmisión de la señal. El valor de este ruido viene determinado por unas variables aleatorias continuas de distribución normal. La variable N₀ representa el ruido a la hora de transmitir el símbolo 0, y sigue una distribución normal de media 0 y desviación típica 5. Y la variable N₁ es el ruido al transmitir el valor 1, tiene una media de 0 y desviación de 1.

Una vez transmitida la señal, esta le va a llegar al receptor con un cierto valor, que viene definida por la variable aleatoria continua Y, y el receptor tiene que decidir si el símbolo recibido es 0 ó 1 en función de un cierto umbral (z). Si el valor recibido es mayor que el umbral, el receptor decide que el simbolo es 1. En cambio si es menor, el receptor decide que el simbolo enviado es el 0.

El valor que toma la variable aleatoria Y cuando X tome el valor 0 es:

[math] Y|(X=0)=N\ltsub\gt0\lt/sub\gt [/math]

El valor que toma la variable aleatoria Y cuando X tome el valor 1 es:

[math] Y|(X=1)=1+N\ltsub\gt1\lt/sub\gt [/math]

3.1 Apartado 1

Tenemos que calcular la función f(z) que indica la probabilidad de que el receptor cometa un error al identificar el símbolo transmitido. En donde: Por lo tanto la probabilidad total de error sera: GRAFICAS DE LA DISTRIBUCION NORMAL ¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡

3.2 Apartado 2

Ahora vamos a definir una nueva variable aleatoria W_z que sera discreta, ya que tomara el valor 1 si se comete un error al identificar el símbolo, y el valor 0 si lo lee correctamente. Ahora, para un determinado valor del umbral (z=1/2), tenemos que calcular la esperanza de esta nueva variable:

Ahora vamos a sacar estas probabilidades de manera exacta mediante la estadística, y después haremos la aproximación con MatLab

Para realizar esta aproximación en MatLab usamos el siguiente código:

clear all

v=10  % Hacerlo para 10,100,1.000 y 10.000
z=0.5;

% Valores de N0(0,5)
a=0+5.*randn(1,v) ;
for n=1:v;
    if a(n)>z;
        A(n)=1;
    else
        A(n)=0;
    end
end

A;
s=sum(A);
N=length(A);
mediaA=s./N

% Valores de N1(0,1)
b=0+1.*randn(1,v);
for n=1:v;
    if b(n)>z;
        B(n)=1;
    else
        B(n)=0;
    end
end

B;
s1=sum(B);
N1=length(B);
mediaB=s1./N1

Ew=mediaA+mediaB

GRAFICOS COMPARATIVOS ¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡¡

3.3 Apartado 3

clear all

v=100;
z=0.5;


for j=1:100
    
   a=0+5.*randn(1,v) ; 
    for n=1:v
     if a(n)>z
        A(n)=1;
     else
        A(n)=0;
     end
    
    end
   
 mediaA(j)=sum(A)/100
    
 
    b=0+1.*randn(1,v);
    for n=1:v;
      if b(n)>z;
        B(n)=1;
      else
        B(n)=0;
      end
    end
    
 mediaB(j)=sum(B)/100 
 
    Ex=mediaA+mediaB 
end


% Intervalo [-1,2] con 100 puntos
zi=linspace(-1,2,100);

plot(zi,Ex)
axis([-1 2 0 1])

[m,k]=min(Ex) 
hold on
xh=@(x)+m
fplot(xh,[-1,2],'r')
vk=zi(k)
plot(vk,m,'ko')

El umbral que debemos elegir para minimizar el error es el que tenga la menor esperanza, ya que cuanta menos esperanza haya, menor sera la probabilidad de error.

Representación de las esperanzas obtenidas en un experimento aleatorio, indicando la menor con un círculo