Tutorial Simulink 2

Complementos de Procesado de Señal y Comunicaciones Máster en Sistemas Multimedia

Introducción a Matlab y Simulink Javier Ramírez Pérez de Inestrosa Dpto. Teoría de la Señal, Telemática y Comunicaciones Universidad de Granada Email: [email protected] Este tutorial se puede obtener en: http://www.ugr.es/~javierrp

¿Qué es Matlab? 

MATLAB es un lenguaje de alto nivel para realizar cálculos cientifico-técnicos.



Integra las herramientas de cálculo necesarias con otras de visualización así como, un entorno de programación de fácil uso.

Introducción a Matlab y Simulink Javier Ramírez

2

Aplicaciones típicas 

     

Cálculo matemático Desarrollo de algoritmos Adquisición de datos Modelado, simulación y prototipado Análisis de datos y visualización Gráficos Desarrollo de aplicaciones e interfaces gráficas de usuario (GUI) Introducción a Matlab y Simulink Javier Ramírez

3

Más sobre MatLab 

MatLab significa “MATrix LABoratory”



El tipo básico de datos es el vector que no requiere ser dimensionado.



Proporciona unos paquetes de extensión (“toolboxes”) para aplicaciones específicas



Estos paquetes incluyen librerías de funciones MatLab (M-files) que extienden las posibilidades de MatLab para resolver problemas específicos Introducción a Matlab y Simulink Javier Ramírez

4

El entorno de Matlab


5

Sintaxis 

Algunos ejemplos sencillos

Entrada

Salida

Comentarios

2+3 7-5 34*212 1234/5786 2^5

ans = 5 ans = 2 Los resultados son los esperados. ans = 7208 ans = 0.2173 Nótese que al resultado se le da el nombre ans. ans = 32

a = sqrt(2)

a = 1.4142

Se puede escoger el nombre de la variable.


6

b = a, pi, 2 + 3i

b = 1.4142 ans = 3.1416 ans = 2.0000 + 3.0000i

Se pueden introducir varios comandos en una sola línea. Pi, i, y j son constantes.

c = sin(pi) eps

c = 1.2246e-016 ans = 2.2204e-016

"eps" es el limite actual de precisión. No se puede operar con números inferiores a eps.

d= [1 2 3 4 5 6 7 8 9 d = 1 2 3 4 5 6 7 8 9 ] e=123456789 e = [1:9] f=123456789 f = 1:9

Definición de vectores. "d", "e", son "f" vectores. Son iguales. El operador “:” se utiliza para formar vectores; cuenta desde el número inicial al final de uno en uno.

g = 0:2:10 f(3) f(2:7) f(:)

Otros usos de “:”. Se utiliza para acceder a parte o la totalidad de los datos de un vector o matriz.

g = 0 2 4 6 8 10 ans = 3 ans = 2 3 4 5 6 7 123456789


7

(nada) ans = 1 2 3

Un punto y coma ";" evita que se visualice la salida. Una coma simple " ' " calcula la traspuesta de una matriz, o en el caso de vectores, intercambia entre vectores fila y columna.

h * h' h .* h h + h

ans = 14 ans = 1 4 9 ans = 2 4 6

Operaciones con vectores. * es la multiplicación matricial. Las dimensiones deben ser las apropiadas. " .* " es la multiplicación componente a componente.

g = [ 1 2 3; 4 5 6; 7 8 9]

g = 1 2 3 4 5 6 7 8 9

Construcción de matrices.

ans = 6 ans = 7 8 9 g = 1 2 3 4 5 4 7 8 9

Accediendo a los elementos de la matriz. ":" se utiliza para acceder a una fila completa.

h = [1 2 3]; h'

g(2,3) g(3,:) g(2,3) = 4


8

Entrada

g^2

g .^ 2

Salida

Comentarios

ans =

30 36 42 66 81 96 Multiplica la matriz por ella misma. 102 126 150 ans = 1 4 9 Eleva al cuadrado cada elemento de la 16 25 36 matriz. 49 64 81


9

Control de la salida 

El comando format 

format compact 



format long 



Muestra únicamente cinco dígitos.

“;” al final del comando. 



Muestra los 15 dígitos que se utilizan en el cálculo.

format short 



Controla el espaciado de líneas.

No visualizar salida:

help format 

Más información. Introducción a Matlab y Simulink Javier Ramírez

10

Más sobre matrices 

Funciones incluidas en MatLab

Entrada

Salida

rand(2)

ans = 0.9501 0.2311 ans = 0.8913 0.7621

rand(2,3) zeros(2) ones(2)

ans = 0 0 ans = 1 1

Comentarios 0.6068 0.4860 0.4565 0.8214 0.0185 0.4447

0 0 1 1

Genera una matriz de números aleatorios entre 0 y 1 Genera una matriz 2x2 de ceros o unos.

eye(2)

ans = 1 0 0 1

Matriz identidad I.

hilb(3)

ans = 1.0000 0.5000 0.3333 0.5000 0.3333 0.2500 0.3333 0.2500 0.2000

Matriz de Hilbert 3x3.


11


Concatenación 



Generar nuevas matrices a partir de otras creadas previamente Por ejemplo: 

Sea la matriz a:

>> a = [1 2; 3 4] a = 1 2 3 4 Introducción a Matlab y Simulink Javier Ramírez

12

Más sobre matrices - concatenación Entrada

Salida

[a, a, a]

ans = 1 2 1 2 1 2 3 4 3 4 3 4

[a; a; a]

ans = 1 3 1 3 1 3

2 4 2 4 2 4

