Captura de vídeo en GUI Matlab

Con este programa se pretende mostrar la forma de capturar vídeo e imágenes en una interfaz gráfica.
El campo tag del axes que presenta el vídeo es axes1 y del axes que contiene la imagen capturada es axes2. El código en la sección de inicialización del programa es el siguiente:
  1. function video_gui_OpeningFcn(hObject, eventdata, handles, varargin)
  2. movegui(hObject,'center')
  3. set(handles.axes1,'XTick',[ ],'YTick',[ ])
  4. set(handles.axes2,'XTick',[ ],'YTick',[ ])
  5. imaqreset
  6. handles.output = hObject;
  7. guidata(hObject, handles);
Con la herramienta Object browser podemos ver el campo tag de cada elemento de la GUI. (Ver Figura siguiente).
El botón de inicio de la GUI abre un programa que de forma automática detecta las cámaras de vídeo (Ver Figura siguiente).
En la función de apertura contiene el código del Listado siguiente.
  1. function sel_camera_OpeningFcn(hObject, eventdata, handles, varargin)
  2. imaqreset;
  3. set(handles.ok_b,'Enable','off')
  4. hw=imaqhwinfo('winvideo');
  5. handles.cam=hw;
  6. set(handles.lista_camaras,'String',{hw.DeviceInfo.DeviceName})
  7. handles.output = hObject;
  8. guidata(hObject, handles);
Con el pop up menu seleccionamos qué cámara usar. Dentro de la función asociada a este elemento (mostrada en el Listado siguiente) está la programación que extrae las principales características de la cámara seleccionada.
  1. function lista_camaras_Callback(hObject, eventdata, handles)
  2. pos=get(handles.lista_camaras,'Value');
  3. hw=handles.cam;
  4. id=hw.DeviceIDs{pos};
  5. set(handles.id_camara,'String',id)
  6. formatos=hw.DeviceInfo(pos).SupportedFormats;
  7. set(handles.formatos,'String',formatos)
  8. list_f = [formatos{1:end}];
  9. si=strfind(list_f,'RGB24_320x240');
  10. if isempty(si)
  11. es_web_ext=0;% Laptop: YUY2
  12. else
  13. es_web_ext=1;% External: RGB
  14. end
  15. handles.es_web_ext=es_web_ext;
  16. handles.id=id;
  17. guidata(hObject, handles)
  18. set(handles.ok_b,'Enable','on')
Al presionar el botón Ok se exportan las características de la cámara seleccionada. El Listado siguiente contiene el código de la función del botón. Note que se usan variables globales para exportar esta información a la primera GUI.
  1. function ok_b_Callback(hObject, eventdata, handles)
  2. global id es_web_ext
  3. es_web_ext = handles.es_web_ext;
  4. id = handles.id;
  5. close (handles.camara)
Regresando a la programación del la primera GUI, una vez seleccionada la fuente de vídeo, el botón de inicio muestra la secuencia de imágenes en el primer axes de la GUI (Ver Listado siguiente).
  1. function inicio_b_Callback(hObject, eventdata, handles)
  2. set(handles.inicio_b,'Enable','off')
  3. start(handles.vidobj);
  4. vidRes = get(handles.vidobj,'VideoResolution');
  5. nBands = get(handles.vidobj,'NumberOfBands');
  6. hImage = image(zeros(vidRes(2), vidRes(1), nBands),'Parent',handles.axes1);
  7. preview(handles.vidobj,hImage);
  8. guidata(hObject, handles);
Con el botón de captura, cuya programación se muestra en el Listado siguiente, obtenemos la imagen a partir del vídeo.
  1. function captura_Callback(hObject, eventdata, handles)
  2. try
  3. rgb = getsnapshot(handles.vidobj);
  4. if handles.es_web_ext==0
  5. rgb = ycbcr2rgb(rgb);
  6. end
  7. image(rgb,'Parent',handles.axes2);
  8. axes(handles.axes2)
  9. axis off
  10. catch
  11. disp('No hay imagen para mostrar')
  12. end
Finalmente, el botón de parada detiene la proyección de vídeo y permitirá seleccionar una nueva fuente de vídeo. (Ver Listado siguiente.)
  1. function parar_b_Callback(hObject, eventdata, handles)
  2. set(handles.inicio_b,'Enable','on')
  3. stoppreview(handles.vidobj)

