Guia Laboratorio PDS UTP

EXPERIENCIAS DE LABORATORIO

Asignatura : Ciclo : Software : Experiencias : Docente : Ciclo académico :

Procesamiento Digital de Señales IX MATLAB 09 Pedro Freddy Huamaní Navarrete 2011 - III

Carrera de Ingeniería Electrónica FACULTAD DE MECATRÓNICA

INGENIERÍA

ELECTRÓNICA

Y

JULIO 2011

Laboratorio de Procesamiento Digital de Señales

2011-III

EXPERIENCIAS DE LABORATORIO PARA LA ASIGNATURA DE PROCESAMIENTO DIGITAL DE SEÑALES EXPERIENCIA DE LABORATORIO 1ra Experiencia de Laboratorio 2da Experiencia de Laboratorio 3ra Experiencia de Laboratorio 4ta Experiencia de Laboratorio 5ta Experiencia de Laboratorio EXAMEN PARCIAL 6ta Experiencia de Laboratorio

7ma Experiencia de Laboratorio

8va Experiencia de Laboratorio 9na Experiencia de Laboratorio

SEMANA DE CLASES

Semana 01 Semana 02 Semana 03 Semana 04 Semana 05 Semana 06 Semana 07 Semana 08 Semana 09 -------------------Semana 11 Semana 12 Semana 13 Semana 14 Semana 15 Semana 16 Semana 17 Semana 18

TEMA DE LABORATORIO Introducción al Toolbox Signa Processing del Matlab Gráfico de señales periódicas y no periódicas. periódicas. Convolución. Cambio de la tasa de muestreo: Decimación e Interpolación. Transformada discreta de Fourier: Directa e Inversa Transformada discreta de Fourier corta en el tiempo. Espectrograma. -------------------Transformada Z. Diagrama de polos y ceros. Estructura de filtros digitales. Filtros digitales no recursivos. Pasa-Bajo, Pasa-Alto, Pasa-Banda y Rechaza-Banda. Filtros digitales recursivos. Transformación de frecuencia. Introducción al filtrado adaptivo: algoritmo Least Mean Square.

EXAMEN FINAL

--------------------

--------------------

EXAMEN SUSTITUTORIO

--------------------

--------------------

MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo Completo – FIEM FIEM – UTP UTP

Página 2


2011-III

ÍNDICE

1. Experiencia de Laboratorio Nº 01: Introducción al Toolbox Signal Processing del Matlab. Pág. 03 2. Experiencia de Laboratorio Nº 02: Gráfico de señales periódicas y no periódicas. Convolución. Convolución. Pág. 10 3. Experiencia de Laboratorio Nº 03: Cambio de la tasa de muestreo: Decimación e Interpolación. Pág. 19 4. Experiencia Experiencia de Laboratorio Laboratorio Nº Nº 04: Transformada discreta de Fourier. Directa e Inversa.

Pág. 24

5. Experiencia Experiencia de Laboratorio Laboratorio Nº Nº 05: Transformada discreta de Fourier corta en el tiempo. Espectrograma. Espectrograma.

Pág. 32

6. Experiencia Experiencia de Laboratorio Laboratorio Nº 06: 06: Transformada Z. Diagrama de Polos y Ceros. Ceros. Estructura de Filtros Digitales. Pág. 39 7. Experiencia Experiencia de Laboratorio Laboratorio Nº 07: 07: Filtros digitales digitales no recursivos. recursivos. Pasa-Bajo, Pasa-Bajo, Pasa Alto, Pasa-Banda y Rechaza-Ba Rechaza-Banda. nda. Pág. 47 8. Experiencia de Laboratorio Nº 08: Filtros digitales recursivos. Transformació Transformación n de frecuencias. Pág. 59 9. Experiencia Experiencia de Laboratorio Laboratorio Nº Nº 09: Introducción Introducción al filtrado adaptativo. Algoritmo Least Least Mean Square


Pág. 65

Página 3


2011-III

EXPERIENCIA DE LABORATORIO Nº 01

INTRODUCCIÓN AL TOOLBOX SIGNAL PROCESSING DEL MATLAB

1. OBJETIVOS: 1.1 Dar una introducción al toolbox de procesamiento de señales del software Matlab. 1.2 Realizar operaciones con vectores para representar señales digitales en función al tiempo. 1.3 Realizar operaciones con matrices para representar imágenes digitales de distintas dimensiones. 2.- FUNDAMENTO TEÓRICO: 2.1

Toolbox de Procesamiento de Señales del Matlab El Matlab, en su versión R2007b, cuenta con el toolbox Signal Processing versión 6.8 de Agosto del 2007, el cual contiene una serie de comandos y/o funciones agrupados según la tabla de contenidos mostrado en la f igura 1.1. >>help signal

Figura 1.1 Funciones agrupadas pertenecientes al Toolbox Signal Processing del Matlab. Para conocer los comandos y/o funciones pertenecientes a cada grupo, basta con dar un click sobre uno de ellos. Por ejemplo, al dar click sobre la palabra WINDOWS, se obtiene una lista de comandos que es mostrada en la figura 1.2. Dicha lista de comandos corresponde a los diferentes tipos de ventanas que son utilizadas en el diseño de filtros digitales. MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo – FIEM – UTP

Página 4


2011-III

Figura 1.2 Lista de comandos del contenido de WINDOWS del Toolbox Signal Processing del Matlab. 2.2

Representación de señales digitales en función al tiempo Las señales digitales, en su mayoría, son representadas en función al tiempo. Por lo tanto, su representación en el entorno del Matlab se realiza a través de vectores filas o vectores columnas. En un caso particular, cuando la señal digital haya sido digitalizada con dos o más números de canales, se optará por utilizar matrices. Donde, cada fila o columna, representará a un canal en particular de dicha señal a analizar. Por ejemplo, a continuación se muestra un vector conformado por números aleatorios entre -16 y 16, que a su vez corresponde a una señal de audio ruidoso con 12000 muestras. >> t = linspace( 0 , 1 , 12000 ); >> V1 = 32 * rand( 1 , 12000 ) - 16; >> V1 = round( V1 ); >> plot( t, V1 ) >> axis( [ 0.38 0.40 -20 20 ] ) >> grid >> sound( V1 , 8000 )

MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo – FIEM – UTP

%generando el vector de números %aleatorios. %redondeando los elementos del %vector V1. %graficando el vector de ruido %en la figura 1.3 %limitando la visualización del %vector de ruido de la figura 1.3 %reproduciendo el audio ruidoso

Página 5


2011-III

20

15

10

5

0

-5

-10

-15

-20 0.38

0.382 0.384 0.386 0.388

0.39

0.392 0.394 0.396 0.398

0.4

Figura 1.3 Parte del gráfico de la señal de audio ruidoso. Asimismo, también es posible representar las señales de audio, utilizando vectores. A continuación se muestra un ejemplo de la forma como se realiza la lectura de un archivo de audio de extensión WAV, a través del software Matlab. %el alumno debe de copiar %cualquier archivo de extensión %WAV a la carpeta de trabajo del %Matlab. >> figure( 3 ), plot( V2 ) %Ver figura 1.4 >> grid >> V2( 10 ) %verificando el valor de la %décima muestra de dicha señal. >> V2 = wavread( ‘ ejemplo.wav ‘ );

0.3

0.2

0.1

0

-0.1

-0.2

-0.3

-0.4 0

0.5

1

1.5

2

2.5

3

3.5

4

4.5 x 10

4

Figura 1.4 Señal de audio leído desde el Matlab


Página 6


2.3

2011-III

Representación de imágenes digitales Las imágenes digitales son representadas numéricamente a partir de matrices o arreglos, donde cada elemento o componente de la matriz o del arreglo, representa a un pixel con una tonalidad de gris o color adecuado. Por ejemplo, a continuación, se representa una matriz conformada por números aleatorios, que a su vez corresponde a una imagen de ruido con 128 filas y 128 columnas, y donde sus elementos se encuentran en el intervalo de 0 a 255. Asimismo, todo elemento o pixel de la matriz con valor igual a 0 representa al color negro, mientras que los pixeles con valores iguales a 255, representan al color blanco. Por otro lado, los pixeles con valores intermedios, corresponden a una intensidad de gris iniciando en el color negro y finalizando en el color blanco. >> IM = 255 * rand( 128 ,128 ); >> IM = round( IM ); >> figure( 1 ) >> colormap( gray ( 256 ))

%genera una matriz con componentes %aleatorios entre 0 y 255. %configura la imagen a tonalidades de %gris %Ver figura 1.5 %verificando el contenido de la cuarta

>> image( IM ) >> IM( 4 , 3) %fila y tercera columna

20

40

60

80

100

120 20

40

60

80

100

120

Figura 1.5 Imagen de tonos de gris generada de números aleatorios en el intervalo de 0 a 255. El Matlab también cuenta con un toolbox exclusivamente para el procesamiento de imágenes, llamado imagen processing. En este toolbox es posible encontrar variedades de funciones y/o comandos que son utilizados en el procesamiento espacial o frecuencial de una imagen digital. Asimismo, es posible cargar variables correspondientes a diversas imágenes con tamaños de 256x256 y 128x128 pixeles. MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo – FIEM – UTP

Página 7


>> load imdemos

2011-III

%cargar distintas variables propias del %Matlab que corresponden a imágenes %digitales en formato de grises.

>> figure( 2 )%ver figura 1.6 >> colormap( gray ( 256 ) ) >> subplot( 2 , 1 , 1) , image( trees) >> subplot( 2 , 1 , 2) , image( saturn )

20

20

40

40

60

60

80

80

100

100

120

120 20

40

60

80

100 120

20

40

60

80

100 120

Figura 1.6 Ejemplo de imágenes en tonos de gris. a) Imagen TREES. b) Imagen SATURN Asimismo, también es posible representar las imágenes en color, utilizando tres matrices. A continuación se muestra un ejemplo de la forma como se realiza la lectura de un archivo de imágenes de extensión JPEG, a través del software Matlab. >> IM2 = imread( ‘ ejemplo.jpg ‘ );

>> figure( 2 )%ver figura 1.7 >> image( IM2 ) >> IM2( 5 , 15 , : ) >> size( IM2 )

%el alumno debe de copiar cualquier %archivo de extensión JPEG o BMP, %a la carpeta de trabajo del Matlab. %verificando el color contenido en el %pixel de la quinta fila y quinceava %columna. %Esto muestra la presencia de tres %matrices de 600x800 pixeles.

ans = 600 800

3


Página 8


2011-III

100

200

300

400

500

600

100

200

300

400

500

600

700

800

Figura 1.7 Archivo de Imagen leído desde el Matlab Finalmente, también es posible representar un video digital en el entorno del Matlab, pero en este caso utilizando arreglos multidimensionales tal como se muestra en el siguiente ejemplo. >> load mri >> whos

%cargando un arreglo

Name Attributes

Size

Bytes

Class

D Map siz

4-D 89x3 1x3

442368 2136 24

uint8 double double

>> size( D ) ans = 128 128

%27 cuadros de 128x128 1

27

>> colormap(gray( 256 ) ) %ver figura 1.8 >> subplot( 3 , 2 , 1 ), image( D( : , : , 1 , 1 ) ) >> subplot( 3 , 2 , 2 ), image( D( : , : , 1 , 2 ) ) >> subplot( 3 , 2 , 3 ), image( D( : , : , 1 , 3 ) ) >> subplot( 3 , 2 , 4 ), image( D( : , : , 1 , 4 ) ) >> subplot( 3 , 2 , 5 ), image( D( : , : , 1 , 5 ) ) >> subplot( 3 , 2 , 6 ), image( D( : , : , 1 , 6 ) )


Página 9


20 40 60 80 100 120

2011-III

20 40 60 80 100 120 20

40

60

80 100 120

