{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Programmieraufgabe 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# some setup\n",
    "%matplotlib inline\n",
    "import numpy as np \n",
    "import matplotlib.pyplot as plt "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Trigonometrische Interpolation und Schnelle Fourier-Transformation\n",
    "\n",
    "Gegeben seien eine $2 \\pi$-periodische Funktion $f$ und $n$ äquidistante Stützstellen auf $[0,2 \\pi)$, d.h.\n",
    "\n",
    " \\begin{equation}\n",
    " (x_0,y_0),\\dots,(x_{n-1},y_{n-1})\n",
    " \\end{equation}\n",
    " mit\n",
    " \\begin{align}\n",
    " x_k &= 2 \\pi \\frac{k}{n} \\\\\n",
    " y_k &= f(x_k) \\ ,\n",
    " \\end{align}\n",
    " für $k \\in \\{0,\\dots,n-1\\}$.\n",
    " \n",
    " In der Vorlesung wurde gezeigt, dass das trigonometrische Polynom\n",
    " \\begin{equation}\n",
    " p(x) = \\frac{1}{n} \\sum^{n-1}_{l=0} \\beta_l \\, e^{i l x}\n",
    " \\end{equation}\n",
    "die Interpolationsaufgabe $p(x_k)=y_k$ genau dann löst, wenn gilt:\n",
    "\\begin{equation}\n",
    " \\beta_k = \\sum^{n-1}_{l=0} y_l \\, e^{-i k x_l} = \\sum^{n-1}_{l=0} y_l \\, e^{-i 2\\pi \\frac{kl}{n}} =: \\bigl(A \\, y \\bigr)_k\n",
    " \\end{equation}\n",
    "Beachten Sie, dass der Normierungsfaktor $\\frac{1}{n}$ gegenüber der Vorlesung nach $p(x)$ geschoben wurde, um mit der Konvention von _numpy_ übereinzustimmen.\n",
    "\n",
    "Den Vektor $\\beta$ nennt man auch die _Diskrete Fourier-Transformierte_ des Vektors $y$.\n",
    "\n",
    "\n",
    "a) Implementieren Sie eine Funktionen, die die Diskrete Fourier-Transformierte eines Vektors berechnet, indem sie die Fourier-Koeffizientenmatrix $A$ aufstellt und anschließend eine Matrix-Vektor-Multiplikation durchführt. Verfahren Sie analog für die inverse DFT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# a)\n",
    "def calcDFT(y):\n",
    "    n = len(y)\n",
    "    # Die imaginäre Einheit in python ist 'j', die zugehörige Zahlendarstellung 0.+1.j oder kurz 1j\n",
    "    # 'fromfunction' erzeugt ein Array durch Funktionsauswertung, 'lambda' ist die syntax für eine anonyme Funktion\n",
    "    \n",
    "    # Fourier-Matrix\n",
    "    ### INSERT CODE ###\n",
    "    \n",
    "    # Matrix-Vektor-Multiplikation\n",
    "    return ### INSERT CODE ###\n",
    "\n",
    "def calcInvDFT(y):\n",
    "    n = len(y)\n",
    "    \n",
    "    # inverse Fourier-Matrix\n",
    "    ### INSERT CODE ###\n",
    "    \n",
    "    # Matrix-Vektor-Multiplikation\n",
    "    return ### INSERT CODE ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) Implementieren Sie (von Hand) die in der Vorlesung gezeigte _Schnelle Fourier-Transformation (FFT)_ nach _Cooley und Tukey_ in eine Funktion, die als Eingabe einen Vektor $y$ nimmt, und dessen Diskrete Fourier-Transformierte zurückgibt. Erstellen Sie eine zweite Funktion, die die inverse Transformation übernimmt. Beachten Sie, dass die Länge $n$ des Eingangsvektors eine Zweierpotenz sein muss."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# b)\n",
    "def calcFFT(y):\n",
    "    # Rekursiver FFT-Algorithmus\n",
    "    n = len(y)\n",
    "    if n==1:\n",
    "        return y\n",
    "    else:\n",
    "        # Prüft in jeder Instanz, ob n (noch) gerade ist\n",
    "        ### INSERT CODE ###\n",
    "        \n",
    "        # Auftrennen in gerade und ungerade Indizes\n",
    "        # y[a:b:s] bedeutet Index-Range von a bis b-1 in Schrittweite s, [::] ist kurz für [Anfang:Ende:1]\n",
    "        ### INSERT CODE ###\n",
    "        \n",
    "        # Zusammenfügen der Teilergebnisse (durch numpy's broadcasting spart man sich for-Schleifen)\n",
    "        ### INSERT CODE ###\n",
    "                \n",
    "        return c\n",
    "    \n",
    "def calcInvFFT(y):\n",
    "    n = len(y)\n",
    "    if n==1:\n",
    "        return y\n",
    "    else:\n",
    "        \n",
    "        ### INSERT CODE ###\n",
    "        \n",
    "        return c/2 # Normierung!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c) Machen Sie sich mit dem Modul _numpy.fft_ vertraut, insbesondere mit den Routinen _fft_ und _ifft_.\n",
    "\n",
    "d) Betrachten Sie nun die Funktion $f(t) = \\cos(10 t) + \\cos(5 t)$. Erstellen Sie ein Array $y_k = f(x_k)$ für $n=16$ und wenden Sie mit jeder der Methoden aus Teilaufgabe a)-c) die Diskrete Fourier Transformation an.\n",
    "Gilt jeweils $\\text{iDFT}(\\text{DFT}(y)) = y$ ?  \n",
    "Vergleichen Sie die Methoden untereinander. Erhalten Sie jeweils das gleiche Ergebnis oder gibt es numerische Abweichungen? Ermitteln Sie den Laufzeit-Unterschied für verschiedene (große) $n$ mittels _%timeit_ ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# d)\n",
    "\n",
    "n = 16\n",
    "\n",
    "# Erzeuge Stützstellen (z.B. np.linspace)\n",
    "### INSERT CODE ###\n",
    "\n",
    "# Gib die DFT aus und prüfe die Inverse\n",
    "print('Naive Lösung:')\n",
    "### INSERT CODE ###\n",
    "print()\n",
    "\n",
    "print('Eigene FFT:')\n",
    "### INSERT CODE ###\n",
    "print()\n",
    "\n",
    "print('Numpy FFT:')\n",
    "### INSERT CODE ###\n",
    "print()\n",
    "\n",
    "# Vergleiche calcFFT und np.fft.fft:\n",
    "print('FFT = FFT?')\n",
    "### INSERT CODE ###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# immernoch d)\n",
    "n=2048\n",
    "\n",
    "# Erzeuge Stützstellen\n",
    "### COPY CODE FROM ABOVE ###\n",
    "\n",
    "# Prüfe Laufzeit (siehe %timeit in der ipython Dokumentation)\n",
    "print(\"Zeit für naive DFT:\")\n",
    "### INSERT CODE ###\n",
    "print()\n",
    "print(\"Zeit für eigene FFT:\")\n",
    "### INSERT CODE ###\n",
    "print()\n",
    "print(\"Zeit für numpy FFT:\")\n",
    "### INSERT CODE ###"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "e ) Konstruieren Sie aus den $\\beta_k$ die Koeffizienten $a_k,b_k$ der reellen Darstellung $q(t)$ des Interpolationspolynoms, siehe dazu S. 169 f. des [alternativen Vorlsungsskripts](http://www.ins.uni-bonn.de/teaching/vorlesungen/AlMaWS13/script.pdf) (beschränken Sie sich auf den Fall gerader $n$). Plotten Sie schließlich die Funktion $f(t)$ sowie das resultierende Interpolationspolynom $q(t)$ für verschiedene $n$ auf ganz $[0,2 \\pi)$. Was fällt Ihnen auf? Wie sieht es aus für z.B. $g(t) = 2\\, \\, \\bigl| \\frac{t}{\\pi}-1\\,\\bigr| -1 $?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# e)\n",
    "n = 16\n",
    "\n",
    "# Funktion f\n",
    "f = lambda t:\n",
    "# Alternative Funktion\n",
    "#f = lambda t:\n",
    "\n",
    "# berechne FFT (mit n sampling points)\n",
    "### INSERT CODE ###\n",
    "\n",
    "\n",
    "m = int(n/2)\n",
    "# Berechne a_k und b_k\n",
    "### INSERT CODE ###\n",
    "\n",
    "# konstruiere Interpolationspolynom\n",
    "q = lambda t:\n",
    "# or def q(t):\n",
    "### INSERT CODE ###\n",
    "\n",
    "# vorhergehende Funktion lässt sich nicht ohne weiteres auf ein numpy-array broadcasten. Erzeuge eine Universal-function:\n",
    "Uq = np.frompyfunc(q,1,1)\n",
    "\n",
    "\n",
    "\n",
    "#konstruiere hochaufgelöstes Gitter für pyplot (mit N sampling points)\n",
    "N = 1000\n",
    "r = np.linspace(0,2*np.pi,N,endpoint=False)\n",
    "# Werte f und q auf diesem Gitter aus\n",
    "yf = f(r)\n",
    "yq = Uq(r)\n",
    "\n",
    "# Plotte Funktion und Interpolierende mit pyplot\n",
    "plt.plot( r, yf, '-b')\n",
    "plt.plot( r, yq, '-r')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "f) Fügen Sie den Daten der Funktion $f$ für $n=1024$ Rauschen hinzu, mittels\n",
    "\\begin{equation}\n",
    "\\bar{y}_k = y_k + X_k \\ , \\quad X_k \\sim U(-\\delta,\\delta) \\, .\n",
    "\\end{equation}\n",
    "Berechnen Sie die Fourier-Transformierte und filtern Sie das Rauschen heraus, indem Sie alle Einträge mit $|\\bar{\\beta}_k| < \\epsilon \\cdot n$ für ein $\\epsilon>0$ auf Null setzen. Führen Sie anschließend die Rücktransformation durch.\n",
    "Plotten Sie die verrauschten sowie die bereinigten Daten für verschiedene $\\delta$ und $\\epsilon$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# f)\n",
    "#parameter\n",
    "delta = 0.3\n",
    "epsilon = 0.1\n",
    "\n",
    "\n",
    "# konstruiere hochaufgelöstes Gitter\n",
    "n = 1024\n",
    "r = np.linspace(0,2*np.pi,n,endpoint=False)\n",
    "\n",
    "# Funktion f\n",
    "### COPY CODE FROM ABOVE ###\n",
    "\n",
    "# konstruiere verrauschtes Signal\n",
    "### INSERT CODE ###\n",
    "yf =\n",
    "\n",
    "# konstruiere entrauschtes Signal\n",
    "### INSERT CODE ###\n",
    "yfC = \n",
    "\n",
    "\n",
    "# Plotte beide Signale\n",
    "plt.plot( r, yf, '-b')\n",
    "plt.show()\n",
    "plt.plot( r, yfC, '-r')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}