GUI de Matlab para encontrar el punto de corte de dos funciones.

Dadas dos funciones, vamos a encontrar el punto de corte a través de una GUI de Matlab. La interfaz gráfica genera ambas funciones al presionar un botón. Luego, se calcula el punto de corte, se asigna su coordenada a un cuadro de texto estático. La siguiente gráfica muestra el entorno de esta GUI.
Programamos que el segundo botón esté desactivado cuando inicia a correr el programa. Esto lo hacemos en la parte correspondiente del archivo .m OpeningFcn.
  1. % --- Executes just before corte_graficas is made visible.
  2. function corte_graficas_OpeningFcn(hObject, eventdata, handles, varargin)
  3. set(handles.pushbutton2,'Enable','off')
  4. handles.output = hObject;
  5. guidata(hObject, handles);
El código que se ejecuta el presionar el primer botón se muestra a continuación. Se generan dos funciones y se colocan sus valores en handles para poder exportar estos valores a otra función del archivo .m.
  1. % --- Executes on button press in pushbutton1.
  2. function pushbutton1_Callback(hObject, eventdata, handles)
  3. x = 0:0.01:10;
  4. y1 = x .^2 + 2;
  5. y2 = x .^3 ;
  6. handles.x=x;
  7. handles.y1=y1;
  8. handles.y2=y2;
  9. handles.axes1;
  10. plot(x, y1, x, y2);
  11. ylim([0,10])
  12. set(handles.pushbutton2,'Enable','on')
  13. guidata(hObject, handles);
Finalmente, el botón etiquetado como Corte ejecuta el código que encuentra el punto de corte de las dos gráficas y presenta la coordenada en el cuadro estático de texto.
  1. % --- Executes on button press in pushbutton2.
  2. function pushbutton2_Callback(hObject, eventdata, handles)
  3. idx = find(handles.y1 - handles.y2 < eps, 1); %// Index of coordinate in array
  4. px = handles.x(idx);
  5. py = handles.y1(idx);
  6. set(handles.text1,'String',['x= ',num2str(px),' y= ',num2str(py)])
  7. handles.axes1;
  8. hold on
  9. plot(px, py, 'ro', 'MarkerSize', 18)
  10. hold off

Cosenoidal con varias fases en Matlab.

Graficar en Matlab la señal z\left( t \right) = 8\cos \left( {2\pi t + \theta } \right) para valores de \theta  = \left[ { - {\pi  \over 2}, - \pi , - {{3\pi } \over 2}} \right] en el intervalo 0 \le t \le 8 .
  1. t=-1:1/100:1;
  2. cc=['r','g','b'];
  3. w0=2*pi;
  4. desv_tmp=zeros(1,3);
  5. for k=1:3
  6. theta=[-pi/2,-pi,-1.5*pi];
  7. desv_tmp(k)=theta(k)/w0;
  8. z=8*cos(2*pi*t + theta(k));
  9. plot(t,z,'LineWidth',2,'Color',cc(k))
  10. hold on
  11. grid on
  12. end
  13. hold off

Gráfica del seno y coseno con la identidad de Euler en Matlab.

Graficar en Matlab las siguientes señales:
x\left( t \right) = {e^{j4\pi t}} + {e^{ - j4\pi t}}
y\left( t \right) = {{\left( {{e^{j4\pi t}} - {e^{ - j4\pi t}}} \right)} \over j}
  1. t=-1:1/100:1;
  2. s1=2*cos(4*pi*t);
  3. s2=2*sin(4*pi*t);
  4. subplot(2,1,1)
  5. plot(t,s1,'LineWidth',2),grid on
  6. subplot(2,1,2)
  7. plot(t,s2,'LineWidth',2),grid on

