Diseñando filtros de radio con SymPy
Diseñando filtros de radio con SymPy¶
Andaba probando un SDR y me topé con que en los 10m hay una interferencia bastante fea de las radios FM. Investigué y parece que hay que filtrar estas radios que saturan el aparato SDR.
Fue difícil encontrar buena información en línea acerca de los filtros para radiofrecuencia. Este es un tema que vi en la licenciatura, en cursos de electricidad y en cursos sobre ecuaciones diferenciales y sobre números complejos. Pero no lo recordaba bien, así que busque y no pude encontrar la derivación de las ecuaciones. Así que me puse a modelar los filtros con SymPy para ahorrarme trabajo y la complejidad. Una vez entendí cómo usar SymPy para estas cosas, fue bastante facil analizar diferentes circuitos con diferentes parámetros. La mejor información que encontre en línea fue esta guía para mejorar la recepción de señales de baja frecuencia, esta guía de un laboratorio de electrónica y estas diapositivas de un curso de electrónica.
Aún me falta experimentar con estos filtros con un generador de ruido. Creo que será interesante hacer uno con los muchos diagramas que hay por ahí y medir la atenuación que estos filtros hacen en la señal en el mundo real.
import sympy
i = sympy.I
Circuito RC¶
El típico filtro pasa bajos (por alguna razón que aún no he encontrado no funciona bien en RF)
Vi, Vc, It, Zc = sympy.symbols("Vi Vc It Zc")
R, C, w = sympy.symbols("R C w", real=True)
Zc = 1/(i*w*C)
It = Vi/(R+Zc)
Vc = It * Zc
Vr = It * R
iw = i*w
Hc = (Vc/Vi)
Hr = (Vr/Vi)
Hr
import numpy
import matplotlib.pyplot as plt
_C = 100e-9
_R = 1e3
ws = numpy.arange(1,100000, 100)
Hs = []
for wi in ws:
Hs.append(sympy.re(Hc.evalf(subs={R: _R, C: _C, w: wi})))
La frecuencia donde la función de transferencia es 0.5:
eq = sympy.Eq(sympy.re(Hc),0.5)
eq
w0 = sympy.solve(eq, w)
w0
[-1/(C*R), 1/(C*R)]
Ahora buscamos el punto donde la ganancia cae a 0.5, esto es en -3dB
numpy.log10(0.5)*10
-3.010299956639812
_w0 = w0[1].evalf(subs={C:_C, R:_R})/(2*numpy.pi)
_w0
$\omega_0 = 1.59Khz $
plt.plot(numpy.log10(ws), numpy.log10(numpy.array(Hs, float))*10)
t_w0 = sympy.log(_w0*2*numpy.pi)/numpy.log(10)
plt.axvline(t_w0, alpha=0.5)
plt.axhline(-3, alpha=0.5)
plt.text(t_w0, -3, "1.59KHz\n-3dB", bbox=dict(facecolor='red', alpha=0.1) )
Text(4.00000000000000, -3, '1.59KHz\n-3dB')
Filtro pasabajos LC¶
L1, C1 = sympy.symbols("L1 C1", real=True)
Vi, Vo, It, Zt, Z1, Z2 = sympy.symbols("Vi Vo It Zt Z1 Z2")
Vc1, Vl1, Il1, Ic1 = sympy.symbols("Vc1 Vl1 Il1 Ic1")
w = sympy.symbols("w", real=True)
#It = Vi/Zt
Z1 = i*w*L1
Z2 = 1/(i*w*C1)
Il1 = Ic1
Ic1 = Vi/(Z1+Z2)
Vc1 = Ic1*Z2
Vl1 = Vi*Z1/(Z1+Z2)
Zt = 1/ (1/Z1 + 1/Z2)
Vc1
Vc1.simplify()
Vo = Vc1.simplify()
H = Vo/Vi
H
Estuve un buen rato averiguando porqué esta función no produce una gráfica de H que tenga sentido. Resulta que la función 1/(1-x^2) tiende al infinito en $C_1L_1 \omega^2 = 1$ (positivo por la izquierda y negativo por la derecha).
_C = 200000e-12
_L = 100e-9
frange = numpy.sqrt(1/(_C*_L)) * 1.5
ws = numpy.arange(1, frange, frange/1e3)
Hs = []
for wi in ws:
Hs.append(sympy.re(H).evalf(subs={C1: _C, L1: _L, w: wi}))
plt.plot(ws, Hs)
[<matplotlib.lines.Line2D at 0x7ff39d8a6cb0>]
La frecuencia en donde está este polo es: $ \omega = \frac{1}{\sqrt{LC}} $
f0 = 1/numpy.sqrt(_L*_C) / 2*numpy.pi
print(f0/1e3, "KHz")
11107.207345395916 KHz
Buscando rápidamente en internet no encontré nada que explicara cómo se supone que funciona esto en la vida real. Obviamente el voltaje no puede tender a infinito, y el filtro no tiene transferencia totalmente de 0 en las frecuencias arriba de $\omega_0$
Notch Filter (RLC)¶
Ya que tenemos un mayor entendimiento de cómo usar SymPy, vamos a probar estos otros filtros más interesantes:
R1, L1, C1 = sympy.symbols("L1 R1 C1", real=True)
I, Vi, Vo, Z1, Z2, Z3 = sympy.symbols("I Vi Vo Z1 Z2 Z3")
w = sympy.symbols("w", real=True)
Z1 = R1
Z2 = i*w*L1
Z3 = 1/(i*w*C1)
I = Vi/(Z1+Z2+Z3)
Vo = I*(Z2+Z3)
Vo
H = Vo/Vi
H_1 = sympy.re(H).simplify()
H_1
De las guías que mencioné arriba, vi que $\frac{1}{\sqrt{L C}}$ es el centro del filtro.
$ $
subs = {
R1: 47,
C1: 22e-12,
L1: 120e-9
}
fmax = 1e12
ws = numpy.exp(numpy.arange(1, numpy.log(fmax), 0.1))
Hs = []
for wi in ws:
subs[w] = wi
Hs.append(sympy.re(H_1.evalf(subs=subs)))
plt.plot(numpy.log(ws), 10*numpy.log10(numpy.array(Hs, float)))
w0 = float(1/numpy.sqrt(subs[C1]*subs[L1]))
plt.axvline(numpy.log(w0),alpha=0.5)
plt.text(numpy.log(w0), 0,
str(numpy.round(w0/(2*numpy.pi) / 1e6, 2)) + "MHz")
plt.fill_betweenx([-4,0], numpy.log(87e6*2*numpy.pi), numpy.log(108e6*2*numpy.pi),alpha=0.5,color="red")
plt.text(numpy.log(90e6*2*numpy.pi), -5, "Radio FM")
plt.fill_betweenx([-4,0], numpy.log(26e6*2*numpy.pi), numpy.log(29e6*2*numpy.pi),alpha=0.5,color="blue")
plt.text(numpy.log(26e6*2*numpy.pi), -5, "10m y CB")
plt.xlim(17,22)
(17.0, 22.0)
Bandpass Filter (RLC)¶
Este filtro tiene L1 y C1 en paralelo, y es lo opuesto (más o menos) al filtro anterior que los tiene en serie.
R1, L1, C1 = sympy.symbols("L1 R1 C1", real=True)
I, Vi, Vo, Z1, Z2, Z3 = sympy.symbols("I Vi Vo Z1 Z2 Z3")
w = sympy.symbols("w", real=True)
Z1 = R1
Z2 = i*w*L1
Z3 = 1/(i*w*C1)
I = Vi/(Z1+(1/(1/Z2+1/Z3)))
Vo = I*(1/(1/Z2+1/Z3))
Vo
H = Vo/Vi
H_1 = sympy.re(H).simplify()
H_1
subs = {
R1: 5,
C1: 1000e-12,
L1: 33e-9
}
fmax = 1e12
ws = numpy.exp(numpy.arange(1, numpy.log(fmax), 0.1))
Hs = []
for wi in ws:
subs[w] = wi
Hs.append(sympy.re(H.evalf(subs=subs)))
plt.plot(numpy.log(ws), 10*numpy.log10(numpy.array(Hs, float)))
w0 = float(1/numpy.sqrt(subs[C1]*subs[L1]))
plt.axvline(numpy.log(w0),alpha=0.5)
plt.text(numpy.log(w0), 0,
str(numpy.round(w0/(2*numpy.pi) / 1e6, 2)) + "MHz")
plt.fill_betweenx([-4,0], numpy.log(87e6*2*numpy.pi), numpy.log(108e6*2*numpy.pi),alpha=0.5,color="red")
plt.text(numpy.log(90e6*2*numpy.pi), -5, "Radio FM")
plt.fill_betweenx([-4,0], numpy.log(26e6*2*numpy.pi), numpy.log(29e6*2*numpy.pi),alpha=0.5,color="blue")
plt.text(numpy.log(26e6*2*numpy.pi), -5, "10m band")
plt.xlim(17,22)
plt.ylim(-10,0)
(-10.0, 0.0)
Conclusión¶
Fue divertido aprender un poco más de SymPy con estos circuitos prácticos, aunque por mucho que este haga el trabajo por uno, uno no puede dejar de tener un entendimiento fundamental de la matemática, la física/electrónica y también de la programación. SymPy no hace toda la tarea, pero al menos sí más de la mitad, diría yo. Dejo todo el código a la vista para mostrar cómo reproducir los cálculos.
Luego tocará experimentar.