20 40 60 80 100 120

20

40

60

80

100 120

20

40

60

80

100 120

20

40

60

80

100 120

20 40 60 80 100 120 20

40

60

80 100 120

20 40 60 80 100 120

20 40 60 80 100 120 20

40

60

80 100 120

Figura 1.8 Representación de las 06 primeras imágenes del arreglo multidimensional.

3.- EJERCICIOS POR SOLUCIONAR: 3.1 Leer un archivo de audio y presentarlo en un vector. Posteriormente, generar un vector aleatorio de la misma dimensión para ser sumado al vector de audio original. Reproducir y captar la diferencia entre ambos vectores. 3.2 Leer un archivo de imagen y presentarlo en una matriz. Posteriormente, generar una matriz aleatoria de la misma dimensión para ser sumado a la matriz de imagen original. Visualizar y captar la diferencia entre ambas matrices. 3.3 Generar aleatoriamente, un arreglo multidimensional conformado por 40 matrices de 64x64 pixeles, y a tonos de gris. Finalmente, presentar en una sola ventana las últimas 8 matrices.


Página 10


2011-III


GRÁFICO DE SEÑALES PERIÓDICAS Y NO PERIÓDICAS. CONVOLUCIÓN.

1. OBJETIVOS: 1.1 Graficar señales periódicas discretas en el tiempo. 1.2 Graficar señales no periódicas discretas en el tiempo. 1.3 Realizar la operación de convolución entre señales discretas. 2. FUNDAMENTO TEÓRICO: 2.1

Funciones discretas Entre las más importantes encontramos: 0, n  k  [n  k ]   El impulso unitario: 1, n  k

0, n  k 1, n  k

 [n  k ]  

El escalón unitario:

2.2

Señales periódicas Las señales periódicas son aquellas señales que muestran periodicidad respecto del tiempo, esto quiere decir que describen ciclos repetitivos. Ver la figura 1.1. Por lo tanto, se cumple la siguiente expresión matemática: x (t) = x (t + T) = x (t + nT),

con “n” como número entero.

4 3 2 1 0 -1 -2 -3 -4

0

0.1

0.2

Figura 1.1


0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Ejemplo de señal periódica.

Página 11


2.3

2011-III

Señales No Periódicas. Las señales no periódicas son aquellas señales en donde no se muestra periodicidad respecto del tiempo, esto quiere decir que no describen ciclos repetitivos. Por lo tanto, no se cumple la expresión matemática anteriormente planteada. A continuación se muestra un ejemplo en la figura 1.2. 0.4 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4

0

0.5

1

1.5

2

2.5 x 10

Figura 1.2 2.4

4

Ejemplo de señal no periódica.

Operación de convolución. La operación de convoluciónes conmutativa y se realiza sobre dos señales discretas y finitas En caso de una señal de entrada x[n], el resultado de la convolución está dada por la siguiente expresión matemática.

x [n]

y [n]

h [n]

y[n]  x[n]  h[n]  h[n]  x[n] 

y[n] 

 x[k ] h[n  k ] k 

3. EJERCICIOS SOLUCIONADOS: 3.1

Graficando señales “contínuas” y discretas. Impulso y Escalón Unitario.

Para representar una señal en forma continua se hace uso del comando o función “plot”. Esta se encarga de unir los puntos dando una apariencia de continuidad. Por otro lado, para graficar una señal discreta, se utiliza el comando o función “stem”, que se encarga de graficar mediante impulsos, la señal a representar. Continuación un ejemplo, ver la figura 1.3. >> n = [ 0 1 2 3 4 5 6 ]; MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo – FIEM – UTP

%tiempo discreto Página 12


2011-III

>> x = [ 5 9 3 -4 0 8 7 ]; %señal discreta >> figure(1) >> subplot(1 , 2 , 1 ) , plot ( n , x) %gráfico “continuo” >> title ( ‘Señal Contínua’ ), xlabel (‘tiempo’ ) >> subplot(1 , 2, 2 ) , stem ( n , x) %gráfico discreto >> title ( ‘Señal Discreta ), xlabel (‘tiempo’ ) Señal Contínua

Señal Discreta

10

10

8

8

6

6

4

4

2

2

0

0

-2

-2

-4

0

2

4

6

-4

0

2

tiempo

Figura 1.3

4

6

tiempo

Ejemplo de señal contínua y discreta.

Una señal impulso: 2 [n - 1] >> n1 = [ 0 : 7]; >> x1 = 2 * [ 0 1 0 0 0 0 0 0]; Una señal escalón: -5 [ n] >> n2 = [ -20 : 1 : 20]; >> x2 = -5 * [ zeros( 1, 20 ) ones(1,21) ]; Una señal de ruido entre 0 y 1: r[n] >> r3 = rand ( 1, 1000 ); >> n3 = 0 :1: 999; >> subplot(3,1,1), stem ( n1 , x1 ) >> subplot(3,1,2), stem( n2 , x2 ) >> subplot(3,1,3), stem( n3 , r3 )


%ver figura1.4

Página 13


2011-III

2

1

0 0

1

2

3

4

5

6

7

0 -2 -4 -6 -20

-15

-10

-5

0

5

10

15

20

1

0.5

0 0

Figura 1.4 3.2

100

200

300

400

500

600

700

800

900

1000

Ejemplo de señal impulso, escalón y ruido.

Graficando una señal sinusoidal. Para representar una señal seno o coseno en el Matlab, debe de indicarse la variable temporal señalando el tiempo de duración de la onda. Por ejemplo, para graficar una señal seno de frecuencia igual a 3 Hz, amplitud igual a 2 voltios y fase igual a 90 º, se aplica el siguiente procedimiento. Discretizando una señal Senoidal continua, para luego graficarla en el dominio del tiempo discreto. x(t) = A * sin (2*pi*f *t + fase ) Para discretizar, reemplazamos t por nT en la expresión anterior. x[nT] = sin (2*pi*f*nT + fase) Donde: T es el periodo de muestreo o 1/Fs x[ n] = sin(2*pi*f*n / Fs + fase) >> Fs = 100; >> n = 0:Fs-1; >> fase = 90; >> A = 2; >> F = 3; >> x = A * sin ( 2*pi* F*n / Fs + fase*pi/180); >> stem ( n , x ,’r’ )


%frecuencia de muestreo

%frecuencia fundamental: % Fs> 2*F % Ver figura 1.5

Página 14


2011-III

2

1.5

1

0.5

0

-0.5

-1

-1.5

-2 0

10

20

Figura 1.5

30

40

50

60

70

80

90

100

Ejemplo de señal sinusoidal.

Sumando un ruido a la señal seno: >> r = length(x); >> R = randn(1,r); >> xR = 2*x + R; >> subplot( 1 , 2 , 1) , plot( n , xR) >> subplot( 1 , 2 , 2) , stem( n , xR) Cuando no se cumple con el teorema de muestreo, se tiene una representación equivocada de la señal discreta. Por ejemplo, a continuación se grafica una onda seno con frecuencia fundamental igual a 20 Hz y frecuencia de muestreo igual a 30 Hz. En este caso no se cumple la relación de tener una Fs> 2 * Fo. >> Fs = 30; >> F = 20; >> n = 0:Fs-1; >> Fase = 90; >> x = sin ( 2*pi* F*n / Fs + fase*pi/180); >> stem ( n , x ,’b’ ) >> hold on >> plot( n , x , ‘r ’ ) >> hold off


% utilizado para congelar la % figura y %volver a graficar % sobre ella. Ver figura 1.6

Página 15


2011-III

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

0

5

10

Figura 1.6

15

20

25

30

Señal sinusoidal distorsionada.

También es posible graficar otros tipos de señales periódicas tal como es el caso de la onda cuadrada. Por lo tanto, para graficar la onda cuadrada, un tren de pulsos, o una modulación por ancho de pulso, se utiliza el comando o función del Matlab denominado: “SQUARE”.

>> help square >> Fs = 1000; >> t = linspace( 0 ,1 , Fs ); >> x = 1 + square( 2 * pi * t * 4 , 20 ); >> plot( t , x )

%ver figura 1.7

2 1.8 1.6 1.4 1.2

1 0.8 0.6 0.4 0.2

0 0

Figura 1.7 3.3

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Señal tren de pulsos obtenido de la onda cuadrada.

Graficando una señal no periódica. Las últimas versiones del Matlab cuentan con comandos o funciones que permiten representar y posteriormente graficar señales de no periódicas, tal es el caso de la señal de electrocardiograma. Para ello, se utiliza la función o comando “ECG”, que permitirá graficar un latido cardiaco mostrando las ondas P, Q, R, S y T. Considerando que el latido corresponde a una persona sana, entonces obtenemos el tiempo de duración para un latido cardiaco: 70 lat / min. >> help ecg >> x = ecg(1000);


% considerando 1000 muestras Página 16


>> Lat = 60 / 70; >> t = linspace( 0 , Lat , 1000 ); >> plot( t , x ), grid >> text(0.35,0.7, 'Complejo QRS')

2011-III

% tiempo de duración de un latido % 70 latidos por minuto % ver figura 1.8

1 0.8 Complejo QRS 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1

0

0.1

0.2

Figura 1.8

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Señal sinusoidal distorsionada.

Asimismo, para graficar la función SINC, se utiliza el comando o función “SINC”. Tal como es mostrado en la figura 1.8. f (t )  sin c(  ) 

sen(  )

 

>> theta=linspace(-10,10,100); >> w =sinc( theta ); >> subplot( 1 , 2 , 1) , >> plot(t,W) >> subplot( 1 , 2 , 2) , >> stem(t,W)

3.4

%Forma Contínua. Figura 1.9 %Forma Discreta. Figura 1.9

Convolución. Para realizar la convolución entre dos señales finitas, o secuencias, habrá que definir cada una de ellas en un vector, y luego utilizar el comando “CONV”. Por

ejemplo, hacer la convolución entre x[n], h 1[n] y h2[n].


Página 17


1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

-0.2

-0.2

-0.4 -10

-5

0

5

Figura 1.9

-0.4 -10

10

-5

0

5

2011-III

10

Función SINC. En forma continua y discreta.

x[n]   [n]  4 [n  1]  5 [n  3]  6 [n  4]

h1[n]   [n]  2 [n  5] h2 [n]  0.5 [n  7]  8 [n  8] 17 [n  9]

x [n]

y [n]

h2 [n]

h1[n]

>> x = [ 1 4 0 5 6 ]; >> h1 = [ 1zeros(1,4) -2 ]; >> h2 = [ zeros(1,7) -0.5 8 17 ]; >> y1 = conv( x , h1 ); % o también conv( h1 , x ) >> y = conv( y1 , h2 ); % o también conv( h2 , y1 ) >> stem( 0:length(y) -1 , y ) % ver figura 1.10 150

100

50

0

-50

-100

-150

-200

-250

-300 0

2

4

Figura 1.10


6

8

10

12

14

16

18

Resultado de la convolución.

Página 18


2011-III

4. EJERCICIOS POR SOLUCIONAR: 4.1

Graficar 750 mili segundos de una señal triangular. Dicha señal deberá tener una frecuencia igual a 12 Hz, una amplitud igual a 1.5 voltios y un nivel DC igual a -0.75 voltios.

4.2

Graficar 10 latidos cardiacos, uno a continuación del otro, de tal forma que el primer, quinto y octavo latido tengan una duración de 0.9 segundos, mientras que los latidos restantes tengan una duración de 0.7 segundos.

4.3

Del diagrama de bloques mostrado (figura 1.11), obtener la señal de salida. Para ello, se plantea dos filtros digitales representados en el tiempo discreto h1[n] y h2[n], así como también una señal de ruido representada por r[n]. Considerar:

x [ n ] =  [n] - 2  [n-2] h1 [ n ] =  [n] - 8  [n-1] + 3  [n-2 ] h2 [ n ] = 2  [n - 1] + 2  [n - 3] r[n]= señal de ruído

r [ n]

3  [n]

+

h1[ n ]

+

X

y [ n]

x [ n]

h2[ n ]

Figura 1.11 Diagrama de bloques por analizar.


Página 19


2011-III


CAMBIO DE LA TASA DE MUESTREO: DECIMACIÓN E INTERPOLACIÓN

1. OBJETIVOS: 1.1 Graficar señales decimadas en el tiempo. 1.2 Graficar señales interpoladas en el tiempo. 1.3 Analizar y graficar cambios de la tasa de muestreo por un factor no entero. 2. FUNDAMENTO TEÓRICO: 2.1

Decimación Es una operación encargada de disminuir la frecuencia de muestreo por un factor entero denominado M.

Pasa-bajo Fcorte = / M

2.2

Interpolación Es una operación encargada de aumentar la frecuencia de muestreo por un factor entero denominado L.

↑L

2.3

M

Pasa-bajo Fcorte = / L

Cambio de la tasa de muestreo por un factor no entero Cuando se utiliza la operación de decimación e interpolación a la vez, con la finalidad de cambiar la tasa de muestreo por un factor no entero: L / M

↑L


Pasa-bajo Fcorte = / L

Pasa-bajo Fcorte = / M

M

Página 20


2011-III


Seguidamente se realiza la decimación e interpolación de una señal periódica en ausencia del filtro pasabajo mostrado en los diagramas de bloques anteriores. >> Fs = 100; >> t = linspace( 0 , 1 , Fs ); >> x = sin( 2 * pi * 1 * t) + 2 * cos( 2 * pi * 3 * t ); >> subplot( 1 , 3 , 1 ) , stem ( t , x ), grid, title(‘ Señal Original ’) >> M = 4; >> L = 3; >> xd = downsample( x , M ); >> Fs1 = Fs / M ; >> t1 = linspace( 0 , 1 , Fs1 ); >> subplot( 1 , 3 , 2 ) , stem( t1 , xd ), grid, title(‘ Señal Decimada ’) >> xi = upsample( x , L ); >> Fs2 = Fs * L ; >> t2 = linspace( 0 , 1 , Fs2 ); >> subplot( 1 , 3 , 3 ) , stem( t2 , xi ) , gr id, title(‘ Señal Interpolada ’) Señal Original

Señal Interpolada

3

3

2

2

2

1

1

1

0

0

0

-1

-1

-1

-2

-2

-2

-3

0

0.5

Figura 3.1 3.2

Señal Decimada

3

1

-3

0

0.5

1

-3

0

0.5

1

Resultado de la decimación e interpolación.

Luego, se muestra el desarrollo de la operación de decimación completa utilizando el filtro pasa-bajo señalado en el diagrama de bloques. Sea una secuencia x[n]=2sin (2* *f*n/Fs). Representemos la versión decimada por 2, considerando una Fs = 100 muestras/seg. y una frecuencia fundamental igual a 5 Hz. >> Fs = 100; >> n = 0 : Fs-1; >> f = 5; >> x = 2 * sin(2*pi*f*n/Fs );


Página 21


>> stem( n , x ) >> help decimate >> x2 = decimate(x,2); >> Fs = Fs / 2; >> stem(0:Fs-1 , x2 )

2011-III

%ver figura 3.2

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

0

Figura 3.2

5

10

15

20

25

30

35

40

45

50

Versión decimada de la onda seno

Se puede observar una característica importante. En la figura 3.1 se percibe la presencia de 100 muestras representando 5 ciclos por segundo. En cuanto que en la figura 3.2, se aprecia la presencia de solo 50 muestras también mostrando 5 ciclos por segundo. 3.3

Luego, se muestra el desarrollo de la operación de interpolación completa utilizando el filtro pasabajo señalado en el diagrama de bloques. Seguidamente interpolamos por 2 para retornar al número de muestras inicial de esta señal. >> Fs = Fs*2; >> xx = interp(x2,2); >> stem(0:Fs-1 , xx )

%ver figura 3.3 2

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

Figura 3.3


0

10

20

30

40

50

60

70

80

90

100

Versión interpolada de la onda seno

Página 22


2011-III

De esta manera, se recupera la Fs de muestreo inicial. A continuación, es posible observar la diferencia entre la señal original y la manipulada por una operación de decimación e interpolación (ver figura 3.4). >> plot(0:Fs-1,x,'r',0:Fs-1,xx,'b') 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5

0

Figura 3.4 3.4

10

20

30

40

50

60

70

80

90

100

Diferencia de señales decimada e interpolada

Asimismo, para lograr el cambio de la Frecuencia de Muestreo, Fs, por un número fraccionario de veces, se procede a realizar ambas operaciones a la vez, tal como lo muestra la siguiente figura.

x [ n ]

 L

Fs * L

Fs * L

Fs

xT [ n ]

 M M

Por ejemplo, si se desea una Fs_Final = 300 muestras/seg, a partir de una Fs_Inicial = 400 muestras/seg, se deberá de realizar las operaciones de decimación e interpolación una seguida de la otra.

Fs_Incial

3

4

Fs_Final

>> Fs = 400; >> t = linspace( 0 , 1 , Fs ); >> x = cos( 2 * pi * t * 0.5 ) + cos( 2 * pi * t * 1.5 ) >> subplot(1,2,1), stem( t , x ) >>help resample MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo – FIEM – UTP

Página 23


>> L = 3; >> M = 4; >> xr =resample(x,L,M); >> Fs_n = Fs * L / M; >> t1 = linspace( 0 , 1 , Fs_n ); >> subplot(1,2,2), stem( t1 , xr )

2011-III

%ver figura 3.5

2

2.5 2

1.5

1.5 1 1 0.5

0.5

0

0 -0.5

-0.5

-1 -1 -1.5 -1.5

-2

-2

-2.5 0

Figura 3.5

0.5

1

0

0.5

1

Señal re-muestreada por un factor no entero.


Con ayuda del comando o la función “ECG”, cargar una señal de electrocardiograma con 1000 muestras por segundo; luego, re-muestrear dicha señal, de tal forma que tres latidos cardiacos continuos se encuentren muestreados a 1200 muestras por segundo.

4.2

Implementar una señal de tono único de frecuencia 3 KHz, amplitud igual a 5 voltios y frecuencia de muestreo de 40 KHz. Cambiar la tasa de muestreo a 20 KHz y posteriormente a 5 KHz. ¿Qué cambios se nota al realizar estas operaciones?. Reproducir dichas se ñales con ayuda del comando “SOUND”, y mostrar sus comentarios y observaciones.

4.3

Implementar una señal multitono conformado por algunas notas musicales, tal como se indica a continuación. Utilizar una frecuencia de muestreo de 10 KHz. SEÑAL =[DO, DO, DO, FA, LA, DO, DO, DO, FA, LA, FA, FA, MI, RE ]; Cambiar la tasa de muestreo de tal forma que el nuevo periodo de muestreo sea de 125 microsegundos. Además, las notas FA y LA deberán tener el doble de duración en el tiempo respecto a las otras notas musicales.


Página 24


2011-III


TRANSFORMADA DISCRETA DE FOURIER. DIRECTA E INVERSA

1. OBJETIVOS: 1.1 Obtener la transformada discreta de Fourier directa e inversa. 1.2 Representar el espectro de frecuencia de una señal periódica. 1.3 Representar el espectro de frecuencia de una señal no periódica. 2. FUNDAMENTO TEÓRICO: 2.1

Transformada Discreta de Fourier (TDF) La TDF es la herramienta principal del Procesamiento Digital de Señales. El Toolbox de SignalProcessingcuenta con un comando o función que nos ayuda a calcular la Transformada Discreta de Fourier “FFT”. Seguidamente se muestra la expresión para el cálculo de la Transformada discreta de Fourier directa e inversa.

X (k ) 

N 1

 x[n] e

-2 j n N

k

n0

x[n] 

1

N 1

X (k ) e  N

2 j n N

k

k 0

k  0, 1, 2, ..., N 1 2.2

Representación del algoritmo de la TDF para el entorno del Matlab A continuación se muestra el algoritmo de la TDF para ser ejecutado en el entorno del Matlab. De la misma manera, es posible adaptar este algoritmo a la sintaxis de cualquier otro software de programación, de tal forma que pueda ser ejecutado sin problema alguno. N = 1024; %definir un valor para N x = 10 * rand( 1 , N) %definir una señal con muestras discretas for k = 0:N-1 a = x(1) * exp( -2 * pi * j * k * 0 / N); for n = 1 : N-1 a = x( n + 1 ) * exp( -2 * pi * j * n * k / N ) + a;


Página 25


2011-III

end TX( k + 1 ) = a; end Comparando este resultado con el del comando o función “FFT”.

>> TX1 = fft( x , N ); >>[ TX’ TX1‘ ]

2.3

Transformada Rápida de Fourier (FFT) La transformada rápida de Fourier es un algoritmo para agilizar el cálculo de la transformada discreta de Fourier, en el cual disminuye el número de operaciones de sumas y multiplicaciones entre números complejos.


Transformada Discreta de Fourier de una señal periódica. Sea una señal coseno x[n], con una frecuencia de muestreo Fs=100 Hz y una frecuencia fundamental de 20 Hz. x[n] = cos(2 *20*n/Fs). A continuación se grafica la señal x[n] en el tiempo para un segundo de duración, Seguidamente, se obtiene la Transformada Discreta de Fourier (DFT o FFT) utilizando una cantidad de muestras N = 16 y N = 512. Finalmente, se grafica el módulo y fase, y la parte Real e Imaginaria de su espectro. >> Fs = 100; >> n=0:Fs-1; >> x = cos(2*pi*n*20/Fs); >> stem(n,x)

%ver figura 4.1

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

10

Figura 4.1


20

30

40

50

60

70

80

90

100

Señal discreta con Fs = 100 Hz.

Página 26


2011-III

>> X_16 = fft(x,16); >> X_512 = fft(x,512); Luego, se obtiene el módulo y fase de la TDF y se grafica respecto a su eje correspondiente. >> mX_16 = abs (X_16); >> fX_16 = angle (X_16);

% módulo para N=16 puntos % fase para N=16 puntos

>> mX_512 = abs (X_512); >> fX_512 = angle (X_512); >> figure(1) >> f_16 = linspace(0,Fs,16); >> f_512 = linspace(0,Fs,512);

% módulo para N=512 puntos % fase para N=512 puntos

>>subplot(2,1,1), stem(f_16 , mX_16) >>subplot(2,1,2), stem(f_16 , fX_16)

%ver figura 4.2

8 6 4 2 0

0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

4 2 0 -2 -4

Figura 4.2

Transformada Discreta de Fourier. Módulo y Fase (N=16).

>> figure(2) >> subplot(2,1,1), stem(f_512 , mX_512) >> subplot(2,1,2), stem(f_512 , fX_512)


%ver figura 4.3

Página 27


2011-III

50 40 30 20 10 0

0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

4 2 0 -2 -4

Figura 4.3

Transformada Discreta de Fourier. Módulo y Fase (N=512).

De igual manera se puede graficar la parte real e imaginaria del espectro de la señal. >> rX_16 = real (X_16); >> iX_16 = imag (X_16); >> rX_512 = real (X_512); >> iX_512 = imag (X_512); >> subplot(2,1,1), stem(f_512 , rX_512) >> subplot(2,1,2), stem(f_512 , iX_512)

%ver figura 4.4

60 40 20 0 -20

0

10

20

30

40

50

60

70

80

90

100

0

10

20

30

40

50

60

70

80

90

100

40 20 0 -20 -40

Figura 4.4

Transformada Discreta de Fourier. Real e Imaginaria (N=512).