Errata del libro Señales y Sistemas de Oppenheim.

  • Ecuación 1.39 es exp(-j0.5t) (página 20).
  • Ecuación 1.123 es x[n].
  • Ecuación 2.14 (página 86).
  • Ecuación 2.21 (página 89).
  • Ecuación 2.30 el h lleva sombrero (página 96).
  • Ecuación 3.61 falta 't' en el exponente.
  • Ecuación 3.124. El argumentos de H es jkwo. No con la exponencial.
  • Problema 3.8. La señal x(t) es real e impar.
  • Página 305. Debe ser mayúscula la x.
  • Página 310, ejemplo 4.13: la última ecuación es -1, no -2.
  • Ecuación 4.45 (página 313) fata 1 sobre 2Pi.
  • Ec 4.49 exponente w0 falta.
  • Luego de la ecuación 4.83 el exponente de la segunda exponencial
  • Ejemplo 5.2. El menos de la seg exponencial
  • Ejemplo 5.7 Entre paréntesis va: y otros múltiplos impares de Pi.
  • Ecuación 5.47 (página 380) el límite inferior es -Inf.
  • Ejemplo 5.9. Para 0, “si n es impar”.
  • Figura 5.17. No es omega 0, es omega c.
  • Ecuación 9.90 es aR. 
  • Figura 9.24 es aR (ver etiqueta también), no R/a.
  • Página 687, luego de 9.94 “...poles at s=-1+-j3”.
  • Ecuación 7.18 es Xc mayúscula.
  • Página 916, ecuación A.53 no es A12, es A21.

Modulación AM.

  1. #AM modulation
  2. # Funcion SINC
  3. import numpy as np
  4. import pylab as plt
  5. from scipy.fftpack import fft
  6.  
  7. N = 500 # number of sample points
  8. dt = 1. / 1000 # sample spacing
  9.  
  10. npp=9
  11. nN=N*npp
  12. t = np.linspace(-N*dt, N*dt, nN)
  13.  
  14. y=np.sinc(150*t)*np.cos(2*np.pi*400*t)
  15. #y=np.abs(np.sinc(t))
  16. plt.figure(1)
  17. plt.plot(t,y, linewidth=2)
  18. plt.grid(True), plt.xlabel('Tiempo')
  19. plt.title(u'x(t)=sinc(150t)*cos(2 $\pi$ (400)t)')
  20. plt.xlim([-6./200,6./200])
  21.  
  22. # FFT
  23. yf = fft(y)
  24. tf = npp*np.linspace(-1./(4.*dt), 1./(4.*dt), nN)
  25. spectrum = 1./nN * np.concatenate([np.abs(yf[nN/2:nN]), np.abs(yf[0:nN/2])])
  26.  
  27. #figure1 = plt.figure(4, (10, 5))
  28. plt.figure(2)
  29. plt.plot(tf, spectrum, linewidth=2)
  30. plt.grid()
  31. plt.xlim([-700,700])
  32. plt.title(u'Espectro de magnitud |X(j$\omega$)|')
  33. plt.xlabel('Frecuencia [Hz]')
  34. plt.ylabel('Magnitud |X(j$\omega$)|')
  35.  

Transformada de Fourier de la función Sinc.

  1. import numpy as np
  2. import pylab as plt
  3. from scipy.fftpack import fft
  4.  
  5. N = 500 # number of sample points
  6. dt = 1. / 1000 # sample spacing
  7.  
  8. npp=3
  9. nN=N*npp
  10. t = np.linspace(-N*dt, N*dt, nN)
  11.  
  12. y=np.sinc(200*t)
  13. #y=np.abs(np.sinc(t))
  14. plt.figure(1)
  15. plt.plot(t,y)
  16. plt.grid(True), plt.xlabel('Tiempo')
  17. plt.title(u'Sinc(200t)')
  18. plt.xlim([-6./200,6./200])
  19.  
  20. # FFT
  21. yf = fft(y)
  22. tf = npp*np.linspace(-1./(4.*dt), 1./(4.*dt), nN)
  23. spectrum = 1./nN * np.concatenate([np.abs(yf[nN/2:nN]), np.abs(yf[0:nN/2])])
  24.  
  25. #figure1 = plt.figure(4, (10, 5))
  26. plt.figure(2)
  27. plt.plot(tf, spectrum, '-')
  28. plt.grid()
  29. plt.xlim([-300,300])
  30. plt.title(u'Espectro de magnitud |X(j$\omega$)|')
  31. plt.xlabel('Frecuencia [Hz]')
  32. plt.ylabel('Magnitud |X(j$\omega$)|')

