Publicado en por Edison Del Rosario

4.3 FIR – Respuesta de frecuencia y superposición

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]
..

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) = 3

y 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.

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]
..

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

# 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)

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]
..

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.

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]
..

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 t

Aplicando 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)

[ ejercicio bk ] [ algoritmo_bk ] [ ejercicio_xn ] [ algoritmo_xn ]