Apuntes de Procesamiento Digital de Señales

UNIVERSIDAD NACIONAL DE SAN AGUSTIN DE AREQUIPA FACULTAD DE INGENIERIA DE PRODUCCION Y SERVICIOS DEPARTAMENTO ACADEMICO DE INGENIERIA ELECTRONICA ESCUELA PROFESIONAL DE INGENIERIA DE TELECOMUNICACIONES

Apuntes de Aula Wildor Ferrel Serruto

Curso: Procesamiento Digital de Señales

Profesor: Wildor Ferrel Serruto

Procesamiento D Digital d de S Señales 1. Introducción Def inición El Procesamiento Digital de Señales (PDS) es la disciplina que estudia los fundamentos matemáticos y algorítmicos del tratamiento de señales y de la información que contienen las señales utilizando un sistema electrónico digital como por ejemplo, un computador, un DSP, un FPGA.  Procesamiento: Realización de operaciones de acuerdo a un algoritmo para transformar los datos o extraer información de ellos.  Digital: Sistema electrónico digital como un computador, un DSP, un FPGA.  Señal: Magnitud variable por medio de la cual se transmite información. Hacer procesamiento digital de señales significa Procesamiento Realizar operaciones o transformaciones Digital mediante un computador u otro circuito electrónico digital de Señales sobre funciones del tiempo y/o del espacio.

Clasif icación d de llas sseñales Por el tipo de función y el tipo de variable   

Señal analógica – Función Función continua de variable continua. Señal de tiempo discreto – Función Función continua de variable discreta. Señal digital – Función Función discreta de variable discreta.

Por el número de dimensiones:   

Señal unidimensional – Ejemplo: Ejemplo: señal de voz, st  Señal bidimensional – Ejemplo: Ejemplo: imagen, s x, y  Señal multidimensional – Ejemplo: Ejemplo: Señal de video en blanco y negro v x, y, t  , señal de video a color u x, y, t   r  x, y, t  g  x, y, t  b x, y, t 

Fundamentos d del P PDS ((Bases m matemáticas) 





Modelado de señales analógicas: o Transformada de Fourier, transformada de Laplace, filtros analógicos. Operaciones de convolución y correlación. Modelado de señales y sistemas de tiempo discreto: o La Transformada Z, la transformada de Fourier de tiempo discreto (DTFT), la transformada discreta de Fourier (DFT). o Operaciones de convolución y correlación. o Estructuras básicas de sistemas de tiempo discreto. Algoritmos de procesamiento digital de señales. o Reducción o incremento de la tasa de muestreo. o Transformada rápida de Fourier (FFT). o Diseño de filtros digitales. 1



Implementación d del P PD S Por software, en computadores de propósito general (Ejemplo: en una PC)  Por hardware (Ejemplo: con FPGAs)  Por software más hardware específico para PDS (Ejemplo: con un DSP TMS320C) Un procesador digital de señal (DSP-Digital Signal Processor), es un dispositivo que implementa hardware especializado para acelerar la ejecución de los algoritmos al goritmos de procesamiento digital de señales. 

Venta jas   

Inmunidad a ruido (mayor precisión). Implementación por software (mayor flexibilidad). Realización de funciones que no son posibles en procesamiento analógico de señales

Desventa ja 

En algunas aplicaciones la desventaja pueden ser el mayor costo y/o el procesamiento lento.

Comparación del modelado de un sistema de tiempo continuo y un sistema d de ttiempo d discreto Sistema de tiempo continuo

x( x(t )

y( y(t )

Filtro Analógico

x( x(t )

y( y(t )

t

t

Ecuación diferencial y t

c1

d y t 

d t Función de transferencia H s

c N

d N y t

d 0 x t

d t Y s X s

d 1

d x t 

d t

d 1 s  d M s M

d 0

1 c 1 s  c N s

2

N

d M

d M x t d t



Respuesta a una señal de entrada y t

x

 h t  d 

Respuesta en frecuencia H  j  

0



Sistema de tiempo discreto

x[n]

y[n]

Filtro Digital

x[n]

y[n]

0

0

n

n

Ecuación en diferencias lineales con coeficientes constantes y n a 1 y n 1  a N y n N b 0 x n b1 x n 1 Función de transferencia 1  M Y  z  b0  b1 z   b M z  H  z   X  z  1  a1 z 1   a N z  N Respuesta a una señal de entrada 



y n

x m h n m m

3

 b M x n M



Respuesta en frecuencia H  e j 

0

2



Competencias e especí f del c curso fi cas d 







Capacidad de analizar y especificar los parámetros fundamentales de un sistema de procesamiento digital de señales. Capacidad de analizar e identificar los principales componentes de un sistema de procesamiento digital de señales. Capacidad de aplicar el modelamiento matemático y algorítmico al procesamiento de señales. Comprensión y dominio de los conceptos básicos sobre las señales de tiempo discreto, los sistemas de tiempo discreto, las transformadas relacionadas y su aplicación en la resolución de problemas de procesamiento de señales.

Conocimientos p previos Fundamentos de cálculo. Algebra de números complejos. Análisis de Fourier de tiempo continuo. Filtros analógicos. Programación en MATLAB.

Bibliograf í ía [1] OPPENHEIM A. V., SCHAFER R.W. Segunda Edición. Tratamiento de Señales en Tiempo Discreto. Prentice Hall Iberia, Madrid, 2000 [2] PROAKIS J. G., MANOLAKIS D. G. Tercera Edición. Tratamiento Digital de Señales. Principios, algoritmos y aplicaciones. Prentice Hall, Madrid, 1998 [3] LI TAN. Digital Signal Processing. Fundamentals and Applications. Elsevier, DeVry University, Decatur, Georgia, 2008 [4] MARIÑO ACEBAL J. B. Segunda Edición. Tratamiento Digital de la Señal. Una introducción experimental Alfaomega, México, 1999

4



2. Conversor IIdeal d de T Tiempo C Continuo a a Tiempo D Discreto Diagrama d de B Bloques d de u un ssistema d de P Procesamiento D Digital d de S Señales Analógicas x sens(t )

x(t )

Filtro Antisolapamiento

Sensor

x[n]

xc(t )

y[n] DAC

DSP

ADC

yc(t ) Filtro de yr (t )

T

Reconstrucción

T

Procesamiento D Digital d de S Señales A Analógicas

x[n]

xc(t )

Sistema de Tiempo Discreto

C/D

y[n] D/C

T

T

Modelo d del C Conversor C C/D xc(t )

x s(t ) ×

Conversor de Tren de Impulsos Ponderados a Secuencia

s(t )

Dominio T Temporal

 t n T

s t n

x s

t

x c

yc(t )

t s t

5

x[n]


x s t


 t n T

x c t n

x c nT  t n T

x s t n

x s t

x n

 t n T

n

Dominio F Frecuencial S j

X s j

X s j

X s j

2

  k 2  T

T k 1 2

X c j

1

T k

1

T k

S j

X c j 

X c j

k

2

T

 s

;

2

T

 k  s

Señal c con E Espectro n no L Limitado X C  j 

1

-S

0

S



X S  j 

 

-2S

-S

0

S

6

2 S





Filtro A Antisolapamiento xc(t )

Filtro xa(t ) Antisolapamiento

Sistema de Tiempo Discreto

x[n] C/D

y[n]

yc(t ) D/C

( ) H aa j T

T

Ejemplo 1 La frecuencia de muestreo es 40 kHz. El filtro anti-solapamiento usado es un filtro de Butterworth pasa-bajas con frecuencia de corte de 8 kHz. El nivel de solapamiento en la frecuencia de corte debe ser 1%. Determine el orden del filtro antisolapamiento. Para el filtro de Butterworth la magnitud de la respuesta en frecuencia es:

| | 

   

H  j 

1

-S

C

0



S

El nivel de solapamiento en la frecuencia de corte es:

      √         



n 1 2 3 4 Nivel de solapamiento (%) 34.30 8.82 2.21 0.55

7



3. Señales d de T Tiempo D Discreto Secuencias b básicas Impulso unitario

 [n] 1

Escalón unitario

0

n

u[n] 1 ... n

0

Secuencia exponencial n x n  A  , n ; A,   R; A,   0 0    1

  1

x[n]

n

0 x[n]

Secuencia senoidal

0

n

x n  Asen 0 n  , n

x[n]

0

n

8



Operaciones b básicas Sean x1 n y x2  n dos secuencias, Adición : y n  x1 n  x2  n

Multiplicación de dos secuencias: y n  x1 n x2  n

Multiplicación por un escalar : y n   x 2  n

Producto interno : y  x1  n  x2  n

  x1  n x 2   n n

Desplazamiento : nd  0 y n  x1  n  nd 

x1[n]

0

n x1[n - nd ]

n

0 x1[n + nd ]

n

0

Ejemplo 2 Dadas las secuencias x1 n y x2 n , escribir un programa en MATLAB para determinar y graficar: a) La secuencia x3  n que es la suma de las secuencias x1 n , x2 n . b) La secuencia x4  n como el producto de las secuencias x1 n , x2 n . 9



c) La secuencia x5  n obtenida mediante el desplazamiento de x1 n a la izquierda en 2 unidades de tiempo. x1[n]

x2[n]

4

4

3

3

2

2

1

1

0

0

n

%************************************************ % Programa Matlab para el Ejercicio 1 %************************************************ % Reinicializar el ambiente clear; clf; % Generar las secuencias x1=[0 0 0 1 2 3 4 3 2 1 0 0 0]; x2=[0 0 4 3 2 1 0 0 0 0 0 0 0]; % a) Adicion x3=x1+x2; % b) Multiplicacion x4=x1.*x2; % c) Dezplazamiento a la izquierda en 2 posiciones x5=zeros(1,13); nd=2; x5(1:13-nd)=x1(nd+1:13); %for m=1:13-nd % x5(m)=x1(m+2); %end; % Graficar x1,x2,x3 subplot(3,1,1); stem([-6:6],x1); ylabel('x1'); xlabel('n'); subplot(3,1,2); stem([-6:6],x2); ylabel('x2'); xlabel('n'); subplot(3,1,3); stem([-6:6],x3); title('Adicion'); ylabel('x3'); xlabel('n'); pause; % Graficar x1,x2,x4 subplot(3,1,3); stem([-6:6],x4); 10

n



title('Multiplicacion'); ylabel('x4'); xlabel('n'); pause; % Graficar x1,x5 subplot(2,1,1); stem([-6:6],x1); ylabel('x1'); xlabel('n'); subplot(2,1,2); stem([-6:6],x5); title('Dezplazamiento'); ylabel('x5'); xlabel('n');

R elaciones iimportantes 

La secuencia escalón unitario se expresa a través de la secuencia impulso unitario de la siguiente forma : n

u n

 m

m 

A su vez, la secuencia impulso unitario se expresa a través de la secuencia escalón unitario en la forma :

 n 

u n

u n 1

Toda secuencia puede ser expresada como una suma ponderada de impulsos unitarios : x m  n

x n

m

m

Clasif icación d de llas ssecuencias p por ssu e extensión Secuencia de extensión finita :  n1   / x n  0, n  n1 y  n2   / x n  0, n  n2

n2

n1 Secuencia de extensión infinita : a) Secuencia derecha

 n1   / x n  0, n  n1 y  n2   / x n  0, n  n2 ... n

n1

b) Secuencia izquierda  n1   / x n  0, n  n1 y  n2   / x n  0, n  n2

... 11

n2

n



c) Secuencia bilateral  n1   / x n  0, n  n1 y  n2   / x n  0, n  n2 ...

...

n

Secuencia p periódica xn es periódica

  N / xn  xn  N n

Periodicidad d de lla ssecuencia c cosenoidal En tiempo continuo, la función cos  0 t es periódica para cualquier valor real de la frecuencia. El periodo es: T

2

0

En tiempo discreto, la relación  0 n = cos  0 n N se cumple si  0 N 2 k , donde N y k son enteros Por lo tanto, la secuencia cos  0 n es periódica si y sólo si cos

2  k

0

para algún N y

k enteros;

Por ejemplo, la secuencia

N

caso contrario, cos  0 n no es periódica. cos

3 4

n

es periódica con periodo

N

8.

En cambio, la

secuencia cos 3 n no es periódica La secuencia compleja Ce j 0 n es periódica si e j 0 n N = e j 0 n . Esto se cumple si  0 N 2 k , donde N y k son enteros. Las exponenciales complejas con frecuencias  0 y  0 2  r son iguales.

Ba jas y y a altas f f recuencias En tiempo continuo, la función cos  0 t oscila más rápidamente a medida que aumenta la frecuencia. En tiempo discreto, puesto que las secuencias cos  0 n y cos  0 2  r n son iguales, las frecuencias  0 y  0 2  r son equivalentes. Para la secuencia cos  0 n cuando la frecuencia aumenta de 0 a  aumentan también las oscilaciones. Sin embargo, cuando la frecuencia aumenta de  a 2 las oscilaciones se hacen más lentas. Las frecuencias en la vecindad de  0 2  k se llaman bajas frecuencias, mientras que las frecuencias en la vecindad de  0 2  k  se dice que son altas frecuencias.

12



4. Sistemas d de ttiempo d discreto Un sistema de tiempo discreto es una transformación que hace corresponder a cada secuencia de entrada una secuencia de salida. y n  Tx n y n x n T Ejemplo. Sistema de diferencia regresiva y n  x n  x n  1 x[n]

-2 -1 0 1 2

y[n]

n

-1 0

n

Tipos d de ssistemas Sistema causal

La salida depende de valores pasados y/o del valor presente de la entrada. Sistema no causal

La salida depende de valores futuros. Ejemplo. Sistema de diferencia progresiva: y n  x n  1  x n Sistema estable.

A toda secuencia limitada de entrada le corresponde una secuencia limitada de salida.  L x x n  Lx , n   Ly

y n  Ly ,  n

Sistema inestable

