{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9113be64-b859-4517-8bae-59c82152acea",
   "metadata": {},
   "source": [
    "Notebook Jupyter      UTC MT12 \n",
    "$$ \\text{TP} 4 \\;: \\; \\textbf{Transformée de Fourier Discrète}$$\n",
    "###### Cette cellule est une cellule Markdown"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fea8551e-9785-46c3-b26e-5061c2e8c1b5",
   "metadata": {},
   "source": [
    "### L'objectif principal de ce TP est de manipuler la transformée de Fourier discrète."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "607f0445-f856-4429-96a1-53c61669cbb8",
   "metadata": {},
   "source": [
    "#### Plan \n",
    "[1. Implémentation de la matrice de transformée de Fourier discrète (TFD ou DFT)](#1.-Implémentation-de-la-matrice-de-transformée-de-Fourier-discrète-(-TFD-ou-DFT-))\n",
    "\n",
    "[2. FFT et \"débruitage\" d'un signal](#2.-FFT-et-débruitage-d-'-un-signal)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ddd4c27-6d9e-4bad-881c-a1f8c33cfc57",
   "metadata": {},
   "source": [
    "# 1. Implémentation de la matrice de transformée de Fourier discrète (TFD ou DFT)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a4149e1-185e-46c2-b651-409a27f76cfd",
   "metadata": {},
   "source": [
    "## 1.1. Rappels sur la transformée de Fourier discrète"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebdc6469-3f5e-4e4d-b91d-3b8c5badbc48",
   "metadata": {},
   "source": [
    "Pour un vecteur $\\textbf{y}=(y_0,...,y_{N-1})^\\top \\in \\mathbb{R}^N$, on rappelle la définition de la transformée de Fourier discrète (~TFD ou DFT) d'ordre $N$ ($N\\geq 2$) :$$Y_n = \\frac{1}{N} \\sum_{k=0}^{N-1} y_k\\, \\omega_N^{-nk},\\qquad n=0,...,N-1,$$ où $\\omega_N$ désigne le nombre complexe suivant :$$\\omega_N = e^{2i\\frac\\pi N}.$$ \n",
    "On note\n",
    "\n",
    "\\begin{align}\n",
    "\\textbf{Y}& =(Y_0,...,Y_{N-1})^\\top, \\quad  \\text{le vecteur de} \\; \\mathbb{R}^N,\\\\\n",
    "\\mathcal{F}_N,& \\quad \\text{l'application linéaire de $\\mathbb{R}^N$ dans $\\mathbb{R}^N$ telle que}  \\; Y_n = (\\mathcal{F}_N \\textbf{y})_n,  \\text{ pour }  n \\in \\{0,\\ldots,N-1\\}.\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76b84378-0e88-4051-863e-7360665cc1b0",
   "metadata": {},
   "source": [
    "## 1.2 Questions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2c4869f-20b0-4b6a-9779-9b44ab29c66e",
   "metadata": {},
   "source": [
    "### Q1.a. Justifier sur une feuille l'écriture matricielle suivante de la TFD \\begin{align*} \\textbf{Y}  = S_N \\, \\textbf{y},\\quad \\text{où} \\;\\; & \\quad S_N = \\frac{1}{N}\\,\\, \\overline{\\Omega}_N,\\\\ \\text{avec} & \\quad  \\Omega_N : \\text{matrice carrée de  dim$=N$ de terme général} \\\\ & \\quad (\\Omega_N)_{nk} = \\text{e}^{2i\\pi \\frac{nk}{N}}  = \\omega_N^{nk}, \\;0 \\leq n,k\\leq N-1.&\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a1b7054-a710-4244-b88a-4f880ff5c3d6",
   "metadata": {},
   "source": [
    "Le produit matriciel $S_N \\, \\textbf{y}$ fournit un vecteur de dimension $N$ dont la n-ième composante est définie par:\n",
    "\\begin{eqnarray*}\n",
    "(S_N\\mathbf{y})_n& =  \\frac{1}{N} \\sum_{k=0}^{N-1} \\text{e}^{-2i\\pi \\frac{nk}{N}}y_k= \\frac{1}{N} \\sum_{k=0}^{N-1} \\omega_N^{-nk} y_k  =Y_n \n",
    "\\end{eqnarray*}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "439bc71d-860e-4509-99df-e85c021f3278",
   "metadata": {},
   "source": [
    "### Q1.b. Justifier sur une feuille que l'inverse de la matrice $S_N$ vérifie $$ (S_N)^{-1} = \\Omega_N.$$ Autrement dit, la transformée de Fourier discrète inverse (TFDI) a pour représentation matricielle la matrice $\\Omega_N$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c229ca88-c5c2-4105-a83f-36d4c4bb0500",
   "metadata": {},
   "source": [
    "### Q2 Implémentation naïve de la TFD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "db02c84a-4493-47a3-ae97-0828942e9180",
   "metadata": {},
   "source": [
    "#### Q2.a. Pour un entier strictement positif $N$,  définir une fonction $\\textbf{TFD1}$ qui implémente la matrice $S_N$ en utilisant deux boucles for."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "65407227-fe84-4f29-a235-e495d0e05f48",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cmath as cm #importe les complexes\n",
    "import scipy as sc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8f5bb4f-e1b6-44a5-a361-5ed1a422ad18",
   "metadata": {},
   "source": [
    "La fonction z.conjugate() calcule le conjugué du nombre complexe z, cette fonction s’applique aussi aux matrices."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ab0cf09-3895-4897-a9e2-6a91d157734f",
   "metadata": {},
   "source": [
    "#### Q2.b Indiquer pourquoi le code de la fonction $\\textbf{TFD1}$ n'est pas très performant. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6e9d768-b664-4c31-b354-8f179fd2fac7",
   "metadata": {},
   "source": [
    "### Q3 Implémentation plus performante de la TFD que la méthode naïve. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca22a9f9-bf0a-411b-abc9-875706792da6",
   "metadata": {},
   "source": [
    "#####  Pour cette question, on considère le vecteur $\\bf{w}\\in\\mathbb{R}^N$ dont le $n$-ième composante est définie par $$w_n= \\sqrt{\\frac{2\\pi}{N}}\\, n\\, e^{i\\pi/4},\\quad n=0,...,N-1$$  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2d486977-381c-4e77-b5d1-39f085e601fc",
   "metadata": {},
   "source": [
    "#### Q3.a. Effectuer sur une feuille  le produit  $$\\bf{w} \\bf{w}^T = \\begin{pmatrix} 0\\\\ \\vdots\\\\ \\sqrt{\\frac{2\\pi}{N}}\\, (N-1)\\, e^{i\\pi/4} \\end{pmatrix} \\begin{pmatrix} 0 & \\ldots & \\sqrt{\\frac{2\\pi}{N}}\\, (N-1)\\, e^{i\\pi/4} \\end{pmatrix} \\; $$ \n",
    "On pourra simplement calculer les éléments $(\\bf{w} \\bf{w}^T)_{nk}$ de cette matrice."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e922717e-7b07-44aa-a23a-a1318d29b38d",
   "metadata": {},
   "source": [
    "#### Q3.b. En déduire sur une feuille  que $\\exp(\\bf{w}\\bf{w}^T)=\\Omega_N $  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5982b4f-5da5-4dbc-afe8-c351ac4bec76",
   "metadata": {},
   "source": [
    "#### Q3.c. Définir une fonction $\\textbf{TFD2}$ qui prend en entrée $N\\geq 1$, qui utilise  $\\Omega_N = \\exp(\\bf{w} \\bf{w}^T)$ et qui renvoie la matrice $S_N$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7bf6a9c8-00a7-4826-8108-c9655a50c5ea",
   "metadata": {},
   "source": [
    "##### Q3.d. Comparer les temps d'exécution des fonctions $\\textbf{TFD1}$ et $\\textbf{TFD2}$ pour $N=100$, $200$, $300$ et $400$. Pour ce faire on importera la library \"time\", puis on utilisera la fonction \"time.time()\" préalablement et à l'issue du l'éxecution de chacune des fonctions $\\textbf{TFD1}$ et $\\textbf{TFD2}$,  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63b3d1cc-8502-4261-aaf4-033c55dc07be",
   "metadata": {},
   "source": [
    "##### Q3.e. Que se passe-t-il lorsqu'on lance la fonction $\\textbf{TDF2}$ pour $N=50000$ ? "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90fa62b9-5dab-4f5b-a2b9-fbb4d8ef2ae0",
   "metadata": {},
   "source": [
    "# 2. FFT et \"débruitage\" d'un signal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21f91444-f6df-4bfa-bf37-776373661f96",
   "metadata": {},
   "source": [
    "Dans la pratique pour calculer efficacement la TFD d'un signal, on utilise l'algorithme de FFT. Dans la suite on va considérer cette manière très performante de calculer la TFD. L'algorithme de FFT est  implémenté  sous le nom \\texttt{fft}."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "695562be-dbc4-404c-8793-c5285e278dd4",
   "metadata": {},
   "source": [
    "## 2.1. Questions sur signal simple"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "968183bf-720d-48d9-96ae-d0fe28a5c54b",
   "metadata": {},
   "source": [
    "#### Soit un signal simple, défini par la fonction $T$-périodique (signal) $f$ avec $T=1$,  $$f(x) = \\sin (2\\pi f_0 x) + \\sin(2 \\pi f_1 x), \\quad \\forall x \\in [0,T],$$ et  pour lequel les fréquences $f_0=50$ et $f_1=120$  sont connues."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "34f63cba-82e6-4fcc-b72d-1d61ce1fec48",
   "metadata": {},
   "source": [
    "#### Q6  Représenter graphiquement le signal $f$ en utilisant $N=1000$ points. On notera $y$ le vecteur $y=f(t)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c958a21c-27a6-4150-86be-96cfa24b0f57",
   "metadata": {},
   "source": [
    "#### Q7 TFD du vecteur $y$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97bb5295-d1e4-4a4f-b16e-8f407182323b",
   "metadata": {},
   "source": [
    "##### Q7.a. Calcul de la TFD de $y$ par la commande np.fft.fft(y). On la notera $\\text{fhat}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "411fa1ea-e561-4e2e-8674-792fa94d928d",
   "metadata": {},
   "source": [
    "##### Q7.b. La TFD du signal $f$ calculée à partir de $y$, est le vecteur $Y$ constitué des nombres complexes $Y_n$, $n=0:N-1$. \n",
    "\n",
    "Tracer le le spectre d'amplitude de $\\text{fhat}$ pour toutes les fréquences. Que constate-t-on ?   \n",
    "\n",
    "*remarque* On représente l'importance des fréquences du spectre du signal par le **spectre d'amplitude** qui correspond au module de la composante $n$ de $\\text{fhat}$ en fonction de la fréquence $n/T$ (pour tout $n=0,\\ldots,N-1$). \n",
    "\n",
    "Par ailleurs, par symétrie, voir la Section 3.2 du Chapitre 3 du cours, on représente généralement le spectre de puissance ou d'amplitude uniquement pour la première moitié des fréquences. \n",
    " "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eb7abdb2-48e2-477d-bea3-7a08c122bb14",
   "metadata": {},
   "source": [
    "## 2.2 Questions: débruitage du signal bruité "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64b60ba1-f63c-41a7-a9eb-dadfe69d2cfd",
   "metadata": {},
   "source": [
    "### Q8 Ajouter un bruit Gaussien centré d'écart-type $2.5$ au signal simple $f$. On notera $funclean$ ce signal bruité. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13be839e-8d07-4c1e-bcf4-299598a9aa22",
   "metadata": {},
   "source": [
    "### Q9  Représenter sur la même fenêtre graphique le graphe de $f$ et celui de $funclean$. \n",
    "\n",
    "Dans le suite, on s'intéresse à débruiter $funclean$  pour retrouver $f$. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a494ded-77b7-4474-b5ae-c340195a0870",
   "metadata": {},
   "source": [
    "### Q10 Calculer $funcleanhat$ la DFT de $funclean$, puis représenter sur une même fenêtre graphique les spectres d'amplitude de $fhat$ et $funcleanhat$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd728c58-bffa-49c3-9e64-a041ab5ccb0d",
   "metadata": {},
   "source": [
    "### Q11 Débruitage de $funclean$\n",
    "Afin de débruiter le signal $funclean$, on va mettre à $0$ les fréquences d'amplitude dont le module est\n",
    "inférieur à $r > 0$, un paramètre. En vous aidant du spectre d’amplitude de $funcleanhat$ proposer une\n",
    "valeur de $r$ et effectuez le débruitage.  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "80635cbd-e1a7-457a-b667-5cbb14272ad5",
   "metadata": {},
   "source": [
    "**Remarque sur le débruitage de $funclean$ proposer à la question Q11** \n",
    "\n",
    "L'opération décrite à la question précédente revient à appliquer un filtre au vecteur $funcleanhat$. Notons $H=(H_n)_n$ les composantes de ce filtre. Dans le domaine fréquentiel la définition de $H=(H_n)$ est la suivante\n",
    "\\begin{align*}\n",
    "H_n = \\left\\{\n",
    "    \\begin{array}{ll}\n",
    "        1 & \\mbox{si } |funcleanhat| \\geq r,\\\\\n",
    "        0 & \\mbox{sinon.}\n",
    "    \\end{array}\n",
    "\\right.\n",
    "\\end{align*}\n",
    "où $r$ désigne la valeur du paramètre choisie. \n",
    "\n",
    "On considère alors le vecteur $fcleanhat= (fcleanhat_n)$ avec\n",
    "\\begin{align*}\n",
    "fcleanhat = funcleanhat \\otimes H,\n",
    "\\end{align*}\n",
    "où $\\otimes$ désigne le \\href{https://fr.wikipedia.org/wiki/Produit_matriciel_de_Hadamard}{produit de Hadamard}. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "670dc7a3-a2cf-4414-9ca0-c92ec6425f7d",
   "metadata": {},
   "source": [
    "### Q12  Retour au domaine temporel. \n",
    "\n",
    "Utiliser l'inverse de la TFD \"np.fft.ifft\" de la librairie \"numpy\" pour revenir dans le domaine temporel et déterminer \n",
    "$fclean$.\n",
    "\n",
    "Puis sur une même fenêtre graphique, tracer $f$ et $fclean$. \n",
    "Le filtrage précédent a-t-il permus de débruiter $funclean$ ? Si tel n'est pas le cas, considérer une nouvelle valeur de $r$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6e188dbd-5a07-40b4-9904-bc20b5b03396",
   "metadata": {},
   "source": [
    "**Remarque** Il est important de comprendre que si la définition du filtre $H$ est naturelle (ou en tout cas intuitive) dans le domaine fréquentiel, sa définition dans le domaine temporel est beaucoup moins aisée. En particulier, noter qu'un filtre passe-bas, passe-haut ou passe bande est souvent (toujours) défini dans le domaine fréquentiel."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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
}