• Author / Uploaded
• Jorge Ríos

Citation preview

Departamento de Electrónica Facultad de Ciencias Exactas, Ingeniería y Agrimensura Universidad Nacional de Rosario

Teoría de Sistemas y Señales

Trabajo Práctico Nº 3 Análisis Frecuencial de Señales Autor: Dr. Juan Carlos Gómez Ing. Silvina Ferradal Ing. Franco Del Colle Dr. Lucas Terissi Dr. Ariel Bayá Srta. Marianela Parodi Dra. Mónica Larese

Mayo de 2011

Trabajo Práctico Nº 3 Análisis Frecuencial de Señales

1. Objetivos El objetivo de este Trabajo Práctico es estudiar el uso de la Transformada de Fourier en Tiempo Discreto (DTFT) para el análisis espectral de señales, a la vez que la implementación en Matlab de algoritmos de cómputo de la DTFT, y el estudio de la performance de diferentes tipos de ventanas empleadas por los algoritmos.

2. Introducción En este Trabajo Práctico se estudiarán las propiedades fundamentales de la Transformada de Fourier en Tiempo Discreto en su uso para el análisis espectral de señales. Para algunas clases de señales (nominalmente aquellas que son de longitud finita, o aquellas que tienen una transformada Z del tipo racional) es posible implementar algoritmos que calculen exactamente la DTFT. Para estas señales se requerirá la programación de funciones Matlab para el cómputo de la DTFT. Esto se hace en la Subsección 3.1.1.. Para el caso de señales de longitud infinita, o que no tengan una Transformada Z racional , no es posible implementar algoritmos para el cómputo exacto de la DTFT. En estos casos se debe recurrir al uso de funciones ventana que trunquen la señal a un número finito de muestras, y permitan la implementación de algoritmos que aproximen la DTFT de la señal. El estudio de varios tipos de ventanas utilizadas en la práctica se realiza en la Subsección 3.1.2.. En la Subsección 3.1.3., se analizan las características de los algoritmos Matlab para el cálculo de la Transformada de Fourier en Tiempo Discreto (genéricamente denominados Transformada Discreta de Fourier (DFT)), en su versión conocida como Transformada Rápida de Fourier (FFT: Fast Fourier Transform). En esta misma Subsección se propone un ejemplo de aplicación para la detección de una señal senoidal inmersa en ruido. Finalmente, en la Subsección 3.2 se presenta el enunciado de los problemas que deberán ser resueltos y entregados en el informe.

3. Desarrollo del Trabajo Práctico 3.1. Problemas Propuestos de Simulación durante la Sesión de Laboratorio 3.1.1. Transformada de Fourier en Tiempo Discreto (DTFT) El objetivo de esta primera subsección es estudiar las propiedades de la Transformada de Fourier en Tiempo Discreto (DTFT: Discrete-Time Fourier Transform), a la vez que desarrollar e implementar algoritmos en Matlab [4] para el cálculo exacto de la misma. Consideraremos dos casos diferentes: señales de longitud finita, para las cuales la DTFT puede calcularse exactamente a partir de la definición, y señales de longitud infinita pero del tipo exponencial, que poseen una Transformada Z racional, para las cuales puede calcularse una expresión analítica de la DTFT que permite también para este caso su cálculo exacto, que de lo contrario sería imposible [1]. Recordemos que para una señal en tiempo discreto x(n) , la DTFT viene dada por [2]:

X (ω ) =

∞

∑ x ( n)e

− jωn

,

−π ≤ ω ≤ π

(1)

n = −∞

que resulta una función a valores complejos de la variable ω , que es periódica, con período fundamental 2π . En el contexto de Matlab, el cómputo de la DTFT basándose en la definición (1) presenta básicamente dos inconvenientes: i. No es posible usar la expresión (1) para el cálculo de la DTFT cuando la señal es de longitud infinita, ya que Matlab sólo puede manipular vectores de longitud finita. Hay sin embargo una excepción a esto, que es el caso en que pueda obtenerse una expresión analítica de la DTFT, que pueda usarse para el cómputo. Un ejemplo de este caso son las señales exponenciales discretas de la forma a TeSyS - TP No. 3

n

