added notebook for Laplace transform

GDET3 Laplace-Transformation.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib notebook\n",
+    "\n",
+    "from ipywidgets import interact, interactive, fixed\n",
+    "import ipywidgets as widgets\n",
+    "from IPython.display import clear_output, display, HTML\n",
+    "from IPython.display import Markdown as md\n",
+    "\n",
+    "from src.laplace_plot import pzPoint, pzPlot\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Polstellen bei j , -j (1-fach)\n",
+    "pp1 = 1j\n",
+    "ord_pp1 = 1\n",
+    "# Nullstellen bei 0 (1-fach)\n",
+    "pn1 = 0\n",
+    "ord_pn1 = 1\n",
+    "\n",
+    "pzp = pzPlot(np.array([pp1]),\n",
+    "             np.array([pn1]),\n",
+    "             np.array([ord_pp1]), \n",
+    "             np.array([ord_pn1]))\n",
+    "\n",
+    "action_type = interactive(pzp.update_action,action=widgets.Dropdown(options=list(pzp.action_types.keys()),value=\"Hinzufügen\", description='Modus'))\n",
+    "display(action_type);\n",
+    "point_type = interactive(pzp.update_mode,mode=widgets.Dropdown(options=list(pzp.mode_types.keys()), value=\"Polstelle\", description='Typ'))\n",
+    "display(point_type)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

src/funcs.py

@@ -144,3 +144,107 @@ def addDirac(array, dirac):
         val_str = str(val)
         new_array[i] = val_str + '*dirac'
     return new_array
+
+
+def calc_t_from_p(t=np.linspace(-6, 6, num=1024), H0=1, pp=np.array([]), pz=np.array([]), ord_p=np.array([]),
+                  ord_z=np.array([])):
+    Polstellen = np.array([])
+    GradPolstellen = np.array([])
+    for c, point in enumerate(pp):
+        Polstellen = np.append(Polstellen, point)
+        GradPolstellen = np.append(GradPolstellen, ord_p[c])
+        if np.abs(np.imag(point)) > 0:
+            Polstellen = np.append(Polstellen, np.conj(point))
+            GradPolstellen = np.append(GradPolstellen, ord_p[c])
+
+    Nullstellen = np.array([])
+    GradNullstellen = np.array([])
+    for c, point in enumerate(pz):
+        Nullstellen = np.append(Nullstellen, point)
+        GradNullstellen = np.append(GradNullstellen, ord_z[c])
+        if np.abs(np.imag(point)) > 0:
+            Nullstellen = np.append(Nullstellen, np.conj(point))
+            GradNullstellen = np.append(GradNullstellen, ord_z[c])
+
+    GradPolstellen = GradPolstellen.astype(int)
+    GradNullstellen = GradNullstellen.astype(int)
+
+    if np.sum(GradPolstellen) >= np.sum(GradNullstellen):
+
+        # TODO - t_min, t_max global ; IstKausalerPol Matrix anpassen; Dirac plotten sollte funktionieren
+        t_min = -6
+        t_max = 6
+        IstKausalerPol = np.array([True for x in Polstellen])
+
+        b = np.zeros((np.sum(GradPolstellen) + 1, 1))
+        A = np.zeros((len(b), len(b)), dtype=complex)
+
+        col = 0
+        tmp = np.array([])
+        for i in range(0, len(Polstellen)):
+            tmp = np.append(tmp, np.array(Polstellen[i] * np.ones(GradPolstellen[i])))
+
+        p = np.poly(tmp)
+        A[-len(p):, col] = p
+
+        for i in range(0, len(Polstellen)):
+            for j in range(1, GradPolstellen[i] + 1):
+                tmp = Polstellen[i] * np.ones(GradPolstellen[i] - j)
+                for k in range(0, len(Polstellen)):
+                    if i != k:
+                        tmp = np.append(tmp, np.array(Polstellen[k] * np.ones(GradPolstellen[k])))
+                p = np.poly(tmp)
+
+                col = col + 1
+                try:
+                    A[-len(p):, col] = p
+                except TypeError:
+                    A[-int(p):, col] = p
+
+        ####
+
+        tmp = []
+        for i in range(0, len(Nullstellen)):
+            tmp = np.append(tmp, np.array(Nullstellen[i] * np.ones(GradNullstellen[i])))
+
+        p = H0 * np.poly(tmp)
+
+        try:
+            b[-len(p):, 0] = p.transpose()
+        except AttributeError:
+            b[-int(p):, 0] = p
+
+        try:
+            A_inv = np.around(np.linalg.inv(A), 3)
+        except np.linalg.LinAlgError:
+            # Not invertible.
+            pass
+
+        a = np.dot(A_inv, b)
+
+        ####
+
+        h = np.zeros(t.shape)
+
+        row = 0
+
+        for i in range(0, len(Polstellen)):
+            for j in range(1, GradPolstellen[i] + 1):
+                row = row + 1
+                amplitude = a[row]
+                if j > 1:
+                    amplitude = amplitude / np.prod(range(1, j))
+
+                tmp = amplitude * np.power(t, j - 1) * np.exp(Polstellen[i] * t)
+                if IstKausalerPol[i]:
+                    tmp[np.where(t < 0)] = 0
+                else:
+                    tmp[np.where(t >= 0)] = 0
+                    tmp = -tmp
+
+                h = h + tmp
+
+        return h
+    else:
+        print("Zaehlergrad größer als Nennergrad. Keine Partialbruchzerlegung möglich!")
+        return np.array([np.NaN for x in t])