Espectro unilateral de un suma de senoides en Python.

  1. import numpy as np
  2. #from scipy import signal
  3. from scipy.fftpack import fft
  4. import matplotlib.pyplot as plt
  5.  
  6. N = 1000 # number of sample points
  7. dt = 1. / 500 # sample spacing
  8.  
  9. frequency1 = 50.
  10. frequency2 = 150.
  11. nN=N*4
  12. t = np.linspace(0.0, N*dt, nN)
  13. s1 = 0.8*np.sin(2*np.pi * frequency1 * t)
  14. s2 = 0.4* np.sin(2*np.pi * frequency2 * t)
  15. y = s1 + s2
  16.  
  17. plt.figure(1)
  18. plt.plot(t, s1)
  19. plt.grid()
  20. plt.title('Senoide de '+str(frequency1)+' Hz')
  21. plt.ylabel('Amplitud')
  22. plt.xlabel('Tiempo [s]')
  23. plt.axis([0, 8/frequency2, -1.5, 1.5])
  24.  
  25. plt.figure(2)
  26. plt.plot(t, s2)
  27. plt.grid()
  28. plt.title('Senoide de '+str(frequency2)+' Hz')
  29. plt.ylabel('Amplitud')
  30. plt.xlabel('Tiempo [s]')
  31. plt.axis([0, 8/frequency2, -1.5, 1.5])
  32.  
  33.  
  34. plt.figure(3)
  35. plt.plot(t, s1+s2)
  36. plt.grid()
  37. plt.title('Suma de senoides de '+str(frequency1)+' Hz'+' y '+str(frequency2)+' Hz')
  38. plt.ylabel('Amplitud')
  39. plt.xlabel('Tiempo [s]')
  40. plt.axis([0, 8/frequency2, -1.5, 1.5])
  41. # FFT
  42. yf = fft(y)
  43. tf = np.linspace(.0, 1./(2.*dt), N/2)
  44. spectrum = 2./nN * np.abs(yf[0:N/2])
  45.  
  46. #figure1 = plt.figure(4, (10, 5))
  47. plt.figure(4)
  48. plt.plot(tf, spectrum, '-')
  49. plt.grid()
  50. plt.title(u'Espectro de magnitud |X(j$\omega$)|')
  51. plt.xlabel('Frequencia [Hz]')
  52. plt.ylabel('Magnitud |X(j$\omega$)|')

Espectro bilateral de un suma de senoides en Python.

  1. import numpy as np
  2. #from scipy import signal
  3. from scipy.fftpack import fft,fftshift
  4. import matplotlib.pyplot as plt
  5.  
  6. N = 1000 # number of sample points
  7. Fs=500.
  8. dt = 1. / Fs # sample spacing
  9.  
  10. frequency1 = 50.
  11. frequency2 = 150.
  12. npp=4
  13. nN=N*npp
  14. t = np.linspace(0.0, N*dt, nN)
  15. s1 = 0.8*np.sin(2*np.pi * frequency1 * t)
  16. s2 = 0.4* np.sin(2*np.pi * frequency2 * t)
  17. y = s1 + s2
  18.  
  19. plt.figure(1)
  20. plt.plot(t, s1)
  21. plt.grid()
  22. plt.title('Senoide de '+str(frequency1)+' Hz')
  23. plt.ylabel('Amplitud')
  24. plt.xlabel('Tiempo [s]')
  25. plt.axis([0, 8/frequency2, -1.5, 1.5])
  26.  
  27. plt.figure(2)
  28. plt.plot(t, s2)
  29. plt.grid()
  30. plt.title('Senoide de '+str(frequency2)+' Hz')
  31. plt.ylabel('Amplitud')
  32. plt.xlabel('Tiempo [s]')
  33. plt.axis([0, 8/frequency2, -1.5, 1.5])
  34.  
  35.  
  36. plt.figure(3)
  37. plt.plot(t, s1+s2)
  38. plt.grid()
  39. plt.title('Suma de senoides de '+str(frequency1)+' Hz'+' y '+str(frequency2)+' Hz')
  40. plt.ylabel('Amplitud')
  41. plt.xlabel('Tiempo [s]')
  42. plt.axis([0, 8/frequency2, -1.5, 1.5])
  43. # FFT
  44. yf = fft(y)
  45. tf = npp*np.linspace( -1./(2.*dt), 1./(2.*dt), nN)
  46. spectrum = 2./nN * np.abs(fftshift(yf))
  47.  
  48. #figure1 = plt.figure(4, (10, 5))
  49. plt.figure(4)
  50. plt.plot(tf, spectrum, '-')
  51. plt.grid()
  52. plt.title(u'Espectro de magnitud |X(j$\omega$)|')
  53. plt.xlabel('Frequencia [Hz]')
  54. plt.ylabel('Magnitud |X(j$\omega$)|')
  55. #plt.axis([-200, 200, 0., 1.5])