μ (n ) , que tienen una DTFT racional en la variable e jω .

2

X (ω ) es una función de la variable continua ω (que puede asumir infinitos valores en el intervalo [ −π ,π ) ). Para el cálculo con Matlab, esto representa un problema de

ii. La DTFT

muestreo en frecuencia. Lo mejor que se puede hacer es realizar el cálculo en un número finito de frecuencias, lo suficientemente grande como para obtener una aproximación suave de la DTFT. Generalmente se toman frecuencias equiespaciadas en el intervalo [ −π ,π ) , con lo que la expresión (1) resulta:

X (ω k ) =

∞

∑ x ( n)e

n = −∞

− jω k n

∞

∑ x ( n )e

=

− j ( 2πk / N )n

k = 0,", N − 1

,

(2)

n = −∞

Sin embargo, esta última expresión todavía no es computable, salvo para el caso en que la señal sea de longitud finita. En este caso, asumiendo que la señal sea de longitud L , la expresión (2) resulta L −1

L −1

n =0

n =0

X (ω k ) = ∑ x(n)e − jω k n = ∑ x(n)e − j (2πk / N )n ,

k = 0,", N − 1

(3)

Es fácil ver que la expresión (3) corresponde a la definición de la Transformada Discreta de Fourier (DFT) con N puntos, que para el caso N = L puede ser calculada eficientemente usando alguno de los algoritmos de la Transformada Rápida de Fourier (FFT: Fast Fourier Transform). A continuación se proponen algunos problemas que requieren la programación de funciones Matlab, para el cálculo de la DTFT . Problema 1: a. Escriba una función Matlab que compute exactamente la DTFT de una señal x(n ) de longitud finita L . • Los argumentos de entrada de la función deben ser: i. x : vector de longitud L conteniendo las muestras de la señal x(n ) .

[

•

• •

ii. N: número de frecuencias equiespaciadas en el intervalo − π , π ) donde la DTFT es evaluada. Los argumentos de salida de la función deben ser: i. X: vector (complejo) conteniendo los valores de la DTFT X (ω ) . ii. ω : vector conteniendo las frecuencias (en rad/seg) donde la DTFT es evaluada. La función deberá así mismo graficar los espectros de amplitud y fase de la señal, y el espectro de densidad de energía. La función debe incluir un "help" que pueda invocarse con el comando >> help function_name

b. Considere la señal

⎧(0,9 )n x(n ) = ⎨ ⎩0

0 ≤ n ≤ 100 c.o.c.

y use la función desarrollada en el apartado a. para calcular la DTFT. c.

Trabajo previo: Verifique el resultado calculando una expresión analítica de la DTFT y comparando la gráfica correspondiente con la obtenida con la función Matlab. ‰

TeSyS - TP No. 3

3

Problema 2: En este problema consideraremos señales de longitud infinita, pero que poseen una DTFT del tipo racional en la variable e

jω

. Esto es equivalente a decir que la señal tenga una Transformada

Z

−1

del tipo racional en la variable z (o z ). En efecto, recordando la definición de la Transformada Z de una señal

X (z ) =

∞

∑ x(n)z

−n

,

(4)

n = −∞

puede verse que si la Región de Convergencia (RDC) de la Transformada Z incluye a la circunferencia de radio unitario, entonces la DTFT X (ω ) de la señal puede calcularse evaluando la transformada Z jω

sobre la circunferencia unitaria (es decir para z = e ). a. Escriba una función Matlab que compute exactamente la DTFT de una señal de longitud infinita, cuya transformada Z es una función racional de la variable z . • Los argumentos de entrada de la función deben ser i. num: vector conteniendo los coeficientes del polinomio numerador de la Transformada Z de la señal, en orden decreciente de potencias de z . ii. den: vector conteniendo los coeficientes del polinomio denominador de la Transformada Z de la señal, en orden decreciente de potencias de z . iii. N: número de frecuencias equiespaciadas en el intervalo − π , π ) , donde la Transformada de Fourier es evaluada.

