GRUPO DE PROCESADO DIGITAL DE SEÑALES
PRÁCTICAS DE INTRODUCCIÓN AL PROCESADO DIGITAL DE SEÑALES.
Marcelino Martínez / Emilio Soria / José V. Francés / Curso 2009-2010
[email protected] [email protected] [email protected]
LABORATORIO: Práctica 1.- Introducción al programa MATLAB. Práctica 2.- Muestreo de Señales. Práctica 3.- Cuantificación. Práctica 4.- Convolución y ecuaciones en diferencias Práctica 5.- Sesión 1. Mini-proyecto de IPDS Práctica 6.- La transformada Z Práctica 7.- Sesión 2. Mini-proyecto de IPDS Práctica 8.- Respuesta en frecuencia de los sistemas discretos Práctica 9.- Sesión 3. Mini-proyecto de IPDS Práctica 10.- Práctica-resumen. Cada una de las sesiones de laboratorio tiene una duración de 3 horas. OBJETIVOS: Los objetivos del módulo de laboratorio de Introducción al procesado digital de señales son que el alumno lleve a la práctica los conceptos adquiridos en el módulo teórico, implementado algoritmos y construyendo sistemas de procesado que sirvan de nexo de unión con la infinidad de aplicaciones prácticas que esta disciplina tiene en la actualidad. BIBLIOGRAFÍA BÁSICA: • • • • • •
"Tratamiento de Señales en Tiempo Discreto".Oppenheim, A.V,; Schafer, R.W.;Buck J.R.Prentice Hall, 2000. ISBN: 84-205-2987-7 "Tratamiento Digital de Señales: Principios, Algorimtos y Aplicaciones".Proakis, J.G.; Manolakis, D.G. Prentice Hall, 1998. ISBN:84-8322-000-8 "Introduction to Signal Processing". Orfanidis, S.J. Prentice-Hall. 1996. ISBN: 0-13-240334-X “Señales y Sistemas". Oppenheim, A.V.; Willsky, A.S.; Young, I.T.Prentice-Hall. 1997. ISBN: 970-170116-X "Tratamiento Digital de Señales: Problemas y Ejercicios Resueltos". Soria, E., Martínez M., Francés, JV., Camps, G., Pearson Prentice-Hall. 2003. ISBN: 84-205-3559-1X “Digital Signal Processing Using Matlab”. Ingle, V., Proakis J.G. Brooks/Cole Thomson Learning. 2000. ISBN:0-534-37174-4
MÉTODO DE EVALUACIÓN: Módulo de laboratorio: Se realizará un examen, para los alumnos que asistan al laboratorio, en la sesión número 10, cuyo valor será el 50% de la nota de laboratorio. El 50% restante se obtendrá de la evaluación continua de cada una de las sesiones (10 %) y de la realización de un mini-proyecto (40%). Habrá un examen de laboratorio, el día de la convocatoria oficial, al que podrán presentarse los alumnos que no hayan asistido a las sesiones de laboratorio o que deseen mejorar su nota que constituirá el 100% de la nota de laboratorio.
NOTA: Es imprescindible la preparación de las prácticas antes de acudir al laboratorio.
2
PRÁCTICA I: INTRODUCCIÓN AL MATLAB. Esta primera sesión, a diferencia de las demás, es una práctica completamente guiada para que el alumno se familiarice con el programa Matlab que se utilizará para el resto de prácticas. En primer lugar es conveniente revisar alguno de los manuales sobre este programa que se encuentran en la página web de la asignatura (www.uv.es/~martsobm) Desde el Administrador de archivos de Windows crea una carpeta en d:\IPDS\minombre en la que se almacenarán todos los programas que escribamos. Ej d:\IPDS\Martínez Al iniciar Matlab veremos la ventana de comando y el prompt (>>). Matlab está esperando a que se introduzcan órdenes.
Desde el prompt Matlab funciona como una calculadora realizando las operaciones que se le indique y mostrando el resultado, si introducimos 3+4*8 nos muestra el resultado, 35 que se almacena en una variable generada automáticamente llamada ans que contiene el resultado (answer) de la última operación. >> 3+4*8 ans = 35 Además de estas operaciones básicas, Matlab tiene infinidad de funciones y comandos que permiten realizar tareas específicas de todo tipo. Estas funciones se encuentran agrupadas formando librerías (Toolboxes). De las mútiples librerías existentes nos centraremos en la librería de Signal Processing. Un ejemplo de comandos de Matlab son pwd y cd-. Podemos obtener información de cualquier comando utilizando la instrucción help. Así si escribimos en el prompt >> help pwd PWD Show (print) current working directory. PWD displays the current working directory. S = PWD returns the current directory in the string S. See also CD. Obtenemos información de este comando, que nos devuelve el directorio actual, y nos sugiere otros comandos relacionados. Si lo ejecutamos >> pwd ans = 3
c:\users\marsel •
Utiliza el comando chdir para situarte en el directorio que has creado anteriormente
>>cd(‘c:\ipds\martinez’) •
Verifique el directorio actual con el comando pwd.
Matlab (MATrix LABoratory) está optimizado para trabajar con matrices de datos por lo que prestaremos especial atención a este tipo de datos y sus operaciones. Definicion de variables. Podemos definir variables asignándoles un valor: >> x=5, y=7; x = 5 Cuando una asignación acaba con ‘;’ no se muestra el resultado por pantalla pero la variable esta definida y puede utilizarse en operaciones posteriores >> z=x+y z = 12 Hay una serie de variables que Matlab tiene definidas por defecto como son: pi, j, i, inf(valor máximo representable) , NaN (el resultado no es un número ej 0/0), etc. Estas variables conservan su valor original si el usuario no las modifica; es decir si escribimos >> pi ans = 3.1416 Pero si realizamos una nueva asignación >> pi=33 pi = 33 dejará de tener su valor original. Las variables i, j originalmente representan a la unidad imaginaria Para declarar vectores procedemos de igual forma pero los elementos de un vector se escribirán entre corchetes separados por comas o espacios en blanco para un vector fila y separados por ‘;’ para un vector columna >> x1=[1 3,4], x2=[3;5;7] x1 = 1
3
4
x2 = 3 4
5 7 Se puede transformar un vector fila en uno columna y viceversa utilizando el operador traspuesta “ ’ ”. >> y1=x1' y1 = 1 3 4 >> y2=x2' y2 = 3
5
7
Se puede acceder a un determinado elemento de un vector con el nombre del mismo e indicando entre paréntesis el índice del elemento al que accedemos >> x2(2) ans = 5 EL ÍNDICE DEL PRIMER ELEMENTO DE UN VECTOR ES EL 1. SI SE INTENTA ACCEDER A UN ELEMENTO CON UN INDICE 0 O NEGATIVO SE PRODUCE UN ERROR.(no confundir con el lenguaje C) >> x2(0) ??? Subscript indices must either be real positive integers or logicals. Para definir una Matriz procedemos de forma similar a un vector, separando los elementos de una fila por comas o espacios en blanco y las columnas ‘ ; ’ >> A=[1 2 5;3,6,4;1+j,3-j 1] A = 1.0000 2.0000 5.0000 3.0000 6.0000 4.0000 1.0000 + 1.0000i 3.0000 - 1.0000i 1.0000 Vemos también como se pueden definir elementos complejos, con los que Matlab opera haciendo operaciones exactamente igual que con los reales. Para acceder a los elementos de una matriz procedemos de la misma forma que para un vector pero en este caso hemos de indicar con un índice el número de fila y con otro el de columna. >> A(3,2) ans = 3.0000 - 1.0000i >> b=A(2,2) b = 6 5
También se pueden generar matrices uniendo vectores >> P=[x1;x2';[7 9 5]] P = 1 3 4 3 5 7 7 9 5 Obsérvese que como x2 lo habíamos definido como un vector columna será necesario trasponerlo, de lo contrario se producirá un error >>P=[x1;x2;[7 9 5]] ??? Error using ==> vertcat All rows in the bracketed expression must have the same number of columns. Se puede eliminar un elemento de un vector asignándolo a la matriz vacía [], >> y1 y1 = 1 3 4 >> y1(2)=[] y1 = 1 4 El elemento indicado se elimina y desplaza al resto. Se pueden realizar operaciones aritméticas sobre los datos con los operadores +, -, *, /, ^ (potencia), para matrices además se puede calcular la inversa (inv), determinante (det), etc. Cuando estas operaciones se realizan sobre matrices es necesario que las dimensiones sean las correctas >> x1*x2 ans = 46 >> x2*x1 ans = 3 5 7
9 15 21
12 20 28
>> x1+x2 ??? Error using ==> + Matrix dimensions must agree. La división entre matrices está definida como el producto por la inversa. >> B=A' 6
B = 1.0000 2.0000 5.0000
3.0000 6.0000 4.0000
1.0000 - 1.0000i 3.0000 + 1.0000i 1.0000
>> A*B ans = 30.0000 35.0000 12.0000 - 1.0000i
35.0000 61.0000 25.0000 - 3.0000i
12.0000 + 1.0000i 25.0000 + 3.0000i 13.0000
0.2545 + 0.3091i 0.1000 - 0.1545i -0.0909
-0.2000 - 0.6000i 0.1000 + 0.3000i 0
0.4636 - 0.0636i 0.8364 - 0.4000i 1.5818 + 0.9273i
0.1000 + 0.3000i 0.2000 + 0.6000i -0.6000 - 1.8000i
0.4636 - 0.0636i 0.8364 - 0.4000i 1.5818 + 0.9273i
0.1000 + 0.3000i 0.2000 + 0.6000i -0.6000 - 1.8000i
>> inv(A) ans = -0.1636 - 0.1273i -0.1000 + 0.0636i 0.2727 >> det(A) ans = 11.0000 -33.0000i >> B/A ans = -0.1909 - 0.2091i -0.1091 + 0.4000i -0.9455 - 0.3818i >> B*inv(A) ans = -0.1909 - 0.2091i -0.1091 + 0.4000i -0.9455 - 0.3818i
En ocasiones nos interesa realizar operaciones entre elementos o matrices pero no como están definidos habitualmente sino indicando que queremos que se realicen elemento a elemento esto se lleva a cabo precediendo el operador por un punto así tendremos los operadores ./, .*, .^ >> A.^2 ans = 1.0000 9.0000 0 + 2.0000i
4.0000 36.0000 8.0000 - 6.0000i 7
25.0000 16.0000 1.0000
>> A^2 ans = 12.0000 + 5.0000i 29.0000 - 5.0000i 18.0000 25.0000 + 4.0000i 54.0000 - 4.0000i 43.0000 11.0000 - 1.0000i 23.0000 - 5.0000i 18.0000 + 1.0000i En la primera operación elevamos cada uno de los elementos al cuadrado mientras que en la segunda estamos calculando A*A. El Operador “:” Un operador muy empleado en los programas es “:”, que genera un secuencia de valores entre el valor inicial y el final con el incremento especificado. Su formado es, inicio:incremento:fin. Si no se especifica el incremento se considera de valor 1. El incremento no necesita ser entero, >> pares=0:2:9 pares = 0
2
4
6
8
7
9
>> impares=1:2:10 impares = 1 3 >> x=-1:0.25:1
5
x = -1.0000 -0.7500 -0.5000 0.5000 0.7500 1.0000
-0.2500
0
0.2500
También se puede utilizar este comando para acceder a los elementos de una matriz. En este caso si únicamente se indica el operador “:” indica que se toman todos los elementos de la fila o columna >> A(:,1) ans = 1.0000 3.0000 1.0000 + 1.0000i >> A(1:2,2:2) ans =
•
2 6 Verifique qué operaciones realizan los siguientes comandos:
>> B(:,1)=[1;9;2]; >>B(1:2:3,:)=[] 8
>>B=magic(4) >>B([1 3],1) Todas las variables que hemos ido declarando se almacenan en el espacio de trabajo de Matlab (workspace) la memoria del ordenador. Podemos ver las variables definidas utilizando los comandos who y whos, cuya única diferencia está en la información que proporcionan >> whos Name
Size
A (complex) B K P Tam ans b impares inpares k n numPages pares pi x x1 x2 y y1 y2 z
Bytes Class 3x3
144
2x4 1x7 3x3 1x1 3x1 1x1 1x5 1x5 1x7 1x11 1x1 1x5 3x3 1x9 1x3 3x1 1x1 2x1 1x3 1x1
64 56 72 8 24 8 40 40 56 88 8 40 72 72 24 24 8 16 24 8
double double double double double double double double double double double double double double double double double double double double
double array array array array array array array array array array array array array array array array array array array array array
Las variables pueden almacenarse en un fichero con la instrucción save >> save variables Las variables se pueden borrar con la instrucción clear. La sentencia siguiente borra las varibles indicadas >>clear x y z Si queremos borrar todas las variables >>clear Si visualizamos las variables, observamos que no hay ninguna >>whos Podemos recuperar las variables almacenadas con la instrucción load >>load variables >>whos Comprobamos que las variables se encuentran de nuevo en el espacio de entorno. LAS VARIBLES TIENEN EL MISMO NOMBRE QUE TENÍAN CUANDO SE ALMACENARON 9
Scripts. Hasta ahora nos hemos limitado a introducir órdenes en la línea de comandos, que no resulta operativo cuando hemos de realizar un gran número de operaciones en determinado orden. La alternativa es escribir un programa, que Matlab se denomina script, y que no son más que una secuencia de comandos almacenados en un fichero de texto que al ser llamado desde Matlab se ejecutará línea a línea en el orden que aparece en el programa. . Los ficheros de texto se pueden crear con cualquier editor de texto, si bien, Matlab lleva su propio editor integrado. •
Abra un nuevo script (File --> New -->Mfile) y escriba el siguiente código
%Mi primer script. Esto es un comentario %Vamos a realizar un programa para representar gráficamente las %ecuaciones y=x, y=2x e y=x/3+0.5 clear close all
%es conveniente borrar todas las variables al empezar %Cierra todas las ventanas gráficas
x=0:0.1:1; %Definimos un vector de valores de x y=x; %Calculamos cada una de las ecuaciones y1=2*x; %Cuando hacemos operaciones sobre vectores Matlab las realiza %automáticamente sobre todos los elementos no es necesario hacerlo %elemento a elemento como ocurre en C y2=x/3+0.5 clf %borra la pantalla grafica plot(x,y,'r') pause %Espera a que pulsemos una tecla para continuar plot(x,y,'r',x,y1,'g',x,y3,'b') pause plot(x,y,'ro',x,y1,'g+',x,y3,'b:') title('Graficas de y, y1 e y2') xlabel('x') ylabel('Amplitud') grid
Almacene el fichero en el directorio que creó al iniciar la sesión con el nombre primero.m (todos los scripts de Matlab tiene la extensión *.m) Compruebe que se encuentra en el directorio en el que ha almacenado el fichero con la instrucción pwd Para ejecutar dicho fichero escriba en el prompt de Matlab el nombre del fichero sin extensión >> primero
• • •
Corrija el error que se ha producido y vuelva a ejecutar el programa Haga la misma representación en el intervalo [-2:2] a con una distancia entre muestra de 0.1 Muestra por pantalla las muestras pares de y2.
>>y2(2:2:end)
•
%También se puede usar y2(2:2:length(y2))
Dibuje las muestras pares de y2
>> plot(x(2:2:end),y2(2:2:end)) Otra forma de superponer 2 gráficas en una misma figura, como se ha hecho en el ejemplo anterior, es utilizando el comando hold. 10
•
Copie las siguientes líneas de código en el programa anterior
figure %crea una nueva ventana gráfica plot(x,y,'r--') hold on plot(x,y1,'g.-') plot(x,y2,'b.') grid on %dibuja una malla de fondo hold off legend('y','y1','y2') xlabel('x') ylabel('Amplitud') print figura1 -djpeg %Almacena la figura como imagen en el fichero %figura1.jpg
En los programas anteriores hemos aprendido a poner leyendas en los ejes y títulos con las instrucciones xlabel, ylabel y title. Es posible representar varias gráficas en una misma figura dividiendo la figura en diferentes subventanas para ellos se utiliza el comando subplot(XYN). X: número de divisiones horizontales de la figura. Y: número de divisiones verticales de la figura N: número de gráfica dentro de la figura empezando a contar desde la esquina superior izda. Veamos que esto es muy sencillo. Las siguientes figuras muestran los parámetros de llamada del comando subplot para cada una de ellas.
11
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6 0.5
0.6 121
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.5
0 0
0.5
1
221
0.5
0 0
0 0 1
0.5
1
1
0.5
0.5
0 0
122
0.5
0.4
0 0 1
1
0 0
0.5
0.5
1
0 0
0.2
0.4
0.6
0.8
1
0.6
0.8
1
0.6
0.8
1
1
212
0.5
0.5
0 0 1
1
222
0.5
0.2
0.4
211
0.5
0 0
1
1
223
211
0.2
0.4
1
224
0.5
0.5
0 0
1
1
223
0.5
0.5
1
0 0
224
0.5
1
En primer lugar especificaremos el suplot en el que queremos dibujar y posteriormente llamaremos a alguna de las funciones de dibujo como plot. •
Copie el siguiente código que le permite hacer una representación en varias ventanas:
%Genere una sinuoside de periodo 8 con 32 muestras N=32 n=0:N-1 s=sin(2*pi*n/8) %Representela con plot y stem subplot(211), plot(n,s) subplot(212), stem(n,s) %Otro tipo de representacion muy usada en PDS disp('Pulse una tecla...') %Saca un mensaje de texto por pantalla pause %Dibujemos el espectro S=fft(s); %Calcula la Transformada de Fourier H=abs(S); %Calcula el modulo de un numero complejo (valor absoluto si es real) fre=2*n/N; clf %Borra la ventana grafica activa subplot(211); plot(fre,H) %Modulo title('Modulo en escala lineal') xlabel('Frecuencia') ylabel('|H(\omega)|') subplot(212); plot(fre,20*log10(H)) %Modulo en dB
12
# !&" ! !!&" !# !
"
#!
#"
$!
$"
%!
%"
# !&" ! !!&" !# ! %!
"
#!
&'()*'+,-+,./0*0+*1-,0* #" $! $"
%!
!"#
$ 23,/),-/10
$"#
%
!"#
$ 23,/),-/10
$"#
%
%"
'!
456!74
$# $! # ! !
456!746(97
%!! ! !%!! !8!! !
Una característica de los scripts es que las variables que se definen durante su ejecución permanecen en el espacio de entorno una vez ha finalizado la ejecución del programa. A continuación veremos las funciones, que no son más que una porción de código que realiza una tarea determinada, exactamente igual que en lenguaje C. Funciones: Una función en Matlab es un segmento de código como el que se indica a continuación: %Funcion que calcula el cociente y el resto de una division entera %[coc,resto]=div_entera(x,y) function [coc,resto]=div_entera(x,y) %Cualquier funcion debe empezar con %la palabra function %[coc,resto] son las variables que devuelve la funcion % (x,y) son los parametros con los que se llama a la funcion coc=x/y; %este seria el cociente pero no entero coc=fix(coc) %nos quedamos con la parte entera utilizanod la funcion de %Matlab fix resto=x-y*coc; %calculamos el resto de la division %Ahora deberemos guardar el fichero con el nombre dado a la funcion en %este caso div_entera.m
Si el número de variables devueltas es solo 1 no es necesario especificarlos entre paréntesis. Podemos devolver tantas variables como consideremos oportunas. Asimismo tampoco hay restricciones en cuanto al número de variables que se le pueden pasar a una función. Una vez almacenada la función podemos utilizarla de la siguiente forma: 13
>> [a,b]=div_entera(15,4);
Los parámetros de llamada pueden ser variables definidas en el espacio de entorno o valores numérico introducidos directamente. Tanto las variables de llamada como los parámetros devueltos no tienen que coindicir con las variables utilizadas en la declaración de la función ya que únicamente se copian los valores. En este caso las variables cociente y resto se copian en las variables a y b respectivamente. La principal diferencia con un script es que las variables internas utilizadas en la función (coc y resto en este caso), no permanecen en el espacio de entorno al retornar de la función. La forma más sencilla de escribir una función es partir de un script y una vez verificado su funcionamiento transformarlo en un función. Las primeras líneas de comentarios que se incluyen en una función son las que nos muestra la instrucción help. Si ejecutamos: >> help div_entera Funcion que calcula el cociente y el resto de una division entera [coc,resto]=div_entera(x,y)
Otro ejemplo de función: Vamos a escribir una función que nos permita representar gráficamente la función en un intervalo elegido por el usuario %Funcion para representar graficamente la funcion % % x^2-1 %f(x)=--------------% x^2+1 % %En el intervalo [x1,x2] especificado por el usuario con N puntos %La funcion devuelve los valores utilizados para la representacion %[x,y]=representa(x1,x2,N) % function [x,y]=representa(x1,x2,N) x=linspace(x1,x2,N); %Genera un vector N datos linealmente espaciados entre x1 y x2 y=(x.^2-1)./(x.^2+1); %Observese la utilizacion del operador ./ y .^ para que las operaciones se %hagan elemento a elemento ya que dado que x es un vector Matlab intentaria %realizar operaciones vectoriales produciendo un error ya que las %dimensiones no serian correctas plot(x,y,'b') %Embellecemos la grafica grid xlabel('X') ylabel('Amplitud') texto=['Representacion en el intervalo [' num2str(x1) ' , ' num2str(x2) ']']; title(texto) %Matlab puede trabajar con variables de tipo texto, que se especifican %entre comillas simples. %Si en un vector de texto tenemos dos cadenas ej: t=['ho' 'la']; %Matlab las concatena directamente almacenado en la variable t='hola' %La funcion num2str(x1) conviente el numero x1 en una cadena de texto.
Si queremos representar la gráfica en el intervalo [-5,6] con 1000 puntos solo hemos de escribir en el prompt de Matlab >> [x,y]=representa(-5,6,1000);
14
12+32425.67-8592592,9-5.23:6,89;!<9=9!> % $&' $&!
)*+,-./0
$&" $ $ !$ !$&" !$&! !$&' !% !!
!"
!#
$ (
#
"
!
Si no necesitamos las variables devueltas podríamos hacer la llamada >> representa(-5,6,1000);
de esta forma únicamente tendríamos la representación, pero no las variables x e y devueltas por la función. LAS
FUNCIONES PUEDEN SER LLAMADAS DESDE UN SCRIPT U OTRA FUNCIÓN, SOLO HEMOS DE CONOCER LOS PARÁMETROS DE LLAMADA Y LOS PARAMETROS DEVUELTOS.
Las librerías de Matlab están formadas por múltiples funciones que realizan tareas específicas. Nosotros podemos crearnos nuestras propias funciones para reutilizar código en las prácticas. •
Transforma en una función el primer script, de manera que la representación de las tres gráficas se pueda realizar en un intervalo especificado por el usuario. El programa deberá devolver el vector X utilizado y las valores de y, y1 e y2.
Algunas funciones útiles: 1. [M,N] = size(X) %Determina las dimensiones de una varible X 2. A=randn(M,N) %Genera un matriz de números aleatorios de M filas y N columnas, si N o M es la unidad será un vector. 3. N=length(x) %Determina el número de elementos (longitud) de un vector. 4. A=ones(N,M) %Genera una matriz de dimensiones NxM con todos los elementos iguales a 1. 5. A=zeros(N,M) %Genera una matriz de dimensiones NxM con todos los elementos iguales a 0. Ejemplo de utilización >> [fi,co]=size(B) fi = 2 co = 15
4 Todas las funciones de Matlab tienen un conjunto de parámetros de llamada, que se indican entre paréntesis y un conjunto de parámetros devueltos que se indican entre corchetes cuando son más de uno, de esta forma pueden utilizarse y recoger sus resultados. Exactamente igual que las funciones definidas por el usuario. Otras funciones relacionadas con el procesado que utilizaremos en las siguientes prácticas son: zplane, filter, roots, poly, residuez, conv, xcorr, freqz, impz, etc Sentencias de control de flujo: escritura de programas Para controlar la ejecución de un programa necesitamos sentencias condicionales y sentencias de repetición. Dado que en estas sentencias se realizan operaciones lógicas verenos cuáles son estos operadores en Matlab. Operadores relacionales Igual: Distinto Menor que Mayor que Menor o igual que Mayor o igual que
== --> eq(A,B) ~= --> ne(A,B) < --> lt(A,B) > --> gt(A,B) <= --> le(A,B) >= --> ge(A,B)
Operadores Lógicos. AND & --> and OR | --> or NOT ~ --> not OR EXCLUSIVO --> xor Cierto si alguno de los elementos del vector no es cero Cierto si ninguno de los elementos del vector es cero
any(X) all(X)
EL RESULTADO DE LOS OPERADORES LÓGICOS ES CIERTO (1) O FALSO (0). CUANDO HACEMOS OPERACIONES LÓGICAS 0 ES CONSIDERADO COMO FALSO Y CUAQUIER OTRO VALOR COMO VERDADERO. Ejemplo: >> x=1:5 x = 1
2
3
4
5
0
1
1
1
>> b=x>2 b = 0
>> c=gt(x,2) c =
16
0
0
1
1
1
>>
Veamos ahora las sentencias de control de flujo mediante programas ejemplo: Sentencia condicional if A > B disp([num2str(A) ' elseif A < B disp([num2str(A) ' elseif A == B disp([num2str(A) ' else disp('Situacion no end
es mayor que ' num2str(B)]) es menor que ' num2str(B)]) es igual que ' num2str(B)]) esperada')
%Cuando tenemos varias sentencias condicionales pueden sustituirse por la %sentencia switch switch (rem(n,4)) %La funcion rem() calcula el resto de la division entera case 0 disp('El resto es 0') case 1 disp('El resto es 1') case 2 disp('El resto es 2') otherwise disp('El resto no es ni 0 ni 1 ni 2') end
Bucles: %Bucle for %Repite la ejecucion un numero fijo de veces x(1)=0.3; x(2)=0.5*x(1); for n = 3:32 x(n)=0.5*x(n-1)+.3*x(n-2) end %Bucle while %Repite un el bucle mientras la condicion sea cierta. a = 0; fa = -Inf; b = 3; fb = Inf; error=1e-10; while b-a > error*b x = (a+b)/2; fx = x^3-2*x-5; if sign(fx) == sign(fa) a = x; fa = fx; else b = x; fb = fx; end end %La solucion es x = 2.0946
Ejemplo de utilización de sentencias condicionales: Escribe un programa que pida un número por teclado. Mientras este número sea distinto de la unidad hará lo siguiente: si el número es par lo dividirá por 2 y si es impar lo multiplicará por 3 y se sumará 1. El programa debe visualizar por pantalla la secuencia generada. Por ejemplo si el número introducido es el 5 la secuencia mostrada será: 5 16 8 4 2 1. Solución: clear
17
close all N=input('Introduce un numero entero '); n=1; %n es indice para almacenar los elementos de la secuencia en un vector n=1; secuencia(n)=N; %En esta variable almacenamos la secuencia generada while(N~=1) %tambien podemos poner (ne(N,1)) [coc,resto]=div_entera(N,2); if(resto==0) N=N/2; else N=3*N+1; end n=n+1; %Incrementamos el indice de la secuencia secuencia(n)=N; %Almacenamos el nuevo valor end secuencia %Sacamos la secuencia por pantalla
Ejercicio: Transforme el programa anterior en un función a la que se le pase como parámetro N y devuelva la secuencia generada. Solución: %Funcion sin ayuda function [secuencia]=gen_secu(N) n=1; %n es indice para almacenar los elementos de la secuencia en un vector n=1; secuencia(n)=N; %En esta variable almacenamos la secuencia generada while(N~=1) %tambien podemos poner (ne(N,1)) [coc,resto]=div_entera(N,2); if(resto==0) N=N/2; else N=3*N+1; end n=n+1; %Incrementamos el indice de la secuencia secuencia(n)=N; %Almacenamos el nuevo valor end %Guardamos en fichero con el nombre gen_secu.m
Para llamar a la función con N=5, escribiríamos: >> y =gen_secu(5)
OBSERVE QUE EN EN EL INTERIOR DE LA FUNCIÓN NO PUEDE PONERSE LA INSTRUCCIÓN CLEAR YA QUE EN ESE CASO LAS VARIABLES QUE SE LE PASAN A LA FUNCIÓN SERÍAN BORRADAS Y SE PRODUCIRÍA UN ERROR. Como N es el parámetro de llamada a la función también se ha suprimido la instrucción input pues de lo contrario el valor de N que se emplearía no sería aquel con el que se llama a la función sino con el que el usuario introdujese. ENTREGA A TU PROFESOR LOS SIGUIENTE EJERCICIOS ANTES DE LA PRÓXIMA SESIÓN DE LABORATORIO 1.-Una
secuencia
exponencial compleja esta definida como: con fd la frecuencia digital Genera una secuencia de este tipo de 48 muestras con una amplitud de 2, una frecuencia digital de 0.06 y fase inicial 0. Representa la parte real y la parte imaginaria de estas secuencia en gráficas distintas dentro de la misma figura (Utiliza el comando subplot). 18
2.-Una secuencia Chirp está definada como: con n=0:Npuntos-1 a) Escribe una función con el siguiente prototipo: x=GenChirp(A,r,fd,incremento,Npuntos,fase_inicial) b) Usando la función anterior, genera y dibuja la siguiente señal.(r=1)
c) Genere y escuche la secuencia siguente . Puedes escuchar la secuencia generada con la instrucción sound(x,8000) Puedes explicar en qué consiste la función chirp.
19
PRÁCTICA II: MUESTREO En esta práctica se estudiarán las consecuencias del teorema de muestreo, haciendo especial hincapié en los efectos que se pueden producir si se muestrea una señal con la frecuencia demasiado baja. 1) Genera dos períodos correspondientes al muestreo de una sinusoide (coseno) analógica de frecuencia 200 Hz, amplitud 1 y frecuencia de muestreo 1 kHz. 2) Realiza la misma operación que en el apartado 1) pero ahora la sinusoide a muestrear es de 1’2 kHz (el mismo número de puntos). 3) Representa en una misma figura las señales obtenidas en los apartados anteriores, ¿qué ocurre?, ¿qué consecuencias sacas de las gráficas?. 4) Se estudiará el efecto del muestreo sobre el espectro de la señal. Genera un segundo de una sinusoide discreta de frecuencia 100 Hz y una amplitud de 1 con un período de muestreo de 1 ms. Representa el espectro de la señal usando: plot([0:N-1]*Fm/N,abs(fft(y))) %N número de muestras de la señal almacenada en y xlabel(‘Frecuencia (Hz)’) % Fm: frecuencia de muestreo ylabel(‘Módulo’) %y: señal de la que queremos ver su espectro title(‘Espectro’) Comenta los resultados. 5) Realiza lo mismo pero ahora la señal a muestrear es la suma de cuatro sinusoides de amplitud uno y frecuencias 100, 200, 600 y 2100 Hz. Comenta los resultados. 6) Realiza la misma operación sustituyendo la de 2100 Hz por la de 1900 Hz. Comenta los resultados. 7) Repite los apartados 5) y 6) con una señal seno, ¿qué ocurre? 8) Genera una señal cuadrada de 1000 puntos con frecuencia de 150 Hz y la muestreas a 1000 Hz. Representa el espectro como lo has hecho en el apartado 4 y explica el resultado.(utiliza la instrucción square) APLICACIÓN PRÁCTICA I • Genera 88.000 muestras de una sinusoide de frecuencia F=300 Hz muestreada a 44100Hz y escúchala con la instrucción sound. • Repite el apartado anterior sustituyendo la sinusoide por una onda cuadrada. ¿Qué diferencias aprecias? APLICACIÓN PRÁCTICA II • a)
b)
Se dispone de dos ficheros de audio, ce44100.wav muestreados a frecuencias 44.1 kHz y 8kHz respectivamente. Lee el fichero ce44100.wav con Matlab utlizando la función wavread, cuyo formato es: [y,Fs,bits]=wavread(file.wav), y escuchalo con la función sound. Vamos a simular un muestreo de la señal original a 22050Hz y a 11025Hz de la siguiente forma: y22050=y(1:2:end) y y11025=y(1:4:end); escucha estas dos nuevas señales y comenta los resultados. Genera 3 segundos de una sinusoide de 300Hz muestreada a 8 kHz, y escúchala. A continuación reprodúcela pero no utilizando la frecuencia a la que se ha muestreado sino otra, por ejemplo a 20kHz y a 4 kHz. ¿Ha cambiado la frecuencia que se escucha ?. ¿Qué relación hay entre la frecuencia analógica de la señal original y la reconstruida con las nuevas frecuencias ?
20
PRÁCTICA III: CUANTIFICACIÓN. En esta práctica analizaremos los efectos de la cuantificación en la conversión analógica digital así como en la implementación de sistemas discretos, para ello vamos a utilizar la función quanti cuyos argumentos son los siguientes: xq=quanti(x,N,m) siendo x la señal original, N el número de bits (2N niveles), ±m es el rango de entrada y xq la salida cuantificada. La cuantificación se realiza por truncamiento 1. Sin utilizar la función quanti, determina cuáles serán los valores obtenidos al cuantificar por redondeo las 5 primeras muestras de una sinusoide de amplitud 3V, frecuencia 20 Hz y frecuencia de muestreo 100Hz, si se emplea un conversor bipolar de intervalo de entrada ±3 V de 4 bits. Repite el proceso si el conversor cuantifica por truncamiento. 2. Efecto sobre la precisión. Se estudiará el efecto de la cuantificación sobre la precisión de un sistema digital. Se implementará el sistema dado por la siguiente ecuación en diferencias:
este sistema determina la raíz cuadrada de A con la condición inicial de x(0)=1. Determina los resultados obtenidos al calcular la raíz de 2 cuando lo implementamos en un dispositivo cuya salida está entre ±5 y usamos un conversor de 4, 8 y 12 bits. a. Implementa la ecuación en diferencias sin cuantificar y obtén las 15 primeras muestras de la salida. b. Cuantifica la salida del apartado a) y represéntala gráficamente junto a la original. 3. Efecto sobre las formas de onda. Se estudiará el efecto de la cuantificación sobre el aspecto una sinusoide. Para ello genera tres períodos de una sinusoide de amplitud 4, frecuencia 50 y frecuencia de muestreo 1000 Hz y cuantifica ésta suponiendo que el parámetro m=5. Prueba para valores de 4 a 16 bits. Representa la señal de error de cuantificación en cada caso (señal formada por el error en cada una de las muestras). 4. Efecto de la cuantificación sobre una señal de audio. Para ello vamos a utilizar el fichero p44100.wav. Lee dicho fichero con la instrucción [y,Fs,bits]=wavread('p44100.wav'); y escucha su contenido con sound(y,Fs). Repite el proceso cuantificando la señal a 8, 4, 3 y 2 bits. (Utiliza como valor de m el valor máximo en valor absoluto de y). ¿A qué se debe el ruido que se escucha ?, ¿ se podría resolver este problema aumentando la frecuencia de muestreo ?. 5. Repite el apartado 4, cuantificando la señal a 8 bits pero aumentando el valor de intervalo de entrada (m) al conversor (15 ó 20 veces superior al empleado en el apartado anterior). Escucha la señal obtenida y comenta los resultados
21
PRÁCTICA IV: CONVOLUCIÓN Y ECUACIONES EN DIFERENCIAS 1.
Un sistema digital sencillo para eliminar altas frecuencias es el conocido como promediado móvil (Moving Average) definido por la siguiente ecuación en diferencias:
a. b. c. d.
Calcula la respuesta impulsional de este sistema para N=4. Genera una señal de ruido gaussiano de 1000 puntos. Utiliza la instrucción randn. Calcula la salida del sistema (a) para la entrada (b) mediante la convolución. Dibuja la entrada y la salida en una misma gráfica con distinto color. ¿Qué diferencias observas entre ambas señales. e. Determina la salida del sistema para la misma entrada variando N entre 4 y 12 utilizando la convolución (conv). Mostrando en cada caso la señal de entrada y la de salida y en otro subplot el espectro de la señal de entrada y de salida. Qué diferencias observas entre la señal de entrada y la señal de salida y entre sus espectros a medida que aumenta N. 2. Un sistema que elimina las oscilaciones de baja frecuencia en un electrocardiograma (oscilaciones de la línea base) puede ser implementado a partir del siguiente diagrama de bloques:
Donde α es un parámetro variable (0<α<1). a. Escribe las ecuaciones en diferencias que te permiten calcularla salida del sistema a partir de la entrada. b. Calcula la salida a partir de las ecuaciones en diferencias utilizando como entrada la señal contenida en el fichero ecg.mat (Frecuencia de muestreo 1kHz) y visualízalas. Dibuja en otra figura el espectro de la señal de entrada y salida superpuesta con distinto color. c. Repite el apartado anterior utilizando la instrucción filter para calcular la salida tomando para , , d. ¿Qué efecto tiene sobre la señal la modificación del parámetro α? 3. CORRELACIÓN. En este apartado se va a usar la autocorrelación para encontrar posibles periodicidades dentro de una señal. Utilizaremos la serie de datos sobre la aparición de manchas en la superficie solar durante los años 1700 hasta el 1995 (fichero sunspots.dat). Cada una de las muestras se ha tomado con un período de 1 mes. Calcula la autocorelación de la señal e interpreta los máximos de la serie obtenida; ¿existe algún periodo en la señal que se pueda destacar?. El fichero radar.mat contiene la señal original (original) enviadas por el sistema detector y el eco recibido tras localizar un objeto (echo). Usa la correlación cruzada entre las dos señales y, sabiendo que se han usado ondas electromagnéticas (velocidad luz=300000 km/s) y una frecuencia de muestreo de 100 kHz determina la distancia a la que se encuentra el objeto. 22
Sesión 1. Mini Proyecto de Introducción al Procesado Digital de Señales.
23
PRÁCTICA VI: TRANSFORMADA Z. 1.-Un sistema digital tiene como respuesta impulsional la serie infinita h(n) definida como:
a) Obtener la ecuación en diferencias del sistema. b) Calcular la función de transferencia (H(z)) del mismo. c) Comprueba que el sistema obtenido tiene esta respuesta impulsional. Utiliza la función impz(). d) Haciendo uso de la instrucción filter(), determina la salida del sistema cuando la entrada es: i) ii) 2.-Se tiene un sistema que ante la entrada , ¿Qué salida,
, se obtendría si
produce una salida ? Compruébalo con
Matlab. La función residuez, permite realizar una descomposición en fracciones simples. Utilizando esta función determina la expresión general de la salida ante la entrada escalón, e indica cuál es el término transitorio y el estacionario. 3. Los sistemas digitales definidos por las dos ecuaciones en diferencias siguientes se conectan en cascada.:
y(n) = x(n) − x(n −1) + 0.9 ⋅ y ( n −1) Utilizando la instrucción filter, determina la salida intermedia y la salida total del sistema, para la dos posibles conexiones en cascada de estos sistemas. La señal de entrada se encuentra almacenada en el fichero ecg50.mat. Comenta los resultados obtenidos. ¿Qué hace cada uno de estos filtros ?
24
Sesión 2. Mini Proyecto de Introducción al Procesado Digital de Señales
25
PRÁCTICA VII: RESPUESTA EN FRECUENCIA DE LOS SISTEMAS DISCRETOS Primera Parte Descarga el fichero pezdemo-v101.zip de http://www.uv.es/~martsobm, que se encuentra en la practica 8 y descomprímelo en tu carpeta. Ejecuta la demo desde Matlab (versión 6.5) colocándote en la carpeta en la que has descomprimido los ficheros y ejecutando el script pezdemo.m Esta demostración ilustra la relación entre el diagrama de polos y ceros de un sistema con su respuesta impulsional y su respuesta en frecuencia. 1) Coloca un cero (o un par complejos conjugados) en ω=0 y observa qué ocurre cuando mueves (pincha con el ratón sobre él y arrastra) dicho cero variando su módulo, manteniendo el ángulo. 2) Repite el apartado anterior variando el ángulo del cero. ¿ qué le ocurre a la respuesta en frecuencia?, ¿Cuántos términos tiene la respuesta impulsional en los apartados 1 y 2) ? 3) Coloca un polo (o un par complejos conjugados) en w=0 y observa qué ocurre cuando mueves (pincha con el ratón sobre él y arrastra) dicho polo variando su módulo, manteniendo el ángulo, ¿Si el polo tiene modulo mayor que la circunferencia unidad que valores toma la respuesta impulsional ? 4) Repite el apartado anterior pero variando el ángulo.(ω) 5) Cuántos términos tiene la respuesta impulsional en los apartados 3 y 4. 6) Diseña un sistema que tenga 2 ceros sobre la circunferencia unidad con ángulo y dos polos con el mismo ángulo pero con modulo algo menor para asegurar la estabilidad. Observa la respuesta del filtro obtenido. ¿De qué tipo de filtro se trata?. Varia la distancia del polo a la circunferencia unidad y observa la respuesta del sistema. ¿Qué efecto tiene la distancia del polo a la circunferencia unidad?. ¿Qué ocurre si el cero no esta sobre la circunferencia unidad ? 7) Utiliza la opción PZ que coloca un polo y cero que es su recíproco conjugado , ¿qué forma tiene la respuesta en frecuencia en módulo? Segunda Parte 1) Determina la respuesta en frecuencia (módulo) de los siguientes sistemas sin utilizar Matlab, posteriormente verifica los resultados con la función freqz
a) Justifica dicha respuesta de acuerdo con la posición de los polos y los ceros de la función de transferencia, utiliza la función zplane b) ¿Cuál es la relación entre la posición del cero y el orden del promediado ? c) Representa la respuesta en módulo en escala lineal y en dB para los dos sistemas anteriores 2) Dibuja la respuesta en frecuencia en módulo, en escala lineal y en dB, y el diagrama de polos y ceros para el siguiente sistema:
3) Responde a las siguientes preguntas:
26
a) Si un sistema tiene un cero en una determinada frecuencia, que forma tendrá su respuesta en frecuencia en dicha frecuencia. Qué ocurre si el cero se encuentra exactamente sobre la circunferencia unidad. b) Cuál es la forma más sencilla de que un sistema que diseñamos amplifique una determinada frecuencia que se le aplique a su entrada. 4) Un “notch filter” se define como aquel sistema que presenta APROXIMADAMENTE la siguiente respuesta en frecuencia (en módulo)
Implementa el sistema de orden 2 mediante ecuación en diferencias con wo= π/2. Comprueba el funcionamiento considerando como entrada la señal con wa= 0, π/4, π/2, 3π/4. ¿Se podría implementar con un sistema de orden 1? 5) Implementa un sistema de orden 2 mediante ecuación en diferencias que responda a las siguientes características: a) Elimina las frecuencias de 0 y fm/2. b) La ganancia se mantiene cercana a uno en el resto de las frecuencias. c) La ganancia para fm/4 es de 1.1 Determina la respuesta en frecuencia del sistema anterior utilizando la instrucción freqz
27
Sesión 3. Mini Proyecto de Introducción al Procesado Digital de Señales
28