Coeficientes de la serie de Fourier conforme T tiende al infinito.

  1. import numpy as np
  2. import pylab as pl
  3. T1=1
  4. T=4*T1
  5. w0=2*np.pi/T
  6. k=np.arange(-40*w0,40*w0,w0)
  7. ak=(2.0*T1)*np.sinc(k)
  8. pl.subplot(4,1,1)
  9. pl.stem(k,ak*T)
  10. pl.hold(True)
  11. ke=np.arange(-40*w0,40*w0,w0/10)
  12. ake=(2.0*T1)*np.sinc(ke)
  13. pl.plot(ke,ake*T,'r--')
  14. pl.hold(False)
  15. pl.grid(True)
  16. pl.xlim([-15,15])
  17. T=8*T1
  18. w0=2*np.pi/T
  19. k=np.arange(-40*w0,40*w0,w0)
  20. ak=(2.0*T1)*np.sinc(k)
  21. pl.subplot(4,1,2)
  22. pl.stem(k,ak*T)
  23. pl.hold(True)
  24. ke=np.arange(-40*w0,40*w0,w0/10)
  25. ake=(2.0*T1)*np.sinc(ke)
  26. pl.plot(ke,ake*T,'r--')
  27. pl.hold(False)
  28. pl.grid(True)
  29. pl.xlim([-15,15])
  30. T=16*T1
  31. w0=2*np.pi/T
  32. k=np.arange(-40*w0,40*w0,w0)
  33. ak=(2.0*T1)*np.sinc(k)
  34. pl.subplot(4,1,3)
  35. pl.stem(k,ak*T)
  36. pl.hold(True)
  37. ke=np.arange(-40*w0,40*w0,w0/10)
  38. ake=(2.0*T1)*np.sinc(ke)
  39. pl.plot(ke,ake*T,'r--')
  40. pl.hold(False)
  41. pl.grid(True)
  42. pl.xlim([-15,15])
  43. T=32*T1
  44. w0=2*np.pi/T
  45. k=np.arange(-80*w0,80*w0,w0)
  46. ak=(2.0*T1)*np.sinc(k)
  47. pl.subplot(4,1,4)
  48. pl.stem(k,ak*T)
  49. pl.hold(True)
  50. ke=np.arange(-40*w0,40*w0,w0/10)
  51. ake=(2.0*T1)*np.sinc(ke)
  52. pl.plot(ke,ake*T,'r--')
  53. pl.hold(False)
  54. pl.grid(True)
  55. pl.xlim([-15,15])

Función sinc en Python.

  1. # Funcion SINC
  2. import numpy as np
  3. import pylab as pl
  4.  
  5. t=np.linspace(-10,10,1000)
  6. s=np.sinc(t)
  7. s_mod=np.abs(np.sinc(t))
  8. pl.figure(1)
  9. pl.plot(t,s)
  10. pl.grid(True), pl.xlabel('Tiempo')
  11.  
  12. pl.figure(2)
  13. pl.plot(t,s_mod)
  14. pl.grid(True), pl.xlabel('Tiempo')

Serie de Fourier de una onda triangular.

  1. import pylab
  2. import numpy as np
  3. #Reconstruccion
  4. A=1.
  5. a0=1./2
  6. t=np.linspace(-3,3,10000)
  7. s_per=0
  8. for k in range(-30,30):
  9. #ak=((-1)^k-1)*A/(np.pi^2*k^2)
  10. if k!=0:
  11. s_per=s_per+((-1)**k-1)*A/(np.pi**2*k**2)*np.exp(1j*k*2*np.pi*t)
  12. s_per=s_per+a0
  13. pylab.plot(t,s_per), pylab.show()
  14. pylab.grid(True)

