{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# EE350: Lab 2 ‒  Transformée de Fourier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "On vous fournis des fonctions de plot pour verifier vos résultats dans le fichier `utils.py`. \n",
    "\n",
    "**Veuillez ne pas modifier le script utils.py.**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "import urllib\n",
    "import os\n",
    "import numpy as np\n",
    "from scipy.signal import convolve2d, get_window\n",
    "from utils import *\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 Transformée de Fourier discrète (TFD)\n",
    "\n",
    "La TFD $X\\left(f\\right)$ du signal $x\\left(n \\right)$ de longueur $N$ est donnée par:\n",
    "$$\n",
    "X\\left(f\\right) = \\sum_{n = n_0}^{n_0 + N - 1} x\\left(n \\right) W_{N}^{-kn}\n",
    "$$\n",
    "\n",
    "et sa transformée inverse par:\n",
    "\n",
    "$$\n",
    "x\\left(n \\right) = \\frac{1}{N} \\sum_{k = 0}^{N - 1} X\\left(k\\right)  W_{N}^{kn}\n",
    "$$\n",
    "\n",
    "où $W_{N} = \\text{exp}\\left(j\\frac{2\\pi}{N}\\right)$.\n",
    "\n",
    "### 1.1 Base de la TFD\n",
    "\n",
    "La première équation ci-dessus montre comment un signal $x$ est projetée sur un ensemble de fonctions de base. Quelles sont ces fonctions de base?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 TFD d'un vecteur de cette base\n",
    "\n",
    "Le signal $x_w[n]$ est défini comme:\n",
    "\n",
    "$$\n",
    "x_w[n] = w[n] \\cdot x[n]\n",
    "$$ \n",
    "\n",
    "$$x[n]=\\cos \\left( \\frac{2 \\pi f n}{f_s} \\right)$$\n",
    "\n",
    "où $w[n]$ est une fenêtre rectangulaire de taille $N=100$ et , avec $f = 90$ Hz et $f_s = 1000$ Hz.\n",
    "\n",
    "Affichez la fonction $x[n]$,le signal $x_w[n]$ fenêtré par un rectangle et la TFD de $x_w[n]$ (seulement son module).\n",
    "\n",
    "- Décrivez la transformée obtenue. Combien de coefficients non-nuls la transformée de Fourier possède-t-elle? Pourquoi ?\n",
    "\n",
    "\n",
    "**Conseil**: La fenêtre peut être générée en utilisant la commande `get_window(...)` du package `scipy.signal`. Vous pouvez trouver des détails sur son utilisation dans la documentation en ligne."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_xn_1(n, xn) # n: discrete time index: np.ndarray, xn: signal: np.ndarray"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_visualize_dft(s: np.ndarray, n_pts: int, fs: float, win_type: str):\n",
    "    \"\"\" Compute DFT of signal s using window \"win_type\"\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    s : np.ndarray (n,)\n",
    "        Signal\n",
    "    n_pts : _type_\n",
    "        Size of the window to apply\n",
    "    fs : float\n",
    "        Sampling frequency\n",
    "    win_type : str\n",
    "        Type of window to use\n",
    "    \n",
    "    \"\"\"\n",
    "    #TODO: Complete this function\n",
    "    # Generate window\n",
    "    # Apply window to signal\n",
    "    # Compute DFT and its shifted version\n",
    "    \n",
    "    \n",
    "    \n",
    "    # Visualize results\n",
    "    # xw: windowed signal: np.ndarray, xwf: DFT of windowed signal: np.ndarray, xwf_shift: shifted DFT of windowed signal: np.ndarray, \n",
    "    # wn: window: np.ndarray, fs: sampling frequency: float, df: frequency bins: np.ndarray, n_pts: size of the window: int, win_type: type of window: str\n",
    "    plot_xwn_dft_1(xw, xwf, xwf_shift, wn, fs, df, n_pts, win_type) \n",
    "    return xwf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute DFT of signal and plot it\n",
    "xwf = compute_visualize_dft(s=xn, n_pts=100, fs=fs, win_type='boxcar')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.3 Longueur de la fenêtre\n",
    "\n",
    "Considérez maintenant une fenêtre rectangulaire avec une taille différente ($N = 75$). Calculez et affichez le signal et sa transformée de Fourier en utilisant la fonction `compute_visualize_dft`. \n",
    "- Quel est le nombre de coeficients non-nuls de la transformée de Fourier? Expliquez la différence avec le cas précédent."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.4 Type de fenêtre\n",
    "\n",
    "Répétez les deux étapes précédentes avec les paramètres suivants :\n",
    "\n",
    "- Fenêtre triangulaire de longueur $N = 100$\n",
    "- Fenêtre triangulaire de longueur $N = 75$\n",
    "- Fenêtre de Blackman-Harris de longueur $N = 100$\n",
    "- Fenêtre de Blackman-Harris de longueur $N = 75$\n",
    "\n",
    "Pour chaque cas, tracez le signal fenêtré $x_w(n)$ ainsi que sa TFD en utilisant la fonction `compute_visualize_dft`.\n",
    "- Quels sont le nombre de coeficients non-nuls dans chaque cas ? Comparez les résultats avec la fenêtre rectangulaire. Comparez les spectres obtenus et commentez ? Concluez sur ce qu'est un bon fenêtrage."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Transformée de Fourier discrète 2D et convolution\n",
    "\n",
    "Dans cette section, nous allons revoir la relation importante entre la transformée de Fourier et l'opération de convolution. Nous illustrerons ces concepts en 2D. \n",
    "\n",
    "Considérons un système linéaire $h$ et invariant qui prends en entrée une image en nuances de gris $y$ et produit une nouvelle image en nuance de gris $v$.\n",
    "\n",
    "L'entrée et la sortie du système sont liées par l'opération de convolution, désignée par $*$. La sortie est la convolution de l'entrée $y$ avec la réponse impulsionnelle du système $h$, c'est-à-dire $v = y * h$. En 2D, l'opérateur de convolution est défini pour chaque point de sortie $(n, m)$ comme :\n",
    "\n",
    "$\n",
    "\\begin{aligned}\n",
    "    v\\left(n, m\\right)  & = \\left(y * h\\right)\\left( n, m \\right) \\\\\n",
    "                        & = \\sum_{o=-\\infty}^{+\\infty}\\sum_{p=-\\infty}^{+\\infty} h\\left(o, p\\right) y\\left(n-o, m-p\\right)\n",
    "\\end{aligned}\n",
    "$\n",
    "\n",
    "Cette équation donne l'expression de chaque pixel $(n,m)$ de l'image de sortie par rapport à l'image d'entrée et au système. La transformée de Fourier et l'opération de convolution sont liées par une propriété importante, qui est la suivante \n",
    "\n",
    "$\n",
    "\\begin{aligned}\n",
    "V\\left(k, l\\right) & = \\mathcal{F}\\left(v\\right)\\\\\n",
    "                & = \\mathcal{F}\\left(y * h \\right) \\\\\n",
    "                & = \\mathcal{F}\\left(y\\right)\\mathcal{F}\\left(h\\right) \\\\\n",
    "                & =   Y\\left(k, l\\right)H\\left(k, l\\right)\n",
    "\\end{aligned}\n",
    "$\n",
    "\n",
    "où $k, l$ sont les fréquences dans le plan 2D, $\\mathcal{F}$ est l'opérateur de transformée de Fourier et $V\\left(k, l\\right), Y\\left(k, l\\right), H\\left(k, l\\right)$ sont les transformées de Fourier de $v\\ , y\\ , h$ respectivement.\n",
    "\n",
    "Dans le domaine de Fourier, la convolution se transforme en une simple multiplication des transformées de Fourier $Y(k, l)$ et $H(k, l)$. La fonction $H(k, l)$ est définie comme la fonction de transfert du système $h$. Dans cet exercice, nous allons montrer que les calculs de $v\\left(n, m\\right)$ avec ces deux solutions sont équivalents.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### 2.1. Domaine temporel\n",
    "\n",
    "Commençons par définir notre système $h(n,m)$. C'est un filtre gaussien 2D de variance $\\sigma^2 = 5$ et de taille $(N_1, N_2) = (21, 21)$. \n",
    "\n",
    "$$  h[n, m] = \\frac{1}{2 \\pi \\sigma^2} e^{- \\frac{1}{2 \\sigma^2} [(n-\\frac{N_1-1}{2})^2 + (m-\\frac{N_2-1}{2})^2] }  $$\n",
    "\n",
    "Ecrivez une fonction `gkernel` qui prend en entrée le nombre de points du filtre et la variance, et sort la réponse impulsionnelle 2D `h`. En utilisant la définition de l'opération de convolution, quelle est l'expression mathématique de la convolution avec notre noyau? Plus précisément, quelles sont les limites des sommes? Pourquoi ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gkernel(N: int,  var: float):\n",
    "    \"\"\" Createa a 2D Gaussian filter with given variance.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    N : int\n",
    "        Size of the filter\n",
    "    var : float\n",
    "        Variance of the filter\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    gfilter : np.ndarray (N, N)\n",
    "        Gaussian filter\n",
    "    \"\"\"\n",
    "    \n",
    "    # Create empty filter\n",
    "    gfilter = np.nonzero((N, N))\n",
    "    \n",
    "    # TODO\n",
    "    \n",
    "    return  gfilter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_gkernel(h, var) # h: Gaussian kernel: np.ndarray, var: variance of the kernel: float"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Maintenant que nous avons notre système, nous avons besoin d'une entrée. Nous allons utiliser une image libre de droits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define Image URL\n",
    "url = \"https://live.staticflickr.com/8674/16504233985_9f1060624e_q_d.jpg\"\n",
    "local = \"astronaut.jpg\"\n",
    "\n",
    "# Check if download is needed\n",
    "if not os.path.exists(local):\n",
    "    print(\"Download image ...\")\n",
    "    urllib.request.urlretrieve(url, local)\n",
    "    \n",
    "# Load and display image\n",
    "y = mpimg.imread(local)\n",
    "# Convert to gray scales\n",
    "y = 0.299 * y[..., 0] + 0.587 * y[..., 1] + 0.114 * y[..., 2]\n",
    "\n",
    "plt.imshow(y, cmap=\"gray\")\n",
    "plt.title(\"Commander Eileen Collins (NASA)\")\n",
    "plt.axis(\"off\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enfin, calculons la sortie du système $v_t$. Utilisez la fonction `convolve2d(...)` du paquet `scipy.signal` avec l'argument `mode='full'` pour calculer la convolution 2D entre l'image d'entrée $y$ et la réponse impulsionnelle du système $h$. Tracez-la en niveaux de gris avec `matplotlib.pyplot`. \n",
    "Quelle est la taille $(M1, M2)$ de l'image? Quelle est la taille $(L1, L2)$ de l'image de sortie ?\n",
    "\n",
    "A partir de la définition de l'opération de convolution :\n",
    "- Quel est l'effet principal de la transformation sur l'image ? Expliquez.\n",
    "- Expliquez la taille de l'image de sortie. Quelle est la relation mathématique entre les dimensions d'entrée et de sortie ?\n",
    "- Du point de vue de la complexité de calcul, quelle est la manière la plus efficace d'implémenter l'opération de convolution ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "plot_2_1(vt) # vt: output image: np.ndarray\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Domaine fréquentiel\n",
    "\n",
    "Maintenant, calculons la sortie $v_f$ dans le domaine de Fourier. Calculez les transformations de Fourier 2D de l'entrée $y$ et de la réponse impulsionnelle du filtre $h$ respectivement. Vous devez spécifier le nombre de points dans la transformée de Fourier finale, qui est $(L_1, L_2)$ dans notre cas. Tracez la valeur absolue de chaque transformée de Fourier en échelle logarithmique. Ensuite, calculez le signal de sortie $v_f$ et tracez le résultat à côté du précédent $v_t$.\n",
    "\n",
    "A partir de votre code: \n",
    "- Expliquez la forme de la valeur absolue des transformées de Fourier.\n",
    "- Comparez les résultats de $v_t$ et $v_f$.\n",
    "- Expliquez pourquoi nous devons spécifier un nombre de points pour les calculs de la transformée de Fourier."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_2_2(vt, vf, Y, H, V) \n",
    "# vt: output image (temporal): np.ndarray, vf: output image (freq): np.ndarray, \n",
    "# Y: TF of input image: np.ndarray, H: TF of kernel: np.ndarray, V: TF of output image: np.ndarray\n",
    "diff = np.mean(np.abs(vf) - vt)\n",
    "print(\"Difference: {:.2f}\".format(diff))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Enfin, calculons le signal $v_f$ avec les transformées de Fourier calculées sur une grille de taille $(M_1, M_2)$. Répétez les opérations du point précédent avec une grille de taille $(M_1, M_2)$ au lieu de $(L_1, L_2)$ et tracez les résultats.\n",
    "\n",
    "- Expliquez les différences entre la nouvelle et la précédente transformée de Fourier. Quelle est l'influence de la taille de la grille ?\n",
    "- Comparez les résultats de $v$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Y :np.ndarray: TF 2D de l'entrée y, H : np.ndarray: TF 2D du filtre h, V : np.ndarray: Produit fréquentiel Y*H, \n",
    "# vt : np.ndarray Résultat convolution temporelle, vf : np.ndarray: Résultat convolution fréquentielle (ifft de V)\n",
    "plot_2_2_2(Y, H, V, vt, vf) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Densités spectrales d'énergie et de puissance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.1 Signal à énergie finie\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Considérons le signal\n",
    "$$\n",
    "x(t) = e^{-|t|}, \\quad t \\in \\mathbb{R}.\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "1. Vérifier que $x(t)$ est un signal à énergie finie.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "2. En utilisant la définition, calculer la fonction d’autocorrélation\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Déterminer la densité spectrale d’énergie $\\Phi_x(f)$ en appliquant la transformée de Fourier à la fonction d'autocorrelation.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Signal à puissance moyenne finie\n",
    "\n",
    "Considérons maintenant le signal périodique\n",
    "$$\n",
    "x(t) = \\cos(2\\pi f_0 t), \\quad t \\in \\mathbb{R}.\n",
    "$$\n",
    "\n",
    "1. Montrer que ce signal est de puissance moyenne finie mais d’énergie infinie.  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "2. Calculer la fonction d’autocorrélation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "3. En déduire la densité spectrale de puissance $\\Phi_x(f)$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<center style=\"color:blue\">Votre réponse ici</center>"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