Página 28


2011-III

Obteniendo y graficando la Transformada Discreta de Fourier Inversa IFFT. >> ix = ifft(X_512 , 512); >> rix = real(ix); >> stem( 0:Fs-1 , rix(1:Fs) )

%ver figura 4.5

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

0

Figura 4.5

10

20

30

40

50

60

70

80

90

100

Transformada Discreta de Fourier Inversa

Graficando la FFT con N=1024 para la suma de tres señales cosenos de frecuencias: 10 Hz, 30 Hz, 43 Hz (Fs=120). La presentación se muestra desde –Fs/2 a Fs/2. >> Fs = 120; >> n = 0:Fs-1; >> x = cos(2*pi*n*10/Fs) + 2*cos(2*pi*n*30/Fs) + 4*cos(2*pi*n*43/Fs); >> tX = fft(x,1024); >> tX = fftshift( tX ); >> mtX = abs(tX); >> f = linspace( -Fs/2 ,Fs/2 ,1024); >> plot( f, mtX), grid %ver figura 4.6 250

200

150

100

50

0 -60

Figura 4.6 MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo – FIEM – UTP

-40

-20

0

20

40

60

Módulo de la TDF de una señal periódica

Página 29


2011-III

Graficando la FFT de una Señal Coseno sumado a un ruido >> Fs=100; >> n = 0:Fs-1; >> x = 2*cos(2*pi*n*10/Fs); >> r = randn(1,Fs) / 1.5 ; >> xr = x + r; >> figure(1) , plot(xr) >> tXR = fft ( xr , 1024 ); >> mtXR = abs (tXR); >> f = linspace(0,Fs,1024); >> plot(f,mtXR) Graficando la FFT con N = 1024 para una Señal Cuadrada de F=10 Hz. >> Fs=200; >> n = 0:Fs-1; >> x = square(2*pi*n*10/Fs); >> figure(1) , stem(x)

%ver figura 4.7

100 90 80 70 60 50 40 30 20 10 0

Figura 4.7

0

10

20

30

40

50

60

70

80

90

100

Módulo de la TDF de una señal periódica con ruido

>> tX = fft ( x , 1024 ); >> mtX = abs (tX); >> mtX = fftshift( mtX ); >> f = linspace( -Fs/2 ,Fs/2 ,1024); >> plot(f,mtX)


%ver figura 4.8

Página 30


2011-III

140

120

100

80

60

40

20

0 -100

Figura 4.8 3.2

-80

-60

-40

-20

0

20

40

60

80

100

Módulo de la TDF de una señal periódica (cuadrada)

Transformada Discreta de Fourier de una señal no periódica. Para mostrar el espectro de frecuencia en módulo de una señal no periódica, se procede a abrir un archivo de audio de la siguiente forma. >> [ x , Fs , nB ] = wavread(‘chimes.wav’); >> % nB : Representa el número de bits utilizado para codificar cada muestra >> % Fs : Representa la frecuencia de muestreo. >> tam = length( x ); >> x1 = x( : , 1 ); >> tX = fft ( x1 , tam ); >> mtX = abs (tX); >> mtX = fftshift( mtX ); >> f = linspace( -Fs/2 , Fs/2 , tam ); >> plot( f , mtX )

%tomando solamente un canal de “x”

%ver figura 4.9

150

100

50

0 -1.5

-1

-0.5

0

0.5

1

1.5 4

x 10

Figura 4.9

Módulo de la TDF de una señal no periódica (audio)


Página 31


2011-III


Seguidamente se muestra una señal x[n] digitalizada con un periodo de muestreo de 125 mSeg y 8 bits por muestra. Si esta señal se encuentra compuesta por la suma de senos y cosenos, se solicita graficar el módulo de su transformada discreta de Fourier. Considere N=256. 3



x[n]  2

sen(

m1

 n 2

4

m) 

 m0

cos(

 (2n  1) 2

m

 2

)

4.2

Grafique el módulo de la transformada discreta de Fourier, del producto de dos ondas periódicas. La primera onda corresponde a una señal cuadrada de frecuencia igual a 250 Hz, amplitud igual a 2.5 voltios y nivel de continua de 2.5 voltios. Y, la segunda onda corresponde a una onda coseno de frecuencia igual 1 KHz y amplitud igual a 2 voltios. Su respuesta deberá ser presentada en el intervalo de -Fs/2 a Fs/2.

4.3

Grafique el módulo de la transformada discreta de Fourier, del producto de una onda periódica y otra no periódica. La onda periódica debe corresponder a una onda coseno de 10 KHz. Elija usted la frecuencia de muestreo y el valor de N adecuado.


Página 32


2011-III


TRANSFORMADA DISCRETA DE FOURIER CORTA EN EL TIEMPO. ESPECTROGRAMA

1. OBJETIVOS: 1.1 Obtener la transformada discreta de Fourier corta en el tiempo (STFT). 1.2 Representar el espectrograma de una señal periódica. 1.3 Representar el espectrograma de una señal no periódica.

2. FUNDAMENTO TEÓRICO: 2.1

Transformada discreta de Fourier corta en el tiempo La transformada discreta de Fourier corta en el tiempo, es una representación de tiempo versus frecuencia, donde se visualiza las variaciones de la frecuencia y de la intensidad, a través de colores, delaseñal que se está representando a lo largo del tiempo. A continuación, la expresión matemática correspondiente. 

STFT  x[n]   X (m,  ) 

 x[n] w[n  m]e

 j  n

n

Donde: x[n] : X(m,) : w[n] :

2.2

señal discreta en el tiempo Transformada de Fourier corta en el tiempo de x[n] ventana para segmentar la señal x[n]

Obtención de la STFT La STFT se obtiene calculando la Transformada Discreta de Fourier a segmentos de la señal por analizar. El tamaño del segmento será delimitado por la ventana w[n], quién definirá diferentes niveles de resolución de la STFT. Cuando se utiliza ventanas estrechas, mejora la resolución en el tiempo pero empeora la de frecuencia. Y, por lo contrario, si se utiliza ventanas anchas, mejora la resolución de frecuencia empeorando la de tiempo. A continuación, en la figura 5.1, se muestra un ejemplo de representación de diferentes niveles de resolución de la STFT, donde se aprecia un aumento de la resolución a nivel de tiempo, pero disminución a nivel de la frecuencia, y viceversa. El espectrograma se obtiene a partir del cuadrado del módulo de la STFT.


Página 33


2011-III

Frecuencia

Frecuencia

Tiempo

Figura 5.1

Tiempo

Comparación de resoluciones resoluciones de la STFT

El software Matlab posee un comando de nombre “SPECTROGRAM”, que permite obtener el espectrograma de una señal a partir de ciertos parámetros a considerar. >> help spectrogram Por ejemplo, a continuación se representa el espectrograma en 2D y 3D de una señal coseno con 1200 mili segundos de duración. Para dicha representación se utiliza primero, una ventana hamming ancha de 200 muestras; luego,una ventana hamming estrecha de 50 muestras.Esto puede ser observado en las figuras 5.2 y 5.3, respectivamente. Para la representación en 3D, recurrir a la opción ROTATE 3D de la barra de herramientas de la ventana figura. >> Fs = 1000; >> t = linspace( 0 , 1.5 , 1.5*Fs ); >> x = cos( 2 * pi * t * 180 ); >> NFFT = 200; % número de muestras de la TDF. >> WINDOW = hamming(NFFT); %tipo de ventana a utilizar. >> NOVERLAP = 180; % número de muestras sobrepuestas. sobrepuestas. >> spectrogram( x ,WINDOW,NOVERLAP,NFFT,Fs); >> Fs = 1000; >> t = linspace( 0 , 1.5 , 1.5*Fs ); >> x = cos( 2 * pi * t * 180 ); >> NFFT = 100; % número de muestras de la TDF. >> WINDOW = hamming(NFFT); %tipo de ventana a utilizar. >> NOVERLAP = 80; % número de muestras sobrepuestas. sobrepuestas. >> spectrogram( x ,WINDOW,NOVERLAP,NFFT,Fs);


Página 34


(a)

Figura 5.2

(b)

Comparación de espectrogramas espectrogr amas en 2D. a) Utilizando una ventana ancha. b) Utilizando una ventana estrecha.

(a)

Figura 5.3

2011-III

(b)

Comparación de espectrogramas en 3D. a) Utilizando una ventana ancha. b). Utilizando una ventana estrecha.


Espectrograma Espectrograma de una señal periódica Considerando como señal periódica, una onda cuadrada de 1 segundo de duración, frecuencia igual a 100 Hz, amplitud igual 2 voltios y nivel DC igual 1.5 voltios. Asimismo, suponiendo un periodo de muestreo de 0.5 mili segundos. >> Ts = 0.5/1000; >> Fs = 1 / Ts; >> t = linspace( 0 , 1 , Fs ); >> x = 1.5 1.5 +square( 2* pi * t * 100 ); >> NFFT = 512; >> WINDOW = hamming(NFFT); hamming(NFFT);


Página 35


2011-III

>> NOVERLAP = 500; >> spectrogram( spectrogram( x ,WINDOW,NOVERLAP,NFFT,Fs); ,WINDOW,NOVERLAP,NFFT,Fs); %Ver figura 5.4

Figura 5.4 3.2

Representación Represent ación del espectrograma espectrogr ama de una señal periódica

Espectrograma Espectrograma de una señal no periódica Para representar el espectrograma de una señal no periódica, se procede a abrir un archivo de extensión WAV, tal como se realizó en la experiencia de laboratorio 01. >> [ V2 , Fs ] = wavread( ‘ ejemplo.wav ‘ );

>> size( V2 ) >> x = V2( : , 1 );

>> figure( 1 ), plot( V2 ) >> grid >> NFFT = 512; >> WINDOW = hamming(NFFT); hamming(NFFT); >> NOVERLAP = 500;

%el alumno debe de copiar %cualquier archivo de %extensiónWAV %extensiónWAV a la carpeta %de trabajo del Matlab. %en caso de ser una señal %tipo estéreo, se toma solo un %canal %Ver figura 5.5

>> figure(2) >> spectrogram( spectrogram( x ,WINDOW,NOVERLAP,NFFT,Fs); ,WINDOW,NOVERLAP,NFFT,Fs); %Ver figura 5.6


Página 36


Figura 5.5

Figura 5.6

2011-III

Representación temporal de una señal no periódica

Representación del espectrograma en 3D de una señal no periódica


Realizar una modulación por amplitud de doble banda lateral, y luego representar su respectivo espectrograma. Considere como señal modulante la suma de tres señales sinusoidales de amplitudes 1, 4 y 8 voltios, así como con frecuencias iguales a 100, 400 y 800 Hz, respectivamente. Asimismo, considere una señal de portadora igual a 10 KHz y amplitud igual a 1 voltio. Deberá de elegir el tamaño y tipo de ventana para el cálculo del espectrograma, así como también la cantidad de muestras sobrepuestas y la frecuencia de muestreo correspondiente. El diagrama de bloques mostrado a continuación ayuda a interpretar el proceso de la modulación por amplitud con doble banda lateral.


Página 37


Señal Modulante

+

X

2011-III

Señal AM

10 Señal Portadora

4.2

Representar el espectrograma de una señal no periódica cuya frecuencia disminuye linealmente, mientras transcurreel tiempo. Deberá de elegir el tipo y tamaño de la ventana a utilizar, así como también el número de muestras que se sobrepondrán al momento de calcular la STFT. Elija una señal con 1.2 segundos de duración y con una amplitud de 1.8 voltios.

4.3

Grabar y digitalizar su propia voz con las siguientes palabras: “casa” y “caza”. Luego, con ayuda del espectrograma, señale y comente si existe alguna diferencia entre ambas pronunciaciones. Deberá de elegir adecuadamente, el tipo y tamaño de ventana, la cantidad de muestras y la amplitud de la señal digital. De preferencia, realice ambas digitalizaciones con un único tiempo de duración y frecuencia de muestreo.


