Filtros FIR Recursivos e Interpolados
Teoría y Aplicaciones a la electrocardiografía
Por Mariano Llamedo Soria
Resumen
En este notebook se presenta un filtro FIR de fase lineal tipo peine y recursivo, que a partir de la interpolación de su respuesta al impulso resulta muy útil para el procesamiento de señales biológicas como el electrocardiograma (ECG). La respuesta en frecuencia tipo peine posee ceros de transmisión tanto en continua (DC ó 0 Hz) como en los múltiplos enteros de la frecuencia de línea eléctrica (50/60 Hz). El cero en DC es fundamental para las interferencias de muy baja frecuencia (<1 Hz), mientras que el resto mitigará cualquier acoplamiento de señal de línea. En este documento se analizan varios aspectos teóricos relevantes del filtro y se finaliza analizando una implementación y uso en ECG real proveniente de una prueba de esfuerzo. Se trata de un registro con notorias interferencias de baja frecuencia y de línea eléctrica de 50 Hz.
Funciones de PyTC2 utilizadas:
Análisis de la respuesta en frecuencia: plot_plantilla, analyze_sys, group_delay
Funciones de diseño digital fir_design_pm, herrmann_lp_fir_order, DC_PWL_removal_recursive_filter
De presentación algebraica: print_latex, a_equal_b_latex_s, print_subtitle, simplify_n_monic_z, expr_simb_expr
Parte I: Análisis Teórico
En la primera parte se realizará un análisis teórico de los conceptos teóricos que fundamentan:
la comprensión de los filtros FIR de fase lineal,
la implementación recursiva de algunos filtros de fase lineal, y
el sobremuestreo o interpolación de la respuesta al impulso de un filtro.
Filtros FIR de fase lineal
En los notebooks anteriores revisamos los filtros de respuesta finita, o FIR, y sus principales virtudes:
Estabilidad, ya que todos sus polos se encuentran en el origen.
Fase lineal, debido a la simetría de la respuesta al impulso.
Retardo constante, como consecuencia de la linealidad de fase.
Tales ventajas van en detrimento del costo computacional por muestra de salida calculada. Es decir, de la cantidad de multiplicaciones y sumas, como también de la cantidad de unidades de memoria necesarias para su implementación. En este notebook se evaluará la posibilidad de implementar ciertos FIR de manera recursiva, sin afectar los beneficios antedichos. Como se verá, el costo computacional por muestra de salida disminuirá considerablemente.
Recordando que la respuesta en frecuencia de un FIR de respuesta simétrica (tipos I o II) se expresan de forma polar como:
donde la fase es
siendo \(M\) el eje de simetría temporal de la respuesta al impulso, coincidiendo con la muestra media para un filtro de \(N\) coeficientes, y sobretodo, con el retardo del filtro constante para toda frecuencia. Luego su respuesta de fase cero será
estando relacionada con la respuesta de módulo mediante \( \left| H(\Omega) \right| = \left| H_R(\Omega) \right| \).
Luego \(Q(\Omega)\) será una función par para estos dos tipos de FIR, dependiendo del largo \(nn\), como se ve en la siguiente tabla.
Tipo |
\( Q(\Omega) \) |
\(H_R(0)\) |
\(H_R(\pi)\) |
Transferencias posibles |
|---|---|---|---|---|
I |
1 |
\(A(\Omega)\) |
\(A(\Omega)\) |
LP-H-BP |
II |
\( \cos{(\Omega/2)} \) |
\(A(\Omega)\) |
0 |
LP-BP |
Vemos que el tipo I será el más versátil, ya que permite diseñar cualquier respuesta, mientras que el tipo II, al forzar un 0 en Nyquist restringe las posibilidades.
En principio se rediseñará un filtro FIR promediador (pasabajos), para convertirlo en un filtro FIR recursivo. Luego se usará para resolver un problema concreto: diseñar un FIR pasa-altos con un cero en \(\Omega=0\) (DC). Por último, mediante el concepto de interpolación de la respuesta al impulso, se convertirá en un filtro peine particularmente útil para el preprocesamiento del ECG. Se presentará la comparación con otros FIR diseñados mediante el resto de metodologías vistas hasta ahora.
Un FIR de fase lineal para remover continua
Un problema fácil de explicar, pero difícil de lograr como tantos en la Ingeniería. En esta sección se retoma un artículo de 2008 de Rick Lyons en dsprelated.com donde trató este tema desde una perspectiva muy interesante, ya que convierte al filtro FIR promediador de \(N\) muestras
cuya transferencia es
en una versión recursiva de la misma transferencia
Si bien a primera vista cuesta entender la equivalencia entre ambas expresiones, conviene pensar el cálculo del promedio entre dos muestras consecutivas. La sumatoria en \(y(n)\) e \(y(n-1)\) tendrá \(N-2\) términos comunes, salvo el primer y último término, es decir \(x(n)\) e \(x(n-N)\). Luego aplicando transformada z llegamos a la consecuente transferencia recursiva
Del mismo modo uno pensaría que esta transferencia no es estable, salvo porque se da una cancelación polo-cero para \(\Omega = 0\) para todo N, resultando \(T(z) = T_R(z)\). De esta forma, un promedio de N muestras pasa de costar N sumas a tan solo 2 sumas, para todo N, como se ve en el siguiente esquema
\(T_R(z)\) - Filtro promediador recursivo ó CIC
Esto se ve claramente en los siguientes ejemplos
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt
from pytc2.sistemas_lineales import group_delay, analyze_sys, plot_plantilla
from pytc2.general import print_latex, print_subtitle, a_equal_b_latex_s, simplify_n_monic_z, expr_simb_expr
import warnings
fig_font_size = 12
plt.rcParams.update({'font.size':fig_font_size})
nn = 5
fs = 2
# Promediador de orden nn
# TransferFunction los asume en potencias decrecientes
num = 1/(nn)* np.ones(nn)
den = np.zeros(nn)
den[0] = 1
my_df = sig.TransferFunction(num, den, dt=1/fs)
# Promediador de orden N recursivo
# TransferFunction los asume en potencias decrecientes
num_R = 1/nn*np.concatenate([[1], np.zeros(nn-1),[-1]])
den_R = np.zeros(nn+1)
den_R[:2] = np.array([1, -1])
my_df_R = sig.TransferFunction(num_R, den_R, dt=1/fs)
plt.close('all')
# Visualizamos su respuesta en frecuencia
figaxes = analyze_sys([my_df, my_df_R], ['típico', 'recursivo'], title_suffix = ': FIR Promediador', annotations=False, fs=fs)
Este filtro pasabajos recursivo, también se lo conoce en otros ámbitos como cascaded integrator comb (CIC), como sugiere su respuesta de módulo.
Filtros de Complementariedad de Demora
Lo interesante de esta transferencia, es que se la puede utilizar para lograr un pasa-altos, mediante el concepto de complementaridad de demora, es decir
esto significa que para un filtro \(T(z)\) de \(N\) coeficientes, cuya demora es \(\frac{N-1}{2}\) ó \(M\) como se definió más arriba para el tipo I. Entonces los módulos de ambos filtros deben respetar que
Es decir, que se complementen en módulo. Luego, el filtro que buscamos surge como el complementario a nuestro pasabajos recursivo, es decir
Según comenta Lyons, conviene que N sea par para que la implementación del filtro sea aún más eficiente, siendo el producto \(1/N\) una mera rotación de bits. Esto conlleva que la demora \(M\) no sea entera, sin embargo, puede solventarse cascadeando una cantidad par (\(C\)) de CIC’s, es decir
entonces para un caso sencillo como \(C = 2\), se llega a
Cuyo diagrama en bloques se ve a continuación
\(T_C(z)\) para \(C = 2\).
Que de forma más explícita se puede expresar como
Luego, reordenando y reutilizando algunas estructuras en pos de la eficiencia computacional se puede llegar a la estructura que propone Lyons en otro artículo en dsprelated
\(T_C(z)\) para \(C = 2\). Modificado de la Fig. 1 en dsprelated
Como se puede anticipar, generalizarlo para otros valores de \(C\) da lugar a coeficientes más espaciados, producto de aumentar la demora. Es decir que a medida que escala \(N\) y \(C\), los coeficientes del numerador y denominador se hacen más dispersos, es decir aumenta la proporción de coeficientes nulos. Este hecho es muy significativo para la eficiencia computacional, y cuando se interpole esta respuesta, se ampliará aún más la escasez o proporción de coeficientes nulos (sparsity en inglés).
Por el momento se realizará un análisis numérico de la respuesta de \(T_2(z)\) para \(N = 5\).
from pytc2.sistemas_lineales import pzmap
nn = 5
fs = 2
# T_2(z)
#
num = np.zeros(2*nn+1)
num[0] = -1
num[nn-1] = nn**2
num[nn] = 2-2* nn**2
num[nn+1] = nn**2
num[2*nn] = -1
num = (1/nn**2) * num
den = np.array([1, -2, 1])
w, h_iir = sig.freqz(num, den, fs = fs)
mag = 20.0 * np.log10(np.abs(h_iir) + np.nanmin(np.abs(h_iir)))
# resguardo el cero en DC para que funcione unwrap
phase = np.angle(h_iir)
phase = np.unwrap(phase[1:])
this_lbl = '$T_2(z)$'
plt.figure(figsize=(12, 5))
plt.plot(w, mag, label= this_lbl )
plt.title(f'Respuesta en Frecuencia de Módulo')
plt.ylabel('Magnitud (dB)')
ax = plt.gca()
plt.legend()
plt.grid(True)
plt.figure(figsize=(12, 5))
plt.plot(w[1:], phase, label= this_lbl) # Bode phase plot
plt.title(f'Respuesta en Frecuencia de Fase')
plt.ylabel('Fase (rad)')
ax = plt.gca()
plt.legend()
plt.grid(True)
plt.figure(figsize=(12, 5))
gd_fir = group_delay(w[1:] * np.pi, phase)
# Para órdenes grandes
plt.plot(w[2:], gd_fir[1:], label=this_lbl ) # Bode phase plot
plt.ylim((0,5))
plt.title(f'Respuesta en Frecuencia de Retardo')
plt.ylabel('Retardo (s)')
ax = plt.gca()
plt.legend()
plt.grid(True)
/home/mariano/pytc2_dev/lib/python3.12/site-packages/scipy/signal/_filter_design.py:482: RuntimeWarning: divide by zero encountered in divide
h = (npp_polyval(zm1, b, tensor=False) /
/home/mariano/pytc2_dev/lib/python3.12/site-packages/scipy/signal/_filter_design.py:482: RuntimeWarning: invalid value encountered in divide
h = (npp_polyval(zm1, b, tensor=False) /
Del análisis de la respuesta en frecuencia se nota claramente la influencia del cero en DC en la respuesta de módulo, y las mismas ventajas de la fase lineal que se analizaron anteriormente.
Para enfatizar esto, visualizamos el diagrama de polos y ceros
my_df = sig.TransferFunction(num, den, dt=1/fs)
# Visualizamos su respuesta en frecuencia
_ = pzmap(my_df, digital=True, fs= fs)
Lamentablemente no se observa con claridad el factor de multiplicidad de ceros y polos. Recurrimos a un análisis pormenorizado.
roots = np.roots(num)
vals = np.hstack((
np.abs(roots).reshape(-1,1),
np.angle(roots).reshape(-1,1)
))
print_subtitle('Ceros')
for mag, ang in vals:
print(f"{mag:0.3f} {ang:0.3f}")
roots = np.roots(den)
vals = np.hstack((
np.abs(roots).reshape(-1,1),
np.angle(roots).reshape(-1,1)
))
print_subtitle('Polos')
for mag, ang in vals:
print(f"{mag:0.3f} {ang:0.3f}")
Ceros
2.618 3.142
2.387 1.750
2.387 -1.750
1.000 0.000
1.000 0.000
1.000 -0.000
1.000 0.000
0.419 1.750
0.419 -1.750
0.382 3.142
Polos
1.000 0.000
1.000 0.000
Como se puede observar, hay 2 polos y 4 ceros, lo que resulta en 2 ceros en la respuesta de módulo de \(T_2(z)\).
Análisis simbólico
Dado que a partir de este punto la complejidad algebraica escala, surge la necesidad de continuar mediante un análisis simbólico, comenzamos reproduciendo los mismos resultados que se presentaron más arriba
import sympy as sp
from IPython.display import display
from pytc2.general import print_latex, print_subtitle, a_equal_b_latex_s
z = sp.symbols('z', complex=True)
N, C = sp.symbols('N, C', real=True, positive=True)
# moving average
Tma = 1/N * (1-z**(-N))/(1-z**(-1))
# delay line of (N-1)/2
Tdl = z**(-(N-1)/2)
num, den = (Tdl - Tma).as_numer_denom()
num = (sp.expand(num/(N*z**(N+1)))).powsimp()
den = (sp.expand(den/(N*z**(N+1)))).powsimp()
Tdc_removal = num/den
print_latex(a_equal_b_latex_s('T_{H}(z)', Tdc_removal ))
# Según Rick Lyons, este sistema sería muy bueno para implementarse
# con N múltiplo de 2**NN, dado que el escalado por N sería simplemente
# una rotación a nivel de bits de NN veces a la derecha, y su implementación
# no necesitaría de multiplicaciones. Sin embargo esta elección impone un
# retardo no entero. Por esta razón se opta por poner dos (incluso cuatro)
# sistemas idénticos en cascada.
# Probamos primero con 2 moving average
Tdc_removal_2 = z**-(N-1) - Tma**2
# emprolijamos la expresion a mano
num, den = Tdc_removal_2.as_numer_denom()
num = (sp.expand(num/(N**2*z**(2*N+2))).powsimp())
den = (sp.expand(den/(N**2*z**(2*N+2))).powsimp())
Tdc_removal_2 = num/den
#display(Tdc_removal_2)
print_latex(a_equal_b_latex_s('T_2(z)', Tdc_removal_2 ))
Si bien no es la misma expresión, se puede verificar que estas expresiones son equivalentes a las presentadas más arriba. Esto se debe a los algoritmos de simplificación que usa SymPy. Ahora sí se probará más promediadores en cascada (C = 4 y 6) sólo para tener una idea de cómo escala la expresión de \(T_C(z)\). Recordar que solamente tiene sentido la elección de \(C\) pares, debido al retardo no entero para \(N\) pares.
# Con 4 promediadores
Tdc_removal_C = z**-(N-1)*C - Tma**C
cc = sp.Rational(4)
# emprolijamos la expresion
num, den = (Tdc_removal_C.subs(C, cc)).as_numer_denom()
num = (sp.expand(num/(N**cc*z**(cc*N+cc)))).powsimp()
den = (sp.expand(den/(N**cc*z**(cc*N+cc)))).powsimp()
Tdc_removal_cc = num/den
print_latex(a_equal_b_latex_s(f'T_{cc}(z)', Tdc_removal_cc ))
# Con 6 promediadores
cc = sp.Rational(6)
Tdc_removal_cc = Tdc_removal_C.subs(C, cc)
# emprolijamos la expresion
num, den = Tdc_removal_cc.as_numer_denom()
num = (sp.expand(num/(N**cc*z**(cc*N+cc)))).powsimp()
den = (sp.expand(den/(N**cc*z**(cc*N+cc)))).powsimp()
Tdc_removal_cc = num/den
print_latex(a_equal_b_latex_s(f'T_{cc}(z)', Tdc_removal_cc ))
Análisis numérico de la respuesta en frecuencia
Ahora que se tiene una buena idea de la expresión matemática de las transferencias pasa-alto de \(T_C(z)\), se procederá a analizar la respuesta en frecuencia numérica. Para ello se usará una implementación que se incluye en PyTC2 mediante una clase de filtros, un tanto más complejos, llamada DC_PWL_removal_recursive_filter. Las transferencias pasa-alto \(T_{C}(z)\) analizadas más arriba, pueden implementarse mediante esta clase, omitiendo los parámetros upsample y fpwl, que explicarán más abajo.
# importamos la clase donde se implementa el filtro
from pytc2.filtros_digitales import DC_PWL_removal_recursive_filter
# cantidad de muestras promediadas
NN = [8, 16, 32]
# cantidad de promediadores
cc = 2
fig = plt.figure(figsize=(12, 5))
for nn in NN:
# algunos parámetros se explicarán más abajo
t_filter = DC_PWL_removal_recursive_filter(samp_avg = nn, upsample = 1, cant_ma = cc, fpwl = 0.2 )
w_theoretical, H_theoretical = t_filter.frequency_response(bTeorica=True)
# normalization
w_theoretical = w_theoretical / np.pi
# Gráficos comparativos
plt.plot(w_theoretical, 20*np.log10(np.abs(H_theoretical) + 1e-5), label=f'$T_{C}(z)$: N={nn} C={cc}')
plt.title(f'Respuesta en Frecuencia de Módulo')
plt.ylabel('Magnitud (dB)')
ax = plt.gca()
plt.legend()
plt.grid(True)
En la figura anterior se representa la excelente respuesta de módulo de \(T_2(z)\) en toda la banda. Se nota claramente la presencia de un cero de transmisión en DC y un bajo ripple u ondulación de la respuesta, especialmente en baja frecuencia. Ahora se continuará analizando en detalle ambos aspectos.
ax.set_ylim( [-.6, .1 ])
ax.set_xlim( [0, .21 ])
fig
En este caso se analiza la rapidez o pendiente con la que se transiciona a la banda de paso. Tomando como referencia el nivel de \(T_2(\Omega_C)= -0.5 \) dB, se nota claramente como aumenta la pendiente con el aumento de \(N\). En general, puede verse que el ripple máximo parece estar holgadamente por encima de -0.5 dB para todos los \(N\) analizados. A continuación se analiza la banda de stop o detenida.
ax.set_ylim( [-60, -30 ])
ax.set_xlim( [0.001, .006 ])
fig
Se toma como referencia para el final de la banda de stop la frecuencia \(\Omega_S\), que se define como \(T_2(\Omega_S)= -60 \) dB. En la siguiente tabla se resumen algunas características de la respuesta de magnitud. Se observa cómo aumenta la pendiente de transición y en consecuencia el ancho de banda (BW) del filtro de manera pronunciada.
N |
\(\Omega_S\) |
\(\Omega_C\) |
BW |
|---|---|---|---|
8 |
0.0044 |
0.2 |
0.8 |
16 |
0.0023 |
0.1 |
0.9 |
32 |
0.0011 |
0.05 |
0.95 |
El rendimiento de la respuesta de módulo analizado, junto con su fase lineal y el bajo costo computacional por muestra, dan cuenta de las excelentes características de este filtro. Recordar que la implementación de las respuestas anteriores toman \(2.C\) sumas y 1 producto, para cualquier cantidad de muestras \(N\). Sin embargo, la cantidad de posiciones de memoria necesarias sí escala con \(N\), y en consecuencia también crece el retardo del filtro.
Siguiendo la metodología propuesta por Herrmann et.al. [1], calculamos la cantidad de coeficientes, que es proporcional al costo computacional, para las transiciones de los filtros presentados más arriba, con un ripple de 0.5 dB y una atenuación de -60 dB.
from pytc2.filtros_digitales import herrmann_lp_fir_order
omega_S = [0.0044, 0.0023, 0.0011]
omega_C = [0.2, 0.1, 0.05]
for ii, (wc, ws) in enumerate(zip(omega_C, omega_S)):
cant_coef = herrmann_lp_fir_order([ws*np.pi, wc*np.pi], [0.5, 60], ripple_in_db=True)
cant_coef = np.ceil(cant_coef[0])
print(f'N:{NN[ii]} Costo:{cant_coef}')
N:8 Costo:22
N:16 Costo:45
N:32 Costo:90
En comparativa, un filtro pasa-altos con una transición normalizada (0.001-0.05) correspondiente con \(T_2(z)\), de \(N = 32\), pero diseñado con el método de Parks-McClellan (PM), requiere al menos 90 coeficientes, es decir 90 multiplicaciones, sumas y posiciones de memoria RAM.
from pytc2.filtros_digitales import fir_design_pm
### Rediseñamos para N = 32
# algunos parámetros se explicarán más abajo
t_filter = DC_PWL_removal_recursive_filter(samp_avg = 32, upsample = 1, cant_ma = cc, fpwl = 0.2 )
w_theoretical, H_theoretical = t_filter.frequency_response(bTeorica=True)
# normalization
w_theoretical = w_theoretical / np.pi
# Parámetros comunes del filtro Pasa Bajos
order = 88
fs = 2.0
wpass = 0.05
wstop = 0.0011
ripple = 0.5 # dB
attenuation = 60 # dB
lgrid = 16
cant_coef, band_edges, desired, weights = herrmann_lp_fir_order([wstop*np.pi, wpass*np.pi], [ripple, attenuation], ripple_in_db=True)
b_pm, _, _ = fir_design_pm(cant_coef, band_edges, desired[::-1], fs = 2*np.pi, weight= weights, grid_density=lgrid, filter_type='highpass')
# Comparar la respuesta en frecuencia de ambos filtros
w, h_pm = sig.freqz(b_pm, worN=w_theoretical, fs=2)
# Graficar ambas respuestas
plt.figure(figsize=(12, 5))
# Gráficos comparativos
plt.plot(w_theoretical, 20 * np.log10(np.abs(H_theoretical) + 1e-4), label=f'$T_2(z)$: N=32 C=2')
plt.plot(w_theoretical, 20 * np.log10(np.abs(h_pm)), label='PM ')
plot_plantilla(filter_type = 'highpass' , fpass = wpass, ripple = ripple , fstop = wstop, attenuation = attenuation, fs = fs)
plt.title('Comparación con filtro diseñado mediante Parks-McClellan (PM)')
plt.xlabel('Frecuencia Normalizada (×π rad/sample)')
plt.ylabel('Amplitud [dB]')
plt.legend()
plt.grid()
ax = plt.gca()
ax.set_xlim( [0, .1 ])
plt.show()
Convergencia exitosa!
Como se observa, la cantidad de coeficientes estimada fue subestimada, y es muy posible que el algoritmo de PM tenga problemas para ubicar ceros en torno a DC. Con esta comparativa el filtro propuesto cobra aún más interés. Sin embargo, en el último artículo de Rick Lyons, a partir de la transferencia \(T_C(z)\), se propone interpolar su respuesta de forma tal de convertirla en una transferencia tipo peine. De esta manera podría aprovecharse sus ceros de transmisión para anular interferencias típicas del electrocardiograma (ECG) como
el movimiento de muy baja frecuencia de la línea de base, provocado por factores tales como la respiración, actividad física, sudoración, degradación de los electrodos, entre otros.
la frecuencia de 50/60 Hz de la línea eléctrica, y todas sus armónicas.
El efecto de interpolar por \(U\) la respuesta al impulso de un filtro, visto desde la perspectiva frecuencial, es equivalente a comprimir la respuesta en frecuencia \(U\) veces. Esto implica para \(T_2(z)\)
Como la clase DC_PWL_removal_recursive_filter ya implementa el factor de sobremuestreo \(U\), su efecto se verá más claro al analizar la respuesta en frecuencia de diferentes filtros
# cantidad de promediadores
cc = 2
uu = 2
# algunos parámetros se explicarán más abajo
t_filter = DC_PWL_removal_recursive_filter(samp_avg = nn, upsample = uu, cant_ma = cc )
w_theoretical, H_theoretical = t_filter.frequency_response(bTeorica=True)
# normalization
w_theoretical = w_theoretical / np.pi
fig = plt.figure(figsize=(12, 5))
plt.plot(w_theoretical, 20*np.log10(np.abs(H_theoretical) + 1e-4), label=f'$T_{C}(z)$: U={uu} N={nn} C={cc}')
plt.title(f'Respuesta en Frecuencia de Módulo')
plt.ylabel('Magnitud (dB)')
ax = plt.gca()
ax.set_ylim( [-80, 5 ])
plt.legend()
plt.grid(True)
El factor de sobremuestreo \(U=2\) podría entenderse como la redefinición de la frecuencia de muestreo original, escalada \(U\) veces $\( f_{Su} = U . f_S \)$
De esta manera puede observarse que el 0 en DC de la transferencia \(T_C(z)\) tiene infinitas réplicas en \(k.f_S\), y como
para \(U=2\) simplemente coincide la primer réplica (k=1) con la nueva \(\tfrac{f_{Su}}{2}\), es decir con la nueva frecuencia de Nyquist. Esto se hará aún más evidente para mayores factores de sobremuestreo.
# cantidad de promediadores
cc = 2
uu = 4
# algunos parámetros se explicarán más abajo
t_filter = DC_PWL_removal_recursive_filter(samp_avg = nn, upsample = uu, cant_ma = cc )
w_theoretical, H_theoretical = t_filter.frequency_response(bTeorica=True)
# normalization
w_theoretical = w_theoretical / np.pi
fig = plt.figure(figsize=(12, 5))
plt.plot(w_theoretical, 20*np.log10(np.abs(H_theoretical) + 1e-4), label=f'$T_{C}(z)$: U={uu} N={nn} C={cc}')
plt.title(f'Respuesta en Frecuencia de Módulo')
plt.ylabel('Magnitud (dB)')
ax = plt.gca()
ax.set_ylim( [-80, 5 ])
plt.legend()
plt.grid(True)
# cantidad de promediadores
cc = 2
uu = 8
# algunos parámetros se explicarán más abajo
t_filter = DC_PWL_removal_recursive_filter(samp_avg = nn, upsample = uu, cant_ma = cc )
w_theoretical, H_theoretical = t_filter.frequency_response(bTeorica=True)
# normalization
w_theoretical = w_theoretical / np.pi
fig = plt.figure(figsize=(12, 5))
plt.plot(w_theoretical, 20*np.log10(np.abs(H_theoretical) + 1e-4), label=f'$T_{C}(z)$: U={uu} N={nn} C={cc}')
plt.title(f'Respuesta en Frecuencia de Módulo')
plt.ylabel('Magnitud (dB)')
ax = plt.gca()
ax.set_ylim( [-80, 5 ])
plt.legend()
plt.grid(True)
Se observa que es conveniente que el factor \(U\) sea \(2^n\) múltiplo, de forma tal que siempre haya un cero en Nyquist, y que los ceros de la transferencia peine ocurran en
siendo
De esta manera se obtiene un filtro para la remoción de interferencias típicas del registro de bio-potenciales, de forma compacta
luego para conocer su estructura, podemos recurrir al análisis simbólico nuevamente
# Factor de sobremuestreo
U = sp.symbols('U', real=True, positive=True)
Tdc_removal_C = z**-(N-1)*C - Tma**C
# Con 2 promediadores
Tdc_removal_cc = Tdc_removal_C.subs({C: 2, z:z**U})
# ordenar en potencias de z decrecientes
k, num, den = simplify_n_monic_z(Tdc_removal_cc, poly_var = z, bSortAscending = False)
Tdc_removal_cc = num/den
print_latex(a_equal_b_latex_s(f'T_{{2,U}}(z)', sp.latex(k) + '\cdot' + sp.latex(Tdc_removal_cc) ) )
# Con 4 promediadores
Tdc_removal_cc = Tdc_removal_C.subs({C: 4, z:z**U})
# ordenar en potencias de z decrecientes
k, num, den = simplify_n_monic_z(Tdc_removal_cc, poly_var = z, bSortAscending = False)
Tdc_removal_cc = num/den
print_latex(a_equal_b_latex_s(f'T_{{4,U}}(z)', sp.latex(k) + '\cdot' + sp.latex(Tdc_removal_cc) ) )
Como se aprecia en las expresiones, el hecho de aumentar \(U\) no aumenta la cantidad de coeficientes diferentes de cero, es decir la cantidad de operaciones matemáticas necesarias por muestra de salida. En el mismo sentido, aumentar el factor de interpolación solo aumenta lo dispersos que estarán dichos coeficientes, aunque también la demora total del filtro. Por el contrario, y como es lógico, aumentar la cantidad de promediadores en cascada, \(C\), sí aumenta el costo computacional, como se puede apreciar en ambas expresiones anteriores.
Referencias
[1] Herrmann O., Rabiner L.R., and Chan D S K.: Practical design rules for optimum finite impulse response lowpass digital filters, Bell System Techical Journal, vol. 52 (July-August), 1973.
[2] Lars Wanhammar, Tapio Saramäki - Digital Filters Using MATLAB-Springer (2020)
[3] Rick Lyons. Linear-phase DC Removal Filter. March 30, 2008.
[4] Rick Lyons. 60-Hz Noise and Baseline Drift Reduction in ECG Signal Processing. January 23, 2021.
Parte II: Aplicaciones a la electrocardiografía
Como se explicó anteriormente, se trata de un filtro de muy bajo costo computacional, demora constante y una respuesta de módulo muy prometedora. A partir de este momento, lo pondremos a prueba con señales de ECG registradas durante una prueba de ergometría.
La señal electrocardiográfica (ECG)
El electrocardiograma de superficie (ECG) es una medición no invasiva del potencial eléctrico producido por el dipolo eléctrico resultante en cada célula cardíaca. Por lo general, para la medición del ECG se colocan electrodos en las extremidades del cuerpo y sobre la superficie del torso, como se aprecia en la figura. Debido a que es una técnica no invasiva y su costo no es elevado, es ampliamente utilizada en la clínica médica para conocer la funcionalidad del corazón. Habitualmente se miden entre una y doce señales de ECG dependiendo de la patología bajo análisis.
La señal de ECG también es afectada por los campos eléctricos generados por todas las otras fuentes biológicas del organismo bajo estudio (ruido muscular), como también por otros factores como el movimiento de los electrodos sobre la piel y el acople de la línea eléctrica (50/60 Hz) entre otros. Una de las primeras tareas previas al análisis del ECG es el filtrado de interferencias. En el siguiente ejemplo, procederemos al filtrado de un ECG registrado durante una prueba de esfuerzo o ergometría. En dichos registros, las interferencias causadas por la actividad muscular y la respiración son muy exageradas dado al esfuerzo al que es sometido el paciente. Entonces se procede filtrar la señal de ECG con una transferencia \(T_{C,U}(z)\), con \(C = 2\) y \(U = 20\), de manera que se elimine la interferencia de línea eléctrica de 50 Hz al mismo tiempo que otras interferencias.
Con la señal de ECG filtrada y la original, se realizará un análisis cualitativo respecto de cómo el filtro resulta:
eficaz para mitigar las interferencias, e
inocuo en los tramos de ECG limpio.
Estos dos aspectos determinarán el correcto funcionamiento del filtro diseñado. A continuación se analiza la eficacia en los diferentes tramos ruidosos de la señal.
import scipy.io as sio
mat_struct = sio.loadmat('ecg.mat')
ecg_one_lead = mat_struct['ecg_lead']
ecg_one_lead = ecg_one_lead.flatten().astype(np.float64)
cant_muestras = len(ecg_one_lead)
fs = 1000 # Hz (Normalizamos a fs/2 = f_nyq)
# frec de la interferencia
f_pwl_interference = 50 # Hz
t_filter = DC_PWL_removal_recursive_filter(fs = fs, fpwl = f_pwl_interference, samp_avg = 64, cant_ma = 2)
ECG_f_rl_fin = t_filter.process(ecg_one_lead)
demora_rl = t_filter.demora
# Nota: se obtuvo esta performance en una notebook estandard con:
# HP 255 15.6 inch G10 Notebook PC (A82ZVUA#ABA)
# 8BA5
# AMD Ryzen 7 7730U with Radeon Graphics 64KiB BIOS 512KiB L1 caché 4MiB L2 caché 16MiB L3 caché
# 16GiB Memoria de sistema SODIMM DDR4 Síncrono Unbuffered (Unregistered) 3200 MHz (0,3 ns)
# 1129116 muestras de ECG a fs = 1kHz
# %timeit ECG_f_rl_fin = t_filter.process(ecg_one_lead)
# 1.18 s ± 10.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
regs_interes = (
np.array([10, 10.2]) *60*fs, # minutos a muestras
np.array([12, 12.4]) *60*fs, # minutos a muestras
np.array([15, 15.2]) *60*fs, # minutos a muestras
np.array([18, 18.2]) *60*fs, # minutos a muestras
)
for ii in regs_interes:
# intervalo limitado de 0 a cant_muestras
zoom_region = np.arange(np.max([0, ii[0]]), np.min([cant_muestras-1, ii[1]]), dtype='uint')
plt.figure(figsize=(12, 6))
plt.plot(zoom_region, ecg_one_lead[zoom_region], label='ECG original', linewidth=2)
# FIR con corrección de demora
plt.plot(zoom_region, ECG_f_rl_fin[zoom_region+demora_rl], label='Filtrado')
plt.title(f'ECG con ruido desde {ii[0]:3.0f} a {ii[1]:3.0f}')
plt.ylabel('Adimensional')
plt.xlabel('Muestras (#)')
axes_hdl = plt.gca()
axes_hdl.legend()
plt.show()
De los tramos con ruido anteriores se puede ver claramente que el filtro es eficiente en mitigar las interferencias de baja frecuencia que se manifiestan en todos los tramos visualizados. Esto se debe al excelente desempeño del filtro en la banda de atenuación entre 0 y 0.5 Hz.
A continuación se analizarán algunos tramos limpios de ECG, para analizar cualquier distorsión que introduzca el filtro de manera indeseada.
regs_interes = (
[4000, 5500], # muestras
[10e3, 11e3], # muestras
)
for ii in regs_interes:
# intervalo limitado de 0 a cant_muestras
zoom_region = np.arange(np.max([0, ii[0]]), np.min([cant_muestras-1, ii[1]]), dtype='uint')
plt.figure(figsize=(12, 6))
plt.plot(zoom_region, ecg_one_lead[zoom_region], label='ECG raw', linewidth=2)
# FIR con corrección de demora
plt.plot(zoom_region, ECG_f_rl_fin[zoom_region+demora_rl], label='Filtrada')
plt.title(f'ECG limpio desde {ii[0]:3.0f} a {ii[1]:3.0f}')
plt.ylabel('Adimensional')
plt.xlabel('Muestras (#)')
axes_hdl = plt.gca()
axes_hdl.legend()
plt.show()
En los últimos dos gráficos se observa que el filtro es inocuo, es decir, que en la banda de paso la señal de ECG no es alterada de manera evidente, sobretodo cuando la señal original no presenta interferencias. Esta característica es difícil de conseguir con filtros lineales, sobretodo cuando las bandas de transición son estrechas como en este caso.
Discusión y Conclusiones
En este documento se presentó un filtro digital con algunas ventajas teóricas:
respuesta al impulso finita, esto garantiza la estabilidad
respuesta de fase lineal, o retardo constante, es decir sin distorsión de retardo
implementación recursiva, o muy bajo costo computacional.
Este filtro también demostró ser muy práctico con señales de ECG real, demostrando eficiencia e inocuidad tanto en la banda pasante como en la de rechazo.
Un aspecto negativo que podría ser señalado para este filtro es la imposibilidad de mitigar interferencias en la banda de 35-50 Hz. Si bien la interferencia en esta banda es bastante común en la práctica clínica, también es cierto que podría subsanarse con un filtro pasabajo de bajo costo computacional y fase lineal, de forma tal que no arruine todas las ventajas que se enumeraron.
En resumen, se trata de un filtro cuyas ventajas superan por mucho las desventajas y que debería ser considerado como una alternativa confiable para el preprocesamiento de señales con interferencias de muy baja frecuencia y de senoidales, como la de la línea eléctrica de 50/60 Hz.