[a, zeros(2); zeros(2), a']

ans = 1 3 0 0

2 4 0 0

0 0 1 2

0 0 3 4


13


Programación 

Se pueden construir matrices mediante programación

for i=1:10, for j=1:10, t(i,j) = i/j; end end 



No se produciría salida puesto que la única línea que podría generar salida (t(i,j) =i/j;) termina en “;” Sin el “;”, Matlab escribiría la matriz t 100 veces!! Introducción a Matlab y Simulink Javier Ramírez

14

Operaciones con matrices 

+, -, *, y / 



Definen operaciones con matrices.

Debemos distinguir: 



“.*”: Multiplicación componente a componente. “*” Multiplicación matricial. Introducción a Matlab y Simulink Javier Ramírez

15

Escalares 

Un escalar es un número.



Matlab los almacena como matrices 1x1



Todas las operaciones entre escalares y matrices se realizan componente a componente salvo: 

La potencia (“^”). Introducción a Matlab y Simulink Javier Ramírez

16

Escalares 

a = 1 3

Ejemplos

2 4

Entrada

Salida

Comentarios

b=2

b=2

Define b como un escalar.

a + b

ans = 3 4 5 6

La suma se hace componente a componente.

a * b

ans = 2 4 6 8

Igual que la multiplicación.

a ^ b

ans = 7 10 15 22

Potencia matricial - a*a

a .^ b

ans = 1 4 9 16

Potencia componente a componente.


17

Vectores 

Un vector es una matriz de una sola fila o columna Entrada

Salida

Comentarios

v = [1 2 3] u = [3 2 1]

v = 1 2 3 u = 3 2 1

Define 2 vectores.

v * u

Error

Las dimensiones no coinciden.

v * u'

ans = 10

Al tomar la traspuesta se corrige el error.

dot(v,u)

ans = 10

Producto escalar (idéntico al anterior).

cross(v,u)

ans = -4 8 -4

El producto vectorial sólo se emplea con vectores en 3 dimensiones.


18

Matrices 

Matlab tiene numerosas funciones predefinidas (help matfun).

Entrada

Salida

Comentarios

k = [16 2 3; 5 11 10; 9 7 6]

k = 16 2 3 5 11 10 9 7 6

Define una matriz.

trace(k)

ans = 33

Traza de una matriz

rank(k)

ans = 3

Rango de una matriz.

det(k)

ans = -136

Determinante de una matriz


19

Matrices Entrada

Salida ans =

inv(k)

[vec,val] = eig(k)

Comentarios 0.0294 -0.0662 0.0956 -0.4412 -0.5074 1.0662 0.4706 0.6912 -1.2206

vec = -0.4712 -0.4975 -0.0621 -0.6884 0.8282 -0.6379 -0.5514 0.2581 0.7676 val = 22.4319 0 0 0 11.1136 0 0 0 -0.5455


Inversa de una matriz Vectores propios y autovalores de una matriz. Las columnas de "vec" contienen los vectores propios; las entradas de la diagonal de "val" son los autovalores.

20

Variables en el espacio de trabajo 

whos  Lista las variables definidas en el entorno. >> whos Name

a b c

Size

Bytes

Class

100x1 100x100 1x1

800 80000 8

double array double array double array

Grand total is 10101 elements using 80808 bytes 

clear  Borra variables del entorno. Introducción a Matlab y Simulink Javier Ramírez

21

Resolución de sistemas de ecuaciones 



Una de las principales aplicaciones de las matrices es la representación de sistemas de ecuaciones lineales. Si a es una matriz de coeficientes, x es un vector columna que contiene las incógnitas y b los términos constantes, la ecuación a x =b

representa el correspondiente sistema de ecuaciones. Introducción a Matlab y Simulink Javier Ramírez

22

Resolviendo ecuaciones 

Para resolver el sistema en MatLab  



x = a \ b x es igual a la inversa de a por b

Ejemplo 



a = [1 2 3; 4 5 6; 7 8 10]; b = [1 1 1]'; Solución: x = -1 1 0 Introducción a Matlab y Simulink Javier Ramírez

23

Salvar y recuperar datos 

Los datos de la sesión se pierden al salir de MatLab.



Para salvar la sesión (entrada y salida) 

Diary(‘session.txt’); 





Guarda los comandos introducidos en la sesión.

Diary ;

Para salvar una o varias matrices  

save datos.mat (guarda todas las variables) save datos.mat x (sólo guarda x) Introducción a Matlab y Simulink Javier Ramírez

24

Salvar y recuperar matrices 

save sesion 



save fichero X 



Salva sólo la variable X

load sesion 



Salva todas las variables en el archivo binario “sesion.mat”.

Recupera los datos previamente salvados

Si los ficheros se pueden salvar en formato texto (ascii). Pueden verse con un editor de textos. Introducción a Matlab y Simulink Javier Ramírez

25

Gráficos 

  

El comando básico es: plot plot(y); plot(x,y); plot(x,y,’b+’,x,z,’gx’);  



color (b,g) blue,green Marcador (+,x)

Personalización del gráfico: 

title, xlabel, ylabel, legend, grid. Introducción a Matlab y Simulink Javier Ramírez

26

Ejemplo



Mes = 1:12; T_Gr = [-2 0 2 4 8 12 14 14 12 8 4 0]; T_Ma = [-4 -2 0 2 6 14 18 18 16 8 2 -2];



plot(Mes, T_Gr, 'bo', Mes, T_Ma, 'rv');



xlabel('Mes'); ylabel('Temperatura (°C)'); title('Temperaturas minimas en Granada y Madrid'); legend('Granada','Madrid'); grid;

 

 

 


27

Resultado:


28

subplot  



   

 

     

Mes = 1:12 T_Gr = [-2 0 2 4 8 12 14 14 12 8 4 0]; T_Ma = [-4 -2 0 2 6 14 18 18 16 8 2 -2]; subplot(2,1,1); plot(Mes, T_Gr, 'bo-'); xlabel('Mes'); ylabel('Temperatura (°C)'); title('Temperaturas minimas en Granada'); grid; subplot(2,1,2); plot(Mes, T_Ma, 'rv-'); xlabel('Mes'); ylabel('Temperatura (°C)'); title('Temperaturas minimas en Madrid'); grid; Introducción a Matlab y Simulink Javier Ramírez

29

Resultado


30

Gráficos tridimensionales   

[x,y] = meshgrid(-2:.2:2, -2:.2:2); z = x .* exp(-x.^2 - y.^2); mesh(z);


31

Gráficos tridimensionales  

   

Z = X.*exp(-X.^2-Y.^2); contour3(X,Y,Z,30) surface(X,Y,Z,'EdgeColor',[.8 .8.8],'FaceColor','none') grid off 0.5 0.4 view(-15,25) 0.3 colormap cool 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4 -0.5 2 1 0 -1 -2

-2

-1.5

-1

-0.5

0


0.5

1

1.5

2

32

Programación 

Ficheros de comandos (scripts)   



Secuencias de comandos. Al invocarlos se ejecutan en el entorno. Las variables creadas son globales.

Ficheros de función    

Permiten definir funciones propias. Variables locales. La información se pasa como parámetros. Se pueden definir subfunciones. Introducción a Matlab y Simulink Javier Ramírez

33

Un ejemplo de función function y = media (x) % Valor medio de x. % Para vectores, media(x) devuelve el valor medio. % Para matrices, media(x) es un vector fila % que contiene el valor medio de cada columna. [m,n] = size(x); if m == 1 m = n; end y = sum(x)/m;


34

Argumentos de funciones 

nargin y nargout 



Número de argumentos de entrada y salida con los que se llama a la función.

Ejemplo: function c = testarg1(a,b) if (nargin == 1) c = a.^2; elseif (nargin == 2) c = a + b; end Introducción a Matlab y Simulink Javier Ramírez

35

Subfunciones function [media,mediana] = estadistica(u) % Función principal % ESTADISTICA Calcula la media y la % mediana utilizando funciones internas. n = length(u); media = mean(u,n); mediana = median(u,n); function a = mean(v,n) % Calcula la media. a = sum(v)/n;

% Subfunción

function m = median(v,n) % Calcula la mediana. w = sort(v); if rem(n,2) == 1 m = w((n+1)/2); else m = (w(n/2)+w(n/2+1))/2; end

% Subfunción


36

Bifurcaciones Condición 1

Condición true

Sentencias

false

Condición

false

Condición 2

true Bloque 1

false

Bloque 2

false

true Bloque 1 true Bloque 2

Bloque 3


37

if Sentencia if

Bifurcación multiple

if condición sentencias end

if condición1 bloque1 elseif condición2 bloque2 elseif condición3 bloque3 else bloque4 end

Ejemplo if rem(a,2) == 0 disp('a is par') b = a/2; end


38

switch switch expresion case case1, bloque1 case {case2, case3, ...} bloque2 ... otherwise, switch valor bloque3 case -1 end disp('negativo'); case 0 disp('cero'); case 1 disp('positivo'); otherwise disp('otro'); end Introducción a Matlab y Simulink Javier Ramírez

39

Bucles

false

Condición

Sentencias

true Sentencias

Condición

true

false


40

for for i = 1:n sentencias end

for i = vector sentencias end

for i = n:-0.2:1 sentencias end

for i = 1:m for j = 1:n sentencias end end


41

while

{continue, break}

while condición sentencias end

Sentencia break Hace que termine la ejecución

Sentencia continue Hace que se pase inmediatamente a la siguiente iteración del bucle for o while

n = 1; while prod(1:n) < 1e100 n = n + 1; end


42

Toolboxes 

Toolboxes para DSP y comunicaciones:   

   

Communications Toolbox Filter Design Toolbox Image Processing Toolbox Signal Processing Toolbox Statistics Toolbox System Identification Toolbox Wavelet Toolbox


43

Procesamiento de señales y comunicaciones 

Matlab dispone de unas librerías para tratamiento digital de señales. 

Signal Processing Toolbox



Communications Toolbox


44

Filtrado de señales 

y= filter(b,a,x); 

Filtra la secuencia x con el filtro descrito por b y a.


45

Respuesta en frecuencia [H,w] = freqz(b,a,N);



Calcula N puntos de la respuesta en frecuencia del filtro definido por b y a. Para el filtro: y(n)= 0.8·y(n-1) + x(n)



20

Magnitude (dB)



Phase (degrees)



freqz(1,[1 -0.8],256);

10 0 -10 0

0.2 0.4 0.6 0.8 Normalized Frequency ( rad/sample)

1

0.2 0.4 0.6 0.8 Normalized Frequency ( rad/sample)

1

0 -20 -40 -60 0


46

Diseño de filtros 

B = FIR1(N,Wn);  



B = FIR1(N,Wn,'high'); 



Filtro paso alta.

Wn = [W1 W2]; B = FIR1(N,Wn,'bandpass'); 



Filtro FIR paso baja de orden N. Wn es la frecuencia de corte normalizada (0
Filtro paso banda.

Wn = [W1 W2]; B = FIR1(N,Wn,'stop'); 

Filtro rechaza banda (notch). Introducción a Matlab y Simulink Javier Ramírez

47

Ejemplo: Filtros FIR paso banda % Respuesta en frecuencia

for i=1:L B = FIR1(N(i),Wn,'bandpass'); [H1,W]= freqz(B,1,NFFT); H(:,i)= H1; end

-40 -60

% Diseño. % Respuesta.

% Visualizacion de la respuesta en frecuencia. subplot(2,1,1); plot(W/pi,20*log10(abs(H))); xlabel('w en unidades de \pi'); ylabel('|H|'); legend('N=8','N=16','N=32','N=64'); grid; subplot(2,1,2); plot(W/pi,unwrap(angle(H))); xlabel('w en unidades de \pi'); ylabel('Fase de H (rads)'); legend('N=8','N=16','N=32','N=64'); grid;

N=8 N=16 N=32 N=64

-20

|H|

NFFT= 256; L = length(N); H = zeros(NFFT,L);

0

% Coeficientes del filtro % Frecuencias de corte.

-80 0

Fase de H (rads)

% Parametros. N= [8 16 32 64]; Wn = [0.4 0.6];

0.2

0.4 0.6 w en unidades de 

0.8

10

1 N=8 N=16 N=32 N=64

0 -10 -20 -30 0

0.2

0.4 0.6 w en unidades de 


0.8

1

48

Diseño y análisis de filtros (fdatool)


49

Signal Processing Tool (SPTool)


50



51



52



53

Diezmado e interpolación Señal original

2 1.5



Y = DECIMATE(X,R)

1 0.5 0

Ejemplo: Decimación en un factor 4. t = 0:.00025:1; % Vector de tiempos x = sin(2*pi*30*t) + sin(2*pi*60*t); y = decimate(x,4);

-0.5 -1 -1.5 -2 0

20

40

80

100

120

20

25

30

Señal decimada

2

stem(x(1:120)); %Original axis([0 120 -2 2]) title('Señal original') figure stem(y(1:30)); %Decimada title('Señal decimada')

60

1.5 1 0.5 0 -0.5 -1 -1.5 -2 0

5


10

15

54

Diezmado e interpolación Señal original

2 1.5



Y = INTERP(X,R)

1 0.5 0 -0.5

t = 0:0.001:1; % Time vector x = sin(2*pi*30*t) + sin(2*pi*60*t); y = interp(x,4); stem(x(1:30)); title('Señal original'); figure stem(y(1:120)); title('Señal interpolada');

-1 -1.5 -2 0

5

10

15

20

25

30

80

100

120

Señal interpolada

2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 0

20


40

60

55

Densidad de potencia espectral periodogram Fs = 1000; t = 0:1/Fs:.3; % Una señal coseno de 200Hz más ruido x = cos(2*pi*t*200)+randn(size(t)); periodogram(x,[],'twosided',512,Fs); % Se usa la ventana por defecto 

pwelch Fs = 1000; t = 0:1/Fs:.296; % Una señal coseno de 200 Hz más ruido x = cos(2*pi*t*200)+randn(size(t)); pwelch(x,[],[],[],Fs,'twosided'); % Ventana por defecto, solapamiento y NFFT. 


56

Ejemplos Periodogram PSD Estimate

-15

-18

-20

-20

-25 -30 -35 -40 -45 -50 0

Welch PSD Estimate

-16

Power Spectral Density (dB/Hz)


-10

-22 -24 -26 -28 -30 -32

100

200

300

400 500 600 Frequency (Hz)

700

800

900

-34 0

100

200

300

400 500 600 Frequency (Hz)


700

800

900

57

Estimación espectral paramétrica 

Método de covarianza pcov(X,ORDER,NFFT,Fs) randn('state',1); x = randn(100,1); y = filter(1,[1 1/2 1/3 1/4 1/5],x); pcov(y,4,[],1000);



Método de covarianza modificado pmcov(X,ORDER,NFFT,Fs) randn('state',1); x = randn(100,1); y = filter(1,[1 1/2 1/3 1/4 1/5],x); pmcov(y,4,[],1000); Introducción a Matlab y Simulink Javier Ramírez

58

Ejemplos Covariance PSD Estimate

-24

-26



-26

-28

-30

-32

-34

-36

-38 0

100

200 300 Frequency (Hz)

400

-30

-32

-34

-38 0

500

Burg PSD Estimate

100


400

500

400

500

Yule-Walker PSD Estimate

-24

-26


-26


-28

-36

-24

-28

-30

-32

-34

-36

-38 0

Modified Covariance PSD Estimate

-24

-28

-30

-32

-34

-36

100


400

500

-38 0

100



59

Modelado. Predicción lineal 

Predictor lineal: 

Estima la muestra siguiente utilizando N muestras anteriores. x(n)

Predictor lineal



xˆ (n)  



N

i 1

xˆ (n)

e(n)  x(n)  xˆ(n)

A(i  1) x(n  i)

[A,E] = LPC(X,N) 

Calcula:  

Los coeficientes del predictor, A. La varianza del error, E.

Minimizando J= E{|e(n)|2}


60

Comunicaciones 

Funciones principales:   



  

Generación de señales aleatorias Análisis de errores Codificación de la fuente (escalar, diferencial) Codificación para el control de errores (convolucional, codificación lineal de bloques) Modulación y demodulación (analógica y digital) Filtrado mediante filtros especiales Aritmética en cuerpos de Galois Introducción a Matlab y Simulink Javier Ramírez

61

Modulación analógica y digital


62

Funciones utilizadas frecuentemente 

randint 

Genera una matriz de números aleatorios enteros uniformemente distribuidos.  

 

out out out out

= = = =

randint; randint(N); randint(M,N); randint(M,N,rg);

% % % %

Un solo valor [0,1] Matriz cuadrada NN Matriz MN Se especifica el rango.


63

Modulación/demodulación digital


64

Modulación/demodulación digital 

“Mapping”+modulación/demodulación 

dmodce/ddemodce y y y y y y y y y



= = = = = = = = =

dmodce(x,Fd,Fs,'method/nomap'...); dmodce(x,Fd,Fs,'ask',M); dmodce(x,Fd,Fs,'fsk',M,tone); dmodce(x,Fd,Fs,'msk'); dmodce(x,Fd,Fs,'psk',M); dmodce(x,Fd,Fs,'qask',M); dmodce(x,Fd,Fs,'qask/arb',inphase,quadr); dmodce(x,Fd,Fs,'qask/cir',numsig,amp,phs); dmodce(x,Fd,[Fs initphase],...);

Sólo “mapping” 

modmap/demodmap modmap('method',...); y = modmap(x,Fd,Fs,'ask',M); y = modmap(x,Fd,Fs,'fsk',M,tone); y = modmap(x,Fd,Fs,'msk'); y = modmap(x,Fd,Fs,'psk',M); y = modmap(x,Fd,Fs,'qask',M); y = modmap(x,Fd,Fs,'qask/arb',inphase,quadr); y = modmap(x,Fd,Fs,'qask/cir',numsig,amp,phs);


65

Ejemplo M = 16; Fd = 1; Fs = 3; x = randint(100,1,M);

% % % % % %

Modulación M-aria. Se muestrea el mensaje original a una muestra por segundo. La señal modulada se muestrea a una frecuencia de 3 muestras por segundo. Mensaje digital aleatorio.

% Modulación M-ary PSK y = dmodce(x,Fd,Fs,'psk',M); % Adición de ruido Gaussiano. ynoisy = y + .04*randn(300,1) + .04*j*randn(300,1); % Diagrama de dispersión a partir de las observaciones ruidosas. scatterplot(ynoisy,1,0,'b.'); % Demodulación para recuperar el mensaje z = ddemodce(ynoisy,Fd,Fs,'psk',M); s = symerr(x,z) % Comprobar la tasa de error de los símbolos. Introducción a Matlab y Simulink Javier Ramírez

66

Gráfico de dispersión (ruido) Scatter plot 1 0.8 0.6

Quadrature

0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 -1

-0.5

0 In-Phase

0.5

1


67

Análisis de errores 

biterr 

Calcula el número de bits erróneos y la tasa de error.

[number,ratio] = biterr(x,y); [number,ratio] = biterr(x,y,k); [number,ratio] = biterr(...,flg); [number,ratio,individual] = biterr(...)


68

Ejemplo x = randint(100,100,4);

% Señal original

% Introducción de errores. % Los errores pueden ser 1, 2, o 3 (no 0). % Colocación de los errores errorplace = (rand(100,100) > .9); errorvalue = randint(100,100,[1,3]); % Valor error errors = errorplace.*errorvalue; y = rem(x+errors,4); % Señal y error sumadas mod 4 % Análisis de errores format short [num_bit,ratio_bit] = biterr(x,y,2) [num_sym,ratio_sym] = symerr(x,y) Introducción a Matlab y Simulink Javier Ramírez

69

Ejercicio: 

Estudiar mediante Matlab/Simulink el efecto del ruido en los sistemas de comunicación digital.  

Construir el diagrama de bloques de simulación. Simular el sistema: 

 

Diferentes esquemas de modulación (ASK, PSK, FSK, MSK). Obtener las curvas de error en función de la SNR. Adicionalmente, considérese PSK con M= 4 y 8 símbolos y compárese las curvas de error frente a la SNR. Introducción a Matlab y Simulink Javier Ramírez

70

Introducción a Simulink 



Simulink es una herramienta para modelado, simulación y análisis de sistemas dinámicos. Soporta tanto sistemas lineales como no lineales:   



en tiempo continuo, muestreados, híbridos y sistemas multifrecuencia (contienen sistemas muestreados a diferente frecuencia). Introducción a Matlab y Simulink Javier Ramírez

71

Entorno de trabajo


72

Construcción del diagrama de bloques


73

Resultado de la simulación 1 0.5 0 -0.5 -1 0

500

1000

0.5

1

1500

2000

2500

3000

1 0.5 0 -0.5 -1 0

1.5

2

2.5

3 4

x 10


74

Ejemplo


75

Simulación


76

Análisis/síntesis LPC de la señal de voz


77

Estimación espectral


78

Resultado de la simulación


79

Sistema de comunicación digital

8-PSK (Gray)


80

Procesado de imágenes 

Lectura de imágenes: 

A = imread(filename, fmt) 

 

[X, map] = imread(...) 



Lee una imagen en escala de grises o en color fmt especifica el formato de imagen (BMP, JPEG, PNG, TIFF, etc)

map

Mapa de color

Visualización de la imagen imshow     

imshow(I) imshow(I,[low high]) imshow(RGB) imshow(BW) imshow(X,map)


81

Lectura de imágenes   

A= imread('texto.tif'); imshow(A); size(A); size(A); 500500





[B map]= imread('flor.tif'); imshow(B); size(A); 5005003


82

Transformaciones espaciales 

Modificación del tamaño  



B = imresize(A, scale) B = imresize(A, [mrows ncols])

Rotar una imagen  

B = imrotate(A,angle) B = imrotate(A,angle,method) 



->

Interpolación {‘nearest’}, ‘bilinear’, ‘bicubic’

B = imrotate(A,angle,method,bbox) 



method bbox

->

Bounding box

{‘crop’} ‘loose’

Recortar una imagen  

B= imcrop(A); B= imcrop(A,rect);

Herramienta interactiva


83

Filtrado 2D 

B= imfilter(A,h)

 



Original

N= 5; h = ones(N,N) / (N*N); Af = imfilter(A,h);

Filtrada N= 5 Introducción a Matlab y Simulink Javier Ramírez

Filtrada N= 20 84

Respuesta en frecuencia de filtros 2D [H, [H, [H, [H,

 

 

f1, f1, f1, f1,

f2] f2] f2] f2]

= = = =

freqz2(h, n1, n2) freqz2(h, [n2 n1]) freqz2(h) freqz2(h, f1, f2)

1

1

0.8 Magnitude

Magnitude

0.8 0.6 0.4

0.6 0.4 0.2

0.2

0 1

0 1 0.5

1 0.5

0

1

-0.5 -1

-1

Fx

0.5

0

0

-0.5 Fy

0.5 0

-0.5

Filtro 55

Fy

-0.5 -1

-1


Fx

Filtro 2020 85

Diseño de filtros 2D

1 0.8 0.6

 



h = fwind1(Hd, win) h = fwind1(Hd, win1, win2) h = fwind1(f1, f2, Hd,...)

0.4 0.2 0 1 0.5 0

Ejemplo:      

 

0

-0.5

-0.5

[f1,f2] = freqspace(21,'meshgrid'); 1.5 Hd = ones(21); r = sqrt(f1.^2 + f2.^2); Hd((r<0.1)|(r>0.5)) = 0; 1 colormap(jet(64)); mesh(f1,f2,Hd); Magnitude



1 0.5

h = fwind1(Hd,hamming(21)); freqz2(h);

-1

-1

0.5

0 1 0.5

1 0.5

0

0

-0.5 Fy

-0.5 -1

-1


Fx

86

Ajuste del nivel de intensidad

6000

15000

4000

10000

4000 3000 2000

2000

5000

0

0

0

100

200

1000 0

100

200

0

0

100


200

87

Ecualización del histograma

250

250

250

200

200

200

150

150

150

100

100

100

50

50

50

0

0

100

200

0

0

100

200

0

0

100


200

88

Ecualización del histograma

6000

15000

4000

10000

2000

5000

4000

2000

0

0

50

100

150

200

250

0

0

50

100

150

0

200

200

200

200

100

100

100

0

0

50

100

150

200

250

0

0

50

100

150

200

250

15000

2000

0

0

50

100

150

200

250

0

50

100

150

200

250

0

50

100

150

200

250

2000

10000 1000 0

1000

5000 0

50

100

150

200

250

0

0

50

100

150

200

250

0


89

Filtrado de ruido (Filtro de Wiener)  

J = wiener2(I, [m n], noise) [J, noise] = wiener2(I, [m n]) 

Wiener2 estima la media y la varianza entorno a cada pixel



A continuación crea un filtro pixel a pixel basado en estas estimaciones

v2 es la varianza del ruido


90

Matlab para reconocimiento de patrones 

Asignación de una clase a un vector de características x del objeto a clasificar: 

Ejemplo:  

Fisher Iris dataset:

http://en.wikipedia.org/wiki/Iris_flower_data_set 503 = 150 muestras de flores Iris de tres especies  

Iris setosa, Iris virginica, Iris versicolor 4 características de cada ejemplo:  Longitud y anchura de los pétalos y sépalos Introducción a Matlab y Simulink Javier Ramírez

91

Matlab para reconocimiento de patrones 

Análisis discriminante:   



knn:    



class = classify(sample,training,group) class = classify(sample,training,group,type) class = classify(sample,training,group,type,prior)

Class Class Class Class

= = = =

knnclassify(Sample, knnclassify(Sample, knnclassify(Sample, knnclassify(Sample,

Training, Training, Training, Training,

Group) Group, k) Group, k, distance) Group, k, distance, rule)

Máquinas de vectores de soporte:    

Group = svmclassify(SVMStruct, Sample) Group = svmclassify(SVMStruct, Sample, 'Showplot‘,ShowplotValue) SVMStruct = svmtrain(Training, Group) SVMStruct =svmtrain(..., 'Kernel_Function', Kernel_FunctionValue)


92

Análisis discriminante: Ejemplo 3.8 3.6



 

load fisheriris SL = meas(51:end,1); SW = meas(51:end,2); group = species(51:end);

3.4 3.2 3 SW



Fisher versicolor Fisher virginica

2.8 2.6 2.4 2.2 2 4.5

5

5.5

6

6.5

7

7.5

8

SL



 

h1 = gscatter(SL,SW,group,'rb','v^',[],'off'); set(h1,'LineWidth',2) legend('Fisher versicolor','Fisher virginica','Location','NW')


93

Análisis discriminante: Ejemplo 

Clasificamos una matriz de datos 2-D: [X,Y] = meshgrid(linspace(4.5,8) ,linspace(2,4)); X = X(:); Y = Y(:); [C,err,P,logp,coeff] = classify([X Y],[SL SW], group,'quadratic');



Visualizar la clasificación: hold on; gscatter(X,Y,C,'rb','.',1,'off'); K = coeff(1,2).const; L = coeff(1,2).linear; Q = coeff(1,2).quadratic; f = sprintf('0 = %g+%g*x+%g*y+%g*x^2+%g*x.*y+%g*y.^2', ... K,L,Q(1,1),Q(1,2)+Q(2,1),Q(2,2)); h2 = ezplot(f,[4.5 8 2 4]); set(h2,'Color','m','LineWidth',2) axis([4.5 8 2 4]) xlabel('Sepal Length') ylabel('Sepal Width') title('{\bf Classification with Fisher Training Data}')


94

Análisis discriminante: Ejemplo Classification with Fisher Training Data 4 3.8

Fisher versicolor Fisher virginica

3.6 3.4

SW

3.2

3 2.8 2.6 2.4 2.2

2 4.5

5

5.5

6

6.5

7

7.5

8

SL


95

Selección de características  

[IDX, Z]= rankfeatures(X, Group) [IDX, Z]= rankfeatures(X, Group,'Criterion', CriterionValue) load fisheriris; X= meas(1:100,:); Group= species(1:100); [IDX, Z] = rankfeatures(X',Group); % Selecionamos las variables más discriminativas data= X(:,[IDX(1) IDX(2)]);

% Selección aleatoria de subconjuntos de entrenamiento y test [train, test] = crossvalind('holdOut',Group); cp = classperf(Group); % Entrenamiento de una máquina de vectores de soporte svmStruct = svmtrain(data(train,:),Group(train),'showplot',true); % Añadimos título. title(sprintf('Kernel Function: %s',... func2str(svmStruct.KernelFunction)),... 'interpreter','none'); % Clasificación del conjunto de test classes = svmclassify(svmStruct,data(test,:),'showplot',true); % Evaluación a partir de la tasa de correctas. classperf(cp,classes,test); cp.CorrectRate


96

Selección de características Kernel Function: linear_kernel

1.8

Kernel Function: linear_kernel 4.5

setosa (training) setosa (classified) versicolor (training) versicolor (classified) Support Vectors

1.6

1.4

setosa (training) setosa (classified) versicolor (training) versicolor (classified) Support Vectors

4

1.2 3.5 1

0.8 3 0.6

0.4

2.5

0.2

0

1

1.5

2

2.5

3

3.5

4

4.5

5

data= X(:,[IDX(1) IDX(2)]);

5.5

2

4

4.5

5

5.5

6

6.5

data= X(:,[IDX(3) IDX(4)]);


97

7

Estas transparencias se pueden obtener en: http://www.ugr.es/~javierrp Para cualquier consulta: Javier Ramírez ([email protected]) Dpto. Teoría de la Señal, Telemática y Comunicaciones

Despacho 22 ETSII Introducción a Matlab y Simulink Javier Ramírez

98

!

"

#

$

%

!"

"#

%$"

!

" $ &

#

'"

$

%

%

'

$

$

) (( * , + ) -

(( ) * (( ) + (( ) , ) *

+

#

) ((

./

./

0

1

, -"

"

% -

.

)

-

/ $ / ./ . / $ 0 1

#

$

) (( 4-5 6 7 89 ((

&

) %

&

$

(

((

#

2 ) 3! * +

(( % (( % (( % (( %

2 ) 3!$*$+:4$-$5:6$7$89

*" # (( ) * , +>= (( )

= " @ A

?( ?(

*$+ * +

; ;

) 3! * +9 ) 3!$*$+9 ) 3!$*$+9< -< ) 3!:*:+9

.

(( ((

)

)

&

@ A > !7B C

'

?( ?(

2 2 2$2 2

2 2 2 2

2 2

2

3 +" D

"

(( (( ((

) 3! ?5 B ?*69 ) @ A ) @

A

) - .; - .< - .= - .> - .; - .? - .> - . 56 9?9 ?7 56 @ 9?A ?8

(( ) 3! * +9 (( ) 34 - 59 (( ) %@ $ A ) 4 !+ *7 *6 !7 (( 30$ 9 ) %@ $ A 0) 4-5 ) BBBBB , -

$ /

'

56 56

%

%.

*E +E" )

56 56 $ 56 $

( (( ) 3B'!'!BB9 (( 1 ) (( ) 3 @B"!> A : (( @1$ A -" F

)

( 7D) (( ) 3B'!'!BB9 56 (( 1 ) 56 $ (( ) @B"!> A 56 $ (( 2 ) @B"*> A 56 $ (( +@1$ $2A

@B"!> A9

; $

-.

56 - .

;

( (( ) 3B'!'!BB9 (( 1 ) (( ) @B"!> A (( @1$ A

-.

; @ @?B7 7?9C 9?9: ?BC

56 - .

4"

/

#

56 758 959: 56 & 56

%

/

#

" )

4

(( (( ((

) 3! *9 ) 3! + -9 @ $ A @ ,*AC@ G*,+ ,-A (( 2 @ $ A (( 3 $29 ) 2 @ $ ((

2

56 ?9 56 9?7 ?A 56 56 56

A

@ $2A

* +

56

/

)

(( . ) +B (( 1 ) B"+ (( 3 *$ *9 ) *@. $1 A (( @ *$ *A !C@ G*,!7 ,8BBA

F

* +

AE

"

!

$

"

$

56

'

I

>

(( ) ;@ $ A (( 3 $29 ) 2 @ A (( ) @ > ;@ (( ) @ A );

L@ $!$?!A

5"

@2AAA>

>

I J I

@2A

>

J)

K (( ) ;@ (( ) (( )

56

$ @ " @ "

) A

J'KA J'KA

*B+

25B2

"

((

@ A

56

@ A

56

(( 3 $1$ 9 ) ((

@ A$

J I J

$

1

@ A )

1@ A

56 56

*+

)

" ) (( ) ;@3! *9$ 3! + -9A 56 ;@ ,*AC@ G*,+ ,-A (( 2 @ A

56 ?9 56 9?7 ?A

A

H

)

66 ) ;@I IA 56 (( ) @ ,*AC@ G*,+> ,-A @ ,*AC@ G*,+ ,-A (( * ) >

* +

F

&

(( ) 3! *9 (( ) 3! + -9 (( 2 @ $

((

G

$

H

(( ) ;@3! *9$ 3! + -9A @ ,*AC@ G*,+ ,-A

((

56 56 56 56

L 2 2

22

2 2

6"

;

(( (( 3

@ A $; 9 )

(( 3 (( ((

"

$; 1

$.9 ) )

56 56

@ $+A

1@ ) *B>

@ A

56

A !B@

56 A 56

1

(( 3M $D $. $. 9 ) M ) D ) . )Q Q . )Q Q

@ A

. $

56 @ !A 56 56 56 56 56 56 56 56 56 56 56

((

56

@ A

N @ A ) @

(( (( &

Q

@ $

OOO . *

-

+

.

-

>. >. 5

Q -

2

Q

2.

"

*D

/

+ $

" $ $3+ 49A

$+A

56 56 56 $ 56 $ 56 56

'

N

56

7" (( (( (( / ((

)

F

6E % .@ A

. ;7

56 MMM ( P

(( / @ A (( @I .IA (( / (( ; (( / (( / ;;

((

-

-

-Q Q ; Q

(( ! ) !BC@@ ,!A>@ ,*A>@ ,+AA (( 3M !$D !$. !$. !9 ) M !) D !) . !) . !) (( @ !A

0

$

&

N -

N

2R;72 $ & . 2R;72 2R;@2 '

* .

(

+ #

O !"

L" %

R

N

-

. $

/

$ & "

$ N

$

%

R

"

%

" R

%

R

2%

R

! 1

2 -!

.

&

) 66

R L

$ .)

$ $

% '

'

%

& '

# %

%

R-

.)

)

$

" $ &

%

/ R

(

%

R %

' ; %

E

0

1) / &

0 #

"

$

R) N

$

% )

. )

/

# ) E

1) / &

& / E

& ;

F )/ $

) # ) D $

%

# E $ P

?D

$

E

)

)

E

)F

S

D ),

%

D

%

' /'

0

), ) / F

M

0 "

L'

2 E

R '

)( )( S L

) "

"

)S

N $

&

% N 2

' D

) ,

F

L

)S

2

N

#

) D ) ,

?B

E

5B

E ) ,

'F

L /

# "$ 2??5?2

@

' E

P

N

L)

) "

)/ S L

F

'

"

)

&

F

R S

N ) /

&

0

$ "

N

)P

N

F

F

)/

2

" D

' )

;B

$ .

M

) S

-

;E

&

" )

2

$

)S

N

N ) S

.

#

#

/) L

M

Q

'

=

R

)/

"

) %

S % )S )S

$

N

; !) /

/ N

& !) /

/ N

N

) S #

#

*"

" / %

+

' -

' %

)

1

* '

, 0

. /

! &#

#

'

)

'

B.

%

)/

2 9.

0

2 )% #

2 7.

0

2-

"

.

) / $

@.

0

)" 2T62 26T2 %

& A.

;

)% / &

8.

;

0

)%

&

#

> ' ! :.

R2

2F '

. 1

) %

2

#

-2F #

UF

!

' # R2 2F

U(

' ' / 2

2 2

,

2

-2F U%

# 2.

$ L 0

2

$ 2%

& 2 $ "

RU% R U%

2%

) 2% ! 1 2

R $

RU%

U%

N

%

+

)

%

*

,

.

1 -

0 %

2%

+"

%

2

"

$ @!A D

#) 2%

U/

2

!

-K ' >

>

$ .

#

@*A

/

#

2% R 2

%

R D

#

2%

2

' @+A / @4A -

. N

$ 2V2

$ 2H2 "

2%

2

' 2%

2 @-A D

# @5A 2Q

# 2-

'

#

. #

@6A 2

1

2 (

#

2 2

%

. !. L

#

N 2

2

) N

OTO D & ! %

N

%

R

V- .

- .

-.

9

)

H- .

C @X

-.

W

, D

)

(( ) ;@379$ 3! 4 89A )7C@ G*,4 ,8A (( ((

@ A @ A ,

% *%

R +

&

$

L

$

OTO D & * %

%

$

R

) - .; - .? - . - .;(X 9- . - .;,X - .> - .;- - .5 - ..>( %

,

-.

-.

F

-. ,

( -.

( -.

-.

R

-

*

+.

$ -

#

#

.

$

$ ) *

Y +,

Y +,* - . *

/ !) #) +

* '! 8

0 *

! " # $ %$ & '(# ! (#) '& #

'!

9

!

'# 8

8

:

0 '

23 1 !)! ) 0

7

*

+ 4

5

( (

-

6

1

' %

R

$

*F

+-

F .

) Z - .;Z - .?Z - . Z - .;9X(XP XP- . Z - .;,X X - . Z - .;- - .5 - ..>(

*

)

Y

+,

+,* - . *

+

/ )!) ) #) ) +

*

'!

! : ;!; # *

: ! " # $ %$ & '(# ! (#) '& #

23 1 !)! ) 0

:

9

*

+ 4 ( (

-

5

6

1

3+5

3+4.5

6.5

++5

6.5

OTO D & + L

%

N

R

!

N

N

/ ,

,

%

? Y

% (

F

)

Y

Y

Y 7

*

+,

5

+,* - . *

/ )!) ) #) ) +

+

# ! 7

5

/ ;

=

!

! / /

;!;

/1

#

=) / = 1

/ ' ) /1 >

! " # $ %$ & '(# ! (#) '&

9 #

23 1 !)! ) 0

+ 4 ( (

-

5

6

1

!

<

% /

$ "

$

& -.

BE

" *%

+

$

'!

)

Y

Y

Y

Y 7

*

+,

5

+,* - . *

/ !) #) +

'!

0 ' '!

#

>

! 7 / /

5

/ ;

=

/1 =) /

/ ' ) /1 = > 1

0 *

9

( (

5

6

1 !

*

>

+ 4

-

-

'# 23 1 !)! ) 0

! " # $ %$ & '(# (#) '& !

#

++5

!

<

.

OTO

& !' # , DQP

2 &

$

/

$

)

%"&'

$ ( ) *

%"&'

$ ( ) * P -P

) 2

;

B

Q

A DQP

9

S- .

7 @ ( @B D @9 @7 \ A ( AB D A9 A7 \

S- . -

$

.

(' -

$ ('

%

R.

0 U

!,D"FS[

)

!"# $ " # (!" # )!"#

*

!"# "#

!" #

!"# L

2 =

#

OTO

& *'

! @ ? = 5

A

B

"

)

5L

-.

,, -.

5L

]

,,

]

)( ]

5 "

$ !-

5 $

D . "

^ $

-.

"

) - . ; ]X - . - . ; ( X - . ? X - .> ? - . - . ; ] X1- . -.;] X -. - . ; - . ? ^X 1- .> ? !X1- . K DQP

)

L - . ; ]XL - . L - . ; ( XP - . ? X XP - . ? L - . L - . ; ] X0- . / - . ; ] XP - . / - . ; / - . ? ^X X0- . ? !X0- . & )

$

(

)

9: *

$

+' # +' , -.,"/ 012 +3' # , - ,4/ 2 5' - 2 6' # ,% # & -*2 7' # -+8,4 2 ' % , & -+8,4 /"2

,

(

+' # ' +' , ' -.,"/ +3' # , ' - ,4/ 5' ' - 2 6' # ,% # &' 7' # ' -+8,4 ' % , &' -+8,4

% 1- .

43;:

R

-.

)

*

012 2 -*2 2 /"2

-.

P -D

$ A

B

Q

9 D

)

0 U

L

2 =

2

DQP

$

7 \

,

-] ]

.

@ %

L -. @B K @9 D

A

;

'

N @ AB D A9 K

$ '

8 % A -( 8B D $ 89 K 87 K

'

L -. A.

1;B

>

' ' '

OTO

& +' O , DQP

/

; &

$

$

P -P B

Q

!"# "#

!" #

!"#

=

2

;

A

S- .

7 K

S- . $

!

;A

A K S- . -4 8 M,

>

.

-

.

S- .

-

.

"Q,PKQ3 / ! *

P

.

+O

: K ];B

-.

S- .

-. V- . ?

S- .

] Y

C K W K S- . -"

*

DQP

9 P

@ D

#

)

2

U

" # (!" # )!"#

%"&'

$ ( ) *

)

!"# $

%"&'

$ ( ) *

)

!,D"FS[

-.

.

]; A

H- .

OTO

& 4' E

R

%

N

$

B EEE

9

$ ;E B

%

) -K L-.

Y

]X] SB- .;

?

%

0- .

^X X 9?-!X ?^X( .X ?-( X!?] X] .

) ) _E@ '

"

9.

X ?( S9- .;

X ?(

"

$

/- .

%

) -.

_ A`

'

) /- .

\- . ?

L- . (- .

Y

]X] SB- .;

Y

?

X ?( S9- .;

X ?(

0- .

^X X 9?-!X ?^X( .X ?-( X!?] X] .

] [- .;

P -D

$

B

)

0 U

A Q

X ?B

L

2 =

2

;

DQP

9 \

,

7 (

N

-] ]

. L N ) _E9

WEa

_B % @ D

$

%

A %

-. AB K A9 K A7 K

1- .

8 %

9 EEE '

' -.

$

R

SIMULINK – MATLAB CONTENIDO 1. ELEMENTOS BÁSICOS 2. EL MOTOR DC 3. SUBSISTEMAS 4. ECUACIONES DIFERENCIALES 5. SIMULACIÓN DE SISTEMAS

INTRODUCCIÓN Simulink es una extensión de Matlab utilizado en el modelamiento y simulación de sistemas. Para arrancar Simulink se puede hacer desde el prompt de Matlab digitando el comando >>Simulink o utilizando el icono

. Se abre la ventana

Simulink Library Browser como se indica abajo y se puede diagramar un nuevo modelo activando el botón New Model , o sea el icono

o de

1

Un modelo es un conjunto de bloques que representa un sistema y como archivo tiene extensión *.mdl

1. ELEMENTOS BÁSICOS Los elementos básicos son líneas y bloques. Los bloques están agrupados en: Sources, Links, Discrete, Continuos, Math, etc., tal como aparecen en la ventana anterior. Cada bloque tiene entradas y salida para realizar su interconexión. Por ejemplo, haga clic en Discrete y luego clic en Discrete Transfer Fcn y arrastre el bloque a la ventana en blanco. Si quiere modificar la función de transferencia del bloque haga doble clic en él y digite los coeficientes del numerador y denominador en la nueva ventana que aparece. Para la función 1/(z 2 +2z +4) con tiempo de muestreo de 1 seg, quedaría:

Para realizar el diagrama en bloques de un sistema se hace lo siguiente:

2

Lo primero es arrastrar los bloques a la página en blanco de forma que,

Step es

la función paso o escalón que se obtiene de Sources, Scope es el osciloscopio que se obtiene de Sinks, Transfer Fcn se obtiene de Continuos, Sum y Gain se

obtienen de Math. Modifique los bloques dando doble clic sobre cada uno de ellos para cambiar sus parámetros o valores e interconéctelos.

Lo segundo es cambiar los nombres a los bloques y asignar las variables o señales haciendo doble clic en el lugar en que se van a colocar y salvar el modelo especificándole un nombre, por ejemplo ejem1.mdl

Por último se debe simular el sistema. Para ello se configura la señal de entrada, en este caso la función paso. Dar doble clic y asignar los siguientes parámetros: Step time=0, Inicial value=0, Final value=1, Sample time=0. Para simular el sistema de control se escoge del menú

o el icono

.y luego

se hace doble clic en Scope para ver su respuesta o salida del sistema. Para observar además la entrada se puede colocar otro Scope a la salida de Step y se puede probar para varios pasos variando su amplitud, tiempo de inicio y tiempo de

3

iniciación del paso. Para observar mejor la respuesta se usa el botón Autoscale (binoculares

)

de la ventana del Scope. Si quiere observar mejor la

respuesta o parte de ella se pueden cambiar los parámetros de simulación, Simulation Simulation parameters. Por ejemplo cambiar el Start time y el Stop time y correr nuevamente la simulación.

2. EJEMPLO: MODELAR UN MOTOR DC Un actuador común en sistemas de control es el motor DC. Provee directamente movimiento rotatorio y acoplado con poleas o correas puede proveer movimiento transnacional.

2.1

ECUACIONES DINÁMICAS

El circuito eléctrico de la armadura y el diagrama de cuerpo libre del rotor es mostrado en la figura con sus ecuaciones dinámicas. 4

(1) Leyes de Newton

(2) Leyes de Kirchhoffs

Los parámetros físicos tienen los siguiente valores : Momento de inercia del rotor : J = 0.01kg.m2/sg2 Rata de amortiguamiento del sistema mecánico: b = 0.1 N.m.sg Constante de la fuerza electromotriz: Ke = Kt = 0.01 Nm/Amp Resistencia eléctrica: R = 1 ohm Inductancia eléctrica: L =0.5H Fuente de voltaje de entrada: V Posición angular: θ Se asume que el rotor y el eje son rígidos

5

2.2

MODELADO DEL MOTOR EN VELOCIDAD

2.3 EXTRAER MODELO LINEAL

Para obtener la función de transferencia del motor primero se trasladan los parámetros del motor al modelo creando un archivo en Matlab (*.m) de la siguiente forma:

% VALORES DE LOS PARÁMETROS DEL MOTOR J = 0.01; b = 0.1; Ke = 0.01; 6

Kt = 0.01; R = 1; L = 0.5;

Se ejecuta este archivo y se simula el modelo para una entrada de paso unitario de valor V = 0.01, con los siguientes parámetros de simulación: Stop time = 3 sg. Arranque la simulación y observe la salida (velocidad del motor).

Como segundo paso se debe obtener el modelo lineal de Matlab del motor. Para esto, borre el bloque Scope y cámbielo por Out obtenido de la librería de Signals&Systems. Haga lo mismo para Step cambiándolo por In de esta misma librería. Los bloques In y Out definen la entrada y salida del sistema que le gustaría extraer. Salve este modelo. El sistema quedará así:

Como tercero y último paso, después de desarrollado el modelo y salvarlo por ejemplo con el nombre MotorDcVel.mdl se ejecutan los siguientes comandos:

7

% OBTENER EL MODELO LINEAL DEL SISTEMA [num, den] = linmod('MotorDcVel') Gps = tf(num, den) La respuesta es :

3. SUBSISTEMAS Abra una nueva ventana y arrastre de la librería Signals&Systems el bloque SubSystem , haga doble clic en este bloque, abra el modelo MotorDcVel.mdl (el que tiene In y Out como terminales) cópielo y péguelo en la nueva ventana de subsistema anterior. Cierre ventanas y aparece una nueva con el bloque con los terminales del subsistema creado. Déle el nombre MotorDcVel. Si a este bloque de subsistema se le da doble clic aparece el modelo completo diseñado anteriormente. Otra forma es señalar los bloques de interés, ir a menú Edit --> create Subsytem

3.1 SISTEMA EN LAZO ABIERTO

Al subsistema creado que constituye la planta de un sistema de control se le va a adicionar un controlador y obtendremos la función de transferencia en lazo abierto y lazo cerrado. 8

% CONTROL DE UN MOTOR DC [num, den]=linmod('ControlMotor') Glazo_abierto = tf(num, den)

Respuesta:

3.2

SISTEMA EN LAZO CERRADO

% CONTROL DE UN MOTOR DC [num, den]=linmod('ControlMotor') Glazo_cerrado= tf(num, den)

Respuesta:

9

3.3

SISTEMA DISCRETO

DIAGRAMA EN SIMULINK

PROGRAMA MATLAB

% SISTEMA DISCRETO DISCRETO T=0.1; [num,den]=dlinmod('MotorDigital',T) Glazo_cerradoz=tf(num,den,T)

Respuesta:

10

4. SOLUCIÓN DE ECUACIONES DIFERENCIALES Ejemplo:

Resolver la siguiente ecuación diferencial:

d2y dy  3  6 y  4t  y' '3 y'6 y  4t  y' '  4t  3 y '6 y 2 dt dt Diagrama Simulink:

Respuestas:

11

Ejemplo:

Comprobar la integración por Simulink.

12

5. SIMULACIÓN DE SISTEMAS 5.1 INTERCAMBIO DE MATLAB A SIMULINK

Para utilizar señale de Matlab a Simulink de la librerís Sources se utiliza el bloque From Workspace.

13

Ejemplo: Resolver la ecuación y’’ + y = e t, y’(0) = 0, y(0) =3

El vector [ t x ] se ejecuta en Matlab en el workspace de la siguiente forma:

>> t = 0:0.001:0.999; >> t = t’; >> x = exp(t)

Al ejecutarse Simulink toma los datos entregados por Matlab. No olvidar colocar condición inicial y(0) = 3 en el integrador. 5.2 INTERCAMBIO DE SIMULINK A MATLAB

Para enviar datos de Simulink a Matlab se utiliza de la librería Sinks el bloque To Workspace.

Ejemplo: Resolver la ecuación: f(t) = Mx’’ + Bx’ + Kx,

M=1, B=1, K= 10, F(t) = 5

14

Diagrama Simulink:

En Matlab: >> plot(t,y)

5.3 EJERCICIOS

Ejercicio1:

15

Si la entrada es una señal senoidal, encontrar las salidas referidas a vC y iL.

Ejercicio2:

Para el siguiente problema hallar la variación de h si el caudal normal Q es de 10 lit/min y en t=5 seg se aplica una perturbación de 2 lit/min. El valor de K=10, A= 2 m2.

A

dh  q(t )  K h dt

Diagrama Simulink:

16

EJERCICIO3: LA BOLA MAGNÉTICA

Ecuaciones:

(1)

m

d 2h i 2  mg  h dt 2

(2) L

di  V  iR dt

Valores:

m=0.1 Kg;

g=9.81;

R=2 Ohm;

L=0.02 H;

=0.001

Diagrama simulink:

17

Controlador:

zeros=[-11.5+7.9i, -11.5-7.9i] polos=[0 -1000] ganancia=-3.3057e+004

Planta:

i(0) = 0;

h(0)=0.05;

h’(0)=0

18

19

EJERCICIO4: TANQUE DE AGUA

Ecuación del modelo:

dVol dh A  bV  a h dt dt

Diagrama simulink:

20

Controlador:

Planta:

EJERCICIO5: MOVIMIENTO PARABÓLICO

21

Ecuaciones:

x' '  0 Movimiento unforme y ' '   g Movimiento acelerado

Condiciones iniciales:

Vo=100 m/sg;

 = 30º

22

23

EJERCICIO6: PÉNDULO SIMPLE

Ecuación:

mL ' ' BL 'wsen  0

Valores:

w (peso) = 2;

L (longitud) = 0.6;

B (amortiguación) = 0.08;

Condiciones iniciales: ’(0) = -2 rad/sg; (0) =  /2

24

Diagrama simulink:

25

EJEMPLO: SISTEMA MECANICO

Parámetros:

m1=40; m2=60; k1=400; k2=400; b1=180; b2=220;

Ecuaciones dinámicas:

f (t )  m1

0  m2

dv1  k1 (v1  v2)dt  (v1  v2)b1 dt

dv2  k1 (v2  v1)dt  (v2  v1)b1  k 2 v2dt  b2v2 dt

Ecuaciones de Laplace:

F ( s)  m1sV1 

0  m2sV 2 

K1 (V 1  V 2)  (V 1  V 2)b1 s

k1 k2 (V 2  V 1)  (V 2  V 1)b1  V 2  b2V 2 s s

26

Diagrama simulink:

EJEMPLO: SISTEMA TERMOQUÍMICO

Se desarrolla una reacción termoquímica en donde el reaccionante A se convierte en un producto B.

Velocidad de reacción: r(t)= k c(t) Constante de velocidad de reacción: k = 0,2 min -1

27

Concentración de la entrada: ci(t) Para t= 0; ci(0)=1.25 lbmol/pie3 Volumen de la masa reaccionante: V= 5 litros Flujo de entrada: F= 1 lt/min

Ecuación dinámica:

V

dc(t )  Fci (t )  Fc (t )  KVc (t ) dt

V

dc(t )  Fci (t )  ( F  KV )c(t ) dt

V dc(t ) F  c(t )  ci (t ) F  KV dt F  KV

Constante de tiempo:

 

V F  KV

Ganancia de estado estacionario:

Ke 

F F  KV

Reemplazando valores:  = 2.5 min; Ke = 0.5; Condición inicial de la concentración: c(0)

0 = Fci(0) - Fc(0) - KVc(0) Reemplazando valores: c(0) = 0.625 lbmol/pie 3 28

Programa en Matlab:

%Entrada al paso. Programa pplineal.m function dy=pplineal(t,y) global K X tau dy=(K*X-y)/tau;

% Entrada rampa. Programa rplineal.m function dy=rplineal(t,y) global K r tau dy=(K*r*t-y)/tau;

% Entrada senoidal. Programa splineal.m function dy=splineal(t,y) global K tau A w dy=(K*A*sin(w*t)-y)/tau;

% Programa principal F=1; V=5; K=0.2; ci0=1.25; c0=solve('F*ci0-F*c0-K*V*c0=0'); c0=eval(c0) %Constante de tiempo tau=V/(F+K*V) % tau=2.5 minutos %Ganancia en estado estacionario Ke=F/(F+K*V) % Ke=0.5 29

global R K tau X r A w Rango Inicio Rango=input('Tiempo de simulacion='); Inicio=input('Condiciones iniciales='); N=input('ESCRIBA 1=PASO, 2=RAMPA, 3=SENO: '); disp(' ') switch N case 1 X=input('Valor del paso='); [t,y]=ode45('pplineal',Rango,Inicio); plot(t,y)

case 2 r=input('valor pendiente de la rampa='); [t,y]=ode45('rplineal',Rango,Inicio); plot(t,r*t,t,y/K,'r')

case 3 A=input('Amplitud del seno='); w=input('Frecuencia del seno='); [t,y]=ode45('splineal',Rango,Inicio); disp('Amplitud del perfil de la respuesta') K*A/sqrt(1+(w*tau)^2) disp('Fase de la respuesta respecto a la entrada') atan(-w*tau) plot(t,A*sin(w*t),t,y,'r') end

30

Programa en Simulink:

EJEMPLO: SISTEMA HIDRAULICO

% HIDRAULICO UNA ETAPA C1=3; R1=1; C2=10; R2=2; qi=2; keyboard plot(t,qo) title('HIDRAULICO') grid 31

pause n=1; while n==1 T=input('Entre tiempo: ') delta=input('Entre valor delta: ') i=find(t<=(T+delta)&t>=(T-delta)); tiempo=t(i) caudal_salida=qo(i) n=input('Entre 1 para seguir y 0 para parar: ') end

EJEMPLO: SISTEMA ELÉCTRICO

% CIRCUITO RC DE DOS ETAPAS R1=10e3; R2=20e3; C1=1e-6; C2=10e-6; ei=10; keyboard plot(t,eo) title('CIRCUITO RC') grid 32

pause n=1; while n==1 T=input('Entre tiempo: ') delta=input('Entre valor de delta: ') i=find(t<=(T+delta)&t>=(T-delta)); tiempo=t(i) voltaje_salida=eo(i) n=input('Entre 1 para seguir y 0 para parar: ') end

33

INTRODUCCION A SIMULINK Matlab (Matrix Laboratory) es un sistema basado en matrices para realizar cálculos matemáticos y de ingeniería. Entre las múltiples herramientas que presenta este programa se encuentra Simulink que es una librería de MATLAB que permite la simulación de procesos mediante diagramas de bloques.

1.

Acceso a la librería de bloques de Simulink: Para acceder a la librería de Simulink se debe abrir inicialmente la ventana principal de Matlab (Matlab Command Window). En esta se puede ejecutar el comando “simulink” o hacer clic en el símbolo correspondiente en la barra de herramientas en la parte superior de esta ventana. Al hacer esto aparecerá el listado de las librerías correspondientes a simulink, donde se podrá tener acceso a todos los bloques que brinda esta herramienta. Para abrir una nueva hoja de trabajo se deberá acceder a través de: File à New à Model, o hacer clic en el símbolo de “hoja nueva”.

Figura 1. Simulink Library Browser

La librería “Simulink” contiene los bloques necesarios para simular un sistema mediante técnicas convencionales, las demás librerías son herramientas adicionales que se utilizan para aplicaciones específicas de control avanzado. En la Figura 2 se muestra el contenido de la librería Simulink.

Figura 2. Librería Simulink. En la figura 3 se pueden observar los bloques mas utilizados en la simulación de procesos. Estos se encuentran en hacer clic en el signo (+) de cada librería. Pueden ser utilizados al hacer clic sobre ellos y arrastrándolos sobre la hoja de trabajo. Igualmente en la tabla 1 se pueden observar una descripción mas detalladas de algunas de estas funciones.

Figura 3. Bloques mas utilizados en la librería Simulink

Tabla 1. Detalles de algunos bloques de la librería Simulink. Bloque

Función

Librería

Parámetros requeridos

Constante

Asigna un valor constante a la entrada.

Sources

Valor de la constante.

Entrada escalón

Introduce un escalón de magnitud específica en un tiempo dado.

Sources

Tiempo del escalón Valor inicial del escalón, Valor Final del Escalón

Entrada Rampa

Introduce una rampa en un tiempo especificado.

Sources

Tiempo de la rampa, pendiente

Entrada Senoidal

Introduce una señal senoidal específicada por el usuario.

Sources

Amplitud de la onda, Fase.

Workspace

Almacena datos de la señal que llega al bloque y la convierte en vector. Si se conecta al reloj se almacena el vector tiempo.

Sinks

Nombre y tipo de la variable (Save format: Matrix)

Scope

Grafica la señal que se introduzca con respecto al tiempo.

Sinks

XYgraph

Grafica la entrada superior en el eje x y la inferior en el eje y

Función de Transferencia

Representa la función transferencia a lazo abierto.

de

Contin.

Numerador y Denominador de la FT

Integrador

Integra una señal en función del tiempo

Contin.

Valor inicial desde el cual se va a integrar

Retardo de transporte

Introduce un retardo en el tiempo en el cual aparece la señal.

Contin.

Valor del retardo (Debe ser un número positivo)

Ganancia

Multiplica la señal por cualquier valor de ganancia que se introduzca.

Math

Valor de la ganancia

Sumador

Suma dos señales.

Math

Número de entradas a sumar

Multiplicador

Multiplica dos señales.

Math

Número de entradas a multiplicar

Matlab-Function

Aplica cualquier función matemática conocida por Matlab a la señal.

Functions and Tables

Función a utilizar

PID

Es un controlador donde se puede introducir una parte proporcional, una integral y una derivativa

Blocksets & Toolbox.:

Proporcional: K, Integral: K/Ti, Derivativo: K.Td

Sinks

Entradas

Rango de los ejes

Simulink-extras: Aditional linear Mux

Permite representar dos señales distintas en una misma gráfica.

Signals &System.

Número de entradas

Obs. Colocando el nombre del bloque de interés en el buscador (Simulink library browser), se puede ubicar directamente en la librería de bloques.

Notas sobre el uso de Matlab.

-

Los bloques pueden ser movidos al arrastrase con el botón izquierdo del mouse y pueden ser copiados al hacer clic sobre ellos con el botón izquierdo del mouse y arrastrando la copia creada.

-

Los bloques se deben unir mediante flechas. Esto se logra haciendo clic en la flecha de salida del bloque deseado y conectándola (sin soltar el botón), a la flecha del bloque que se desea unir.

-

Para que las modificaciones en el programa hagan efecto este deberá ser grabado después de realizar los cambios.

-

Al hacer clic con el botón derecho sobre una señal, se podrá obtener una “línea” de esta señal para llevarla o conectarla a un bloque deseado.

-

Al barrer el mouse sobre un grupo de bloques se podrán mover estos a la vez y copiarlos en grupo.

-

El “save format” de los “workspace” debe ser colocado en “array” para poder “graficar” las variables deseadas posteriormente.

2. Ejemplos de Uso de Simulink

Modelos Matemáticos no Linealizados. a)

Simular la siguiente ecuación diferencial y encontrar su respuesta ante una entrada escalón.

d ( x)  − 60  + 5 exp  * X (t ) = 1 dt  T 

Donde, X0 = 0 en T=500

Solución:

-

Seleccionar los bloques necesarios para representar el modelo y llevarlos a la hoja de trabajo. Para ello debe buscar los bloques en las librerías correspondientes (Ver tabla 1 y figuras 2 y 3), seleccionar cada uno haciendo "click" sobre él para marcarlo y arrastrarlo con el "mouse" hasta la ventana.

Figura 4. Bloques necesarios para la representación de la ecuación.

-

Armar el modelo. Las ecuaciones diferenciales pueden representarse en bloques de "Simulink" como función del tiempo sin linealizarlas ni llevarlas al dominio de Laplace. En primer lugar, se debe despejar la derivada temporal para expresarla en función de los demás términos de la ecuación:

d ( x)  − 60  = 1 − 5 exp  * X (t ) dt  T 

(1)

Luego, se debe establecer qué valores en la ecuación son constantes y cuáles son función del tiempo, en este caso, la conversión X depende del tiempo, pero la temperatura T es constante. El primer término del lado derecho de la ecuación (1) se puede representar como un escalón unitario o como una entrada constante. En el segundo término [5exp(-60/T).X(t)], la temperatura se representa, igualmente, como una entrada escalón o como una entrada constante, luego se invierte con un bloque Matlab Function, donde se especifica la función 1/u, se multiplica por una ganancia de –60 con un bloque Gain, y se introduce nuevamente en una Matlab Función para obtener la exponencial, que va a ser multiplicada por 5 con otro bloque Gain, como se muestra en la figura 5. La variable X(t) no se conoce porque es el resultado de integrar el lado derecho de la ecuación, esta variable debe multiplicarse con la exponencial con un bloque producto para formar el segundo término de la ecuación diferencial, luego ambos términos se combinan con un bloque Sum para obtener la ecuación (dX(t)/dt) completa que pasa por un integrador para obtener la variable X(t), que se realimenta al bloque producto. Se debe colocar igualmente un bloque de reloj para que el simulador contabilice el tiempo. Todas las señales que se deseen guardar o ser posteriormente llamadas para graficar, deben ser alimentadas a un bloque workspace (Save format : Array)

Figura 5. Diagrama de bloques de la simulación de la ecuación diferencial. Haciendo clic en cada bloque se pueden cambiar sus parámetros y sus nombres. En este caso se colocan los siguientes: Entrada Escalón (T)

Integrador

Matlab Function (1/T)

Matlab Function Exp(E/KT)

Step Time = 0 Valor Inicial = 500 Valor Final = 500

Inicial Value = 0

Fuction = 1/u

Fuction = exp

De esta forma tenemos:

Figura 6. Diagrama de Bloques con sus parámetros.

-

A continuación se abre el menú Simulation à Simulation Parameters y se modifica el tiempo de parada. En este caso se pondrá 8 seg. Posteriormente se simula la ecuación diferencial apretando el botón o símbolo de “play” en la parte superior de la pantalla y se espera a que la maquina realice el calculo (indicado en la parte inferior derecha de la ventana)

-

Después de realizar la simulación se regresa a la ventana principal de Matlab y se grafican los resultados colocando el comando: plot(T,X) y se obtiene:

Figura 7. Simulación de la respuesta de la ecuación diferencial.

Para observar la respuesta del sistema ante una perturbación se coloca:

Entrada Escalón (T)


Para que el simulador acepte el cambio, el archivo debe ser guardado (no debe aparecer un “asterisco” al lado del nombre del programa en la parte superior de la ventana). Después de simular y graficar, se obtiene:

Figura 8. Respuesta del sistema de la ecuación diferencial ante entrada escalón. Se observa la perturbación del sistema en el tiempo = 4 seg.

Ecuaciones Diferenciales Acopladas b) Simular el siguiente sistema de ecuaciones diferenciales.

dT = A(T − To1) * Fo − B(T − Ta) dt

A = 0.08

B = 0.05

dTa = C (Ta − To 2) * Fp + D(Ta − T ) dt

C = 0.5

D = 0.01

Donde Fo, To1 y To son posibles perturbaciones cuyos valores son: Fo = 0.01 m3/min, To1 = 280K, To2 = 350K

Solución:

-

Los Bloques seleccionados son:

Figura 9. Bloques usados para la representación del sistema.

-

Armar el modelo. Las ecuaciones diferenciales se "escriben en bloques de Simulink" como función del tiempo sin realizar ninguna modificación de las mismas para linealizarlas o llevarlas al dominio de Laplace.

Figura 10. Diagrama de bloques del proceso con sus parámetros.

Los parámetros iniciales de los bloques y simulación serán: Entrada Escalón (T01)

Entrada Escalón (T02)

Entrada Escalón (Fo)

Integrador (T)

Integrador (Ta)

Stop Time



Step Time = 0 Valor Inicial = 0.01 Valor Final = 0.01

Initial Value = 280

Initial Value = 350

200

-

Después de Simular y Graficar usando: Plot(t,[T,Ta]) se obtiene,

Figura 11. Simulación del sistema de ecuaciones diferenciales

Funciones de transferencia a lazo abierto y lazo cerrado c)

Dada la siguiente función de trasferencia obtenga la respuesta del sistema a lazo abierto ante una entrada escalón.

FT =

1 s + 4s + 2 2

Solución:

-



-

Armar el modelo.


Los parámetros iniciales de los bloques y simulación serán:

Entrada Escalón

Función de Trasferencia

Parámetros Simulación


Numerador = [1] Denominador = [1 4 2]

Stop Time = 20

-

Después de Simular y Graficar usando: Plot(T,X) se obtiene,

Figura 14. Respuesta temporal de la función de transferencia.

d)

Compare la respuesta a lazo cerrado obtenida para la función de transferencia del ejemplo anterior con la respuesta obtenida al aumentar 5 veces la ganancia y con la respuesta obtenida al usar un controlador de función de transferencia:

1  Gc1 = K 1 +   Ti 

Donde, K = 1

y Ti = 0.5

Solución:

-



-

Armar el modelo.


Los parámetros de los bloques son:

Entrada Escalón

Ganancia

Mux

PID

Stop Time


Gain = 5

Number of inputs =3

Proporcional = 1 Integral = 1/0.5 Derivative = 0

20

-

Después de Simular y Graficar usando: Plot(T,s) se obtiene,

Figura 17. Respuesta temporal de la función de transferencia.

Tutorial Simulink 2

Recommend Documents