1 files changed, 45 insertions, 42 deletions
diff --git a/buch/papers/ellfilter/python/elliptic2.py b/buch/papers/ellfilter/python/elliptic2.py
index 29c6f47..20a7428 100644
--- a/buch/papers/ellfilter/python/elliptic2.py
+++ b/buch/papers/ellfilter/python/elliptic2.py
@@ -1,5 +1,6 @@
 # %%
 
+import enum
 import matplotlib.pyplot as plt
 import scipy.signal
 import numpy as np
@@ -8,7 +9,9 @@ from matplotlib.patches import Rectangle
 
 import plot_params
 
-def ellip_filter(N):
+N=5
+
+def ellip_filter(N, mode=-1):
 
     order = N
     passband_ripple_db = 3
@@ -26,7 +29,16 @@ def ellip_filter(N):
         fs=None
     )
 
-    w, mag_db, phase = scipy.signal.bode((a, b), w=np.linspace(0*omega_c,2*omega_c, 4000))
+    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)
 
@@ -38,105 +50,96 @@ def ellip_filter(N):
     return w/omega_c, FN2 / epsilon2, mag, a, b
 
 
-plt.figure(figsize=(4,2.5))
+f, axs = plt.subplots(2, 1, figsize=(5,3), sharex=True)
 
-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)
+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)
 
-plt.gca().add_patch(Rectangle(
+axs[0].add_patch(Rectangle(
     (0, 0),
     1, 1,
     fc ='green',
     alpha=0.2,
     lw = 10,
 ))
-plt.gca().add_patch(Rectangle(
+axs[0].add_patch(Rectangle(
     (1, 1),
-    0.01, 1e2-1,
+    0.00992, 1e2-1,
     fc ='orange',
     alpha=0.2,
     lw = 10,
 ))
 
-plt.gca().add_patch(Rectangle(
-    (1.01, 100),
+axs[0].add_patch(Rectangle(
+    (1.00992, 100),
     1, 1e6,
     fc ='red',
     alpha=0.2,
     lw = 10,
 ))
 
-zeros = [0,0.87,1]
+zeros = [0,0.87,0.995]
 poles = [1.01,1.155]
 
 import matplotlib.transforms
-plt.plot( # mark errors as vertical bars
+axs[0].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,
+        axs[0].transData,
+        axs[0].transAxes,
     ),
 )
-plt.plot( # mark errors as vertical bars
+axs[0].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,
+        axs[0].transData,
+        axs[0].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)
+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)
 
-plt.gca().add_patch(Rectangle(
+axs[1].add_patch(Rectangle(
     (0, np.sqrt(2)/2),
     1, 1,
     fc ='green',
     alpha=0.2,
     lw = 10,
 ))
-plt.gca().add_patch(Rectangle(
+axs[1].add_patch(Rectangle(
     (1, 0.1),
-    0.01, np.sqrt(2)/2 - 0.1,
+    0.00992, np.sqrt(2)/2 - 0.1,
     fc ='orange',
     alpha=0.2,
     lw = 10,
 ))
 
-plt.gca().add_patch(Rectangle(
-    (1.01, 0),
+axs[1].add_patch(Rectangle(
+    (1.00992, 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)|$")
+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()