[

• Los argumentos de salida deben ser i. X: ii. w:

vector conteniendo los valores de la DTFT. vector conteniendo las frecuencias donde la DTFT es evaluada.

• La función deberá así mismo graficar los espectros de amplitud y fase de la señal, y el espectro de densidad de energía, en función de la frecuencia ω . • La función debe incluir un "help" que pueda invocarse con el comando >> help function_name NB: Para la implementación de esta función puede hacer uso de la función de Matlab polyval(u1,u2) que permite evaluar el polinomio, cuyos coeficentes en orden decreciente de potencias se ingresan como primer argumento u1, en el segundo argumento u2. b. Considere la señal

⎛1⎞ x(n ) = ⎜ ⎟ ⎝2⎠

n−2

⎛7⎞ ⎝8⎠

n

μ (n − 2) + ⎜ ⎟ μ (n )

,

y use la función desarrollada en el apartado a. para calcular la DTFT. c.

Trabajo Previo: Verifique el resultado calculando una expresión analítica de la DTFT

d. Compare la gráfica correspondiente con la obtenida con la función Matlab. ‰ 3.1.2.

Cómputo de la DTFT usando ventanas

En esta Subsección estudiaremos el uso de ventanas para el cómputo de la DTFT de señales de longitud infinita, que no necesariamente tienen una transformada Z racional. Como se vió en la Subsección anterior, no es posible usar la expresión (1) para el cómputo de la DTFT de una señal de longitud infinita, excepto para el caso en que la señal tenga una transformada Z racional. La alternativa es entonces truncar la señal a un número finito de muestras usando lo que se denomina una función de ventana, y luego computar la DTFT de la señal truncada que aproximará la DTFT de la señal original, bajo ciertas condiciones. Numerosos tipos diferentes de ventanas han sido propuestos en la literatura para su uso en TeSyS - TP No. 3

4

análisis espectral de señales (ver [1],[2] para más detalles). El tipo más común es la denominada ventana rectangular que sólo trunca la señal a un número finito de muestras, sin afectar el valor de las muestras. En el dominio temporal, la ventana rectangular de longitud L viene dada por:

1≤ n ≤ L c.o.c.

⎧1 w(n ) = ⎨ ⎩0

(5)

Si x(n ) es la señal de longitud infinita, la señal truncada xˆ (n ) se puede obtener como

xˆ (n ) = x(n )w(n ) .

(6)

La DTFT, X (ω ) , de la señal original está entonces relacionada con la DTFT , Xˆ (ω ) , de la correspondiente señal truncada, a través de la convolución (periódica) en el dominio frecuencial

1 Xˆ (ω ) = 2π

π

∫ X (λ )W (ω − λ )dλ .

(7)

−π

Puede verse entonces que al aproximar el espectro X (ω ) , de la señal original, por el espectro, Xˆ (ω ) , de la correspondiente señal truncada, se comete un error que dependerá de la forma de la función de ventana w(n ) (o de su espectro W (ω ) ) a través de la relación (7). Se proponen a continuación algunos problemas donde se estudiarán las características espectrales de algunas de las funciones de ventana disponibles en Matlab [3] (nominalmente las funciones hamming, hanning y boxcar, que implementan las ventanas de Hamming, Hann, y rectangular

respectivamente), y se analizará el error introducido en la aproximación de X (ω ) por Xˆ (ω ) .

Problema 3: Las funciones Matlab hamming, hanning y boxcar, permiten implementar las funciones de ventana de Hamming, Hann, y rectangular respectivamente. Estas funciones están definidas en el dominio temporal como [3]: • Ventana Hamming (hamming)

⎧ ⎛ 2πn ⎞ ⎪0.54 − 0.46cos⎜ ⎟ w(n + 1) = ⎨ ⎝ L −1 ⎠ ⎪⎩0 •

c.o.c.

Ventana de Hann (hanning)

⎧1 ⎛ ⎛ 2πn ⎞ ⎞ ⎪ ⎜1 − cos⎜ ⎟ ⎟⎟ w(n ) = ⎨ 2 ⎜⎝ ⎝ L +1⎠⎠ ⎪0 ⎩ •

0 ≤ n ≤ L −1

1≤ n ≤ L c.o.c.

Ventana Rectangular (boxcar)

⎧1 w(n ) = ⎨ ⎩0 a.

1≤ n ≤ L c.o.c.

Trabajo Previo: Calcule la DTFT de la ventana rectangular en forma teórica y determine los cruces por ceros del espectro de amplitud.

b. Utilice la función Matlab implementada en el Problema 1 para computar el espectro de amplitud de la ventana rectangular. En una misma gráfica muestre superpuestos el espectro para tres longitudes L distintas, L = 20, 31, 81. Use una escala lineal para el espectro de amplitud.

TeSyS - TP No. 3

5

c.

Mida la amplitud y el ancho del lóbulo principal de la DTFT. Analice los puntos de cruces por cero y la amplitud de los lóbulos laterales. Además, verifique la ubicación de los puntos de cruces por cero a partir del cálculo analítico de la DTFT de un pulso rectangular. NOTA: Recuerde que puede usar el comando ginput para extraer valores de una gráfica en Matlab.

d. Considere la siguiente señal de longitud infinita

⎛π ⎞ x(n) = cos⎜ n ⎟ ⎝2 ⎠ Utilice una ventana rectangular de longitud L=45 para computar una aproximación del espectro de la señal y graficarlo. e.

Sea la señal

x(n) = cos(0.5π .n ) + cos(0.51π .n )

Use la ventana rectangular para computar y graficar una aproximación del espectro de la señal x(n). Determine la relación entre la resolución en frecuencia y la longitud L de la ventana. A partir de este resultado, calcule la longitud L mínima de modo que se distingan las componentes armónicas de la señal x(n). ‰ 3.1.3.

Transformada Discreta de Fourier (DFT)

En esta sección se estudiará el uso de la Transformada Discreta de Fourier en el análisis espectral de señales. Particularmente se considerará la implementación de la DFT en un algoritmo del tipo FFT (Fast Fourier Transform) disponible en Matlab (la función fft). Como un problema de aplicación se empleará la DFT para determinar la amplitud y frecuencia de una señal senoidal inmersa en ruido. Recordemos que para una señal en tiempo discreto x( n) , de longitud finita L , la Transformada de Fourier en Tiempo Discreto viene dada por: L −1

X (ω ) = ∑ x(n)e − jωn ,

0 ≤ ω ≤ 2π

n =0

Si se muestrea

X (ω ) en frecuencias equiespaciadas ω k =

las muestras resultan:

2πk , k = 0,1,", N − 1 , y donde N ≥ L , N

− j 2 kπn ⎛ 2πk ⎞ L −1 N X (k ) = X ⎜ ⎟ = ∑ x ( n )e N ⎝ ⎠ n =0

o equivalentemente N −1

X ( k ) = ∑ x ( n)e

− j 2 kπn

N

,

k = 0,1,", N − 1

n =0

ya que x(n) = 0 para L ≤ n ≤ N-1. La última ecuación constituye la definición de la denominada Transformada Discreta de Fourier con N puntos ( N -point DFT) [2]. Existen diversas formas eficientes de implementar computacionalmente la DFT. Una de las más difundidas es la denominada Transformada Rápida de Fourier (FFT: Fast Fourier Transform), que permite computar la DFT con un tiempo de ejecución reducido. En Matlab, la FFT está implementada en la función fft. Esta función emplea un algoritmo radix-2 si la longitud de la secuencia es una potencia entera de 2 (ver [2], Sección 9.3.3 en Cap 9, para mayores detalles), y un algoritmo más lento si no lo es. Dada una señal almacenada en el vector x, el siguiente comando >> y = fft(x);

TeSyS - TP No. 3

6

computará la DFT con N=length(x) puntos de la secuencia x. El número de puntos de la DFT puede especificarse independientemente de la longitud de la secuencia mediante un argumento de entrada adicional de la función. Por ejemplo, los comandos para computar la FFT con N=100 puntos serían: >> N = 100; >> y = fft(x,N); Si N es menor que la longitud de x, entonces el vector x es completado con ceros hasta la longitud N. En caso de ser mayor, la secuencia es truncada a una longitud igual a N. Nota: La función fft no genera el vector de frecuencias, que debe ser generado independientemente. La función devuelve los valores de la DFT que corresponderían a frecuencias entre 0 y 2π.

El retorno al dominio temporal puede realizarse mediante el cómputo de la IDFT, utilizando para ello la función Matlab ifft, de la siguiente manera: >> x = real(ifft(y)); donde el vector y corresponde a la DFT de la señal x, y x es la señal temporal original. Cabe aclarar que se toma la parte real de la señal obtenida mediante la función ifft, dado que la señal resultante puede presentar valores complejos con parte imaginaria muy pequeña debido a errores númericos. También puede especificarse el número N de puntos considerados en el cálculo de la IDFT, el cual debe coincidir con el valor de N usado previamente en la función fft. Si se omite el valor N, la función ifft utilizará por defecto N=length(y). Nota: La función ifft de Matlab asume que el espectro se encuentra definido entre 0 y 2π. Los siguientes ejemplos ilustran el uso de la funciones fft e ifft. Ejemplo 1: Supongamos que se desea calcular la DFT de la señal en tiempo discreto x(n) que se obtiene al muestrear la señal:

x(t ) = 4 + 3sin(2π 5t)

con una frecuencia Fs = 100 Hz. La siguiente sucesión de comandos Matlab permite generar la señal x(n) y calcular y graficar la amplitud y fase de la DFT de x(n) . t=0:0.01:1; x=4+3*sin(5*pi*2*t); y=fft(x,10000); magnitud=abs(y); fase=unwrap(angle(y)); frec=[0:2*pi/10000:2*pi]; frec=frec(:,1:10000); subplot(211) plot(frec,magnitud);xlabel('frecuencia w [rad/seg]'); ylabel('Magnitud'); title('Espectro de x(n)');axis([0 2*pi 0 200]);grid; subplot(212) plot(frec,fase*180/pi);xlabel('frecuencia w [rad/seg]'); ylabel('Fase [grados]'); axis([0 2*pi -4e4 2e4]);grid; La salida gráfica correspondiente está representada en la Figura 3.1.1.

TeSyS - TP No. 3

7

Espectro de x(n) 200

Magnitud

150

100

50

0

0

1

2

3

4

5

6

4

5

6

frecuencia w [rad/seg] 4

2

x 10

Fase [grados]

1 0 -1 -2 -3 -4

0

1

2

3

frecuencia w [rad/seg]

Figura 3.1.1: Salida gráfica correspondiente al Ejemplo 1. Ejemplo 2: La siguiente secuencia de comandos permite obtener nuevamente la señal temporal a partir del cálculo de la IDFT, visualizando en una misma gráfica la señal original y el resultado obtenido a partir de la IDFT: x_new=real(ifft(y,10000)); figure, plot(t,x,'r-o',t,x_new(1:length(x)),'b-x'); title('Señal original vs. resultado de la IDFT'); xlabel('tiempo [segundos]'); ylabel('Amplitud'); legend('Señal original','Resultado de la IDFT'); grid; La salida gráfica correspondiente se visualiza en la Figura 3.1.2.

Figura 3.1.2: Salida gráfica correspondiente al Ejemplo 2. TeSyS - TP No. 3

8

Problema 4: El archivo de datos tp3_1.mat contiene un vector x que consiste de las muestras de una señal senoidal inmersa en ruido, muestreada con una frecuencia de Fs = 1 KHz. Se asume que el ruido es blanco Gaussiano. a.

Use la función fft de Matlab para computar y graficar la DTFT de la señal x, usando una ventana rectangular de longitud apropiada. Grafique la amplitud de la DTFT en por unidad, y la fase en grados, con el eje de abscisas calibrado en frecuencias continuas.

b. Basándose en las gráficas obtenidas determine la amplitud y frecuencia de la señal senoidal. Justifique teóricamente los cálculos realizados. c.

Repita los apartados a. y b., pero usando ahora una ventana de Hann. Compare los resultados con los obtenidos en a. y b.

d. Repita los apartados a. y b. pero usando una longitud de ventana igual a la mitad de la usada en a. y b. Compare los resultados obtenidos. ‰

3.2. Problemas a Presentar en el Informe Problema 5: Procesamiento de señales de electrocardiograma (ECG) El objetivo de este problema consiste en demostrar algunas de las tareas de pre-procesamiento requeridas por señales electrocardiográficas provenientes de la actividad eléctrica del corazón, antes de que la señal pueda ser usada para diagnóstico automático. Se pretende además determinar el ritmo cardíaco del paciente. Introducción El ciclo cardíaco comprende una contracción secuencial de las aurículas y los ventrículos del corazón. La actividad eléctrica combinada de las diferentes células miocárdicas produce corrientes eléctricas que se difunden a través de los fluidos corporales. Estas corrientes son tan altas que producen diferencias de potencial detectables entre distintas partes de la superficie del cuerpo. La señal registrada como la diferencia entre dos potenciales sobre la superficie corporal se llama “ECG lead” (señal de ECG). La señal usada en este problema fue registrada por medio de cinco electrodos colocados sobre los brazos y piernas de un voluntario, los cuales miden las diferencias de potencial entre los brazos izquierdo y derecho, y entre el pie izquierdo y el brazo izquierdo. El quinto electrodo va a tierra y se colocó sobre la pierna derecha del paciente.

Figura. 3.2.1. Definición de ondas, segmentos e intervalos en una señal normal de ECG. El patrón regular de picos producido por el corazón se repite para cada latido (ver Figura 3.2.1). Este patrón de picos potenciales se llama electrocardiograma o ECG. El pequeño pico inicial, P, marca la contracción de las aurículas. Los picos más grandes que siguen a P, el complejo QRS, son una superposición de la relajación de las aurículas y la contracción ventricular. Los potenciales eléctricos que surgen de la relajación auricular son, sin embargo, mucho menores que los potenciales originados por la contracción ventricular. Por esta razón, se suele considerar que el complejo QRS marca la contracción ventricular. A continuación del complejo QRS sigue la onda T asociada con la relajación ventricular. TeSyS - TP No. 3

9

El pulso Normalmente el corazón de un adulto se despolariza 60 a 90 veces por minuto. Una tasa de despolarización inferior se llama bradicardia, mientras que una superior se denomina taquicardia. El pulso normal de un recién nacido es mucho más alto que el de un adulto. El pulso se puede calcular como la cantidad de intervalos R-R contenidos en un minuto. Este número se denota BPM (Beats Per Minute – Latidos por minuto). Eliminación de señales provenientes del movimiento muscular La señal ecg.mat está muestreada a una frecuencia de 900 Hz y sus valores se encuentran expresados en volts. El factor de amplificación en las señales del electrodo es de 200. a.

Grafique la señal medida con las unidades apropiadas en los ejes.

El paciente por lo general se mueve un poco, lo cual agrega una señal de baja frecuencia proveniente de los músculos. Esta señal se superpone con la señal de ECG y dificulta el estudio de sus detalles. b. En el dominio de la transformada de Fourier, elimine la señal de movimiento muscular anulando las muestras correspondientes a las frecuencias por debajo de 0,8 Hz (no deberá anular la frecuencia F = 0 Hz correspondiente a la componente de continua). Calcule la transformada inversa de Fourier para volver al dominio temporal. Observación: Tenga en cuenta que la función ifft de Matlab asume que el espectro se encuentra definido entre 0 y 2π. c. Grafique en una misma figura (utilizando el comando subplot) la señal original y la nueva señal resultante del filtrado ecg1. Explique sus diferencias. Eliminación de interferencia de línea Frecuentemente el ECG se encuentra contaminado con ruido de línea, que se superpone con la señal de ECG. d. Anule la/s muestra/s correspondiente/s a la frecuencia de 60 Hz. e. Grafique en una misma figura la señal ecg1 y la señal ecg2 resultante del nuevo filtrado. Explique sus diferencias. Determinación del ritmo cardíaco f. Con la señal ya filtrada calcule el ritmo cardíaco del paciente y determine si sufre de taquicardia, bradicardia o se encuentra en situación normal. Problema 6: a.

Para cada una de las siguientes ventanas, i. Ventana Rectangular (boxcar) ii. Ventana Hann (hanning) iii. Ventana Hamming (hamming) Utilice la función Matlab implementada en el Problema 1 para computar los espectros de amplitud. En una misma gráfica muestre superpuestos los espectros para tres longitudes L distintas, L = 20,31,81. Use una escala lineal para el espectro de amplitud.

b. Para cada una de las ventanas en el punto a, mida la amplitud y el ancho del lóbulo principal de la DTFT. Analice los puntos de cruces por cero y la amplitud de los lóbulos laterales. En el caso de la ventana rectangular, verifique la ubicación de los puntos de cruces por cero a partir del cálculo analítico de la DTFT de un pulso rectangular. NOTA: Recuerde que puede usar el comando ginput para extraer valores de una gráfica en Matlab. c.

Muestre superpuestos, en una misma gráfica, los espectros de amplitud de las ventanas rectangular, Hann y Hamming para una longitud L=45.

d. La existencia de los lóbulos laterales en el espectro de la ventana, provoca que la potencia de la señal truncada usando la ventana se distribuya en todo el rango de frecuencias. Compare la altura de los lóbulos laterales para las tres ventanas consideradas en el apartado c.

TeSyS - TP No. 3

10

e.

Considere la siguiente señal de longitud infinita

⎛π ⎞ x(n) = cos⎜ n ⎟ ⎝2 ⎠ Utilice una ventana rectangular de longitud L=45 para computar una aproximación del espectro de la señal y graficarlo. Muestre, en la misma gráfica, el espectro teórico de la señal. Compare ambos espectros. f.

El ancho del lóbulo principal del espectro de la ventana es un indicador de la resolución en frecuencia que se obtiene cuando se usa la ventana para analizar el espectro de una señal. En otras palabras, el ancho del lóbulo principal está relacionado con la mínima separación en frecuencia entre dos componentes armónicos de la señal, tal que el espectro presente picos bien diferenciados correspondientes a cada armónico. Esto puede verse fácilmente si se usa la ventana para computar el espectro de una señal compuesta por dos armónicos de frecuencias próximas. Sea la señal

x(n) = cos(0.5π .n ) + cos(0.51π .n )

Use la ventana rectangular del apartado e. para computar y graficar una aproximación del espectro de la señal x(n). Superponga en la gráfica el espectro teórico de la señal. Determine la relación entre la resolución en frecuencia y la longitud L de la ventana. A partir de este resultado, calcule la longitud L mínima de modo que se distingan las componentes armónicas de la señal x(n). g.

Analice qué ocurriría en caso de utilizar las ventanas de Hann y de Hamming para truncar la señal del apartado f. Calcule, para estas dos ventanas, la longitud L mínima necesaria para distinguir ambas componentes de la señal.

Problema 7: Análisis frecuencial utilizando ventanas El uso de la Transformada de Fourier es de fundamental importancia para el análisis frecuencial de señales. Sin embargo, si se calcula la transformada de Fourier de una porción de larga duración de la señal, se pierde la información temporal de su contenido frecuencial, es decir, se conoce el contenido armónico de la señal pero no se puede especificar en qué momento está presente cada armónico. En las aplicaciones en las cuales se requiere conocer la información frecuencial en función del tiempo de la señal, se suele realizar un análisis de Fourier de tiempo corto. Para ello se divide la señal en segmentos de corta duración denominados frames. Un frame es una porción de L muestras sucesivas de una señal de longitud M, que se obtiene a partir de aplicar una ventana a la señal y desplazarla temporalmente. En la Figura 3.2.2, se divide una señal de longitud M, en N frames de longitud L. De esta manera, calculando la Transformada de Fourier a cada frame se puede conocer la evolución del contenido armónico de la señal en el tiempo. Esta técnica se denomina Transformada de Fourier en Tiempo Corto (Short-Time Fourier Transform).

Señal

1

2

k

(k-1)*L + 1

...

N-1

N

k*L

Figura. 3.2.2. División en frames de una señal. TeSyS - TP No. 3

11

Sistema de Marcación por Tonos En el sistema de marcación por tonos utilizado en telefonía, también llamado sistema multifrecuencial o DTMF (Dual-Tone Multi-Frequency), cuando un usuario pulsa en el teclado de su teléfono la tecla correspondiente a un dígito que quiere marcar, se envían dos tonos de distinta frecuencia (uno por columna y otro por fila de acuerdo a la Tabla 3.2.1). Tabla 3.2.1: Frecuencias asociadas a cada dígito del sistema DTMF. FH [Hz]

FL [Hz]

697 770 852 941

1209

1336

1477

1633

1 4 7 *

2 5 8 0

3 6 9 #

A B C D

Así, la señal generada al presionar un dígito, tendrá dos componentes y una forma genérica:

x ( t ) = A sen ( 2π FL t ) + A sen ( 2π FH t )

FS ⎯⎯ →

x ( n ) = A sen ( 2π f L n ) + A sen ( 2π f H n )

Luego, la central telefónica detecta las frecuencias contenidas en la señal y determina el dígito que se marcó. En los casos en que las frecuencias FL y FH difieren de sus valores nominales, indicados en la Tabla 1, en ±1.8% , la central telefónica descarta el dígito enviado. El objetivo de este problema es analizar el contenido frecuencial de una señal asociada al marcado de un número telefónico, por medio de la Transformada de Fourier en Tiempo Corto, para identificar el número telefónico marcado, determinando cada uno de los dígitos que lo componen. a.

Use la función Matlab wavread para cargar el archivo “tonos.wav” que contiene la señal de audio a analizar. Sintaxis: [Y, Fs, Nbits] = wavread('tonos.wav'); Donde Y son las muestras de la señal, Fs la frecuencia muestreo y Nbits la cantidad de bits por muestra utilizados. b. Con la ayuda de las funciones fft y fftshift de Matlab compute la DTFT de la señal Y completa y grafique el espectro de amplitud en función de la frecuencia en tiempo continuo asociada, en el rango entre [-Fs/2 , Fs/2]. Indique si es posible o no determinar el número marcado a partir del espectro de amplitud de la señal completa. c. Realice una función que tome como argumentos de entrada dos frecuencias, F1 y F2 (F1 < F2), y como argumento de salida la variable digito, correspondiente al dígito identificado. Para ello compare F1 con cada una de las frecuencias FL y F2 con cada una de las frecuencias FH. En caso de que alguna de las frecuencias estuviera fuera del rango (FL ± 1.8% y FH ± 1.8%), la variable de salida debe ser digito=[]. d. Realice un script que: i. Divida la señal en frames con una duración de 450mseg, tal que la longitud del frame resulta L = (450e-3)*Fs muestras. ii. Para cada frame, calcule la DTFT, determine automáticamente las dos frecuencias de mayor amplitud y calcule las frecuencias en tiempo continuo asociadas (en Hz) (para ello puede utilizar la función max de Matlab). Utilice la función del ítem c. para determinar el dígito correspondiente a este frame. iii. Una vez procesados todos los frames, muestre en pantalla el número identificado. e. Suponga ahora que la señal se envía a través de un canal de transmisión con una respuesta al impulso h(n), como se indica en la Figura 3.2.3.

TeSyS - TP No. 3

12

y(n)

h(n)

yfiltrada(n)

Figura 3.2.3. Canal de transmisión. El canal de transmisión es modelado con una estructura FIR (Finite Impulse Response) de fase lineal y de longitud 80. Obtenga la respuesta al impulso del canal a partir del comando h=fir1(80,0.325) y filtre la señal contenida en el archivo “tonos.wav”. Luego, determine nuevamente el número telefónico marcado a partir de la señal filtrada. Compare los resultados con los obtenidos en el ítem d. y extraiga conclusiones. Ayuda: grafique la respuesta en frecuencia del canal.

4. Referencias [1] Burrus, C. & McClellan, J. & Oppenheim, A. & Parks, T. & Schafer, R. & Schuessler, H. (1994). Computer-based Exercises for Signal Processing using Matlab, Prentice Hall, Englewood Cliffs, New Jersey. [2] Proakis, J. & Manolakis, D. (1992). Digital Signal Processing: Principles, Algorithms and Applications. 2nd Edition, Macmillan Publishing Company, New York. [3] The Math Works, Inc. (1996). Signal Processing Toolbox - User's Guide, V.4. The Math Works, Inc., 24 Prime Park Way, Natick, MA 01760-1500. [4] The MathWorks, Inc. (1997). Matlab, The Language of Technical Computing – Getting Started with Matlab, Version 5., 24 Prime Park Way, Natick, MA 01760-1500.

TeSyS - TP No. 3

13