{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Kolokwium poprawkowe  \n",
    "27.01.2020\n",
    "\n",
    "### Proszę zmienić nazwę tego notebooka na następującą postać: AS_Kolokwium_Poprawkowe_numer_indeksu.ipynb, wpisując we właściwym miejscu własny numer indeksu.\n",
    "\n",
    "_Rozwiązania proszę odesłać na adres jarekz@fuw.edu.pl_\n",
    "\n",
    "\n",
    "Importy i funkcje, które przydadzą się do rozwiązania zadania, treści zadań znajdują się poniżej."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pylab as plt\n",
    "import scipy.signal as ss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def charkterystyki(a,b,f,T,Fs):\n",
    "    t = np.arange(-T, T, 1/Fs)\n",
    "    \n",
    "    f, h = ss.freqz(b, a, worN= f, fs = Fs) \n",
    "    m = np.abs(h)\n",
    "    \n",
    "    f, grupowe = ss.group_delay((b,a), w=f, fs = Fs)\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    plt.subplot(2,2,1)\n",
    "    plt.title('moduł transmitancji')\n",
    "    plt.plot(f,20*np.log10(m))\n",
    "    plt.ylabel('[dB]')\n",
    "    plt.ylim([-20, 1])\n",
    "    plt.grid(True)\n",
    "    \n",
    "    plt.subplot(2,2,3)\n",
    "    plt.title('opóźnienie grupowe')\n",
    "    plt.plot(f,grupowe)\n",
    "    plt.ylabel('próbki')\n",
    "    plt.xlabel('Częstość [Hz]')\n",
    "    plt.grid(True)\n",
    "    plt.ylim([0, np.max(grupowe)+1])\n",
    "    \n",
    "    plt.subplot(2,2,2)\n",
    "    plt.title('odpowiedź impulsowa')\n",
    "    x = np.zeros(len(t))\n",
    "    x[len(t)//2] = 1\n",
    "    y = ss.lfilter(b,a,x)\n",
    "    plt.plot(t, x)\n",
    "    plt.plot(t, y)\n",
    "    plt.xlim([-T/2,T])\n",
    "    plt.grid(True)\n",
    "    \n",
    "    plt.subplot(2,2,4)\n",
    "    plt.title('odpowiedź schodkowa')\n",
    "    s = np.zeros(len(t))\n",
    "    s[len(t)//2:] = 1\n",
    "    ys = ss.lfilter(b,a,s) \n",
    "    plt.plot(t, s)\n",
    "    plt.plot(t, ys)\n",
    "    plt.xlim([-T/2,T])\n",
    "    plt.xlabel('Czas [s]')\n",
    "    plt.grid(True)\n",
    "    \n",
    "    fig.subplots_adjust(hspace=.5)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def parametryAR(x,p):\n",
    "    '''funkcja estymująca parametry modelu AR \n",
    "    argumenty:\n",
    "    x - sygnał\n",
    "    p - rząd modelu\n",
    "    f. zwraca:\n",
    "    a - wektor współczynników modelu\n",
    "    epsilon - estymowana wariancja szumu\n",
    "\n",
    "    funkcja wymaga zaimportowania modułu numpy as np\n",
    "    '''\n",
    "    N = len(x)\n",
    "    ak = np.correlate(x,x,mode='full')\n",
    "    norm_ak = np.hstack((np.arange(1,N+1,1),np.arange(N-1,0,-1)))\n",
    "    ak=ak/norm_ak\n",
    "    R=ak[N-1:]\n",
    "    RL  = R[1:1+p]\n",
    "    RP = np.zeros((p,p))\n",
    "    for i in range(p):\n",
    "        aa = ak[N-1-i:N-1-i+p]\n",
    "        RP[i,:] = aa\n",
    "    a=np.linalg.solve(RP,RL)\n",
    "    sigma2 = (ak[N-1] - np.sum(a*ak[N:N+p]))\n",
    "    return a, sigma2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def TFRPlot(TFR, t_mapy, f_mapy, sig, Fs=128,title =''):\n",
    "    '''\n",
    "    Funkcja do rysowania map czas-częstość z sygnałem zaprezentowanym poniżej\n",
    "    TFR - mapa czas-częstość (time-freqyency representation\n",
    "    t_mapy, f_mapy - wektory reprezentujące osie czasu i częstości\n",
    "    sig - sygnał do wyrysowania pod mapą (np. ten, z którego powstała mapa)\n",
    "    Fs - częstość próbkowania sygnału \n",
    "    title - tytuł do wyświetlenia ponad mapą\n",
    "    '''\n",
    "    df = f_mapy[1]-f_mapy[0]\n",
    "    dt = t_mapy[1]-t_mapy[0]\n",
    "    t = np.arange(0,len(sig))/Fs\n",
    "    \n",
    "    plt.figure()\n",
    "    sygAxes = plt.axes([0.05, 0.05, 0.8, 0.1])\n",
    "    tfAxes = plt.axes([0.05, 0.15, 0.8, 0.8])\n",
    "    sygAxes.plot(t,sig)\n",
    "    plt.xlim((t_mapy.min(), t_mapy.max()))\n",
    "    #plt.setp(sygAxes, yticklabels=[])\n",
    "    tfAxes.imshow(TFR,aspect='auto',origin='lower',interpolation='nearest', \n",
    "                  extent=(t_mapy.min()-dt/2,t_mapy.max()+dt/2,f_mapy.min()-df/2,f_mapy.max()+df/2))\n",
    "    plt.setp(tfAxes,xticklabels=[])\n",
    "    plt.title(title)\n",
    "    plt.show()\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zadanie 1.  Procesy AR (5pkt)\n",
    "\n",
    "#### Proszę wczytać plik \n",
    "https://www.fuw.edu.pl/~jarekz/KOLOKWIUM_AS/zad1_syg.bin \n",
    "\n",
    "Jest to sygnał jednokanałowy, liczby są zapisane w formacie ‘<f’, próbkowanie 128Hz. Proszę wczytać sygnał i wykreślić jego przebieg."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę napisać funkcję tworzącą wykres AIC. (1pkt)\n",
    "Dla przypomnienia, kryterium Akaike (AIC):\n",
    "\n",
    "  $\\mathrm{AIC}(p)= \\frac{2p}{N} +\\ln(V) $\n",
    "\n",
    "$p$ - liczba parametrów modelu, $N$ - liczba próbek sygnału,$V$ - wariancja szumu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę wybrać optymalny rząd na podstawie wykresu i uzasadnić swój wybór (1pkt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Miejsce na komentarz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę wyestymować współczynniki modelu AR dla wybranego rzędu, a następnie obliczyć widmo modelu AR (1pkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę skomentować otrzymane wyniki (1pkt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zadanie 2. Filtrowanie sygnałów (6pkt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę wczytać plik \n",
    "https://www.fuw.edu.pl/~jarekz/KOLOKWIUM_AS/zad2_syg.bin\n",
    "\n",
    "Jest to sygnał jednokanałowy, liczby są zapisane w formacie ‘float64’, próbkowanie 128Hz. Proszę wykreślić jego przebieg. Jest to sygnał złożony z 2 gaborów o częstości 10Hz i skali 1s, delty oraz szumu."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Skonstruuj dwa projekty filtra, o zadanych parametrach (2pkt):\n",
    "- Filtr Chebyshev II dla pasma przenoszenia 9-11Hz, z pasmem przejściowym o długości 1Hz, tak by minimalne tłumienie w paśmie zaporowym wynosiło 25 dB.\n",
    "\n",
    "- Filtr pasmowo przepustowy Butterworth’a 3 rzędu z wartościami odcięcia 8 i 12 Hz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Wraz z projektem filtra proszę wyrysować jego charakterystyki (0,5pkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę zastosować filtr do wczytanego sygnału (1pkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę skomentować wynik (3pkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zadanie 3. Mapy czas częstość (5pkt)\n",
    "Przydatna funkcja:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę uzupełnić funkcje spektrogram (2pkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spektrogram(x, okno, trans , Fs):\n",
    "    Nx = len(x)\n",
    "    No = len(okno)\n",
    "    okno = #unormuj okno\n",
    "    pozycje_okna = np.arange(0, Nx, trans)                        \n",
    "    t = pozycje_okna/Fs\n",
    "    N_trans = len(pozycje_okna)    \n",
    "    f = np.fft.rfftfreq(No, 1/Fs)\n",
    "    P = np.zeros((len(f),N_trans))\n",
    "    z = np.zeros(int(No/2))\n",
    "    sig = np.concatenate((z,x,z))\n",
    "    for i,poz in enumerate(pozycje_okna):\n",
    "        s = sig[poz:...]\n",
    "        s = s*okno \n",
    "        S = ...\n",
    "        P_tmp = ...   \n",
    "        P_tmp = P_tmp.real/Fs \n",
    "        if len(s)%2 ==0: # dokładamy moc z ujemnej części widma \n",
    "            P_tmp[1:-1] *=2\n",
    "        else:\n",
    "            P_tmp[1:] *=2            \n",
    "        P[:,i] = P_tmp\n",
    "    return t,f,P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  Proszę wczytać plik o nazwie \n",
    "https://www.fuw.edu.pl/~jarekz/KOLOKWIUM_AS/zad3_syg.bin\n",
    "\n",
    "#### Jest to sygnał jednokanałowy, liczby są zapisane w formacie ‘float’, próbkowanie 128Hz. Proszę wykreślić jego przebieg."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę wywołać spektrogram dla sygnału 'zad3_syg.bin' dla różnych długości okien i skomentować jaki mają wpływ na rozdzielczość czasową i przestrzenną (2pkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Zadanie 4. Własności okienek i periodogram (5pkt)\n",
    "\n",
    "#### Proszę wczytać plik \n",
    "\n",
    "https://www.fuw.edu.pl/~jarekz/KOLOKWIUM_AS/zad4_syg.bin\n",
    "\n",
    "#### Jest to sygnał jednokanałowy, liczby są zapisane w formacie  ‘<f’, próbkowanie 100Hz. Proszę wykreślić jego przebieg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Proszę wykreślić periodogram okienka i zokienkowanego sygnału dla okna prostokątnego i hamming’a w skali decybelowej; widma okienek wyestymować dla `nfft = 1024` (2pkt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Skomentowac różnice pomiędzy otrzymanymi widmami w oparciu i tw. o splocie (3pkt) \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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}