diff options
Diffstat (limited to 'buch/papers/ellfilter/python/elliptic.py')
| -rw-r--r-- | buch/papers/ellfilter/python/elliptic.py | 162 |
1 files changed, 101 insertions, 61 deletions
diff --git a/buch/papers/ellfilter/python/elliptic.py b/buch/papers/ellfilter/python/elliptic.py index 9f209e9..b3336a1 100644 --- a/buch/papers/ellfilter/python/elliptic.py +++ b/buch/papers/ellfilter/python/elliptic.py @@ -5,19 +5,62 @@ import scipy.special import scipyx as spx import numpy as np import matplotlib.pyplot as plt -import matplotlib from matplotlib.patches import Rectangle -matplotlib.rcParams.update({ - "pgf.texsystem": "pdflatex", - 'font.family': 'serif', - 'font.size': 9, - 'text.usetex': True, - 'pgf.rcfonts': False, -}) +import plot_params def last_color(): - plt.gca().lines[-1].get_color() + return plt.gca().lines[-1].get_color() + +# define elliptic functions + +def ell_int(k): + """ Calculate K(k) """ + m = k**2 + return scipy.special.ellipk(m) + +def sn(z, k): + return spx.ellipj(z, k**2)[0] + +def cn(z, k): + return spx.ellipj(z, k**2)[1] + +def dn(z, k): + return spx.ellipj(z, k**2)[2] + +def cd(z, k): + sn, cn, dn, ph = spx.ellipj(z, k**2) + return cn / dn + +# https://mathworld.wolfram.com/JacobiEllipticFunctions.html eq 3-8 + +def sn_inv(z, k): + m = k**2 + return scipy.special.ellipkinc(np.arcsin(z), m) + +def cn_inv(z, k): + m = k**2 + return scipy.special.ellipkinc(np.arccos(z), m) + +def dn_inv(z, k): + m = k**2 + x = np.sqrt((1-z**2) / k**2) + return scipy.special.ellipkinc(np.arcsin(x), m) + +def cd_inv(z, k): + m = k**2 + x = np.sqrt(((m - 1) * z**2) / (m*z**2 - 1)) + return scipy.special.ellipkinc(np.arccos(x), m) + + +k = 0.8 +z = 0.5 + +assert np.allclose(sn_inv(sn(z ,k), k), z) +assert np.allclose(cn_inv(cn(z ,k), k), z) +assert np.allclose(dn_inv(dn(z ,k), k), z) +assert np.allclose(cd_inv(cd(z ,k), k), z) + # %% Buttwerworth filter F_N plot @@ -37,7 +80,7 @@ plt.gca().add_patch(Rectangle( plt.gca().add_patch(Rectangle( (1, 1), 0.5, 1, - fc ='green', + fc ='orange', alpha=0.2, lw = 10, )) @@ -47,7 +90,8 @@ plt.grid() plt.xlabel("$w$") plt.ylabel("$F^2_N(w)$") plt.legend() -plt.savefig("F_N_butterworth.pdf") +plt.tight_layout() +plt.savefig("F_N_butterworth.pgf") plt.show() # %% Cheychev filter F_N plot @@ -69,7 +113,7 @@ plt.gca().add_patch(Rectangle( plt.gca().add_patch(Rectangle( (1, 1), 0.5, 1, - fc ='green', + fc ='orange', alpha=0.2, lw = 10, )) @@ -79,57 +123,10 @@ plt.grid() plt.xlabel("$w$") plt.ylabel("$F^2_N(w)$") plt.legend() -plt.savefig("F_N_chebychev.pdf") +plt.tight_layout() +plt.savefig("F_N_chebychev.pgf") plt.show() -# %% define elliptic functions - -def ell_int(k): - """ Calculate K(k) """ - m = k**2 - return scipy.special.ellipk(m) - -def sn(z, k): - return spx.ellipj(z, k**2)[0] - -def cn(z, k): - return spx.ellipj(z, k**2)[1] - -def dn(z, k): - return spx.ellipj(z, k**2)[2] - -def cd(z, k): - sn, cn, dn, ph = spx.ellipj(z, k**2) - return cn / dn - -# https://mathworld.wolfram.com/JacobiEllipticFunctions.html eq 3-8 - -def sn_inv(z, k): - m = k**2 - return scipy.special.ellipkinc(np.arcsin(z), m) - -def cn_inv(z, k): - m = k**2 - return scipy.special.ellipkinc(np.arccos(z), m) - -def dn_inv(z, k): - m = k**2 - x = np.sqrt((1-z**2) / k**2) - return scipy.special.ellipkinc(np.arcsin(x), m) - -def cd_inv(z, k): - m = k**2 - x = np.sqrt(((m - 1) * z**2) / (m*z**2 - 1)) - return scipy.special.ellipkinc(np.arccos(x), m) - - -k = 0.8 -z = 0.5 - -assert np.allclose(sn_inv(sn(z ,k), k), z) -assert np.allclose(cn_inv(cn(z ,k), k), z) -assert np.allclose(dn_inv(dn(z ,k), k), z) -assert np.allclose(cd_inv(cd(z ,k), k), z) # %% plot arcsin @@ -314,3 +311,46 @@ for n in (1,2,3,4): plt.plot(omega, np.abs(G)) plt.grid() plt.show() + + + + +# %% + + +k = np.concatenate(([0.00001,0.0001,0.001], np.linspace(0,1,101)[1:-1], [0.999,0.9999, 0.99999]), axis=0) +K = ell_int(k) +K_prime = ell_int(np.sqrt(1-k**2)) + + +f, axs = plt.subplots(1,2, figsize=(5,2.5)) +axs[0].plot(k, K, linewidth=0.1) +axs[0].text(k[30], K[30]+0.1, f"$K$") +axs[0].plot(k, K_prime, linewidth=0.1) +axs[0].text(k[30], K_prime[30]+0.1, f"$K^\prime$") +axs[0].set_xlim([0,1]) +axs[0].set_ylim([0,4]) +axs[0].set_xlabel("$k$") + +axs[1].axvline(x=np.pi/2, color="gray", linewidth=0.5) +axs[1].axhline(y=np.pi/2, color="gray", linewidth=0.5) +axs[1].text(0.1, np.pi/2 + 0.1, "$\pi/2$") +axs[1].text(np.pi/2+0.1, 0.1, "$\pi/2$") +axs[1].plot(K, K_prime, linewidth=1) + +k = np.array([0.1,0.2,0.4,0.6,0.9,0.99]) +K = ell_int(k) +K_prime = ell_int(np.sqrt(1-k**2)) + +axs[1].plot(K, K_prime, '.', color=last_color(), markersize=2) +for x, y, n in zip(K, K_prime, k): + axs[1].text(x+0.1, y+0.1, f"$k={n:.2f}$", rotation_mode="anchor") +axs[1].set_ylabel("$K^\prime$") +axs[1].set_xlabel("$K$") +axs[1].set_xlim([0,6]) +axs[1].set_ylim([0,5]) +plt.tight_layout() +plt.savefig("k.pgf") +plt.show() + +print(K[0], K[-1]) |