UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
PRÁCTICAS DE DIGITAL SIGNAL PROCESSING - DSP PROFESOR: Esp., Ing., RAMIRO PERDOMO PERDOMO RIVERA
CURSO B-2012
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de Procesamiento Digital de Señales
Práctica 1 REPASO DE SEÑALES Y SISTEMAS Esta práctica pretende introducir los conceptos necesarios para comenzar a trabajar con DSP en MATLAB. Se supone un conocimiento básico del manejo de vectores y matrices en MATLAB, MATLAB, adquirido en semestres anteriores.
1 DSP en MATLAB 1.1 Introducción Objetivo: reforzar los conceptos de las clases de teoría sobre DSP, comenzar a introducir señales señales en MATLAB, hacer operaciones operaciones con ellas ellas y representar representar graficas. Implementar Implementar sistemas sistemas muy sencillos. Más en particular, en esta práctica se introducirán los conceptos de vector de valores y vector base de tiempos y se utilizaran para representar algunas señales básicas. b ásicas. Después, Desp ués, se s e exploraran las operaciones elementales con señales (suma, producto, desplazamientos...). Por último, se introducirá el concepto de sistema comparándolo con el concepto de "función de MATLAB". La práctica propone ejercicios de dificultad creciente. Los primeros ejercicios de cada tipo se explican muy detalladamente en el enunciado, los siguientes deben ser resueltos por los alumnos utilizando los conocimientos de teoría y los que se han ido obteniendo en ejercicios anteriores. Duración: Una o dos sesiones.
1.2 Señales en MATLAB 1.2.1 Limitaciones de MATLAB Antess de nada Ante nada,, nóte nótese se que que en MATL MATLAB AB (com (comoo en cual cualqu quie ierr otro otro simu simula lado dor r informático) informático),, es imposible introducir introducir funciones de variable variable continua ya que se trabaja trabaja con conjuntos discretos de valores (vectores y matrices). Por tanto, para poder trabajar con señales señales continu continuas as vamos vamos a necesi necesitar tar tomar tomar un con conjun junto to "represe "representa ntati tivo" vo" de valores. Eso significa que nuestra señales van a ser conjuntos de valores tornados de la verdadera señales continua en instantes determinados, generalmente se tratara de inst instant antes es equ equie iespa spaci ciad ados os (mue (muest stre reoo unif unifor orme me). ). No vamo vamoss a estud estudia iarr aquí aquí las las propiedades y consecuencias consecue ncias del proceso de muestreo. Solamente comentaremos que el espacio temporal entre muestras deberá ser lo "suficientemente pequeño".
Figura 1.- Proceso de muestreo
1.2.2 Representación de una Señales Básica En este apartado vamos a representar una señal es muy simple. Vamos a empezar por una sinusoide, sinusoide, esto es: x(t) = sen(w0t). sen(w0t). Vamos a tomar el valor de frecuencia frecuencia angular de 2π w0 = Tc/4 rad/seg. Esto corresponde con un periodo de T0 =
frecuencia hertziana de
f o
=
1 T o
=
ω o
2π
=
ω o
= 8 seg y con una
0.125 Hz
Intentaremos iniciar con este apartado una buena costumbre. Vamos a almacenar los valores de x(t) en un vector x al que Llamaremos "vector de valores". En paralelo crearemos otro vector con los valores de la variable t. Este vector (t) lo llamaremos "vector de base de tiempos" o, simplemente, "base de tiempos". Nótese que el contenido de t(l) es un instante de tiempo (en segundos) y x(l) es el valor (voltios) que toma la señal en ese instante (lo mismo para t(2) y x(2)... ). Nótese también que la señales en MATLAB es el conjunto de los dos vectores (t y x).
Vamos a representar la señal de interés entre -10 y 10 segundos. Vamos a tomar como incremento 0.1. En estas condiciones la base de tiempos la podemos crear así: t = -10:0.1:10; Y el vector de valores: wo = pi/4; x = sin(wo*t); % Nótese el uso de la notación: % "seno de un vector''. Para hacer una representación grafica usaremos la función plot. plot(t,x); Y tendríamos que obtener una grafica similar a esta:
Figura 2.- Grafica de la función senoidal. Pregunta 1.- Utilizando la herramienta “lupa" de las graficas de MATLAB, intente medir el periodo de esta señales. Respuesta:
1.2.6 Representación de Otras Señales En este apartado se les pide que representen otras señales sencillas y no se les da el problema resuelto como en el apartado anterior. Vamos a representar todas entre —5 y 10 segundos con incrementos de 0.05. Fíjense que al tratarse de representaciones con las mismas características temporales, pueden usar la misma base de tiempos para todos los casos. A continuación se indican las señales a representar y los resultados que se deberían obtener. x(t) = u(t) (señales "escalón unitario'5).
Figura 6.- Escalón unitario.
Pregunta 2.- Las gráficas de sinusoides amortiguadas incluyen representación de la envolvente de la curva (en línea roja discontinua). ¿Como haríamos para incluir esta representación?
1.6 Operaciones Simples con Señales Ahora vamos a hacer algunas Operaciones simples con una señal definida por su grafica. Vamos a comenzar introduciendo en MATLAB la señal x(t) dada por este dibujo:
Figura 16.- Señal definida por una grafica. En principio, nos vale cualquier base de tiempos que abarque desde el —2 hasta el 2. Vamos a darle un poco de margen por ambos lados y la hacemos de —6 a 6 con incrementos de 0.05. t = -6:0.05:6; Hay que calcular la ecuación para la recta que va entre t = -2 y t = 0. Esta recta tendrá la ecuación típica: r(t) = mt+b donde: • m es la pendiente. Al tratarse de una recta creciente debe resultar positiva. El valor de m se puede calcular como la tangente del ángulo que forma la recta con 1 0 1 y1 y 0 m x1 x0 0 ( 2 ) 2 el eje horizontal (eje t). En este caso: • b es el término independiente. Se puede calcular conociendo m y cualquier punto de la recta, por ejemplo r(-2) = 06 o r(0) = 1 (con este último obtenemos una ecuación extremadamente simple que resulta b = 1). De hecho, el término independiente de una recta siempre es igual al punto en que corta al eje vertical (cuidado con esta propiedad porque, a veces, no es posible ver el punto de forma tan trivial como en este caso). =
− −
−
=
−
=
−
Por tanto, x(t) entre —2 y 0 (y solo entre —2 y 0) es igual a la expresión: t/2+1. Por supuesto, x(t) es igual a 1 entre 0 y 2 e igual a cero en el resto de intervalos. Veamos como crear el vector de señal. L « length (t); % Averiguar la longitud x = zeros(1,L); % Primero todo ceros pi = f ind(t=-2) ; % Busco el -2 p2 = find(t=0) ; % Busco el 0 p6 = find(t==2) ; % Busco el 2 % Con el incremento elegido; -2, 0 y 2 estarán % en la base de tiempos x(pl:p2) = t(pl:p2)/2+l; % Parte de recta creciente x(p2:p6) = 1; % Parte constante % x en t=0 lo hemos calculado con dos fórmulas % dando el mismo resultado
Figura 17.- Resultado de introducir x(t) en MATLAB. Una vez que tenemos la señal x(t) en MATLAB vamos a realizar algunas operaciones simples con:
a) Traslación: calcular y representar xi(t) = x(t-6). Para este apartado veremos dos soluciones:
EJECUTE LAS DOS Y FÍJENSE BIEN EN LA DIFERENCIA DE MÉTODO. Solución 1: Sabemos que la señal es la misma pero todos los puntos se trasladaran de tiempo t a tiempo t+6. Esto es: basta con sumar 6 a la base de tiempos. tl = t + 6; xl = x;
Solución 2: Movemos el vector de valores 6 unidades de tiempo hacia delante (y hacemos crecer la base de tiempos). t__aux = 0.05:0.05:6; % Tres unidades de tiempo % No empieza en cero porque lo vamos a % añadir por la derecha L6 = length(t_aux); % Averiguar cuantos valores son los 6 segundos x__aux = zeros (1,L6) ; % Ceros para poner por la izda xl = [x__aux x] ; % Los valores son los mismos % (con ceros por delante) tl = [t max(t)+t_aux]; % Crear nueva base de tiempos
Figura 18.- Señal desplazada 6 segundos.
b) Escalado : calcular y representar x2(t) = x(2t). Sabemos que la serial es la misma pero "comprimida" (a la mitad). Por tanto, el punto situado en t pasara a t/2. Esto es: basta con dividir por 2 la base de tiempos.
t2 = t/2; % Tres unidades de tiempo x2 = x;
Figura 19.- Señal comprimida a la mitad. Nótese que haciendo esto hemos dividido por dos el "incremento temporal'' de la base de tiempos (los tiempos ahora están separados 0.025 segundos). c) Reflexión: calcular y representar x6(t) = x(-t). Ahora se trata de que el punto situado en t pasará a -t. La primera idea es cambiar de signo la base de tiempos. Eso es correcto pero no es suficiente porque tanto la base de tiempos como el vector de valores están en orden inverso al que debieran. Eso debemos resolverlo invirtiendo ambos vectores. t6 = -t; % Base de tiempos invertida x6 = x; % Los valores son los mismos final = length(t6); t6 = t6(final:-!:!); x6 = x6(final:-!:!); % Invertir
Figura 20.- Señal reflejada. d) Sumar: sumar x(t) y Xj(t) (la señal del apartado a). Para sumar dos señales no hay mas que sumar los valores pero DEBEMOS TENER IA MISMA BASE DE TIEMPO PARA AMBAS SEÑALES Si conservamos las variables xl y tl tenemos mucho hecho (se refiere a las variables del apartado a, solución 2). Viendo que tl es la misma base de tiempos de x pero extendida 6 unidades de tiempo no hace falta mas que extender igualmente la base de tiempos de x (y añadir el número adecuado de ceros a los valores de x). Lextra = length(t1) - length(t); %Diferencia de longitudes temporales (en número de valores) t = t1 % Extendemos la base de tiempo de x x = [x zeros(1,Lextra)]; % Añadimos ceros a x t4 = t; x4 = x+x1; % Calcular la señal suma (x4)
Figura 21 Señal suma Pregunta 6 Si les da tiempo, generar la señal x(t) de la figura y calculen y representen las señales: x1 = x(t-2), x2 = x(t/6), x6 = x(-t), x4 = x(t) - x1(t)
1.4 Concepto de Sistema Ahora vamos a realizar en MATLAB un par de sistemas sencillos. Hemos definido un sistema como cualquier ente capaz de transformar señales, es decir, que realiza operaciones sobre las señales. Para afianzar la idea de que un sistema transforma señales vamos a emplear el concepto de “función de MATLAB”, para crear nuestros primeros sistemas. Un sistema lo vamos a definir como una función que podemos crear con el editor de MATLAB ( o con cualquier otro editor) utilizando un encabezamiento del tipo: Function [y, ty] = nombre(x, tx) % % Instrucciones que generen “y” y ty” % ( a partir de x y de tx) % %Estas instrucciones se graban en “Nombre.m” Obsérvese que hemos querido reforzar la idea de que una señal en MATLAB es una base de tiempos y un vector de valores. La función recibe como argumentos los dos vectores de la señal x(t) (entrada) y devuelve como resultados los dos vectores de la señal y(t) (salida). Por ejemplo, para implementar el sistema que eleva al cuadrado (y(t)= [x(t)] 2 podemos usar la función: function [y,ty] = Cuadrado (x , tx) ty = tx; % La base de tiempos es la misma y = x.A2; % Los valores se elevan al cuadrado
Si conservamos las variables t y x del apartado anterior podemos probar esta función: [y, ty] = Cuadrado (x , t) ; % Fíjense en la diferencia entre el nombre % de la variable t y el del argumento tx Obtenemos la siguiente señal (recordemos que habíamos extendido la base de tiempos hasta t=6):
Figura 22.- Señal del apartado 4 elevada al cuadrado.
Ejercicio: crear una función que desplace la señal de entrada to unidades de tiempo hacia la derecha (es decir: implementen el sistema: y(t) = x(t-to)). Nótese que la función deberá recibir al menos to como parámetro de entrada adicional, con lo que el encabezado será: function [y,ty] = Desplazamiento (x,tx,tO) Nota.- Probablemente, también será necesario conocer el incremento temporal que existe entre los valores de la base de tiempos de entrada. Podemos añadirlo como parámetro pero ¿,Podemos deducirlo a partir de tx?
Bibliografía [1] C. Burrus, J. McClelland, A. Oppenheim, T. Parks, R. Schafer, and H. Scuessler, editors. Computer-based exemses for Signal Processing using MATLAB. Prentice-Hall, 1994 [2] John Eaton. Octave Documentation. GNU, 1997 (En fotocopiadora). [6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing. Prentice-Hall, 1989. [4] A.V. Oppenheim, A.S. Willsky, and I.T. Young. Signals and Systems. Prentice-Hall, 1986.
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales Práctica 2 OBJETIVO
Esta práctica pretende reforzar los conocimientos adquiridos en clase de teoría sobre la convolución. La convolución se puede ver como una operación entre señales que además es el fundamento de todos los sistemas lineales e invariantes (ya que todo sistema LTI se puede formular como una convolución). Al realizar esta práctica se encontraran con una serie de ejercicios de cálculo de convolución que tendrán que realizar con ayuda del computador. Lo más conveniente (como se les irán recomendando a lo largo del texto), es comprobar los resultados de esos ejercicios con los resultados que pueden obtener analíticamente.
2 La Integral de Convolución en Tiempo Continuo 2.1 Introducción Objetivo: reforzar los conceptos teóricos sobre convolución, comprender las limitaciones del computador para trabajar con las señales y su convolución, comparar los resultados del computador con los analíticos. La práctica propone ejercicios sobre convolución de dificultad creciente. Los primeros ejercicios se explican muy detalladamente en el enunciado, los siguientes deben ser resueltos por los alumnos utilizando los conocimientos de teoría y los que se van ido obteniendo de ejercicios anteriores. Duración: Una o dos sesiones.
3.2 Convoluciones de Señales Sencillas 3.2.1 Convolución de Señales Finitas Vamos a calcular la convolución entre dos señales x(t) y h(t) que, en este caso, van a ser iguales. x y h van a ser (ambas) un pulso rectangular centrado en t=0, de amplitud 1 y que
se extenderá desde t=-l a t=l. Si dibujamos esa señal antes de intentar introducir en MATLAB, k grafica seria:
Vamos a ver ahora (aunque esta parte pertenece a la práctica anterior) como generada en MATLAB: t = -5:0.1:5; % Base de tiempos (común) L = length (t); % Longitud pi = find(t==-l); p2 = find(t==l); % Localizar los puntos x = zeros(1,L); x(pl:p2) = 1; % Pulso entre -1 y 1 h = x; % h es igual a x Después de esto podemos representar x (o h) con plot(tpc) obteniendo:
Y ahora implemente elaborando un programa en MATLAB que realice la función "Convolución" con la instrucción Function para calcular y(t)=x(t)*h(t). Nótese que el resultado es una señal que nos viene dada con un vector de valores y una nueva base de tiempos. Si representamos el resultado con un plot, obtendremos:
Pregunta 1.- Comprueben analíticamente el resultado de esta convolución calculándola como se vió en su asignatura Señales y Sistemas.
Como segundo ejercicio, deben realizar la convolución de las siguientes señales: --x(t): pulso triangular centrado en t=0, altura 1, entre t=-2 y t=2. - h(t): pulso rectangular de altura 1, entre t=0 y t=2. Las señales descritas son (gráficamente):
Nótese que para generar el triángulo es necesario calcular las ecuaciones de dos rectas (recuerden como se hizo en la practica anterior). Los resultados en MATLAB deberán ser
Pregunta 2.- Comprueben analíticamente el resultado de esta convolución calculándola como hemos visto en clase de teoría.
2.6.2 Convolución con Señales Infinitas La convolución donde intervienen señales infinitas en tiempo son imposibles de calcular en MATLAB ya que no se puede tener un vector de valores de tamaño infinito. Sin embargo, vamos a hacer algunos ejemplos donde veremos que si es posible hacer cálculos aproximados. Empezaremos por una señal que podríamos llamar "cuasi-fmita" como es el caso de 3 32 t h(t ) e u (t ) 2 "cuasi-finita" porque a partir de t~2 sus valores son prácticamente nulos). Lo llamamos h(t) porque la intención es interpretarlo como un sistema que va a recibir varias entradas, esto es: vamos a hacer convolución de esta señal con varias señales x(t) diferentes. −
=
Como curiosidad interesante, esta señal es la respuesta al impulso de un sistema definido por un circuito lineal, se trata del siguiente circuito RC (donde el producto RC, llamado constante de tiempo, debe valer 2/6 segundos):
Vamos a empezar haciendo la convolución de h(t) con el pulso del apartado anterior. Esto es:
Para resolver este problema deben decidir hasta que valor de t van a tomar valores de la exponencial (debe ser un valor donde h(t) ya sea muy pequeño). Los resultados en MATLAB deben ser
Pregunta 6.- Comprueben analíticamente el resultado de esta convolución calculándola como hemos visto en clase de teoría.
Pregunta 4.- Este ejercicio también se puede resolver analíticamente utilizando la teoría de circuitos que estudiaron en la primera parte de la carrera. Se trata de un circuito RC en régimen transitorio. La entrada, x(t), pasa de 0 a 1 en el instante t=0 provocando que i(t) deje de ser nula y que el condensador se empiece a cargar. Si no cambiara nada, la situación en tiempo infinito seria de tensión de salida 1 y corriente 0 (condensador a tope de carga). Sin embargo en t=2, con el condensador a media carga, x(t) baja a 0 generando un segundo transitorio. En ese segundo transitorio el condensador se descargara llegando a tensión nula en el infinito (descarga total). Comprueben que si resuelve el circuito usando las leyes de Kirchoff (y las ecuaciones diferenciales que generan), el resultado para y(t) es el mismo. Ahora vamos a cambiar la señal de entrada por x(t)=u(t). Nótese que vamos a aumentar considerablemente la complicación porque se trata de dos señales infinitas y, además, la x(t) no tiende a cero como si hace la h(t). No queda más remedio que truncar la señal u(t) y, aunque escojamos un valor grande, realmente vamos a convertir x(t) en un pulso finito (como el del ejercicio anterior). A continuación, presentamos los resultados de truncar u(t) en t=5 (final de nuestra base de tiempos). Véase lo parecidos que son los resultados con los del ejercicio anterior.
Para interpretar estas gráficas (intentando saber que es correcto y que no) podemos recordar el circuito de la figura 6. Si la entrada es u(t) tendremos un único transitorio que acaba con la carga completa del condensador y el potencial de salida a 1. Es decir: y(t) alcanza asintoticamente el valor 1 y nunca decrece. El decrecimiento a partir de t=5 que se ve en la figura 9 se debe a la aproximación de cortar u(t) a partir de t=5. Si suponemos que u(t) se trunca hasta tiempo to, pueden dibujar las graficas de h(t) y x(t- t) y comprobar que no hay diferencia en el resultado de convolución hasta t=to. (a partir de ese punto comienzan los errores que podemos llamar "efectos de borde" ya que se deben a que toda señal simulada debe tener un final o “borde”). Pregunta 5.- Dibujar las señales h(T) y x(t- t) como se propone en el párrafo anterior (consideren x(t) truncada en un instante genérico t=to).
Pregunta 6.- Calcular analíticamente la convolución de x(t) y h(t) (sin truncar ninguna señal) y compruebe que la grafica de MATLAB solo es correcta hasta t=5. Pregunta 7.- Calcular y(t) utilizando teoría de circuitos. E1 resultado debe ser idéntico al de la pregunta 6. 2.6.6 Convolución con Señales Periódicas Vamos a terminar la práctica con un ejemplo de señales periódicas. Mantendremos h(t) igual a la del apartado anterior (señal exponencial que además es respuesta al impulso de un circuito RC) y haremos x(t) = cos(2ftt). Cuando una de las dos señales de una convolución es periódica, el resultado es periódico con el mismo periodo (no vamos a ver el caso en que las dos señales no 2π T o 1 π 2 sean periódicas). En este caso, x(t) es periódica con periodo fundamental seg con lo que la salida y(t) también deberá serlo. Nótese que si interpretamos la convolución como el calculo del potencial de salida de un circuito, este caso es un ejemplo de análisis en régimen senoidal permanente (que podría resolverse utilizando fasores). El problema de resolver esta convolución utilizando MATLAB es muy parecido al del ejercicio anterior. Esto es: allí donde trunquemos el coseno van a aparecer problemas ya que es una señal infinita y nos estamos quedando con un trozo. Tomando el coseno entre t=-5 y t=5 se obtiene (en MATLAB): =
=
En la gráfica anterior solo podemos fiarnos de los resultados entre t= -4 y t=4 (aproximadamente). En este rango vemos que la señal de salida es periódica de periodo 1. Es más: la salida es también un coseno, esa es una propiedad de los sistemas LTI que todavía no hemos visto: la respuesta a un coseno siempre es otro coseno (donde puede cambiar la amplitud y la fase). Pregunta 9.- Utilicen la técnica de los fasores para calcular la salida del circuito de la figura 6. Entre t= -4 y t=4 debería ser igual a la gráfica de MATLAB. Pregunta 10.- ¿De donde hemos sacado el intervalo de valores fiables: t= -4 y t=4? Consejo: dibujar x(t-i) y h(i) y ver en que rango de t el dibujo es equivalente al que se tendría con un coseno infinito. Nótese que es aproximado y que suponer que va de -4 a 4 es equivalente a suponer que la exponencial ya esta muy debilitada en t=l.
Página
Bibliografía [1] C. Burrus, J. McClelland, A. Oppenheim, T. Parks, R. Schafer, and H. Scuessler, editors. Computer-based exemses for Signal Processing using MATLAB. Prentice-Hall, 1994 [2] John Eaton. Octave Documentation. GNU, 1997 (En fotocopiadora). [6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing. Prentice-Hall, 1989. [4] A.V. Oppenheim, A.S. Willsky, and I.T. Young. Signals and Systems. Prentice-Hall, 1986.
Página
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales Práctica 3
MUESTREO INTERPOLACIÓN Y ALIASING OBJETIVO: Observar el muestreo de una señal y el efecto conocido como “Alising”.
MATERIAL Y EQUIPO: Computadora con MATLAB
FUNDAMENTOS TEORÍCOS Si la frecuencia más alta contenida en una señal analógica xa(t) es Fmax y la señal se muestrea a una tasa Fs >2 Fmax, entonces xa(t) se puede recuperar totalmente a partir de sus muestras mediante la siguiente función:
Si el criterio no es satisfecho, existirán frecuencias cuyo muestreo coincide con otras
(el llamado aliasing).
Página
DESARROLLO: 1. Simular el muestreo de la señal y(t) a una frecuencia de muestreo de 10 Hz y graficarla.
2. Recuperar la señal analógica y(t) utilizando las muestras obtenidas en el punto 1 y aplicando la sumatoria de funciones de interpolación.
3. Simular el muestreo de una onda sinusoidal pura de 300 Hz. Utilizar una frecuencia muestreo de 800 Hz.
de
4. Simular el muestreo de las ondas sinusoidales puras cuyas frecuencias se indican en los incisos a, b, c, d y e a una frecuencia de muestreo de 800 Hz. a. 125 Hz b. 215 Hz c. 305 Hz d. 395 Hz e. 500 Hz Observar los cambios en las formas de onda de la señal muestreada y la señal recuperada. Graficar en la misma ventana la señal original, la señal muestreada y la señal recuperada. 5. Repetir el punto anterior pero ahora con un periodo de muestreo de 1 ms y con las frecuencias de a. 7525 Hz b. 7650 Hz c. 7775 Hz d. 7900 Hz Observar los cambios en las formas de onda de la señal muestreada y la señal recuperada. Graficar en la misma ventana la señal original, la señal muestreada y la señal recuperada.
Página
6. Encontrar tres señales diferentes que tengan la misma representación discreta. 7. Establecer conclusiones
PROGRAMA DE APOYO %Practica 2. Procesamiento Digital de Senales %Muestreo e interpolacion %Limpiar variables, funciones, ventana de comandos y figuras clear; clc; clf %Constantes f1=4.7 %Armónico de mayor frecuencia de la señal f2=2 %frecuencia de la señal M=3
%Limite superior de la amplitud para la grafica m=-3
%Armónico
de
menor
%Limite inferior de la
%amplitud para la grafica a=-1 %Tiempo de inicio de la señal b=1 %Tiempo final de la señal fs=10 %frecuencia de muestreo ab=b-a %Intervalo de la señal %Crear el vector de tiempo cuasicontinuo fc=100*f1 %Frecuencia usada para visualizar la señal como analógica tc=1/fc; t=[a:tc:b]; %Crea el vector de tiempo de señal muestreada ts=1/fs; td=[a:ts:b]; %Señal a interpolar yc=sin(2*pi*f1*t)+sin(2*pi*f2*t); %Señal muestreada yd=sin(2*pi*f1*td)+sin(2*pi*f2*td); %Graficar ambas señales subplot(2,1,1) stem(td,yd,'o') title('Señal muestreada') axis([a,b,m,M]) subplot(2,1,2) %Dos cuadros de figura en vertical plot(t,yc) title('Señal original') axis([a,b,m,M]) hold on %Rutina para recuperar la señal original a partir de la muestreada [xm,xn]=meshgrid(t,td);
Página
aux=(pi/ts)*(xm-xn)+eps; yi=yd*(sin(aux)./(aux)); %Grafica de la señal recuperada plot(t,yi,'r') title('Señal recuperada') axis([a,b,m,M]) hold off
CONCLUSIONES Bibliografía [1] C. Burrus, J. McClelland, A. Oppenheim, T. Parks, R. Schafer, and H. Scuessler, editors. Computer-based exemses for Signal Processing using MATLAB. Prentice-Hall, 1994 [2] John Eaton. Octave Documentation. GNU, 1997 (En fotocopiadora). [6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing. Prentice-Hall, 1989. [4] A.V. Oppenheim, A.S. Willsky, and I.T. Young. Signals and Systems. Prentice-Hall, 1986.
Página
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales Práctica 4 CUANTIZACIÓN OBJETIVO: Analizar la relación entre el ruido de cuantización, la frecuencia de muestreo y el paso de cuantización.
MATERIAL Y EQUIPO: Computadora con MATLAB
INTRODUCCIÓN: La cuantización se refiere al proceso en el que una señal analógica se aproxima a una señal que puede tomar solamente un número finito de valores. La señal digital resultado de la cuantificación es diferente a la señal analógica que la originó debido a lo que se conoce como error de cuantificación. El error de cuantificación se interpreta como un ruido añadido a la señal tras el proceso de decodificación digital. La cuantificación no tendrá ninguna consecuencia si el ruido añadido con la cuantificación se mantiene por debajo del ruido presente en la señal analógica original.
DESARROLLO: Hacer un programa en MATLAB que simule la conversión analógica a digital y que recupere la señal original a partir de la señal digital con la sumatoria de
Página
funciones de interpolación. (Analizar el programa que aparece al final de este documento). Trabajar con la señal definida por:
3. Muestrear la señal a una frecuencia 3 veces mayor a la del armónico de mayor frecuencia y luego digitalizar con una resolución de 1 unidad. Graficar el muestreo de la señal y su digitalización. 4. Reconstruir la señal analógica a partir de sus muestras digitalizadas y compararla con la señal original. Usar una frecuencia de muestreo 3 veces mayor a la del armónico mas alto. 5. Hacer una grafica de la magnitud del error de cuantización. Hacer una resta entre la seña muestreada la señal digitalizada. 6. Aumentar la frecuencia de muestreo 20 veces más y observar si se reduce la magnitud del error. Se debe incrementar la frecuencia de muestreo a 60 veces mayor que el armónico mayor. 7. Regresar ahora a la frecuencia de muestreo anterior (3 veces mayor a la del armónico de mayor frecuencia) y mejorar la resolución en un factor de 20, es decir usar una resolución de 0.05 unidades. Observar lo que sucede con la magnitud del error. 8. Simular un la digitalización de un convertidor A/D de 3 bits con niveles entre 0 y 5 Volt. 9. Anotar conclusiones
CONCLUSIONES Bibliografía [1] C. Burrus, J. McClelland, A. Oppenheim, T. Parks, R. Schafer, and H. Scuessler, editors. Computer-based exemses for Signal Processing using MATLAB. Prentice-Hall, 1994 [2] John Eaton. Octave Documentation. GNU, 1997 (En fotocopiadora). [6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing. Prentice-Hall, 1989. [4] A.V. Oppenheim, A.S. Willsky, and I.T. Young. Signals and Systems. Prentice-Hall, 1986.
Página
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales Práctica 5 OPERACIÓNES CON SEÑALES OBJETIVO: Realizar operaciones con señales digitales de audio.
MATERIAL Y EQUIPO: Computadora con MATLAB Audifonos
INTRODUCCIÓN: En muchas ocasiones es necesario considerar señales que son el resultado de una pequeña transformación de otra señal. Un tipo importante de este tipo de transformación en el tiempo es el corrimiento en el tiempo donde la señal original es desplazada en el eje del tiempo, ya sea para atrasarla o adelantarla; para una señal discreta equivale a un corrimiento en la variable independiente n y se representa como x[n-no] (“no” es el corrimiento). Otro es la inversión en el tiempo donde la señal es vista como una reflexión en el en n=0 y se representa como x[-n]. Además de estas operaciones están las suma, resta y multiplicación entre señales.
DESARROLLO: 1. Guardar los archivos de texto con el sonido digitalizado en la carpeta
“work” de MATLAB.
2. Abrir uno de los archivos para verificar que se trata de datos numéricos.
Enseguida se muestra el contenido del archivo
3. Cargar los datos del sonido digitalizado (archivo de texto) al ‘workspace’
de MATLAB. Usar la función load. Guardar en una variable con el mismo nombre . Se muestra el ejemplo con el archivo datop.
Después de hacer esto debe aparecer la variable en la “workspace”
ventana
4. Generar el vector ‘np’, ‘nd’, ‘nde’, ‘ns’ correspondiente al tiempo de las
variables. Deben tener el mismo número de elementos que los que tiene la variable. La variable datop tiene 14001 elementos. Esto se puede ver en la ventana “workspace”
Entonces el vector de tiempo debe tener 14001 elementos y se genera con la siguiente instrucción Después de esto debe aparecer una nueva variable llamada “np” en el “workspace”
5. Graficar cada una de las señales usando la función stem (usar puntos en la
grafica). Enseguida se muestra el ejemplo para la gráfica de datop vs np.
Con esto se obtiene la siguiente gráfica
6. Reproducir las señales como sonido a una frecuencia de muestreo de
8000 Hz utilizando la función de MATLAB ‘sound(señal, Fs)’. Fs es la frecuencia de muestreo, en este caso Fs=8000. 7. Reproducir las mismas señales con una frecuencia de muestreo de 5000 y
observar la diferencia en el sonido. Anotar lo que sucede y dar una explicación. 8. Reproducir las mismas señales con una frecuencia de muestreo de 13000 y observar la diferencia en el sonido. Anotar lo que sucede y dar una explicación. 9. Obtener la suma de las cuatro señales y guardarla en una variable con nombre datosuma. 10. Graficar la señal datosuma y reproducirla como sonido. 11. Realizar las operaciones necesarias para concatenar las cuatro señales en el mismo orden que en (1) dándoles un espaciamiento de 2000 muestras y guardarlas en la señal ‘datosconcatenados’. Graficar y reproducir. 12. Invertir la señal ‘datosconcatenados’ y guardalos en ‘invertida’. 13. Graficar y reproducir la señal ‘invertida’
CONCLUSIONES
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales PRÁCTICA 6 RESPUESTA DEL SISTEMA OBJETIVO Obtener la respuesta de un sistema discreto por medio de la convolución y la evaluación directa de las ecuaciones de diferencias.
MATERIAL Y EQUIPO: Computadora con MATLAB
INTRODUCCIÓN: La convolución discreta es una operación matemática entre dos señales discretas (el símbolo de la operación es un * asterisco) que tiene gran importancia en procesamiento digital de señales ya que la respuesta y[n] de un sistema lineal e invariante en el tiempo se puede obtener a partir de su respuesta al impulso h[n] y la señal de entrada x[n].
La convolución discreta está definida por la sumatoria
Un sistema discreto lineal e invariante en el tiempo puede ser descrito matemáticamente por una ecuación de diferencias de coeficientes constantes de la forma,
Esta ecuación describe la relación entre la señal de salida y la señal de entrada. Se puede ver que están involucrados los valores de instantes actuales así como de instantes anteriores tanto para la señal de entrada como para la de salida. De esta expresión se puede despejar el valor actual de la señal de salida y[n] con lo que se obtiene,
Donde se puede ver que la ecuación de diferencias da una forma recursiva para obtener la salida actual del sistema utilizando los valores de la señal de entrada previos así como el actual y también los valores previos de la misma señal de salida. MATLAB cuenta con una función que evalúa este tipo de ecuación de diferencias dada una cierta señal de entrada. La función se denomina FILTER. FILTER filtro digital. Y = FILTER(B,A,X) filtra los datos del vector X con el filtro descrito por los vectores A y B. La función FILTER evalúa la siguiente ecuación de diferencias. a(1)*y(n) + a(2)*y(n-1) + a(3)*y(n-2) + ... = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-2) + ... Donde B=[b(1) b(2) …], A=[a(1) a(2) …], X= señal de entrada
DESARROLLO: 1. Hacer una función que calcule la convolución de dos funciones. 2. Obtener la convolución de x[n] y y[n] definidas a continuación.
3. La respuesta al impulso de un sistema discreto lineal e invariante en el tiempo es h[n]. ¿Qué respuesta tendrá este sistema si se le aplica la señal x[n] definida a continuación? (Utilizar la convolución).
Solución
n
4. Hacer una función que genere una secuencia exponencial del tipo a . n (ejemplo 0.7 ). 5. Obtener la respuesta de un sistema discreto lineal e invariante cuando la entrada es x[n]. La respuesta al impulso del sistema es h[n]. Graficar la señal de entrada, la señal de salida y la respuesta al impulso del sistema.
8. Un sistema esta descrito por la ecuación de diferencias donde x[n] es la entrada y y[n] es la salida. Calcular y graficar la respuesta al impulso y al escalón unitario. 9. Se tiene un diferenciador digital descrito por: Obtener la respuesta de este sistema para las siguientes entradas Pulso rectangular Pulso triangular pulso sinusoidal
CONCLUSIONES: Bibliografía [1] C. Burrus, J. McClelland, A. Oppenheim, T. Parks, R. Schafer, and H. Scuessler, editors. Computer-based exemses for Signal Processing using MATLAB. Prentice-Hall, 1994 [2] John Eaton. Octave Documentation. GNU, 1997 (En fotocopiadora). [6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing. Prentice-Hall, 1989. [4] A.V. Oppenheim, A.S. Willsky, and I.T. Young. Signals and Systems. Prentice-Hall, 1986.
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales Práctica 7 DESARROLLO EN FRACCIONES PARCIALES Y GRAFICAS DE POLOS Y CEROS OBJETIVO: Obtener el desarrollo de fracciones parciales y las graficas de polos y ceros con ayuda de MATLAB.
MATERIAL Y EQUIPO: Computadora con MATLAB
INTRODUCCIÓN: Una operación que debe hacerse con frecuencia en el cálculo de la respuesta del sistema es el desarrollo en fracciones parciales, para que permita el uso de tablas de transformadas. Para esta operación se cuenta con la función de MATLAB “residuez”. Las gráficas de polos ceros de la función de transferencia son muy importantes ya que nos dan información sobre la dinámica del sistema. MATLAB cuenta con la función “ zplane” que permite realizar este tipo de gráficas de forma simple introduciendo los coeficientes de la función de transferencia
DESARROLLO: Hacer el desarrollo en fracciones parciales y las gráficas de polos y ceros de las siguientes funciones.
2
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales Práctica 8 RESPUESTA EN FRECUENCIA OBJETIVO: Aplicar las funciones de MATLAB para obtener la respuesta en frecuencia de sistemas discretos.
MATERIAL Y EQUIPO: Computadora con MATLAB
INTRODUCCIÓN: Una de las características de un sistema lineal e invariante en el tiempo es que la respuesta en estado estacionario del sistema a una entrada sinusoidal es otra señal de tipo sinusoidal; la diferencia entre estas es solamente de magnitud y fase. Si conociéramos la forma en que afecta el sistema a una entrada sinusoidal de cualquier frecuencia podríamos determinar su respuesta a cualquier señal de entrada ya que todas las señales se pueden considerar como una combinación lineal de señales sinusoidales. Por eso resulta conveniente caracterizar a los sistemas con su respuesta en frecuencia, es decir con la información sobre el cambio que produce en la magnitud y en la fase de las señales sinusoidales de entrada antes de llevarlas a la salida.
DESARROLLO: Obtener la respuesta en frecuencia y la respuesta al impulso de los siguientes sistemas. Utilizar preferentemente la función FVTOOL
a)
b)
c)
CONCLUSIONES: Bibliografía [1] C. Burrus, J. McClelland, A. Oppenheim, T. Parks, R. Schafer, and H. Scuessler, editors. Computer-based exemses for Signal Processing using MATLAB. Prentice-Hall, 1994 [2] John Eaton. Octave Documentation. GNU, 1997 (En fotocopiadora). [6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing. Prentice-Hall, 1989. [4] A.V. Oppenheim, A.S. Willsky, and I.T. Young. Signals and Systems. Prentice-Hall, 1986.
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales Práctica 9 ANALISIS ESPECTRAL OBJETIVO: Utilizar la transformada rápida de Fourier para obtener la descomposición de una señal discreta en sus armónicos.
MATERIAL Y EQUIPO: Computadora con MATLAB
INTRODUCCIÓN: El análisis espectral se refiere al proceso de descomposición de una señal en sus componentes de frecuencia. Con este análisis se obtiene de cada componente de frecuencia una magnitud y una fase que representan lo que conocemos como transformada de Fourier. Para el caso de señales discretas se tiene la correspondiente transformada de Fourier de tiempo discreto (DTFT), la cual es una representación de la misma señal pero en el domino de la frecuencia discreta. Con la DTFT se obtiene una función continua de la frecuencia discreta q ue s e pue de obtener directamente de la expresión matemática que la define,
(1) Si la señal a analizar fuera de duración infinita seria imposible evaluar numéricamente la expresión anterior debido a las limitaciones de memoria en los equipos de cómputo. Sin embargo si la señal es de duración finita entonces la ecuación se puede calcular para cualquier valor de frecuencia. Desafortunadamente la cantidad de operaciones y los requerimientos de memoria aumentan de forma exponencial con el número de muestras. Para sobrepasar este inconveniente se definió una nueva transformada denominada transformada discreta de Fourier DFT que equivale al desarrollo en series de Fourier
para la señal a analizar. Para esto se supone que la señal representa solo un periodo de una señal ficticia de la cual se calcula su serie. Además se desarrollo la trasformada rápida de Fourier FFT que calcula la DFT mediante un la algoritmo que realiza las operaciones de forma eficiente. Existe una rutina en MATLAB “fft” que calcula la FFT y que puede usarse para analizar espectralmente una secuencia de duración finita. La función regresa el mismo número de datos que los que se ingresan. Si se utiliza una secuencia de 10 elementos la fft regresa 10 datos que representan los componentes de frecuencia de la secuencia original espaciados por 2π/N radianes, donde N representa el número de muestras de la señal. El orden en que la función “fft” de MATLAB entrega los componentes de frecuencia es diferente a como estamos acostumbrados a graficarlos (frecuencias positivas a la derecha y negativas a la izquierda), nos da los componentes de f recuencia negativa a la derecha después del ultimo componente de frecuencia positiva. Para ordenar los componentes de frecuencia en el orden acostumbrado se utiliza la función “fftshift”. Otras funciones comunes en el análisis espectral son “abs” y “angle” la primera para obtener el valor absoluto de una señal compleja y la segunda para la fase.
DESARROLLO: 1. Graficar el contenido espectral de la señal
x[n]
=
(
cos 0.1
n
)
usando la transformada rápida de Fourier. 2.
Generar una señal cuyas frecuencias espectral junto con la secuencia.
sea: 0.9*π graficar su contenido Y
3. Generar una señal con tres tonos diferentes sobrepuestos cuyas frecuencias son 0.1π, 0.3π 0.7π y graficar su contenido espectral junto con la secuencia. En este caso la expresión matemática para la señal a analizar es una suma de tres funciones coseno a las frecuencias mencionadas.. 4. Analizar el contenido espectral de la señal impulso unitario localizado en n=100 con 0
6. Analizar el contenido espectral de las siguientes señales -20
pulso sinusoidal
7.
un pulso rectangular con inicio en n=100 y fin en 150 con 0
CONCLUSIONES: Bibliografía [1] C. Burrus, J. McClelland, A. Oppenheim, T. Parks, R. Schafer, and H. Scuessler, editors. Computer-based exemses for Signal Processing using MATLAB. Prentice-Hall, 1994 [2] John Eaton. Octave Documentation. GNU, 1997 (En fotocopiadora). [6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing. Prentice-Hall, 1989. [4] A.V. Oppenheim, A.S. Willsky, and I.T. Young. Signals and Systems. Prentice-Hall, 1986.
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de
Procesamiento Digital de Señales Práctica 10 FILTROS DIGITALES DE RESPUESTA INFINITA AL IMPULSO IIR OBJETIVO: Comparar la respuesta en frecuencia de los filtros de respuesta infinita al impulso.
MATERIAL Y EQUIPO: Computadora con MATLAB
INTRODUCCIÓN: Un filtro electrónico es un elemento que discrimina una determinada frecuencia o gama de frecuencias de una señal que pasa a través de él, pudiendo modificar tanto su amplitud como su fase. Existen diferentes tipos de filtros, según la respuesta en frecuencia que se desee de estos, así que se puede requerir de filtros pasa bajas, pasa altas, pasa banda o rechazo de banda. Cualquiera de las respuestas en frecuencia deseadas se puede obtener con diferentes tipos de filtros que tiene características especiales en las bandas de paso, de transición o de rechazo. Los filtros más comunes son los que se muestran en la figura 1. El filtraje se puede obtener de forma analógica o digital.
Figura 1. Gráficas de respuesta en frecuencia de los principales tipos de filtros donde se muestra sus diferencias en las bandas de paso, de rechazo y de transición. DESARROLLO: Comparar la respuesta en frecuencia de los diferentes tipos de filtros digitales IIR. 1. Para un filtro pasa bajas de quinto orden, una frecuencia de corte de 0.5 radianes y rizo de 1 dB (donde aplique). 2. Para un filtro pasa banda de orden 10, frecuencias de corte en 0.3 y 0.6 radianes y rizo de 1 dB. 3. Comparar la respuesta en frecuencia del filtro Chebyshev tipo 1 pasa bajas para una frecuencia de 0.3 rad y rizos de 1 dB, 3 dB y 6 dB. 4. Comparar la respuesta en frecuencia del filtro Chebyshev tipo 2 pasa altas para una frecuencia de corte de 0.3 rad y rizos de 1 dB, 3 dB y 6 dB. 5. Comparar la respuesta en frecuencia del filtro Elíptico pasa bajas para una frecuencia de corte de 0.3 rad y rizos de 1 dB, 3 dB y 6 dB en banda de paso y de rechazo.
CONCLUSIONES: Bibliografía [1] C. Burrus, J. McClelland, A. Oppenheim, T. Parks, R. Schafer, and H. Scuessler, editors. Computer-based exemses for Signal Processing using MATLAB. Prentice-Hall, 1994 [2] John Eaton. Octave Documentation. GNU, 1997 (En fotocopiadora). [6] A.V. Oppenheim and R.W. Schafer, Discrete-Time Signal Processing. Prentice-Hall, 1989. [4] A.V. Oppenheim, A.S. Willsky, and I.T. Young. Signals and Systems. Prentice-Hall, 1986.
UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA, PROGRAMA DE ELECTRÓNICA
Laboratorio de Procesamiento Digital de Señales
Práctica 10 FILTROS DIGITALES DE RESPUESTA FINITA AL IMPULSO FIR OBJETIVO: Comparar la respuesta en frecuencia de los diferentes tipos de filtros de respuesta finita al impulso.
MATERIAL Y EQUIPO: Computadora con MATLAB
INTRODUCCIÓN: Los filtros digitales de respuesta finita al impulso (FIR por sus siglas en ingles) son filtros que obtiene la señal de salida usando solamente valores de la secuencia de entrada y ningún valor de la señal de salida, por eso su respuesta a un impulso en la entrada tiene una longitud finita y por lo tanto siempre es estable. Sin embargo tienen la desventaja de necesitar un orden mayor respecto a los filtros IIR para cumplir las mismas características de respuesta en frecuencia.
DESARROLLO: Comparar la respuesta en frecuencia de los filtros digitales que aplican el método de la ventana, la aproximación por mínimos cuadrados y muestreo en frecuencia. 1. Para un filtro pasa bajas de orden 15, una frecuencia de corte de 0.5 radianes. 2. Para un filtro pasa banda de orden 50, frecuencias de corte en 0.3 y 0.6 radianes. 3. Comparar la respuesta en frecuencia del filtro Chebyshev tipo 1 pasa bajas para una frecuencia de 0.3 rad y rizos de 1 dB, 3 dB y 6 dB. 4. Comparar la respuesta en frecuencia del filtro Chebyshev tipo 2 pasa altas para una frecuencia de corte de 0.3 rad y rizos de 1 dB, 3 dB y 6 dB. 5. Comparar la respuesta en frecuencia del filtro Elíptico pasa bajas para una frecuencia de corte de 0.3 rad y rizos de 1 dB, 3 dB y 6 dB en banda de paso y de rechazo.
A continuación se dan algunos ejemplos para que se realicen la práctica 10
Fundamentación Teórica El filtro digital es la implementación en hardware o software de una ecuación en diferencias con una entrada digital.
Respuesta al impulso (convolución) del filtro:
Ejemplo: Respuesta al impulso de un filtro con coeficientes a(1)=1, a(2)=-0.9, b(1)=1 con MatLab. n=0:49; %señal impulso imp = [1; zeros(49,1)]; %coeficientes del filtro b=1; a=[1 -0.9]; %respuesta al impulso h = filter(b,a,imp); stem(n,h)
Respuesta en frecuencia del filtro:
Ejemplo: % coeficientes del filtro [b,a] = cheby1(12,0.5,200/500); %respuesta en frecuencia [h,f] = freqz(b,a,256,1000); %grafica de la magnitud mag=abs(h); subplot(121) plot(f,m) %grafica de la fase fase=unwrap(f*180/pi); subplot(122) plot(f,fase)
Los filtros digitales tienen: • • • •
Alta inmunidad al ruido Alta precisión, limitada por los errores de redondeo en la aritmética empleada Fácil modificación de las características del filtro Muy bajo costo
Los filtros se clasifican en filtros FIR (Respuesta impulsional finita) y filtros IIR (Respuesta impulsional Infinita)
FILTROS FIR Un filtro FIR de orden M tiene la siguiente función de transferencia:
Tiene como función de transferencia:
• •
•
La secuencia b(k) son los coeficientes del filtro Es no recursivo, o sea, la salida depende solamente de las entradas y no de las salidas pasadas La función de transferencia sólo tiene ceros, excelente estabilidad.
FILTROS IIR Tiene como ecuación en diferencias:
Tiene como función de transferencia:
•
• • •
Es recursivo, o sea, que su salida además de las entradas depende de las salidas pasadas. Tiene polos y ceros, tiene problemas de estabilidad. La fase no es lineal con la frecuencia El orden del filtro es mucho menor que un filtro FIR para la misma aplicación
DISEÑO DE FILTROS DIGITALES El diseño consiste en obtener los coeficientes del filtro para conseguir unos requerimientos específicos. Su implementación obedece en escoger y aplicar a una estructura particular del filtro esos coeficientes. Los filtros se normalizan a la frecuencia de Nyquist, o sea, a la frecuencia de muestreo dividida por dos:
Por ejemplo, para filtrar 30 Hz con un filtro pasabajas y fs =100 Hz con un Butterworth de orden 5: [b,a] = butter(5,30/50) = butter(5,0.6) Para convertir la frecuencia normalizada a frecuencia angular se debe multiplicar por π. Una especificación más rigurosa podría ser riple en la banda de paso (passband Rs), atenuación en la banda de rechazo (stopband-Rp) o en la banda de transición (ws-wp), etc.
1. DISEÑO DE FILTROS IIR Hay varios métodos:
a) IIR USANDO FILTROS ANÁLOGOS Filtro Butterworth Comprende diseños de filtros pasabajo, pasa banda, pasa alto, y banda rechazo. La respuesta en magnitud es plana en la banda de paso. Filtro pasabajo de orden n con frecuencia de corte wn
[b,a] = butter(n,wn) [b,a] = butter(n,wn, ‘ftype’) [z,p,k] = butter(n,wn) [z,p,k] = butter(n,wn, ‘ftype’) Tipo= high, low, bandpass, bandstop (wn=[w1 w2]) La frecuencia de corte es donde la magnitud del filtro es Ω=1 Para encontrar el orden y la frecuencia de corte dadas las especificaciones: [n,Wn] = buttord(Wp,Ws,Rp,Rs)
en
Ejemplo: Pasa banda Se quiere un filtro pasabanda de 100 a 2000 Hz, la banda stop arranca en 500 Hz, la frecuencia frecuencia de muestreo es de 10 KHz, al menos 1 dB de riple en la banda de paso y al menos 60 dB de atenuación en la banda stop. [n,Wn] = buttord([1000 2000]/5000,[500 2500]/5000,1,60) %n = 12, Wn = 0.1951 0.4080 [b,a] = butter(n,Wn); [sos,g] = tf2sos(b,a); Hd = dfilt.df2tsos(sos,g); h = fvtool(Hd); set(h,'Analysis','freq')
Filtro Chebyshev Tipo I Minimiza la diferencia entre el ideal y la respuesta de frecuencia actual sobre la banda banda de de paso paso incorpor incorporand ando o un un equir equiriple iple de Rp dB en la banda banda de paso. paso. La respuesta en la banda rechazo es plana (maximally flat). La transición de la banda de paso a la banda de rechazo es más rápida que en el de Butterworth.
[z,p,k] = cheby1(n,R,Wp) : Filtro pasa bajo z,p,k] = cheby1(n,R,Wp,'ftype') [b,a] = cheby1(n,R,Wp) : Filtro pasa bajo [b,a] = cheby1(n,R,Wp,'ftype') El orden del filtro es n con frecuencia de corte en la banda de paso normalizada en Wp y R dB de riple pico a pico en la banda de paso.
Ejemplo: Pasa bajo Para una frecuencia de muestreo de 1000 Hz diseñar un filtro Chebyshev Chebyshev Tipo I con 0.5 dB de riple en la banda de paso y una frecuencia de borde en la banda de paso de 300 Hz. [z,p,k] = cheby1(9,0.5,300/500); [sos,g] = zp2sos(z,p,k); Hd = dfilt.df2tsos(sos,g); h = fvtool(Hd) set(h,'Analysis','freq')
%Respuesta en frecuencia freqz(b,a,512,1000)
Filtro Chebyshev Tipo II Minimiza la diferencia con el filtro ideal en la banda stop incorporando un equiriple de Rs dB en la banda stop. La respuesta en la banda de paso es plana (Maximally flat).
[z,p,k] = cheby2(n,R,Wst) [z,p,k] = cheby2(n,R,Wst,'ftype') [b,a] = cheby2(n,R,Wst) [b,a] = cheby2(n,R,Wst,'ftype')
Ejemplo: Pasa bajo Para una frecuencia de muestreo de 1000 Hz diseñe un filtro pasa bajo Chebyshev II con atenuación en la banda stop de 20 dB debajo de la banda de paso y una frecuencia de borde en banda stop de 300 Hz. [z,p,k] = cheby2(9,20,300/500); [sos,g] = zp2sos(z,p,k); % Convert to SOS form Hd = dfilt.df2tsos(sos,g); % Create a dfilt object h = fvtool(Hd); % Plot magnitude response set(h,'Analysis','freq') % Display frequency response
Filtro Elíptico Es un filtro equiriple tanto en la banda de paso como en la banda de rechazo. Riple en la banda de paso Rp, riple en la banda stop Rs. Minimiza el ancho de la transición.
[z,p,k] = ellip(n,Rp,Rs,Wp) [z,p,k] = ellip(n,Rp,Rs,Wp,'ftype') [b,a] = ellip(n,Rp,Rs,Wp) [b,a] = ellip(n,Rp,Rs,Wp,'ftype') Wp frecuencia normalizada en banda de paso, Rp riple en dB en la banda de paso, Rs riple en dB en la banda rechazo.
Ejemplo: Pasa bajo Diseñar un filtro pasa bajo Elíptico de orden 6 con fp=300 Hz, 3 dB en la banda de paso y 50 dB de atenuación en la banda rechazo.
[z,p,k] = ellip(6,3,50,300/500); [sos,g] = zp2sos(z,p,k); % Convert to SOS form Hd = dfilt.df2tsos(sos,g); % Create a dfilt object h = fvtool(Hd) % Plot magnitude response set(h,'Analysis')
RESUMEN .
Ejemplos:
b) DISEÑO DE IIR EN FORMA DIRECTA Se diseña en forma directa especificando la respuesta en frecuencia. El método encuentra la transformada inversa FFT y la resuelve utilizando la ecuación Yule – Walker. [b,a] = yulewalk(n,f,m) La frecuencia f es un vector de 0 a 1, donde 1 representa la frecuencia de Nyquist. La magnitud m es un vector que contiene la respuesta de la magnitud deseada en los puntos de f.
Ejemplo: Filtro multibanda de orden 10 m = [0 0 1 1 0 0 1 f = [0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1]; [b,a] = yulewalk(10,f,m); [h,w] = freqz(b,a,128); Plot(f,m,w/pi,abs(h))
Ejemplo: Filtro pasa bajo de orden 8 f = [0 0.6 0.6 1]; m = [1 1 0 0]; [b,a] = yulewalk(8,f,m); [h,w] = freqz(b,a,128);
1
0 0];
plot(f,m,w/pi,abs(h),'--') legend('Ideal','Diseño yulewalk ') title('Comparación de la respuesta en frecuencia')
2. DISEÑO DE FILTROS FIR Los filtros FIR tienen las siguientes ventajas: -
Tienen fase lineal Son siempre estables Eficientes realizaciones en hardware Transientes de duración finita
La principal desventaja es que el orden del filtro para una similar respuesta es mucho más alto que la de un filtro IIR.
Respuesta en frecuencia:
A(w) es real si:
Simetría par: K1 = 0
Si N es impar: Tipo I
Si N es par: Tipo II
Simetría Impar:
Si N es Par: Tipo IV
MÉTODOS DE DISEÑO
a) MÉTODO DE VENTANEO La truncación de una secuencia en el dominio del tiempo causa el fenómeno de Gibbs como una discontinuidad en el dominio de la frecuencia. Si hd(n) es la secuencia del prototipo ideal, si se trunca por una ventana rectangular:
La multiplicación de dos funciones en el dominio del tiempo corresponde a la convolución de sus Transformadas de Fourier,
Ventana triangular de Bartlett:
Ventanas de coseno generalizado: Rectangular, Hanning, Hamming y Blackman Ventana Rectangular Hanning Hamming Blackman Ventana Kaiser:
a 1.0 0.5 0.54 0.42
b 0.0 0.5 0.46 0.5
c 0.0 0.0 0.0 0.8
La ventana rectangular es la mas simple pero presenta el fenómeno de Gobbs. La de Bartlett o triangular reduce el overshoot pero dispersa considerablemente la banda de transición. Las hanning, Hamming y Blackman proveen una truncación mas suave y una respuesta de frecuencia mejor. La mejor puede ser la de kaiser que permite ajustar el compromiso entre overshoot y banda de transición con el parámetro β. Considérese un filtro digital ideal pasa bajo con frecuencia de corte wo en rad/seg. Este filtro tiene una magnitud igual a 1 en frecuencias menores a wo y magnitud igual a 0 en frecuencias entre wo y π. Su respuesta al impulso es:
Este filtro no es implementable porque su respuesta al impulso es infinito y no causal. Para crear una respuesta de duración finita hay que truncarlo por medio de una ventana. Ejemplo: Filtro pasabajo de orden 51 con frecuencia de corte wo=0.4π rad/seg, tiene coeficientes: %longitud 51 ventana rectangular n=-25:25; %coeficientes b=0.4*sinc(0.4*n); fvtool(b,1)
Se presenta el efecto Gibbs en la banda de paso, esta distorsión se disminuye si se aplica una ventana tipo haming:
%longitud 51 ventana Hamming n=-25:25; %coeficientes b=0.4*sinc(0.4*n); b = b.*hamming(51)'; fvtool(b,1)
Las Funciones fir1 y fir2 de Matlab basan el proceso en el uso de ventanas. Dado el orden del filtro y la descripción del filtro deseado, estas funciones retornan la Transformada de Fourier Inversa del filtro con ventana. La multiplicación de una ventana en el dominio del tiempo causa una convolución en el dominio de la frecuencia. VENTANAS En diseño de filtros digitales se escoge una ventana para amortiguar los efectos Gibbs que resulta de la truncación de una señal infinita en el tiempo.
EJEMPLOS: Ventana Rectangular Ventana rectangular de longitud de 50 %ventana rectangular n = 50; w = rectwin(n); wintool
Ventana Bartlett %ventana de Bartlett de 8 n = 8; w = bartlett(n); wintool
%w =[0;0.2857;0.5714;0.8571;0.8571;0.5714;0.2857;0]
Ventanas cosenoidales Las ventanas Blackman, Flat Top, Hamming, Hann, and rectangular son casos especiales de de ventanas de coseno. Son combinaciones de secuencias senoidales con frecuencias de 0, 2π/N-1, 4π/N-1 donde N es el número de puntos de la ventana. ind = (0:n-1)' *2 *pi/(n-1); w = A - B* cos(ind) + C* cos(2*ind); Para N=64:
Ventana Kaiser La ventana de Kaiser es una aproximación para maximizar la energía del lóbulo principal frente a los lóbulos laterales. El parámetro β controla este peso del lóbulo principal. n = 50; w1 = kaiser(n,1); w2 = kaiser(n,4); w3 = kaiser(n,9); [W1,f] = freqz(w1/sum(w1),1,512,2); [W2,f] = freqz(w2/sum(w2),1,512,2); [W3,f] = freqz(w3/sum(w3),1,512,2); plot(f,20*log10(abs([W1 W2 W3]))); grid; legend('beta = 1','beta = 4','beta = 9',3)
Para encontrar el orden del filtro se utiliza la función: [n,Wn,beta,ftype] = kaiserord(f,a,dev) f: frecuencia de corte de las bandas a: amplitud deseada en las bandas dev: especifica el riple de la banda de paso y la atenuación de la banda stop, en forma de ganancia absoluta no en dB Para calcular los coeficientes del filtro se usa la función fir1: b = fir1(n,Wn,kaiser(n+1,beta),ftype,'noscale')
ftype: es ‘high’ para pasa alto y ‘stop’ para banda stop. Para pasa banda puede ser ‘dc-0’ so la primera banda es banda stop o ‘dc-1’si es pasa banda. Algoritmo
Donde α es la atenuación en banda stop en dB
Donde n es el orden del filtro y Δw es la región de transición más pequeña. EJEMPLOS: Diseñe un filtro pasabajo con banda de paso definida de 0 a 1 KHz y banda stop de 1500 Hz a 4 KHz. El riple de la banda de paso es del 5% y la atenuación en la banda stop de 40 dB. %orden del filtro kaiser fs = 8000; fc = [1000 1500]; mag = [1 0]; %5% =0.05, -40dB =20 log(0.01) dev = [0.05 0.01]; [n,Wn,beta,ftype] = kaiserord(fc,mag,dev,fs); h = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); freqz(h) %n = 36, Wn = 0.3125, ftype =low, beta = 3.3953
Ejemplo: Diseñar un filtro pasa banda de longitud impar, usar fir1 %filtro de longitud impar (orden par) fs = 8000; fc = [1000 1300 2210 2410]; mag = [0 1 0]; dev = [0.01 0.05 0.01]; [n,Wn,beta,ftype] = kaiserord(fc,mag,dev,fs); n = n + rem(n,2); h = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); [H,f] = freqz(h,1,1024,fs); plot(f,abs(H)), grid on %n = 90, beta = 3.3953, ftype = DC-0
Ejemplo: Diseñar un filtro pasabajo con corte en la banda de paso de 1500 Hz y corte en banda stop de 2000 Hz, riple en banda de paso de 0.01 y en banda stop de 0.1, frecuencia de muestreo de 8000 Hz. [n,Wn,beta,ftype] = kaiserord([1500 2000],[1 0],[0.01 0.1],8000); b = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale'); %n = 36, Wn = 0.4375, ftype = low USO DE LA FUNCIÓN: fir1 Para diseño de filtros estándar pasa bajo, pasa alto, pasa banda y banda stop se utiliza fir1
Ejemplo: n = 50; Wn = 0.4; %por defecto es pasa bajo con ventana de hamming b = fir1(n,Wn); En general: b = fir1(n,Wn,'ftype',window) •
• •
Wn es un número entre 0 y 1, donde 1 corresponde a la frecuencia de Nyquist. Si Wn = [ w1 w2] fir1 retorna un filtro pasa banda con w1
El tipo de filtro se programa con ftype: ‘high’ filtro pasa alto ‘stop’ para un banda stop, si wn=[w1 w2] es el rango de frecuencias de la banda stop ‘dc-1’ la primera banda del filtro multibanda es pasa banda ‘dc-0’ la primera banda del filtro multibanda es un banda stop Si h(n) es la Transformada de Fourier Inversa de la respuesta en frecuencia ideal y w(n) es la ventana, los coeficientes del filtro son: • •
• •
B(n)= w(n) h(n),
1≤n≤N
EJEMPLOS Diseñe un filtro FIR pasa banda de orden 48 con pasa banda de 0.35 ≤ w ≤ 0.65 % uso de fir1 b = fir1(48,[0.35 0.65]); freqz(b,1,512)
Diseñe un filtro FIR pasa alto para que atenúe las frecuencias posteiores a fc = 0.48 y ventana Chebyshev de 30 dB de riple. %pasa alto b = fir1(34,0.48,'high',chebwin(35,30)); freqz(b,1,512)
USO DE LA FUNCIÓN: fir2 La función fir2 también diseña filtros FIR ventaneados pero con respuesta de frecuencia lineal arbitraria. En contraste con fir1 que solamente diseña filtros estándar. n = 50; f = [0 .4 .5 1]; m = [1 1 0 0]; b = fir2(n,f,m); En general: b = fir2(n,f,m,window) Ejemplo: Filtro pasa bajo de orden 30. %uso de fir2 f = [0 0.6 0.6 1]; m = [1 1 0 0]; b = fir2(30,f,m); [h,w] = freqz(b,1,128); plot(f,m,w/pi,abs(h)) legend('Ideal','diseño con fir2') title('comparación de magnitudes en respuesta en frecuencia')
b) MÉTODO DE MULTIBANDA CON BANDAS DE TRANSICIÓN LEAST SQUARED ERROR Las especificaciones de los filtros están dada generalmente en el dominio de la frecuencia y puesto que la energía de la señal está relacionada con el cuadrado de la señal, el criterio de error cuadrático es el indicado. Un método es considerar la integral del cuadrado del error (Integral squared error).
Se trata es de minimizar este error. La DFT es:
La IDFT es:
Diferenciador Un diferenciador ideal puede ser un filtro de fase lineal FIR. La respuesta en frecuencia del diferenciador es:
Los coeficientes del filtro son:
Banda de transición: Otra forma es minimizar el error aplicándolo a la banda de transición:
Uso de pesos: Introducir una función de peso W(w) al error.
APROXIMACIÓN DE CHEBYSHEV Minimiza el valor del máximo error y tienen un equiriple en el comportamiento de la respuesta en frecuencia.
El algoritmo de Remes Exchange desarrollado por Parks y McClellan, hace que la función de error tome valores de ±δ para un conjunto de r+1 frecuencias fm , m=1,…, r+1.
Función: firls (Least Squares) Minimiza el error cuadrático medio entre la respuesta de frecuencia deseada y la respuesta de frecuencia actual.
b = firls(n,f,a) b(k) = b(n+2-k)
f: es un vector de parejas de puntos de frecuencias entre 0 y 1, donde 1 corresponde a la frecuencia de Nyquist. a: Es un vector que contiene la amplitud deseada en los puntos especificados de f En forma general: b = firls(n,f,a,w,'ftype') w: es el peso en la banda de frecuencia ftype: ‘hilbert’ para filtros de fase lineal con simetría impar (tipo III y tipo IV) b(k) = - b(n+2-k). Usa la Transformada de Hilbert “differentiator’ usa técnica especial de pesos en las bandas Ejemplo: Diseñe un filtro pasa bajo de orden 255 con banda de transición: %uso firls b = firls(255,[0 0.25 0.3 1],[1 1],[ 1 1 0 0]); Usando fdatool
Diseñar un filtro antisimétrico de orden 24 %filtro antisimétrico de orden 24 F = [0 0.3 0.4 0.6 0.7 0.9]; A = [0 1 0 0 0.5 0.5]; b = firls(24,F,A,'hilbert'); for i=1:2:6, plot([F(i) F(i+1)],[A(i) A(i+1)],'r'), hold on end [H,f] = freqz(b,1,512,2); plot(f,abs(H)), grid on, hold off legend('Ideal','Diseño firls ')
Usando fdatool:
En resumen esta function function diseña filtros de fase lineal tipo I, II, III, IV. Los tipo I y II son por defecto para n par e impar respectivamente. Los ‘hilbert’ y ‘differentiator’ producen tipo III (n par) y IV (n impar) Función: firpm (Parks-McClellan) Esta function function implementa implementa el algorit algoritmo mo de Parks-McClell Parks-McClellan an que usa la teoría teoría de Remez y la aproximación de Chebyshev. Minimizan el máximo error entre la respuesta de frecuencia deseada y la actual. Son filtros equiriple. %uso de firpm n = 20; % orden del filtro filtro f = [0 0.4 0.4 0.5 0.5 1]; 1]; % bord bordes es de la band bandaa de de frec frecue uenc nciias a = [1 1 0 0]; % amplitudes deseadas b = firpm(n,f,a); bb = firls(n,f,a); fvtool(b,1,bb,1) legend('firpm','firls')
Nótese que con firls presenta mejor respuesta en la banda de paso y en la banda stop, pero firpm es mejor en la banda de transición. Mediante trazos de líneas firls y firpm pueden realizar filtros, pasa bajo, pasa alto, banda stop, psasa banda, multibanda. %tipos de filtros con firls o firpm %pasa banda f = [0 0.3 0.4 0.7 0.8 1]; a = [0 0 1 1 0 0]; %pasa alto a = [0 0 1 1]; f = [0 0.3 0.4 0.5 0.8 1]; %banda stop a = [1 1 0 0 1 1]; %multibanda f = [0 0.1 0.15 0.25 0.3 0.4 0.45 0.55 0.6 0.7 0.75 0.85 0.9 1]; a = [1 1 0 0 1 1 0 0 1 1 0 0 1 1];
Uso del vector peso Tanto firls como firpm permiten minimizar el error en ciertas bandas de frecuencias especificando un vector peso. Por ejemplo, un filtro equiriple pasa bajo con riple 10 veces menor en la banda stop que en la banda de paso es: %uso de pesos en las bandas n = 20; % orden del fitro f = [0 0.4 0.5 1]; % bordes de la bandas de frecuencias a = [1 1 amplitudes deseadas w = [1 10]; % vector peso b = firpm(n,f,a,w);
0
0];
%
fvtool(b,1)
Filtros antisimétricos/Transformada de Hilbert. Para filtros de simetría impar Tipo III (orden par) o Tipo IV (orden impar) %filtros antisimétricos b = firpm(21,[0.05 1],[1 1],'h'); % Pasa alto Hilbert bb = firpm(20,[0.05 0.95],[1 1],'h'); % pasa banda Hilbert fvtool(b,1,bb,1) legend('pasa alto', 'pasa banda')
La señal analítica de una señal compleja x está compuesta por la parte real de x y la Transformada de Hilbert como la parte imaginaria (orden par).
%señal analítica fs = 1000; % frecuencia de muestreo t = (0:1/fs:2)'; % vector tiempo x = sin(2*pi*300*t); % señal senoidal de 300 Hz xh = filter(bb,1,x); % Transformada Hilbert de x xd = [zeros(10,1); x(1:length(x)-10)]; % 10 muestras de retardo xa = xd + j*xh; % señal analítica
Diferenciadores La diferenciación de una señal en el dominio del tiempo es equivalente a la multiplicación de su transformada de Fourier por una función rampa. Diferenciar una señal es pasarla por un filtro H(w)= jw %diferenciador b = firpm(21,[0 1],[0 pi],'d'); bb = firpm(20,[0 0.9],[0 0.9*pi],'d'); fvtool(b,1,bb,1) legend('diferenciador orden impar','diferenciador orden par')
c) MÉTODO DEL MÍNIMO CUADRADO RESTRINGIDO Constrained Least Squares (CLS) La restricción consiste en que no es necesario definir las bandas de transición. Permite considerar umbrales más altos o más bajos con riple máximo permisible. Función fircls1 Se usa para el diseño de pasa bajo o pasa alto. b = fircls1(n,wo,dp,ds) b = fircls1(n,wo,dp,ds,'high')
wo: es la frecuencia de corte normalizada dp: riple en la banda de paso ds: riple en la banda stop %uso de fircls1 n = 55;wo = 0.3; dp = 0.02; ds = 0.008; b = fircls1(n,wo,dp,ds); fvtool(b)
Usando pesos en las bandas: b = fircls1(n,wo,dp,ds,wp,ws,k) Donde k es la relación entre el error en la banda de paso y el error en la banda stop
%Uso de pesos en fircls1 n = 55; wo = 0.3; dp = 0.02; ds = 0.004; wp = 0.28; ws = 0.32; k = 10; h = fircls1(n,wo,dp,ds,wp,ws,k); fvtool(h,1)
Función fircls Se usa para filros FIR multibanda en el cual se pueden especificar la cantidad máxima de riple en cada banda. b = fircls(n,f,amp,up,lo) f: vector de frecuencias de transición amp: magnitudes up, lo: define los límites alto y bajo para cada banda
%uso de fircls n=150; f=[0 0.4 1]; a=[1 0]; up=[1.02 0.01]; lo =[0.98 -0.01]; %con despliegue de las bandas b = fircls(n,f,a,up,lo,'both'); fvtool(b)
n = 129; f = [0 0.3 0.5 0.7 0.9 1]; a = [0 0.5 0 1 0]; up = [0.005 0.51 0.03 1.02 0.05]; lo = [-0.005 0.49 -0.03 0.98 -0.05]; h = fircls(n,f,a,up,lo); fvtool(h,1)
IMPLEMENTACIÓN DE FILTROS DIGITALES Después de haber diseñado el filtro, la función de transferencia o ecuación en diferencias del filtro debe implementarse y hay varias maneras de hacerlo. Estas maneras de realizarlas son representadas en diagramas en bloques y son llamadas estructuras del filtro y se escoge aquella que sea menos sensitiva a los errores en los coeficientes o que minimice el ruido introducido por la cuantización de los coeficientes. Función dfilt, filter Procedimiento para la implementación: 1. Se generan los coeficientes del filtro IIR o FIR por el método de diseño apropiado 2. Crear el objeto del filtro de los coeficientes usando dfilt 3. Aplicar al objeto del filtro al dato x usando la función filter Por ejemplo, diseñar, implementar como estructura forma directa traspuesta II y aplicar a un filtro Butterworth %implementación de filtros [b,a] = butter(5,0.4); Hd = dfilt.df2t(b,a); % Forma directa II transpuesta fs=100;
t=0:1/fs:1; x=sin(2*pi*t*3)+0.25*sin(2*pi*t*40); y=filter(Hd,x); plot(t,x,t,y) legend('x','y')
ESTRUCTURA DE FILTROS 1. IMPLEMENTACIÓN DE FILTROS IIR a) Forma Directa I: dfilt.df1
%forma directa I [b,a] = butter(4,.5); Hd = dfilt.df1(b,a); num=get(Hd,'Numerator') %num = 0.0940 0.3759 0.5639 0.3759 0.0940 den=get(Hd,'Denominator') %den = 1.0000 -0.0000 0.4860 -0.0000 0.0177 fvtool(b,a)
b) Forma Directa I sos (second order section): dfilt.df1sos El filtro es implementado en cascada por secciones de Segundo orden (tres coeficientes en el numerador y tres coeficientes en el denominador).
Ejemplo: para un filtro ellíptico , pasa bajode orden 6, Rs=1, Rp=60 dB, Wp=0.4 %forma directa I sos [b,a] = ellip(6,1,60,.4); % coeficientes del filtro % s: coeficientes g: ganancia de la sección
[s,g] = tf2sos(b,a);
% Convierte a SOS
Hd=dfilt.df1sos(s,g); get(Hd,'sosmatrix')
fvtool(Hd)
c) Forma Directa I Transpuesta: dfilt.df1t
Ejemplo: [b,a] = butter(4,.5); Hd = dfilt.df1t(b,a); num=get(Hd,'Numerator') %num = 0.0940 0.3759 0.5639 0.3759 0.0940 den=get(Hd,'Denominator') %den = 1.0000 -0.0000 0.4860 -0.0000 0.0177 fvtool(b,a) d) Forma Directa I Transpuesta sos: dfilt.df1tsos
%forma directa I sos [b,a] = ellip(6,1,60,.4); % coeficientes del filtro % s: coeficientes g: ganancia de la sección [s,g] = tf2sos(b,a); % Convierte a SOS Hd=dfilt.df1tsos(s,g); get(Hd,'sosmatrix')
e) Forma Directa II: dfilt.df2
Ejemplo: [b,a] = butter(4,.5); Hd = dfilt.df2(b,a) coeffs(Hd) realizemdl(Hd) f) Forma Directa II sos: dfilt.df2sos
Ejemplo: [z,p,k] = ellip(6,1,60,.4); [s,g] = zp2sos(z,p,k);
Hd = dfilt.df2sos(s,g) g) Forma Directa II Transpuesta: dfilt.df2t
Ejemplo: [b,a] = butter(4,.5); Hd = dfilt.df2t(b,a) h) Forma Directa II Transpuesta sos: dfilt.df2tsos
[z,p,k] = ellip(6,1,60,.4); [s,g] = zp2sos(z,p,k); Hd = dfilt.df2tsos(s,g)
2. IMPLEMENTACIÓN DE FILTROS FIR a) Forma directa FIR: dfilt.dffir
b = firpm(30,[0 .1 .2 .5]*2,[1 1 0 0]); Hd = dfilt.dffir(b) fvtool(Hd) b) Forma directa transpuesta FIR: dfilt.dffirt
b = firpm(30,[0 .1 .2 .5]*2,[1 1 0 0]); Hd = dfilt.dffirt(b) c) Forma directa simétrica FIR: dfilt.dfsymfir
Ejemplo: Orden Par b = [-0.008 0.06 0.44 0.44 0.06 -0.008]; Hd = dfilt.dfsymfir(b) fvtool(Hd) Ejemplo: Orden Impar b = [-0.01 0.1 0.8 0.1 -0.01]; Hd = dfilt.dfsymfir(b) fvtool(Hd) d) Forma directa antisimétrica FIR: dfilt.dfasymfir
Ejemplo: %orden impar tipo IV Hd = firpm(25,[0 .4 .5 1],[0 0 1 1],'h'); fvtool(Hd) %orden par tipo III h=firpm(44,[0 .3 .4 1],[0 .2 0 0],'differentiator'); fvtool(Hd)
e) Forma lattice ARMA (autoregressive, moving-average): dfilt.latticearma