1 files changed, 152 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..20a7428
--- /dev/null
+++ b/buch/papers/ellfilter/python/elliptic2.py
@@ -0,0 +1,152 @@
+# %%
+
+import enum
+import matplotlib.pyplot as plt
+import scipy.signal
+import numpy as np
+import matplotlib
+from matplotlib.patches import Rectangle
+
+import plot_params
+
+N=5
+
+def ellip_filter(N, mode=-1):
+
+    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
+    )
+
+    if mode == 0:
+        w = np.linspace(0*omega_c,omega_c, 2000)
+    elif mode == 1:
+        w = np.linspace(omega_c,1.00992*omega_c, 2000)
+    elif mode == 2:
+        w = np.linspace(1.00992*omega_c,2*omega_c, 2000)
+    else:
+        w = np.linspace(0*omega_c,2*omega_c, 4000)
+
+    w, mag_db, phase = scipy.signal.bode((a, b), w=w)
+
+    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
+
+
+f, axs = plt.subplots(2, 1, figsize=(5,3), sharex=True)
+
+for mode, c in enumerate(["green", "orange", "red"]):
+    w, FN2, mag, a, b = ellip_filter(N, mode=mode)
+    axs[0].semilogy(w, FN2, label=f"$N={N}, k=0.1$", linewidth=1, color=c)
+
+axs[0].add_patch(Rectangle(
+    (0, 0),
+    1, 1,
+    fc ='green',
+    alpha=0.2,
+    lw = 10,
+))
+axs[0].add_patch(Rectangle(
+    (1, 1),
+    0.00992, 1e2-1,
+    fc ='orange',
+    alpha=0.2,
+    lw = 10,
+))
+
+axs[0].add_patch(Rectangle(
+    (1.00992, 100),
+    1, 1e6,
+    fc ='red',
+    alpha=0.2,
+    lw = 10,
+))
+
+zeros = [0,0.87,0.995]
+poles = [1.01,1.155]
+
+import matplotlib.transforms
+axs[0].plot( # mark errors as vertical bars
+    zeros,
+    np.zeros_like(zeros),
+    "o",
+    mfc='none',
+    color='black',
+    transform=matplotlib.transforms.blended_transform_factory(
+        axs[0].transData,
+        axs[0].transAxes,
+    ),
+)
+axs[0].plot( # mark errors as vertical bars
+    poles,
+    np.ones_like(poles),
+    "x",
+    mfc='none',
+    color='black',
+    transform=matplotlib.transforms.blended_transform_factory(
+        axs[0].transData,
+        axs[0].transAxes,
+    ),
+)
+
+for mode, c in enumerate(["green", "orange", "red"]):
+    w, FN2, mag, a, b = ellip_filter(N, mode=mode)
+    axs[1].plot(w, mag, linewidth=1, color=c)
+
+axs[1].add_patch(Rectangle(
+    (0, np.sqrt(2)/2),
+    1, 1,
+    fc ='green',
+    alpha=0.2,
+    lw = 10,
+))
+axs[1].add_patch(Rectangle(
+    (1, 0.1),
+    0.00992, np.sqrt(2)/2 - 0.1,
+    fc ='orange',
+    alpha=0.2,
+    lw = 10,
+))
+
+axs[1].add_patch(Rectangle(
+    (1.00992, 0),
+    1, 0.1,
+    fc ='red',
+    alpha=0.2,
+    lw = 10,
+))
+
+axs[0].set_xlim([0,2])
+axs[0].set_ylim([1e-4,1e6])
+axs[0].grid()
+axs[0].set_ylabel("$F^2_N(w)$")
+axs[1].grid()
+axs[1].set_ylim([0,1])
+axs[1].set_ylabel("$|H(w)|$")
+plt.tight_layout()
+plt.savefig("elliptic.pgf")
+plt.show()
+
+print("zeros", a)
+print("poles", b)
+
+
+
+