src/laplace_plot.py

 import bisect
-
-import sympy as sp
+import src.funcs as func
 
 from src.ient_plots import *
 
@@ -92,9 +91,9 @@ class pzPlot():
         self.fig.canvas.mpl_connect('button_press_event', onclick)
 
         self.handles['axh'] = plt.subplot(gs[1, 0])
-        # ient_axis(self.handles['axh']);
+        ient_axis(self.handles['axh'])
         self.handles['axH'] = plt.subplot(gs[1, 1])
-        # ient_axis(self.handles['axH']);
+        ient_axis(self.handles['axH'])
 
         plt.subplots_adjust(wspace=.25)
 
@@ -264,29 +263,29 @@ class pzPlot():
         # plt.close(backend.fig)
 
     def update_plot(self):
-        p = sp.symbols('p')
-        f = sp.symbols('f', real=True)
-        t = sp.symbols('t', real=True, nonnegative=True)
-
-        poles = np.array(list((x.p for x in self.pp)))
+        poles = np.array(list((x.p for x in self.pp)), dtype=complex)
         poles_order = np.array(list((x.order for x in self.pp)))
-        zeroes = np.array(list((x.p for x in self.pz)))
+        zeroes = np.array(list((x.p for x in self.pz)), dtype=complex)
         zeroes_order = np.array(list((x.order for x in self.pz)))
-
         H = 1
-        for pz, oz in zip(zeroes, zeroes_order):
-            oz = int(oz)
-            H *= (p - pz) ** oz
-            if pz.imag > 0:
-                H *= (p - pz.conjugate()) ** oz
-
-        for pp, op in zip(poles, poles_order):
-            op = int(op)
-            H *= 1 / (p - pp) ** op
-            if pp.imag > 0:
-                H *= 1 / (p - pp.conjugate()) ** op
-
-        # Inverse transform. In 2018, sympy was not able to handle bilateral Laplace transforms.
-        # Therefore, we calculate the partial fraction and transform each term back separately.
-        # roc_sigmas = self.roc['sigma']#np.sort(np.array(list((x.p for x in self.pp ))).real)
-        # self.H = H
+        length = 1024
+
+        t = np.linspace(-6, 6, num=length)
+        h_t = func.calc_t_from_p(t, H, poles, zeroes, poles_order, zeroes_order)
+        H_f = np.power(t, 2)
+
+        self.handles['axh'].clear()
+        ient_axis(self.handles['axh'])
+        self.handles['axh'].set_xlim(-6, 6)
+        self.handles['axh'].set_ylim(np.amin(np.real(h_t)) - 0.2, np.amax(np.real(h_t)) + 0.2)
+        self.handles['axh'].set_xlabel(r'$\rightarrow \mathrm{t}$')
+        self.handles['axh'].set_ylabel(r'$\uparrow \mathrm{s(t)}$')
+        self.handles['axh'].plot(t, np.real(h_t))
+
+        self.handles['axH'].clear()
+        ient_axis(self.handles['axH'])
+        self.handles['axH'].set_xlim(-6, 6)
+        self.handles['axH'].set_ylim(np.amin(H_f) - 0.2, np.amax(H_f) + 0.2)
+        self.handles['axH'].set_xlabel(r'$\rightarrow \mathrm{w}$')
+        self.handles['axH'].set_ylabel(r'$\uparrow \mathrm{|S(jw)| [dB]}$')
+        self.handles['axH'].plot(t, H_f)