GDET3 Autokorrelationsfunktion.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University\n",
    "%matplotlib widget\n",
    "\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact, interactive\n",
    "\n",
    "import scipy as sp\n",
    "import scipy.special\n",
    "import scipy.signal\n",
    "\n",
    "from ient_nb.ient_plots import *\n",
    "from ient_nb.ient_signals import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "    <img src=\"ient_nb/figures/rwth_ient_logo@2x.png\" style=\"float: right;height: 5em;\">\n",
    "</div>\n",
    "\n",
    "# Auto- und Kreuzkorrelationsfunktion von Energiesignalen\n",
    "## Einleitung\n",
    "\n",
    "### Energiesignal\n",
    "\n",
    "Ein Signal $s(t)$ heißt Energiesignal, wenn seine Signalenergie endlich ist, es gilt also\n",
    "$$E_s = \\int\\limits_{-\\infty}^\\infty |s(t)|^2 \\mathrm{d} t < \\infty \\text{ .}$$  \n",
    "\n",
    "Viele wichtige Signale haben keine endliche Gesamtenergie, z. B. alle periodischen Signale, die Sprungfunktion oder zeitlich nicht begrenzte Zufallssignale. Für Signale mit Dirac-Impulsen ist die Energie nicht definiert.\n",
    "\n",
    "### Kreuzkorrelationsfunktion\n",
    "Für zwei Energiesignale $s(t)$ und $g(t)$ kann die *Kreuzkorrelationsfunktion* berechnet werden. Diese zeigt an, wie ähnlich sich zwei Signale bei unterschiedlichen Verschiebungen $\\tau$ sind.\n",
    "$$\n",
    "\\varphi_{sg}^\\mathrm{E}(\\tau) \n",
    "= \\int\\limits_{-\\infty}^\\infty s^\\ast(t) g(t+\\tau)\\mathrm{d}t \n",
    "\\stackrel{t\\rightarrow -t}{=} \\int\\limits_{-\\infty}^\\infty s(-t)^\\ast g(\\tau-t) \\mathrm{d} t \n",
    "= s^\\ast(-\\tau) \\ast g(\\tau)\n",
    "= \\left[\\varphi_{gs}^\\mathrm{E}(-\\tau)\\right]^\\ast\n",
    "$$\n",
    "Die Kreuzkorrelationsfunktionen $\\varphi_{sg}^\\mathrm{E}(\\tau)$ und $\\varphi_{gs}^\\mathrm{E}(\\tau)$ sind also zueinander zeitgespiegelt und konjugiert-komplex.\n",
    "Für die Berechnung der Kreuzkorrelationsfunktion müssen die beiden Signale nicht zwingend Energiesignale sein, es ist ausreichend, dass das Integral berechnet werden kann. Man beachte hier auch die Ähnlichkeit zur Faltung zweier Signale, bei der Berechnung kann entsprechend vorgegangen werden. \n",
    "### Autokorrelationsfunktion\n",
    "\n",
    "Die Autokorrelationsfunktion entspricht der Kreuzkorrelationsfunktion, wenn $g(t) = s(t)$ gilt. Sie zeigt an, wie ähnlich ein Signal sich selber bei einer zeitlichen Verschiebung $\\tau$ ist\n",
    "$$\n",
    "\\varphi_{ss}^\\mathrm{E}(\\tau) \n",
    "= \\int\\limits_{-\\infty}^\\infty s^\\ast(t) s(t+\\tau)\\mathrm{d}t \\text{.}\n",
    "$$\n",
    "Die Autokorrelationsfunktion besitzt folgende allgemeingültige Eigenschaften:\n",
    "* Die Autokorrelationsfunktion ist immer eine gerade Funktion.\n",
    "* Den maximalen Wert nimmt eine Autokorrelationsfunktion für $\\tau=0$ an, da in diesem Fall größte Ähnlichkeit vorliegt. \n",
    "* Das Maximum der Autokorrelationsfunktion eines Energiesignals ist gleich seiner Energie\n",
    "$$\\varphi_{ss}^\\mathrm{E}(0)=E_s \\text{.}$$\n",
    "* Bei zeitlich begrenzten Signalen hat die Autokorrelationsfunktion die doppelte Breite des Signals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interaktive Demo\n",
    "Die interaktive Demo ermöglicht es, sich die Kreuz- und Autokorrelation für verschiedene Signale anzusehen. \n",
    "$$\n",
    "\\varphi_{sg}^\\mathrm{E}(\\tau) \n",
    "= s^\\ast(-\\tau) \\ast g(\\tau) \n",
    "= g(\\tau) \\ast s^\\ast(-\\tau)\n",
    "= \\int\\limits_{-\\infty}^\\infty g(t) s^\\ast(t-\\tau) \\mathrm{d} t \n",
    "$$\n",
    "\n",
    "Es ist auch möglich, eine eigene Funktion $s_0(t)$ zu definieren. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_0 = lambda t: rect(t/2-1/2)*(-t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def correlation(s, g):\n",
    "    # Correlate s and g numerically\n",
    "    return sp.signal.convolve(np.conj(s(-t)),g(t), mode='full')*deltat\n",
    "\n",
    "fs = 2000 # Samplingrate\n",
    "(t,deltat) = np.linspace(-8,8, 16*fs, retstep=True) # Zeitachse\n",
    "(tau,deltatau) = np.linspace(-16,16, 2*len(t)-1, retstep=True) # Korrelation\n",
    "\n",
    "signal_types = {'Rechteck'           : rect,\n",
    "                'Dreieck'            : tri, \n",
    "                'Sprungfunktion'     : unitstep, \n",
    "                'si-Funktion'        : lambda t: si(t*np.pi), \n",
    "                'Exponentialimpuls'  : lambda t: unitstep(t)*np.exp(-t),\n",
    "                'Gauß-Signal'        : gauss, \n",
    "                'Doppelrechteck'     : lambda t: rect(t*2+0.5)-rect(t*2-0.5),\n",
    "                'Rampe'              : lambda t: t*rect(t-0.5), \n",
    "                'Versch. Rechteck'   : lambda t: -rect(t-0.5),\n",
    "                'Eigene Kreation'    : s_0,\n",
    "               }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wähle Signale für $s(t)$ und $g(t)$ im Drop-Down-Menü aus. Für beide Signale kann auch die Breite $T$ und die Verschiebung $t_0$ angepasst werden. Die jeweiligen Signale werden dann in der nachfolgenden Abbildung angezeigt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "fig0, axs0 = plt.subplots(1, 2, figsize=(8,2))\n",
    "@widgets.interact(s_type=widgets.Dropdown(options=list(signal_types.keys()), description=r'Wähle $s(t)$:'),\n",
    "                  s_T=widgets.FloatSlider(min=0.5, max=4, value=1, step=.1, description=r'Dehnung T', style=ient_wdgtl_style), \n",
    "                  s_t0=widgets.FloatSlider(min=-2, max=2, value=0, step=.1, description=r'Verschiebung $t_0$', style=ient_wdgtl_style), \n",
    "                  g_type=widgets.Dropdown(options=list(signal_types.keys()), description=r'Wähle $g(t)$:'),\n",
    "                  g_T=widgets.FloatSlider(min=0.5, max=4, value=1, step=.1, description=r'Dehnung T', style=ient_wdgtl_style), \n",
    "                  g_t0=widgets.FloatSlider(min=-2, max=2, value=0, step=.1, description=r'Verschiebung $t_0$', style=ient_wdgtl_style))\n",
    "def update_signals(s_type, s_T, s_t0, g_type, g_T, g_t0):\n",
    "    global s, g, phi_sg # reused in second interactive plot\n",
    "    s = lambda t: signal_types[s_type]((t-s_t0)/s_T);\n",
    "    g = lambda t: signal_types[g_type]((t-g_t0)/g_T);\n",
    "    phi_sg = correlation(s, g) # numerical correlation\n",
    "    \n",
    "    if not axs0[0].lines: # plot s(t) and g(t)\n",
    "        ax = axs0[0]; ax.plot(t, s(t), 'rwth');\n",
    "        ax.set_xlabel(r'$\\rightarrow t$'); ax.set_ylabel(r'$\\uparrow s(t)$')\n",
    "        ax.set_xlim([-2.9, 2.9]); ax.set_ylim([-1.19, 1.19]);  ient_axis(ax); ient_grid(ax);\n",
    "        \n",
    "        ax = axs0[1]; ax.plot(t, g(t), 'rwth');\n",
    "        ax.set_xlabel(r'$\\rightarrow t$'); ax.set_ylabel(r'$\\uparrow g(t)$')\n",
    "        ax.set_xlim(axs0[0].get_xlim()); ax.set_ylim(axs0[0].get_ylim()); ient_axis(ax); ient_grid(ax);\n",
    "    else: # update lines\n",
    "        axs0[0].lines[0].set_ydata(s(t)); \n",
    "        axs0[1].lines[0].set_ydata(g(t));\n",
    "    try: # if convolution figure is already opened, update s(t)\n",
    "        if axs[0].lines:\n",
    "            axs[0].lines[0].set_ydata((g(t)));\n",
    "            ient_update_ylim(axs[0], np.concatenate((g(t), s(t))), 0.19, 5); ient_update_ylim(axs[1], phi_sg, 0.19, 5);\n",
    "            update_plot(-2); # update correlation plot\n",
    "    except: pass\n",
    "    ient_update_ylim(axs0[0], s(t), 0.19, 5);  ient_update_ylim(axs0[1], g(t), 0.19, 5);    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In der nachfolgenden Grafik kann für die beiden oben ausgewählten Funktionen nun das Ergebnis der Kreuzkorrelationsfunktion betrachtet werden. Hierfür muss der Schieberegler von links nach rechts geschoben werden. Wenn das Kästchen bei Integrand angeklickt ist, wird die aktuell überlappende Fläche der beiden Funktionen gestrichelt angezeigt. Diese entspricht dem Wert der Kreuzkorrelationsfunktion. Dies kann im unteren Teil der Grafik interaktiv betrachtet werden. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false,
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "fig, axs=plt.subplots(2, 1, figsize=(8,6),) # gridspec_kw = {'width_ratios':[3, 1]}\n",
    "@widgets.interact(tau_shft=widgets.FloatSlider(min=-4, max=4, value=-2, step=.1, description=r'Verschiebung $\\tau$', style=ient_wdgtl_style), \n",
    "                  show_integrand=widgets.Checkbox(value=True, description='Zeige Integrand', style=ient_wdgtl_style))\n",
    "def update_plot(tau_shft, show_integrand=True):\n",
    "    tau_ind = np.where(tau>=tau_shft); tau_ind = tau_ind[0][0]; phi_plot = phi_sg.copy(); phi_plot[tau_ind:] = np.nan; # hide g(t') with t'>t\n",
    "    sg = np.conj(s(t-tau_shft))*g(t) # integrand\n",
    "    \n",
    "    if not axs[0].lines: # Call plot() and decorate axes. Usually, these functions take some processing time\n",
    "        ax = axs[0]; ax.plot(t, g(t), 'rwth', label=r'$g(t)$'); # plot g(t)\n",
    "        ax.plot(t, np.conj(s(t-tau_shft)), 'grun', label=r'$s^\\ast(t-\\tau)$'); # plot s(t-tau)\n",
    "        ax.plot(t, sg, '--', color='orange', lw=1, label=r'$g(t)s^\\ast(t-\\tau)$'); # plot integrand\n",
    "        ient_annotate_xtick(ax, r'$\\tau$', tau_shft, -0.1, 'rwth', 15); # mark t on tau axis\n",
    "        ax.fill_between(t, 0, sg, facecolor=\"none\", hatch=\"//\", edgecolor='k', linewidth=0.0); # hatch common area\n",
    "        ax.set_xlabel(r'$\\rightarrow t$');\n",
    "        ax.set_xlim([-4.2,4.2]); ient_update_ylim(ax, np.concatenate((g(t), s(t))), 0.19, 5);\n",
    "        ax.legend(); ient_grid(ax); ient_axis(ax);\n",
    "        \n",
    "        ax = axs[1]; ax.plot(tau, phi_plot); # plot phi_sg(tau)\n",
    "        ax.plot([tau_shft, tau_shft], [0, phi_sg[tau_ind]], 'ko--', lw=1);\n",
    "        ax.set_xlabel(r'$\\rightarrow \\tau$'); \n",
    "        ax.set_ylabel(r'$\\uparrow \\varphi_{sg}^\\mathrm{E}(\\tau) = \\int g(t) s^\\ast(t-\\tau)\\mathrm{d}t$'); \n",
    "        ax.set_xlim(axs[0].get_xlim()); ient_update_ylim(ax, phi_sg, 0.19, 5); \n",
    "        ient_axis(ax); ient_grid(ax); fig.tight_layout(); \n",
    "        \n",
    "    else: # Replace only xdata and ydata since plt.plot() takes longer time\n",
    "        ax = axs[0]; ax.lines[1].set_ydata(np.conj(s(t-tau_shft))); ax.lines[2].set_ydata(sg); # update signals\n",
    "        ax.texts[0].set_x(tau_shft); ax.lines[3].set_xdata([tau_shft,tau_shft]) # update labels\n",
    "        if ax.collections: ax.collections[0].remove(); # update integrand\n",
    "        if show_integrand: ax.fill_between(t,0,sg, facecolor=\"none\", hatch=\"//\", edgecolor='k', linewidth=0.0);\n",
    "        ax = axs[1]; ax.lines[0].set_ydata(phi_plot); # update signals\n",
    "        ax.lines[1].set_xdata([tau_shft, tau_shft]); ax.lines[1].set_ydata([0, phi_sg[tau_ind]]); # update labels\n",
    "\n",
    "    axs[0].lines[2].set_visible(show_integrand)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aufgaben\n",
    "* Wähle zwei gleiche Funktionen aus. Überprüfe die oben angegebenen Eigenschaften der Autokorrelationsfunktion. Variiere hierbei die Funktionen und betrachte die Autokorrelationsfunktion verschiedener Signale. \n",
    "* Wie ändert sich das Ergebnis, wenn die Breite eines der beiden Signale sich ändert?\n",
    "* Wähle zwei verschiedene Signale aus und beobachte, wo in diesen Fällen die maximale Kreuzkorrelation auftritt. Kannst du erklären, warum?"
   ]
  },
  {
   "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, Ü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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}