Página 38


2011-III


TRANSFORMADA Z. DIAGRAMA DE POLOS Y CEROS. ESTRUCTURA DE FILTROS DIGITALES.

1. OBJETIVOS: 1.1 Obtener laTransformada Z con ayuda del toolbox Symbolic del Matlab. 1.2 Representar el diagrama de polos y ceros de filtros digitales. 1.3 Representar estructuras de filtros digitales con la Transformada Z.


Transformada Z La Transformada Z es una generalización de la Transformada de Fourier y es relacionada como la contraparte de la Transformada de Laplace. La expresión matemática que define a dicha transformada es mostrada a continuación, donde se define como Transformada Z Bilateral. En el caso de considerarse una Transformada Z Unilateral, por ejemplo para el caso de secuencias o señales causales, la sumatoria se inicia en n=0. 

X ( z ) 

 x[n] z

n

m

Seguidamente se muestra algunas funciones discretas y su respectiva Transformada Z. Sea:

x[n]

 X (z )

 [n]

 1

 [n]



Entonces:


1 1  z 1

Página 39


 z  m

 [n  m]  [n  m]



a n  [ n]



a 2.2

nm

2011-III

z m 1  z 1 1 1  az 1 z m

 [n  m] 

1  az 1

Diagrama de Polos y Ceros Un sistema discreto lineal con coeficientes constantes, con entrada x[n] y salida y[n], es representado a través de la siguiente ecuación en diferencia: N

a

k

M

b x[n  k ]

y[n  k ] 

k 0

k

k 0

……………….. (6.1)

Entonces, calculando su respectiva Transformada Z, la expresión anterior se transforma a: N

a

k

M

b

z Y (z)  -k

k 0

k

-k

z X (z)

k 0

De donde se obtiene la función de transferencia, tal como se indica a continuación.

H ( z ) 

Y ( z ) X ( z )

1



2

 M

bo  b1 z  b2 z  ...  b M z 1

2

 N

ao  a1 z  a2 z  ...  a N z

Luego, se utiliza la representación de diagrama de polos y ceros, en el plano Z (ver figura 6.1), para representar los polos y ceros de la función de transferencia anteriormente mostrada. Si los polos permanecen en el interior del círculo unitario (sea a la izquierda o derecha del eje imaginario), se considera un caso estable. Por el contrario, si por lo menos uno de los polos se sitúa en el exterior del círculo unitario, el sistema es inestable.


Página 40


Figura 6.1 2.3

2011-III

Círculo unitario para la representación de polos y ceros

Estructura de Filtros Digitales Para representar la estructura o diagrama de bloques de un filtro digital, se hace uso de tres elementos principales: retardo o delay, sumatoria y amplificadores o ganancias. De la ecuación en diferencia general (ecuación 6.1), es posible distinguir tres casos fundamentales de modelos: Modelo MA (medias móviles): Modelo AR (auto-regresivo): Modelo ARMA :

cuando N = 0 y M ≠ 0 cuando N ≠ 0 y M = 0 cuando N ≠ 0 y M ≠ 0

Para el primer caso, modelo MA, se obtiene su correspondiente Transformada Z y se procede a graficar su estructura o diagrama de bloques utilizando tres elementos principales (ver figura 6.2). H ( z )

Y ( z )

b o

b 1 z

1

X ( z )

X (z)

...

b M z

M

...

Z-1

b1

Z-1

b2

∑

Figura 6.2

2

1

Z-1

b0

b 2 z

bM

∑

...

Y (z)

∑

Estructura de filtros digitales correspondiente al modelo MA.


Página 41


2011-III


Transformada Z en forma literal Con ayuda del toolboxSymbolic del Matlab, es posible obtener la Transformada Z directa o inversa en forma literal. >> help ztrans >> help iztrans >> syms n z >> f1 = 0.2^(n-2); >> F1 = ztrans( f1 ) >> pretty( F1 )

%se definen las variables literales: tiempo y Z. %se obtiene la Transformada Z directa

25 z ------z - 1/5 >> F2 = 10 / ( 1 – 0.6 / z - 0.07 / z / z ) >> f2 = iztrans( F2 ) %se obtiene la Transformada Z inversa >> pretty( f2 ) n n 5 (-1/10) 35 (7/10) ---------- + ---------4 4

3.2

Diagrama de polos y ceros. Dada una función de transferencia en términos de Z, es posible graficar su respectivo diagrama de polos y ceros con la finalidad de comprobar la estabilidad o inestabilidad en el sistema discreto. H ( z )

Y ( z ) X ( z )

2 z 1 0 . 7 z

1

1

0 . 1 z

2

>> Nu = [ 0 2 ]; >> De = [ 1 0.7 0.1 ]; >> zplane( Nu , De ) >> roots( De) ans =

%obteniendo las raíces del denominador

-0.5000 -0.2000


Página 42


Figura 6.3

2011-III

Diagrama de polos y ceros de la función de transferencia H(z) del ejemplo anterior.

>> [ r , p , k ] = residuez( Nu , De )

% Obtener la expansión en % fracciones parciales

r= -6.6667 6.6667 p= -0.5000 -0.2000 k= [] >> dstep( Nu , De , 20 )

Figura 6.4

% para visualizar la respuesta a un % escalón unitario, para los % primeros 20 segundos

Representación temporal de la respuesta a una entrada escalón unitario.


Página 43


3.3

2011-III

Estructura de filtros digitales. El segundo modelo de representación de una ecuación en diferencia con coeficientes constantes, corresponde al modelo auto-regresivo (AR). A continuación se presenta el diagrama de bloques correspondiente y su análisis a través del Matlab (ver figura 6.5). Sea el siguiente ejemplo de función de transferencia, se solicita analizar el diagrama de bloques y posteriormente la respuesta ante una entrada escalón unitario. H ( z ) 

Y ( z ) X ( z )



4 1

2

3

1  0.2 z  0.71 z  0.072 z

4 X (z)

Y (z)

∑

-0.2

Z-1

∑

0.71

Z-1

∑ Z-1 0.072

Figura 6.5

Estructura de filtros digitales correspondiente al modelo AR.

>> Nu = 1; >> De = [ 1 0.2 -0.71 -0.072 ]; >> zplane( Nu , De )

%ver figura 6.6

Para visualizar la respuesta a un escalón unitario, se debe de restringir el número de muestras de la señal de salida, debido a que este tipo de modelo corresponde a un IIR (respuesta infinita al impulso). >> x = [ 1zeros( 1 , 40 ) ]; >> y = filter(Nu , De , x ); >> stem( 0:length(y)-1 , y )


%ver figura 6.7

Página 44


2011-III

>> xlabel(' tiempo ' ) >> ylabel(' Amplitud ')

Figura 6.6

Diagrama de polos y ceros de la función de transferencia H(z) del ejemplo anterior.

Figura 6.7

Representación temporal de la respuesta a una entrada escalón unitario.

4. EJERCICIOS POR SOLUCIONAR: 4.1 Obtenga la Transformada Z Inversa de H(z) haciendo uso de la expansión en fracciones parciales y con ayuda del Toolbox Symbolic del Matlab. Posteriormente, grafique h[n] para el intervalo de tiempo: 0 ≤ n ≤ 12. H ( z ) 

Y ( z ) X ( z )

1




3

4 z  7 z

1  0.4 z   0.72 z   0.142 z   0.0839 z   0.0063 z  1

2

3

4

5

Página 45


2011-III

4.2

Para la pregunta anterior, considere que la entrada X(z) es igual a: 1 + 2z –3. Por lo tanto, se solicita determinar la salida Y(z) haciendo uso de la expansión en fracciones parciales. Luego, con ayuda del Toolbox Symbolic del Matlab, obtenga su respectiva Transformada Z Inversa. Finalmente, represente la señal de salida y[n] en el dominio del tiempo discreto para el intervalo de 0 ≤ n ≤ 25.

4.3

Para la siguiente función de transferencia (modelo ARMA), grafique su correspondiente diagrama de bloques. Asimismo, deberá mostrar el diagrama de polos y ceros, así como la respuesta ante una entrada tipo escalón de amplitud igual a 4. H ( z ) 

Y ( z ) X ( z )

2




4

1  4 z  7 z 







1 2 3 4 1  0.45 z  0.5775 z  0.2206 z  0.0159 z

Página 46


2011-III


FILTROS DIGITALES NO RECURSIVOS: PASA-BAJO, PASA-ALTO, PASA-BANDA Y RECHAZA BANDA.

1. OBJETIVOS: 1.1 Diseñar filtros digitales no recursivos tipo Pasa-Bajo, Pasa-Alto, Pasa-Banda y Rechaza-Banda. 1.2 Realizar operaciones de filtrado digital.


Filtros digitales no recursivos (FIR) Son sistemas discretos que se caracterizan por ser estables. Es decir, todos los polos de su función de transferencia se encuentran en el interior del círculo unitario. Asimismo, poseen una fase lineal y necesitan un orden elevado para aproximarse al tipo de filtro ideal. Seguidamente se muestra un ejemplo de función de transferencia de un filtro digital no recursivo, su respuesta impulsional y su correspondiente diagrama de bloques.

H ( z ) 

Y ( z ) X ( z )

 a  bz 1  cz 2  bz 3  az 4

h[n]  a [n]  b [n  1]  c [n  2]  b [n  3]  a [n  4]

Z-1

Z-1

X (z)

a

b

Z-1

c

b

∑


Z-1

a

Y (z)

Página 47


2.2

Filtro Pasa-Bajo, Pasa-Alto, Pasa-Banda y Recha-Banda Entre los filtros digitales no recursivos más comunes encontramos: Filtro Pasa-bajo: posee una frecuencia de corte y deja pasar frecuencias menores a esta. Filtro Pasa-alto: posee una frecuencia de corte y deja pasar frecuencias mayores a esta. Filtro Pasa-banda: posee dos frecuencias de corte y deja pasar frecuencias comprendidas entre estas dos frecuencias. Filtro Rechaza-banda: posee dos frecuencias de corte y rechaza frecuencias comprendidas entre estas dos frecuencias.

2.3

2011-III

todas las todas las todas las todas las

Diseño de Filtros por Windowing Consiste en la elección de una ventana (hamming, hanning, blackman, rectangular, triangular, káiser, etc.) para multiplicarla con la función SINC, correspondiente a un filtro ideal. Ventana Hamming

 2  n  M  1 

w[n]  0.53836  0.46164 * cos

0  n  M  1 M : tamaño de la ven tan a

Función SINC o filtro ideal en el dominio del tiempo discreto

hi [n]  sin c(n) 

sen( n)

 n

Por lo tanto, para diseñar el filtro digital no recursivo, se procede a multiplicar la ventana y el filtro ideal, en el dominio del tiempo discreto.

h[n]  w[n] * hi [n] 2.4

Operación de filtrado digital La operación de filtrado digital se realiza de muchas maneras. Una de ellas concierne a la operación de convolución entre la señal de entrada y la respuesta impulsionaldel filtro digital no recursivo, tal como lo muestra la figura 7.1.


Página 48


x [n]

Figura 7.1

h [n]

2011-III

y [n]

Representación de la entrada y salida de un filtro digital.

y[n]  x[n]  h[n]  h[n]  x[n] 3. EJERCICIOS SOLUCIONADOS: 3.1

