{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Wavelettransformation zur Rauschreduktion in Audiosignalen\n",
    "\n",
    "Dieser Beleg ist eine Zusammenarbeit von Johannes Müller und Willy Müller. Er wurde für die Hochschule Anhalt, Fachbereich 6 bearbeitet und erstellt. Für den erfolgreichen Abschluss des Moduls Signaltheorie sollte ein tangierendes Thema für einen Beleg ausgesucht werden, welches die Kenntnisse der Studierenden in den gelehrten Themen reflektiert.\n",
    "\n",
    "## Einleitung und Motivation\n",
    "\n",
    "Wir haben uns daher gefragt: Inwiefern kann die Wavelettransformation als Mittel zur Reduktion von Hintergrundgeräuschen in Audiosignalen verwendet werden? In der Zeitdomäne ist das sogennante *Ducking* eine häufig verwendete Methode und in der Frequenzdomäne wird oft das Spektrum des Störsignals von dem des Eingangssignals abgezogen. Mit Wavelets könnte eine Kombination dieser Effekte erreicht werden.\n",
    "\n",
    "Um diese Hypothese zu überprüfen, müssen zuerst die Rahmenbedingungen gesetzt werden.\n",
    "\n",
    "## Analyse der Problemstellung\n",
    "\n",
    "Es soll ein System erstellt werden, dass ein zeitdiskretes Eingangssignal $x(t)$ entgegennimmt und ein Signal $y(t)$ ausgibt. Die Amplitude der Signale ist über das Interval $$-1 \\leq x \\leq 1; -1 \\leq y \\leq 1$$ begrenzt. Die Diskretisierung in der Zeitachse ist auf die übliche Rate von $\\mathsf{sr} = 22050 \\ \\mathrm{Hz}$ festgelegt, so dass Audiofrequenzen bis ca. $10 \\ \\mathrm{kHz}$ abgebildet werden können. Die Dauer des Eingangssignals ist in folgeden Tests auf $T = 5 \\ \\mathrm{s}$ festgelegt, wobei das System dieses Signal in noch kleineren zeitlichen Abschnitten bearbeiten wird und somit auch ein zeitlich unbeschränktes Signal verarbeiten kann.\n",
    "\n",
    "Das Eingangssignal $x(t)$ setzt sich hier beispielsweise aus dem Nutzsignal $source(t)$ und dem weißen Rauschen $r(t)$ zusammen, wobei das Rauschen nur zu bestimmten Segmenten des Nutzsignals addiert wird.\n",
    "Zusätzlich zum Eingangssignal $x(t)$ steht dem System das Referenzsample $r(t)$ zur Verfügung, welches zur Bestimmung des Pegels des Störsignals herangezogen werden kann."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 294060 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Erzeugung des Störsignals\n",
    "\n",
    "import numpy\n",
    "sr = 22050 # sample rate\n",
    "T = 5.0    # seconds\n",
    "t = numpy.linspace(0, T, int(T*sr), endpoint=False) # time variable\n",
    "r = numpy.random.normal(loc=0, scale=1e-2, size=int(T*sr)) # white noise\n",
    "\n",
    "from IPython.display import Audio\n",
    "import librosa\n",
    "import librosa.display\n",
    "Audio(r, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 294060 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Erzeugung des Eingangssignals\n",
    "\n",
    "# Quelle: https://ccrma.stanford.edu/workshops/mir2014/audio/\n",
    "source, _ = librosa.load('audio/T39-piano.wav', sr=sr, duration=T)\n",
    "x = source + r * (numpy.cos(t*2) < 0)\n",
    "Audio(x, rate=sr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zeitliche Unterteilung\n",
    "\n",
    "Das Eingangssignal wird nun in kleinere überlappende Abschnitte unterteilt, welche anschließend zur Rauschminderung weiterverarbeitet werden. Die Länge der Abschnitte so gewählt, dass für die minimalste Frequenz des Signals mindestens eine volle Sinuskurve im Abschnitt enthalten ist. Sie darf aber auch nicht zu lang sein, denn in der Echtzeitverarbeitung müssten alle enthaltenen Werte gepuffert werden, was zusätzliche Latenz herbeibringt. Schließlich sollte sie eine volle Zweierpotenz betragen, um die Abschnitte in optimierten Algorithmen, wie folgend FFT- oder Waveletalgorithmen, effizient verarbeitet zu können. Die Länge der Abschnitte ist hier auf 2048 Sample festgelegt, was einem Zeitinterval von $$\\frac{\\mathsf{win\\_size}}{\\mathsf{sr}} = \\frac{2048}{22500} = 91.02\\ \\mathbf{ms}$$ entspricht und somit alle Kriterien ausreichend erfüllt."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "win_size = 2048\n",
    "stride = win_size >> 1\n",
    "\n",
    "window = numpy.hanning(win_size)\n",
    "\n",
    "def signal_to_frames(signal):\n",
    "    frames = librosa.util.frame(signal, frame_length=win_size, hop_length=stride)\n",
    "    return numpy.swapaxes(frames, 0, 1)\n",
    "\n",
    "def frames_to_signal(frames):\n",
    "    count = len(frames) * stride + win_size\n",
    "    ret = numpy.ndarray(count)\n",
    "    ret.fill(0)\n",
    "    offset = 0\n",
    "    for frame in frames:\n",
    "        ret[offset:offset+len(frame)] += frame * window\n",
    "        offset += stride\n",
    "    return ret[:int(sr*T)]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Die Abschitte überlappen sich zu 50%. Die nachfolgende Grafik zeigt, wie die überlappenden Abschnitte mit Hilfe der Hanning-Fensterfunktion rekombiniert werden."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 60440 characters]",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "frames = signal_to_frames(numpy.ones(stride*8))\n",
    "plt.figure(figsize=(14, 5))\n",
    "plt.xlabel(\"Zeit in s\")\n",
    "plt.ylabel(\"Amplitude\")\n",
    "for i, frame in enumerate(frames):\n",
    "    plt.plot(t[i*stride:i*stride+win_size], frame * window, '--')\n",
    "plt.plot(t[0:stride*9], frames_to_signal(frames));\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Methoden der Rauschkategorisierung\n",
    "\n",
    "Die Amplitude des Rauschens im Eingangssignal ist im betrachteten Zeitraum nicht konstant. Jedoch ist die Energie des Rauschens im gesamten Frequenzspektrum des Signals verteilt und kann damit von dem Nutzsignal unterschieden werden. Wenn der Grenzwert zwischen Störsignal und Nutzsignal bekannt ist, dann kommen unabhängig von der Domäne zwei Arten von Grenzwertfunktionen in Betracht:\n",
    "- Ein harter Filter $f1$, bei dem alle Werte unter einem Grenzwert auf Null gesetzt werden, und\n",
    "- Ein weicher Filter $f2$, bei dem der Bereich außerhalb des Grenzwertes in Richtung Null skaliert wird.\n",
    "\n",
    "Diese Filter sind folgend für Werte $arr$ mit dem Grenzwert $g$ mit $0 < g <= max(arr)$ definiert:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# Harter Filter\n",
    "def f1(arr, g):\n",
    "    return (abs(arr) >= g) * arr\n",
    "\n",
    "# Weicher Filter\n",
    "def f2(arr, g, maximum=1):\n",
    "    return numpy.sign(arr) * numpy.clip(abs(arr) - g, 0, maximum - g) * (1 / (maximum-g))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wahl des Waveletalgorithmus\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['bior1.1', 'bior1.3', 'bior1.5', 'bior2.2', 'bior2.4', 'bior2.6', 'bior2.8', 'bior3.1', 'bior3.3', 'bior3.5', 'bior3.7', 'bior3.9', 'bior4.4', 'bior5.5', 'bior6.8', 'coif1', 'coif2', 'coif3', 'coif4', 'coif5', 'coif6', 'coif7', 'coif8', 'coif9', 'coif10', 'coif11', 'coif12', 'coif13', 'coif14', 'coif15', 'coif16', 'coif17', 'db1', 'db2', 'db3', 'db4', 'db5', 'db6', 'db7', 'db8', 'db9', 'db10', 'db11', 'db12', 'db13', 'db14', 'db15', 'db16', 'db17', 'db18', 'db19', 'db20', 'db21', 'db22', 'db23', 'db24', 'db25', 'db26', 'db27', 'db28', 'db29', 'db30', 'db31', 'db32', 'db33', 'db34', 'db35', 'db36', 'db37', 'db38', 'dmey', 'haar', 'rbio1.1', 'rbio1.3', 'rbio1.5', 'rbio2.2', 'rbio2.4', 'rbio2.6', 'rbio2.8', 'rbio3.1', 'rbio3.3', 'rbio3.5', 'rbio3.7', 'rbio3.9', 'rbio4.4', 'rbio5.5', 'rbio6.8', 'sym2', 'sym3', 'sym4', 'sym5', 'sym6', 'sym7', 'sym8', 'sym9', 'sym10', 'sym11', 'sym12', 'sym13', 'sym14', 'sym15', 'sym16', 'sym17', 'sym18', 'sym19', 'sym20']\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 294060 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 47848 characters]",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 59460 characters]",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pywt\n",
    "\n",
    "print(pywt.wavelist(kind='discrete'))\n",
    "algo = pywt.Wavelet('coif2')\n",
    "r_frames = signal_to_frames(r)\n",
    "x_frames = signal_to_frames(x)\n",
    "\n",
    "def frames_to_coeffs(frames):\n",
    "    maxlevel = pywt.dwt_max_level(win_size, algo)\n",
    "    return pywt.wavedec(frames, algo, level=maxlevel)\n",
    "\n",
    "def coeffs_to_frames(coeffs):\n",
    "    return pywt.waverec(coeffs, algo)\n",
    "\n",
    "x_coeffs = frames_to_coeffs(x_frames)\n",
    "r_coeffs = frames_to_coeffs(r_frames)\n",
    "\n",
    "plt.figure(figsize=(14, 5))\n",
    "for i, coeff in enumerate(r_coeffs):\n",
    "    r_coeffs[i] = numpy.max(numpy.abs(r_coeffs[i]), axis=0)\n",
    "    plt.plot(r_coeffs[i])\n",
    "\n",
    "for i, coeff in enumerate(x_coeffs):\n",
    "    x_coeffs[i] = f2(coeff, r_coeffs[i])\n",
    "    \n",
    "y_frames = coeffs_to_frames(x_coeffs)\n",
    "\n",
    "y = frames_to_signal(y_frames)\n",
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.waveshow(x, sr=sr)\n",
    "librosa.display.waveshow(y, sr=sr)\n",
    "Audio(y, rate=sr)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<librosa.display.AdaptiveWaveplot at 0x7fb8fc380460>"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 53796 characters]",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "yt = f2(x, 0.1)\n",
    "Audio(yt, rate=sr)\n",
    "\n",
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.waveshow(x, sr=sr)\n",
    "\n",
    "librosa.display.waveshow(yt, sr=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 193648 characters]",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24856 characters]",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = librosa.stft(x, n_fft=2**11)\n",
    "Y = f2(abs(X), 0.9, numpy.max(abs(X))) * numpy.exp(1j * numpy.angle(X))\n",
    "\n",
    "y = librosa.istft(Y)\n",
    "y = y * max(x) / max(y)\n",
    "\n",
    "\n",
    "Xdb = librosa.amplitude_to_db(abs(X))\n",
    "Ydb = librosa.amplitude_to_db(abs(Y))\n",
    "\n",
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.specshow(Xdb, sr=sr, x_axis='time', y_axis='hz')\n",
    "plt.figure(figsize=(14, 5))\n",
    "\n",
    "librosa.display.specshow(Ydb, sr=sr, x_axis='time', y_axis='hz');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 52800 characters]",
      "text/plain": [
       "<Figure size 1008x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(14, 5))\n",
    "librosa.display.waveshow(x, sr=sr)\n",
    "librosa.display.waveshow(y, sr=sr);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "                <audio  controls=\"controls\" >\n",
       "                    <source src=\"data:audio/wav;base64,[encoded data, 293608 characters]\" type=\"audio/wav\" />\n",
       "                    Your browser does not support the audio element.\n",
       "                </audio>\n",
       "              "
      ],
      "text/plain": [
       "<IPython.lib.display.Audio object>"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Audio(y, rate=sr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