Por lo menos, a una secuencia limitada de entrada le corresponde una secuencia ilimitada de salida. Ejemplo. Sistema acumulador : n

y n 

 x k 

k 

13



Sistema lineal

Sistema que cumple con el principio de superposición. T ax1  n  bx 2  n  aT x1  n   bT x 2  n

Sistema invariante en el tiempo

Un desplazamiento en el tiempo de la secuencia de entrada produce el mismo desplazamiento de la secuencia de salida.



T x n  y n  T x n  nd

  yn  n  d

R espuesta d de u un ssistema L LTI Respuesta de un sistema al impulso unitario h n  T  n En un sistema LTI, la respuesta a una secuencia de entrada se expresa: 

 x mh n  m

y n 

m

Sistema FIR : h n tiene extensión finita Sistema IIR : h n tiene extensión infinita

Cálculo d de lla ssalida p por c convolución x n

y n



h n y n

x m h n

m

m

Ejemplo 3 La secuencia x n es aplicada a la entrada de un sistema LTI. Encontrar la secuencia de salida, si la respuesta del sistema al impulso unitario es h n . h[n]

x[n] 4 3 2

1

1

0

0

n

Convolución llineal Es la operación efectuada sobre secuencias : 

x1  n x2  n 

 x  m x  n  m

m

En sistemas LTI tenemos : y n  x n h n 14

1

2

n



Propiedades de la convolución   

Conmutativa Asociativa Distributiva

Conexión en serie de sistemas LTI h2  n

h1 n

≡

h1 n h2  n

≡

h1  n  h2  n

Conexión en paralelo de sistemas LTI h1 n

+ h2  n

Convolución p por b bloques Método de solapamiento-suma

Definición de un bloque: x r n

x n rL , 0 n L 1 0,

c .c . 

Secuencia de entrada expresada a través de los bloques xn   xr n  rL r  0

Secuencia de salida: yn  xn* hn ;

Bloque de salida:

    yn    xr n  rL * hn   xr n  rL* hn ; r 0  r 0 

y r n

x r n

L P 1

L



*



yn 



 y n  rL r  0

r

h n 

P

Se observa que y r n rL se solapa con y r 1 n solapamiento se suman.

r 1 L , y que las muestras de

Ecuación e en d dif erencias llineales c con c coef icientes c constantes ((LCCDE) Muchos sistemas lineales e invariantes en el tiempo (sistemas LTI) se describen mediante la ecuación: N

M

 a yn  k    b xn  r  k

r

k  0

r  0

Si hacemos a0  1  y n  b0 x n  b1 x n  1 ...b M x n  M 

 a1 y n  1 ...a N y n  N 

15



Cálculo d de lla ssalida p por r recursión Ejemplo 4 2

Para el sistema descrito por yn 

2

 xn  1 

2  yn  1  yn  2

hallar la salida si xn   n, si las condiciones iniciales son: y 1  0 , y 2  0 Solución: y0  y1  y2  y3  y4 

y5 

2 2 2 2 2 2 2 2 2 2

2 2

0 

2 0 0  0

1 

2

2 0  0  2

0 

2

0 

2 1 

0 

2

0 

2 0 

2

 0 1

2

2



2

2

2

1  0

2

2 2

2

2

ym 

2

y 3 

2 2

 xn  1 

 xm  1 

2

2

 2

yn  2  

yn 

2

2  yn  1  yn

2  ym  1  ym  2

 xn  1 

2  yn  1  yn  2

2

2  y 2  y 1 

2

 x 2 

2 2

0 

2 00  0

Dada la ecuación en diferencias: y n

b 0 x n

b1 x n 1

a1 y n 1

... b M x n M

... a N y n N

(1)

Despejamos y n N : y n N

b0 a N 1

a N

x n

y n

b1 a N a1 a N

x n 1

y n 1

...

...

b M a N

a N 1 a N

x n M

y n

(2)

N 1

Se calculan las muestras de salida a partir de las condiciones iniciales: , y n 0 2 , y n 0 1 , y n 0 , , y n 0 N 1 , y n 0 N , y n 0 N 1 ,                          

(2)

Condiciones 16

(1)



Causalidad d de u un ssistema L LTI 0,

h n

n 0

Estabilidad d de u un ssistema L LTI Un sistema LTI es estable

h n n

Suficiencia

Demostraremos que, si

entonces el sistema es estable

h n n

Sea x n acotada, x n Luego, y n

.

Lx

k =

x k h n k

=

x n

k h k

x n

k

k

Lx

h k

k

k

y n

Lx

h k

Si

h k

entonces y n es acotada. En

h n

n . consecuencia, el sistema es estable. k

Necesidad

Demostraremos que, si

entonces el sistema no es estable.

h n n

Sea x n definida por: x n

h

n

h

n

, si h n

0

si h n

0

0,

x n

es acotada, ya que x n

La salida es: y n

1.

x k h n

k

k

Para n

0 , tenemos:

y 0

x k h k

y 0

k = k

h

h

k

h

k

h 2

k =

h

k

h

k k

k

k

Si

h n n

, entonces

h

.

k

k

Esto significa que y n no es acotada. En consecuencia, el sistema no es estable. 17



Ejemplo 5 La secuencia 1, 0  n  10, x n   c. c. 0,

es aplicada a la entrada de un sistema LTI con respuesta al impulso

 21 n , n  0 h n   .  0, c. c. a) Calcular analíticamente la secuencia de salida por convolución lineal. Sólo los gráficos obtenerlos con un programa en MATLAB. b) Compruebe su resultado con la función conv del MATLAB y muestre los gráficos de las secuencias.

%**************************************************** % Programa para el Ejercicio 2 a) %**************************************************** % Reinicializar el ambiente clear; clf % Formar la secuencia de salida con 20 muestras y0_9=2-(1/2).^[0:9]; y10_19=((2^10)-1)*((1/2).^[10:19]); y=[y0_9 y10_19]; % Graficar la secuencia de salida stem([0:19],y); title('Secuencia de salida'); ylabel('y[n]'); xlabel('n'); %***************************************************** % Programa para el Ejercicio 2 b) %***************************************************** % Reinicializar el ambiente clear; clf % Formar la secuencia de entrada x=[ones(1,10) zeros(1,10)]; % Formar la respuesta al impulso h=[(1/2).^[0:19]]; % Efectuar la convolucion yy=conv(x,h); % Tomar 20 muestras de la secuencia de salida y=yy(1:20); % Graficar x,h,y subplot(3,1,1); stem([0:19],x); grid; title('Secuencia de entrada'); ylabel('x[n]'); xlabel('n'); subplot(3,1,2); 18



stem([0:19],h); grid; title('Respuesta al impulso'); ylabel('h[n]'); xlabel('n'); subplot(3,1,3); stem([0:19],y); grid; title('Secuencia de salida'); ylabel('y[n]'); xlabel('n');

Ejemplo 6 Escriba un programa en MATLAB para efectuar las siguientes tareas : a) Generar una secuencia senoidal x n de 50 muestras con frecuencia 2 25 rad/muestra. b) Obtener una secuencia x1 n adicionando a la secuencia inicial un ruido aleatorio uniformemente distribuido en el intervalo -0,25 a 0,25. c) Obtener la secuencia y n recursivamente filtrando la secuencia x1 n con un sistema promedio móvil de tamaño 5.

%**************************************************** % Programa para el Ejercicio 3 %**************************************************** % Reinicializar el ambiente clear; clf % Formar la secuencia de entrada N=50; x=sin(2*pi*[0:N-1]/25); % Formar la secuencia con ruido aditivo x1=x+0.5*(rand(1,N)-0.5); % Calcular la secuencia de salida para el sistema de % promedio movil de tamano 5 for n=1:N y(n)=mean(x1(max(n-4,1):n)); end; % Graficar las secuencias subplot(3,1,1); stem([0:N-1],x); grid; ylabel('x[n]'); xlabel('n'); subplot(3,1,2); stem([0:N-1],x1); grid; ylabel('x1[n]'); xlabel('n'); subplot(3,1,3); stem([0:N-1],y); grid; ylabel('y[n]'); xlabel('n');

19



Correlación c cruzada Para las secuencias x n e y n , la secuencia r xy n

x k y k

n

r xy n

k

x k

n

y k

k

se llama correlación cruzada de x n e y n . La correlación cruzada de y n e x n es: r yx n

y k x k

n

=

y k k

k

Se cumple que

r xy n

r yx

n x k

n

.

Comparando la expresión de la convolución x n

y n

x k y n

k

con

k

r xy n

x k y k

n

podemos escribir

r xy n

x n

y

n

k

Autocorrelación

Para la secuencia x n la secuencia de autocorrelación es r xx n

x k x k n k

En MATLAB se usa la función xcorr: x=[ 4 2 1 0]; y=[-1 -1 1 1]; Rxy=xcorr(x,y) Ryx=xcorr(y,x) Rxx=xcorr(x) Ryy=xcorr(y)

Correlación de secuencias periódicas

Para las secuencias periódicas x n e y n la correlación cruzada se define en la forma: r xy n

lim M

M

1 2 M 1 k

x k y k n M

La autocorrelación de x n será: r xx n

lim M

M

1 2 M 1 k

x k x k n M

Si las secuencias x n e y n tienen un mismo perido igual a N, el promedio en un intervalo infinito es igual al promedio en un único intervalo mayor o igual al periodo, es decir, siendo M≥N: r xy n

Las secuencias

1 M 1 x k y k M k 0 r xy n

n

r xx n

1 M 1 x k x k M k 0

n

y r xx n son periódicas y tienen el mismo periodo N. 20



Aplicación de la correlación en la determinación de la distancia de blancos sn    xn  D  wn r sx n tiene un pico en n  D

Aplicación de la correlación en la identificación de una señal periódica oculta en una señal con ruido

Sea la secuencia y n de la forma: y n

x n

w n

x n es una secuencia periódica con periodo desconocido N , w n es un ruido aleatorio.

Se desea determinar el periodo N. Para ello calculamos la autocorrelación de y n : r yy n

=

1 M 1 y k y k M k 0

1 M 1

M k 0

x k

w k

1 M 1 = x k x k M k 0

+

1 M 1 w k x k M k 0

= r xx n

r xw n

n =

x k n

w k n

1 M 1 n + x k w k M k 0

n +

1 M 1 w k w k M k 0

r wx n

= n + n =

r ww n

Ejemplo 7 Escriba un programa en MATLAB que realice las siguientes tareas:   con 500 muestras. a) Genere la señal   b) Forme la señal adicionando a un ruido aleatorio, uniformemente distribuido, con amplitud 1 y una media de 0. c) Calcule como la autocorrelación periódica de . d) Halle como la autocorrelación periódica de . e) A partir de , determine el periodo de . f) Forme un tren de impulsos con el periodo encontrado. g) Determine como la correlación cruzada periódica de con el tren de impulsos. h) Determine como la auto-correlación periódica de .

[]     [ ]  [ ]

[] [] [] [] []

[]

 [ ]

[] []

 [ ] []

21



%******************************************** % Aplicacion de la correlacion en la % identificacion de una señal periódica % oculta en una señal con ruido %******************************************** close all; clear all; N=500; n=[0:N-1]; x=0.5*cos((pi/26)*n)+0.5*cos((pi/13)*n); w=2*(rand(1,N)-0.5); y=x+w; Rx=xcorr(x,'biased'); Ry=xcorr(y,'biased'); RRy=xcorr(Ry,'biased'); %Determinamos el periodo [max1,pos1] = max(RRy) [min2,pos2] = min(RRy(pos1+1:length(RRy))) [max3,pos3] = max(RRy(pos1+pos2+1:length(RRy))) periodo=pos2+pos3 % Formamos el tren de impulsos s=zeros(1,length(y)); m=1; while m


plot(nRx,Rx);grid; title('Rx - Autocorrellacion de x '); NRy=length(Ry); nRy=[0:NRy-1]-((NRy-1)/2); subplot(312); plot(nRy,Ry);grid; title('Ry - Autocorrellacion de y '); NRRy=length(RRy); nRRy=[0:NRRy-1]-((NRRy-1)/2); Ninicio=((NRRy-1)/2)-((NRy-1)/2); Nfin=((NRRy-1)/2)+((NRy-1)/2); subplot(313); plot(nRy,RRy(Ninicio:Nfin));grid; title('RRy - Autocorrelacion de Ry'); pause; subplot(311); plot(n,x);grid; title('Secuencia original'); subplot(312); NRys=length(Rys); Ninicio=(NRys-1)/2; plot(Rys(Ninicio:Ninicio+N-1));grid; title('Rys - Correlacion cruzada de y con el tren s'); subplot(313); NRRys=length(RRys); Ninicio=(NRRys-1)/2; plot(RRys(Ninicio:Ninicio+N-1));grid; title('RRys - Autocorrelacion de Rys');

23


Profesor: Wildor Ferrel Serruto Secuencia original

1 0.5 0 -0.5 -1

0

50

100

150

200

250

300

350

400

450

500

300

350

400

450

500

300

350

400

450

500

100

200

300

400

500

100

200

300

400

500

100

200

300

400

500

300

350

400

450

500

300

350

400

450

500

300

350

400

450

500

Secuencia con ruido 2 1 0 -1 -2

0

50

100

150

200

250 Tren de impulsos

1

0.5

0

0

50

100

150

200

250

Rx - Autocorrellacion de x 0.6 0.4 0.2 0 -0.2 -500

-400

-300

-200

-100

0 Ry - Autocorrellacion de y

0.6 0.4 0.2 0 -0.2 -500

10

x 10

-400

-300

-200

-100

-3

0 RRy - Autocorrelacion de Ry

5 0 -5 -500

-400

-300

-200

-100

0

Secuencia original 1 0.5 0 -0.5 -1

0

50

100

150

200

250

Rys - Correlacion cruzada de y con el tren s 2 1 0 -1

0

50

100

150

200

