diff options
Diffstat (limited to '')
-rw-r--r-- | buch/papers/ellfilter/python/elliptic2.py | 152 |
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) + + + + |