Filtro digital no recursivo pasa-bajo. Para diseñar un filtro digital no recursivo pasa-bajo, se inicia señalando el orden del filtro a requerir, luego la frecuencia de corte, la frecuencia de muestreo y el tipo de ventana a utilizar. A continuación se muestra un ejemplo de diseño para un filtro FIR Pasa-Bajo con Fc=2 KHz, Fs = 10 KHz y una ventana haming. >> M = 66; >> Fc = 2000; >> Fs = 10000; >> fcN = Fc / Fs; >> wc = 2*pi*fcN; >> for n=0:M/2; if n==0 fi(n+1) = wc/pi; else fi(n+1) = (sin(wc*n))/(pi*n); end end >> hi = [ fliplr( fi(2:M/2+1)) fi ]; >> w = hamming( M+1 ); >> h = (w').*hi; >> stem( 0:length(h)-1 , h ) >> xlabel('Tiempo discreto') >> grid

% orden del filtro FIR % frecuencia de corte % frecuencia de muestreo % frecuencia normalizada

% filtro ideal en el tiempo discreto % ventana Hamming % filtro real tipo digital % filtro real en el dominio del tiempo % ver Figura 7.2

>> freqz(h , 1 , 1024 , 'whole' , Fs )

% filtro FIR en el dominio de la % frecuencia. Desde 0 a Fs. % ver Figura 7.3

>> zplane( h , 1 )


Página 49


Figura 7.2

Figura 7.3 3.2

2011-III

Representación del filtro digital FIR en el dominio del tiempo

Representación del filtro digital FIR en el dominio de la frecuencia

Filtro digital no recursivo pasa-alto. Para diseñar un filtro digital no recursivo pasa-alto, se inicia señalando el orden del filtro a requerir, luego la frecuencia de corte, la frecuencia de muestreo y el tipo de ventana a utilizar. Se debe partir del diseño de un filtro pasa-bajo, luego se transforma el filtro ideal. A continuación se muestra un ejemplo de diseño para un filtro FIR Pasa-Alto con Fc=3 KHz, Fs = 10 KHz y utilizando una ventana triangular. >> M = 86; >> Fc = 3000; >> Fs = 10000; >> fcN = Fc / Fs; >> wc = 2*pi*fcN;


% orden del filtro FIR % frecuencia de corte % frecuencia de muestreo % frecuencia normalizada

Página 50


>> for n=0:M/2; if n==0 fi(n+1) = 1 - wc/pi; else fi(n+1) = -(sin(wc*n))/(pi*n); end end >> hi = [ fliplr( fi(2:M/2+1)) fi ]; >> w = triang( M+1 ); >> h = (w').*hi; >> stem( 0:length(h)-1 , h ) >> xlabel('Tiempo discreto') >> grid

2011-III

% filtro ideal Pasa-Alto en el tiempo discreto % ventana Triangular % filtro real tipo digital % filtro real en el dominio del tiempo % ver Figura 7.4

>> freqz(h , 1 , 1024 , 'whole' , Fs )

% filtro FIR en el dominio de la % frecuencia. Desde 0 a Fs. % ver Figura 7.5

>> zplane( h , 1 )

Figura 7.4 3.3


Filtro digital no recursivo rechaza-banda. Para diseñar un filtro digital no recursivo rechaza-banda, se inicia señalando el orden del filtro a requerir, luego las frecuencias de corte, la frecuencia de muestreo y el tipo de ventana a utilizar. Se debe partir del diseño de dos filtros pasa-bajo, transformando luego el filtro ideal. A continuación se muestra un ejemplo de diseño para un filtro FIR rechazabanda con Fc 1=2.1 KHz, Fc2=3.5 KHz, Fs = 10 KHz y utilizando una ventana hanning. >> M = 66; >> Fc1 = 2100; >> Fc2 = 3500; >> Fs = 10000; >> fcN1 = Fc1 / Fs;


% orden del filtro FIR % frecuencia de corte 1 % frecuencia de corte 2 % frecuencia de muestreo % frecuencia normalizada 1 Página 51


>> fcN2 = Fc2 / Fs;

Figura 7.5

2011-III

% frecuencia normalizada 2


>> wc1 = 2*pi*fcN1; >> wc2 = 2*pi*fcN2; >> for n=0:M/2; if n==0 fi(n+1) = 1 - (wc2-wc1)/pi; else fi(n+1) = (sin(wc1*n)-sin(wc2*n))/(pi*n); end end >> hi = [ fliplr( fi(2:M/2+1)) fi ]; % filtro ideal en el tiempo discreto >> w = hanning( M+1 ); % ventana Triangular >> h = (w').*hi; % filtro real tipo digital >> stem( 0:length(h)-1 , h ) % filtro real en el dominio del tiempo >> xlabel('Tiempo discreto') % ver Figura 7.6 >> grid >> freqz(h , 1 , 1024 , 'whole' , Fs )

% filtro FIR en el dominio de la % frecuencia % ver Figura 7.7

>>zplane( h , 1 )


Página 52


Figura 7.6

Figura 7.7

2011-III



Del diagrama de polos y ceros del filtro digital es posible graficar la magnitud y fase respecto a la frecuencia. Para ello, se necesita subdividir el plano Z en un número de partes igual a una potencia de 2. Luego, dividir el espectro de frecuencia, desde 0 hasta Fs, por cada una de las partes subdividas sobre la circunferencia unitaria. Una vez realizada la subdivisión, se procede a tomar distancias, para el caso del módulo, o ángulos, para el caso de la fase, de cada punto sobre el círculo unitario a cada cero. El polo no se toma en cuenta para el módulo debido a que la distancia es igual a 1 por tratarse de una circunferencia de radio unitario.


Página 53


2011-III

A continuación se muestra un ejemplo del producto de las distancias para el caso de un filtro FIR de orden 2. >> Fc = 220; >> Fs = 1000; >> orden = 2; >> h = fir1( orden , Fc / (Fs/2), triang(orden + 1)); >> zplane(h,1) H ( z )  0.2077  0.5846 z 1  0.2077 z 2 H ( z )  0.2077(1  2.8146 z 1  z 2 )



Móduloi  20 * log10  H o



 dist cero , P i     dist  polo , P i  

Fasei   cero, P i     polo, P i  i : punto de subdivisión sobre el círculo unitario Donde: Ho: dist( cero , Pi ): dist( polo , Pi ): (cero, Pi ): (polo, Pi ):

Factorización de la función de transferencia. distancia de cada cero a un punto Pi del plano Z. distancia de cada polo a un punto Pi del plano Z. ángulo de cada cero a un punto Pi del plano Z. ángulo de cada polo a un punto Pi del plano Z.

Para el ejemplo anterior, los ceros son: Cero1 = -2.3973 Cero2 = -0.4171 Y los polos se encuentran en el origen, por lo tanto su distancia a cada punto P i sobre el círculo unitario del plano Z, siempre será 1. Entonces: Ho = 0.2077


Página 54




Móduloi  20 * log10  H o



2011-III

 dist cero , P i     dist  polo , P i  

Fasei   cero, P i     polo , P i  i : punto de subdivisión sobre el círculo unitario Para : F  0 Hz d 01  dis ([2.3973,0],[0,0]) d 0 2  dis ([0.4171,0],[0,0]) Para : F  125 Hz d 12 51  dis ([2.3973,0], [0.7071,0.7071]) d 12 5 2  dis ([0.4171,0],[0.7071,0.7071]) Para : F  250 Hz d 25 01  dis ([2.3973,0], [0,1]) d 25 0 2  dis ([0.4171,0],[0,1]) Para : F  375 Hz d 37 51  dis ([2.3973,0], [0.7071,0.7071]) d 37 5 2  dis ([0.4171,0],[ 0.7071,0.7071]) Para : F  500 Hz d 50 01  dis ([2.3973,0], [1,0]) d 50 0 2  dis ([0.4171,0],[1,0]) Módulo0 Hz  20 * log10 0.2077 * d 01 * d 0 2  Módulo125 Hz  20 * log10 0.2077 * d 12 51 * d 12 52  Módulo25 0 Hz  20 * log10 0.2077 * d 25 01 * d 25 0 2  Módulo375 Hz  20 * log10 0.2077 * d 37 51 * d 37 52  Módulo50 0 Hz  20 * log10 0.2077 * d 50 01 * d 50 0 2 


Página 55


2011-III

Lo mismo se realiza para los otros puntos de frecuencia (P i) correspondientes a: 625 Hz, 750 Hz y 875 Hz, aunque no es necesario calcularlos debido a la simetría cuando se grafica el módulo. De manera similar se procede a obtener los ángulos de cada cero y polo hacia cada uno de los 8 puntos ubicados sobre la circunferencia unitaria, con la finalidad de obtener el gráfico de fase.

Figura 7.8 3.4

Representación del filtro digital FIR en plano Z

Operación de filtrado digital. El filtrado de una señal digital con cualquiera de los filtros anteriormente indicados, es posible obtenerlos de muchas maneras. A continuación se señalan dos formas simples y muy utilizadas. En primer lugar, el filtrado digital es posible obtenerlo a través de una convolución entre la señal de entrada discreta en el tiempo y el filtro digital expresado, también, en el dominio del tiempo. Dicha convolución se obtiene utilizando el comando CONV. Seguidamente un ejemplo. Sea una señal periódica compuesta por la suma de tres ondas: coseno de 250 Hz, seno de 350 Hz y coseno de 500 Hz. >> Fs = 4000; >> n = linspace( 0 , 1 , Fs ); >> x = 2*cos( 2*pi*n*250 ) + 5*sin( 2*pi*n*350 ) - cos( 2*pi*n*500); >> TX = fft( x , 4096 ); >> F = linspace( 0 , Fs , 4096 ); >> subplot( 2 , 2, 1) , plot( n, x ) %Señal de entrada en el tiempo >> subplot( 2 , 2, 2) , plot( F, abs(TX) ) %TDF de la señal de entrada


Página 56


2011-III

>> orden = 80; >> w = hamming( orden+1 ); >> Fc = [ 300 405 ] / (Fs/2);

%Ventana hamming. %Frecuencia de corte normalizada

>> h = fir1(orden , Fc , w ); >> y = conv( h , x ); >> nn = linspace( 0 , 1 , length(y) ); >> TY = abs( fft( y , 4096 ) );

%Diseño del filtro utilizando FIR1 %Filtrado digital o convolución solo %para secuencias finitas. %TDF de la señal de salida.

>> subplot( 2 , 2, 3) , plot( nn, y ) >> subplot( 2 , 2, 4) , plot( F, TY )

%Señal de salida en el tiempo %TDF de la señal de salida %Ver figura 7.9

Figura 7.9 a) Zoom de la señal de entrada en el tiempo. b) Representación en frecuencia de la señal de entrada. c) Zoom de la señal de salida en el tiempo. d) Representación en la frecuencia de la señal de salida. En segundo lugar, el filtrado digital, también es posible obtenerlo a través del comando FILTER donde participan tanto la señal de entrada como el filtro digital, en el dominio del tiempo discreto. A continuación, para el ejemplo anteriormente planteado, se realiza la siguiente modificación. >> Nu = h ; >> De = 1; >> y = filter( Nu , De , x); Donde “Nu” representa a los coeficientes del polinomio del numerador de la función de transferen cia del filtro digital diseñado. Asimismo, “De” corresponde


Página 57


2011-III

al denominador que en este caso es igual a “1” por tratarse de un filtro FIR,

cuya función de transferencia es del tipo:

H ( z ) 

Y ( z ) X ( z )

 bo  b1 z 1  b2 z 2  b3 z 3  . . .  b M z  M

4. EJERCICIOS POR SOLUCIONAR: 4.1 A partir del siguiente bosquejo, diseñe el filtro digital correspondiente haciendo uso de la ventana de Blackman. Utilizar Ws = 2000  rad/seg. Luego, grafique el filtro en el dominio del tiempo discreto, el diagrama de polos y ceros, y la representación en frecuencia desde –Fs/2 a Fs/2.

- 3 dB

-41 dB

0

400 440 620 680 1000 rad / se

4.2

