Add new GUI for discrete FT (adapted from Praktikum TI)

5940beec · Christian Rohlfing · 0ccc0710 · 5940beec
Verified Commit 5940beec authored 4 years ago by Christian Rohlfing
--- a/GDET3 Diskrete Fourier-Transformation GUI.ipynb
+++ b/GDET3 Diskrete Fourier-Transformation GUI.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University\n",
+    "%matplotlib widget\n",
+    "\n",
+    "from ipywidgets import interact, interactive, fixed, HBox, VBox\n",
+    "import ipywidgets as widgets\n",
+    "from IPython.display import clear_output, display, HTML\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "from rwth_nb.plots import colors\n",
+    "import rwth_nb.plots.mpl_decorations as rwth_plt\n",
+    "from rwth_nb.misc.signals import *\n",
+    "from rwth_nb.misc.transforms import *\n",
+    "\n",
+    "t,deltat = np.linspace(-10,10,50001,  retstep=True) # t-axis\n",
+    "f,deltaf = np.linspace(-10,10,len(t), retstep=True) # f-axis\n",
+    "kMax = 16 # number of k values in sum for Sa(f)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "<div>\n",
+    "    <img src=\"figures/rwth_ient_logo@2x.png\" style=\"float: right;height: 5em;\">\n",
+    "</div>\n",
+    "\n",
+    "# Demonstrator Diskrete Fourier-Transformation\n",
+    "\n",
+    "Zum Starten: Im Menü: Run <span class=\"fa-chevron-right fa\"></span> Run All Cells auswählen.\n",
+    "\n",
+    "## Einleitung\n",
+    "\n",
+    "$s(t)$ wird ideal mit Abtastrate $r$ abgetastet.\n",
+    "\n",
+    "Diskrete Fourier-Transformation mit $F = \\frac{1}{M}$\n",
+    "$$\n",
+    "S_\\mathrm{d}(k) \n",
+    "= \n",
+    "\\sum_{n=0}^{M-1} s_\\mathrm{d}(n) \\mathrm{e}^{-\\mathrm{j}2\\pi kFn} \\quad k=0,\\ldots,M-1 \\text{ ,}\n",
+    "$$\n",
+    "wobei $s_\\mathrm{d}(n)$ um $M$ periodisches, diskretes Zeitsignal mit $t = n \\cdot \\frac{1}{r}$.\n",
+    "\n",
+    "Für die Rücktransformation in den Zeitbereich gilt dann\n",
+    "$$\n",
+    "s_\\mathrm{d}(n)\n",
+    "=\n",
+    "\\frac{1}{M} \\sum_{k=0}^{M-1}S_\\mathrm{d}(k) \\mathrm{e}^{\\mathrm{j}2\\pi kFn} \\quad n=0,\\ldots,M-1 \\text{ ,}\n",
+    "$$\n",
+    "wobei\n",
+    "$f = k \\cdot \\frac{r}{M}$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "M = 24; # period of discrete signal\n",
+    "r = 4; # sampling rate\n",
+    "\n",
+    "n = np.arange(M) # discrete time axis\n",
+    "sd = tri(n/r) + tri((n-M)/r) # discrete signal sd(n)\n",
+    "\n",
+    "F = 1/M \n",
+    "Sd = np.zeros_like(sd, dtype=np.complex)\n",
+    "for k in np.arange(M): # very slow computation for teaching purpose! Fasten your seatbelts with FFT :)\n",
+    "    Sd[k] = np.sum(sd * np.exp(-2j*np.pi*k*F*n))\n",
+    "#Sd = np.fft.fft(sd)\n",
+    "Sd = np.real(Sd) # discard super small imaginary part (1e-16)\n",
+    "\n",
+    "# Plot\n",
+    "fig,axs = plt.subplots(2,1)\n",
+    "ax = axs[0]; rwth_plt.stem(ax, n, sd); \n",
+    "ax.set_xlabel(r'$\\rightarrow n$'); ax.set_ylabel(r'$\\uparrow s_\\mathrm{d}(n)$')\n",
+    "rwth_plt.axis(ax);\n",
+    "\n",
+    "ax = axs[1]; rwth_plt.stem(ax, np.arange(M), Sd/r); \n",
+    "ax.set_xlabel(r'$\\rightarrow k$'); ax.set_ylabel(r'$\\uparrow S_\\mathrm{d}(k)/r$')\n",
+    "rwth_plt.axis(ax); fig.tight_layout()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Demonstrator"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "source_hidden": true
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def decorate_second_axis(ax):\n",
+    "    ax.spines['top'].set_position(('axes',.9)); ax.spines['right'].set_color('none'); ax.spines['left'].set_color('none'); ax.spines['bottom'].set_color('none');\n",
+    "    ax.xaxis.set_label_coords(.98,1); ax.xaxis.label.set_horizontalalignment('right') \n",
+    "\n",
+    "fig, axs = plt.subplots(2, 1, **rwth_plt.landscape); #fig.tight_layout()\n",
+    "highlight_rect_style = {'facecolor': \"none\", 'edgecolor': 'rwth:orange', 'linestyle':'--', 'linewidth':2.0}\n",
+    "\n",
+    "@widgets.interact(  r = widgets.IntSlider(min=0, max=10, step=2, value=0, description='$r$ [Hz]'),\n",
+    "                  Mbr = widgets.IntSlider(min=0, max=10, step=2, value=0, description='$M/r$'))\n",
+    "def update_plots(r, Mbr):\n",
+    "    M = r*Mbr\n",
+    "    global n, snT, k, SkF\n",
+    "    # Continuous functions s(t) and S(f)\n",
+    "    s = tri\n",
+    "    S = lambda f: si(np.pi*f) ** 2\n",
+    "    \n",
+    "    # Ideal sampling\n",
+    "    Sa_f = np.zeros_like(S(f))\n",
+    "    nT = []; snT = nT; n = nT;\n",
+    "    if r > 0:\n",
+    "        T = 1/r\n",
+    "        \n",
+    "        # Sampled function s(n)\n",
+    "        nT, snT = sample(t, s(t), T); n = nT/T;\n",
+    "        \n",
+    "        # Sampled continuous spectrum S_a(f)\n",
+    "        for k in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements \n",
+    "            Sa_f += S(f-k/T)\n",
+    "        #Sa_f = Sa_f/T # normalize here\n",
+    "        \n",
+    "        # Periodic copies of s(n)\n",
+    "        kF = []; SkF = kF; k = kF;\n",
+    "        sd = np.zeros_like(snT)\n",
+    "        if M > 0:\n",
+    "            F = 1/M # in Hz\n",
+    "            kF, SkF = sample(f, Sa_f, r*F); k = kF/F/r\n",
+    "            for d in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements \n",
+    "                sd += s((n-d*M)*T)\n",
+    "    display('M={}'.format(M))\n",
+    "    \n",
+    "    # Plot\n",
+    "    if not axs[0].lines: # Call plot() and decorate axes. Usually, these functions take some processing time\n",
+    "        ax = axs[0]; ax.set_title('Zeitbereich');\n",
+    "        ax.plot(t, s(t), 'k', label='$s(t)$');\n",
+    "        ax.set_xlabel(r'$\\rightarrow t $ [s]')\n",
+    "        ax.set_xlim([-5,10]); ax.set_ylim([-.09, 1.19]); ax.legend(loc='upper right', bbox_to_anchor=(1.0, 0.8)); rwth_plt.axis(ax);\n",
+    "        \n",
+    "        axn = ax.twiny(); rwth_plt.stem(axn, [1], [1], color='rwth:blue', label=r'$s(n)$'); \n",
+    "        axn.set_xlabel(r'$\\rightarrow n$'); decorate_second_axis(axn);\n",
+    "        axn.legend(loc='upper right', bbox_to_anchor=(1.0, 0.6));\n",
+    "        axn.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)  \n",
+    "               \n",
+    "        ax = axs[1]; ax.set_title('Frequenzbereich');\n",
+    "        ax.plot(f, S(f), 'k', label='$S(f)$')\n",
+    "        ax.plot(f, Sa_f, 'rwth:red', label='$S_\\mathrm{a}(f)/r$')\n",
+    "        ax.set_xlabel(r'$\\rightarrow f $ [Hz]')\n",
+    "        ax.set_xlim([-5,10]);  ax.set_ylim([-.09, 1.19]); ax.legend(loc='upper right', bbox_to_anchor=(1.0, 0.8)); rwth_plt.axis(ax);\n",
+    "        \n",
+    "        axk = ax.twiny(); rwth_plt.stem(axk, [1], [1], label=r'$S_\\mathrm{d}(k)/r$');\n",
+    "        axk.set_xlabel(r'$\\rightarrow k$'); decorate_second_axis(axk);\n",
+    "        axk.legend(loc='upper right', bbox_to_anchor=(1.0, 0.5));\n",
+    "        axk.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style) \n",
+    "        \n",
+    "    else: # Update only the plot lines themselves. Should not take too much processing time\n",
+    "        ax = axs[0] # time domain\n",
+    "        axn = fig.axes[2];\n",
+    "        if r > 0:\n",
+    "            if len(axn.collections) > 1: axn.collections[1].remove()\n",
+    "            if M == 0:\n",
+    "                rwth_plt.stem_set_data(axn.containers[0], n, snT); #axn.containers[0].set_label(r'$s(n)$')\n",
+    "            else:\n",
+    "                rwth_plt.stem_set_data(axn.containers[0], n, sd);  #axn.containers[0].set_label(r'$s_\\mathrm{d}(n)$')\n",
+    "                axn.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)  \n",
+    "                        \n",
+    "        ax = axs[1] # frequency domain\n",
+    "        ax.lines[1].set_ydata(Sa_f)\n",
+    "        axk = fig.axes[3];\n",
+    "        if M > 0 and r > 0:\n",
+    "            if len(axk.collections) > 1: axk.collections[1].remove()\n",
+    "            axk.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)  \n",
+    "            rwth_plt.stem_set_data(axk.containers[0], k, SkF)\n",
+    "            \n",
+    "    axn = fig.axes[2]; axn.set_visible(r > 0)\n",
+    "    axk = fig.axes[3]; axk.set_visible(M > 0 and r > 0)\n",
+    "    if r > 0:\n",
+    "        axn.set_xlim((ax.get_xlim()[0]*r,ax.get_xlim()[1]*r))\n",
+    "        if M > 0:\n",
+    "            axk.set_xlim((ax.get_xlim()[0]*M/r,ax.get_xlim()[1]*M/r))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "---\n",
+    "This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources) (OER). Feel free to use the notebook for your own purposes. The code is licensed under the [MIT license](https://opensource.org/licenses/MIT). \n",
+    "\n",
+    "Please attribute the work as follows: \n",
+    "*Christian Rohlfing, Emin Kosar, Übungsbeispiele zur Vorlesung \"Grundgebiete der Elektrotechnik 3 - Signale und Systeme\"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University."
+   ]
+  }
+ ],
+ "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.8.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
+%% Cell type:code id: tags:
+
+``` python
+# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University
+%matplotlib widget
+
+from ipywidgets import interact, interactive, fixed, HBox, VBox
+import ipywidgets as widgets
+from IPython.display import clear_output, display, HTML
+
+import matplotlib.pyplot as plt
+from rwth_nb.plots import colors
+import rwth_nb.plots.mpl_decorations as rwth_plt
+from rwth_nb.misc.signals import *
+from rwth_nb.misc.transforms import *
+
+t,deltat = np.linspace(-10,10,50001,  retstep=True) # t-axis
+f,deltaf = np.linspace(-10,10,len(t), retstep=True) # f-axis
+kMax = 16 # number of k values in sum for Sa(f)
+```
+
+%% Cell type:markdown id: tags:
+
+<div>
+    <img src="figures/rwth_ient_logo@2x.png" style="float: right;height: 5em;">
+</div>
+
+# Demonstrator Diskrete Fourier-Transformation
+
+Zum Starten: Im Menü: Run <span class="fa-chevron-right fa"></span> Run All Cells auswählen.
+
+## Einleitung
+
+$s(t)$ wird ideal mit Abtastrate $r$ abgetastet.
+
+Diskrete Fourier-Transformation mit $F = \frac{1}{M}$
+$$
+S_\mathrm{d}(k)
+=
+\sum_{n=0}^{M-1} s_\mathrm{d}(n) \mathrm{e}^{-\mathrm{j}2\pi kFn} \quad k=0,\ldots,M-1 \text{ ,}
+$$
+wobei $s_\mathrm{d}(n)$ um $M$ periodisches, diskretes Zeitsignal mit $t = n \cdot \frac{1}{r}$.
+
+Für die Rücktransformation in den Zeitbereich gilt dann
+$$
+s_\mathrm{d}(n)
+=
+\frac{1}{M} \sum_{k=0}^{M-1}S_\mathrm{d}(k) \mathrm{e}^{\mathrm{j}2\pi kFn} \quad n=0,\ldots,M-1 \text{ ,}
+$$
+wobei
+$f = k \cdot \frac{r}{M}$.
+
+%% Cell type:code id: tags:
+
+``` python
+M = 24; # period of discrete signal
+r = 4; # sampling rate
+
+n = np.arange(M) # discrete time axis
+sd = tri(n/r) + tri((n-M)/r) # discrete signal sd(n)
+
+F = 1/M
+Sd = np.zeros_like(sd, dtype=np.complex)
+for k in np.arange(M): # very slow computation for teaching purpose! Fasten your seatbelts with FFT :)
+    Sd[k] = np.sum(sd * np.exp(-2j*np.pi*k*F*n))
+#Sd = np.fft.fft(sd)
+Sd = np.real(Sd) # discard super small imaginary part (1e-16)
+
+# Plot
+fig,axs = plt.subplots(2,1)
+ax = axs[0]; rwth_plt.stem(ax, n, sd);
+ax.set_xlabel(r'$\rightarrow n$'); ax.set_ylabel(r'$\uparrow s_\mathrm{d}(n)$')
+rwth_plt.axis(ax);
+
+ax = axs[1]; rwth_plt.stem(ax, np.arange(M), Sd/r);
+ax.set_xlabel(r'$\rightarrow k$'); ax.set_ylabel(r'$\uparrow S_\mathrm{d}(k)/r$')
+rwth_plt.axis(ax); fig.tight_layout()
+```
+
+%% Cell type:markdown id: tags:
+
+## Demonstrator
+
+%% Cell type:code id: tags:
+
+``` python
+def decorate_second_axis(ax):
+    ax.spines['top'].set_position(('axes',.9)); ax.spines['right'].set_color('none'); ax.spines['left'].set_color('none'); ax.spines['bottom'].set_color('none');
+    ax.xaxis.set_label_coords(.98,1); ax.xaxis.label.set_horizontalalignment('right')
+
+fig, axs = plt.subplots(2, 1, **rwth_plt.landscape); #fig.tight_layout()
+highlight_rect_style = {'facecolor': "none", 'edgecolor': 'rwth:orange', 'linestyle':'--', 'linewidth':2.0}
+
+@widgets.interact(  r = widgets.IntSlider(min=0, max=10, step=2, value=0, description='$r$ [Hz]'),
+                  Mbr = widgets.IntSlider(min=0, max=10, step=2, value=0, description='$M/r$'))
+def update_plots(r, Mbr):
+    M = r*Mbr
+    global n, snT, k, SkF
+    # Continuous functions s(t) and S(f)
+    s = tri
+    S = lambda f: si(np.pi*f) ** 2
+
+    # Ideal sampling
+    Sa_f = np.zeros_like(S(f))
+    nT = []; snT = nT; n = nT;
+    if r > 0:
+        T = 1/r
+
+        # Sampled function s(n)
+        nT, snT = sample(t, s(t), T); n = nT/T;
+
+        # Sampled continuous spectrum S_a(f)
+        for k in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements
+            Sa_f += S(f-k/T)
+        #Sa_f = Sa_f/T # normalize here
+
+        # Periodic copies of s(n)
+        kF = []; SkF = kF; k = kF;
+        sd = np.zeros_like(snT)
+        if M > 0:
+            F = 1/M # in Hz
+            kF, SkF = sample(f, Sa_f, r*F); k = kF/F/r
+            for d in np.arange(-kMax, kMax+1): # evaluate infinite sum only for 2*kMax+1 elements
+                sd += s((n-d*M)*T)
+    display('M={}'.format(M))
+
+    # Plot
+    if not axs[0].lines: # Call plot() and decorate axes. Usually, these functions take some processing time
+        ax = axs[0]; ax.set_title('Zeitbereich');
+        ax.plot(t, s(t), 'k', label='$s(t)$');
+        ax.set_xlabel(r'$\rightarrow t $ [s]')
+        ax.set_xlim([-5,10]); ax.set_ylim([-.09, 1.19]); ax.legend(loc='upper right', bbox_to_anchor=(1.0, 0.8)); rwth_plt.axis(ax);
+
+        axn = ax.twiny(); rwth_plt.stem(axn, [1], [1], color='rwth:blue', label=r'$s(n)$');
+        axn.set_xlabel(r'$\rightarrow n$'); decorate_second_axis(axn);
+        axn.legend(loc='upper right', bbox_to_anchor=(1.0, 0.6));
+        axn.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)
+
+        ax = axs[1]; ax.set_title('Frequenzbereich');
+        ax.plot(f, S(f), 'k', label='$S(f)$')
+        ax.plot(f, Sa_f, 'rwth:red', label='$S_\mathrm{a}(f)/r$')
+        ax.set_xlabel(r'$\rightarrow f $ [Hz]')
+        ax.set_xlim([-5,10]);  ax.set_ylim([-.09, 1.19]); ax.legend(loc='upper right', bbox_to_anchor=(1.0, 0.8)); rwth_plt.axis(ax);
+
+        axk = ax.twiny(); rwth_plt.stem(axk, [1], [1], label=r'$S_\mathrm{d}(k)/r$');
+        axk.set_xlabel(r'$\rightarrow k$'); decorate_second_axis(axk);
+        axk.legend(loc='upper right', bbox_to_anchor=(1.0, 0.5));
+        axk.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)
+
+    else: # Update only the plot lines themselves. Should not take too much processing time
+        ax = axs[0] # time domain
+        axn = fig.axes[2];
+        if r > 0:
+            if len(axn.collections) > 1: axn.collections[1].remove()
+            if M == 0:
+                rwth_plt.stem_set_data(axn.containers[0], n, snT); #axn.containers[0].set_label(r'$s(n)$')
+            else:
+                rwth_plt.stem_set_data(axn.containers[0], n, sd);  #axn.containers[0].set_label(r'$s_\mathrm{d}(n)$')
+                axn.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)
+
+        ax = axs[1] # frequency domain
+        ax.lines[1].set_ydata(Sa_f)
+        axk = fig.axes[3];
+        if M > 0 and r > 0:
+            if len(axk.collections) > 1: axk.collections[1].remove()
+            axk.fill_between([-.1, M-1+.1], -.05, 1.05, **highlight_rect_style)
+            rwth_plt.stem_set_data(axk.containers[0], k, SkF)
+
+    axn = fig.axes[2]; axn.set_visible(r > 0)
+    axk = fig.axes[3]; axk.set_visible(M > 0 and r > 0)
+    if r > 0:
+        axn.set_xlim((ax.get_xlim()[0]*r,ax.get_xlim()[1]*r))
+        if M > 0:
+            axk.set_xlim((ax.get_xlim()[0]*M/r,ax.get_xlim()[1]*M/r))
+```
+
+%% Cell type:markdown id: tags:
+
+---
+This notebook is provided as [Open Educational Resource](https://en.wikipedia.org/wiki/Open_educational_resources) (OER). Feel free to use the notebook for your own purposes. The code is licensed under the [MIT license](https://opensource.org/licenses/MIT).
+
+Please attribute the work as follows:
+*Christian Rohlfing, Emin Kosar, Übungsbeispiele zur Vorlesung "Grundgebiete der Elektrotechnik 3 - Signale und Systeme"*, gehalten von Jens-Rainer Ohm, 2020, Institut für Nachrichtentechnik, RWTH Aachen University.