From 71d369b3d24a42b5abb66e9b26477472e0100b2f Mon Sep 17 00:00:00 2001
From: Nicolas Tobler <nicolas.tobler@ost.ch>
Date: Sun, 21 Aug 2022 18:06:18 +0200
Subject: added pole zero calculation

---
 buch/papers/ellfilter/elliptic.tex        |   3 +-
 buch/papers/ellfilter/python/elliptic2.py |   7 ++-
 buch/papers/ellfilter/python/elliptic3.py | 101 ++++++++++++++++++++++++++++++
 3 files changed, 108 insertions(+), 3 deletions(-)
 create mode 100644 buch/papers/ellfilter/python/elliptic3.py

(limited to 'buch/papers/ellfilter')

diff --git a/buch/papers/ellfilter/elliptic.tex b/buch/papers/ellfilter/elliptic.tex
index fc9d5b6..26d90f1 100644
--- a/buch/papers/ellfilter/elliptic.tex
+++ b/buch/papers/ellfilter/elliptic.tex
@@ -85,9 +85,10 @@ k_1 = k^N \prod_{i=1}^L \sn^4 \Bigg( \frac{2i - 1}{N} K, k \Bigg),
 N = 2L+r.
 \end{equation}
 Die Herleitung ist sehr umfassend und wird in \cite{ellfilter:bib:orfanidis} im Detail angeschaut.
+$k_1$ jedoch gar nicht berechnet werden, um die elliptisch rationale Funktion zu erhalten.
 
 Um ein elliptisches Filter auszulegen, kann die Ordnung $N$ und der Parameter $k$ gewählt werden.
-$k_1$ muss dann mit \eqref{ellfilter:eq:degeqsol} oder mit numerischen Methoden berechnet werden.
+% $k_1$ muss dann mit \eqref{ellfilter:eq:degeqsol} oder mit numerischen Methoden berechnet werden.
 Je kleiner $k$ gewählt wird, desto grösser wird die Dämpfung des Filters im Sperrbereich im Verhältnis zum Durchlassbereich.
 Allerdings verliert das Filter dabei auch an Steilheit.
 Wenn $k$ und $k_1$ bekannt sind, können die Position der Pol- und Nullstellen $p_i$ und $n_i$ in einem Raster konstruiert werden, wie dargestellt in Abbildung \ref{ellfilter:fig:cd2}.
diff --git a/buch/papers/ellfilter/python/elliptic2.py b/buch/papers/ellfilter/python/elliptic2.py
index 6f03ecf..3d9065d 100644
--- a/buch/papers/ellfilter/python/elliptic2.py
+++ b/buch/papers/ellfilter/python/elliptic2.py
@@ -29,6 +29,9 @@ def ellip_filter(N, mode=-1):
         fs=None
     )
 
+    print("poles", a)
+    print("zeros", b)
+
     if mode == 0:
         w = np.linspace(0*omega_c,omega_c, 2000)
     elif mode == 1:
@@ -148,8 +151,8 @@ plt.tight_layout()
 plt.savefig("elliptic.pgf")
 plt.show()
 
-print("zeros", a)
-print("poles", b)
+print("poles", a)
+print("zeros", b)
 
 
 
diff --git a/buch/papers/ellfilter/python/elliptic3.py b/buch/papers/ellfilter/python/elliptic3.py
new file mode 100644
index 0000000..10accbb
--- /dev/null
+++ b/buch/papers/ellfilter/python/elliptic3.py
@@ -0,0 +1,101 @@
+# %%
+
+import matplotlib.pyplot as plt
+import scipy.signal
+import numpy as np
+import matplotlib
+from matplotlib.patches import Rectangle
+import scipy.special
+import scipyx as spx
+
+# import plot_params
+
+def last_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
+
+N = 6
+L = (N//2) * 2
+r = N - L
+
+k = 0.9143
+
+i = np.arange(1, L+1)
+ui = (2*i - 1) / N
+k1 = k**N * np.prod(sn(ui*ell_int(k), k)**4)
+k1 = 0.0165
+k1 = 0.0058
+
+
+kp = np.sqrt(1-k**2)
+k1p = np.sqrt(1-k1**2)
+
+K = ell_int(k)
+Kp = ell_int(kp)
+K1 = ell_int(k1)
+K1p = ell_int(k1p)
+
+# assert np.allclose(Kp*K1*N/K, K1p, rtol=0.001)
+
+zeros = K/N * (np.arange(N)*2 + 1)
+poles = zeros + (1j * Kp)
+# if len(poles) % 2 == 0:
+#     poles = np.delete(poles, len(poles)//2)
+
+
+plt.plot(np.real(zeros), np.imag(zeros), "o")
+plt.plot(np.real(poles), np.imag(poles), "x")
+# plt.plot([0,K1], [0,K1p])
+# plt.plot([0,K], [0,Kp])
+plt.show()
+
+zeros = cd(zeros, k)
+poles = cd(poles, k)
+
+plt.plot(np.real(zeros), np.imag(zeros), "o")
+plt.plot(np.real(poles), np.imag(poles), "x")
+plt.ylim([-0.1,0.1])
+plt.xlim([-2.5,2.5])
+plt.show()
+
+w = np.linspace(0,2, 2000)
+
+def make_RN(w):
+    y = np.prod(w[:, None] - zeros[None], axis=-1) / np.prod(w[:, None] - poles[None], axis=-1)
+    y /= np.prod(1 - zeros) / np.prod(1 - poles)
+    return y
+
+
+RN = make_RN(w)
+
+plt.semilogy(w, np.abs(RN))
+plt.ylim([0.1,1000])
+
+plt.plot(w, np.ones_like(w) / k1)
+
+plt.show()
+
+H = 1 / (1 + RN**2)
+
+plt.semilogy(w, np.abs(H))
+plt.ylim([0.00001,1])
+plt.show()
-- 
cgit v1.2.1