modelisationanalysespectrale/duval/nonParametrique/analyze.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Chargement des données"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib widget\n",
    "import pickle\n",
    "import scipy.signal\n",
    "\n",
    "\n",
    "my_data = np.genfromtxt('../Supelec_2012_SIR_Spectral_Analysis_EA_v001.csv', delimiter=',')\n",
    "# c'est plus rapide de charger des pickle que des csv\n",
    "a=\"\"\"\n",
    "with open(\"data.pickle\", \"wb\") as file:\n",
    "    pickle.dump(my_data, file)\n",
    "\n",
    "with open(\"data.pickle\", \"rb\") as file:\n",
    "    my_data = pickle.load(file)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOy9eYwkWXrY93txVR5VeVRW1tFVfVQf0z0z3Ts7M7uzXA7FNbnLa6W1RF2mJEiEbHBt2NZflmX9YRgrwLZk2BAMGYYByjRMA5JoUict00uRS2lWu5zduXpmu3v6qu7q6qrqOrKOzMo743j+44uIjIzK6hmL1miE7QcUKjMy4h3f+953f18orTXP2rP2rD1rz9oPVzP+TU/gWXvWnrVn7Vn75Nsz4v+sPWvP2rP2Q9ieEf9n7Vl71p61H8L2jPg/a8/as/as/RC2Z8T/WXvWnrVn7YewPSP+z9qz9qw9az+E7RMn/kqpn1VK3VVKrSil/uonPf6z9qw9a8/aswbqk4zzV0qZwD3gp4AN4G3gz2itP/zEJvGsPWvP2rP2rH3ikv9rwIrW+qHWegD8GvBHP+E5PGvP2rP2rP3QN+sTHm8RWE983wC+kLxBKfV14OsAJuarOVVI/GagdSBfTlJYVPibSlzToAwDHQQn339SH+n/yZa89nHuTz+TvBYNbQg/HlmnGn9vdD2GS7TO8Nnk9ZF1pufJGPik7z+ppeczZn5Phd1HrG9kz8c9T3hPEBwbI77O+HtH4JOcv0r1kRw/Ca/U9REYKlCGiQ784+tNtfi5JLySv6XXn8a7xPrSY33kviY+H4NXaq1jYTFmX5VloX3vZPxJzX/sGqNHx8wpPiPpPU/h1InzHYd7Y5556v6Ne/7Y5Mevd+zZetpZO2m+48YJvzf14Z7Wujqmt7h90sT/I5vW+peBXwYoqGn9Bb4MhAjlecP/9vA7cOxz9MxIMxPjpO5RjkPQ6Ry7pgcDuTciHokxjHwu/l3Zx+cxsq7k3FNzjuYCoMzUGmwrnsdT1xb1Z1rxWqLP0bMnzSMJH+U4ct+YZ076PgLTMcgbwWbk+QS8jFxO5ntSnzByT/p568xp/M2t489kc/FaYhhrhrhjWsf2Nrm2GI6JPYqejQmjlegjIiDJdYS6dYxfqbGMnOBQhJsRXsW4aFrH8CQNzxNhFuKnyg5xO4Jj3ML5RnhiFPP49caxvUjjUXKs5Fri66l7R3A8cXZPOndJfNA6MUa4pvT1kTNln9BPAt5JGjIyJ/M4no/M0R4dx5gM4TXm/KfP2gjdso+f83GwjfYoOcf0/cfusyx+x/21tWMdptonbfbZBE4nvi+F18a39EFidEOS14xcDu15YwlV9D/6i1rymQhJomswJBrKsuK/ZB9BpzOWgMfzSMxx3NxjZE/8T/8WzS0973QbQbrBIJ7TyEFPgjbJ4JJzHwzGMpqYKSTWm4ZpEk4nwSwN92jM9FqiucV9hGtKPh+1JOGP5mGWisf6NXK544cxNc/kXibhmHw2Db/k9xECmJhj0OnEfST3O9qfJO5GuBjdZ5aK8RzMUnEsrJJ9JMdMEu0kHCNcT+KM9jz8euNYv9Hv6XGSZyuadxp3xrVxuJre02RLnueg08GYzMe4kWY6yXOT3L/ot+RZSu7/yLX0GRwMRvY/Of8kvNK0Z9yaxjGHcWckifvRXJLrSO7BSfd9VPukJf+3gUtKqWWE6P8C8Gef9kD6QCYll5EDGyF1iOxpYKQ3GBghjOmNPYmJpLWC5PPj5jRCYMJn0pJOcl3jmElaoxlH0I8hVIo5JscZOcipQx0dpOQYaSQbN7dkS8Ju3DzGSkRjkDYptY1oPyminF6HfHdH1pHej3Fzj2Az9nCG/R/bp8T98dxSTD/52ziGmoZREleAWBKPGPNJ8Ev2Gz2bZiQj2mxybk85Z8mW7CsWulL9pc9i8tn4c0qbTTL8NNGOno3HCvc2ElzS2tkIM0jtRRpu6XMdzX8cTNL3pfuNxosY/Tg4ps9FmqAP522P9PNR40e/pYW5p7VPVPLXWnvAfwr8NnAb+HWt9a2P+3yScCU3epw0Pg4AaQ6ZvH6SFHBsDd5QU/goII8jZumxk79FCG5M5oefU1JNdD197SRJLS2VJ4lTWsKJ5pGUDpO/JxFr3Bqi341c7pj5SDkOZqk40m+0jmPjpH8PD1R6XyJJzqhWjhGQk2CdZqRp3EnDJS2NJSX4cfOJJM6IID4dXnZ8f9TPOM0hCYdx96QFjzRuJglgWnBKtjQBj+5L78s4STqJI9E8xwpHCbw4Ju0mGHR67BE4hQLUOCEwucZ4L7zjWuqxcxEJONHcU4JK8pmTfkviXXSOkowweW+swVhDrSeJI0mGPzLfEwSL6LdxprmT2idu89da/xbwWx/vZkZtpOFGJrnpOEktDUhliW0uaLWPScYjCD/mYEVjpSWccZpCuqXvOem+ZAs6HegwItEmpbUkQoxba3qe4yTj5OFIzze59mSfx0xVY5Aw2ou0/T4+TJ3UHMfAM+o7TVTS8xx5rnE0MreTTA/J/RwnKcZwHiMZJ+E2YiJMHO5j/oHU3NNSvR64x+6JcC35Pw3fk7SLk+CV3I9x/pWI4CXnlmT0x4jcCcJSNEdjMo8aPMWGHY07DkdSaxlHfGOcHjO/9FjHJP/UWUqPlZbS036D5PWRccYIlcm+x2lREYEf8WWlGPkxGCTOcVqLH9Hy3GPDHWuf+gzfcRLWOOJ1TIJISXZBqx3/BqP2/FhataxjfccbcwJHTUupyc/jDkgkzScJTtRP8nP60EfPjkgsjnNMakgTzUhKGidlJf/SMEzap2OCmCCqaeky3f+4PYn3JiURpolJ+qBENteTxkn2H0lg6cOaHDtNoCMpLTlu2j5/7JCnHekpU8E4mB5jTmO0syRDSZtF0ucgDeenMe+4n5RJFBiLa+MY6NOY8cj3YmHkWe15sZSb7CtNvMedpTSBBjDC/scR3DRep/E08oMk+x2ngcYWhRPgkMbRpNSfvj/Z0mONI/hRf8nfk/CK/R0fIwjkae3TTfzVeGlhHBKlD9xJnDrdV5Izn0QcxxFOGD1U0cab1eoxxpJsI+pg2pyRIEwx0U0ysZSWkzTRwBAp0uOeRBSSRCVC+OQ84+fHONeT8Is+pxlR+vtJ0tlJxDx5/zjCnLaHJuF4rL9Iuk2sZZyUle4vCZ/k+GkTUXrO40xq49o4opKGR9IxG9+bYjbpOSf3ZpwJbCxT8ryxa34a0U4zpmQbZbpubBIat/7kmE/zaQF4W9sCl3J5FPfHSPkxEx9j+oLRvUzi9tPMJycJLOMYSnr+6cCF5Fyjax8l3CQ1vyRDeNp5Gtc+3cQ/bEkikpSgozYOqSNkip917JHrJzGVpxH79Odj3N+y8Gu1Y5uQnHtyDWlinrz3JBikiW1y3ZEd8SSYJNeZvCc6HDHTSmgU6XWOk9xjmKak6rGOu7S/IS2lJ6I4kog9wnDGOHfTeJGEUzS32OSRGnNk7PAvGRqZPpBjJeUUEY/+n2SjTR7WJH6ncTDWVlJM5tj6UriR/Jwm7uOeS449Ms+UQDWCM2POhfY8gtr+iO08KagkGVHapzWurwiGyXmbF5flt3Z7OGbCB5Dcm2i/TwruSO9fGreSuPJRsEnjyAgzGYMTaZNhevz03sEQv5NCQpr5jIPpuPapJ/4nEZt0FEP6magN1Vh3ZHPguC0Tjku5yf/pz2nunyRSUUszqXEIkZzzOIkhWm/yvnHzVJYVm7ei70ltIE3EkvMdZy5KrmEEaZPMbIxz9Glhe9E4aXNacp3HpNmE4yzZYh9QygyQ3JNkX8kDliQqSYkz7eiMiMq4/k6SnpOHcQR3E+ON9Jky76SZU5KxJsdLjj9iphtD4JNrSrZxWkX62bQwkBZEjml8xcIIvo47Y9Ga0ud4nKCSlqqDR+tj6cJJ60jb1NPn+RhepRjvOLqRvA8YMWsdw92UwJZmUuOEyLSjOq2RJMc6Jmx+TKfvp574fxTBflool7IszMWFY5sSSRzjDtNIeGFKyklvehqZ0oQiltrS45yAUCc
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure()\n",
    "plt.imshow(np.abs(my_data), aspect=\"auto\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Constantes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# paramètres du signal\n",
    "Te = 10**(-6)\n",
    "Fe = 1/Te\n",
    "\n",
    "def make_cos_window(te, Tcut, length):\n",
    "    dT=Tcut/16\n",
    "    len1 = int((Tcut-dT)/te)\n",
    "    len2 = int(dT/te)\n",
    "    len3 = length - len1 - 2*len2\n",
    "    cospart = np.cos(np.linspace(0, np.pi, len2))*0.5 + 0.5\n",
    "    return np.concatenate([\n",
    "        cospart[::-1],\n",
    "        np.ones(len1),\n",
    "        cospart,\n",
    "        np.zeros(len3),\n",
    "    ])\n",
    "\n",
    "# fenêtre cosinus\n",
    "cos_window = make_cos_window(Te, 0.6e-3, 4096)\n",
    "# coefficients pour filtre de chebychev\n",
    "chebyfilt = scipy.signal.cheby2(2, 40, [2e3, 2e5], fs=Fe, btype='bandpass', output='sos')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Aide à l'affichage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_fft(sig, te, absolute=True):\n",
    "    N = len(sig)//2\n",
    "    signal_f = np.fft.fft(sig)[:N]\n",
    "    if absolute:\n",
    "        signal_f = abs(signal_f)\n",
    "    freqs = np.arange(0, N)*0.5/te\n",
    "    return freqs, signal_f\n",
    "\n",
    "def plot_fft(sig, te, absolute=True):\n",
    "    plt.plot(*get_fft(sig, te, absolute))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analyse d'un signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sig_simple(signal):\n",
    "    time = np.arange(0, len(signal))*Te\n",
    "    signal_win = cos_window * signal\n",
    "    signal_filt = scipy.signal.sosfilt(chebyfilt, signal_win)\n",
    "\n",
    "    freqs, sig_tf = get_fft(signal_filt, Te)\n",
    "    \n",
    "    energy_f = np.multiply(sig_tf, sig_tf)\n",
    "    bary_f = sum(freqs * energy_f) / sum(energy_f)\n",
    "    var_f = np.sqrt( sum((freqs - bary_f)**2 * energy_f) / sum(energy_f) )\n",
    "\n",
    "    plt.figure()\n",
    "    plt.plot(time, signal, linewidth=1)\n",
    "    plt.plot(time, signal_win, linewidth=1)\n",
    "    plt.plot(time, signal_filt, linewidth=1)\n",
    "\n",
    "    plt.figure()\n",
    "    plot_fft(signal, Te)\n",
    "    plot_fft(signal_win, Te)\n",
    "    plot_fft(signal_filt, Te)\n",
    "    plt.vlines([bary_f-var_f, bary_f, bary_f+var_f], 0, max(sig_tf), colors=['r', 'g', 'r'])\n",
    "\n",
    "signal = my_data[:,101]\n",
    "sig_simple(signal)    \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analyse de tous les signaux, pour extraire les paramètres significatifs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def extract_params(signal):\n",
    "    # fenetrage\n",
    "    signal_win = cos_window * signal\n",
    "    \n",
    "    # filtrage\n",
    "    signal_filt = scipy.signal.sosfilt(chebyfilt, signal_win)\n",
    "\n",
    "    # fft\n",
    "    freqs, sig_tf = get_fft(signal_filt, Te)\n",
    "    \n",
    "    # moyenne, variance, moments d'ordre 3 à 10 (pondérées)\n",
    "    energy_f = np.multiply(sig_tf, sig_tf)\n",
    "    bary_f = sum(freqs * energy_f) / sum(energy_f)\n",
    "    var_f = sum((freqs - bary_f)**2 * energy_f) / sum(energy_f)\n",
    "    \n",
    "    moments = [sum((freqs - bary_f)**n * energy_f) / sum(energy_f) for n in range(3, 11)]\n",
    "\n",
    "    return [bary_f, var_f] + moments\n",
    "\n",
    "\n",
    "barys_vars_f = []\n",
    "for i in range(my_data.shape[1]):\n",
    "    barys_vars_f.append(extract_params(my_data[:, i]))\n",
    "\n",
    "\n",
    "vars_f, barys_f, mom3, mom4, mom5, mom6, mom7, mom8, mom9, mom10 = zip(*barys_vars_f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.scatter(vars_f, barys_f, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moyenne fréquentielle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.scatter(vars_f, mom3, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moment d'ordre 3\")\n",
    "\n",
    "_=\"\"\"\n",
    "plt.figure()\n",
    "plt.scatter(vars_f, mom4, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moment d'ordre 4\")\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(vars_f, mom5, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moment d'ordre 5\")\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(vars_f, mom6, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moment d'ordre 6\")\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(vars_f, mom7, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moment d'ordre 7\")\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(vars_f, mom8, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moment d'ordre 8\")\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(vars_f, mom9, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moment d'ordre 9\")\n",
    "\n",
    "plt.figure()\n",
    "plt.scatter(vars_f, mom10, marker='.')\n",
    "plt.xlabel(\"Variance fréquentielle\")\n",
    "plt.ylabel(\"Moment d'ordre 10\")\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}