250 RRys - Autocorrelacion de Rys

1 0.5 0 -0.5

0

50

100

150

200

250

24



5. La T Transf ormada Z Z Cálculo p por m medio d de ttransf ormadas La secuencia exponencial compleja: x se llama auto secuencia porque:

z n y n

=

z n

n

A z n

T

h k x n

k

=

k

h k z k

=

h k z n k

k

k

=

H z

z n

k

h k z k

= H z z n H z es el autovalor, se llama función de transferencia del sistema y es la transformada Z de la respuesta al impulso y n

El cálculo de la salida se puede hacer como sigue: x[n]

X ( z )

h[n]

H ( z )



Y ( z )

y[n]

Transformadas usadas en el análisis y diseño de sistemas de tiempo discreto:  Transformada Z  Transformada de Fourier de tiempo discreto (DTFT)  Transformada de Fourier discreta (DFT)  Transformada rápida de Fourier (FFT) x [n]

X Y

h [n]

H

Transformada Z directa :

X  z  



 xn z

n

n  

Transformada Z inversa :

xn 

1





X  z  z n 1dz

j C 2

25

y [n]



Ejemplo 8 x [n]

0

2

1

3 1

2

0

2 3 X  z   x0 z  x1 z  x2 z   x3 z 

 2  z 0  1 z 1  3  z 2  2  z 3 X  z   2  z 1  3  z 2  2  z 3 Se debe indicar la región de convergencia (RC). RC: Plano z completo excepto z = 0

Ejemplo 9

x [n]

Hallar la transformada Z de la secuencia n

 8  xn    un  9 

...

Solución:  8  n X  z      un z n    9 





 8      z n n 0  9  

 8     z 1   n0  9 

X  z  

1 8 1 1  z 9

1 8

1

1  z 9

 

1 1  

  1

n

Im{z} 8 9

1

z

1

RC : z 

|

8

8 9

9

RC

Cálculo d de lla ttransf ormada Z Z iinversa 

 n 

n 0

n



n

0

n



Inspección mas propiedades. Expansión en fracciones parciales. Expansión en serie de potencias.

26

Re{z}



Inspección mas propiedades

Ejemplo 10 Encontrar la secuencia cuya transformada Z es : 1

X z z

1

,

RC : z

1 2

2

Solución

X  z  

1 1

1  z 1 2

1 1

1

1 n 2

1 1 z 2

x n

1 n 1 2

z

1

n

a u n 1

x n n d

u n

az

1

, z

X z z

a

n d

u n 1

Expansión en fracciones parciales

Ejemplo 11 Hallar la transformada z inversa de la función

Im z 

     ||   mediante la expansión en fracciones parciales. Solución

           

|

     



1 2

2

Re z 

RC

                                                   

[]  []  []  [] Ejemplo 12

Hallar la respuesta al impulso del sistema causal con función de transferencia

        Haga la expansión en fracciones parciales mediante la función residuez del MATLAB. 27



Solución Corremos el programa b =[1 -1 -6]; a =[1 -1.5 -1]; [R,P,K]=residuez(b,a)

En la ventana de comandos del matlab obtenemos: R = -0.8000 -4.2000 P = 2.0000 -0.5000 K = 6

      []  [] [] []

Luego, Por tanto,

Expansión en serie de potencias

Encontrar la secuencia cuya transformada Z es: 1

X z

1 az

, RC : z

1

a

Solución 1

az

1 1 a

1

a

1

a

1

z

a

1

1

z a

2

z

2

a

3

z

3



z 2

z a

z

2

a

2

z

2

a

2

z

2

a

3

3

z 3

a

z

3



A partir de la condición de la RC tenemos: a 1 z  1

2 2

3 3

X  z    a z  a z  a z   



 a

,

m1

X  z  



  a z n

n

n  1

x n

anu

n

a

n  m

m m

z

z

1

28

1



Propiedades d de lla R R egión d de C Convergencia Secuencia derecha :

z 0 RC

Si z

z 0

z RC

z 0

z RC

Secuencia izquierda :

z 0 RC

Si z

La RC no puede contener polos Secuencia derecha

pk  RC

RC : z  max pk

Im{z}

Re{z}

... n

n1

RC

Secuencia izquierda

RC : z  min pk

Im{z}

Re{z}

... Secuencia bilateral

...

n2

RC

n Im{z}

...

Re{z}

n

RC 29



Función d de T Transf erencia Transformada Z de su respuesta impulsional H  z   Z hn  A partir de la LCCDE Y  z  b0  b1 z 1    b M z  M H  z    X  z  1  a1 z 1    a N z  N Los polinomios del numerador y del denominador se representan respectivamente

H  z  

B( z) y A( z)

B z 

A z 

Estabilidad 



hn   

n  





n

hn z

n  

z 1



Im{z}

La RC de la función de transferencia de un sistema estable contiene a la circunferencia de radio unitario

1

Re{z}

RC Causalidad La RC de la función de transferencia de un sistema causal es el exterior de una circunferencia y contiene al punto z = . Im{z}

Re{z}

RC

30



Ejemplo 13 Un sistema LTI causal tiene función de transferencia :

H  z  

1

3  4 z 1

 1.5 z  2

1  3.5 z

1) Determine manualmente la RC y la respuesta impulsional. ¿Es estable el sistema? 2) Compruebe los resultados anteriores usando MATLAB Solución

H  z  

1

3  4 z 1

1  3.5 z

 1.5 z  2





1

3 1  43 z





3 z  z 

4 3



1  3 z  1  0.5 z    z  3 z  0.5 1

1

Ceros :

     Polos :    

Im{z}

La RC es :

0.5 4/3

Re{z} 3

RC

No es estable el sistema.

                                       []    []  [] %************************************************ % Determinar el diagrama de polos y ceros % y la respuesta impulsional % de un sistema de tiempo discreto %************************************************ % Reinicializar el ambiente close all; clear all; % Funcion de Transferencia 31



B = [3 -4 0 ]; A = [1 -3.5 1.5]; % Diagrama de polos y ceros ceros=roots(B) polos=roots(A) zplane(B,A); grid; pause; % Respuesta impulsional [h,n]=impz(B,A); subplot(2,1,1); stem(n,h); grid; ylabel('h[n]'); xlabel('n'); title('Respuesta Impulsional segun impz'); % Respuesta impulsional analitica %n=[0:11]; hh=((1/2).^(n))+2*((3).^(n)); subplot(2,1,2); stem(n,hh); grid; ylabel('hh[n]'); xlabel('n'); title('Respuesta Impulsional analitica');

1.5

1

0.5 Imaginary Part 0

-0.5

-1

-1.5

-1

-0.5

0

0.5

1 Real Part

32

1.5

2

2.5

3



ceros = 0 1.3333 polos = 3.0000 0.5000

4

x 10 5

Respuesta Impulsional segun impz

3 h[n] 2 1 0 0 4

2

x 10 5

4

6 8 n Respuesta Impulsional analitica

10

12

10

12

3 hh[n] 2 1 0 0

2

4

6 n

33

8



6. Transf ormada d de F Fourier d de T Tiempo Discreto ((DTFT) A partir de la Transformada Z

 X  z  z e

j

X e

Im{z}

j

z  e j ,

      ω

1 Re{z} 

    x[n]  e    j

j n

X e

n

Esta es la DTFT directa

Existencia d de lla D D T FT Condición suficiente: Si la secuencia es absolutamente sumable, la DTFT existe. Es decir, si

, la DTFT existe.

x n n

La D DTFT iinversa x n

1 2 



X e j

e j n d 



Ejemplo 14 Determinar la DTFT de la secuencia impulso unitario x n Solución: X e j

 n e n

jn

1 X e j 

 [n

1

1

0

Ejemplo 15

 n .

0

n

[]  ( )[]

 Graficar la magnitud y la fase de la DTFT de la secuencia  Solución La secuencia es absolutamente sumable, por tanto, la DTFT existe,

( )          34





close all; clear all; w=[-2*pi:2*pi/200:2*pi]; j=sqrt(-1); X=1./(1-(1/2)*exp((-1)*j*w)); subplot(211); plot(w/pi,abs(X));grid; title('Magnitud'); subplot(212); plot(w/pi,180*angle(X)/pi);grid; ylabel('°');title('Fase'); Magnitud 2 1.8 1.6 1.4 1.2 1 0.8 -2

-1.5

-1

-0.5

0

0.5

1

1.5

2

0.5

1

1.5

2

Fase 30 20 10 °

0 -10 -20 -30 -2

-1.5

-1

-0.5

0

Interpretación d de lla D DTFT

|X ( Z )| Im{ Z }

e j |H (e j )|

Re{ Z }

-

35

0

 



R espuesta e en F Frecuencia

  j

Es la DTFT de la respuesta impulsional : H e

Magnitud de la Respuesta en Frecuencia :

H e j  (Respuesta en amplitud o

ganancia del sistema)

 H e j



Fase de la Respuesta en Frecuencia : desplazamiento de fase del sistema)

 (Respuesta

en fase o

j    H  z  z e j H e A partir de la Función de Transferencia :

Signif icado d de lla R R espuesta e en A Amplitud y y lla R R espuesta e en F Fase. Ejemplo 16 Dado el sistema de tiempo discreto causal con función de transferencia: 0.0042  z

 0.0147  z -2  0.0034  z -3 H  z   -1 -2 -3 -4 1 - 3.1582  z  4.0820  z - 2.5235  z  0.6279  z -1

Obtener experimentalmente el gráfico de la respuesta en amplitud y la respuesta en fase del sistema.

Solución

x[n]

y[n]

H ( z )

%************************************************ % Obtención Experimental de la Respuesta % en amplitud y la Respuesta en fase de un % Sistema de Tiempo Discreto %************************************************ % Reinicializar el ambiente close all; clear all; 36



% Función de Transferencia B = [0.0000 0.0042 0.0147 0.0034 0.0000] A = [1.0000 -3.1582 4.0820 -2.5235 0.6279] % Respuesta impulsional [h,n]=impz(B,A); stem(n,h); grid; ylabel('h[n]'); xlabel('n'); pause; % Obtención Experimental de la Respuesta en magnitud % y la Respuesta en Fase %N=length(n); N=1000; n=1:N-1; W=[0*pi:pi/800:pi]; % Vector de Frecuencias Nw=length(W) for m = 1:Nw x=sin(W(m)*n); % Secuencia de entrada y=filter(B,A,x); % Obtencion de la salida x=x(N-900:N-1); y=y(N-900:N-1); % Salida estable ganancia(m)= max(y)/max(x); xh = hilbert(x); yh = hilbert(y); fase(m)=wrapToPi( angle(yh(800)) - angle(xh(800)) ); end % Con freqz H=freqz(B,A,W); Hdecibelios=20*log10(abs(H)); % Graficos subplot(2,2,1); plot(W/pi,20*log10(ganancia)); grid; xlabel('w/pi'); ylabel('|H(e^jw)|,db'); title('Respuesta en amplitud obtenida experimentalmente'); subplot(2,2,3); plot(W/pi,Hdecibelios); grid; xlabel('w/pi'); ylabel('|H(e^jw)|,db'); title('Con la funcion freqz'); subplot(2,2,2); plot(W/pi,180*fase/pi); grid; xlabel('w/pi'); ylabel('


0.25

0.2

0.15

0.1 h[n] 0.05

0

-0.05

-0.1

0

50

100

150

n Respuesta en amplitud obtenida experimentalmente

Respuesta en fase obtenida experimentalmente

20

200

0 100 w)|,db -20

w),° j
j |H(e

0

-40 -100 -60

-80

-200 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

w/pi

0.5

0.6

0.7

0.8

0.9

1

0.7

0.8

0.9

1

w/pi

Con la funcion freqz

Con la funcion freqz

20

200

0 100 w)|,db -20

w),° j
j |H(e

0

-40 -100 -60

-80

-200 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

w/pi

0.5 w/pi

Ejemplo 17 Para el sistema de promedio móvil de tamaño N 1) Determine analíticamente la respuesta en frecuencia. 2) Grafique la magnitud y la fase para N=5. 3) Compruebe los resultados anteriores usando MATLAB. Solución

h[n]

0.4

1/ N

0

N

38

n

0.6




    hne j

H e


N 1

 j n

n

 e n0

1 N

 j n

N N j  j   2  2  e e 2 e      N 1     j  j  j  e 2  e 2 e 2     

1e 1

 N

 j N

1e

R

 

r

 j

r 0

N

 j

  j

H e

5 sen 1 2 5

sen 2

e

 N sen 1 2

 N

sen 2



e

 j  N 1  2

 j 2

%****************************************** % Respuesta Frecuencial del % Sistema Promedio Móvil N=5 %****************************************** % Reinicializar el ambiente close all; clear all; % 1) Función H hallada analiticamente N=5; w=[(-1)*pi:(2*pi)/800:pi]; j=sqrt(-1); H=(1/N).*((sin(w.*(N/2)))./(sin(w./2)))... .*exp(((-1)*j).*(w./2).*(N-1)); % 2)Gráficos de la solución analítica % Graficar la magnitud de H(e^jw) subplot(2,2,1); plot(w/pi,abs(H)); grid; title('Magnitud calculada analiticamente'); ylabel('|H(e^jw)|'); xlabel('w/pi'); % Graficar la fase de H(e^jw) subplot(2,2,2); plot(w/pi,angle(H)); grid; title('Fase calculada analíticamente'); ylabel('< H(e^jw)'); xlabel('w/pi'); % 2)Comprobacion con freqz B=[1/5 1/5 1/5 1/5 1/5]; A=[1 0 0 0 0]; H=freqz(B,A,w); % Graficar la magnitud de H(e^jw) subplot(2,2,3); plot(w/pi,abs(H)); grid; title('Magnitud (Comprobacion con freqz)'); 39

R1



1

1