Serie de Fourier de una onda cuadrada.

  1. import pylab
  2. import numpy as np
  3. #Con SINC
  4. k=np.arange(-30,30)
  5. ak=(1/8.0)*np.sinc(k/8.0)*np.exp(-1j*k*np.pi/8.0)
  6. m=np.abs(ak)
  7. f=np.angle(ak)
  8. #pylab.figure(2)
  9. #pylab.stem(k,m,markerfmt=" ")
  10. #pylab.show()
  11. #Reconstruccion
  12. t=np.linspace(-2,2,10000)
  13. s_per=0
  14. for k in range(-30,30):
  15. #print s_per
  16. s_per=s_per+(1/8.0)*np.sinc(k/8.0)*np.exp(-1j*k*np.pi/8.0)*np.exp(1j*k*2*np.pi*t)
  17. pylab.plot(t,s_per), pylab.show()

Convolución continua con Python.

  1. #By Mark Wickert from Signals and Systems For Dummies
  2. def expo_conv(t,A,alpha):
  3. import numpy as np
  4. y = np.zeros(len(t))
  5. for k, tk in enumerate(t):
  6. if tk >= 0:
  7. y[k] = A/alpha*(1 - np.exp(-alpha*tk))
  8. return y
  9.  
  10. import numpy as np
  11. import pylab as pl
  12. import ssd
  13. t = np.arange(-4,14,.01)
  14. xc2 = ssd.step(t)
  15. hc2 = ssd.step(t)*np.exp(-1*t)
  16. yc2_num,tyc2 = ssd.conv_integral(xc2,t,hc2,t,('r','r'))
  17. #pl.subplot(211)
  18. pl.plot(tyc2,yc2_num,linewidth=2)
  19. pl.grid(True)
  20. pl.xlabel('Tiempo (t)')
  21. pl.title('y(t)=x(t)*h(t)')
  22. pl.axis([-8,10,-0.5,1.5])
  23. #pl.subplot(212)
  24. #pl.plot(t,expo_conv(t,1,1))

Convolución discreta con Python.

  1. #By Mark Wickert from Signals and Systems For Dummies
  2. import numpy as np
  3. import pylab as pl
  4. import ssd
  5.  
  6. n=np.arange(-4,6)
  7. xd1=2*ssd.drect(n,4)
  8. hd1=1.5*ssd.dimpulse(n)-0.5*ssd.drect(n+1,3)
  9. #hd1=-0.5*ssd.dimpulse(n+1)+1*ssd.dimpulse(n)-0.5*ssd.dimpulse(n-1)
  10. yd1_num,nd1=ssd.conv_sum(xd1,n,hd1,n)
  11. #pl.subplot(1,3,1)
  12. pl.figure(1)
  13. pl.stem(n,xd1)
  14. pl.title('x[n]'), pl.axis([-4,5,-1,3]),pl.grid(True)
  15. #pl.subplot(1,3,2)
  16. pl.figure(2)
  17. pl.stem(n,hd1),pl.grid(True)
  18. pl.title('h[n]')
  19. #pl.subplot(1,3,3)
  20. pl.figure(3)
  21. pl.stem(nd1,yd1_num), pl.axis([-4,5,-2,2]),pl.grid(True)
  22. pl.title('y[n]=x[n]*h[n]')

Descomposición de una secuencia discreta en una combinación lineal de impulsos desplazados ponderados.

  1. import numpy as np
  2. import pylab as pl
  3. from impulse_step_ramp import *
  4.  
  5. #-----------------------------------------------------
  6. def udelta(n):
  7. y=(n==0)*1 #Convert bool to int [...*1]
  8. #i=np.where(y==1)
  9. return y
  10. #-----------------------------------------------------
  11.  
  12. n=np.arange(-10,10)
  13. s=udelta(n+6)-3*udelta(n+5)-4*udelta(n+4)+1*udelta(n+3)-5*udelta(n+2)+7*udelta(n+1)-1*udelta(n)+\
  14. 3*udelta(n-1)+3*udelta(n-2)-5*udelta(n-3)+2*udelta(n-4)-1*udelta(n-5)-1*udelta(n-6)
  15. pl.stem(n,s)
  16. pl.axis([-10,10,-10,10])
  17. pl.grid(True)
  18. pl.xlabel('Tiempo [n]')
  19. for i, txt in enumerate(n):
  20. pl.annotate(s[i], (n[i],s[i]+0.5))
  21. pl.show