Para el diseño de un filtro digital FIR Pasa-Banda se utilizó la ventana de hamming, y una frecuencia de muestreo de 4 KHz. Si la ganancia obtenida en una de sus frecuencias de corte de 500 Hz, fue de -0.843 dB, se solicita graficar el diagrama de polos y ceros de dicho filtros, así como su respuesta en frecuencia tanto de módulo como de fase (para el intervalo de –Fs/2 a Fs/2). H(z) = - 0.0462 + 0.9076 z-2 - 0.0462 z-4

4.3

Diseñar un filtro FIR Pasa-Banda de orden 100, muestreado a 10KHz y con frecuencias de corte igual 2KHz y 3.2 KHz. Utilice la ventana de hamming. Luego, se solicita construir una rutina de programación en el Matlab para graficar el espectro de frecuencia de módulo y fase utilizando el producto de distancias o suma de ángulos a partir del diagrama de polos y ceros. Procure que la atenuación en la banda pasante se sitúe entre 0 y 3 dB y utilice por lo menos 512 muestras.


Página 58


2011-III


FILTROS DIGITALES RECURSIVOS. TRANSFORMACIÓN DE FRECUENCIA.

1. OBJETIVOS: 1.1 Diseñar filtros digitales recursivos tipo Pasa-Bajo, Pasa-Alto, Pasa-Banda y Rechaza-Banda. 1.2 Realizar operaciones de filtrado digital. 2. FUNDAMENTO TEÓRICO: 2.1

Filtros digitales recursivos (IIR) Son sistemas discretos que se caracterizan por presentar polos y ceros dispersos en todo el plano Z. Asimismo, poseen una fase no linealy no necesitan de un orden elevado para aproximarse al tipo de filtro ideal. En la figura 8.1 se muestra un ejemplo de una función de transferencia de un filtro digital recursivo y su correspondiente estructura discreta.

H ( z ) 

x [n]

Y ( z ) X ( z )

1



1

3

Z-1

a, b, c, d , e, f , g  

+

b

y [n]

+f

c

+

+

+ +

Z-1

-g

+ +

-1

+

+

+

+ Z

+

a

,

+ +e

+

2

w[n]

+

+

3

1  ez  fz  gz

+

+

2

a  bz  cz  dz

d

Figura 8.1. Estructura discreta de un filtro IIR.


Página 59


2.2

2011-III

Diseño por aproximación de la función Butterworth Se parte del diseño de filtros pasa-bajos analógicos, para luego transformarlo al dominio discreto. Seguidamente, la figura 8.2 muestra la repuesta en módulo de un filtro pasa-bajo IIR y sus respectivas bandas de frecuencia, donde:

1: 2: f p: f s: f n:

rizado u ondulación de la banda pasante. rizado u ondulación del para banda (banda de rechazo). frecuencia límite de la banda pasante. frecuencia límite del para banda (banda de rechazo). frecuencia normalizada ( fs / fp ).

Para realizar el diseño se normaliza la frecuencia de acuerdo a las especificaciones, se determina el orden del prototipo de filtro pasa-bajo, se obtiene la función de transferencia normalizada, y se finaliza desnormalizando a través de las transformaciones de frecuencia. A continuación se muestran las expresiones matemáticas necesarias para la realización del diseño:

 2  100.1 1 1

0.1  2

 10

log n

 0.1  1

 10

1/ 2

 1    1 

log f n 

H( f ) 

1

Banda de transición

1 - 1 Banda pasante

Para banda

2 f f p

fs

Figura 8.2. Respuesta en módulo de un filtro IIR pasa-bajo. 2.3

Transformaciones de frecuencia y transformación bilineal El diseño de filtros IIR parte de un prototipo pasa-bajo analógico. Por lo tanto, el orden del filtro se duplicará cuando se procede al diseño de filtros pasabanda o elimina-banda.


Página 60


2011-III

También es posible utilizar los comandos del Matlab: LP2LP, LP2BP, LP2BS y LP2HP, para lograr la transformación de frecuencia. Posteriormente, a través de la transformación bilineal (comando BILINEAR), se procede a la digitalización del filtro utilizando la siguiente expresión:

s 

2  z  1 





T s  z  1 

Por ejemplo, sea el filtro pasa-bajo analógico H(s), el cual se digitaliza a través de la Transformación Bilineal utilizando un periodo de muestreo igual a 2 segundos.

1 H ( s) 

2 s  5

 H(z) 

2 2  z  1 

 5 T s  z  1 

3

1

 z 1

3 2 1 1  z 3

>> Nc = [ 2 ]; >> Dc = [ 1 5 ]; >> Ts = 2; >> [ Nd , Dd ] = bilinear( Nc , Dc , 1/Ts ) Nd = 0.3333 0.3333 Dd = 1.0000 0.6667

2.4

Filtrado digital y respuesta impulsional Para determinar la respuesta impulsional, h[n], se procede a colocar como entrada al filtro digital, IIR, la función [n]. También es posible obtenerlo con ayuda del comando “DIMPULSE”.

Para realizar la operación de filtrado, se recurre a la representación en ecuaciones de diferencia, a partir dela estructura del filtro mostrado en la figura 8.1. Donde w[n] denota a una variable auxiliar en el tiempo “n”, y w[n -1], w[n-2], w[n-3] denotan el tiempo pasado de dicha variable auxiliar w[n]. Para el instante n=0, w[n-1] = w[n-2] = w[n-3] = 0. Asimismo, también es posible obtener la operación de filt rado utilizando el comando “FILTER”.

w[n]  x[n]  e * w[n  1]  f * w[n  2]  g * w[n  3] y[n]  a * w[n]  b * w[n 1]  c * w[n  2]  d * w[n  3]


Página 61


2011-III


Filtro digital recursivo pasa-bajo (Butterworth) Se desea diseñar un filtro pasa-bajo digital a partir de la función analógica Butterworth. Para ello, se parte de un filtro analógico con las siguientes características: atenuación de pasa-banda menor a 1dB a una frecuencia de 0.70 Hz, atenuación de para-banda mayor a 30 dB para una frecuencia de 1.45 Hz, y un periodo de muestro igual a 10 mili segundos. Calculamos n y

 2 . Luego, normalizamos la frecuencia.

>> Fp = 0.70; >> Fs = 1.45; >> d1 = -1; >> d2 = -30;

% atenuación en la banda pasante (en dB) % atenuación en la banda de rechazo (en dB)

 2

>> e2 = power( 10 , - 0.1*d1 ) - 1 ; % variable >> Fn = Fs / Fp; >> n = 0.5 * log10( ( 10^(-0.1*d2) - 1 )/( 10^(-0.1*d1) - 1 ) ) / log10( Fn ) >> n = ceil( n ); >> k = 1:n; >> u = -sin( (2*k - 1 ) * pi / 2 / n ); >> w = cos( (2*k - 1 ) * pi / 2 / n ); >> polos = u + w*j; >> num = 1; >> den = real( poly( polos ) ); >> figure(1) >> pzmap( num ,den) % diagrama de polos y ceros >> freqs( num , den ) % respuesta en frecuencia del filtro % normalizado analógico También es posible obtener los polos del filtro normalizado H N(s), a partir de la función BUTTAP y conociendo el orden del filtro. >> [ z , p , k ] = buttap( n ); >> num1 = poly( z ); >> den1 = poly( p ); >> figure(2) >> pzmap( num1 , den1 )

% diagrama de polos y ceros

Se obtiene: H N ( s) 

1 s 6  3.86s 5  7.46s 4  9.14s3  7.46s2  3.86s  1

Se desnormaliza el filtro normalizado H N(s), y se obtiene H D(s). Finalmente se digitaliza dicho filtro.


Página 62


2

H D ( s)  H N ( s w p * v3 )

H (v3 ) 

1 1   v3 2

2*n



2011-III

1 2

Donde: HD(s): representa al filtro analógico desnormalizado. >> e2 = 0.2589; >> n = 6; >> v3 = solve(' 1/( 1 + 0.2589 * power( v3 , 2*6 )) = 0.5 '); >> v3 = 1.1191947453492187293069439198769 >> wp = 2*pi*Fp; >> wo = wp * v3; >> [ num2 , den2] = lp2lp ( num , den , wo); >> pzmap( num2 , den2 )

%diagrama de polos y ceros del filtro % pasa-bajo analógico desnormalizado % respuesta en frecuencia

>> freqs( num2 , den2 )

El filtro analógico pasa-bajo desnormalizado H(s) resultante, será:

H D ( s) 

14226.609962 s  19.02s  180.86s  1090.37s 3  4382.40s 2  11166.62s  14226.61 6

5

4

>> Ts = 0.01; >> Fs = 1 / Ts; >> [ numD , denD ] = bilinear ( num2 , den2 , Fs ); >> zplane( numD , denD ) >> puntos = 1024; >> freqz( numD , denD , puntos , 'whole', Fs )

% digitalizando % ver la figura 8.3

El filtro digital pasa-bajo H(z) resultante, será: 9

H ( z ) 



 3.03 z 2  4.04z -3  3.03 z 4  1.21z -5  0.20 z 6  -1 -2 -3 -4 -5 -6 1 - 5.81z  14.07z - 18.17z  13.20z - 5.12z  0.83z

10 0.20  1.21z


-1

Página 63


2011-III

Figura 8.3.Respuesta en módulo del filtro pasa-bajo digital IIR 3.2

Filtro digital recursivo pasa-alto (Chebyshev I y II) A continuación se muestran las características de frecuencia y atenuación para el diseño de un filtro digital IIR pasa-alto, partiendo de la aproximación de la función analógica Chebyshev tipo I: Rp = 3 dB : Rs = 40 dB : Wp = 1200  rad/seg. : Ws = 700  rad/seg. :

atenuación en la banda pasante. atenuación en la banda de rechazo. frecuencia límite en la banda pasante. frecuencia límite en la banda de rechazo.

>> Rp = 3; >> Rs = 40; >> Fs = 10000; >> Ws = 2 * pi * Fs ; >> wp = 1200*pi / ( Ws/2);%frecuencia digital normalizada >> ws = 700*pi / ( Ws/2); %frecuencia digital normalizada >> [ orden, wc ] = cheb1ord( wp, ws, Rp, Rs); >> [ nu2 , de2 ]= cheby1( orden , Rp , wc , 'high'); >> puntos = 1024; >> freqz( nu2 , de2 , puntos , 'whole', Fs ) % ver la figura 8.4


Página 64


Figura 8.4

2011-III

Respuesta en frecuencia (módulo y fase) del filtro digital pasa alto a partir de la función Chebyshev I

Asimismo, se muestran las características de frecuencia y atenuación para el diseño de un filtro digital IIR pasa-banda, partiendo de la aproximación de la función analógica Chebyshev tipo II: Rp = 3 dB : Rs = 38 dB : Wp = [3400 4200] rad/seg.: Ws = [ 1800 5600 ] rad/seg.:

atenuación en la banda pasante. atenuación en la banda de rechazo. frecuencias límites en banda pasante. frecuencias límites en banda- rechazo.

>> Rp = 3; >> Rs = 38; >> Fs = 10000; >> Ws = 2 * pi * Fs; >> wp = [ 1400*pi 4200*pi ] / ( Ws /2 ); >> ws = [ 800*pi 4850*pi ] / ( Ws /2 ); >> [ orden, wc ] = cheb2ord( wp, ws, Rp, Rs); >> [ nu2 , de2 ]= cheby2( orden , Rs , wc); >> puntos = 1024; >> freqz( nu2 , de2 , puntos , 'whole', Fs )

3.3

%frecuencia digital normalizada %frecuencia digital normalizada

%ver figura 8.5

