{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a6c59f53",
   "metadata": {},
   "source": [
    "# Traitement du Signal - TP8 : Filtrage Numérique\n",
    "\n",
    "On a vu il y'a quelques semaines, juste avant les vacances de la Toussaint, certains filtres idéaux. Pour rappel, les filtres ont pour objectif de sélectionner uniquement certaines fréquences d'un signal tout en éliminant des fréquences qu'on ne souhaite plus avoir. On avait appliqué les filtres idéaux dans le domaine fréquentiel, par simple multiplication.\n",
    "\n",
    "Aujourd'hui, on va s'intéresser aux filtres réalisables. Ce sont des filtres qui eux sont appliqués en temporel, et donc applicables en temps réel (c'est-à-dire qu'on n'a pas besoin d'avoir le signal entier pour le traiter). Ils sont également réalisables car ils sont implémentables physiquement (ex: via un circuit électrique)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b29c49e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Import des librairies\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ac8f0ef",
   "metadata": {},
   "source": [
    "## Exercice 1 : On filtre comment dans le monde réel ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b41291f8",
   "metadata": {},
   "source": [
    "Commençons par les filtres à réponse impulsionnelle finie (RIF). Ces filtres s'appliquent via convolution avec le signal à traiter en temporel, et sont faciles à implémenter. De plus, ils ont l'avantage d'être stables et peuvent avoir une phase linéaire. La conception d'un filtre RIF se base souvent sur la définition d'un filtre idéal (ceux vus en TP5), mais repassés en temporel.\n",
    "\n",
    "Dans un premier temps, on va construire un premier filtre RIF, dont la réponse impulsionnelle est de la forme suivante: \n",
    "\n",
    "\\begin{equation*}\n",
    "    h[n] = \\frac{2 f_c}{F_e} sinc(2 f_c \\frac{n-n_0}{F_e})\n",
    "\\end{equation*}\n",
    "\n",
    "Ici, $h$ est de taille $N$ (qu'on choisit impair), et $n_0$ le centre de la réponse impulsionnelle ($n_0 = \\frac{N-1}{2}$). $F_e$ est la fréquence d'échantillonnage et $f_c$ est la fréquence de coupure du filtre.\n",
    "\n",
    "Développez une fonction qui créera la réponse impulsionnelle $h$ en fonction de la taille $N$, la fréquence de coupure $f_c$ et la fréquence d'échantillonnage $F_e$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a5c3296",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Développement du filtre RIF passe-bas\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1376f1fa",
   "metadata": {},
   "source": [
    "Construisez le filtre de réponse impulsionnelle $h$ avec $N = 1001$, $f_c = 50$ Hz et $F_e = 1000$ Hz. Tracez le résultat obtenu."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d08bcb3d",
   "metadata": {},
   "source": [
    "**_QUESTION :_** Durant combien de temps le filtre évolue-t-il ? \n",
    "\n",
    "*Ca veut dire qu'il faut que vous calculiez l'axe temporel du filtre si vous n'aviez pas compris...*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "416e1c4c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Construction et affichage de h\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3465e19b",
   "metadata": {},
   "source": [
    "**_REPONSE :_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c1e9b99",
   "metadata": {},
   "source": [
    "Calculez la transformée de Fourier de la réponse impulsionnelle et affichez le spectre d'amplitude.\n",
    "\n",
    "*Note : Ici, il n'est pas nécessaire de normaliser le spectre d'amplitude, étant donné qu'on calcule la FFT d'un filtre (déjà normalisé) et non d'un signal*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3776ca9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# FFT de la réponse impulsionnelle et affichage du spectre d'amplitude\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4918d299",
   "metadata": {},
   "source": [
    "**_QUESTIONS :_**\n",
    "- Quelle est la nature du filtre développé ?\n",
    "- Quelle est la fréquence de coupure du filtre (en fonction du spectre d'amplitude) ?\n",
    "- Quelle(s) différence(s) remarquez-vous par rapport à un filtre idéal ?\n",
    "\n",
    "*Note : Pour trouver la fréquence de coupure, utilisez la fonction find_peaks de Scipy*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25e0c237",
   "metadata": {},
   "source": [
    "**_REPONSES :_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d6ffbc5",
   "metadata": {},
   "source": [
    "Afin de corriger les imperfections de ce filtre, on **fenêtre toujours la réponse impulsionnelle avec une fenêtre de Hann** (ou Hamming, ou Blackman). Corriger la réponse impulsionnelle et tracez le spectre d'amplitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5be2f344",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Correction de la réponse impulsionnelle via une fenêtre de Hann, calcul de la FFT et affichage du spectre d'amplitude\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80532ad8",
   "metadata": {},
   "source": [
    "**_QUESTION :_** Est-ce que le fenêtrage de Hann a amélioré le spectre d'amplitude du filtre ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c35d9065",
   "metadata": {},
   "source": [
    "**_REPONSE :_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e051e359",
   "metadata": {},
   "source": [
    "Prenons maintenant un signal $x(t)$, défini par l'équation suivante:\n",
    "\n",
    "\\begin{equation*}\n",
    "    x(t) = 3 \\cdot sin(2 \\pi f_0 t - \\frac{\\pi}{2}) + sin(2 \\pi f_1 t + \\frac{\\pi}{3}) + 0.5 \\cdot sin(2 \\pi f_2 t) + 0.1 \\cdot sin(2 \\pi f_3 t + \\frac{\\pi}{6}) \n",
    "\\end{equation*}\n",
    "\n",
    "Avec $f_0 = 5 Hz$, $f_1 = 40 Hz$, $f_2 = 100 Hz$ et $f_3 = 200 Hz$.\n",
    "\n",
    "Construisez ce signal, évoluant durant 4 secondes avec $F_e$ = 1000 Hz. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37fda796",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Construction de x(t) et affichage du signal\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "789d008e",
   "metadata": {},
   "source": [
    "Construisez maintenant $y(t)$, qui correspond **théoriquement** au filtrage de $x(t)$ par un filtre passe-bas à fréquence de coupure de 50 Hz (c'est-à-dire $x(t)$ sans les sinus coupés par le filtre), et affichez le signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f2f45b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Construction de y(t) et affichage du signal\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f1f2303",
   "metadata": {},
   "source": [
    "Filtrez $x(t)$ par convolution avec la réponse impulsionnelle $h(t)$ (via la fonction *np.convolve*), et gardez uniquement le signal sur les 4 premières secondes (taille de $x(t)$). Affichez le résultat du filtrage ainsi que le signal $y(t)$ sur un même graphique."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54d482da",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Filtrage de x(t) via convolution avec h(t) et affichage sur graphique avec y(t)\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4d67903",
   "metadata": {},
   "source": [
    "**_QUESTION :_** Que remarquez-vous ? Le résultat du filtrage est-il correct ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b68aa5ca",
   "metadata": {},
   "source": [
    "**_REPONSE :_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "161755d4",
   "metadata": {},
   "source": [
    "**_QUESTION :_** A quoi correspond le retard observé sur le signal ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8461a2a",
   "metadata": {},
   "source": [
    "**_REPONSE :_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f126090a",
   "metadata": {},
   "source": [
    "Bravo, vous avez conçu un filtre FIR passe-bas. Avant de passer à la suite, à vous de développer les fonctions pour concevoir :\n",
    "- un filtre passe-bande\n",
    "- un filtre passe-haut\n",
    "- un filtre coupe-bande\n",
    "\n",
    "*Pour vous aider : grâce à la propriété de linéarité de la transformée de Fourier, vous pouvez concevoir les 3 filtres demandés via une combinaison linéaire de filtres passe-bas. Je vous conseille de dessiner les portes des filtres passe-bas en fréquentiel, afin de trouver les bonnes opérations à réaliser.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8c77071",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Fonctions de conception de filtres passe-bande, passe-haut et coupe-bande\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d996900d",
   "metadata": {},
   "source": [
    "Créez une version de chacun de ces filtres et visualisez le spectre d'amplitude après transformée de Fourier pour vérifier que votre code est correct."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e8fa1dc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Création, transformée et affichage du spectre d'amplitude d'un filtre passe-bande/passe-haut/coupe-bande\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4c789dc",
   "metadata": {},
   "source": [
    "On va maintenant étudier l'impact de la taille $N$ de la réponse impulsionnelle du filtre sur les performances du filtrage. Construisez un même filtre passe-bas en modifiant la taille $N$ de la réponse impulsionnelle de $h(t)$ ($N = 51, 1001, 5001$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80744ced",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Construction de h(t) en variant N et affichage\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0a442a9",
   "metadata": {},
   "source": [
    "Filtrez le signal $x(t)$ avec chacunes de vos réponses impulsionnelles et affichez le résultat sur 3 figures/sous-figures différentes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd55febd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Filtrage de x(t) avec h(t) de différentes tailles N et affichage\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "553581d2",
   "metadata": {},
   "source": [
    "**_QUESTION :_** En quoi la taille $N$ de la réponse impulsionnelle $h(t)$ influence le résultat du filtrage de $x(t)$ ?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe35a4d7",
   "metadata": {},
   "source": [
    "**_REPONSE :_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c78bb64c",
   "metadata": {},
   "source": [
    "On va maintenant tracer le diagramme de Bode de chacun des filtres construits. Vous avez ci-dessous une fonction qui trace automatiquement le diagramme en question :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f0d975d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import freqz\n",
    "\n",
    "def plot_bode(h, Fe, titre_figure=None):\n",
    "\n",
    "    \"\"\"\n",
    "    Trace le diagramme de Bode d'un filtre de réponse impulsionnelle h\n",
    "    h : tableau NumPy 1D, la réponse impulsionnelle à tracer\n",
    "    Fe : float, la fréquence d'échantillonnage\n",
    "    titre_figure: str (optionnel), le titre de la figure \n",
    "    \"\"\"\n",
    "\n",
    "    # Calcul de H(e^{jw})\n",
    "    w, H = freqz(h, worN=2048)\n",
    "\n",
    "    # Axe fréquentiel en Hz\n",
    "    f = w * Fe / (2 * np.pi)\n",
    "\n",
    "    # Tracé du diagramme de Bode\n",
    "    plt.figure(figsize=(10, 6))\n",
    "\n",
    "    # Module en dB\n",
    "    plt.subplot(2, 1, 1)\n",
    "    plt.plot(f, 20 * np.log10(np.abs(H) + 1e-12))\n",
    "    plt.xticks(np.arange(0, Fe/2 + 1, 50))\n",
    "    plt.ylabel(\"Gain (dB)\")\n",
    "    if titre_figure is not None:\n",
    "        plt.title(titre_figure)\n",
    "    plt.grid(True)\n",
    "\n",
    "    # Phase en radians\n",
    "    plt.subplot(2, 1, 2)\n",
    "    plt.plot(f, np.unwrap(np.angle(H)))\n",
    "    plt.xticks(np.arange(0, Fe/2 + 1, 50))\n",
    "    plt.xlabel(\"Fréquence (Hz)\")\n",
    "    plt.ylabel(\"Phase (rad)\")\n",
    "    plt.grid(True)\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "733cb468",
   "metadata": {},
   "source": [
    "Tracez le diagramme de Bode pour chacun des filtres (avec la taille $N$ qui varie)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d76a2c5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Affichage des diagrammes de Bode des filtres en fonction de leur taille N\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6205bf21",
   "metadata": {},
   "source": [
    "**_QUESTION :_**  En quoi la taille $N$ de la réponse impulsionnelle $h(t)$ influence le diagramme de Bode des filtres ?\n",
    "\n",
    "*Note : concentrez-vous uniquement sur le diagramme du gain (dB), la phase n'est pas importante ici...*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20a14631",
   "metadata": {},
   "source": [
    "**_REPONSE :_**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1109c65",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Exercice 2 : Filtrage en temps réel ! En direct live !\n",
    "\n",
    "Un des avantages d'un filtre RIF, par rapport aux filtres idéaux, est qu'il peut être appliqué en temps réel. Cela veut dire qu'on n'a pas besoin d'avoir le signal en entier pour le filtrer. \n",
    "\n",
    "Pour rappel, un signal $x[n]$, filtré par convolution avec une réponse impulsionnelle $h[n]$ de taille $N$, donne le signal $y[n]$ suivant :\n",
    "\n",
    "\\begin{equation*}\n",
    "    y[n] = \\sum_{k=0}^{N} h[k] x[n-k]\n",
    "\\end{equation*}\n",
    "\n",
    "On peut donc voir ici que pour calculer $y$ à l'instant $n$, il nous faut la réponse impulsionnelle en entier ainsi que les $N$ dernières valeurs du signal d'entrée $x[n]$. De ce fait, au cours du temps, il n'est pas nécessaire d'avoir en mémoire l'intégralité du signal d'entrée $x$ mais uniquement les $N$ dernières valeurs.\n",
    "\n",
    "Le code ci-dessous reproduit un signal acquis et tracé en temps réel (ralenti) à partir d'une fonction $x(t)$ définie :\n",
    "\n",
    "\\begin{equation*}\n",
    "    x(t) = sin(2 \\pi \\cdot 50 \\cdot t) + 0.5 \\times sin (2 \\pi \\cdot 200 \\cdot t)\n",
    "\\end{equation*}\n",
    "\n",
    "L'acquisition et l'affichage se fait ici en continu, donc vous pouvez interrompre le code à tout moment (CTRL+C ou bouton d'arrêt d'exécution de la cellule)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e7ef539",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulation de l'acquisition et du tracé en temps réel de x(t)\n",
    "def acquisition_signal(t):\n",
    "    \"\"\"\n",
    "    Fonction qui simule l'acquisition du signal x en fonction du temps t\n",
    "    t: float, l'indice temporel\n",
    "    \"\"\"\n",
    "    return np.sin(2*np.pi*50*t) + 0.5*np.sin(2*np.pi*200*t)\n",
    "\n",
    "\n",
    "from IPython.display import clear_output, display\n",
    "\n",
    "Fe = 800\n",
    "T = 0.5           # durée de la fenêtre affichée (s)\n",
    "Nwin = int(Fe*T)  # nombre de points dans la fenêtre\n",
    "\n",
    "# Axe temporel de la fenêtre d'affichage\n",
    "t_win = np.arange(0, Nwin) / Fe\n",
    "\n",
    "# Signal x dans la fenêtre d'affichage\n",
    "x_win = np.zeros(Nwin)\n",
    "\n",
    "# Préparation figure\n",
    "fig, ax1 = plt.subplots(1, 1, figsize=(10, 4), sharex=True)\n",
    "\n",
    "ax1.set_title(\"Acquisition temps réel : x(t)\")\n",
    "ax1.grid(True)\n",
    "\n",
    "line_x, = ax1.plot([], [], label=\"x(t)\")\n",
    "ax1.legend()\n",
    "ax1.set_xlabel(\"Temps (s)\")\n",
    "ax1.set_xlim(0, T)\n",
    "ax1.set_ylim(-2, 2)\n",
    "\n",
    "try:\n",
    "    n = 0\n",
    "    while True:\n",
    "        # Temps absolu courant\n",
    "        t_current = n / Fe\n",
    "\n",
    "        # Acquisition simulée\n",
    "        x_n = acquisition_signal(t_current)\n",
    "\n",
    "        # Mise à jour de la fenêtre\n",
    "        if n < Nwin:\n",
    "            # Phase de \"remplissage\" : on ajoute à la suite\n",
    "            t_win[n] = t_current\n",
    "            x_win[n] = x_n\n",
    "            # Mise à jour des données affichées\n",
    "            line_x.set_data(t_win[:n+1], x_win[:n+1])\n",
    "\n",
    "        else:\n",
    "            # Fenêtre pleine : décalage à gauche + ajout du nouveau point à la fin\n",
    "            t_win[:-1] = t_win[1:]\n",
    "            x_win[:-1] = x_win[1:]\n",
    "            t_win[-1] = t_current\n",
    "            x_win[-1] = x_n\n",
    "            # Mise à jour des données affichées\n",
    "            line_x.set_data(t_win, x_win)\n",
    "\n",
    "        # Limites en x : fenêtre glissante de largeur T\n",
    "        if t_current < T:\n",
    "            ax1.set_xlim(0, T)\n",
    "        else:\n",
    "            ax1.set_xlim(t_current - T, t_current)\n",
    "\n",
    "        # Rafraîchissement dans le notebook\n",
    "        clear_output(wait=True)\n",
    "        display(fig)\n",
    "\n",
    "        n += 1\n",
    "\n",
    "except KeyboardInterrupt:\n",
    "    print(\"\\nAcquisition interrompue par l'utilisateur.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6bdb2279",
   "metadata": {},
   "source": [
    "Dans ce code, l'acquisition est simulée via la fonction *acquisition_x*. On peut mettre à la place un signal réel lu via un capteur, mais ici, on se contentera de simuler la fonction *x(t)*."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a176174c",
   "metadata": {},
   "source": [
    "On va filtrer le signal $x(t)$ par un filtre RIF passe-bas comme vu dans l'exercice 1. Construisez $h[n]$ de taille $N = 201$ avec une fréquence de coupure $f_c = 100$ Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d3877fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Construction et affichage de la réponse impulsionnelle h correspondant à un filtre RIF passe-bas\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "400fe8c2",
   "metadata": {},
   "source": [
    "L'objectif va être donc de filtrer $x[n]$ en temps réel. Le code ci-dessous est adapté pour afficher $x[n]$, le signal d'origine, et $y[n]$, le signal filtré, en temps réel. Pour l'instant, $y[n]$ n'étant pas calculé (ça sera à vous de le faire), $y[n]$ = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dddcf6be",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulation de l'acquisition de x(t), du filtrage en (t) et du tracé en temps réel des deux signaux\n",
    "Fe = 800\n",
    "T = 0.5\n",
    "Nwin = int(Fe*T)\n",
    "\n",
    "# Axe temporel de la fenêtre d'affichage\n",
    "t_win = np.arange(0, Nwin) / Fe\n",
    "\n",
    "# Signaux x, y dans la fenêtre d'affichage\n",
    "x_win = np.zeros(Nwin)\n",
    "y_win = np.zeros(Nwin)\n",
    "\n",
    "# Préparation figure\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6), sharex=True)\n",
    "\n",
    "ax1.set_title(\"Acquisition et filtrage FIR temps réel\")\n",
    "ax1.set_ylabel(\"x(t)\")\n",
    "ax1.grid(True)\n",
    "\n",
    "ax2.set_ylabel(\"y(t)\")\n",
    "ax2.set_xlabel(\"Temps (s)\")\n",
    "ax2.grid(True)\n",
    "\n",
    "line_x, = ax1.plot([], [], label=\"x(t)\")\n",
    "line_y, = ax2.plot([], [], label=\"y(t)\", color=\"orange\")\n",
    "\n",
    "ax1.legend()\n",
    "ax2.legend()\n",
    "\n",
    "ax1.set_xlim(0, T)\n",
    "ax2.set_xlim(0, T)\n",
    "\n",
    "ax1.set_ylim(-2, 2)\n",
    "ax2.set_ylim(-2, 2)\n",
    "\n",
    "try:\n",
    "    n = 0\n",
    "    while True:\n",
    "        # Temps absolu courant\n",
    "        t_current = n / Fe\n",
    "\n",
    "        # Acquisition simulée\n",
    "        x_n = acquisition_signal(t_current)\n",
    "        \n",
    "        # \"Filtrage\" de x(t) stocké dans y(t) \n",
    "        y_n = 0\n",
    "\n",
    "        # Mise à jour de la fenêtre\n",
    "        if n < Nwin:\n",
    "            # Phase de \"remplissage\" : on ajoute à la suite\n",
    "            t_win[n] = t_current\n",
    "            x_win[n] = x_n\n",
    "            y_win[n] = y_n\n",
    "            \n",
    "            # Mise à jour des données affichées\n",
    "            line_x.set_data(t_win[:n+1], x_win[:n+1])\n",
    "            line_y.set_data(t_win[:n+1], y_win[:n+1])\n",
    "\n",
    "        else:\n",
    "            # Fenêtre pleine : décalage à gauche + ajout du nouveau point à la fin\n",
    "            t_win[:-1] = t_win[1:]\n",
    "            x_win[:-1] = x_win[1:]\n",
    "            y_win[:-1] = y_win[1:]\n",
    "            t_win[-1] = t_current\n",
    "            x_win[-1] = x_n\n",
    "            y_win[-1] = y_n\n",
    "            # Mise à jour des données affichées\n",
    "            line_x.set_data(t_win, x_win)\n",
    "            line_y.set_data(t_win, y_win)\n",
    "\n",
    "        # Limites en x : fenêtre glissante de largeur T\n",
    "        if t_current < T:\n",
    "            ax1.set_xlim(0, T)\n",
    "            ax2.set_xlim(0, T)\n",
    "        else:\n",
    "            ax1.set_xlim(t_current - T, t_current)\n",
    "            ax2.set_xlim(t_current - T, t_current)\n",
    "\n",
    "        # Rafraîchissement dans le notebook\n",
    "        clear_output(wait=True)\n",
    "        display(fig)\n",
    "\n",
    "        n += 1\n",
    "\n",
    "except KeyboardInterrupt:\n",
    "    print(\"\\nAcquisition interrompue par l'utilisateur.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56b1d6ec",
   "metadata": {},
   "source": [
    "A vous maintenant de compléter le code ci-dessous, qui permettra ainsi de réaliser le filtrage *en temps-réel* de $x[n]$ pour obtenir $y[n]$. \n",
    "\n",
    "Pour vous guider, il y'a deux étapes à compléter :\n",
    "- *\"A COMPLETER 1\"* : Il faut mettre à jour *buffer_x*, un vecteur de taille $M$ contenant les $M$ dernières valeurs de $x$. Le remplissage se fait dans l'ordre inverse (c'est-à-dire que la dernière valeur du vecteur correspond à la valeur la plus récente de $x$)\n",
    "- *\"A COMPLETER 2\"* : Il faut calculer $y[n]$, en fonction des $M$ dernières valeurs de $x$ et de la réponse impulsionnelle $h$, en suivant l'équation en début d'exercice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92bb26ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Filtrage en temps réel de x[n] via la réponse impulsionnelle h du filtre FIR passe-bas\n",
    "Fe = 800\n",
    "T = 0.5\n",
    "Nwin = int(Fe*T)  # nombre de points dans la fenêtre\n",
    "\n",
    "# Axe temporel de la fenêtre d'affichage\n",
    "t_win = np.arange(0, Nwin) / Fe\n",
    "\n",
    "# Signaux x, y dans la fenêtre d'affichage\n",
    "x_win = np.zeros(Nwin)\n",
    "y_win = np.zeros(Nwin)\n",
    "\n",
    "buffer_x = np.zeros(N) # Vecteur qui stocke les N dernières valeurs de x\n",
    "\n",
    "# Préparation figure\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6), sharex=True)\n",
    "\n",
    "ax1.set_title(\"Acquisition et filtrage FIR temps réel\")\n",
    "ax1.set_ylabel(\"x(t)\")\n",
    "ax1.grid(True)\n",
    "\n",
    "ax2.set_ylabel(\"y(t)\")\n",
    "ax2.set_xlabel(\"Temps (s)\")\n",
    "ax2.grid(True)\n",
    "\n",
    "line_x, = ax1.plot([], [], label=\"x(t)\")\n",
    "line_y, = ax2.plot([], [], label=\"y(t)\", color=\"orange\")\n",
    "\n",
    "ax1.legend()\n",
    "ax2.legend()\n",
    "\n",
    "ax1.set_xlim(0, T)\n",
    "ax2.set_xlim(0, T)\n",
    "\n",
    "ax1.set_ylim(-2, 2)\n",
    "ax2.set_ylim(-2, 2)\n",
    "\n",
    "try:\n",
    "    n = 0\n",
    "    while True:\n",
    "        # Temps absolu courant\n",
    "        t_current = n / Fe\n",
    "\n",
    "        # Acquisition simulée\n",
    "        x_n = acquisition_signal(t_current)\n",
    "        \n",
    "        # A COMPLETER 1\n",
    "        # Mise à jour du buffer x (qui stocke les M dernières valeurs de x)\n",
    "        ...\n",
    "\n",
    "        # A COMPLETER 2\n",
    "        # Calcul de y_n, la valeur de y au temps n, correspondant au filtrage de x par h à l'instant n\n",
    "        y_n = ...\n",
    "\n",
    "        # Mise à jour de la fenêtre\n",
    "        if n < Nwin:\n",
    "            # Phase de \"remplissage\" : on ajoute à la suite\n",
    "            t_win[n] = t_current\n",
    "            x_win[n] = x_n\n",
    "            y_win[n] = y_n\n",
    "            \n",
    "            # Mise à jour des données affichées\n",
    "            line_x.set_data(t_win[:n+1], x_win[:n+1])\n",
    "            line_y.set_data(t_win[:n+1], y_win[:n+1])\n",
    "\n",
    "        else:\n",
    "            # Fenêtre pleine : décalage à gauche + ajout du nouveau point à la fin\n",
    "            t_win[:-1] = t_win[1:]\n",
    "            x_win[:-1] = x_win[1:]\n",
    "            y_win[:-1] = y_win[1:]\n",
    "            t_win[-1] = t_current\n",
    "            x_win[-1] = x_n\n",
    "            y_win[-1] = y_n\n",
    "            # Mise à jour des données affichées\n",
    "            line_x.set_data(t_win, x_win)\n",
    "            line_y.set_data(t_win, y_win)\n",
    "\n",
    "        # Limites en x : fenêtre glissante de largeur T\n",
    "        if t_current < T:\n",
    "            ax1.set_xlim(0, T)\n",
    "            ax2.set_xlim(0, T)\n",
    "        else:\n",
    "            ax1.set_xlim(t_current - T, t_current)\n",
    "            ax2.set_xlim(t_current - T, t_current)\n",
    "\n",
    "        # Rafraîchissement dans le notebook\n",
    "        clear_output(wait=True)\n",
    "        display(fig)\n",
    "\n",
    "        n += 1\n",
    "\n",
    "except KeyboardInterrupt:\n",
    "    print(\"\\nAcquisition interrompue par l'utilisateur.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef5af1a7",
   "metadata": {},
   "source": [
    "Super ! On a notre filtre RIF fonctionnel en temps réel ! Maintenant, on peut faire mieux. Car les filtres RIF ont deux inconvénients non négligeables :\n",
    "- La taille de la réponse impulsionnelle est conséquente (il faut au minimum une centaine de coefficients, voir plus si on veut qu'il soit plus efficace !). De ce fait, le calcul est long (si on veut faire de l'embarqué avec de l'électronique notamment)\n",
    "- Dû à la taille de $h$, le retard est aussi important (ici, plus d'un dixième de seconde). Dans un contexte de simulation comme ici, c'est ok, mais dans des applications plus critiques (ex: traitement de capteur pour un véhicule autonome), le retard n'est clairement pas négligeable\n",
    "\n",
    "Pour palier à ces deux problèmes, on va pouvoir aborder une nouvelle catégorie de filtres : les filtres à impulsion infinie (RII)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1d32d50",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "## Exercice 3 : RII-ra bien qui RII-ra le dernier\n",
    "\n",
    "Comme vu en CM, les filtres RII sont des filtres dits récursifs, c'est-à-dire que le filtrage dépend du signal original mais également du signal précédemment filtré. Un signal $x[n]$, filtré par un filtre RII, donne le signal $y[n]$ suivant :\n",
    "\n",
    "\\begin{equation*}\n",
    "    y[n] = \\sum_{k=0}^{M} b_k x[n-k] - \\sum_{p=1}^{N} a_p y[n-p]\n",
    "\\end{equation*}\n",
    "\n",
    "Le filtre RII ici n'est pas définit par une réponse impulsionnelle $h(t)$ comme pour un filtre RIF, mais par des coefficients :\n",
    "- $M$ coefficients $b_k$, liés au signal d'entrée $x(t)$\n",
    "- $N$ coefficients $a_p$, liés au signal filtré $y(t)$ aux instants précédents.\n",
    "\n",
    "Contrairement au filtre RIF, il va falloir non seulement stocker les $M$ précédentes valeurs de $x(t)$ mais également les $N$ précédentes valeurs de $y(t)$ pour un filtre RIF, soit $N+M$ valeurs stockées. Cependant, au lieu de stocker une centaine de valeurs, voire plus, on peut avoir ici des résultats très satisfaisants avec moins de dix coefficients au total."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e1e48fdc",
   "metadata": {},
   "source": [
    "Avant d'expliquer plus en détail, on va construire d'abord un filtre RII. Pour cela, créez un filtre de Butterworth passe-bas d'ordre 4 via la fonction *scipy.signal.butter* (Allez voir la doc !). Affichez ensuite les coefficients *a* et *b* obtenus."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6528bdd4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Construction du filtre de Butterworth\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa6feae2",
   "metadata": {},
   "source": [
    "Vous pouvez voir ici qu'avec un filtre de Butterworth d'ordre 4, vous avez $N=5$ coefficients $a_p$ ($a_0$, $a_1$, $a_2$, $a_3$ et $a_4$) et $M=5$ coefficients $b_k$ ($b_0$, $b_1$, $b_2$, $b_3$ et $b_4$). Cela veut dire qu'il faudra uniquement stocker les 5 dernières valeurs de $x$ (donc $x[n]$, $x[n-1]$, $x[n-2]$, $x[n-3]$ et $x[n-4]$), et les 5 dernières valeurs de $y$ pour pouvoir calculer $y[n]$ .\n",
    "\n",
    "Notez que le coefficient $a_0$ = 1. Ce coefficient est en réalité à ne pas prendre en compte, car il correspond au coefficient de $y$ à l'instant $n$, donc finalement ce qu'on cherche à calculer."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc1e113a",
   "metadata": {},
   "source": [
    "On va donc reprendre le code final du précédent exercice et l'adapter pour fonctionner avec un filtre RII. A vous de compléter le code aux endroits suivants :\n",
    "\n",
    "- *\"A COMPLETER 1\"* : Il faut initialiser *buffer_x* et *buffer_y*, les vecteurs stockant respectivement les $M$ dernières valeurs de $x$ et les $N$ dernières valeurs de $y$\n",
    "- *\"A COMPLETER 2\"* : Il faut mettre à jour *buffer_x* en décalant toutes les valeurs d'un pas à droite, puis en ajoutant en début de vecteur $x[n]$\n",
    "- *\"A COMPLETER 3\"* : Il faut calculer $y[n]$ en suivant l'équation aux différences décrite en début d'exercice\n",
    "- *\"A COMPLETER 4\"* : Il faut mettre à jour *buffer_y* en décalant toutes les valeurs d'un pas à droite, puis en ajoutant en début de vecteur $y[n]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e814d7a2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Filtrage RII en temps réel de x(t)\n",
    "\n",
    "# Paramètres généraux\n",
    "Fe = 800\n",
    "T = 0.5\n",
    "Nwin = int(Fe*T)  # nombre de points dans la fenêtre\n",
    "\n",
    "# Axe temporel de la fenêtre d'affichage\n",
    "t_win = np.arange(0, Nwin) / Fe\n",
    "\n",
    "# Signaux x, y dans la fenêtre d'affichage\n",
    "x_win = np.zeros(Nwin)\n",
    "y_win = np.zeros(Nwin)\n",
    "\n",
    "# A COMPLETER 1\n",
    "# Création des buffers pour x et y\n",
    "...\n",
    "\n",
    "# -----------------------------\n",
    "# Préparation figure\n",
    "# -----------------------------\n",
    "fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 6), sharex=True)\n",
    "\n",
    "ax1.set_title(\"Acquisition et filtrage Butterworth (IIR) temps réel\")\n",
    "ax1.set_ylabel(\"x(t)\")\n",
    "ax1.grid(True)\n",
    "\n",
    "ax2.set_ylabel(\"y(t)\")\n",
    "ax2.set_xlabel(\"Temps (s)\")\n",
    "ax2.grid(True)\n",
    "\n",
    "line_x, = ax1.plot([], [], label=\"x(t)\")\n",
    "line_y, = ax2.plot([], [], label=\"y(t)\", color=\"orange\")\n",
    "\n",
    "ax1.legend()\n",
    "ax2.legend()\n",
    "\n",
    "ax1.set_xlim(0, T)\n",
    "ax2.set_xlim(0, T)\n",
    "\n",
    "ax1.set_ylim(-2, 2)\n",
    "ax2.set_ylim(-2, 2)\n",
    "\n",
    "try:\n",
    "    n = 0\n",
    "    while True:\n",
    "        # Temps absolu courant\n",
    "        t_current = n / Fe\n",
    "\n",
    "        # Acquisition simulée\n",
    "        x_n = acquisition_signal(t_current)\n",
    "\n",
    "        # A COMPLETER 2\n",
    "        # Mise à jour du buffer de x (via un décalage à droite des valeurs pour insérer x[n] en début de buffer)\n",
    "        ...\n",
    "\n",
    "        # A COMPLETER 3\n",
    "        # Calcul de y[n] en suivant l'équation aux différences\n",
    "        y_n = ...\n",
    "\n",
    "\n",
    "        # A COMPLETER 4\n",
    "        # Mise à jour du buffer de y (via un décalage à droite des valeurs pour insérer y[n] en début de buffer)\n",
    "        ...\n",
    "\n",
    "        # Mise à jour de la fenêtre\n",
    "        if n < Nwin:\n",
    "            # Phase de remplissage\n",
    "            t_win[n] = t_current\n",
    "            x_win[n] = x_n\n",
    "            y_win[n] = y_n\n",
    "\n",
    "            line_x.set_data(t_win[:n+1], x_win[:n+1])\n",
    "            line_y.set_data(t_win[:n+1], y_win[:n+1])\n",
    "        else:\n",
    "            # Fenêtre pleine : on décale à gauche et on ajoute en fin\n",
    "            t_win[:-1] = t_win[1:]\n",
    "            x_win[:-1] = x_win[1:]\n",
    "            y_win[:-1] = y_win[1:]\n",
    "            t_win[-1] = t_current\n",
    "            x_win[-1] = x_n\n",
    "            y_win[-1] = y_n\n",
    "\n",
    "            line_x.set_data(t_win, x_win)\n",
    "            line_y.set_data(t_win, y_win)\n",
    "\n",
    "        # Limites en x : fenêtre glissante de largeur T\n",
    "        if t_current < T:\n",
    "            ax1.set_xlim(0, T)\n",
    "            ax2.set_xlim(0, T)\n",
    "        else:\n",
    "            ax1.set_xlim(t_current - T, t_current)\n",
    "            ax2.set_xlim(t_current - T, t_current)\n",
    "\n",
    "        # Rafraîchissement dans le notebook\n",
    "        clear_output(wait=True)\n",
    "        display(fig)\n",
    "\n",
    "        n += 1\n",
    "\n",
    "except KeyboardInterrupt:\n",
    "    print(\"\\nAcquisition interrompue par l'utilisateur.\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18e3a33e",
   "metadata": {},
   "source": [
    "Réalisez maintenant un filtrage passe-haut de Butterworth (ordre 4, $f_c$ = 100 Hz), et appliquez-le en temps réel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f49f0007",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Filtrage de Butterworth passe-haut de x(t) et affichage en temps réel\n",
    "..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b4b8f80",
   "metadata": {},
   "source": [
    "Pour finir, afin de mieux explorer la documentation de Scipy, effectuez un filtre passe-bas de Tchebychev, puis un filtre passe-bas de Cauer (ou elliptique), et appliquez le filtrage en temps réel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5404030d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Filtrage de Tchebychev passe-bas de x(t) et affichage en temps réel\n",
    "..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77cb466c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# A COMPLETER\n",
    "# Filtrage de Cauer passe-bas de x(t) et affichage en temps réel\n",
    "..."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "deep",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}