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