Operación de filtrado digital Se diseña un filtro pasa-bajo digital Butterworth de orden 4 y frecuencia de corte de 200 Hz. A continuación se procede a filtrar una señal cuadrada con frecuencia principal igual a 100 Hz, y muestreada a 2 KHz. >> Fs = 2000; >> Fc = 200; >> orden = 4; >> puntos = 1000; >> [ Nu , De ] = butter( orden , Fc / (Fs/2) )


Página 65


2011-III

>> freqz( Nu , De , puntos , ‘whole’ , Fs )

>> n = 0: Fs - 1; >> x = square( 2* pi * 100 * n / Fs ); >> subplot( 2 , 1 , 1 ), plot( n , x ) >> axis( [ 200 360 -1.5 1.5 ] ) >> y = filter( Nu , De , x ); >> subplot( 2 , 1 , 2 ), plot( y ) >> axis( [ 200 360 -1.5 1.5 ] )

Figura 8.5

% filtrado de la señal cuadrada % ver la Figura 8.6

Respuesta en frecuencia (módulo y fase) del filtro digital pasa banda a partir de la función Chebyshev II


Diseñar un filtro analógico pasa-alto Butterworth con las siguientes especificaciones: atenuación de pasa-banda hasta 1.8 dB en la frecuencia de 20.8 Hz y atenuación de para-banda de 31 dB en la frecuencia de 12.6 Hz. Posteriormente, digitalizar dicho filtro utilizando un periodo de muestreo de 100 mili segundos. Finalmente, graficar su correspondiente diagrama de polos y ceros, su respuesta impulsional y la representación en frecuencia (magnitud y fase) desde –Fs/2 a Fs/2.


Página 66


(a)

Figura 8.6. 4.2

2011-III

(b)

a) Respuesta en frecuencia del filtro diseñado. b) Señal de entrada y salida al filtro digital diseñado.

A partir del siguiente bosquejo, diseñe el filtro digital correspondiente haciendo uso de la función Chebyshev Tipo II. Utilizar Ws = 2000  rad/seg. Luego, grafique el filtro en el dominio del tiempo discreto, el diagrama de polos y ceros, y la representación en frecuencia desde –Fs/2 a Fs/2.

0 dB - 3 dB

-41 dB

0

400

440

620

680

1000

rad / seg

4.3

Diseñar un filtro digital IIR Pasa-Banda de orden 8, a partir de la función Chebyshev I. Deberá de contar con un rizado menor a 2 dB en la banda pasante.El filtro debe encontrarse muestreado a 10KHz y con frecuencias de corte igual 2KHz y 3.2 KHz.


Página 67


2011-III


INTRODUCCIÓN AL FILTRADO ADAPTATIVO. ALGORITMO LEAST MEAN SQUARE

1. OBJETIVOS: 1.1 Utilizar el algoritmo LMS (least mean square) para diseñar filtros adaptivos. 1.2 Realizar una operación de filtrado adaptativo. 2. FUNDAMENTO TEÓRICO: 2.1

Filtro adaptivo Son sistemas discretos que se auto diseñan a través de un algoritmo recursivo, cambiando sus parámetros de forma dinámica a través del tiempo. Para ello, emplean un criterio de adaptación que se encuentra relacionado con la obtención de una señal de error existente: entre la señal de salida y la deseada. Estos sistemas adaptivos son utilizados en diferentes tipos de aplicaciones: a) Predicción. b) Modelado. c) Eliminación de interferencias y de eco. d) Identificación de sistemas. e) Control adaptivo. A continuación, en la Figura 9.1, se muestra un diagrama de bloques simplificado de un filtro adaptivo.

Figura 9.1.


Diagrama de bloques de un filtro adaptivo.

Página 68


2.2

2011-III

Algoritmo LMS El algoritmo LMS (least mean square), es un algoritmo simple que se utiliza en la teoría de filtrado adaptivo para hallar los coeficientes del filtro que minimice el valor cuadrático medio de la señal error. Esta señal error se define como la diferencia entre la señal deseada y la señal obtenida a la salida del filtro. Su principal desventaja se presenta en la tasa de convergencia lenta. A continuación se muestra un resumen del algoritmo LMS. 1. Se inicializan los pesos o coeficientes del filtro. w[n] 2. Se obtiene la salida del filtro. y [ n ] = w[ n ] * x[ n ] 3. Se calcula el error. e [ n ] = d[ n ] - y[ n ] 4. Se adaptan los coeficientes. w[ n + 1 ] = w[ n ] + u e[ n ] * x[ n ]


Filtrado adaptivo para identificación de sistemas Seguidamente se muestra un ejemplo aplicativo del uso del filtrado adaptivo para el caso de identificación de sistemas. Por lo tanto, se diseña, a priori, un filtro FIR pasa alto de orden 28, frecuencia de corte de 100 Hz y frecuencia de muestreo de 1000 Hz. Luego, se hace pasar una señal de error a través del filtro diseñado dando origen a una señal de salida que se denominará señal deseada d(t). Se seleccionan 41 coeficientes para el filtro adaptivo y se procede a realizar la operación de adaptación. >> x = 0.2*randn(1,600); >> Fc = 100; >> Fs = 1000; >> h = fir1( 40 , Fc / (Fs/2) , 'high' ); >> d = filter( h ,1,x); >> w0 = zeros(1,41); >> mu = 0.8; >> S = initlms(w0,mu); >> [y,e,S] = adaptlms(x,d,S);


%entrada al filtro %Filtro pasa alto FIR %por identificar %Señal deseada %Coeficientes iniciales %del filtro adaptivo %Tasa de aprendizaje %del filtro LMS %inicializa filtro adaptivo %filtrado adaptivo

Página 69


2011-III

Según el diagrama de bloques mostrado en la figura 9.2, cuando las señales “y” y “d” son próximas una de la otra, el filtro adaptivo consiguió identificar al filtro real denotado por “h”.

FILTRO PASA ALTO FIR DE ORDEN 20

+

Error: e(t)

+ -

Ruido: x(t) FILTRO ADAPTIVO LMS

Figura 9.2

Señal deseada: d(t)

Señal de salida: y(t)

Diagrama de bloques de filtrado adaptivo para identificación de sistemas.

>> subplot( 3, 1, 1 ) , plot( 1:600 , y ), grid >> title('Señal de salida del filtro adaptivo') >> subplot( 3, 1, 2 ) , plot( 1:600 , d ) , grid >> title('Señal de salida del filtro a identificar') >> error = d - y; >> subplot( 3, 1, 3 ) , plot( 1:600 , error ) >> title('Señal de error')

Figura 9.3

%ver figura 9.3

Representación gráfica de las señales a) salida del filtro adaptivo. b) salida del filtro identificado. c) señal de error.


Página 70


2011-III

Posteriormente, se muestran los coeficientes del filtro real identificado, y los coeficientes del filtro adaptivo utilizado. >> figure(2) >> hold on >> stem( 0:40 , h , 'or') , grid >> stem( 0:40 , S.coeffs , '*b'); >> hold off >> legend(' Filtro Real ','Filtro Estimado');

Figura 9.4 3.2

%ver figura 9.4

Identificación del filtro pasa alto a través del filtrado adaptivo

Filtrado adaptivo para cancelación de ruido A continuación se muestra un nuevo ejemplo aplicativo del uso del filtrado adaptivo para el caso de cancelación de ruido. Por lo tanto, se plantea, a priori, una señal coseno de 5 Hz con ruido, y muestreada a 80 Hz. Luego, se seleccionan 80 coeficientes para el filtro adaptivo y se procede a filtrar una nueva señal de ruido con la finalidad de cancelar ambos ruidos y tener como señal de error, una señal muy próxima a la señal coseno pero sin ruido. A continuación, en la figura 9.5, se muestra el diagrama de bloques correspondiente:

Coseno + Ruido Señal deseada: d (t)

+

Error e(t):

+

Coseno

Señal de Ruido x (t)

Figura 9.5

FILTRO ADAPTIVO LMS

Señal de salida: y(t)

Diagrama de bloques de filtrado adaptivo para cancelación de ruido.


Página 71


>> Fs = 80; >> t = linspace( 0 , 1, Fs); >> d1 = 2*cos( 2 * pi * 5 * t ); >> d = 0.2*randn( size(t) ) + d1; >> x =0.2*randn( size(t) );

>> w0 = zeros(1,80); >> mu = 0.05;

>> S = initlms( w0 , mu ); >> [y,e,S] = adaptlms(x,d,S); >> subplot(2, 1, 1), plot( d1 ), grid >> title('Señal COSENO sin ruido') >> subplot(2,1,2), plot( e ), grid >> title('Señal de ERROR (cancelación de ruido)')

Figura 9.6

2011-III

%señal limpia %señal deseada %coseno + ruido %señal de ruido %(entrada al filtro %adaptivo) %Coeficientes iniciales %del filtro adaptivo %Tasa de aprendizaje %del filtro LMS

%inicializa filtro adaptivo %filtrado adaptivo

%ver figura 9.6

Señal COSENO sin ruido y señal de ERROR producto de la cancelación de ruido

También es posible observar el error encontrado en la señal COSENO sin ruido y la señal donde se realizó la cancelación del ruido. Para ello, se debe realizar la siguiente operación: >> error = d1 - e; >> figure(2) >> plot( error ) MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo – FIEM – UTP

Página 72


3.1

2011-III

Algoritmo de Filtrado adaptivo Utilizando los comandos del software Matlab, se simplifica el uso del algoritmo de filtrado adaptivo. Por lo tanto, a continuación se representa el algoritmo de filtrado adaptivo LMS utilizando la teoría indicada en la sección 2.2, con la finalidad de entenderla y poder programarla en cualquier otro lenguaje de programación de alto nivel. Se aplica este nuevo ejemplo para el caso de identificación de sistemas. Entonces, como señal deseada se considera la salida del filtro pasa bajo con frecuencia de corte igual a 200 Hz y muestreo de 1000 Hz, previamente diseñado. >> N = 10000; >> Fs = 1000; >> h = fir1( 28 , 200 / (Fs/2) ); >> X = randn(1,N); >> D = filter( h , 1 , X ); >> M = 28; >> mu = 0.02; >> Y = zeros(1,N); >> E = Y; >> w = [1 zeros(1,M-1) ]'; >> pr = 0:-1:-M+1;

%tamaño del vector de entrada

%número de coeficientes del %filtro adaptivo %factor de aprendizaje %salida del filtro adaptivo %re-ordenando las muestras

>> for n = M:N XR = X(n+pr); Y(n) = XR * w; E(n) = D(n) - Y(n); w = w + ( mu*E(n)*XR )'; end

%algoritmo iterativo

>> [ H1 , W ] = freqz( h , 1 , 1024 , Fs ); >> [ H2 , W ] = freqz( w , 1 , 1024 , Fs );

%respuesta en frecuencia

>> mH1 = 20*log10( abs( H1 ) ); >> mH2 = 20*log10( abs( H2 ) ); >> fH1 = unwrap( angle( H1 ) ); >> fH2 = unwrap( angle( H2 ) ); >> subplot(2,1,1) >> plot( W , mH1 , 'r' , W , mH2 , '--b' ), grid %Magnitud >> xlabel(' Frecuencia (Hz) ') %ver figura 9.7-a >> ylabel(' Magnitud (dB) ') >> legend(' Magnitud del Filtro Real ', 'Magnitud del Filtro Estimado' ); >> subplot(2,1,2) >> plot( W , fH1 , 'r' , W , fH2 , '--b' ), grid %Fase >> xlabel(' Frecuencia (Hz) ') %ver figura 9.7-b >> ylabel(' Fase (º) ') >> legend(' Fase del Filtro Real ', 'Fase del Filtro Estimado' ); MSc. Pedro Freddy Huamaní Navarrete Docente Tiempo Completo – FIEM – UTP

Página 73

Guia Laboratorio PDS UTP

Recommend Documents