Para sistemas lineales invariantes en el tiempo, LTI, las entradas senoidales complejas tienen como respuesta otra señal senoidal de la misma frecuencia pero con diferente magnitud y fase. TODOS los sistemas LTI tienen la propiedad "entrada-sinusoide" genera una "salida-sinusoide".
Para un filtro FIR LTI, ante una entrada discreta en tiempo tipo señal exponencial compleja, la salida también es una señal discreta en tiempo, tipo exponencial compleja con diferente amplitud, de la misma frecuencia.
|H(e^{j \omega})| = \sum _{k=0}^{M} b_k e^{-j\omega k} = \sum _{k=0}^{M} h[k] e^{-j\omega k} y[n] = |H(e^{j \omega})| x[n] \text{ ; para } x[n] = A e^{j\varphi}e^{-j\omega n}1. Ejercicio bk
Referencia: McClellan Ejemplo 6.1 y 6.2 p217
Considere un sistema LTI con coeficientes de la ecuación de diferencias bk = [1,2,1]. Al convertirlo en la función respuesta de frecuencias del sistema H(ω) se tiene:
H(e^{j \omega}) = 1 + 2 e^{-j\omega} + e^{-j \omega 2}como bk es simétrico, se obtiene como factor común e(-jω) para simplificar la expresió con un coseno.
= e^{-j \omega} \Big( e^{j \omega} + 2 + e^{-j \omega} \Big) = e^{-j \omega} \Big( 2 +2 \cos(\omega) \Big)donde la magnitud es:
|H(e^{j \omega})|= \Big( 2 +2 \cos(\omega) \Big)y la fase:
\angle H(e^{j \omega})= e^{-j \omega}Al disponer de una señal de entrada con una frecuencia de ω=π/3:
x[n] = 2 e^{j \pi /4} e^{j(\pi /3)n}se evalúa la magnitud de |H|:
|H(e^{j \pi /3})|= \Big( 2 +2 \cos(\pi /3) \Big) = 3y la fase de ∠H:
\angle H(e^{j \pi /3 })= e^{-j \pi/3}siendo y[n] = (H)(x[n])
y[n] = 3 e^{-j\pi /3}\Big( 2 e^{j \pi /4} e^{j(\pi /3)n} \Big) = (3)(2) e^{j \pi /4-j\pi /3} e^{j\pi /3 n} = 6 e^{-j \pi /12} e^{j\pi /3 n} = 6 e^{j \pi /4} e^{j\pi(n-1) /3}la salida del sistema es igual a la entrada, multiplicada por 3 y un cambio de fase a -π/3 que es un retraso de una muestra. Recordando que esto aplica cuando bk es simétrico y la entrada es una exponencial compleja.
2. Algoritmo bk en Python
Primero se analiza la simetría del vector bk, invirtiendo el vector y restando de su forma original, si todos los resulados son ceros, el vector es simétrico.
Acorde con la simetría de bk se procede a la creación de H y luego yn a partir de xn.
Para xn en el algoritmo, se requiere usar la dsp.euler_args_one_term() para dominio del tiempo contínuo t de la Unidad 1.4 Espectro – Suma de sinusoides y la fórmula de Euler.
La función en telg1034.py amplía el uso al dominio de tiempo discreto n al sustituir la variable independiente con las siguientes instrucciones y se obtiene el resultado esperado.
# revisa si unasenal es discreta con n, Unidad 4 Filtros Digitales
esdiscreta = False
varindepend = X.free_symbols
if (n in varindepend) and not(t in varindepend):
X = X.subs(n,t)
esdiscreta = True
# analiza unasenal con t
El resultado con el algoritmo es:
bk_simetrico: True
H_magn : exp(I*w) + 2 + exp(-I*w)
H_magcos: 2*cos(w) + 2
H_fase : -w
H : (2*cos(w) + 2)*exp(-I*w)
x[n]: 2*exp(I*pi/4)*exp(I*pi*n/3)
wk : pi/3
y[n]: 6*exp(I*pi/4)*exp(I*pi*(n - 1)/3)
>>>Instrucciones en Python
Algoritmos y funciones: telg1034.py
# ejemplo 6.1 p217 FIR - Respuesta de frecuencias
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1034 as dsp
# variables continuas
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes
# variables discretas
from telg1034 import n
# INGRESO
bk = [1,2,1]
#bk = [1,-2,4,-2,1]
xn = 2*sym.exp(I*pi/4)*sym.exp(I*(pi/3)*n)
# PROCEDIMIENTO
M = len(bk)-1
# revisa simetria bk
bk = np.array(bk)
H_factor = 0
sumainvertida = np.sum(np.abs(bk-np.flip(bk)))
bk_simetrico = False
if sumainvertida == 0 and M>1:
bk_simetrico = True
H_factor = sym.sympify(M)/2
# respuesta de frecuencia
H = sym.S.Zero
for k in range(0,M+1,1):
H = H + bk[k]*sym.exp(-I*w*(k))
# H cuando bk es simetrico.
H_magn = H ; H_magcos = 0 ; H_fase = 0
if bk_simetrico:
H_magn = sym.expand(H/sym.exp(-I*w*H_factor))
H_magcos = H_magn.rewrite(sym.cos)
H_fase = -w*H_factor
H = H_magcos*sym.exp(I*H_fase)
# señal de entrada x[n]
xn_arg = dsp.euler_args_one_term(xn)
wk = xn_arg['freq_rad']
# señal de salida y[n]
if bk_simetrico:
n_veces = H_fase/w
Ampl = H_magcos.subs(w,wk)*xn_arg['amplitud']
yn = Ampl*sym.exp(I*xn_arg['fase'])*sym.exp(I*wk*(n+n_veces))
else:
yn = xn*H.subs(w,wk)
# SALIDA
print('bk_simetrico:',bk_simetrico)
print('H_magn :',H_magn)
print('H_magcos:',H_magcos)
print('H_fase :',H_fase)
print(' H :',H)
print('x[n]:',xn)
print('wk :',wk)
print('y[n]:',yn)
3. Ejercicio x[n] de varias sinusoides
Referencia: McClellan Ejemplo 6.4 p220
Para el filtro FIR con coeficientes b[k] = [1,2,1] , encuentre la señal de salida cuando la señal de entrada tiene varios componentes.
x[n] = 4 + 3 \cos \Big( \frac{\pi}{3} n - \frac{pi}{2} \Big) + 3 \cos \Big(\frac{7\pi}{8}n \Big)El ejercicio se desarrolla término a término por tener diferentes frecuencias. Se aplica el principio de superposición.
4. Algoritmo X[n] en Python
Para el ejercicio, xn tiene una señal senoidal mas general y de varios componentes. Por lo que el procedimiento del algoritmo anterior se amplía para cada uno de los componentes.
La señal de entrada se debe revisar sus componentes sumas, procediendo de forma semejante a lo realizado en la Unidad 1.2 Señal – Periodo y Desplazamiento en tiempo. La función dsp.x_list_term_Add() para separar los términos suma entrega su resultado tanto para la variable t ó n.
Sin embargo la obtención de los parámetros de cada término suma de xn con la función dsp.cos_args_one_term(), debe ampliar la función de t a n, semejante al ejercicio que precede con las siguientes instrucciones.
# revisa si unasenal es discreta con n, Unidad 4. Filtros digitales
esdiscreta = False
varindepend = x.free_symbols
if (n in varindepend) and not(t in varindepend):
x = x.subs(n,t)
esdiscreta = True
# analiza unasenal con tAplicando H en cada término xn acorde a la simetría de bk, se tiene el siguiente resultado con el algoritmo es:
H_magn : exp(I*w) + 2 + exp(-I*w)
H_magcos: 2*cos(w) + 2
H_fase : -w
H : (2*cos(w) + 2)*exp(-I*w)
x[n]: 3*cos(7*pi*n/8) + 3*cos(pi*n/3 - pi/2) + 4
y[n]: 0.456722804932279*cos(pi*(7*n/8 - 7/8)) + 9*cos(pi*(n - 1)/3 - pi/2) + 16
>>>Instrucciones en Python
# ejemplo 6.4 p220 FIR - Respuesta de frecuencias a varias senoidales
# telg1034 DSP fiec-espol edelros@espol.edu.ec
import numpy as np
import matplotlib.pyplot as plt
import sympy as sym
import telg1034 as dsp
# variables continuas
from telg1034 import t,A,w,f,p,pi,DosPi,I,equivalentes
# variables discretas
from telg1034 import n
# INGRESO
bk = [1,2,1]
#bk = [1,-2,4,-2,1]
#xn = 3*sym.cos(n*pi/3 - pi/2)
xn = 4 + 3*sym.cos(n*pi/3 - pi/2) + 3*sym.cos(7*pi*n/8)
# PROCEDIMIENTO
M = len(bk)-1
# revisa simetria bk
bk = np.array(bk)
H_factor = 0
sumainvertida = np.sum(np.abs(bk-np.flip(bk)))
bk_simetrico = False
if sumainvertida == 0 and M>1:
bk_simetrico = True
H_factor = sym.sympify(M)/2
# respuesta de frecuencia
H = sym.S.Zero
for k in range(0,M+1,1):
H = H + bk[k]*sym.exp(-I*w*(k))
# H cuando bk es simetrico.
H_magn = H ; H_magcos = 0 ; H_fase = 0
if bk_simetrico:
H_magn = sym.expand(H/sym.exp(-I*w*H_factor))
H_magcos = H_magn.rewrite(sym.cos)
H_fase = -w*H_factor
H = H_magcos*sym.exp(I*H_fase)
# señal de entrada x[n]
x_senal = dsp.x_list_term_Add(xn) # términos suma
x_conteo = len(x_senal)
yn = sym.S.Zero
for i in range(0,x_conteo,1):
unasenal = x_senal[i]
cond1 = unasenal.has(sym.cos) or unasenal.has(sym.sin)
cond2 = not(unasenal.has(n))
if cond1 or cond2:
xn_arg = dsp.cos_args_one_term(unasenal)
Ak = xn_arg['amplitud']
wk = xn_arg['freq_rad']
fk = xn_arg['fase']
if bk_simetrico:
H_magnk = dsp._float_is_int(H_magcos.subs(w,wk).evalf())
H_fasek = H_fase.subs(w,wk)
if wk!=sym.S.Zero:
n_veces = H_fasek/wk
yk = H_magnk*Ak*sym.cos(wk*(n+n_veces) + fk)
else:
yk = H_magnk*Ak
else:
yk = H.subs(w,wk)*unasenal
else: # no es senoidal o constante
yk = H*unasenal
print(' ____revisar término:',unasenal.subs(t,n))
# principio de superposición
yn = yn +yk
# SALIDA
print('H_magn :',H_magn)
print('H_magcos:',H_magcos)
print('H_fase :',H_fase)
print(' H :',H)
print('x[n]:',xn)
print('y[n]:',yn)