GDET3 Autokorrelationsfunktion.ipynb 13.2 KB
Newer Older
1
2
3
4
5
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
Christian Rohlfing's avatar
Christian Rohlfing committed
6
   "metadata": {
7
8
9
10
    "hide_input": false,
    "jupyter": {
     "source_hidden": true
    }
Christian Rohlfing's avatar
Christian Rohlfing committed
11
   },
12
13
   "outputs": [],
   "source": [
Christian Rohlfing's avatar
Christian Rohlfing committed
14
    "# Copyright 2020 Institut für Nachrichtentechnik, RWTH Aachen University\n",
15
    "%matplotlib widget\n",
16
17
    "\n",
    "import ipywidgets as widgets\n",
18
    "from ipywidgets import interact, interactive\n",
19
20
21
    "\n",
    "import scipy as sp\n",
    "import scipy.special\n",
22
    "import scipy.signal\n",
23
    "\n",
24
25
    "from ient_nb.ient_plots import *\n",
    "from ient_nb.ient_signals import *"
26
27
28
29
30
31
32
33
34
35
   ]
  },
  {
   "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",
Christian Rohlfing's avatar
Christian Rohlfing committed
36
37
38
    "# Auto- und Kreuzkorrelationsfunktion von Energiesignalen\n",
    "## Einleitung\n",
    "\n",
39
40
    "### Energiesignal\n",
    "\n",
41
42
    "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",
43
    "\n",
44
    "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",
45
    "\n",
46
    "### Kreuzkorrelationsfunktion\n",
47
    "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",
Christian Rohlfing's avatar
Christian Rohlfing committed
48
49
50
51
52
53
    "$$\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",
54
    "$$\n",
55
56
    "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",
57
58
    "### Autokorrelationsfunktion\n",
    "\n",
59
60
61
62
63
64
65
66
67
68
69
    "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."
70
71
72
73
74
75
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
Christian Rohlfing's avatar
Christian Rohlfing committed
76
77
78
79
80
81
82
83
    "## 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",
84
    "\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
85
    "Es ist auch möglich, eine eigene Funktion $s_0(t)$ zu definieren. "
86
87
88
89
90
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
91
   "metadata": {},
92
   "outputs": [],
Christian Rohlfing's avatar
Christian Rohlfing committed
93
94
95
96
97
98
99
100
101
102
103
104
105
   "source": [
    "s_0 = lambda t: rect(t/2-1/2)*(-t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "outputs": [],
106
   "source": [
Christian Rohlfing's avatar
Christian Rohlfing committed
107
    "def correlation(s, g):\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
108
    "    # Correlate s and g numerically\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
109
    "    return sp.signal.convolve(np.conj(s(-t)),g(t), mode='full')*deltat\n",
110
    "\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
111
112
113
    "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",
114
    "\n",
115
    "signal_types = {'Rechteck'           : rect,\n",
116
117
118
119
120
    "                '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",
121
    "                'Doppelrechteck'     : lambda t: rect(t*2+0.5)-rect(t*2-0.5),\n",
122
123
    "                'Rampe'              : lambda t: t*rect(t-0.5), \n",
    "                'Versch. Rechteck'   : lambda t: -rect(t-0.5),\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
124
    "                'Eigene Kreation'    : s_0,\n",
125
    "               }"
126
127
128
   ]
  },
  {
129
   "cell_type": "markdown",
130
131
   "metadata": {},
   "source": [
Christian Rohlfing's avatar
Christian Rohlfing committed
132
    "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."
133
134
135
   ]
  },
  {
136
137
   "cell_type": "code",
   "execution_count": null,
Christian Rohlfing's avatar
Christian Rohlfing committed
138
   "metadata": {
139
140
141
142
    "hide_input": false,
    "jupyter": {
     "source_hidden": true
    }
Christian Rohlfing's avatar
Christian Rohlfing committed
143
   },
144
   "outputs": [],
145
   "source": [
146
    "fig0, axs0 = plt.subplots(1, 2, figsize=(8,2))\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
147
148
149
150
151
152
153
    "@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",
154
    "    global s, g, phi_sg # reused in second interactive plot\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
155
156
    "    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",
157
    "    phi_sg = correlation(s, g) # numerical correlation\n",
158
    "    \n",
Christian Rohlfing's avatar
Christian Rohlfing committed
159
    "    if not axs0[0].lines: # plot s(t) and g(t)\n",
160
161
162
163
    "        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",
Christian Rohlfing's avatar
Christian Rohlfing committed
164
165
    "        ax = axs0[1]; ax.plot(t, g(t), 'rwth');\n",
    "        ax.set_xlabel(r'$\\rightarrow t$'); ax.set_ylabel(r'$\\uparrow g(t)$')\n",
166
    "        ax.set_xlim(axs0[0].get_xlim()); ax.set_ylim(axs0[0].get_ylim()); ient_axis(ax); ient_grid(ax);\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
167
    "    else: # update lines\n",
168
    "        axs0[0].lines[0].set_ydata(s(t)); \n",
Christian Rohlfing's avatar
Christian Rohlfing committed
169
170
171
    "        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",
172
    "            axs[0].lines[0].set_ydata((g(t)));\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
173
    "            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",
174
    "            update_plot(-2); # update correlation plot\n",
175
    "    except: pass\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
176
177
178
    "    ient_update_ylim(axs0[0], s(t), 0.19, 5);  ient_update_ylim(axs0[1], g(t), 0.19, 5);    "
   ]
  },
179
  {
Christian Rohlfing's avatar
Christian Rohlfing committed
180
   "cell_type": "markdown",
181
182
   "metadata": {},
   "source": [
183
    "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. "
184
185
   ]
  },
186
187
188
  {
   "cell_type": "code",
   "execution_count": null,
Christian Rohlfing's avatar
Christian Rohlfing committed
189
   "metadata": {
190
191
192
193
    "hide_input": false,
    "jupyter": {
     "source_hidden": true
    }
Christian Rohlfing's avatar
Christian Rohlfing committed
194
   },
195
196
   "outputs": [],
   "source": [
197
    "fig, axs=plt.subplots(2, 1, figsize=(8,6),) # gridspec_kw = {'width_ratios':[3, 1]}\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
198
    "@widgets.interact(tau_shft=widgets.FloatSlider(min=-4, max=4, value=-2, step=.1, description=r'Verschiebung $\\tau$', style=ient_wdgtl_style), \n",
199
    "                  show_integrand=widgets.Checkbox(value=True, description='Zeige Integrand', style=ient_wdgtl_style))\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
200
201
    "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",
Christian Rohlfing's avatar
Christian Rohlfing committed
202
    "    sg = np.conj(s(t-tau_shft))*g(t) # integrand\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
203
    "    \n",
204
205
206
    "    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",
Christian Rohlfing's avatar
Christian Rohlfing committed
207
    "        ax.plot(t, sg, '--', color='orange', lw=1, label=r'$g(t)s^\\ast(t-\\tau)$'); # plot integrand\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
208
209
210
    "        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",
211
212
    "        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",
Christian Rohlfing's avatar
Christian Rohlfing committed
213
    "        \n",
214
    "        ax = axs[1]; ax.plot(tau, phi_plot); # plot phi_sg(tau)\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
215
    "        ax.plot([tau_shft, tau_shft], [0, phi_sg[tau_ind]], 'ko--', lw=1);\n",
Christian Rohlfing's avatar
Christian Rohlfing committed
216
217
    "        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",
218
219
    "        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",
Christian Rohlfing's avatar
Christian Rohlfing committed
220
221
    "        \n",
    "    else: # Replace only xdata and ydata since plt.plot() takes longer time\n",
222
223
224
    "        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",
Christian Rohlfing's avatar
Christian Rohlfing committed
225
    "        if show_integrand: ax.fill_between(t,0,sg, facecolor=\"none\", hatch=\"//\", edgecolor='k', linewidth=0.0);\n",
226
227
    "        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",
Christian Rohlfing's avatar
Christian Rohlfing committed
228
    "\n",
229
    "    axs[0].lines[2].set_visible(show_integrand)"
Christian Rohlfing's avatar
Christian Rohlfing committed
230
231
232
233
   ]
  },
  {
   "cell_type": "markdown",
234
235
236
237
238
239
240
241
242
243
   "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",
244
245
   "metadata": {},
   "source": [
246
    "---\n",
247
248
249
    "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",
250
    "*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."
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
   ]
  }
 ],
 "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",
Christian Rohlfing's avatar
Christian Rohlfing committed
270
   "version": "3.8.1"
271
272
273
  }
 },
 "nbformat": 4,
274
 "nbformat_minor": 4
275
}