ylabel('|H(e^jw)|'); xlabel('w/pi'); % Graficar la fase de H(e^jw) subplot(2,2,4); plot(w/pi,angle(H)); grid; title('Fase (Comprobacion con freqz)'); ylabel('< H(e^jw)'); xlabel('w/pi'); Magnitud calculada analiticamente

Fase calculada analíticamente

1

3 2

0.8

1 w) j < H(e 0

w)| 0.6 j |H(e 0.4

-1 0.2

-2

0 -1

-0.8

-0.6

-0.4

-0.2

0 w/pi

0.2

0.4

0.6

0.8

-3 -1

1

-0.8

-0.6

Magnitud (Comprobacion con freqz)

-0.4

-0.2

0 w/pi

0.2

0.4

0.6

0.8

1

0.4

0.6

0.8

1

Fase (Comprobacion con freqz)

1

3 2

0.8

1 w) j < H(e 0

w)| 0.6 j |H(e 0.4

-1 0.2

-2

0 -1

-0.8

-0.6

-0.4

-0.2

0 w/pi

0.2

0.4

0.6

0.8

-3 -1

1

-0.8

-0.6

-0.4

-0.2

0 w/pi

0.2

Filtros c clásicos 

Filtro pasa bajas

j

|H (e )|

- 

Filtro pasa altas

0





0





j

|H (e )|

-

40

Curso: Procesamiento Digital de Señales 


Filtro pasa banda

|H (e j )|

- 

Filtro supresor de banda

0





0





|H (e j )|

- Filtros IIdeales C Clásicos 

Filtro pasa bajas j

H (e )

[]   - 

Filtro pasa altas

-c

0

c





0

c





j

H (e )

[]  []   -

-c

Ejemplo 18 1) Generar una secuencia senoidal x[n] de 80 muestras con frecuencia 2/25 rad/muestra. 2) Obtener una secuencia x1[n] adicionando a la secuencia inicial un ruido aleatorio uniformemente distribuido en el intervalo -0,25 a 0,25. 3) Obtener la secuencia y[n] filtrando la secuencia x1[n] con un sistema promedio móvil de tamaño 5 recursivamente usando la función filter. 4) Lo mismo del punto 3 usando la función fft. 41


x [n]


X

Tamaño=84

Y h [n]

y [n] Tamaño=84

H

Tamaño=84

%***************************************** % Filtrado con el Sistema de Promedio % Movil de Tamaño 5 %***************************************** % Reinicializar el ambiente close all; clear all; % 1) Formar la secuencia de entrada N=80; n=[0:N-1]; x=sin(2*pi*n/25); % 2) Formar la secuencia con ruido aditivo x1=x+0.5*(rand(1,N)-0.5); % 3) Calcular la salida usando la función filter B=[1/5 1/5 1/5 1/5 1/5]; A=[1 0 0 0 0 ]; y1=filter(B,A,x1); % 4) Usando la función fft % El tamaño de h, x1 debe ser N+5-1=N+4 h=[1/5 1/5 1/5 1/5 1/5 zeros(1,N-1)]; x2=[x1 zeros(1,4)]; H=fft(h); X2=fft(x2); Y2=H.*X2; y2=ifft(Y2); y2=real(y2); % Graficar las secuencias subplot(4,1,1); stem(n,x); grid; ylabel('x[n]'); xlabel('n'); subplot(4,1,2); stem(n,x1); grid; ylabel('x1[n]'); xlabel('n'); subplot(4,1,3); stem(n,y1); grid; ylabel('y1[n]'); xlabel('n'); title('Salida obtenida empleando la funcion filter'); subplot(4,1,4); stem(n,y2(1:N)); grid; ylabel('y2[n]'); xlabel('n'); title('Salida obtenida empleando la funcion fft'); 42



1 x[n]

0 -1

0

10

20

30

0

10

20

30

0

10

20

30

0

10

20

30

40 n

50

60

70

80

40 n Salida obtenida empleando la funcion filter

50

60

70

80

40 n Salida obtenida empleando la funcion fft

50

60

70

80

50

60

70

80

2 x1[n] 0 -2

2 y1[n]

0

-2

2 y2[n]

0

-2

40 n

DTFT c con f f unciones iimpulso Sea una DTFT representada en la forma: X e j



1

- 0

-

X e

j

x n

   0 1  = X e j 2  

= = = =

1 2  1 2 

1 2  1

e

 0 n

j 0 n

 0

0





   0

e j n d  =



 

 0

 

 0

e

j n

d  =

 

 

 0

e

j n

d 



e e

j 0 n

1 2 

j 0 n

2  1 = c os  0 n  cos

1

e

e

j 0 n

1 2 



 

 0

e j n d  =



=

j 0 n

  

 0

2     0

  

 0

con periodo 2  43

con periodo 2 



7.Análisis e en e el D Dominio T Transf ormado El r retardador iideal La respuesta al impulso del retardador ideal es:

hid n   n  nd 

 ( )    ( )  ( )

 Su respuesta en frecuencia es   cuya magnitud y la fase son:   ;   Se observa que la fase es función lineal de la frecuencia. Filtros c clásicos iideales c con r retardo

El filtro pasa-bajas ideal con retardo tiene respuesta en frecuencia:

()   ||     ) []  (    

y respuesta al impulso

R etardo d de g grupo Es una medida de la linealidad de la fase:

    {( )}

Significado del retardo de grupo xn  sn cos  0n

j

; yn  ? y n

x n

H (e

 

X e

j

)

 S e     H e   e      H e   e     S e     H e   e 

 12 S e j 

 

 12 S e j 



  0

  0

j

1 2

j

j H e j



1 2

0

j H e j



j

  0

0

-0

-

   

Y e j  X e j  H e j

j

 H (e

j

 

j H e j

j

  0



)

-0 -

-0

0

0 0

44











Aproximamos  H e j  en  0 :  H e j   0    0    Aproximamos  H e j en   0 :  H e j    0    0   



   S e     H e   e   

Y e j

1 2

j   0

j 0

   S e     e 

Y e j

j   0

yn  sn   0  e

j  0 

j

j 0

e

 

1 2

j 0

yn  H e

j 0

 12 S e j



 S e j 



  0

j0

1 2

H e



 

  0 

 

H e

j 0  n  0 

0

  0

j0

e

 H e   e    j 0

 e

j

j  0 

 0

j 0

n  0 

 

  0 

 

H e

 sn   0  e j

0

1 2

e

 j0

 

H e

j 0

1 2

e

 j0

 sn     cos n         0

0

0

0

0

Interpretación g gráf ica d de lla r respuesta e en f f recuencia

        ( )            

A

Filtro R R esonador D Digital Filtro pasabanda de banda estrecha. Se puede implementar con dos polos conjugados complejos situados cerca de la circunferencia de radio uno. Filtro resonador con ceros en el origen. H  z  

Para b0

que

1  r  e se

b0

j 0

b0



 z 1  1  r  e j  z 1  , 

0

cumpla

 

H e j

   1

H e

j 0

 1  r  1  r 2  2r  cos2 0  .

45

el

1  r  e

j 0

 e j  1  r  e j  e j  0



coeficiente

b0





debe

ser



 1  r 2  La frecuencia de resonancia es  r  arccos  cos 0  . El ancho de banda de 3 db  2r  es   21  r  .

Ejemplo 19 Hallar la función de transferencia de un resonador digital con ganancia pico unitaria a 50 Hz, un ancho de banda de 3 db de 6 Hz, sabiendo que la frecuencia de muestreo es 300 Hz. clf; clear; f0=50; deltaf=6; fs=300; w0=f0*2*pi/fs deltaw=deltaf*2*pi/fs r=1-0.5*deltaw b0=(1-r)*sqrt(1+r*r-2*r*cos(2*w0)) B=[b0 0 0] j=sqrt(-1); A=conv([1 (-1)*r*exp(j*w0)],[1 (-1)*r*exp((-1)*j*w0)]) W=[0:0.5:150]*2*pi/fs; subplot(1,2,1); H=freqz(B,A,W); plot(W*fs/(2*pi),abs(H));grid; subplot(1,2,2); W=([45:0.1:55]*2*pi)/fs; H=freqz(B,A,W); plot(W*fs/(2*pi),abs(H));grid; 1.4

1.2

1.2

1.1

1

1

0.8

0.9

0.6

0.8

0.4

0.7

0.2

0.6

0 0

50

100

0.545 46 47 48 49 50 51 52 53 54 55

150

Filtro resonador con ceros en z  1 y z  1. j 2 1 1 1  z  1  z  1  e    j H e   b0 H  z   b0 1  r  e j  z 1  1  r  e j  z 1  , 1  r  e j  e j  1  r  e j  e j  0

0

0

0

Ejercicio 1. Para ambos tipos de resonador, graficar con MATLAB la magnitud y la fase en el intervalo de   a  ; para  0   3 , r  0.8 y para  0   3 , r  0.95 .

46



Filtro R R anura Filtro con uno o más cortes profundos en su respuesta en frecuencia. Para implementar un filtro ranura se pueden tomar dos ceros complejos conjugados sobre la circunferencia de radio uno. j j  z 1  H  z   b0  1  e  z 1   1  e Para reducir el ancho de banda del corte se pueden introducir dos polos conjugados complejos j  j 1  1  e  z  1  e  z 1  H  z   b0  1  r  e j  z 1  1  r  e j  z 1   0

 0

0

0

0

0

Ejemplo 20 Para el filtro ranura de dos ceros y dos polos, graficar con MATLAB la respuesta en amplitud en el intervalo de a  para  0  4 , r  0.85 . Solución





close all; clear all; w0=pi/4; r=0.85; j=sqrt(-1); b0=((1-r*exp(j*w0))*(1-r*exp((-1)*j*w0)))/... ((1-exp(j*w0))*(1-exp((-1)*j*w0))); B=b0*conv([1 (-1)*exp(j*w0)],[1 (-1)*exp((-1)*j*w0)]); A=conv([1 (-1)*r*exp(j*w0)],[1 (-1)*r*exp((-1)*j*w0)]); W=[0:pi/400:pi]; H=freqz(B,A,W); plot(W/pi,abs(H));grid; title('Respuesta en amplitud de un filtro ranura'); Respuesta en amplitud de un filtro ranura 1.4

1.2

1

0.8

0.6

0.4

0.2

0

0

0.1

0.2

0.3

0.4

0.5

47

0.6

0.7

0.8

0.9

1



Filtro P Peine Su respuesta en magnitud consiste en una serie de picos regularmente espaciados, cuya figura es semejante a la de un peine. El sistema promedio móvil es un ejemplo sencillo de filtro peine. yn  M 1 1

H  z   M 11

Si

M

 k  0

z  k , H  z   M 11

reemplazamos

  j

H L e

1 M 1

sen M 21  L sen 2 L 

Ejemplo 21

1  z 

1  z 

1

z e

 M 1

, H e j   M 11 L

por

z

sen M 21   sen 2 

e

 j M 2 

tenemos

M

 xn  k  k  0

. H L  z  

 L M 1

1 M 1

1  z

 L

1  z

,

 j M  L 2

    

Para el filtro peine con y graficar con MATLAB la respuesta en magnitud y en fase en el intervalo de a  . Solución

close all; clear all; % Respuesta en frecuencia L=10; M=3; w=[0:(2*pi)/800:pi]; j=sqrt(-1); H=(1/(M+1)).*((sin(w.*(L*((M+1)/2))))./(sin(w.*(L/2))))... .*exp(w.*((-1)*j*M*L*(1/2))); % Grafico de la magnitud de H(e^jw) subplot(2,1,1); plot(w/pi,abs(H)); grid; title('Respuesta en amplitud de un Filtro Peine'); ylabel('|H(e^jw)|'); xlabel('w/pi'); % Grafico de la fase de H(e^jw) subplot(2,1,2); plot(w/pi,angle(H)); grid; title('Respuesta en fase de un Filtro Peine'); ylabel('< H(e^jw)'); xlabel('w/pi');

48


Profesor: Wildor Ferrel Serruto Respuesta en amplitud de un Filtro Peine

1 0.8

| )

0.6

w j e ( H |

0.4 0.2

0

0

0.1

0.2

0.3

0.4

0.5 w/pi

0.6

0.7

0.8

0.9

1

0.7

0.8

0.9

1

Respuesta en fase de un Filtro Peine 3 2 1 )

