{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "changing-skating",
   "metadata": {},
   "source": [
    "# TD1 exo 1 sur le TFD (Corrige Octave)\n",
    "\n",
    "Exécution interactive en ligne ici [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/balaise31/Signal/master?labpath=discret%2Ftd%2FFREQ_code%2Fexo1_tfd_corr_octave.ipynb)\n",
    "\n",
    "|        | Sujet | Corrigé |\n",
    "|--------|-------|---------|\n",
    "|Python  | [sujet python](./exo1_tfd.ipynb) | [corrigé python](./exo1_tfd_corr.ipynb) |\n",
    "|Octave  | [sujet octave](./exo1_tfd_octave.ipynb) | [corrigé octave](./exo1_tfd_corr_octave.ipynb) |\n",
    "\n",
    "\n",
    "\n",
    "Retour au [sujet de TD papier](../FREQ_sujet.ipynb)\n",
    "\n",
    "\n",
    "## TFD avec l'algorithme fft\n",
    "---\n",
    "Là tout est simple car tout est discret et N périodique :\n",
    "\n",
    "- $t= k T_e$ avec N points\n",
    "- $f=k\\Delta_f$ avec N points dans $[0, F_e [$ ouvert à droite !!!\n",
    "\n",
    "Faisons un essai avec un TFD à 4 points. On peut utiliser l'algo FFT de calcul rapide de TFD d'octave/matlab...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "operational-steal",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "clear all;\n",
    "close all;\n",
    "\n",
    "N=4;\n",
    "Fe=100; Te = 1/Fe;\n",
    "%% VOTRE CODE : k= ? et t= ?\n",
    "% déclarant un vecteur ligne k du segment entier [[0, 4[[ \n",
    "% le vecteur temps correspondant\n",
    "\n",
    "\n",
    "\n",
    "%% VOTRE CODE : n= ? et f= Df=?\n",
    "% déclarant un vecteur ligne n du segment entier [[0, 4[[ \n",
    "% le vecteur des 4 fréquences de [0 ; Fe[ <- semi ouvert ! \n",
    "% fréquences espacées de Df (la résolution fréquencielle ou le F0 des SdF)\n",
    "\n",
    "\n",
    "\n",
    "%% VOTRE CODE : delta = @(?) ?\n",
    "% qui crée une fonction anonyme delta(k) qui vaut (k==0) \n",
    "\n",
    "\n",
    "\n",
    "l=2;\n",
    "%% VOTRE CODE : x0 = ?; xl = ?\n",
    "% qui utilise la fonction delta pour calculer x0 = delta(k)\n",
    "% et xl = delta(k) retardée de l\n",
    "\n",
    "\n",
    "%% VOTRE CODE : stem (?,?,'r')\n",
    "% qui affiche x0 et xl avec \"stem\" et \"hold on\" \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5c3df387-977a-45ab-8b79-d167f41a02c5",
   "metadata": {},
   "source": [
    "On peut maintenant faire `help fft` dans la console (File → New Console for Notebook)\n",
    "\n",
    "Et afficher les TFD des signaux : Attention ils **sont complexes** "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82940bdb-637a-4ec5-8da0-991490435167",
   "metadata": {},
   "outputs": [],
   "source": [
    "%% VOTRE CODE : tfd_x0 = ? et tfd_xl = ?\n",
    "\n",
    "\n",
    "% On affiche les vecteurs\n",
    "disp(tfd_x0)\n",
    "disp(tfd_xl)\n",
    "\n",
    "% On fait un graphe du spectre (fréquenciel)\n",
    "stem(f,real(tfd_xl),'bs'); hold on;\n",
    "stem(f,imag(tfd_xl),'go');\n",
    "legend([\"reel\"; \"imaginaire\"])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "vocal-courtesy",
   "metadata": {},
   "source": [
    "## TFD Matricielle\n",
    "---\n",
    "\n",
    "On peut créer une matrice avec en rangées tous les vecteurs de la base fréquencielle $w_0$ à $w_{N-1}$ et en colonne le temps pour $k$ de 0 à $N-1$. Puis utiliser cette matrice pour calculer tous les produits scalaires d'un signal avec les vecteurs $w_n$ : ce sont les coordonnées fréquentielles du signal donc la TFD.\n",
    "\n",
    "### Construction de la matrice des vecteurs fréquenciels $w_n$\n",
    "---\n",
    "Inspirez-vous de ce qui est fait avec 3 échantillons dans [VEC2 TFD matriciel](../..//cours/notebooks/VEC2_bases_frequentielles.ipynb#Base-fr%C3%A9quentielle-orthogonale-:-TFD)\n",
    "\n",
    "Construisez une représentation matricielle contenant chaque vecteur $w_i$ de la base en colonne : \n",
    "  * la matrice `s` avec un échantillonnage presque continu du temps k (en rangées) \n",
    "  * la matrice `W` avec seulement 4 échantillons temporels (en rangées toujours)\n",
    "les colonnes sont liées aux fréquences :"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9d95790-bca5-44eb-a5bd-c9828eb62ab4",
   "metadata": {},
   "outputs": [],
   "source": [
    "%% MODIFIEZ LE CODE ISSU DE VEC2 : pour passer en 4 échantillons\n",
    "N = 3;  \n",
    "t=-0.5:0.01:3.5 ; \n",
    "k = 0:2 ;  \n",
    "\n",
    "Fe= 1; \n",
    "Df = Fe/N; % On peut le nommer f0 pour évoquer les SdF\n",
    "n = 0:2;  % ICI encore\n",
    "f = n*Df;\n",
    "\n",
    "s  = exp(i*2*pi*Df*t'*n);  \n",
    "W =  exp(i*2*pi*Df*k'*n);\n",
    "\n",
    "addpath(\"../../cours/notebooks/utiles\")  \n",
    "\n",
    "\n",
    "%% MODIFIEZ ICI !\n",
    "\n",
    "%avoir=1; %% CHANGEZ le nombre pour sélectionner w_i à afficher\n",
    "avoir=0:2; %% DECOMMENTEZ pour tout afficher\n",
    "\n",
    "affiche3d(t,s(:,(avoir+1)),k,W(:,(avoir+1)));  \n",
    "subplot(2,2,1); title([\"Affichage de w_n pour n =  \" num2str(avoir)])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "arctic-spank",
   "metadata": {},
   "source": [
    "Vérifions que l'on a bien une base orthonormée.\n",
    "\n",
    "Pour cela on calcule le produit scalaire discret avec les vecteurs coordonnées  $<\\!<u,v>\\!> = \\overline{{}^T\\!{V}}. U$\n",
    "\n",
    "La transposée du conjugué est dite la **transposée Hilbertienne** et est notée $V^H=\\overline{{}^T\\!{V}}$.\n",
    "\n",
    "Avec octave on utilise l'apostrophe `V'`=$V^H$\n",
    "\n",
    "Si l'on calcule $W^H.W$ nous obtenons ainsi tous les produits scalaires de la base.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d63d550-cf15-4855-8ccb-2ca3ac3e3d4a",
   "metadata": {},
   "outputs": [],
   "source": [
    "% On extrait deux vecteurs\n",
    "w0 = W(:,1); \n",
    "w2 = W(:,2);\n",
    "\n",
    "arrondi = @(x) round(x*100)/100;\n",
    "\n",
    "%% VOTRE CODE : scal_w0_w2 = ?  norme_w2 = ?\n",
    "%calculant les produits scalaires <w_0,w_2> et <w_2,w_2>\n",
    "\n",
    "\n",
    "\n",
    "arrondi(scal_w0_w2)\n",
    "arrondi(norm_carre_w2)\n",
    "%Tous les produits scalaires\n",
    "scalaires = W'*W;\n",
    "arrondi(scalaires)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc9ca392-e4f4-49a6-8b1c-32fd78fa8445",
   "metadata": {},
   "source": [
    "La matrice $W^H=(\\overline{W_N^{nk}})_{0<n,k<N-1}$, est bien orthogonale (dans $\\mathbb{C}$) avec\n",
    "$$ W. W^H = N . Id$$ mais pas normée.\n",
    "\n",
    "Avec $W^H$ dite la matrice adjointe ou transposée de Hilbert : transposée du conjugué.\n",
    "\n",
    "Et donc $W^{-1}=W^H/N$ \n",
    "\n",
    "Donc on peut calculer une tranformée en fréquence avec $W^H$ et son inverse avec $W/N$ :\n",
    "\n",
    "$s \\overset{TFD}{\\longrightarrow} \\hat{S} = W^H.s \\overset{TFD^{-1}}{\\longrightarrow} s = \\frac{W}{N} \\,.\\, \\hat{S} = \\frac{W}{N} \\,.\\, W^H.s = Id.s = s $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "630d2965-2b07-46bb-9a81-53eb90d0a0f7",
   "metadata": {
    "tags": []
   },
   "source": [
    "### Calcul de la TFD en matriciel\n",
    "\n",
    "On utilise alors la matrice $W^H$ pour passer du temporel vers le fréquentiel :\n",
    "\n",
    "$W^H . s  =  \\hat{S}$\n",
    "\n",
    "Dans la suite, nous recalculons $M=W^H$ directement avec $W^H = \\left({W_N}^{n.k}\\right)_{k,n}$ (maintenant k en rangée et n en colonne pour **transposer**) avec $W_N=e^{-i\\frac{2\\pi}{N}}$ le **conjugué** de la racine Nème de l'unité. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "enormous-housing",
   "metadata": {},
   "outputs": [],
   "source": [
    "N=4;\n",
    "W_N = exp(-i*2*pi/N); % on prend directement le conjugué\n",
    "% cette fois-ci k est horizontal et n vertical\n",
    "% car on veut directement la transposée\n",
    "k=0:(N-1);\n",
    "n = k;\n",
    "\n",
    "M_tfd = W_N.^(n'*k);\n",
    "\n",
    "%% VOTRE CODE : verif que M = W^H noté W' en octave\n",
    "% utilisez la fonction sum() pour avouter les valeurs absolues \n",
    "%  des erreurs d'une matrice.\n",
    "\n",
    "\n",
    "\n",
    "M_tfdi = M_tfd'/N; % transformée inverse par simple transposé de hilbert\n",
    "tfdi_de_tfd = arrondi(M_tfdi * M_tfd) % On doit avoir la matrice identité\n",
    "\n",
    "\n",
    "%% VOTRE CODE : x0 = , xl = , x = \n",
    "% avec x le signal (1,-1,1,-1). en octave [1,2] est une matrice ligne et [1; 2] colonne\n",
    "% attention le vecteur signal doit être vertical !\n",
    "% transpose(x) est différent de x' !\n",
    "\n",
    "\n",
    "\n",
    "disp([\"  x0\", \"  xl\", \"   x\"])\n",
    "disp(arrondi([ x0, xl, x]))\n",
    "\n",
    "% VOTRE CODE : tfd_x0 =  ,  tfd_xl =  ,tfd_x =   ,  \n",
    "% calculer matriciellement les transformées \n",
    "\n",
    "\n",
    "\n",
    "disp([\"tfd(x0\", \" xl\", \"  x )\"])\n",
    "disp(arrondi([ tfd_x0, tfd_xl, tfd_x]))\n",
    "\n",
    "%% VOTRE CODE : fft_x0, fft_xl, fft_x\n",
    "% calculer avec l'algo fft de octave  les transformées \n",
    "\n",
    "\n",
    "\n",
    "disp([\"fft(x0\", \" xl\", \"  x )\"])\n",
    "disp(arrondi([ fft_x0, fft_xl, fft_x]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12ffc20c-936c-45a3-8775-e1be56dd995d",
   "metadata": {},
   "source": [
    "On peut constater que la base fréquencielle n'est pas normée (vecteurs de norme $\\sqrt{N}$, car $<< \\vec{w_n}, \\vec{w_n}>> = N$), la norme du signal mesurée en fréquenciel $\\sqrt{<\\hat{S},\\hat{S}>}$ est amplifié. Coefficients $\\sqrt{N}$ fois trop grands et norme du vecteur fréquenciel $\\sqrt{N}$ fois trop grande. On ne divise pas par N par convention permettant d'économiser du temps de calcul. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3cb4b03d-3381-4690-9a3d-e47bfdaac3ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "norme = @(x) sqrt(x'*x);\n",
    "norme(fft_x0)/norme(x0)\n",
    "norme(fft_xl)/norme(xl)\n",
    "norme(fft_x)/norme(x)\n",
    "\n",
    "sqrt(N)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55ece506-06e8-4961-b7bb-b31fac890760",
   "metadata": {},
   "source": [
    "L'idéal serait de diviser par N d'abord pour être normé et avoir : \n",
    "\n",
    "$s \\overset{\\text{TFD normée}}{\\longrightarrow} \\hat{S} = \\frac{W^H}{N}.s  \\overset{\\text{TFD norm}^{-1}}{\\longrightarrow} s =  W. \\hat{S}$ \n",
    "\n",
    "Mais on a une TFD non normée qu'il faut compenser en réciproque en divisant par N :\n",
    "\n",
    "$s \\overset{\\text{TFD}}{\\longrightarrow} \\hat{S} = {W^H}.s  \\overset{\\text{TFD}^{-1}\\text{compensée} }{\\longrightarrow} s =  \\frac{W}{N}. \\hat{S}$ \n",
    "\n",
    "> Les coefficients de la TFD sont donc laissé N fois trop grand ce qui fait \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9b4ee311-3b15-4c0c-b861-4764ea5cad65",
   "metadata": {},
   "outputs": [],
   "source": [
    "norme(fft_x0/N)/norme(x0)\n",
    "norme(fft_xl/N)/norme(xl)\n",
    "norme(fft_x/N)/norme(x)\n"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "4.2.2"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}