1 files changed, 149 insertions, 0 deletions
diff --git a/buch/papers/ellfilter/python/elliptic2.py b/buch/papers/ellfilter/python/elliptic2.py
new file mode 100644
index 0000000..29c6f47
--- /dev/null
+++ b/buch/papers/ellfilter/python/elliptic2.py
@@ -0,0 +1,149 @@
+# %%
+
+import matplotlib.pyplot as plt
+import scipy.signal
+import numpy as np
+import matplotlib
+from matplotlib.patches import Rectangle
+
+import plot_params
+
+def ellip_filter(N):
+
+    order = N
+    passband_ripple_db = 3
+    stopband_attenuation_db = 20
+    omega_c = 1
+
+    a, b = scipy.signal.ellip(
+        order,
+        passband_ripple_db,
+        stopband_attenuation_db,
+        omega_c,
+        btype='low',
+        analog=True,
+        output='ba',
+        fs=None
+    )
+
+    w, mag_db, phase = scipy.signal.bode((a, b), w=np.linspace(0*omega_c,2*omega_c, 4000))
+
+    mag = 10**(mag_db/20)
+
+    passband_ripple = 10**(-passband_ripple_db/20)
+    epsilon2 = (1/passband_ripple)**2 - 1
+
+    FN2 = ((1/mag**2) - 1)
+
+    return w/omega_c, FN2 / epsilon2, mag, a, b
+
+
+plt.figure(figsize=(4,2.5))
+
+for N in [5]:
+    w, FN2, mag, a, b = ellip_filter(N)
+    plt.semilogy(w, FN2, label=f"$N={N}, k=0.1$", linewidth=1)
+
+plt.gca().add_patch(Rectangle(
+    (0, 0),
+    1, 1,
+    fc ='green',
+    alpha=0.2,
+    lw = 10,
+))
+plt.gca().add_patch(Rectangle(
+    (1, 1),
+    0.01, 1e2-1,
+    fc ='orange',
+    alpha=0.2,
+    lw = 10,
+))
+
+plt.gca().add_patch(Rectangle(
+    (1.01, 100),
+    1, 1e6,
+    fc ='red',
+    alpha=0.2,
+    lw = 10,
+))
+
+zeros = [0,0.87,1]
+poles = [1.01,1.155]
+
+import matplotlib.transforms
+plt.plot( # mark errors as vertical bars
+    zeros,
+    np.zeros_like(zeros),
+    "o",
+    mfc='none',
+    color='black',
+    transform=matplotlib.transforms.blended_transform_factory(
+        plt.gca().transData,
+        plt.gca().transAxes,
+    ),
+)
+plt.plot( # mark errors as vertical bars
+    poles,
+    np.ones_like(poles),
+    "x",
+    mfc='none',
+    color='black',
+    transform=matplotlib.transforms.blended_transform_factory(
+        plt.gca().transData,
+        plt.gca().transAxes,
+    ),
+)
+
+plt.xlim([0,2])
+plt.ylim([1e-4,1e6])
+plt.grid()
+plt.xlabel("$w$")
+plt.ylabel("$F^2_N(w)$")
+plt.legend()
+plt.tight_layout()
+plt.savefig("F_N_elliptic.pgf")
+plt.show()
+
+
+
+plt.figure(figsize=(4,2.5))
+plt.plot(w, mag, linewidth=1)
+
+plt.gca().add_patch(Rectangle(
+    (0, np.sqrt(2)/2),
+    1, 1,
+    fc ='green',
+    alpha=0.2,
+    lw = 10,
+))
+plt.gca().add_patch(Rectangle(
+    (1, 0.1),
+    0.01, np.sqrt(2)/2 - 0.1,
+    fc ='orange',
+    alpha=0.2,
+    lw = 10,
+))
+
+plt.gca().add_patch(Rectangle(
+    (1.01, 0),
+    1, 0.1,
+    fc ='red',
+    alpha=0.2,
+    lw = 10,
+))
+
+plt.grid()
+plt.xlim([0,2])
+plt.ylim([0,1])
+plt.xlabel("$w$")
+plt.ylabel("$|H(w)|$")
+plt.tight_layout()
+plt.savefig("elliptic.pgf")
+plt.show()
+
+print("zeros", a)
+print("poles", b)
+
+
+
+