1
Tarea #1 ELO-313 Procesamiento digital de señales con aplicaciones Felipe Ávila Cárcamo 201121002-6 I. Evaluación de propiedades de sistemas Para el presente problema, se estudian 3 archivos de MATLAB precompilados bbox1.p, bbox2.p y bbox3.p, los cuales son 3 sistemas desconocidos o cajas negras, y de las cuales se quieren hacer pruebas para determinar propiedades características, como linealidad, invarianza en el tiempo, causalidad y estabilidad. Linealidad: Para probar la linealidad de los sistemas estudiados, se sigue la sugerencia dada de excitar los sistemas con ruido blanco, para así asegurar excitaciones ricas en frecuencia. Para esto se toman dos entradas, llamadas x1 (t) y x2 (t), generadas con el comando de MATLAB wgn(m, n, p) siendo sus argumentos m, n la dimensión del vector (o matriz), y p la potencia en dB con respecto a 1[W ]. A cada entrada se le asocia una salida de la forma yn = bboxN (xn ), con n = 1, 2 y N = 1, 2, 3. A su vez se emplea una tercera entrada denominada y3 = bboxN (x3 ), con x3 (t) = a · x1 (t) + b · x2 (t), entonces si se tiene que y3 6= a · y1 + b · y2 para algún a, b, entonces el sistema se considera no lineal, o sea: bboxN (a · x1 (t) + b · x2 (t)) 6= a · bboxN (x1 (t)) + b · bboxN (x2 (t)) ⇒ bboxN (·) Es no lineal
(1)
En las siguientes imágenes, se observarán en azul el resultado de a · y1 + b · y2 y en rojo el resultado de y3 . Respuestas de bbox1() a ruido blanco 100 ay1 + by2 y3 50
Amplitud
0
−50
−100
−150
0
0.2
0.4
0.6
0.8
1
Tiempo [s]
(a) Respuestas obtenidas
(b) Respuestas con zoom
Figura 1: Prueba de linealidad sobre bbox1() con 2 entradas de ruido blanco. Se observa que no hay diferencia alguna (no se aprecia la gráfica en azul dado que coincide con la respuesta en rojo), aunque esta prueba no es suficiente para afirmar linealidad de bbox1().
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
2
4
4
Respuestas de bbox2() a ruido blanco
x 10
ay1 + by2 y3 3 2
Amplitud
1 0 −1 −2 −3 −4
0
0.2
0.4
0.6
0.8
1
Tiempo [s]
(b) Respuestas con zoom
(a) Respuestas obtenidas
Figura 2: Prueba de linealidad sobre bbox2() con 2 entradas de ruido blanco. Tal como en el caso de bbox1(), ambas respuestas de bbox2() son similares, pero no basta para afirmar linealidad del sistema. Respuestas de bbox3() a ruido blanco
Respuestas con zoom de bbox3() a ruido blanco
40
15 ay1 + by2 y3
ay1 + by2 y3
30 10 20 5 Amplitud
Amplitud
10 0
0
−10 −5 −20 −10 −30 −40
0
0.2
0.4
0.6
0.8
1
−15
4
4.2
4.4
4.6
Tiempo [s]
(a) Respuestas obtenidas
4.8 5 Tiempo [s]
5.2
5.4
5.6
5.8 −5
x 10
(b) Respuestas con zoom
Figura 3: Prueba de linealidad sobre bbox3() con 2 entradas de ruido blanco En este caso, claramente se aprecia en las respuestas con el zoom hecho que existen discrepancias, por lo tanto el sistema bbox3() es un sistema no lineal, pues basta un caso en el cual exista discrepancia entre las respuestas analizadas para hacer la afirmación. Para los sistemas bbox1() y bbox2() no bastaba con esta prueba para afirmar linealidad, pues las respuestas estudiadas deben ser iguales para cualquier tipo de señal en el tiempo. Invarianza en el tiempo: Para hacer la prueba de invarianza, se debe demostrar el supuesto que para cualquier parámetro τ ∈ R que: bboxN (x(t)) = y(t) ⇒ bboxN (x(t − τ )) = y(t − τ )
(2)
O sea, que un corrimiento (hacia la izquierda o derecha) en la entrada del sistema, genera la misma salida corrida la misma cantidad de tiempo que la original. Para esto se usan entradas sinusoidales en los sistemas, y se obtienen las siguientes figuras, siendo en azul la respuesta y(t − 0.2), siendo y(t) = bboxN (x(t)) la salida antre la entrada sinusoidal sin corrimiento, y en rojo la respuesta de la sinusoide con corrimiento hacia la derecha de 0.2[s] bboxN (x(t − 0.2)):
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
3
Prueba de invarianza para bbox1 6 y(t − 0.2) bbox1(x(t − 0.2)) 4
Amplitud
2
0
−2
−4
−6
0
0.2
0.4
0.6
0.8
1
Tiempo [s]
Figura 4: Respuestas obtenidas para bbox1().
Prueba de invarianza para bbox2 50 y(t − 0.2) bbox2(x(t − 0.2)) 40
Amplitud
30
20
10
0
−10 0.1
0.2
0.3
0.4
0.5 Tiempo [s]
0.6
0.7
0.8
0.9
Figura 5: Respuestas obtenidas para bbox2().
Prueba de invarianza para bbox3 1 y(t − 0.2) bbox3(x(t − 0.2))
0.8 0.6 0.4
Amplitud
0.2 0 −0.2 −0.4 −0.6 −0.8 −1
0
0.2
0.4
0.6
0.8
1
Tiempo [s]
Figura 6: Respuestas obtenidas para bbox3(). Claramente, para las cajas negras bbox1 y bbox3 no se aprecian los dos gráficos pues al estar superpuestos y ser exactamente iguales, sólo se observa un color, por lo que esta prueba no es suficiente para decir que son variantes en el tiempo (y tampoco
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
4
invariantes, pues la demostración debe ser para cualquier corrimiento, y para cualquier tipo de entrada). Por otra parte, las respuestas obtenidas para bbox2 presentan diferencias, las cuales basta para afirmar que el sistema es variante en el tiempo. Causalidad y estabilidad. En esta sección, se evaluará la causalidad y estabilidad de los sistemas. Para esto, se excitarán los sistemas con respuestas a impulsos, siendo el impulso generado como [zeros(1, 29), 1, zeros(31, 60)]. Los resultados obtenidos fueron los siguientes: Respuestas de bbox1() a impulso
Respuestas de bbox2() a impulso
1.5
1 0.8 0.6 0.4
1 Amplitud
Amplitud
0.2 0 −0.2 0.5
−0.4 −0.6 −0.8
0
0
10
20
30 Muestra
40
50
−1
60
0
(a) Respuesta a impulso bbox1
10
20
30 Muestra
40
50
60
(b) Respuesta a impulso bbox2 Respuestas de bbox3() a impulso
1 0.8 0.6 0.4
Amplitud
0.2 0 −0.2 −0.4 −0.6 −0.8 −1
0
10
20
30 Muestra
40
50
60
(c) Respuesta a impulso bbox3
Figura 7: Respuestas a impulso de los sistemas dados. Comenzando el análisis de los gráficos mostrados con el sistema bbox1 dado que es el único sistema invariante y lineal, se observa que la respuesta tiene forma de una progresión geométrica (analizando los datos) desde el instante en que se aplica √ el impulso (muestra [k] = 30). Siendo en un principio igual a 2 y disminuyendo con razón √ de 0.74 por cada muestra. Por lo que se puede decir que su respuesta a impulso en este caso queda descrita por h1 = 2(0.74)n−30 µ(30 − n), con µ() la función escalón unitario. Para el caso la condicón de estabilidad BIBO puede ser probada mediante: ∞ X
|hN [n]| < ∞
(3)
n=−∞
Y dado que la progresión geométrica en este caso es de razón estrictamente menor a 1, entonces la sumatoria presentada es acotada, por lo que el sistema es estable. Cabe recalcar como el sistema eso sí no es causal, dado que la respuesta entregada aparece en instantes anteriores al cual el impulso es recién aplicado. Basta este contraejemplo para asegurar que bbox1 entonces es no causal. En resumen entonces bbox1 hasta ahora es lineal, invariante, estable en sentido BIBO pero no es causal, siendo un ejemplo claro de la independencia de estas propiedades en sistemas. Por otra parte para bbox2(), tal como en el ejemplo anterior se aprecia que existe respuesta en instantes anteriores al cual se aplica el impulso, por lo que basta como contra ejemplo para asegurar que el sistema no es causal. A su vez, dado que es un sistema variante en el tiempo, no basta la respuesta a impulso para asegurar su estabilidad. Por último, para bbox3() presenta una respuesta a impulso igual a cero en todo instante para las muestras adquiridas. Dado que el sistema es no lineal, no se puede concluir ni estabilidad o causalidad con el gráfico obtenido.
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
5
II. Inversión de sistemas Sea un sistema y = SA (x) dado por y[n] =
1 y[n 2
− 1] − x[n].
Parte A: Se pide obtener una expresión para el sistema inverso y = SB (x) de modo tal que δ[n] = SB (SA (δ[n])), donde δ[n] representa el impulso en tiempo discreto. Siendo un sistema lineal, se quiere entonces encontrar un SB (x) tal que anule el efecto de SA (x) para cualquier x[n]. Para obtener el sistema se usa la herramienta de la transformada Z, que en el caso es útil, dado que como el sistema SA (x) es lineal, basta encontrar su respuesta a impulso con la transformada Z despejando el cociente Y (z)/X(z), para luego así invertir la expresión y volver a componer la ecuación de diferencias, siendo la que represente entonces al sistema SB (x). Entonces, aplicando la transformada Z se tiene: 1 y[n − 1] − x[n]/Z 2 1 −1 Y (z) = z Y (z) − X(z) 2 Y (z) 1 = HA (z) = 1 −1 X(z) z −1 2 y[n] =
(4) (5) (6)
Por lo tanto y como se dijo anteriormente, se desea entonces un HB (z) tal que es (HA (z))−1 , entonces: HB (z) =
Y (z) 1 1 −1 = z −1 − 1 ⇒ Y (z) = X(z) z −1 X(z) 2 2
1 x[n − 1] − x[n] 2 Siendo entonces la última ecuación presentada el sistema inverso SB (x) buscado. ∴ y[n] =
(7) (8)
Parte B: Se pide crear funciones en MATLAB que implementen los sistemas SA (x) y SB (x) para cualquier entrada x[n], y posterior obtener las respuestas de impulso de SA , SB y SB (SA ). Para el desarrollo de las funciones, se considera únicamente entradas causales, y que para ambos sistemas SA y SB se inicializa la salida en un principio, siendo para las muestras n = 0 y[0] = x[0] = 0, y como x es causal, para las muestras n < 0 se tiene salidas nulas. Las funciones se pueden observar en el anexo, sección II. Con estas funciones, graficando las respuestas a impulso se obtuvieron los siguientes resultados:
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
6
Respuesta a impulso del sistema SB (x) 0.5
Amplitud
0
−0.5
−1
0
5
(a) Sistema SA (x)
10
15 Muestra
20
25
30
(b) Sistema SB (x) Respuesta a impulso del sistema SB (SA (x))
1 0.9 0.8 0.7
Amplitud
0.6 0.5 0.4 0.3 0.2 0.1 0
0
5
10
15 Muestra
20
25
30
(c) Sistema SB (SA (x))
Figura 8: Respuestas a impulso de los sistemas pedidos. En los casos, el impulso se consideró sin ceros a la izquierda del instante en que se dispara, dado que los sistemas son causales. Observando el tercer gráfico de la figura (8c) se aprecia que se obtiene la entrada original, demostrando así que el sistema SB es el inverso de SA , logrando recuperar el impulso de tiempo discreto. III. Integrales y derivadas en tiempo discreto Se planea estudiar los siguientes sistemas en tiempo discreto: d • La derivada dado por el sistema y = S1 (x) donde y(t) = ´ t dt x(t). • La integral dada por el sistem y = S2 (x) donde y(t) = x(τ )dτ . −∞ Parte A: En esta sección se busca formular sistemas en tiempo discreto que aproximen los formulados en tiempo continuo presentados en el enunciado S1 y S2 . Se comenzará con el análisis para la derivada, la cual se puede escribir aproximadamente en tiempo discreto haciendo uso del tiempo de muestreo T como: y(kT ) =
dx(t) x(kT ) − x((k − 1)T ) ≈ dt t=kT T
(9)
Haciendo uso de la notación discreta en que n = kT (dado que T se supone fijo, y k es simplemente el número de la muestra), se tiene entonces el sistema como: x[n] − x[n − 1] y[n] = (10) T Ahora bien, para la integral se usa la aproximación mediante la serie de Riemann, quedando expresado como: y(kT ) =
k X
x(uT )T
(11)
x[u]T
(12)
u=−∞
Usando nuevamente la notación discreta, se tiene entonces que: y[n] =
n X u=−∞
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
7
Ahora bien, trabajando la expresión: y[n] = T
n−1 X
! x[u] + x[n]
= y[n − 1] + x[n]T
(13)
u=−∞
Siendo la ecuación (13) la resultante de la aproximación de la integral en sistema de tiempo discreto, teniendo así las ecuaciones de diferencias para los sistemas S1 y S2 . Cabe destacar en primera instancias que las ecuaciones obtenidas corresponden mutuamente a sus inversas, lo cual es válido ya que la derivada y la integral son operaciones inversas de por sí. Otro comentario importante de los resultados obtenidos es que no son únicos, dado que por ejemplo en el caso se empleo para la integral la aproximación por series de Riemann o sumas de áreas rectangulares, cuando también existe por ejemplo mediante regla trapezoidal, y a su vez la derivada en vez de haberla trabajado con datos adyacentes, haberla realizado como la diferencia de sí misma con una muestra anterior, pero no lo antecesora, si no dos o tres muestras atras, y dividida la diferencia por el factor correspondiente del periodo de muestreo entre ambas muestras. Ahora bien para analizar la estabilidad, basta analizar las transformadas Z de las ecuaciones (10) y (13). Procediendo para la derivada se obtiene: x[n] − x[n − 1] Y (z) 1 − z −1 y[n] = /Z ⇒ = (14) T X(z) T Se observa que el polo se encuentra en z = 0, siendo dentro del círculo unitario en el plano complejo, por lo que es estable el sistema. A su vez, se puede considerar la estabilidad vista en (3) donde se tiene que ante un impulso, la respuesta de la salida sería simplemente { T1 , T1 }, por lo que la suma de esto es acotado, dado que el periodo de muestreo T 6= 0. Ahora bien, para la integral se tiene: Y (z) T y[n] = y[n − 1] + x[n]T /Z ⇒ = (15) X(z) 1 − z −1 Claramente en este caso, se posee un polo en z = 1, esto es sobre el circulo unitario en el plano complejo, por lo que se considera el sistema como inestable. Esto se pudo haber deducido de igual considerando que ante una entrada acotada como un escalón, la integral nos arroja una rampa que claramente no es acotada y diverge. Parte B: Para esta sección se pide generar una señal x(t) = sin(2 · π · 100 · t), duración de un segundo y aplicarle los sistemas de derivada discreta e integral discreta. Para escoger el periodo de muestreo de la señal, se decide que si se quiere tener 100 puntos por periodo, y el periodo dura 10 [ms], entonces con un periodo de 100 [µs] se logra lo buscado. Se diseñan funciones en MATLAB que implementan los sistemas discretos de S1 y S2 , y que tienen de argumento la señal a derivar (o integrar) y el periodo de muestreo. A su vez se considera que y[0] = x[0] pues se toma en cuenta que la sinusoide de entrada a los sistemas es causal (igual a cero para todo tiempo menor a cero). En el anexo del presente documento aparece en detalle el código usado. Los gráficos obtenidos son los siguientes:
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
8
Sinusoide de 100[Hz] x(t)
−3
Derivada discretada de x(t)
1
800
3.5
0.8 600
X: 0.005 Y: 0.003182
3
X: 0.01 Y: 627.9
Integración discreta de x(t)
x 10
0.6 400
2.5
200
2
0.4
0
Amplitud
Amplitud
Amplitud
0.2
0
1.5
−0.2 −200
1
−400
0.5
−600
0
−0.4
−0.6
−0.8
−1
0
0.01
0.02 0.03 Tiempo [s]
0.04
0.05
−800
0
0.01
0.02 0.03 Tiempo [s]
0.04
0.05
−0.5
0
0.01
0.02 0.03 Tiempo [s]
0.04
0.05
Figura 9: Gráfica de los primeros 50[ms] de la señal sinusoidal, su derivada discreta e integral discreta. Analizando los resultados, se puede afirmar que lo obtenido es muy cercano a lo esperado para los sistemas originales. Para la derivada el resultado teórico es: d {sin(2 · π · 100 · t)} = 200πcos(2 · π · 100 · t) (16) dt Siendo entonces la derivada un coseno de amplitud 200π ≈ 628.31, y dando la amplitud para la derivada discreta muy cercano a este valor (627.9), y cumpliendo a su vez con ser un coseno claramente. Por otra parte, la integral (dado que se considera causal la entrada, entonces como la integral es continua, empieza en cero) teóricamente es: ˆ t ˆ t −cos(200πτ ) τ =t sin(2 · π · 100 · τ )dτ = sin(2 · π · 100 · τ )dτ = = 0.00159(1 − cos(200πt)) (17) 200π τ =0 −∞ 0 Siendo entonces siempre el valor mayor o igual a cero, y a su vez teniendo amplitud igual a 0.00318 , resultado que es muy similar al obtenido y que se ve en la figura (9). El resultado de la integral de igual manera es intuitivo, dado que como se considera entrada causal, se parte integrando el lóbulo positivo de la sinusoide llegando a un máximo, y luego disminuyendo hasta llegar a cero, para volver a ser positivo, teniendo entonces la integral un valor medio mayor a cero. Parte C: Ahora, considerando una entrada como x[n] = δ[n] − δ[n − 5] entre los tiempos de muestreos −10 ≤ n ≤ 20 se desea graficar la señal x[n] junto a la respuesta a los sistemas discretos de derivador e integrador. Para esto se usa las mismas funciones en MATLAB definidas con anterioridad y que se encuentran en el anexo. Se obtiene el siguiente resultado:
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
9
Derivada discretada de x(t)
Integración discreta de x(t) 1
0.8
0.8
0.9
0.6
0.6
0.8
0.4
0.4
0.7
0.2
0.2
0.6
0
Amplitud
1
Amplitud
Amplitud
Entrada x[n] = δ[n] − δ[n − 5] 1
0
0.5
−0.2
−0.2
0.4
−0.4
−0.4
0.3
−0.6
−0.6
0.2
−0.8
−0.8
0.1
−1 −10
−5
0
5 Muestra
10
15
−1 −10
20
−5
0
5 Muestra
10
15
0 −10
20
−5
0
5 Muestra
10
15
20
Figura 10: delta Parte D: Tal como la sección anterior, se considera una entrada en este caso como x[n] = µ[n] − µ[n − 11] entre los tiempos de muestreos −10 ≤ n ≤ 20 se desea graficar la señal x[n] junto a la respuesta a los sistemas discretos de derivador e integrador. Se toma en cuenta entonces que la entrada a los sistemas es 1 para los 0 ≤ n ≤ 10, y para cualquier otra muestra la entrada es cero. Se obtiene el siguiente resultado:
Entrada x[n] = µ[n] − µ[n − 11]
Derivada discretada de x(t)
1
1
0.9
0.8
0.8
0.6
0.7
0.4
0.6
0.2
Integración discreta de x(t) 12
10
0.5
Amplitud
Amplitud
Amplitud
8
0
0.4
−0.2
0.3
−0.4
0.2
−0.6
0.1
−0.8
6
4
2
0 −10
−5
0
5 Muestra
10
15
20
−1 −10
−5
0
5 Muestra
10
15
20
0 −10
−5
0
5 Muestra
10
15
20
Figura 11: escalon La salida del sistema del derivador discreto es la que se espera correspondiendo a un impulso en el momento en que se aplica cada escalón y conservando el signo. Por otro parte la salida del sistema del integrador discreto corresponde claramente UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
10
a una rampa que comienza a crecer cuando comienza el efecto del escalón y termina de crecer cuando se aplica el otro escalón de polaridad inversa, manteniendose constante, dado que deja de integrar valor alguno. Parte E: Para esta sección, se carga en MATLAB mediante el comando audioread() los archivos de audio music.wav y speech.wav, y se emplean de entrada al derivador e integrador discreto con los cuales se ha estado trabajando, y luego se describe cómo afecta esto en la señal de audio. Lo que se obtuvo fue interesante, dado que en primer lugar se aprecia una distorsión en ambas salidas (para derivador e integrador), pero aún se puede reconocer los sonidos, siendo entonces que las salidas de los sistemas obtenidos no son irreconocibles. Por otra parte se notó que la intensidad cambió de manera considerable. Para el derivador, lo obtenido tenía mayor volumen, o sea era más intensa la salida y lo que se explica a que el derivador discreto divide la diferencia por el tiempo de muestreo, entonces dado que el tiempo de muestreo es pequeño, entonces existe una amplificación en las muestras claro. Para el integrador, se notó que la intensidad disminuyó de manera abrupta, a lo cual se le encuentra explicación de manera análoga a lo descrito con anterioridad. Aparte de lo descrito con anterioridad, también se recalca que lo obtenido de la derivación de las señales de audio fueron más ruidosas, con alta ponderación de los sonidos agudos y poco efecto de los tonos graves. Por otra parte lo obtenido con el integrador se notó que los graves estaban más presentes, mientras que los tonos agudos disminuyeron y el ruido no se sintió de manera intensa. La explicación tiene que ver simplemente a que la derivación es una operación que favorece los cambios abruptos, considerándose entonces en término de filtraje un filtro pasa alto en términos frecuenciales y siendo el motivo que se acentúen los tonos agudos y el ruido (dado que el ruido es una señal aleatoria con grandes cambios en el tiempo), mientras que el integrador el fenómeno es similar pero para bajas frecuencias, actuando como filtro pasa bajo. Parte F: Ahora se busca, en función de las aproximaciones realizadas con anterioridad, encontrar la ecuación de diferencias de la doble derivada y = S1 (S1 (x)), la doble integral y = S2 (S2 (x)), la derivada de la integral y = S1 (S2 (x)) y la integral de la derivada y = S2 (S1 (x)). Dado que tanto la derivada S1 e integral S2 son sistemas LTI, lo pedido se puede obtener con facilidad usando las respuesta a impulso, considerando que sistemas en cascada es lo mismo que tener las funciones de impulso respectiva multiplicándose. Para la doble derivada se tiene ya su respuesta a impulso en (14). Por lo tanto se tendrá que: Hdd (z) = Hd (z)Hd (z) =
1 − 2z −1 + z −2 T2
(18)
Pasando entonces a la ecuación por diferencias se tiene entonces: 1 (x[n] − 2x[n − 1] + x[n − 2]) T2 Ahora bien, la respuesta a impulso de la doble integral se obtiene como: ydd [n] =
T2 1 − 2z −1 z −2 Por lo tanto, el sistema descrito en ecuación de diferencias para la doble integral es: Hii (z) = Hi (z)Hi (z) =
yii [n] = T 2 x[n] + 2yii [n − 1] − yii [n − 2]
(19)
(20)
(21)
Finalmente, la derivada de la integral y la integral de la derivada, dado que tanto derivada con integral son operaciones inversas y que se puede apreciar con que la transformada Z de las ecuaciones de diferencias vistas en (14) y 15) son recíprocas, entonces poner los sistemas en cascada anularían mutuamente sus efectos (pues es la multiplicación de sus respuestas a impulsos) sin importar el orden (dado que la multiplicación es conmutativo), darán como salida la señal de entrada, es decir: yid = x[n], ydi = x[n]
(22)
Ahora bien, si los sistemas mostrados anteriormente fuesen filtros, para el sistema más simple visto en (22) son filtros pasa todo, dado que la salida es igual a la entrada en todo instante y sin atenuación en ninguna frecuencia, entonces por términos de estabilidad BIBO, para cualquier entrada acotada tendrá una salida acotada, el sistema es estable. A su vez, el sistema derivador doble dado que conserva polos y ceros del sistema derivador discreto simple pero ahora con multiplicidad dos, tiene el polo en z = 0 por lo que su región de convergencia contiene el círculo unitario y es estable, y ahora atenuará aún más los efectos de baja frecuencia y aumentará las frecuencias altas siendo un sistema entonces pasa alto. Para el sistema de integral doble, se tiene entonces claramente que su región de convergencia ROC no contiene el círculo unitario, pues conserva el mismo polo que el sistema de integral discreto, sólo que ahora con el polo doble en z = 1. Para analizar el filtro, dado que el sistema como se dijo, conserva los polos y ceros del sistema integrador discreto, pero ahora con multiplicidad UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
11
dos, se puede concluir que será del tipo pasa bajo, y lo cual se nota en el efecto del T 2 que entonces, mientras más pequeño el periodo de muestreo, más se acentua el efecto de la entrada, siendo entonces como un efecto promedio de los valores pasados.
IV. Estimación de tiempo de retardo en radares Sea un radar que transmite una señal xa (t) y recibe una señal ya (t) donde: ya (t) = axa (t − td ) + va (t)
(23)
Donde va (t) es ruido aleatorio. Las señales en tiempo discreto son: y[n] = ax[n − D] + v[n]
(24)
Parte A: Se desea saber como poder calcular el retardo D usando la correlación cruzada ryx [m]. Primero, la correlación cruzada entre dos señales se define como: ∞ ryx [m] =
X
y[n]x[n − m]
(25)
n=−∞
Entonces al reemplazar y[n] en la ecuación anterior se tiene: ∞ X
ryx [m] =
(ax[n − D] + v[n])x[n − m] =
n=−∞
∞ X
ax[n − D]x[n − m] +
n=−∞
∞ X
v[n]x[n − m]
(26)
n=−∞
De esta ecuación se puede apreciar que la correlación en el caso está dada por dos sumatorias, una correlacionando la entrada con la señal aleatoria y otra correlacionando la señal de entrada consigo mismo pero retardada. La información que se busca se encuentra en la sumatoria formada por la entrada multiplicada por ella misma adelantada m muestras, y se puede demostrar que la sumatoria claramente será máxima, cuando m = D, o sea cuando x[n] se multiplique por ella misma en cada muestra, sin desfase o retardos, por lo que en cercanías de m = D esta sumatoria será máxima. El problema que presenta lo expuesto anteriormente, es que se tiene datos determinísticos (la sumatoria de las entrada con retardo multiplicada por si mismo adelantada m muestras) y datos aleatorios (la sumatoria que incluye la entrada adelantada m muestras por el ruido), por lo que existiría un grado de inexactitud en la obtención del valor de D, fuertemente ligado a la energía que posea el ruido asociado al sistema, a su vez al factor de ponderación a usado. Entonces a mayor a, es mayor la precisión de la estimación del retardo D. Para analizar si es posible usar la correlación ryy [m] para encontrar el retardo, se procede como antes, trabajando la ecuación como sigue. ryy [m] =
∞ X
y[n]y[n − m]
(27)
(ax[n − D] + v[n])(ax[n − D − m] + v[n − m])
(28)
n=−∞ ∞
=
X
n=−∞ ∞
=
X
a2 x[n − D]x[n − D − m] +
n=−∞
∞ X
ax[n − D]v[n − m] +
n=−∞
∞ X
av[n]x[n − D − m] + rvv
(29)
n=−∞
En este caso, se tiene un valor determinístico que es la correlación entre la señal con desfase consigo mismo adelantada m muestras, y los demas valores aleatorios, pero si se desea maximizar el valor de la correlación ryy que es de la manera en que se procedió anteriormente, se requerirá que m = 0, por lo que el retardo D no sería posible de obtener directamente. Por último, cabe recalcar que para el cálculo de D mediante el uso de la correlación cruzada ryx en términos de que si es óptimo en el sentido de mínimos cuadrados, se considera primero que el error en mínimos cuadrados es descrito por: Els (f ) =
n X (yk − f (xk ))2 k=1
n
(30)
DOnde yk es la señal misma de salida y f (xk la predicción realizada, y dado que en este caso hay presencia de ruido en el sistema, lo óptimo es usar algún método que minimice este error entre la señal misma, y la predicción. En el sentido propuesto, el método consiste en encontrar el retardo D óptimo de la entrada x[n] para que así el error muestra a muestra de la salida propiamente tal y la predicción realizada sea mínima.
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
12
Parte B: Sea una secuencia x[n] conocida como secuencia de Barker, dada por: x[n] = +1; +1; +1; +1; +1; 1; 1; +1; +1; 1; +1; 1; +1
(31)
2
y un ruido v[n] Guassiano de media cero y varianza σ = 0.01. Se escribe un código en MATLAB que genera una secuencia y[n], para muestras entre 0 ≤ n ≤ 199, con a = 0.9 y D = 20. Para esto, se crea una función sisr ad() que aplica el sistema discreto del radar descrito anteriormente (códigos se encuentran en el anexo). Ahora bien la gráfica obtenida es la siguiente Entrada x[n] sec. Barker
Salida y[n] del sist.con entrada x[n]
1
1.5
0.8 0.6
1
0.4 0.2
0.5
0 −0.2
0
−0.4 −0.6
−0.5
−0.8 −1
0
20
40
60
80
100 120 Muestras
140
160
180
−1
200
0
20
40
60
80
100 120 Muestras
140
160
180
200
Figura 12: Gráfico de entrada (secuencia de Barker) a la izquierda y salida a la derecha. Parte C: Ahora se busca obtener la gráfica de las correlaciones ryy [m] y correlación cruzada ry x[m]. Se crea una función que nos calcula la correlación y que se encuentra en el anexo adjunto al presente informe. Las gráficas obtenidas son las siguientes: Correlación ryy[m]
Correlación cruzada rxy[m]
14
12
12
10
10
8
8 6 6 4 4 2
2
0
0 −2
0
20
40
60
80
100 Lag m
120
140
160
180
−2
200
0
20
40
60
80
100 Lag m
120
140
160
180
200
Figura 13: Gráfico de autocorrelación a la izquierda, y a la derecha la correlación cruzada Para determinar la SNR, se determina el valor con la siguiente fórmula: SN R(y) = 10log10 Siendo la potencia definida como: Py =
n X
Pseñal Pruido
y[k]2
(32)
(33)
k=1
Se crea una función en MATLAB llamada SN Rdb() que se encuentra en el anexo que permite obtener rápidamente lo buscado. En el casose obtiene que para el ruido de varianza σ 2 = 0.01 da: SN R(y) = 7.8211[dB]
(34)
Parte D: Usando los códigos mencionados anteriormente, se tiene entonces primero para un σ 2 = 0.1 una salida como la siguiente:
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
13
Entrada x[n] sec. Barker
Salida y[n] del sist.con entrada x[n]
1
2
0.8 1.5 0.6 1
0.4 0.2
0.5
0 0
−0.2 −0.4
−0.5
−0.6 −1 −0.8 −1
0
20
40
60
80
100 120 Muestras
140
160
180
200
−1.5
0
20
40
60
80
100 120 Muestras
140
160
180
200
Figura 14: Gráfico de entrada (secuencia de Barker) a la izquierda y salida a la derecha para σ 2 = 0.1.
Correlación ryy[m]
Correlación cruzada rxy[m]
30
12
25
10 8
20
6 15 4 10 2 5
0
0
−5
−2
0
20
40
60
80
100 Lag m
120
140
160
180
200
−4
0
20
40
60
80
100 Lag m
120
140
160
180
200
Figura 15: Gráfico de autocorrelación a la izquierda, y a la derecha la correlación cruzada para σ 2 = 0.1. En el caso, la SNR es: SN R(y) = 1.9863[dB] Para un σ 2 = 0.1
(35)
El valor claramente disminuye, pues al tener una mayor varianza, la señal de ruido tendrá más potencia pues su variación con respecto en la media es mayor. Ahora bien para un σ 2 = 1 se tiene: Entrada x[n] sec. Barker
Salida y[n] del sist.con entrada x[n]
1
3
0.8 2 0.6 1
0.4 0.2
0
0 −1
−0.2 −0.4
−2
−0.6 −3 −0.8 −1
0
20
40
60
80
100 120 Muestras
140
160
180
200
−4
0
20
40
60
80
100 120 Muestras
140
160
180
200
Figura 16: Gráfico de entrada (secuencia de Barker) a la izquierda y salida a la derecha para σ 2 = 1.
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
14
Correlación ryy[m]
Correlación cruzada rxy[m]
250
10
200
5
150 0 100 −5 50 −10
0
−50
0
20
40
60
80
100 Lag m
120
140
160
180
200
−15
0
20
40
60
80
100 Lag m
120
140
160
180
200
Figura 17: Gráfico de autocorrelación a la izquierda, y a la derecha la correlación cruzada para σ 2 = 1. En el caso, la SNR es: SN R(y) = 0.135[dB] Para un σ 2 = 1
(36)
En el caso, disminuyó notablemente la SNR, dando a entender que la potencia de la salida es comparable a la del ruido mismo, o sea el ruido afecta en gran medida a la salida, pues sus variaciones son mucho mayores que los casos anteriores. También se resalta, que en las correlaciones cruzadas vistas anteriores, en todas el máximo se encontró para m = 20, y a su vez para la autocorrelación el máximo estaba en m = 0, que coincide con el análisis hecho en la parte A de la presente sección. A su vez, se ve que a mayor varianza, el máximo disminuía en magnitud respecto a los otros valores, y lo que coincide con el hecho que al incluir más variabilidad en la señal de ruido, más impreciso se vuelve el modo mostrado para calcular el retardo. Parte E: Se pide hacer todo lo mencionado anteriormente, pero para una secuencia Barker repetida 10 veces. Entonces se tienen los siguientes resultados usando los códigos mostrados anteriormente (simplemente ahora hay que modificar la entrada x[n]). Para el σ 2 = 0.01 se tiene: Entrada x[n] sec. Barker
Salida y[n] del sist.con entrada x[n]
1
1.5
0.8 1
0.6 0.4
0.5
0.2 0
0
−0.2 −0.5
−0.4 −0.6
−1
−0.8 −1
0
20
40
60
80
100 120 Muestras
140
160
180
200
−1.5
0
20
40
60
80
100 120 Muestras
140
160
180
200
Figura 18: Gráfico de entrada (secuencia de Barker repetida 10 veces) a la izquierda y salida a la derecha para σ 2 = 0.01.
Correlación ryy[m]
Correlación cruzada rxy[m]
120
120
100
100
80
80
60
60
40
40
20
20
0
0
−20
0
20
40
60
80
100 Lag m
120
140
160
180
200
−20
0
20
40
60
80
100 Lag m
120
140
160
180
200
Figura 19: Gráfico de autocorrelación a la izquierda, y a la derecha la correlación cruzada para σ 2 = 0.01.
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
15
La SNR en el caso es de: SN R(y) = 17.6277[dB] Para un σ 2 = 0.01
(37)
2
Por otro lado los gráficos obtenidos para σ = 0.1 son: Correlación ryy[m]
Correlación cruzada rxy[m]
120
120
100
100
80
80
60
60
40
40
20
20
0
0
−20
0
20
40
60
80
100 Lag m
120
140
160
180
−20
200
0
20
40
60
80
100 Lag m
120
140
160
180
200
Figura 20: Gráfico de entrada (secuencia de Barker repetida 10 veces) a la izquierda y salida a la derecha para σ 2 = 0.1.
Correlación ryy[m]
Correlación cruzada rxy[m]
140
120
120
100
100
80
80 60 60 40 40 20
20 0
0
−20
−20
0
20
40
60
80
100 Lag m
120
140
160
180
200
0
20
40
60
80
100 Lag m
120
140
160
180
200
Figura 21: Gráfico de autocorrelación a la izquierda, y a la derecha la correlación cruzada para σ 2 = 0.1. La SNR en el caso es de: SN R(y) = 7.985[dB] Para un σ 2 = 0.1
(38)
Y por último los gráficos obtenidos para σ 2 = 1 son los siguientes: Entrada x[n] sec. Barker
Salida y[n] del sist.con entrada x[n]
1
4
0.8 3 0.6 2
0.4 0.2
1
0 0
−0.2 −0.4
−1
−0.6 −2 −0.8 −1
0
20
40
60
80
100 120 Muestras
140
160
180
200
−3
0
20
40
60
80
100 120 Muestras
140
160
180
200
Figura 22: Gráfico de entrada (secuencia de Barker repetida 10 veces) a la izquierda y salida a la derecha para σ 2 = 1.
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
16
Correlación ryy[m]
Correlación cruzada rxy[m]
300
120
250
100
200
80
150
60
100
40
50
20
0
0
−50
0
20
40
60
80
100 Lag m
120
140
160
180
200
−20
0
20
40
60
80
100 Lag m
120
140
160
180
200
Figura 23: Gráfico de autocorrelación a la izquierda, y a la derecha la correlación cruzada para σ 2 = 1. La SNR en el caso es de: SN R(y) = 1.8841[dB] Para un σ 2 = 1
(39)
Con respecto a todo lo mostrado anteriormente, se observa que al igual que ante, en las correlaciones cruzadas el máximo se encontraba en D = 20, y no sólo eso, puesto que las SNR obtenidas fueron mayores respecto a las SNR obtenidas para las mismas varianzas pero con una secuencia de Barker no repetida. La explicación se debe a que en los casos en realidad lo que se hizo fue hacer exhaustivamente el mismo experimento de antes de la única secuencia de Barker, pero en este caso se hacía repetidas veces (10 precisamente), lo cual elimina la aleatoriedad presente en el experimento en el caso original. Esto permite hacer un mejor cálculo de D, pues ahora representa una extrapolación del análisis presentado con solo una secuencia de Barker, y lo cual explica los máximos locales factores del retardo D presentes en los gráficos de las correlaciones cruzadas. Parte F: Se busca encontrar la relación matemática entre rxx y ryy . Para esto primero se revisa la ecuación (27) donde se tiene que: ryy =
∞ X
∞ X
a2 x[n − D]x[n − D + m] +
ax[n − D]v[n + m] +
v[n]ax[n − D + m] + rvv
(40)
n=−∞
n=−∞
n=−∞
∞ X
Si se hace un cambio de variable i = n − D en la primera sumatoria y segunda sumatoria se tiene entonces que: ryy =
∞ X
a2 x[i]x[i − m] +
i=−∞
∞ X
∞ X
ax[i]v[i + D − m] +
i=−∞
av[n]x[n − D − m] + rvv
(41)
n=−∞
Entonces se tiene que: ryy = a2 rxx + arxv [m − D] + arvx [m + D] + rvv
(42) 2
Por lo tanto, se puede ver la relación entre ryy y rxx , siendo el último multiplicado por a y con sumandos aleatorios. Por lo tanto si se considera el ruido de baja potencia (o varianza), entonces ryy será muy similar pero escalado con un valor de a2 . Se hacen pruebas para ruido de varianza de σ 2 = 0.01 y un a = 2, y por entrada la secuencia Barker sin repetición, obteniendo los siguientes resultados de gráficos: Correlación ryy[m]
Correlación cruzada rxx[m]
60
14
X= 0 Y= 53.5093
50
40
10
30
8
20
6
10
4
0
2
−10
0
20
X= 0 Y= 13
12
40
60
80
100 Lag m
120
140
160
180
200
0
0
20
40
60
80
100 Lag m
120
140
160
180
200
Figura 24: Con a = 2 y σ 2 = 0.01 gráfico de autocorrelaciones, a la izquierda de ry y[m] y a la derecha de rxx [m] Se ve claramente que son muy parecidas, sólo el ruido presente en la autocorrelación de la salida, y a su vez se ve en el máximo un escalamiento, siendo el máximo de la autocorrelación de la salida 53, y el de la entrada 13, considerando entonces que está presente la amplificación cercana a a2 = 4, por lo que el análisis teórico propuesto es válido. UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
17
Anexo: Códigos Usados en MATLAB Sección I: Evaluación de propiedades de sistemas.
1 2 3 4 5 6
t=0:1e−6:1; x1=wgn(1,length(t),0); x2=wgn(1,length(t),0); a=2; b=10; x3=a*x1+b*x2;
7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
%Linealidad %bbox1 y1=bbox1(x1); y2=bbox1(x2); y3=bbox1(x3); figure plot(t,a*y1+b*y2); hold on plot(t,y3,'r'); hold off h = legend('$ay_1+by_2$','$y_3$'); set(h,'Interpreter','latex'); grid on; ylabel('Amplitud'); xlabel('Tiempo [s]'); title('Respuestas de bbox1() a ruido blanco');
24 25 26 27 28 29 30 31 32 33 34 35 36
figure plot(t(1,40:60),a*y1(1,40:60)+b*y2(1,40:60)),'b−−o'; hold on plot(t(1,40:60),y3(1,40:60),'r−−o'); hold off h = legend('$ay_1+by_2$','$y_3$'); set(h,'Interpreter','latex'); grid on; ylabel('Amplitud'); xlabel('Tiempo [s]'); xlim([t(1,40) t(1,60)]); title('Respuestas con zoom de bbox1() a ruido blanco');
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52
%bbox2 y1=bbox2(x1); y2=bbox2(x2); y3=bbox2(x3); figure plot(t,a*y1+b*y2); hold on plot(t,y3,'r'); hold off h = legend('$ay_1+by_2$','$y_3$'); set(h,'Interpreter','latex'); grid on; ylabel('Amplitud'); xlabel('Tiempo [s]'); title('Respuestas de bbox2() a ruido blanco');
53 54 55 56 57 58 59 60 61 62 63 64 65
figure plot(t(1,40:60),a*y1(1,40:60)+b*y2(1,40:60),'b−−o'); hold on plot(t(1,40:60),y3(1,40:60),'r−−o'); hold off h = legend('$ay_1+by_2$','$y_3$'); set(h,'Interpreter','latex'); grid on; ylabel('Amplitud'); xlabel('Tiempo [s]'); xlim([t(1,40) t(1,60)]); title('Respuestas con zoom de bbox2() a ruido blanco');
66
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
18
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81
%bbox3 y1=bbox3(x1); y2=bbox3(x2); y3=bbox3(x3); figure plot(t,a*y1+b*y2); hold on plot(t,y3,'r'); hold off h = legend('$ay_1+by_2$','$y_3$'); set(h,'Interpreter','latex'); grid on; ylabel('Amplitud'); xlabel('Tiempo [s]'); title('Respuestas de bbox3() a ruido blanco');
82 83 84 85 86 87 88 89 90 91 92 93 94
figure plot(t(1,40:60),a*y1(1,40:60)+b*y2(1,40:60),'b−−o'); hold on plot(t(1,40:60),y3(1,40:60),'r−−o'); hold off h = legend('$ay_1+by_2$','$y_3$'); set(h,'Interpreter','latex'); grid on; ylabel('Amplitud'); xlabel('Tiempo [s]'); xlim([t(1,40) t(1,60)]); title('Respuestas con zoom de bbox3() a ruido blanco');
95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
%Invarianza x1=sin(2*pi*30*t); x2=[zeros(1,200000),sin(2*pi*30*t(1,200001:1000001))]; y1=bbox1(x1); y2=bbox2(x1); y3=bbox3(x1); %bbox1 figure plot(t,[zeros(1,200000),y1(1,1:800001)]); hold on plot(t,bbox1(x2),'r'); hold off h = legend('$y(t−0.2)$','$bbox1(x(t−0.2))$'); set(h,'Interpreter','latex'); grid on; ylabel('Amplitud'); xlabel('Tiempo [s]'); title('Prueba de invarianza para bbox1');
114 115 116 117 118 119 120 121 122 123 124 125 126 127
%bbox2 figure plot(t,[zeros(1,200000),y2(1,1:800001)]); hold on plot(t,bbox2(x2),'r'); hold off h = legend('$y(t−0.2)$','$bbox2(x(t−0.2))$'); set(h,'Interpreter','latex'); grid on; ylabel('Amplitud'); xlabel('Tiempo [s]'); title('Prueba de invarianza para bbox2'); xlim([0.1 0.9]);
128 129 130 131 132 133 134 135 136 137
%bbox3 figure plot(t,[zeros(1,200000),y3(1,1:800001)]); hold on plot(t,bbox3(x2),'r'); hold off h = legend('$y(t−0.2)$','$bbox3(x(t−0.2))$'); set(h,'Interpreter','latex'); grid on;
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
19
138 139 140
ylabel('Amplitud'); xlabel('Tiempo [s]'); title('Prueba de invarianza para bbox3');
141 142 143 144 145 146
%Respuesta a impulso impulso=[zeros(1,29),1,zeros(1,30)]; y1imp=bbox1(impulso); y2imp=bbox2(impulso); y3imp=bbox3(impulso);
147 148 149 150 151 152 153
figure stem(y1imp); grid on ylabel('Amplitud'); xlabel('Muestra'); title('Respuestas de bbox1() a impulso');
154 155 156 157 158 159 160
figure stem(y2imp); grid on ylabel('Amplitud'); xlabel('Muestra'); title('Respuestas de bbox2() a impulso');
161 162 163 164 165 166 167
figure stem(y3imp); grid on ylabel('Amplitud'); xlabel('Muestra'); title('Respuestas de bbox3() a impulso');
Sección II: Sistemas inversos. Para los sistemas SA y SB se crearon las siguientes funciones: 1 2 3 4 5
1 2 3 4 5
function y = S_A(x) y(1)=−x(1) for n= 2:length(x) y(n)=0.5*y(n−1)−x(n); end
function y = S_B(x) y(1)=−x(1) for n= 2:length(x) y(n)=0.5*x(n−1)−x(n); end
Generación de gráficos y otros: 1 2 3 4 5 6 7 8
imp=[1,zeros(1,29)]; figure stem(S_A(imp)); grid on xlabel('Muestra'); ylabel('Amplitud'); h = title('Respuesta a impulso del sistema $S_A(x)$'); set(h,'Interpreter','latex');
9 10 11 12 13 14 15 16
figure stem(S_B(imp)); grid on xlabel('Muestra'); ylabel('Amplitud'); h = title('Respuesta a impulso del sistema $S_B(x)$'); set(h,'Interpreter','latex');
17
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
20
18 19 20 21 22 23 24
figure stem(S_B(S_A(imp))); grid on xlabel('Muestra'); ylabel('Amplitud'); h = title('Respuesta a impulso del sistema $S_B(S_A(x))$'); set(h,'Interpreter','latex');
Sección III: Derivadas e integrales discretas. Función que implemente el sistema para el derivador en tiempo discreto: 1 2 3 4 5
function y = ddis(x,T) y(1)=x(1); for n= 2:length(x) y(n)=(1/T)*(x(n)−x(n−1)); end
Función que implemente el sistema para el integrador en tiempo discreto: 1 2 3 4 5
function y = idis(x,T) y(1)=x(1); for n= 2:length(x) y(n)=y(n−1)+ T*x(n); end
Generación de gráficos: 1 2 3 4 5 6 7 8 9 10 11
T=1e−4; t=0:T:1; x=sin(200*pi*t); %Parte b figure subplot(1,3,1); plot(t(1:500),x(1:500)); grid on title('Sinusoide de 100[Hz] x(t)'); xlabel('Tiempo [s]'); ylabel('Amplitud');
12 13 14 15 16 17 18
subplot(1,3,2); plot(t(1:500),ddis(x(1:500),T)); grid on title('Derivada discretada de x(t)'); xlabel('Tiempo [s]'); ylabel('Amplitud');
19 20 21 22 23 24 25
subplot(1,3,3); plot(t(1:500),idis(x(1:500),T)); grid on title('Integración discreta de x(t)'); xlabel('Tiempo [s]'); ylabel('Amplitud');
26 27 28 29 30 31 32 33 34 35 36 37 38 39
%Parte c T=1; t=−10:T:20; x=zeros(1,length(t)); x(1,11)=1; x(1,16)=−1; figure subplot(1,3,1); stem(t,x); grid on h = title('Entrada $x[n]=\delta[n] − \delta[n−5]$'); set(h,'Interpreter','latex'); xlabel('Muestra');
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
21
40
ylabel('Amplitud');
41 42 43 44 45 46 47
subplot(1,3,2); stem(t,ddis(x,T)); grid on title('Derivada discretada de x(t)'); xlabel('Muestra'); ylabel('Amplitud');
48 49 50 51 52 53 54
subplot(1,3,3); stem(t,idis(x,T)); grid on title('Integración discreta de x(t)'); xlabel('Muestra'); ylabel('Amplitud');
55 56 57 58 59 60 61 62 63 64 65 66 67 68
%parte D T=1; t=−10:T:20; x=zeros(1,length(t)); x(1,11:21)=1; figure subplot(1,3,1); stem(t,x); grid on h = title('Entrada $x[n]=\mu[n] − \mu[n−11]$'); set(h,'Interpreter','latex'); xlabel('Muestra'); ylabel('Amplitud');
69 70 71 72 73 74 75
subplot(1,3,2); stem(t,ddis(x,T)); grid on title('Derivada discretada de x(t)'); xlabel('Muestra'); ylabel('Amplitud');
76 77 78 79 80 81 82
subplot(1,3,3); stem(t,idis(x,T)); grid on title('Integración discreta de x(t)'); xlabel('Muestra'); ylabel('Amplitud');
83 84 85 86 87 88 89 90
%parte E [x1,Fs1]=audioread('music.wav'); [x2,Fs2]=audioread('speech.wav'); y11=ddis(x1,1/Fs1); y12=idis(x1,1/Fs1); y21=ddis(x2,1/Fs2); y22=idis(x2,1/Fs2);
Sección IV: Retardo en radares. Función que implementa el sistema del radar discreto: 1 2 3 4 5
function y=sis_rad(x,v,D,a) y(1:D)=v(1,1:20); for n=21:length(x) y(n)=a*x(n−D)+v(1,n); end
Función para la correlación cruzada: 1 2 3 4
function rxy=corr_cruz(x,y) rxy=zeros(1,length(x)); for m=1:length(x) for n=1:length(x)
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
22
if n+m−1<=length(x) rxy(m)=rxy(m)+y(n)*x(n+m−1); end
5 6 7
end
8 9
end
Función para obtener la SNR de dos señales: 1 2 3 4 5 6 7 8 9 10
function val=snrdb(senal,ruido) potsen=0; potruido=0; for n=1:length(senal) potsen=potsen+(abs(senal(n)))^2; end for n=1:length(ruido) potruido=potruido+(abs(ruido(n)))^2; end val=10*log10(potsen/potruido);
Función para crear los gráficos de esta sección: 1 2 3 4 5 6 7 8 9 10 11 12
function y=parte4graf(x,y,rxy,ryy,muestra) figure subplot(1,2,1) stem(muestra,x); title('Entrada x[n] sec. Barker'); xlabel('Muestras'); grid on subplot(1,2,2) stem(muestra,y); title('Salida y[n] del sist.con entrada x[n]'); xlabel('Muestras'); grid on
13 14 15 16 17 18 19 20 21 22 23 24
figure subplot(1,2,1) stem(muestra,ryy); xlabel('Lag m'); grid on title('Correlación ryy[m]'); subplot(1,2,2) stem(muestra,rxy); xlabel('Lag m'); title('Correlación cruzada rxy[m]'); grid on
Script para generar los gráficos necesarios de cada sección: 1 2 3 4 5 6 7 8 9 10 11
muestra=0:1:199; D=20; a=0.9; x=[1 1 1 1 1 −1 −1 1 1 −1 1 −1 1,zeros(1,187)]; %parte b v=sqrt(0.01)*randn(1,200); y=sis_rad(x,v,D,a); ryy=corr_cruz(y,y); rxy=corr_cruz(y,x); parte4graf(x,y,rxy,ryy,muestra); snr4c=snrdb(y,v)
12 13 14 15 16 17 18 19
%parte d sigma^2=0.1 v=sqrt(0.1)*randn(1,200); y=sis_rad(x,v,D,a); ryy=corr_cruz(y,y); rxy=corr_cruz(y,x); parte4graf(x,y,rxy,ryy,muestra); snr4d1=snrdb(y,v)
20
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones
23
21 22 23 24 25 26 27
%parte d sigma^2=1 v=sqrt(1)*randn(1,200); y=sis_rad(x,v,D,a); ryy=corr_cruz(y,y); rxy=corr_cruz(y,x); parte4graf(x,y,rxy,ryy,muestra); snr4d2=snrdb(y,v)
28 29 30 31 32 33 34 35 36 37 38
%parte e x=[x(1,1:13) x(1,1:13)]; x=[x(1,1:26) x(1,1:26) x(1,1:26) x(1,1:26) x(1,1:26),zeros(1,70)]; %parte e sigma^2=0.01 v=sqrt(0.01)*randn(1,200); y=sis_rad(x,v,D,a); ryy=corr_cruz(y,y); rxy=corr_cruz(y,x); parte4graf(x,y,rxy,ryy,muestra); snr4e=snrdb(y,v)
39 40 41 42 43 44 45 46
%parte e sigma^2=0.1 v=sqrt(0.1)*randn(1,200); y=sis_rad(x,v,D,a); ryy=corr_cruz(y,y); rxy=corr_cruz(y,x); parte4graf(x,y,rxy,ryy,muestra); snr4e1=snrdb(y,v)
47 48 49 50 51 52 53 54
%parte e sigma^2=1 v=sqrt(1)*randn(1,200); y=sis_rad(x,v,D,a); ryy=corr_cruz(y,y); rxy=corr_cruz(y,x); parte4graf(x,y,rxy,ryy,muestra); snr4e2=snrdb(y,v)
UTFSM - ELO-313 - 2016
Tarea #1 Procesamiento digital de señales con aplicaciones