Senoidal, exponencial compleja y señal amortiguada con Python.

Señal coseno.
  1. import numpy as np
  2. import pylab as pl
  3.  
  4. t=np.linspace(0,1,100)
  5. f0=5
  6. s=np.cos(2*np.pi*f0*t-np.pi/2)
  7. pl.plot(t,s)
  8. pl.grid('on')
  9. pl.xlabel('Tiempo [s]')
  10. pl.ylabel('Amplitud')
  11. pl.title('Senoidal')
  12. pl.show()
Exponencial continua real.
  1. import numpy as np
  2. import pylab as pl
  3. n=np.arange(-15,15)
  4. e=np.exp(-0.1*n)
  5. pl.plot(n,e,'r')
  6. pl.grid(True)
  7. pl.show()
Relación de Euler para seno y coseno.
  1. import numpy as np
  2. import pylab as pl
  3.  
  4. t=np.linspace(0,0.5,100)
  5. coseno_s=0.5*(np.exp(1j*4*np.pi*t)+np.exp(1j*4*np.pi*t))
  6. seno_s=0.5*((np.exp(1j*4*np.pi*t)+np.exp(1j*4*np.pi*t)))/1j
  7. pl.subplot(1,2,1)
  8. pl.plot(t,coseno_s)
  9. pl.grid('on')
  10. pl.xlabel('Tiempo [s]')
  11. pl.ylabel('Amplitud')
  12. pl.title('Coseno')
  13. pl.subplot(1,2,2)
  14. pl.plot(t,seno_s)
  15. pl.grid('on')
  16. pl.xlabel('Tiempo [s]')
  17. #pl.ylabel('Amplitud')
  18. pl.title('Seno')
  19. pl.show()
Seno amortiguado.
  1. c=np.sin(2*np.pi*n/15)
  2. s_a=e*c
  3. #pl.stem(n,s_a)
  4. pl.plot(n,e,n,-e)
  5. pl.plot(n,s_a,'r')
  6. pl.grid(True)
  7. pl.show()

Simulación de un sistema para transmisión de datos binarios en Python


Uno de los mejores libros que he revisado este año es Signals and Systems For Dummies de Mark Wickert. En una de sus aplicaciones titulada Simulate the System for Transmitting Binary Data in Python añadí una leve corrección al código en la parte de gráficas y llamado a funciones de la librería ssd.py.

  1. import numpy as np
  2. import pylab as pl
  3. import ssd
  4. from scipy import signal
  5.  
  6. r, b, data0 = ssd.BPSK_tx(100000,10,1.5,0,'src')
  7. # 100000 symbols, Ns=10, Df = 1.5*Rb, 0dB down
  8. r = ssd.cpx_AWGN(r,100,10) # EsN0=100dB, Ns=10
  9. Pr,f = ssd.my_psd(r,2**10,Fs=10)
  10. pl.figure(1)
  11. pl.plot(f,10*np.log10(Pr))
  12. z = signal.lfilter(b,1,r)# b is the SRC filter
  13. #pl.figure(2)
  14. ssd.eye_plot(np.real(z[2000:6000]),20)#20 samp wind
  15. yI,yQ = ssd.scatter(z[2000:6000],10,0)
  16. pl.figure(3)
  17. pl.plot(yI,yQ,'.')
  18. pl.axis('equal')
  19. r, b, data0 = ssd.BPSK_tx(100000,10,1.5,0,'src')
  20. r = ssd.cpx_AWGN(r,20,10) # EsN0=20dB, Ns=10
  21. Pr,f = ssd.my_psd(r,2**10,Fs=10)
  22. pl.figure(4)
  23. pl.plot(f,10*np.log10(Pr))
  24. z = signal.lfilter(b,1,r)
  25. #pl.figure(5)
  26. ssd.eye_plot(np.real(z[2000:6000]),20)#20 samp wind
  27. #ssd.scatter_plot(z[2000:6000]*np.exp(1j*np.pi/5),10,0)
  28. yI,yQ = ssd.scatter(z[2000:6000]*np.exp(1j*np.pi/5),10,0)
  29. pl.figure(6)
  30. pl.plot(yI,yQ,'.')
  31. pl.axis('equal')
Related Posts Plugin for WordPress, Blogger...