{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analiza wyników pomiaru amplitudy sygnału generatora\n",
    "\n",
    "Archiwa\n",
    "[pomiar_sterowany_komputerem.tar.gz](pomiar_sterowany_komputerem.tar.gz) i\n",
    "[pomiar_sterowany_komputerem.zip](pomiar_sterowany_komputerem.zip) zawierają\n",
    "komplet plików do interaktywnego uruchomienia przedstawionych przykładów.\n",
    "Polecamy Pythona w wersji\n",
    "[Anaconda](https://www.anaconda.com/products/distribution). \n",
    "\n",
    "Na komputerach w pracowniach Python Anaconda jest zapisany w katalogu\n",
    "/work/anaconda3-2024. Aby skorzystać z tej wersji należy wykonać następującą\n",
    "konfigurację w wierszu poleceń.\n",
    "\n",
    "```\n",
    "eval \"$(/work/anaconda3-2024/bin/conda shell.bash hook)\"\n",
    "```\n",
    "\n",
    "Konfiguracja nie jest trwała i obowiązuje do końca pracy powłoki. W przypadku\n",
    "jednoczesnego korzystania z kilku egzemplarzy powłoki konfigurację należy\n",
    "wykonać w każdej z nich."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import json, os, scipy.optimize\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.image as mpimg\n",
    "from matplotlib.backends.backend_pdf import PdfPages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Czytam dane zarejestrowane w czasie pomiaru."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('generator.json', 'r') as f:\n",
    "    results = json.load(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Liczba zarejestrowanych pomiarów."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "84"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Zarejestrowane przebiegi to ciągi napięć $(u_i)$. Do punktów $(t_i, u_i\\pm \\frac{1}{2}u\\_err_i)$ dopasowuję funkcję \n",
    "\n",
    "   $$a\\sin(2\\pi f t+\\phi)+b$$\n",
    "\n",
    "gdzie $f$ jest nominalną częstotliwością mierzonego sygnału, a $a$, $\\phi$, $b$ parametrami. Obliczam wartości $a$, $\\phi$, $b$ i odpowiadające im niepewności $\\Delta a$, $\\Delta\\phi$, $\\Delta b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_sine(dt, ydat, yerr, freq):\n",
    "    # Współrzędne dobieram tak, że początek układu znajduje się na środku\n",
    "    # ekranu oscyloskopu i pokrywa się z punktem wyzwalania.\n",
    "    xdat = np.linspace(-len(ydat)/2*dt, len(ydat)/2*dt, len(ydat),\n",
    "                       endpoint=False)\n",
    "\n",
    "    # Definicja doposowywanej funkcji.\n",
    "    def f(x, a, ph, b):\n",
    "        return a*np.sin(2*np.pi*freq*x+ph)+b\n",
    "\n",
    "    # Kandydat początkowy f(x)=(Vpp/2)sin(2πfx).\n",
    "    p0 = (max(ydat)-min(ydat))/2., 0., 0.\n",
    "    # Dopasowanie metodą najmniejszych kwadratów.\n",
    "    p1, pcov = scipy.optimize.curve_fit(f, xdat, ydat, p0, sigma=yerr)\n",
    "    return *p0, *p1, *np.sqrt(np.diag(pcov))\n",
    "\n",
    "def get_sine(dt, n, freq, a, ph, b):\n",
    "    # Współrzędne dobieram tak, że początek układu znajduje się na środku\n",
    "    # ekranu oscyloskopu i pokrywa się z punktem wyzwalania.\n",
    "    x = np.linspace(-n/2*dt, n/2*dt, n, endpoint=False)\n",
    "    return a*np.sin(2*np.pi*freq*x+ph)+b\n",
    "\n",
    "for e in results:\n",
    "    # Czytam dane z poszczególnych kanałów pomiarowych i z logicznego kanału\n",
    "    # math.\n",
    "    for src in ['ch1', 'ch2', 'ch3', 'ch4', 'math']:\n",
    "        if src in e:\n",
    "            c = e[src]\n",
    "            # Sprawdzam czy oscyloskop przesłał dane.\n",
    "            if not c['bytes']:\n",
    "                continue\n",
    "\n",
    "            # Dopasowuję sinusoidę a*sin(2*pi*freq+phi)+b\n",
    "            (c['first_guess_a'], c['first_guess_ph'], c['first_guess_b'],\n",
    "                c['fitted_a'], c['fitted_ph'], c['fitted_b'],\n",
    "                c['fitted_a_err'], c['fitted_ph_err'], c['fitted_b_err']\n",
    "                ) = fit_sine(c['xincrement'], c['u'], c['u_err'], e['freq'])\n",
    "\n",
    "            # Chcę aby współczynnik a nie był ujemny.\n",
    "            if c['fitted_a']<0:\n",
    "                c['fitted_a'] *= -1\n",
    "                c['fitted_ph'] += np.pi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tworzę dziennik [generator.pdf](generator.pdf) zawierający wykresy wszystkich pomiarów."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "with PdfPages('generator.pdf') as pages:\n",
    "    for e in results:\n",
    "        # Obliczam ilość zarejestrowanych kanałów.\n",
    "        n = 0\n",
    "        for src in ['ch1', 'ch2', 'ch3', 'ch4', 'math']:\n",
    "            if src in e:\n",
    "                c = e[src]\n",
    "                if c['bytes']:\n",
    "                    bytes_cnt = len(c['bytes'])\n",
    "                    n += 1\n",
    "        if n == 0:\n",
    "            continue\n",
    "\n",
    "        fig, ax = plt.subplots(n+1)\n",
    "\n",
    "        # Widok ekranu oscyloskopu.\n",
    "        ax[0].set_xlim(0, bytes_cnt)\n",
    "        ax[0].set_ylim(0, 255)\n",
    "        ax[0].set_xticks([int(e) for e in np.linspace(0, bytes_cnt-1, num=13,\n",
    "            endpoint=True)])\n",
    "        for idx, label in enumerate(ax[0].xaxis.get_ticklabels()):\n",
    "            label.set_visible(False)\n",
    "        ax[0].set_yticks([int(e) for e in np.linspace(0, 255, num=11,\n",
    "            endpoint=True)])\n",
    "        for idx, label in enumerate(ax[0].yaxis.get_ticklabels()):\n",
    "            label.set_visible(False)\n",
    "        for src in ['ch1', 'ch2', 'ch3', 'ch4', 'math']:\n",
    "            if src in e:\n",
    "                c = e[src]\n",
    "                # Anaconda /work/anaconda3 w OKWF nie rozumie jeszcze\n",
    "                # napisów f'...' Stosuję starszą funkcję format.\n",
    "                ax[0].plot(c['bytes'], label='{:s} {:.2e} V/div'.format(\n",
    "                    src.upper(), c['scale']))\n",
    "        ax[0].set_title('{}\\nf = {:.2e} Hz, {:.2e} s/div'.format(\n",
    "            e['time'], e['freq'], e['time_scale']))\n",
    "        ax[0].legend(loc='upper right')\n",
    "        ax[0].grid()\n",
    "\n",
    "        # Sygnały w poszczególnych kanałach.\n",
    "        ax_idx = 1\n",
    "        for src in ['ch1', 'ch2', 'ch3', 'ch4', 'math']:\n",
    "            if src in e:\n",
    "                c = e[src]\n",
    "                t = np.linspace(-len(c['u'])/2*c['xincrement'],\n",
    "                    len(c['u'])/2*c['xincrement'], num=len(c['u']),\n",
    "                    endpoint=False)\n",
    "                ax[ax_idx].set_xlim(0, len(c['u']))\n",
    "                ax[ax_idx].set_xticks([int(e) for e in np.linspace(0,\n",
    "                    len(c['u'])-1, num=13, endpoint=True)])\n",
    "                for idx, label in enumerate(ax[ax_idx].xaxis.get_ticklabels()):\n",
    "                    label.set_visible(False)\n",
    "                ax[ax_idx].plot(c['u'], label=src.upper())\n",
    "                ax[ax_idx].plot(get_sine(c['xincrement'], len(c['bytes']),\n",
    "                    e['freq'], c['first_guess_a'], c['first_guess_ph'],\n",
    "                    c['first_guess_b']), label='first guess')\n",
    "                ax[ax_idx].plot(get_sine(c['xincrement'], len(c['bytes']),\n",
    "                    e['freq'], c['fitted_a'], c['fitted_ph'], c['fitted_b']),\n",
    "                    label='fitted sine')\n",
    "                ax[ax_idx].set_ylabel('U [V]')\n",
    "                ax[ax_idx].legend(loc='upper right')\n",
    "                ax[ax_idx].grid()\n",
    "                ax_idx += 1\n",
    "\n",
    "        pages.savefig(fig)\n",
    "        plt.close(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usuwam pomiary, które zostały powtórzone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "results2 = []\n",
    "prev = None\n",
    "for curr in results:\n",
    "    if prev!=None and curr['freq']!=prev['freq']:\n",
    "        results2.append(prev)\n",
    "    prev = curr\n",
    "if prev!=None:\n",
    "    results2.append(prev)\n",
    "results = results2\n",
    "results2 = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Liczba pomiarów po odrzuceniu pomiarów, które zostały powtórzone."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "84"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  Wykres zmian napięcia w funkcji częstotliwości"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22920 characters]\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "freq = [e['freq'] for e in results]\n",
    "attenuation_osc = [e['vpp1']/16 for e in results]\n",
    "attenuation_fit = [e['ch1']['fitted_a']/8 for e in results]\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.set_xscale('log')\n",
    "ax.set_yscale('log')\n",
    "ax.set_xlim(freq[0], freq[-1])\n",
    "ax.plot(freq, attenuation_osc, label='oscilloscope')\n",
    "ax.plot(freq, attenuation_fit, label='fitted sine')\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "ax.set_title('Attenuation')\n",
    "ax.set_xlabel('frequency [Hz]')\n",
    "ax.set_ylabel('ratio')\n",
    "\n",
    "plt.show()\n",
    "plt.close(fig)"
   ]
  }
 ],
 "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