w j e ( H <

0 -1 -2 -3

0

0.1

0.2

0.3

0.4

0.5 w/pi

0.6

Ejercicio 2. Para el sistema con función de transferencia H  z  

1. 2. 3. 4.

16

1 4

1  z

4

1  z

Dibujar el diagrama de polos y ceros y determinar la ROC Hallar la respuesta impulsional Hallar y graficar la magnitud de la respuesta en frecuencia ¿Qué tipo de filtro es?

Sistema iinverso Para un sistema LTI su sistema inverso es aquel que al conectársele en cascada da por resultado un sistema identidad H  z   Gz   1 G z  

1

H  z 

Ver ejemplo 5.4 y ejemplo 5.5 de Oppenheim

Sistema p pasa ttodo Es el sistema cuya respuesta en frecuencia tiene magnitud constante igual a la unidad.

H e j   1 Ejemplo, H ap  z  

  j

H ap e

z 1  a 1

1  a  z

e  j  a 1 a  e

 j

=e

 j



1 a  e 1 a  e

j

 j

En un sistema pasa-todo cada polo está emparejado con un cero conjugado recíproco. z 1  a 1 Si a es real: H ap  z   1 . Cero: z  . Polo: z  a 1  a  z a 49


Profesor: Wildor Ferrel Serruto Im{z}

Re{z} 1

Sistema d de f f ase m mí nima Un sistema estable y causal se llama sistema de fase mínima si su sistema inverso también es causal y estable. Los polos y los ceros de un sistema de fase mínima deben estar dentro de la circunferencia unitaria.

Descomposición e en p pasa ttodo y y f f ase m mí nima Un sistema estable y causal siempre se puede expresar como la conexión en cascada de un sistema de fase mínima y un sistema pasa todo.

H  z   H min  z   H ap  z  Ejemplo 22 H  z  

1

1  2  z

1  0.2 z  1  0.7 z  Para el sistema: Hallar la descomposición en pasa todo y fase mínima 1

1

Im{z}

Re{z} 0.2

0.7

1

2

Im{z}

0.5 0.2

50

0.7

Re{z}

1

2

Curso: Procesamiento Digital de Señales 1

H  z  

1  2  z

1  0.2 z   1  0.7 z   1

H ap  z  

1

z 1  a



Profesor: Wildor Ferrel Serruto 1

1  0.5  z

1

1

   1  0.2 z  1  0.7 z  =

1  0.5  z 1  0.5  z

1

1

1



1  2  z

1

1  0.5  z

;

1

1  a  z

1  0.5  z 

1 1  1  2  z     H  z     2 1  0.2 z 1  1  0.7 z 1   2  1  0.5  z 1 1

1

H  z     2

1  0.5  z

1  0.2 z   1  0.7 z   1

1



z 1  0.5 1

1  0.5  z

Sistema c con f f ase llineal Un sistema LTI se dice que tiene fase lineal si: H e j

j  j  = H e   e

Un sistema LTI es de fase lineal generalizada si:  j      j H e j = Ae  e j Ae  - Es una función real. Para que un sistema causal de función de transferencia racional tenga fase lineal generalizada su respuesta al impulso de tamaño N debe ser real, FIR, simétrica ( h n h N 1 n ) o anti-simétrica ( h n h N 1 n ). N impar N par h n simétrica Tipo I Tipo II h n anti-simétrica Tipo III Tipo IV

x [n]

x [n]

n

0

n

0

Tipo I

Tipo II

x [n]

x [n]

n

0

n

0

Tipo III

Tipo IV

Para el filtro de tipo II, tenemos N

H e

j

=

2

1

h n e n 0

jn

+

N 1 n

h n e jn

;

m N 1 n ; n N 1 m

N 2

51

Curso: Procesamiento Digital de Señales N

H e j

H e j

= =

2

N

1

jn

h n e

+

=

2

1

m 0

N

N

2

1

jn

h n e

2

+

m 0

jn

j N 1 n

e

2

=

e

j

N 1 2

n

N 1

j

e

2

n

N 1

e

j

2

n 0 N

H e

1

h n

n 0 j

h N 1 m

h m e j N 1 m N

e

; h m

1

1

h n

m e j N 1 m

h N 1

n 0

n 0 N

2


2

=

1

2 h n cos 

N 1

N 1

n

2

n 0

e

j

2

N 1

H e

j

= H r

e

j

e

j

2

N impar h n simétrica Tipo I h n anti-simétrica Tipo III Tipo I

H r e

N 1

h

2

N

II

2

 

j o H i e

j

j 0

H e

N 1

+

2

N par Tipo II Tipo IV

2 h n cos 

n 0

N 1 2

1

n

h

2 h n cos 

2

n 0

j

N 1 2

N 1

+  2hn

2

n0

N 2

N 1

 

H e

1

 hn

0

n

0

0

n

0

2

n

n0

N 1

III

2

2 h n sen

j



N 1



N 1

n 0 N

IV

2

2

1

2 h n sen

j n 0

2

Conclusiones:

Los filtros de tipo II no pueden ser pasa altas ni supresores de banda.  Los filtros de tipo III sólo pueden ser filtros pasa banda.  Los filtros de tipo IV no pueden ser pasa bajas ni supresores de banda.  Los filtros de tipo I puede implementar cualquier filtro. Es el único tipo capaz de realizar filtros supresores de banda. Los ceros de la función de transferencia de un sistema de fase lineal se presentan en pares conjugados recíprocos. 

52



8. R R elación e entre lla D DTFT y y lla C CTFT x[n]

xc(t )

1

C/D X c

j

X e

x[n]

xc(t )

j

C/D

T

2

3

(1) x n

x c

nT

X c

j

T

(2)

x c t x n

1

e jt d 

X c j

2  1

X c j

2 

e jnT d  ,

T , 



x n

x n

1

X c j

2 

 T

e j n

2 k 1 

1 2  k

2 k 1 

x n x n



2  k



1  2  

1 T

 1 X c j T T

k

T

 2  k T

2  k ,



x n



e j

n



d  2  La relación entre las transformadas es: 

X e j



1

X c j

T k

54



e j n d 



X e j

T

2  k , d 

e j  2  k n d 

(3) 1

d 

, d 

e j n d  ,



 2  k 1 X c j T T

X c j



d  T



1

X e j



2 k T

d 



Interpretación g gráf ica Para k

Para k

0:

1:

Para k 1 :

1

X 0 e j

T

X c j

1

X 1 e j

T 1

X 1 e j

T

 T

X c j

X c j



X c j 



;



T

,

 T



 2  T

 2  T



1

-S

0

-M X  e j

X 1 e j



S

M





1

X 1 e j

X 0 e j

T …

…

-MT 0

-2

2

MT



X S  j 

1

T …

…

-S

-M

0

S

M



X S  j  X e j T  X c j 

T X e j  T 0

55



T

 c.c .



T



Estructura d del c conversor D D/C iideal X c j

X e jT H r j Conversor de secuencia en tren de impulsos ponderados

x[n]

x s(t )

Filtro de reconstrucción

H r  j 

T

H r  j  T







T

T



Relación entre la señal reconstruida y la secuencia

x n e jTn

X c j

H r j

n

x c t

x n  t nT

* h r t

n

hr t  

sen T t   T

t

xc t  



 xn

n  

56

sen  T t  n   T t  n 

xc(t )



Conversor D D/C p práctico x[n]

xa(t )

Retenedor de orden cero

X a j

X e j T

x[n]

xc(t ) xa(t )

0

0

n p (t )

t

1

0 x a t

x n

T

t

p t nT

n

x n P j e jnT

X a j n

X a j

P j

x n

e jnT

n

X a j

P j X e jT

P j

p t

P j

T

e jt dt

e jt dt

0

P j

T e jt j 0

P j

1

P j

1

j 1 j

e jT

e

T j 2e

1 T j 2

e

T T j j 2e 2

sen P j

57

2

T 2



e

j

T 2



P j X e jT

X a j

sen2T   j T 2 X a  j   2  e  X e jT 



sen2T   j T 2  X e jT X a  j   T T  e





2

X e

jT



T 1 T

sen2T 

X c  j   X e

 e j  X a  j  T 2

2

jT

 H  j  r

T 1 T

2

sen

T 2



e

j T 2

 X a  j  H r  j

 2T j   e  X a  j  T    T  X c  j   sen T  2  0 c.c. T 2





 2T j  e  H rr  j    sen2T   0

T 2

x[n]

Retenedor de orden cero

xa(t )

  T     T c.c.

Filtro de reconstrucción H rr  j  

X c j

X a j

X e j T

T

H rr  j

58

xc(t )



Ejemplo 23 Graficar la respuesta en amplitud del filtro de reconstrucción, conectado en cascada con el retenedor de orden cero, para la frecuencia de muestreo de 8 kHz. T T

 2 j  e  H rr  j    sen2T   0

2

  T     T c.c.

% Reinicializar el ambiente close all; clear all; % Hallamos la función Hrr fs=8000; T=1/fs; W=[(-1)*2*pi/T:(2*pi/T)/400:2*pi/T]; N=length(W); j=sqrt(-1); for i=1:1:N if abs(W(i))
-6000

-4000

-2000

0 f

59

2000

4000

6000

8000



9. C Cambio d de lla T Tasa d de M Muestreo R educción d de lla f f recuencia d de m muestreo e en u un f f actor e entero Sistema compresor x n

x c

x d n

x c nT

x d n

x nM

X e j



Si M

1

X c j

1 MT r

kM , 0

X d e j

X d e

x c nMT

T k

X d e j

r i

nT

j

X c



2 k T

j

 2 r MT

i M 1 ,

k

1 M 1 1 M i 0 T k j 1 M 1 X e M i 0

2 , X d e

j

1 2

X e j

X c



j

 2 i 2 kM MT

2  i M

X e

j

 2

X e

j

 2  2

Ver Figura 4.21 y Figura 4.22 de Oppenheim

Sistema diezmador Filtro pasabajas Ganancia = 1 Corte = /M

60

X d e j



Incremento d de lla f f recuencia d de m muestreo e en u un f f actor e entero x n

x c

nT

x i n

x c nT

x c n

T L

Sistema expansor

x k  n

x e n

k L

k

X e e j

x k n

X e e j

k L

 n

k L

e j n

k

x k k

X e e j

 n

e j n

n

x k e j k L k

X e e j

X e j L

Ver Figura 4.25 de Oppenheim

Sistema interpolador Filtro pasabajas Ganancia = L Corte = /L

61



Cambio d de lla ttasa d de m muestreo e en u un f f actor n no e entero Filtro pasa-bajas Ganancia=L Corte= min(π/L, π/M)

Procesamiento d digital d de sseñales a analógicas Sistema de tiempo discreto

H e j X c

j

X e j

Y e j

Y r j

T Y e j T ,

Y r j 

T H e jT X e j T ,

Y r j 

T H e j T

Y r j

H e j T X c j ,

 T

1 T



X c j

H eff j 

Y r j

 T

T T  T

 T



 T

,

 T





 T

H e j T 0

62

 T

 



T ,  T



Ejemplo 24 Dado el sistema

Sistema de tiempo H e j X c

Y e j

X e j

j

donde,



X c j 

Y r j



1

-2 1000 H e j 

0

2 1000 1 0





c

c.c.

¿Cuál es la frecuencia mínima para que al muestrear la entrada no se produzca solapamiento? 

Si  c y r t

2

, ¿Cuál es la frecuencia mínima de muestreo que asegura que

x c t

Para  c

Y e j ,

 2

? y una frecuencia de muestreo igual a 2000 Hz, graficar

Y r j

, y hallar

H eff j 

63

X e j ,



Aplicación d del d diezmado e en lla c conversión A A/D

Sistema de tiempo discreto

Filtro antisolapamiento





X c j 

1

- N

 N

0



H aa j 

1

-M1

- N

0

 N

M1

Filtro antisolapamiento exacto Corte=/M

Filtro antisolapamiento simple

64









X c j 

1

-4 N

0

- N

 N

4 N



 N

4 N



0

 N

4 N



0

 N T

2



2



H aa j 

1

-4 N

0

- N

X a j 

1

-4 N

- N X e j  ˆ

1

T

-2

 M

H e j

1

0

-2

X d e j 1

T

-2

-

T

0

 65

MT

2





10.

Cuantización

Cuantización por redondeo: Unipolar ∆ - Tamaño de paso de cuantización

Bipolar:

,- ,- ,66





f e 

    



0

Potencia del ruido (varianza):

    ∫   

e

  √                          Error RMS de cuantización:

;

;

D - Intervalo dinámico L – Número de niveles

Razón señal ruido de cuantización:

 La SQNR mejora en 6 dB por cada bit adicional.  La SQNR se reduce si el intervalo dinámico excede la rms de la señal.

Ejemplo 25 Una señal varía en el intervalo -2.5V a +2.5V. ¿Cuál debe ser la cantidad de bits para que el error rms de cuantización sea menor que 5mV? D=5V;

  √   √    √   √      ;

; B=9 bits.

67


11.    


Estructuras d de S Sistemas d de T Tiempo Discreto

Forma directa I Forma directa II Forma serial Forma paralela Directa I, II

 

Serial

Paralela

Inmunidad a los efectos de aritmética de precisión finita Facilidad de cálculo

Símbolos

x 2[n]

x 1[n]

a

x[n]

x 1[n]+ x2[n]

a a x [n]

-1

x[n]

z

x [n-1]

Ejemplo 26 Para el sistema con función de transferencia

dar la representación en la forma directa I, forma directa II, forma serial y forma paralela. Solución

y[n] = b0 x[n]+ b1 x[n-1]+ b2 x[n-2] - a1 y[n-1] – a2 y[n-2] y[n] = 0.5 x[n]-0.5 x[n-2] – 1.3 y[n-1] – 0.36 y[n-2] 68



Forma directa I

y[n] = 0.5 x[n]-0.5 x[n-2] – 1.3 y[n-1] – 0.36 y[n-2]

y[n]

x[n]

x[n]

y[n]

0.5 -1

-1

z

z

-1.3

0

-1

-1

z

-0.36

-0.5

z

Forma directa II

x[n]

y[n]

0.5 -1

z

-1

z

-1.3

0 -1

-0.36

z

-1

z

w[n]

x[n]

-0.5

0.5 -1

z -1.3

0 -1

-0.36

z -0.5

   xn  1.3wn  1  0.36wn  2

wn

   0.5wn  0.5wn  2

y n

69

y[n]



Forma serial

w1[n]

x[n]

y1[n]

0.5

w2[n] -

-1

z -0.4

1

z -0.9

-0.5

1

    y n  0.5w n  0.5w n  1 w n  y n  0.9w n  1 yn  w n  w n  1 w1 n  x n  0.4 w1 n  1 1

1

2

1

1

2

2

2

Forma paralela

-1.39

w2[n]

x[n]

2.1

y2[n]

-1

z -0.4 w3[n]

-0.21 -

z -0.9

   1.39 xn w n  xn  0.4w n  1 y n  2.1w n y1 n 2

2

2

2

   xn  0.9w n  1

w3 n

3

  yn  y n   y n   y n  y3 n  0.21w3 n 1

2

3

70

y3[n]

y[n]

y[n]



Forma d directa II y n

b 0 x n

b1 x n 1

a1 y n 1

...

... b M x n M

a N y n N

x[n]

,

H z

b0

b1 z 1

...

b M z M

1

a1 z 1

...

a N z N

y[n]

b0 -1

z -1

z -a1

b1

-1

-1

z

z

-a N-1

bM-1

-1

-1

z

z

-a N

bM

v n

b 0 x n

b1 x n

y n

a1 y n

V z

b0

1

1

...

...

b M x n

a N y n

... b M z M

b1 z 1 1

Y z 1 a1 z

1

N

...

a N z

N

M

v n

X z

V z

Forma d directa III x[n]

y[n]

b0 -1

z -a1

b1 -1

z

-a N-1

bM-1 z -1

-a N w n

y n

a1 w n

b0 w n

1

...

b1 w n

a N w n

1

...

1

W z 1 a1 z

1

...

a N z

N

N

bM x n

b M w n

M

Y z

X z

71

b0

b1 z 1

... b M z M

W z



Forma sserial M 1

1 g k z H z

b0

k 1 N 1 k 1

N s

1 b1 k z 1

b 2 k z 2

k 1

1 a1 k z 1

a 2 k z 2

b0

M 2

* 1 h k z 1 1 h k z 1

k 1

1 c k z

H z

1

1

N 2 k 1

* 1 1 d k z 1 1 d k k z

b0

y 0 n w k n

y k n y n

x n a 1 k w k n

w k n b 0 y N

s

1

a 2 k w k n

b1 k w k n 1

2

y k 1 n b 2 k w k n 2

n

72



Forma p paralela M N

H z

C k z

k

k 0

N 1

N 2

N 1

A k

N 2

B k

N 2

* B k

k 11

c k z 1

k 11

d k z 1

k 11

* d k z 1

N

M N

H z

C k z

k

k 0

N s k 1

b 0 k 1

b1 k z 1

a1 k z 1

a 2 k z 2

b01

b11

b02

b12

b03

b13

w k n y k n

a1 k w k n b 0 k w k n M N

y n

C k x n k 0

1

a 2 k w k n

2

x n

b1 k w k n 1 N s

y k n k 1

73



R epresentación de sistemas de tiempo discreto en el espacio de estados Ejemplo 27 Para el sistema con función de transferencia,

Hallar la representación en el espacio de estados en la forma:

,,--  ,, --  ,,--

Solución Las ecuaciones de la forma directa II son:

,-,-  ,-  ,,- ,- ,,- ,- ,- ,- ,,-  ,-  ,-  ,,- ,- ,- ,,      [,,]0 1[ ]0 -   ,- 1,,-, -[,],,    0  1 01   , -  ;

;

En matlab escribimos: b=[0.5 0 -0.5] a=[1 1.3 0.36] [A,B,C,D]=tf2ss(b,a)

74


12.


La S Serie d de F Fourier D Discreta

Una secuencia periódica x~ n con periodo N se puede expresar a través de la serie de Fourier discreta (DFS): ~ n x

1 N

2

N 1

j kn ~ X k e N

k 0 ~ n x

La DFS es la descomposición de en una suma de N exponenciales complejas armónicamente relacionadas. Los valores de los coeficientes de la DFS se calculan con la fórmula: ~ X k

N 1

~ n x

e

j

2 N

kn

n 0

Ejemplo 28 ~ Hallar los coeficientes de la DFS de la secuencia periódica x n esta secuencia como una serie de Fourier Discreta

y expresar

1

... Solución

... 0

3

6

 ,-  ∑̃,-  

n

̃ ,-  ̃,-  ̃,-  ,- ,,- -  ,-      ,-  ̃,-   ∑ ,- ̃,-   ,- ,-   ,-  /  ̃,-    . /    . ̃,-    .   /  ./          ̃,-     () ̃,-    .  /   ./           ̃,-     () /    . /  ̃,-     .      ̃,-             ̃,-               75



̃,-    ()  () La D DFS y y lla D DTFT

Para la secuencia periódica x~ n con periodo N los coeficientes de la DFS son: ~ X k

N 1

~ n x

j

e

2 N

kn

n 0

Definimos la secuencia de extensión finita x~ n , 0 n N 1 x n 0, c .c . La DTFT de la secuencia

X e j = X e j =

x n

es:

N 1

x n e jn =

x n e jn

n N 1

n 0

~ n e jn x

n 0

Se observa que se cumple

~ X k

X e j



Ejemplo 29

2 N

k

1

Hallar la DTFT de la secuencia

0 j X e 

x n e jn

=

=

2

n

1 e j 0  1 e j 1  1  e j

n

Si comparamos esta DTFT con los coeficientes de la DFS del ejemplo anterior vemos que se cumple la relación:  j 2 k ~ X k   1  e 3  X e j 





 23 k





La c convolución p periódica ~ Dadas las secuencias periódicas x , x~ 2 n ambas con periodo N, se define 1 n la convolución periódica como la operación:

~ n x 1

x~ 2 n

N 1

x~1 m

x~ 2 n

m

m 0

Ejemplo 30 Hallar la convolución periódica de las secuencias x~1 m y x~ 2 m ~ n x 3

N 1

~ m x 1

m 0

~ n x 2

m

76



~ m x 1

4

4

3

…

4

3

2

3

2

2

1

1

…

m

0 ~ m x 2

…

…

1

m

0 ~ x 2

…

m

1

…

0

m ~

x 2

…

1

m

…

1

m

0 ~ 0 x 3 ~ x 3 1 ~ 2 x 3

~ 3 x 3 ~ x 3 4

7 5

3 6 9

Ver Tabla 8.1 de Oppenheim

77


13.


La T Transf ormada D Discreta d de F Fourier (DFT)

Dada una secuencia de extensión finita de longitud N x n

n

0 Para encontrar

x n

N -1

a partir de su DTFT, deben ser suficientes N muestras de

X e j en las frecuencias  0 ,  1 , ...,  N 1 . Puesto que la DTFT está dada por la expresión X e j = x n e jn n

tenemos N ecuaciones con N incógnitas

X e

j 0

= x 0

j 0

x 1 e

x 2 e

j 0 2

... x N 1 e

j 0 N 1

… j j N 1 j N 1 2 j N 1 N 1 ... x N 1 e X e N 1 = x 0 x 1 e x 2 e Este sistema de ecuaciones tiene solución única si sólo si las frecuencias  0 ,  1 , ...,  N 1 son distintas. 2 j k , k 0 ,1,..., N 1 , las N muestras de X e Si escogemos  se llaman N

Transformada Discreta de Fourier (DFT) N 1

x n

X k

e

j

2 kn N

, 0

k N 1

n 0

c .c .

0,

y su transformada inversa es: 1

x n

N

N 1

X k e

j

2

N

kn

, 0

n N 1

k 0

c .c .

0,

Usando el factor Twiddle

   N 1

X k

n 0

kn , 0 x n W N

k N 1 c .c .

0,

y su transformada inversa es: 78



N 1

1

x n

N

k 0

X k W N kn , 0 n N 1 c .c .

0,

Operación n módulo N n

Re siduo de n N Re siduo de

N

,

n

N ,

N

n

0

n

0

La DFT y la DFS Si definimos x~ n

x n

r N , entonces

X k

~

X k

, 0 k N 1 .

r

Esto significa que X k corresponde a un periodo de la DFS de una secuencia periódica x~ n obtenida por la repetición de x n . También, podemos escribir que, si x~ n

x

, entonces

n N

La DFT y la DTFT X e j

2 N



X k

k

,

0

0,

k N 1 c .c .

La DFT y la Transformada Z X z

X k

j

2

z e N

k

, 0

0,

c.c .

Desplazamiento C Circular x n

n

0 1 2 3 4 x

n

k N 1

3 5

n

0 1 2 3 4 79

~

X k

X

k N



La c convolución c circular Dadas las secuencias x1 n , circular como la operación: x 1 n

N x 2 n

de extensión finita, se define la convolución

x 2 n N 1

x 1 m

m 0

x 2

n

, 0 n N 1

m N

Ejemplo

Hallar la convolución circular de las secuencias

x1 n

, 5

.

x 2 m 4

5

x 1 m

x 2 n

3 2

1

1

m

0 5

x 2

m

0

m 5 4 3 2

1

m

0 x 2

1

m 5

5 4

x1 m

3

5 x 2 m

2 12

1 10

m

0

9 8

x 2

2

m 5

6

5 4 3 2 1

0

0

m

80

m



La c convolución c circular y y lla c convolución llineal Si

x1 n

o

x 2 n

x1

tiene extensión M N , entonces n N x 2 n x1 n x 2 n para M 1 n N 1 .

Ejemplo

La convolución lineal de las secuencias del ejemplo anterior es: x1 n

x 2 n

12 9 9

6 5 3 1

n

0 Si

x1 n

tiene extensión N 1 y x1 n

tiene extensión N 2 , entonces x1 n x 2 n si N N 1 N 2 1 .

x 2 n

N x 2 n

Ejemplo

La convolución circular de longitud 7 de las secuencias convolución lineal. 5

x 2 n

es igual a la

x 2 n 3 2

1

0

,

4

7

x1 n

x1 n

1

n

0

Ver Tabla 8.2 de Oppenheim

81

n


14.


Comple jidad d del c cálculo d de lla ssalida

1. Por convolución x n 

y n



h n 

y n

x m h n

m

m

Sean las secuencias x

n

,

h n

Operaciones  +

de extensión finita de longitud M

1+2+...+M-1+M+M-1+...+2+1 0+1+...+M-2+M-1+M-2+...+1+0

Total M (M-1)

2. Mediante la DFT de longitud 2M-1

x[n]

X [k ]

h[n]



H [k ] M 1

X k

x n

e

Y [k ]

j

2

N

kn

,

k

0 ,1,..., N 1

n 0

N 2 M 1

Por cada muestra de

X k

Operaciones  +

Total M M-1

Operaciones  +

Total M (2M-1) (M-1) (2M-1)

Para una DFT

Para el producto

X k H k

Operaciones  +

Total 2M-1 0

82

y[n]



Para la IDFT y n

Por cada muestra de

2 M 2

1 2 M 1

Y k e

j

2 N

kn

,

0

n

2 M

2

k 0

y n

Operaciones  +

Total 2M-1 2M-2

Operaciones  +

Total 2 (2M-1) (2M-2) (2M-1)

Para toda la IDFT

Total mediante DFT Operaciones  +

2M(2M-1) + (2M-1) + (2M-1) 2(M-1) (2M-1) + (2M-2) (2M-1)

Total Aprox. 2  8M  8(M-1)2

La T Transf ormada R R ápida d de F Fourier ((FFT) Mediante diezmado en el tiempo j

N 1

X k

x n

e

2

N

k n

,0

n 0 j

X k

x n e

2

N

2

1

X k

x 2 r e r 0

N 1

k n

n par

N

k

j

x n e

2

N

k n

n impar

j

2

N

N

k 2 r

2

1

j

x 2 r 1 e r 0

83

2

N

k 2 r 1


N 2

j

1

X k

2 N

k r j

2

x 2 r e

e

2

N

k

2

1

x 2 r 1 e

j

2 N

k r

2

r 0

r 0       

DFT de muestras pares

DFT de muestras impares

G

k

N

e

2

N

k

H

k

2

G

k

N 2

N 2

Usando la notación W N

X k

N

    

j

X k


e

j

2 N

k W N H

k

N 2

DF T de N/2 puntos

DFDFT T de de N/2 puntos N/2 puntos

84



Número de operaciones en cada etapa Operaciones Total N  + N Número de etapas:

log 2 N

Total Operaciones por FFT: Operaciones  + Estructura Mariposa

Total N log 2 N N log 2 N

85


r

W N

N 2

e

2

j

N

r

N 2


r W N

Número de operaciones: Operaciones  +

Total (N/2) log 2 N log 2 N

Complejidad del cálculo de la salida Operaciones 2 M  3

+

2

N

Total log log 2 2 M

M 5

2 M

3 2 M log 2 2 M

6 M 1

3 log 2 M log 2 M

Ejemplo

Para N=256 Operaciones  +

Cálculo por convolución M = 65 65536 536 (M-1)2 = 65025

Cálculo con FFT M 5 3 log 2 M = 7424 6 M 1 log 2 M = 13824

Cálculo d de lla ssalida m mediante lla D D FT

x[ x[n] h[n]

X [k ] H [k ]



1 86

Y [k ]

y[ y[n]

2

3



Filtrado F FIR p y ssuma c con D DFT por ssolapamiento y L

L

x1 n

L

1 P 1 –

ceros x 2 n

P 1 1 –

ceros x 3 n

y 1 n

P 1 1 puntos –

que se suman y 2 n

1 puntos P 1 –

que se suman y 3 n

87



Filtrado F FIR p y a almacenamiento c con D DFT por ssolapamiento y

L

L

L

x1 n

P-1

P 1 1 –

ceros x 2 n

x 3 n

y 1 n

1 P 1 –

puntos que se desprecian

y 2 n

1 P 1 –


y 3 n P 1 1 –


88


15.


Diseño d de F Filtros D Digitales

El objetivo es determinar los coeficientes de la ecuación en diferencias de forma que el filtro cumpla un esquema de tolerancias

y n  b0 x n  b1 x n  1...b M xn  M 

 a1 y n  1... .. . a N y n  N  Al diseñar un filtro pasa bajas se puede exigir que la curva de la magnitud de la respuesta en frecuencia cumpla las especificaciones dadas en el esquema j j 

|H (e H (e )| 1+ 1+  p 1 1- 1 p r 0

H e j 

 p

r





 p

r





dB

0

- p

-r

      ⁄      ⁄  

       (  )⁄ 89



Métodos d de D Diseño d de F Filtros Filtros IIR A partir de filtros analógicos  Invarianza impulsional  Transformación bilineal  Aproximación de derivadas Método directo  Minimización del error cuadrático medio Filtros FIR  Por enventanado  Por optimización  Por muestreo de la respuesta frecuencial

Diseño d de u un f f iltro d digital a a p partir d de u un f f iltro a analógico

Esquema de Tolerancias del Filtro Digital

Paso 1 Esquema de Tolerancias del Filtro Analógico

Paso 2 Función de Transferencia del Filtro Analógico H a( s)

Paso 3 Función de Transferencia del Filtro Digital H ( z )

Paso 4 Ecuación en Diferencias

90



Filtro d de B Butterworth

|H (j)|, dB 0



r



- p

-r

  10   10  1  log   10  1    N    r   2  log   p    r

 C 

 p 10

 p 1

   2 N 10 10  1     p

sk  C

 sk 



  N  1  k  2

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

|H a( j)| 1

H a  s  



1 2

N C

N

  s  s

k



k 1

0 91

 C





Filtro d de C Chebyshev

|H (j)|, dB  0

r



    10 10  1  1 cosh    10  1    N    1   r  cosh   p    r

 p 10

- p

-r  p

  10 10

1

   1  1    2 1  N

1

1

a  12   N  

1  N

b  12   N  

sk  a   p  cos  k  j  b   p  sen k  k 

   N  1  k  2 2 N

; k  0,1,..., N  1

K :

 1 , N par  2 1    H a 0    1, N impar  

|H a( j)| 1 1- p

0

H a  s  

K N

  s  s

k



k 1

92

 p





Cálculo del Coeficiente K de la Función de Transferencia

K :

|H a( j)|

 1 , N par  H a 0    1   2 1, N impar

1 1- p

0 H a  s  

 p

 1 , N par  H a  j    0   1   2 1, N impar



K N

  s  s

k



 1 , N par  H a 0   N   1   2   sk  1, N impar K

k 1

k 1 N  1    sk , N par  2 1   k 1 K   N    sk , N impar  k 1

Invarianza IImpulsional Se basa en el muestreo de la respuesta impulsional del filtro analógico

h n

T h a nT

Paso 1 Al pasar de las especificaciones del filtro digital a las especificaciones del filtro analógico basta hacer : 



T

Relación entre los polos del filtro analógico y los polos del filtro digital N

H a s

ha t  

A k

k 1 s N

 k 1

s k

Ak  e sk t  u t  93



N

h n

T h a nT

T

A k e

s k nT

u nT

k 1

N

T

A k e

s k T n

u n

k 1 N

H z

T

A k

k 1 1 e s k T z 1

Los polos del filtro digital son :

p k

Si

e

s k T

e

T Re S k

Re sk   0 

e

jT Im S k

pk  1

En consecuencia, si el filtro analógico es estable, entonces el filtro digital será también estable. Como Ha( j) nunca es perfectamente limitada en frecuencia, se produce superposición espectral. Por ello, este método no se aplica a filtros pasa altas ni supresores de banda.

Paso 3

H a s

B a s

H z

A a s Residue()

H a s

h a t

N

H z

A k e

A k

T

k 1 s s k

s k t

A z

Residuez()

N A k

N

B z

k 1 1 e

N

h n

u t

T

A k k 1

k 1

94

e

s k T

s k T n

z 1

u n



Ejemplo 31 A partir de un filtro analógico de Butterworth diseñar un filtro digital IIR pasa bajas, por Invarianza Impulsional, con las especificaciones :

H e j  ,

dB

0,2  0,3 





0 -1

-15

 p

1

 r  15

H a  j 

 p

 p 

0.2

 r 

T

0.3

 r

T



0

  p

  r H a  s  

 N C

N  5.89

N

  s  s 

N  6

k

k 1

s1, 6 

0.703 T

 C  e j 180 75 

s 2 ,5  95

0.703

0.703 T

T e j 180  45 


s3, 4 

0.703 T


e j 18015  6

 0.703    T   H a  s    s  s1  ...   s  s6 

H a  s  

A1

 s  s  1

 ... 

A6

 s  s  6

A1, 6  0.144  j 0.249  T A2 ,5   1.07 T A3, 4  0.927  j1.61 T H z

T

A1 1 z 1 z 1

s1, 6T

 0.649  j 0.524

s 2 , 5T

 0.535  j 0.290

z 1, 6  e z 2 ,5  e

s3 , 4T

z 3, 4  e

...

A 6 1 z 6 z 1

 0.499  j 0.092 %**************************************************** % Diseño de un filtro digital a partir de un filtro % analogico de Butterworth por Invarianza Impulsional %**************************************************** % Reinicializar el ambiente clear;clf; ap=1; ar=15; wp=0.2*pi; wr=0.3*pi; % PASO 1 T=1; Wp=wp/T; Wr=wr/T; % PASO 2 % Hallamos la orden del filtro N=(log10(((10^(ar/10))-1)/((10^(ap/10))-1)))... /(2*log10(Wr/Wp)); N=ceil(N); N=N % Hallamos la frecuencia de media potencia Wc=(Wp)/(((10^(ap/10))-1)^(1/(2*N))); % Hallamos los polos del filtro analogico theta= (pi*(N+1)+2*pi*[0:N-1])./(2*N); j=sqrt(-1); sk=Wc*cos(theta)+j*Wc*sin(theta); % Funcion de transferencia del filtro analogico 96



B=[ Wc^(N)]; A=poly(sk); % PASO 3 % Expandimos en fracciones parciales [R,P,K] = residue(B,A); % Hallamos los polos del filtro digital pk=exp(P*T); [b,a] = residuez(R,pk,K); % PASO 4 b=real(b) a=real(a) % Graficar la magnitud de la respuesta en frecuencia W=[0:pi/50:0.35*pi]; H=freqz(b,a,W); Hdecibelios=20*log10(abs(H)); plot(W/pi,Hdecibelios); grid; title('Respuesta en frecuencia'); ylabel('|H(e^jw)|, db'); xlabel('w/pi'); Respuesta en frecuencia

0

-5 b d , | )-10 w j ^ e ( H |

-15

-20

-25

0

0.05

0.1

0.15

w/pi

0.2

0.25

0.3

0.35

N = 6 b = 0.0000 0.0006 0.0101 0.0161 0.0041 0.0001 0 a = 1.0000 -3.3635 5.0684 -4.2759 2.1066 -0.5706 0.0661 yn = 0.0000*x n + 0.0006*x n-1 + 0.0101*x n-2 + + 0.0161*xn-3 + 0.0041*xn-4 + 0.0001*x n-5 + + 0*xn-6 + + 3.3635*yn-1 - 5.0684*yn-2 + 4.2759*y n-3 + - 2.1066*yn-4 + 0.5706*yn-5 - 0.0661*y n-6

97



Transf ormación B Bilineal Para obtener la función de transferencia del filtro digital a partir de la función de transferencia del filtro analógico se aplica la transformación

H  z   H a  s  s  2 1 z 1 T

1 z 1

Correspondencia entre el plano S y el plano Z

1 

2 T

Plano S

Plano Z

El eje imaginario del plano S es mapeado en la circunferencia unitaria del plano Z :  j

z  e j

s   2 T

1 e

1  e  j

 j T 2 tan 2 

j  j T 2 tan  2  Paso1

2



T

2

  tan

|H (e j )|

|H a( j)|

0

0

 98



2





La transformación bilineal inversa es :

1

T s 2

1

T s 2

z

A cada polo y a cada cero de Ha( s) aplicamos esta transformación y encontramos los polos y los ceros de H( z) N

  z  z  k

H  z   A 

k 1 N

  z  p  k

k 1

El coeficiente A se calcula de la condición :

H e j 

  0

 H a  j 

0

Cálculo del coeficiente A al usar el Filtro Analógico de Chebyshev de tipo I

H e j 

  0

 H a  j 

 H 1  H a 0

0

 1 , N par  2 H 1  H a 0   1   1, N impar N

 1  z 

 1 , N par  k 1 2   H 1  A  N   1    1  pk  1, N impar k

k 1

N  1  pk    1  , N par  k N 1 2  1   1  z k    k 1 A   N   1  pk   k 1 , N impar  N   1  z k   k 1

99



Ejemplo 32 A partir de un filtro analógico de Chebyshev tipo I diseñar un filtro digital IIR pasa bajas, por Transformación Bilineal, con las especificaciones :

H e j  ,

dB

0,2  0,3 





0 -1

-15

 p  1

 0.2  0.650  p  tan   T  2   0.3  1.02 2  r  T tan   T  2  2 T

 r  15

H a  j 

 p  r



0

  p

  r

100




H a  s  


N  3.02

K N

  s  s

k

N  4



k 1

s1, 4  s2,3 

K 

0.044 T 4

 0.219  j 0.265

0.044

T

4 T H a  s    s  s1    s  s2    s  s3    s  s4 

 0.091  j 0.640 T

H z

H a s

s s k

2

1 2 1 z

s

T

1 z

z 1

T z

2 T

T s 2 k

z 1

2

s k

1

1

s

1

T

2

1 z 1

2

T

1 z 1

T z

z 1 1

T 2

s k z 1

z 1 T s 2 k

z 1 z 1 0.044 T 4

H z 16

4

1

T 2

s k z 1

T 4 k 1 0.044

H z

T 4 16

s 2 k

z 1 z 1 4

T 4 k 1

1

z 1

z 1

T 2

s k z 1

Los polos del filtro digital son :

Los ceros del filtro digital son :

T

p k

T

s 2 k

1

T

1

T s 2 k

2

s k

z 1 , 2 , 3 , 4

101

z 1

1

T s 2

1

T s 2

z 1



%****************************************** % Diseno de un filtro digital % a partir del filtro de Chebyshev I % por transformacion bilineal %****************************************** % Reinicializar el ambiente clear;clf; ap=1; ar=15; wp=0.2*pi; wr=0.3*pi; % PASO 1 Wp=2*tan(wp/2); Wr=2*tan(wr/2); % PASO 2 % Hallamos el orden del filtro N=acosh(sqrt((10^(ar/10)-1)... / (10^(ap/10)-1))) / acosh(Wr/Wp); N=ceil(N); N=N % Hallamos los polos del filtro analogico epsilon=sqrt(10^(ap/10)-1); j=sqrt(-1); angulos=( pi*(N+1) + [0:N-1]*2*pi )./(2*N); alfa=(epsilon^(-1)) + sqrt(1+epsilon^(-2)); aa=(1/2)*( (alfa^(1/N)) - (alfa^(-1/N)) ); bb=(1/2)*( (alfa^(1/N)) + (alfa^(-1/N)) ); sk=aa*Wp*cos(angulos)+j*bb*Wp*sin(angulos); % Hallamos la constante multiplicativa k preliminar if rem(N,2)==1 K=1; else K=1/sqrt(1+epsilon^2); end % PASO 3 % Hallamos los polos y zeros del filtro digital pk=(1+sk/2)./(1-sk/2); zk=-1*ones(1,N); A=K*prod(ones(1,N)-pk)/prod(ones(1,N)-zk); % PASO 4 % Coeficientes del filtro digital b=A*poly(zk); a=poly(pk); a=real(a) b=real(b) 102



% GRAFICO W=[0:pi/50:0.45*pi]; H=freqz(b,a,W); Hdecibelios=20*log10(abs(H)); plot(W/pi,Hdecibelios); grid; title('Respuesta en frecuencia'); ylabel('|H(e^jw)|, db'); xlabel('w/pi'); Respuesta en frecuencia 0 -5 -10 -15

b d , | ) -20 w j ^ e ( -25 H |

-30 -35 -40 -45 0

N b a y

0.05

0.1

0.15

0.2

0.25 /pi

0.3

= 4 = 0.0018 0.0073 0.0110 0.0073 = 1.0000 -3.0543 3.8290 -2.2925 n = 0.0018*xn+0.0073*xn-1+0.0110*x n-2+ +0.0073*xn-3+0.0018*xn-4+ +3.0543*yn-1-3.8290*yn-2+ +2.2925*yn-3-0.5507*yn-4

103

0.35

0.0018 0.5507

0.4

0.45



Método d de E Enventanado Un filtro pasa bajas ideal tiene respuesta en frecuencia:

H i (e j ) 1

-

-C

La respuesta impulsional es :

0

C





hi [n] 1/4

hi n  

sen C n 

...

...

n n

0

hi [n]

1/4

...

... 0 w [n]

n

1

Ventana Rectangular

hn  hi n wn

h [n]

0

n

0

n

1/4

104



Operación en el Dominio Temporal

hn  hi n wn

hi n  

sen  C n  n

1,  N 21  n  N 21 wn   0, c.c.

Ventana Rectangular

Respuesta al impulso real

 sen C n   ,  N 21  n  N 21 hn   n  0, c.c.

Operación en el Dominio Frecuencial j j 

j j 

H i (e )

W (e W (e )

1

-

-C

0

C





0

-









j j 

H (e H (e )

H e j   H i e j  W e j 

W e

j



sen  N 2  sen  2  -

0

DTFT de la Ventana  |W (e |W (e j j )|

A -

0 

 

105



 - Ancho del lóbulo principal

A - Amplitud máxima de los lóbulos secundarios, relativa al lóbulo principal En la ventana rectangular tenemos : A= -13 dB   4 N

Aproximación de la Respuesta Frecuencial Ideal j j

j j

|H (e )|

|H i (e )|

1+ 1-

C

-

 W (e

j(( -  ) j

)

Cuando aumenta el tamaño de la ventana se reduce la banda de transición; sin embargo, la amplitud del risado no disminuye (efecto de Gibbs).  = -21 dB  - Error máximo entre H(e j ) y Hi (e j ) para cualquier valor de N es

Ejemplos de ventanas

•

w [n]

Ventana Rectangular

1

wn   1, 0  n  N  1

0

•

N n

w [n]

Ventana de Bartlett o Triangular

1

 2 N n11 , wn   2  n 1 2 ,  N 1 

0  n  N 21 N 1 2

 n  N  1 106

0

N n



w [n] 1

•

Ventana de Hanning

0

wn  12  12 cos

2 n 1 N 1

N

n

N

n

, 0  n  N  1 w [n] 1

•

Ventana de Hamming

0

wn  0,54  0,46 cos

  , 0  n  N  1 2 n N 1

w [n]

•

1 Ventana de Blackman

wn   0,42  0,5 cos



2  n 1 N 1

  0,08 cos

4  n 1 N 1

,

0

N

n

w [n]

•

1

Ventana de Kaiser

0 wn  

2 I 0    1  n      



   N  1 2,

I 0   

 , 0  n  N  1

I 0 . es función de Bessel

107

=4

N

n



Tabla comparativa

Ventana

A, dB

, dB

de Kaiser

p/ Kaiser

Rectangular

-13

4/N

-21

0

1,81/N

Bartlett

-25

8/N

-25

1,33

2,37/N

Hanning

-31

8/N

-44

3,86

5,01/N

Hamming

-41

8/N

-53

4,86

6,27/N

Blackman

-57

12/N

-74

7,04

9,19/N

Procedimiento de diseño por enventanado buscando utilizar la ventana más simple 1. Dado el esquema de tolerancias, se especifica la magnitud de la respuesta en frecuencia ideal j

|H i (e )| 1+ p 1- p r 0 2. 3. 4. 5. 6.

 p

r





Hallar hi [n] Escoger una ventana con error máximo inferior a las tolerancias Calcular N y redondear hacia arriba Obtener el filtro h [n] = hi [n] wi [n] Verificar si las tolerancias son satisfechas: o Si no satisface, aumentar N e ir al paso 5 hasta que satisfaga o Si satisface, reducir N e ir al paso 5 mientras satisface

108



Ejemplo 33 Diseñar un filtro pasa-bajas que cumpla las especificaciones j

|H i (e )| 1 + 0,02 1 - 0,02 0,01 0

0,2

0,3





1. Trazamos la magnitud de la respuesta en frecuencia ideal j

|H i (e  )| 1 + 0,02 1 - 0,02 0,01 0 2. Hallamos hi [n] :

0,2

0,3

 C   p   r 2



hi n 



sen0.25 n  n

3. Escogemos una ventana con error máximo inferior a las tolerancias :

min  r ,  p  0,01 20 log 0,01  40dB Hanning, Hamming o Blackman 4. Calculamos N:

  8 N    r   p ;

N 

8  r   p

N  80

%******************************************** % Diseño de un filtro digital por Enventanado %******************************************** % Reinicializar el ambiente clear;clf; % PASO 1 deltap=0.02; deltar=0.01; wp=0.2*pi; wr=0.3*pi; % PASO 2 wc=(wp+wr)/2; % PASO 3 delta=min(deltar,deltap); deltadB=20*log10(delta) 109



% PASO 4 N=8*pi/(wr-wp); N=ceil(N); if rem(N,2)==1 %Para filtros LP no es necesario N=N else N=N+1 end %N=N-20 % PASO 5 % Formamos la Ventana de Hanning w = hanning(N); w=w'; n=[0:N-1]; stem(n,w); grid; title('Ventana de Hanning'); ylabel('w[n]'); xlabel('n'); % Respuesta Impulsional Ideal pause; hi=(sin(wc*(n-(N-1)/2)))./((n-(N-1)/2)*pi); hi(((N-1)/2)+1)=wc/pi; stem(n,hi); grid; title('Respuesta Impulsional Ideal'); ylabel('hi[n]'); xlabel('n'); % Respuesta Impulsional del Filtro pause; h=hi.*w; stem(n,h);grid; title('Respuesta Impulsional del Filtro'); ylabel('h[n]'); xlabel('n'); % Formamos los coeficientes de la Ecuación b=h; a=zeros(1,N); a(1)=1; % Graficamos la magnitud de la % respuesta en frecuencia W=[0:pi/400:0.4*pi]; H=freqz(b,a,W); Hdecibelios=20*log10(abs(H)); pause; plot(W/pi,Hdecibelios); grid; title('Respuesta en frecuencia'); ylabel('|H(e^jw)|, db'); xlabel('w/pi');

110



Ventana de Hanning

1 0.9 0.8 0.7 0.6

] n [ w 0.5

0.4 0.3 0.2 0.1 0

0

10

20

30 n Respuesta Impulsional Ideal

40

50

60

0.25 0.2

0.15

] n [ i h

0.1 0.05 0

-0.05

-0.1

0

10

20

30 n

40

50

60

40

50

60

Respuesta Impulsional del Filtro 0.25

0.2

0.15 ] n [ h

0.1

0.05

0

-0.05

0

10

20

30 n

111



Respuesta en frecuencia 20

0

-20 b d , | ) w j ^ e ( H |

-40

-60

-80

-100

0

0.05

0.1

0.15

0.2 w/pi

0.25

0.3

0.35

0.4

b = 0.0000 -0.0001 0.0000 0.0003 0.0008 0.0008 0.0000 -0.0015 -0.0028 -0.0025 0.0000 0.0039 0.0066 0.0056 0.0000 -0.0079 -0.0131 -0.0108 0.0000 0.0147

0.0243 0.0202 0.0000 -0.0283 -0.0483 -0.0422 0.0000 0.0733 0.1575 0.2245 0.2500 0.2245 0.1575 0.0733 0.0000 -0.0422 -0.0483 -0.0283 0.0000 0.0202

0.0243 0.0147 0.0000 -0.0108 -0.0131 -0.0079 0.0000 0.0056 0.0066 0.0039 0.0000 -0.0025 -0.0028 -0.0015 0.0000 0.0008 0.0008 0.0003 0.0000 -0.0001

0.0000

Dif erencias e entre F Filtros F FIR e IIIR    

Los filtros IIR producen en general distorsión de fase, es decir la fase no es lineal con la frecuencia. Los filtros FIR son de fase lineal. El orden de un filtro IIR es mucho menor que el de un filtro FIR para una misma aplicación. Los filtros FIR son siempre estables.

112



Diseño m mediante lla v ventana d de K K aiser wn  

2 I 0    1  n      



I 0   

   N  1 2,

 , 0  n  N  1

I 0 . es función de Bessel

Para calcular el parámetro de forma β hacemos: – Parámetro de ondulación. - Atenuación real en la banda de supresión.

(  )               {                %****************************************************** % Diseño de un filtro pasa-bajas empleando la ventana de Kaiser

%****************************************************** close all; clear all; %PASO 1 wp=0.2*pi wr=0.3*pi deltap=0.02 deltar=0.01 delta=min(deltap,deltar) %PASO2 wc=(wp+wr)/2 %PASO3 deltaw=wr-wp A=(-1)*20*log10(delta) if A<21 beta=0 elseif A<=50 beta=0.5842*((A-21)^(0.4))+0.07886*(A-21) else beta=0.1102*(A-8.7) end N=((A-8)/(2.285*deltaw))+1 N=ceil(N) if rem(N,2)==1 N=N else N=N+1 end alfa=(N-1)/2 %PASO 4 wkaiser=kaiser(N,beta) w=wkaiser' 113



n=0:N-1; hi=(sin(wc*(n-alfa)))./(pi*(n-alfa)) hi(((N-1)/2)+1)=wc/pi; x=length(hi) y=length(w) h=hi.*w % Formamos los coeficientes de la Ecuación b=h; a=zeros(1,N); a(1)=1; % Graficamos la magnitud de la % respuesta en frecuencia W=[0:pi/400:0.4*pi]; H=freqz(b,a,W); Hdecibelios=20*log10(abs(H)); plot(W/pi,Hdecibelios); grid; title('Respuesta en frecuencia'); ylabel('|H(e^jw)|, db'); xlabel('w/pi'); Respuesta en frecuencia 10

0

-10

-20

b d , | )

-30

w j

e ( H | -40

-50

-60

-70

-80 0

0.05

0.1

0.15

0.2 w/pi

0.25

0.3

0.35

0.4

%****************************************************** %Diseño de un filtro pasa-altas empleando la ventana de Kaiser

%****************************************************** close all; clear all; %PASO 1 wp=0.5*pi wr=0.35*pi deltap=0.021 deltar=0.021 delta=min(deltap,deltar) %PASO2 wc=(wp+wr)/2 %PASO3 deltaw=wp-wr A=(-1)*20*log10(delta) if A<21 beta=0 114



elseif A<=50 beta=0.5842*((A-21)^(0.4))+0.07886*(A-21) else beta=0.1102*(A-8.7) end N=((A-8)/(2.285*deltaw))+1 N=ceil(N) if rem(N,2)==1 N=N else N=N+1 end alfa=(N-1)/2 %PASO 4 wkaiser=kaiser(N,beta) w=wkaiser' n=0:N-1; hi=(sin(wc*(n-alfa)))./(pi*(n-alfa)) %hi(((N-1)/2)+1)=wc/pi; hi=(-1)*hi hi(((N-1)/2)+1)=1-(wc/pi) h=hi.*w % Formamos los coeficientes de la Ecuación b=h; a=zeros(1,N); a(1)=1; % Graficamos la magnitud de la % respuesta en frecuencia W=[0.2*pi:pi/400:pi]; H=freqz(b,a,W); Hdecibelios=20*log10(abs(H)); plot(W/pi,Hdecibelios); grid; title('Respuesta en frecuencia'); ylabel('|H(e^jw)|, db'); xlabel('w/pi');

Respuesta en frecuencia 10

0

-10

-20

-30 b d , | ) w -40 j e ( H |

-50

-60

-70

-80

-90 0.2

0.3

0.4

0.5

0.6 w/pi

115

0.7

0.8

0.9

1



Método d de m muestreo e en f f recuencia

 , -

Se obtiene muestreando la respuesta en frecuencia deseada en puntos equiespaciados en el dominio de la frecuencia.

Ejemplo 34

 

Diseñar un filtro pasa-bajas con frecuencia de corte Solución Tomamos 10 muestras ( ) dentro de un periodo de la respuesta en frecuencia ideal:

 |,-|  *+ ,-    , ,        ,-    ,-,- ,- ,-  ,-  ,- ,-,-,-,-,- Para que la respuesta al impulso simétrica con respecto a .

sea real,

debe ser conjugada

, ,

Ahora calculamos la DFT inversa: clear;clf; j=sqrt(-1); H0=1; H1= exp((-1)*j*0.9*pi); H2= exp((-1)*j*1.8*pi); H3=0; H4=0; H5=0; H6=0; H7=0; H8=conj(H2); H9=conj(H1); H=[H0 H1 H2 H3 H4 H5 H6 H7 H8 H9]; h=ifft(H) W=[0:2*pi/100:2*pi]; H=freqz(h,1,W); plot(W/pi,abs(H));grid; 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0

0.2

0.4

0.6

0.8

1

116

1.2

1.4

1.6

1.8

2



El resultado es:

h = [0.0716 -0.0794 -0.1000 0.1558 0.4520 ... 0.4520 0.1558 -0.1000 -0.0794 0.0716] Cerca a la frecuencia de corte hay una sobre-elongación; para reducirla podemos hacer

,- 

clear;clf; j=sqrt(-1); H0=1; H1= exp((-1)*j*0.9*pi); H2= 0.5*exp((-1)*j*1.8*pi); H3=0; H4=0; H5=0; H6=0; H7=0; H8=conj(H2); H9=conj(H1); H=[H0 H1 H2 H3 H4 H5 H6 H7 H8 H9]; h=ifft(H) W=[0:2*pi/100:2*pi]; H=freqz(h,1,W); plot(W/pi,abs(H));grid; 1.4 1.2 1 0.8 0.6 0.4 0.2 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

El resultado es:

h = [-0.0093 -0.0485 0.0000 0.1867 0.3711 ... 0.3711 0.1867 0 -0.0485 -0.0093]

117

2



%Diseño de filtro pasabajas N=20; wc=0.8*pi; gana1=zeros(1,(N/2)+1); for i=1:(N/2) if 2*(i-1)*pi/N
118



Filtro A Adaptativo Cancelación de ruido con filtro adaptativo

Señal y Ruido

,-,-,-

ADC

+ +

-

Ruido

ADC

 , -

 , Filtro Adaptativo

,-,-

Algoritmo LMS

  ,-,%---------------------------------% Filtrado Adaptativo % Algoritmo LMS %---------------------------------close all; clear all; [s, fs] = wavread('arch_voz.wav'); s=s'; N=length(s) t=0:1:N-1; t=t/fs; x=0.8*(rand(1,N)-0.5); n=filter([0 0 0 0 0 0.5],1,x); d=s+n; mu=0.01; Nw=31; w=zeros(1,Nw); y=zeros(1,N); e=y; for m=Nw+1:1:N-1 sum=0; for i=1:1:Nw sum=sum+w(i)*x(m-i+1); end 119

 , -

DAC



y(m)=sum; e(m)=d(m)-y(m); for i=1:1:Nw w(i)=w(i)+2*mu*e(m)*x(m-i+1); end end subplot(2,1,1); plot(t,s);grid;title('Señal original'); subplot(2,1,2); plot(t,d);grid;title('Señal original + Ruido'); pause; subplot(2,1,1); plot(t,s);grid;title('Señal original'); subplot(2,1,2); plot(t,e);grid;title('Señal restaurada'); wavwrite(d,'arch_voz_ruid.wav'); wavwrite(e,'arch_voz_rest.wav'); Señal original 0.15 0.1 0.05 0 -0.05 -0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.5

0.6

0.7

0.5

0.6

0.7

Señal original + Ruido 0.3 0.2 0.1 0 -0.1 -0.2 -0.3 -0.4

0

0.1

0.2

0.3

0.4

Señal restaurada 0.2 0.1 0 -0.1 -0.2

0

0.1

0.2

0.3

0.4

120



ANEXO Tablas de Transformadas tomadas de OPPENHEIM A. V., SCHAFER R.W. Tratamiento de Señales en Tiempo Discreto.

121

PROCESAMIENTO DIGITAL DE SEÑALES - TABLAS DE TRANSFORMADAS

Transformada Z

121


122


DTFT

123


124


125


DFS

126

Apuntes de Procesamiento Digital de Señales

Recommend Documents