{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "hide_input": false
   },
   "outputs": [],
   "source": [
    "# Copyright 2019 Institut für Nachrichtentechnik, RWTH Aachen University\n",
    "%matplotlib notebook\n",
    "\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import clear_output, display, HTML\n",
    "\n",
    "import scipy.special as special # erfc\n",
    "from scipy import signal # convolution\n",
    "\n",
    "from ient_nb.ient_plots import *\n",
    "from ient_nb.ient_signals import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "    <img src=\"figures/rwth_ient_logo@2x.png\" style=\"float: right;height: 5em;\">\n",
    "</div>\n",
    "\n",
    "# Amplitudenabtastung\n",
    "\n",
    "## Übersicht\n",
    "\n",
    "![Blockdiagramm](figures/amplitudenabtastung_block_diagram.png)\n",
    "\n",
    "Signal nach Korrelationsfilter $h(t)$: $$y(t)=g(t)+n_\\mathrm{e}(t)=m(t)\\ast h(t) + n(t) \\ast h(t)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hilfsfunktionen\n",
    "\n",
    "Faltung sowie Konvertierung String zu Bits und umgekehrt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def convolution(s, h):\n",
    "    # Convolve s and h numerically\n",
    "    g = signal.convolve(s, h, mode='full')*deltat; g = g[0:len(s)];\n",
    "    return g\n",
    "\n",
    "def str2bits(s):\n",
    "    return np.array([int(x) for x in list(''.join(format(ord(i),'b').zfill(8) for i in s))])\n",
    "\n",
    "def bits2str(b):\n",
    "    try: # decode string\n",
    "        s=''.join(np.array([chr(i) for i in np.packbits(b)]))\n",
    "    except: # decoding failed\n",
    "        s = '-'\n",
    "    return s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sender\n",
    "\n",
    "Streng geheime Nachricht wird in binäres Signal $a_n$ codiert:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = str2bits('Institut für Nachrichtentechnik')\n",
    "La = len(a) # Anzahl an Bits\n",
    "a[0:10] # Ausgabe der ersten 10 Werte"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameter und Achsen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameter\n",
    "fs = 1000 # Samplingrate\n",
    "T=1 # Breite des Trägersignals in Sekunden\n",
    "Tmax = La*T # Gesamtlänge des Sendesignals\n",
    "\n",
    "# Achsen\n",
    "(t,deltat) = np.linspace(0, La*T, La*T*fs, retstep=True) # Zeitachse\n",
    "nT = np.linspace(1,La*T,La*T) # nT-Achse\n",
    "nT_idxs = (nT*fs).astype(int); nT_idxs[-1] = nT_idxs[-1]-1 # nT in Samples"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trägersignal $s(t)=\\mathrm{rect}\\left(\\frac{t}{T}-\\frac{1}{2}\\right)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Trägersignal\n",
    "s = lambda t: rect(t/T-0.5)\n",
    "#s = lambda t: rect(2*(t-0.25))-rect(2*(t-0.75))\n",
    "    \n",
    "# Sendeenergie wenn a_n=1\n",
    "Es = T #  s(t) Rechteck der Breite T!\n",
    "    \n",
    "# Plot Trägersignal s(t)\n",
    "fig,ax = plt.subplots(1,1); plt.plot(t/T, s(t), 'rwth'); \n",
    "ax.set_xlabel(r'$\\rightarrow t/T$'); ax.set_ylabel(r'$\\uparrow s(t)$');\n",
    "ax.set_xlim([0,2]); ax.set_ylim([-0.1,1.1]); ient_axis(ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sendesignal $m(t)=\\sum_{n=-\\infty}^\\infty a_n s(t-nT)$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = np.zeros(len(t)) # Sendesignal\n",
    "for l in range(len(a)):\n",
    "    m = m + a[l]*s(t-l*T)\n",
    "    \n",
    "# Plot\n",
    "fig,ax = plt.subplots(1,1); plt.plot(t/T, m, 'rwth');\n",
    "ax.set_xlabel(r'$\\rightarrow t/T$'); ax.set_ylabel(r'$\\uparrow m(t)$');\n",
    "ax.set_xlim([-0.1,11.9]); ax.set_ylim(0,1.09); ient_axis(ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kanal\n",
    "\n",
    "$n(t)$ Gauß-verteiltes additives weißes Rauschen mit Leistungsdichte $N_0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_channel(N0, l=len(m)):\n",
    "     # Konvertierung Leistungsdichte in Leistung (im Diskreten). Einheit Leistungsdichte: Ws=W/Hz, Einheit Leistung: W\n",
    "    deltaf = 1/fs\n",
    "    sigma_n_squared = N0/deltaf # Varianz von n(t)\n",
    "    \n",
    "    # n(t) Werte sind Gauß-verteilt\n",
    "    n = np.random.normal(0, np.sqrt(sigma_n_squared), l)\n",
    "    \n",
    "    return n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Augenblicksleistung des gefilterten Signals $n_\\mathrm{e}(t)=n(t)\\ast g(t)$: \n",
    "$$N=\\mathcal{E}\\{n_\\mathrm{e}^2(T)\\} = N_0 \\int\\limits_{-\\infty}^\\infty |h(t)|^2 \\mathrm{d}t$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_N(N0,h):\n",
    "    return N0*sum(h**2)*deltat\n",
    "\n",
    "def calculate_N0(N,h):\n",
    "    return N/(sum(h**2)*deltat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Empfänger\n",
    "\n",
    "Filter $h(t) = s(T-t)$ mit $s(t)=\\mathrm{rect}\\left(\\frac{t}{T}-\\frac{1}{2}\\right)$ und \n",
    "$$y(t)=g(t)+n_\\mathrm{e}(t)=m(t)\\ast h(t) + n(t) \\ast h(t)$$\n",
    "\n",
    "Zunächst wird $g(t) = m(t) \\ast h(t)$ statisch berechnet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Korrelationsfilter\n",
    "h = s(T-t)\n",
    "g = convolution(m,h); # Faltung: g(t) = m(t) * h(t)\n",
    "\n",
    "# Plot\n",
    "fig,ax = plt.subplots(1,1)\n",
    "plt.plot(t/T, m, 'rwth'); plt.plot(t/T, g, 'grun')\n",
    "ax.set_ylabel(r'$\\uparrow m(t),g(t)$'); ax.set_xlabel(r'$\\rightarrow t/T$');\n",
    "ax.set_xlim([-0.1,11.9]); ax.set_ylim(0,1.09); plt.legend({r'$m(t)$',r'$g(t)$'},loc='upper right'); ient_axis(ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Da später die Störleistungsdichte $N_0$ des Kanals variiert werden soll, wird die Faltung $n_\\mathrm{e}(t) = n(t) \\ast h(t)$ in einer Funktion berechnet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def simulate_receiver(n,h,g,C):\n",
    "    # Faltung mit Filter h(t)\n",
    "    # eigentlich y(t) = [m(t)+n(t)]*h(t); allerdings ist m(t) determiniert und g(t) vorher berechnet worden\n",
    "    ne = convolution(n, h);\n",
    "    y = g + ne\n",
    "    \n",
    "    # Abtaster\n",
    "    yT = y[nT_idxs]; gT = np.round(g[nT_idxs]); neT = ne[nT_idxs];\n",
    "    \n",
    "    # Entscheidungsstufe\n",
    "    a_dec = np.array([0 if yTn < C else 1  for yTn in yT]) # Entscheidung: y(nT) > C: an=1, sonst: an=0\n",
    "    str_dec = bits2str(a_dec) # Decodiere String aus Bits\n",
    "    \n",
    "    return (str_dec, a_dec, y, ne, yT,gT,neT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sendeleistung und -energie:\n",
    "\n",
    "- Augenblicksleistung $S_\\mathrm{a} = \\mathcal{E}\\{g^2(T)\\}=g^2(T)$\n",
    "- Energie (wenn $a_n=1$, ansonsten $E_s=0$!): $E_s = \\int\\limits_{-\\infty}^{\\infty}|s(t)|^2\\mathrm{d}t$\n",
    "- Energie pro gesendetem Bit: $E_b = \\frac{E_s}{2}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Sa = 1 # Augenblickssendeleistung: Sa = g(T)^2=1 wenn a_n=1\n",
    "\n",
    "Eb = Es/2 # Energie pro gesendetem Bit\n",
    "C = np.sqrt(Sa)/2 # Entscheidungsschwelle\n",
    "\n",
    "# Die gleichen Werte, diesmal nur mit numerischer Ungenauigkeit\n",
    "gT = g[nT_idxs];\n",
    "Sa_measured = np.mean(gT[a==1]**2) \n",
    "Es_measured = sum(np.abs(s(t))**2)*deltat # Integral -> Summe im Diskreten; dt -> deltat (Abstand zwischen zwei Samplen)\n",
    "Eb_measured = Es_measured/2\n",
    "C_measured = np.sqrt(Sa_measured)/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistiken\n",
    "\n",
    "Bitfehlerwahrscheinlichkeit\n",
    "$$ P_b = \\frac{1}{2}\\mathrm{erfc}\\left(\\sqrt{\\frac{E_b}{4N_0}}\\right) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_prob(Eb,N0):\n",
    "    Pb =  0.5*special.erfc(np.sqrt(Eb/(4*N0)))\n",
    "    return Pb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interaktive Demo: Unipolare Binärübertragung"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(fig,axs) = plt.subplots(2, 2, figsize=(8,6))\n",
    "\n",
    "interact_manual.opts['manual_name'] = 'Update'\n",
    "@interact_manual(N0 = widgets.FloatSlider(min=0.001, max=.1, step=0.001, value=0.01, description='$N_0$'))\n",
    "def update_signals(N0=0.01):\n",
    "    # Simulation\n",
    "    n = simulate_channel(N0) # Simuliere Kanal\n",
    "    str_dec, a_dec, y, ne, yT,gT,neT = simulate_receiver(n,h,g,C) # Simuliere Empfänger\n",
    "    \n",
    "    # Plot\n",
    "    if not axs[0,1].lines: # Plot\n",
    "        ax=axs[0,0]; ax.plot(t/T, m, 'rwth');\n",
    "        ax.set_ylabel(r'$\\uparrow m(t)$'); ax.set_xlabel(r'$\\rightarrow t/T$'); ax.set_ylim(0,1.09);\n",
    "\n",
    "        ax=axs[0,1]; ax.plot(t/T, m+n, 'rwth');\n",
    "        ax.set_ylabel(r'$\\uparrow m(t)+n(t)$'); ax.set_xlabel(r'$\\rightarrow t/T$');\n",
    "        \n",
    "        ax=axs[1,0]; ax.plot(t/T, g, 'rwth'); \n",
    "        ax.plot(nT/T, gT, color='rwth', marker='o', lw=0)\n",
    "        ax.set_ylabel(r'$\\uparrow g(t)$'); ax.set_xlabel(r'$\\rightarrow t/T$'); ax.set_ylim(-1,2.09);\n",
    "\n",
    "        ax=axs[1,1]; ax.plot(t/T, y, 'rwth'); \n",
    "        ax.plot(nT/T, yT, color='rwth', marker='o', lw=0)\n",
    "        ax.set_ylabel(r'$\\uparrow y(t)$'); ax.set_xlabel(r'$\\rightarrow t/T$'); ax.set_ylim(-1,2.09);\n",
    "        ax.plot(ax.get_xlim(), [C, C], 'k--',lw=1); ax.text(0, C+0.1, r'$C$', fontsize=12, horizontalalignment='left',verticalalignment='bottom',backgroundcolor='w'); \n",
    "    \n",
    "        [ax.set_xlim([-0.1,11.9]) and ient_axis(ax) for ax in axs.flatten()];\n",
    "    else: # update lines\n",
    "        ax = axs[0,1]; ax.lines[0].set_ydata(m+n)        \n",
    "        ax = axs[1,1]; ax.lines[0].set_ydata(y); ax.lines[1].set_ydata(yT);\n",
    "        \n",
    "    # Leistungen, S/N Verhältnis\n",
    "    N = calculate_N(N0,h); SN = Sa/N;\n",
    "    Pb = calculate_prob(Eb, N0);\n",
    "    N_measured = np.mean(neT**2); SN_measured = Sa_measured/N_measured;  \n",
    "    Pb_measured = np.mean(a_dec != a);\n",
    "        \n",
    "    # Output\n",
    "    display(HTML(\n",
    "       \"\"\"<table><tr><td>TX a</td>{}</tr><tr><td>RX a</td>{}</tr></table><br />Empfangen: {}<br />\n",
    "            Sa/N={:.2f} dB (berechnet: {:.2f} dB), log10(Pb)={:.4f} (berechnet: {:.4f})\"\"\".format(\n",
    "               '<td>{}</td>'.format('</td><td>'.join(str(tmp) for tmp in a)),\n",
    "               '<td>{}</td>'.format('</td><td>'.join(str(tmp) for tmp in a_dec)), str_dec, \n",
    "           10*np.log10(SN_measured), 10*np.log10(SN), np.log10(Pb_measured+np.finfo(float).eps), np.log10(Pb), \n",
    "       )\n",
    "    ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wahrscheinlichkeitsdichtefunktion\n",
    "\n",
    "Mit Verteilungsdichtefunktion von $n_\\mathrm{e}(nT)$ \n",
    "$$p_{n_\\mathrm{e}}(x) = \\frac{1}{\\sqrt{2\\pi N}} \\exp\\left[ - \\frac{x^2}{2N} \\right]$$\n",
    "und VDF von $g(nT)$\n",
    "$$p_g(t) = \\mathrm{Prob}\\{a_n=0\\} \\cdot \\delta(x) + \\mathrm{Prob}\\{a_n=1\\} \\cdot \\delta\\left(x-\\sqrt{S_a}\\right)$$\n",
    "folgt VDF von $y(nT) = g(nT) + n_\\mathrm{e}(nT)$ \n",
    "$$p_y(x) = p_g(x) \\ast p_{n_\\mathrm{e}}(x) \n",
    "= \\frac{1}{\\sqrt{2 \\pi N}} \\left( \n",
    "\\mathrm{Prob}\\{a_n=0\\} \\cdot \\exp\\left[ -\\frac{x^2 }{2N} \\right] + \n",
    "\\mathrm{Prob}\\{a_n=1\\} \\cdot \\exp\\left[ -\\frac{\\left(x-\\sqrt{S_\\mathrm{a}}\\right)^2}{2N} \\right] \\right)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_pdf(x,a,N,mean0,mean1):\n",
    "    # a_n = 0\n",
    "    pa0 = np.sum(a==0)/len(a); \n",
    "    py0 = 1/np.sqrt(2*np.pi*N)*np.exp(-(x-mean0)**2/(2*N))\n",
    "    # a_n = 1\n",
    "    pa1 = np.sum(a==1)/len(a);\n",
    "    py1 = 1/np.sqrt(2*np.pi*N)*np.exp(-(x-mean1)**2/(2*N))\n",
    "    # Gesamtwahrscheinlichkeit\n",
    "    py = pa0*py0 + pa1*py1\n",
    "    return py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Demo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(logSN,deltaSN) = np.linspace(-5,25,1000, retstep=True) # logSN axis\n",
    "(fig_pdf,axs_pdf) = plt.subplots(2, 1, figsize=(8,8))\n",
    "np.seterr(divide = 'ignore') \n",
    "@interact(N0 = widgets.FloatSlider(min=0.003, max=.5, step=0.001, value=0.005, \n",
    "                                   description=r'$N_0$', readout_format='.3f', continuous_update=False))\n",
    "def update_signals_pdf(N0=0.01):\n",
    "    # Simulation\n",
    "    n = simulate_channel(N0)  # Simuliere Kanal\n",
    "    str_dec, a_dec, y, ne, yT,gT,neT = simulate_receiver(n,h,g,C)  # Simuliere Empfänger\n",
    "        \n",
    "    # Leistungen, S/N Verhältnis\n",
    "    N = calculate_N(N0,h); SN = Sa/N;\n",
    "    Pb = calculate_prob(Eb, N0)\n",
    "    N_measured = np.mean(neT**2); SN_measured = Sa_measured/N_measured;\n",
    "    Pb_measured = np.mean(a_dec != a)\n",
    "    \n",
    "    # Verteilungsdichte p_y(x)\n",
    "    py_measured, bins = np.histogram(yT,bins=100,range=(-0.5,1.5),density=True)  # gemessen mit Histogramm\n",
    "    x = (bins[:-1] + bins[1:]) / 2 \n",
    "    py = calculate_pdf(x,a,N,0,np.sqrt(Sa))  # berechnet\n",
    "    \n",
    "    # Plot\n",
    "    if not axs_pdf[0].lines:\n",
    "        ax = axs_pdf[0]; #ax.plot(x, py_measured, 'rwth'); \n",
    "        ient_stem(ax, x, py_measured, 'rwth', markerfmt=\" \")\n",
    "        ax.plot(x, py, 'grun')\n",
    "        ax.set_xlabel(r'$\\rightarrow x$'); ax.set_ylabel(r'$\\uparrow p_y(x)$'); ient_axis(ax);\n",
    "                \n",
    "        ax = axs_pdf[1];\n",
    "        ax.semilogy(logSN, 0.5*special.erfc(np.sqrt(10**(logSN/10)/8)), color=\"rwth\")\n",
    "        ax.semilogy(10*np.log10(SN_measured), Pb_measured , marker='o', markersize=5, color=\"rot\");\n",
    "        ax.set_xlabel(r'$\\rightarrow E_b/N_0$ [dB]'); ax.set_ylabel(r'$\\rightarrow P_{b}$');\n",
    "        ax.text(ax.get_xlim()[1]-1, .24, r'$0{,}5$', fontsize=12, \n",
    "                horizontalalignment='right',verticalalignment='top',backgroundcolor='w');\n",
    "        ax.plot(ax.get_xlim(), [0.5, 0.5], 'k--',lw=1); ax.grid(); ax.set_ylim([1e-10,2]);\n",
    "    else: # update lines\n",
    "        ax = axs_pdf[0]; \n",
    "        ient_stem_set_data(ax.containers[0], x, py_measured)\n",
    "        ax.lines[-1].set_data(x, py); \n",
    "        ax = axs_pdf[1]; ax.lines[1].set_xdata(10*np.log10(SN_measured+eps)); ax.lines[1].set_ydata(Pb_measured);        \n",
    "        \n",
    "    # Output\n",
    "    display(HTML('gemessen: log10(Pb)={:.4f} (berechnet: {:.4f})'.format(\n",
    "        np.log10(Pb_measured),np.log10(Pb))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 \"Informationsübertragung\"*, gehalten von Jens-Rainer Ohm, 2019, 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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}