IUE Rayleigh-Verteilung.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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 scipy import signal # convolution, \n",
    "\n",
    "from ient_nb.ient_plots import *\n",
    "from ient_nb.ient_signals import *\n",
    "from ient_nb.ient_audio import *\n",
    "\n",
    "# Gaussian PDF\n",
    "gaussian_pdf = lambda x, ms, sigmas: 1/(sigmas*np.sqrt(2*np.pi)) * np.exp(-(x-ms)**2/(2*sigmas**2))\n",
    "\n",
    "fs = 22050; # sampling rate\n",
    "(t, deltat) = np.linspace(-10, 10, 20*fs, retstep=True) # t axis in seconds\n",
    "Amax = 15 # limit of histogram (the total range would be from -Amax to Amax)\n",
    "(x, deltax) = np.linspace(-Amax, Amax, 4001, retstep=True) # x axis (corresponds to amplitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div>\n",
    "    <img src=\"figures/rwth_ient_logo@2x.png\" style=\"float: right;height: 5em;\">\n",
    "</div>\n",
    "\n",
    "# Rayleigh-Verteilung\n",
    "\n",
    "## Einleitung\n",
    "\n",
    "Ein Gauß-verteilert Rauschprozess mit Leistungsdichte $N_0$ wird über einen Bandpass $h_\\mathrm{BP}(t)$ mit Mittenfrequenz $f_0$ und Bandbreite $f_\\Delta$ übertragen.\n",
    "\n",
    "![Blockdiagramm](figures/rayleigh_block_diagram.png)\n",
    "\n",
    "Die Tiefpass-Hüllkurve $|n_\\mathrm{T}|(t)$ ist in diesem Fall Rayleigh-verteilt."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rauschprozess $n_0(t)$\n",
    "\n",
    "Gauß-verteiles Rauschen $n_0(t)$ mit Leistungsdichte $N_0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Sample one realization of n0(t)\n",
    "N0 = 0.0004\n",
    "deltaf = 1/fs\n",
    "sigma_n0_squared = N0/deltaf # Variance of n_0(t)\n",
    "sigma_n0 = np.sqrt(sigma_n0_squared)\n",
    "n0 = np.random.normal(0, sigma_n0, len(t))\n",
    "\n",
    "# Estimate p_n0(x) with histogram\n",
    "pn0_hist, bins = np.histogram(n0, bins=50, range=(-Amax, Amax), density=True)\n",
    "x0_hist = (bins[:-1] + bins[1:])/2 # x-axis\n",
    "\n",
    "# Plot time domain signal and both calculated and measured PDFs\n",
    "fig, axs = plt.subplots(2, 1); \n",
    "ax = axs[0]\n",
    "ax.plot(1000*t, n0, 'rwth'); \n",
    "ax.set_xlabel(r'$\\rightarrow t$ [ms]', bbox=ient_wbbox); ax.set_ylabel(r'$\\uparrow {}^k n_0(t)$'); \n",
    "ax.set_xlim([-21,21]); ax.set_ylim([-11,11]); ient_axis(ax);\n",
    "\n",
    "ax = axs[1]\n",
    "ient_stem(ax, x0_hist, pn0_hist, 'black-50', markerfmt=\" \", label='Geschätzt');\n",
    "ax.plot(x, gaussian_pdf(x, 0, sigma_n0), 'rwth', label='Berechnet'); \n",
    "ax.set_xlabel(r'$\\rightarrow x$'); ax.set_ylabel(r'$\\uparrow p_{n_0}(x)$'); ax.legend(); ient_axis(ax);\n",
    "\n",
    "ient_audio_play(n0, fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bandpassrauschen $n(t)$\n",
    "\n",
    "Erzeuge Filterkoeffizienten eines Bandpasses $h_\\mathrm{BP}(t)$ mit Mittenfrequenz $f_0 = 1000\\,\\mathrm{Hz}$ und Bandbreite $f_\\Delta=200\\,\\mathrm{Hz}$ (Butterworth-Filter 4. Ordnung)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def butter_bandpass(f0, fdelta, fs, order=5):\n",
    "    nyq = 0.5 * fs # Nyquist frequency\n",
    "    normal_min = (f0-fdelta/2) / nyq # normalized \n",
    "    normal_max = (f0+fdelta/2) / nyq # normalized \n",
    "    b, a = signal.butter(order, [normal_min,normal_max], btype='bandpass', analog=False)\n",
    "    return b, a\n",
    "\n",
    "# Filter requirements.\n",
    "order = 4 # order :)\n",
    "f0 = 1000  # center frequency\n",
    "fdelta = 200 # band width\n",
    "\n",
    "b, a = butter_bandpass(f0, fdelta, fs, order) # generate filter coefficients\n",
    "_, H = signal.freqz(b, a, worN=int(len(t)/2+1)) # compute H(f)=H(z=e^(j 2 pi f)) out of b, a\n",
    "\n",
    "f = np.linspace(0, fs/2, len(H));\n",
    "\n",
    "# Plot\n",
    "fig,ax = plt.subplots(1,1); ax.plot(f, np.abs(H), 'rwth'); \n",
    "ax.axvline(f0, color='k', linestyle='--',lw=0.5); # cutoff frequency\n",
    "ax.set_xlabel(r'$\\rightarrow f$ [Hz]'); ax.set_ylabel(r'$\\uparrow |H_\\mathrm{BP}(f)|$'); \n",
    "ax.set_xlim([0,5500]); ax.set_ylim([-.25,1.19]); ient_axis(ax); "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$n(t)$ ist das sog. \"Bandpassrauschen\"\n",
    "\n",
    "![Blockdiagramm](figures/rayleigh_block_diagram.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def butter_filter(s, b, a):\n",
    "    g = signal.lfilter(b, a, s)\n",
    "    return g\n",
    "\n",
    "# Filter n0(t) to yield n(t)\n",
    "n = butter_filter(n0, b, a)\n",
    "\n",
    "# Plot time domain signals\n",
    "fig, axs = plt.subplots(2, 1); \n",
    "ax = axs[0]; ax.plot(1000*t, n0, 'rwth'); \n",
    "ax.set_xlabel(r'$\\rightarrow t$ [ms]', bbox=ient_wbbox); ax.set_ylabel(r'$\\uparrow {}^k n_0(t)$'); \n",
    "ax.set_xlim([0,110]); ax.set_ylim([-11,11]); ient_axis(ax);\n",
    "\n",
    "ax = axs[1]; ax.plot(1000*t, n, 'rwth'); \n",
    "ax.set_xlabel(r'$\\rightarrow t$ [ms]', bbox=ient_wbbox); ax.set_ylabel(r'$\\uparrow {}^k n(t)$'); \n",
    "ax.set_xlim([0,110]); ax.set_ylim([-1.9,1.9]); ient_axis(ax);\n",
    "\n",
    "ient_audio_play(n0, fs, r'${}^k n_0(t)$:');\n",
    "ient_audio_play(n, fs, r'${}^k n(t)$:');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Das Bandpassrauschen $n(t)$ ist ebenfalls Gauß-verteilt mit Varianz $\\sigma_n^2=2N_0 f_\\Delta$ (Formel (9.24))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate standard deviation and variance\n",
    "sigma_n_squared = 2*N0*fdelta # Equation (9.24)\n",
    "sigma_n = np.sqrt(sigma_n_squared)\n",
    "\n",
    "print(\"Varianz berechnet: {:.5f}, gemessen: {:.5f}\".format(sigma_n_squared, np.var(n)))\n",
    "\n",
    "# Calculate PDF\n",
    "pn = gaussian_pdf(x, 0, sigma_n)\n",
    "\n",
    "# Estimate PDF\n",
    "pn_hist,bins = np.histogram(n, bins=50, range=(-5,5), density=True)\n",
    "x_hist = (bins[:-1] + bins[1:]) / 2\n",
    "\n",
    "# Plot\n",
    "fig,axs = plt.subplots(2,1); \n",
    "ax = axs[0]; \n",
    "ient_stem(ax, x0_hist, pn0_hist, 'black-50', markerfmt=\" \", label='Geschätzt');\n",
    "ax.plot(x, gaussian_pdf(x, 0, sigma_n0), 'rwth', label='Berechnet'); \n",
    "ax.set_xlabel(r'$\\rightarrow x$'); ax.set_ylabel(r'$\\uparrow p_{n_0}(x)$'); \n",
    "ax.legend(); ient_axis(ax);\n",
    "\n",
    "ax = axs[1]; \n",
    "ient_stem(ax, x_hist, pn_hist, 'black-50', markerfmt=\" \");\n",
    "ax.plot(x, pn, 'rwth'); \n",
    "ax.set_xlabel(r'$\\rightarrow x$'); ax.set_ylabel(r'$\\uparrow p_n(x)$'); ient_axis(ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tiefpass-Hüllkurve $|n_\\mathrm{T}(t)|$\n",
    "\n",
    "Das (komplexwertige) äquivalente Tiefpasssignal (bezüglich Mittenfrequenz $f_0$) lässt sich über die analytische Komponente ${}^k n_+(t)$ berechnen:\n",
    "$${}^k n_\\mathrm{T}(t) = 2 {}^k n_+(t) \\mathrm{e}^{-\\mathrm{j}2\\pi f_0 t}$$\n",
    "mit analytischer Komponente $s_+(t)$ eines Signals $s(t)$:  $\\mathrm{Re}\\{s_+(t)\\} = \\frac{1}{2} s(t)$ \n",
    "und $\\mathrm{Im}\\{s_+(t)\\} = \\frac{1}{2} s(t) \\ast \\frac{1}{\\pi t}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_eq_lp(s, f0): # s: BP signal, f0: center frequency\n",
    "    # signal.hilbert is NOT calculating the Hilbert transform, \n",
    "    # but the analytic component s_plus(t) using the Hilbert transform\n",
    "    s_plus = 0.5*signal.hilbert(s) \n",
    "    sT = 2*s_plus*np.exp(-2*1J*np.pi*f0*t)\n",
    "    return sT\n",
    "\n",
    "# Equivalent lowpass signal of n(t)\n",
    "nT = calculate_eq_lp(n, f0)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(1, 1); \n",
    "ax.plot(t*1000, n, 'rwth', label=r'${}^k n(t)$'); \n",
    "ax.plot(t*1000, np.abs(nT), 'grun', label=r'$|{}^k n_\\mathrm{T}(t)|$'); \n",
    "ax.set_xlim([0,110]); ax.set_xlabel(r'$\\rightarrow t$ [ms]', bbox=ient_wbbox); \n",
    "ax.legend(); ient_axis(ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Einhüllende $|{}^k n_\\mathrm{T}(t)|$ ist Rayleigh-verteilt mit Verteilungsdichtefunktion (Formel (9.27)):\n",
    "$$p_{|n_\\mathrm{T}|}(x) = \\epsilon(x) \\frac{x}{\\sigma_n^2} \\mathrm{e}^{- \\frac{x^2}{2\\sigma_n^2}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotte ${}^k n(t)$, Einhüllende $|{}^k n_\\mathrm{T}(t)|$ sowie Verteilungsdichtefunktionen (Abbildung 9.10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate PDF\n",
    "rayleigh_pdf = lambda x, sigma: unitstep(x)*x/sigma**2 * np.exp(-x**2/(2*sigma**2))\n",
    "pnTabs = rayleigh_pdf(x, sigma_n)\n",
    "\n",
    "# Estimate PDF of |n_T(t)| with histogram\n",
    "pnTabs_hist, bins_pnTabs = np.histogram(np.abs(nT), bins=50, range=(0,5), density=True)\n",
    "xT_hist = (bins_pnTabs[:-1] + bins_pnTabs[1:]) / 2\n",
    "\n",
    "# Plot\n",
    "fig, axs = plt.subplots(1, 2, figsize=(8,4), gridspec_kw = {'width_ratios':[2, 1]}); \n",
    "ax = axs[0];\n",
    "ax.plot(t*1000, n, 'rwth', label=r'${}^k n(t)$'); \n",
    "ax.plot(t*1000, np.abs(nT), 'grun', label=r'$|{}^k n_\\mathrm{T}(t)|$'); \n",
    "ax.plot(t*1000, -np.abs(nT), 'grun', lw=0.5, linestyle='--'); \n",
    "ax.set_xlabel(r'$\\rightarrow t$ [ms]', bbox=ient_wbbox); \n",
    "ax.set_xlim([0,220]); ax.legend(); ient_axis(ax);\n",
    "\n",
    "ax = axs[1]; \n",
    "ax.plot(pn, x, 'rwth', label=r'$p_n(x)$'); \n",
    "ax.plot(pnTabs, x, 'grun', label=r'$p_{|n_\\mathrm{T}|}(x)$'); \n",
    "ax.set_ylabel(r'$\\uparrow x$'); \n",
    "ax.set_ylim(axs[0].get_ylim()); ax.legend(); ient_axis(ax);\n",
    "ax.text(0, sigma_n, r'$\\sigma_n$', color='grun', fontsize=12, bbox=ient_wbbox);\n",
    "ax.plot([0, np.max(pnTabs)], [sigma_n, sigma_n], '--', color='grun', lw=0.5);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotte Verteilungsdichtefunktionen erneut, diesmal mit Schätzungen"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "fig,axs = plt.subplots(2,1); \n",
    "ax = axs[0]; \n",
    "ient_stem(ax, x_hist, pn_hist, 'black-50', markerfmt=\" \", label='Geschätzt');\n",
    "ax.plot(x, pn, 'rwth'); \n",
    "ax.set_ylabel(r'$\\uparrow p_n(x)$'); ax.set_xlabel(r'$\\rightarrow x$'); \n",
    "ax.set_xlim([-5,5]); ax.set_ylim([0,1.05]); ient_axis(ax);\n",
    "\n",
    "ax = axs[1]; \n",
    "ient_stem(ax, xT_hist, pnTabs_hist, 'black-50', markerfmt=\" \", label='Geschätzt');\n",
    "ax.plot(x, pnTabs, 'grun', label='Berechnet'); \n",
    "ax.set_xlabel(r'$\\rightarrow x$'); ax.set_ylabel(r'$\\uparrow p_{|n_\\mathrm{T}|}(x)$'); \n",
    "ax.set_xlim([-5,5]); ax.set_ylim([0,1.75]); ax.legend(); ient_axis(ax);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Verteilungsfunktion der Rayleigh-Verteilung\n",
    "$$P_{|n_\\mathrm{T}|}(x) = \\int\\limits_{-\\infty}^x  p_{|n_\\mathrm{T}|}(\\xi) \\mathrm{d}\\xi = \\big[1- \\mathrm{e}^{- \\frac{x^2}{2\\sigma_n^2}}\\big] \\epsilon(x)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate PF\n",
    "rayleigh_pf = lambda x, sigma: unitstep(x)*(1-np.exp(-x**2/(2*sigma**2)))\n",
    "PnTabs = rayleigh_pf(x, sigma_n)\n",
    "\n",
    "# Estimate PF\n",
    "PnTabs_hist = np.cumsum(pnTabs_hist)*np.mean(np.diff(bins_pnTabs))\n",
    "\n",
    "# Plot\n",
    "fig,axs = plt.subplots(2, 1); \n",
    "ax = axs[0]\n",
    "ient_stem(ax, xT_hist, PnTabs_hist, 'black-50', markerfmt=\" \", label='Geschätzt');\n",
    "ax.plot(x, PnTabs, 'grun', label='Berechnet'); \n",
    "ax.set_xlabel(r'$\\rightarrow x$'); ax.set_ylabel(r'$\\uparrow P_{|n_\\mathrm{T}|}(x)$'); \n",
    "ax.set_xlim([-5,5]); ax.set_ylim([0,1.19]); ax.legend(); ient_axis(ax);\n",
    "\n",
    "ax = axs[1]    \n",
    "ient_stem(ax, xT_hist, pnTabs_hist, 'black-50', markerfmt=\" \", label='Geschätzt');\n",
    "ax.plot(x, pnTabs, 'grun', label='Berechnet'); \n",
    "ax.set_xlabel(r'$\\rightarrow x$'); ax.set_ylabel(r'$\\uparrow p_{|n_\\mathrm{T}|}(x) = \\frac{\\mathrm{d}}{\\mathrm{d} x}P_{|n_\\mathrm{T}|}(x)$'); \n",
    "ax.set_xlim([-5,5]); ax.set_ylim([0,1.8]); ient_axis(ax);\n",
    "ax.text(sigma_n, -0.01, r'$\\sigma_n$', color='rot', fontsize=12, bbox=ient_wbbox, verticalalignment='top', horizontalalignment='center');\n",
    "ax.plot([sigma_n, sigma_n], [0, np.max(pnTabs)], '--', color='rot', lw=1, );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mittelwert $m_{|n_\\mathrm{T}|}$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_nTabs = np.sqrt(np.pi/2)*sigma_n\n",
    "\n",
    "(m_nTabs, np.mean(np.abs(nT)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Leistung $L_{|n_\\mathrm{T}|}^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L_n = sigma_n_squared # m_n = 0\n",
    "L_nTabs = 2*L_n\n",
    "\n",
    "(L_nTabs, np.mean(np.abs(nT)**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Streuung $\\sigma_{|n_\\mathrm{T}|}^2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma_nTabs_squared = sigma_n_squared * (2-np.pi/2)\n",
    "\n",
    "(sigma_nTabs_squared, np.var(np.abs(nT)))"
   ]
  },
  {
   "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
}