{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div align=\"center\">\n",
    "    \n",
    "# ΠΑΝΕΠΙΣΤΗΜΙΟ ΚΡΗΤΗΣ\n",
    "# Τμήμα Επιστήμης Υπολογιστών\n",
    "\n",
    "## HY-215: Εφαρμοσμένα Μαθηματικά για Μηχανικούς\n",
    "## Εαρινό Εξάμηνο 2024-2025\n",
    "### Διδάσκων: Γ. Καφεντζής\n",
    "    \n",
    "## Ανάθεση: 19-5-2025\n",
    "## Παράδοση (ηλεκτρονικά, όπως αναφέρεται στο τέλος του αρχείου): \n",
    "## Τρίτη, 3-6-2025, 23:59:59\n",
    "<h1><center>5η Σειρά Εργαστηρίου</center></h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Σε κάθε άσκηση `Python`, εισάγουμε πρώτα τις βιβλιοθήκες που θα χρησιμοποιήσουμε.\n",
    "\n",
    "Επιλέξτε το παρακάτω κελί κώδικα και πατήστε Shift + Enter, ώστε να το εκτελέσετε. \n",
    "\n",
    "Το ίδιο θα κάνετε με οποιοδήποτε κελί κώδικα. Αν τροποποιήσετε κάποιο κελί κώδικα, θα πρέπει να το ξανατρέξετε για να δείτε τα αποτελέσματά του."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.signal import butter, freqz, lfilter, residue\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.io.wavfile as wav\n",
    "import numpy as np\n",
    "from scipy import signal\n",
    "from IPython.display import Audio\n",
    "from utils import plot_laplace_frequency_response"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1η Άσκηση\n",
    "\n",
    "Έχετε ήδη καταλάβει ότι το Ανάπτυγμα σε Μερικά Κλάσματα είναι μια σωτήρια 😃 τεχνική όσον αφορά αντίστροφους μετασχηματισμούς, τόσο Fourier όσο και Laplace. Φυσικά, η ``Python`` δε θα μπορούσε να μην έχει την αντίστοιχη συνάρτηση για να μας διευκολύνει, και αυτή είναι η ``residue`` της $\\texttt{SciPy}$ βιβλιοθήκης. \n",
    "\n",
    "Έστω λοιπόν ότι έχουμε μια ρητή συνάρτηση του $j2\\pi f$ ή του $s$. Η σύνταξή της εντολής για να μας κάνει το ανάπτυγμα είναι η ακόλουθη:\n",
    "\n",
    "$$\\textrm{r, p, k = residue(b, a)}$$\n",
    "\n",
    "με \n",
    "\n",
    "- $\\textrm{b}$ να είναι ένα διάνυσμα που περιέχει τους συντελεστές του πολυωνύμου (του $s$ ή $j2\\pi f$) του αριθμητή σε φθίνουσα σειρά βαθμού (πρώτα ο συντελεστής του όρου με το μεγαλύτερο βαθμό, κλπ, κλπ, και τελευταία η σταθερά), και \n",
    "- $\\textrm{a}$ να είναι ένα διάνυσμα που περιέχει τους συντελεστές του πολυωνύμου (του $s$ ή $j2\\pi f$) του παρονομαστή, επίσης σε φθίνουσα σειρά βαθμού (πρώτα ο συντελεστής του όρου με το μεγαλύτερο βαθμό, κλπ, κλπ, και τελευταία η σταθερά).\n",
    "\n",
    "Η συνάρτηση επιστρέφει τρία διανύσματα: \n",
    "- το διάνυσμα $\\textrm{r}$ περιέχει τους αριθμητές των μερικών κλασμάτων που προκύπτουν, \n",
    "- το διάνυσμα $\\textrm{p}$ περιέχει τους πόλους των απλών κλασμάτων, δηλ. τις ρίζες των παρονομαστών των μερικών κλασμάτων, και \n",
    "- το διάνυσμα $\\textrm{k}$ περιέχει συντελεστές όρων που προκύπτουν κατά τη διαίρεση πολυωνύμων στην περίπτωση που η ρητή συνάρτηση έχει βαθμό πολυωνύμου αριθμητή $M$ μεγαλύτερο ή ίσο με τον αντίστοιχο του παρονομαστή, $N$, δηλ. $M \\geq N$.\n",
    "\n",
    "Ας δούμε ένα παράδειγμα. Έστω:\n",
    "\n",
    "$$X(s) = \\frac{s^2 + 2s + 1}{s^3 -3s +2}$$\n",
    "\n",
    "Ο αριθμητής έχει συντελεστές πολυωνύμου τους $[1, 2, 1]$ και ο παρονομαστής $[1, 0, -3, 2]$ αντίστοιχα, με το $0$ στη δεύτερη θέση του διανύσματος να βρίσκεται εκεί γιατί ο δευτεροβάθμιος όρος δεν υπάρχει στο πολυώνυμο του παρονομαστή. Έτσι, ο αντίστοιχος συντελεστής του είναι $0$. \n",
    "\n",
    "Εκτελούμε λοιπόν:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "b = [1, 2, 1]\n",
    "a = [1, 0, -3, 2]\n",
    "\n",
    "r, p, k = residue(b, a)\n",
    "\n",
    "print(f\"r = {r} \\np = {p} \\nk = {k}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Από αυτό καταλαβαίνουμε αμέσως ότι επειδή το διάνυσμα $k$ είναι άδειο, δεν υπάρχουν όροι από τη διαίρεση των πολυωνύμων - λογικό κι αναμενόμενο, καθώς ο αριθμητής μας είναι δευτεροβάθμιος και ο παρονομαστής τριτοβάθμιος - δεν υπήρχε λόγος να κάνουμε διαίρεση των πολυωνύμων. Από τα υπόλοιπα, παρατηρούμε ότι μπορούμε να γράψουμε το αποτέλεσμα του αναπτύγματος ως\n",
    "\n",
    "$$X(s) = \\frac{s^2 + 2s + 1}{s^3 - 3s + 2} = \\frac{1/9}{s+2} + \\frac{4/3}{s-1} + \\frac{8/9}{s-1}  = \\frac{1}{9}\\left(\\frac{1}{s+2}\\right) + \\frac{4}{3}\\left(\\frac{1}{s-1}\\right) + \\frac{8}{9}\\left(\\frac{1}{s-1}\\right)$$\n",
    "\n",
    "Δυστυχώς, η ``Python`` μας αφήνει εδώ. 😃 Στη συνέχεια, ανάλογα με το πεδίο σύγκλισης, μπορούμε εμείς να γυρίσουμε τα σήματά μας πίσω στο χρόνο μέσω των γνωστών πινάκων μετασχηματισμού Laplace. \n",
    "\n",
    "---\n",
    "\n",
    "Αναπτύξτε σε μερικά κλάσματα τους παρακάτω μετασχηματισμούς:\n",
    "\n",
    "(α') $$X(s) = \\frac{24}{s(s+8)}$$\n",
    "(β') $$X(s) = \\frac{20}{s^2+4s+3}$$\n",
    "(γ') $$X(s) = \\frac{s^2}{s^2-4s+4}$$\n",
    "(δ') $$X(s) = \\frac{s^4 + 2s^3 + 5s + 10}{s^2 - 1}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Απαντήσεις:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (α')\n",
    "# INSERT CODE HERE\n",
    "# Πρέπει να βρείτε\n",
    "# r = [ 3. -3.] , p = [ 0. -8.] , k = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (β')\n",
    "# INSERT CODE HERE\n",
    "# Πρέπει να βρείτε\n",
    "# r = [ 10. -10.] , p = [-1. -3.] , k = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (γ')\n",
    "# INSERT CODE HERE\n",
    "# Πρέπει να βρείτε\n",
    "# r = [4. 4.] , p = [2. 2.] , k = [1.]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (δ')\n",
    "# INSERT CODE HERE\n",
    "# Πρέπει να βρείτε\n",
    "# r = [-2.  9.] , p = [-1.  1.] , k = [1. 2. 1.]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2η Άσκηση\n",
    "\n",
    "Εργάζεστε σε μια από τις πρώτες εταιρίες κινητής τηλεφωνίας, και το πόστο σας είναι $\\textrm{μηχανικός σχεδίασης ϕίλτρων}$.\n",
    "\n",
    "Ο προϊστάμενός σας συγκαλεί σύσκεψη στην οποία αποφασίζεται ότι εσείς πρέπει να σχεδιάσετε και να αναπτύξετε ένα ΓΧΑ σύστημα με απόκριση συχνότητας $H(f)$ για εφαρμογές επικοινωνίας ϕωνής, το οποίο θα αποκόπτει τις συχνότητες μεγαλύτερες από κάποια δοθείσα $f_c$ (η οποία λέγεται, όπως ξέρετε, συχνότητα αποκοπής - cut-off frequency) ενώ ϑα κρατά όσο γίνεται ανέπαφες τις συχνότητες μικρότερες από $f_c$. Τέτοια συστήματα γνωρίζετε ότι ονομάζονται $\\textbf{ϕίλτρα επιλογής συχνοτήτων}$, και για αυτήν την άσκηση θα τα αποκαλούμε απλά ως ϕίλτρα.\n",
    "\n",
    "Ο προϊστάμενός σας, που δε γνωρίζει θεωρία σημάτων και συστημάτων (ενώ εσείς πλέον είστε ξεφτέρια! 😃), σας παραδίδει την απόκριση συχνότητας $H(f)$ που θέλει να ϕτιάξετε, στο παρακάτω σχήμα, και σας αναφέρει ότι το ζητούμενο $f_c$ ισούται με $f_c = 3000$ Hz, αφού το ϕίλτρο θα ενσωματωθεί σε στρατιωτικά ασύρματα τηλεφωνικά συστήματα, όπου το εύρος ζώνης επικοινωνίας - και το κόστος λειτουργίας - είναι περιορισμένο, καθώς σημασία έχει να είναι απλά καταληπτή η επικοινωνία και όχι η ποιότητά της."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"LP.png\" alt=\"Lowpass\" width=\"600\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(α')** Γράψτε στο παρακάτω κελί (κάνοντας double click πάνω του) την εξήγηση που θα δώσετε στον προϊστάμενό σας σχετικά με το αδύνατο της πραγματοποίησης του παραπάνω φίλτρου!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Απάντηση:\n",
    "\n",
    "(γράψτε εδώ)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(β')** Αφού τον πείσατε για την ορθότητα του παραπάνω ερωτήματος, σας αναθέτει να υλοποιήσετε ένα φίλτρο που να πλησιάζει όσο γίνεται αυτό που σας ζήτησε αρχικά, και να είναι υλοποιήσιμο. \n",
    "Στην προσπάθειά σας αυτή, αναζητείτε μια συνάρτηση που να ταιριάζει στις \"προδιαγραφές\" σας. Ένας μαθηματικός φίλος σας αναφέρει ότι έχει υπόψη του μια συνάρτηση η οποία να πλησιάζει το ζητούμενο φίλτρο σας, και την οποία σχεδιάζει πρόχειρα στο χαρτί, όπως στο παρακάτω σχήμα. \n",
    "\n",
    "<img src=\"butter2.png\" alt=\"Butterworth\" width=\"600\"/>\n",
    "\n",
    "Η συνάρτηση ονομάζεται\n",
    "\n",
    "### $$\\text{συνάρτηση Butterworth}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Μη έχοντας καλύτερη εναλλακτική, του ζητάτε να σας δώσει τη μαθηματική περιγραφή της συνάρτησης. Σας δίνει μια περιγραφή στο χώρο της συχνότητας που βρήκε σε κάποιο μαθηματικό εγχειρίδιο, ως\n",
    "\n",
    "$$\\large |H(f)|^2 = \\displaystyle \\frac{1}{1 + \\left(\\frac{j2\\pi f}{j2\\pi f_c}\\right)^{2N}}$$\n",
    "\n",
    "με $N$ την **τάξη** της συνάρτησης, όπως σας ανέφερε. \n",
    "\n",
    "Προσέξτε ότι σας δίνει μια περίγραφή του **μέτρου** του φίλτρου **στο τετράγωνο** (δεν ήταν και ιδιαίτερα βοηθητικός, μάλλον βαριόταν 😄)...\n",
    "\n",
    "$\\rightarrow$ Μετατρέψτε τη συνάρτηση αυτή στο χώρο του μετασχ. Laplace (βρείτε δηλ. το $|H(s)|^2$) θέτοντας $s = j2\\pi f$. \n",
    "\n",
    "**Γράψτε την απάντησή σας παρακάτω:**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Απάντηση:\n",
    "\n",
    "$$\\large |H(s)|^2 = ?$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(γ')** Θέλετε να μελετήσετε τη συμπεριφορά του φίλτρου - όπως το ονομάζετε πλέον - Butterworth, για να το κατανοήσετε καλύτερα. \n",
    "\n",
    "$\\rightarrow$ Βρείτε τους πόλους του $|H(s)|^2$ στο $s$-επίπεδο.\n",
    "\n",
    "**Γράψτε την απάντησή σας παρακάτω:**\n",
    "\n",
    "$\\underline{Βοήθεια:}$\n",
    "Πολλαπλασιάστε την απάντησή σας στο προηγ. ερώτημα σε αριθμητή και παρονομαστή με το $(j2\\pi f_c)^{2N}$ και βρείτε τις ρίζες του παρονομαστή, χρησιμοποιώντας τη μεθοδολογία του De Moivre από την πρώτη κιόλας διάλεξη του μαθήματος. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Απάντηση: \n",
    "\n",
    "Οι πόλοι $s_k$ είναι οι $$\\large s_k = ?$$\n",
    "\n",
    "για $k=0,1,\\cdots,2N-1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(δ')** Γνωρίζετε από τη θεωρία σας ότι επειδή το φίλτρο σας είναι πραγματικό σήμα στο χρόνο (δηλ. $h(t) \\in \\Re$), θα ισχύουν οι γνωστές συμμετρίες στο χώρο του Fourier, δηλ.\n",
    "\n",
    "$$H(-f) = H^*(f)$$\n",
    "\n",
    "Με χρήση αυτής, μπορούμε να γράψουμε\n",
    "\n",
    "$$|H(f)|^2 = H(f)H^*(f) = H(f)H(-f) = H(s)H(-s)\\Big|_{s=j2\\pi f}$$\n",
    "\n",
    "Άρα οι πόλοι που βρήκατε νωρίτερα στο χαρτί σας ανήκουν στο σύστημα $$|H(s)|^2 = H(s)H(-s)$$ \n",
    "\n",
    "Καταλαβαίνετε ότι κάποιοι πόλοι από αυτούς θα ανήκουν στο $H(s)$ και οι υπόλοιποι στο $H(-s)$, δεδομένου ότι έχουμε γινόμενο συναρτήσεων μεταφοράς. Το σύνολο των πόλων που έχετε βρει θα είναι η τομή του συνόλου των πόλων των δυο συστημάτων, $H(s)$ και $H(-s)$.\n",
    "\n",
    "Έχοντας βρει τους πόλους από το προηγούμενο ερώτημα, σχεδιάστε το διάγραμμα πόλων-μηδενικών του $|H(s)|^2$ σε ένα χαρτί σας (χωρίς να ασχοληθείτε με τα μηδενικά) κι επιλέξτε από τους πόλους που σχεδιάσατε ένα υποσύνολο πόλων που θα ανήκει στο $H(s)$, ώστε το σύστημα που θα προκύψει από αυτούς να είναι $\\textbf{ευσταθές και αιτιατό}$. Προσέξτε ότι αν $s_p$ είναι ένας πόλος του $H(s)$, τότε ο $-s_p$ είναι $\\textbf{αναγκαστικά}$ πόλος του $H(-s)$. Επίσης, ένας πόλος $s_p$ θα ανήκει είτε στο ένα σύστημα είτε στο άλλο, και όχι και στα δυο ταυτόχρονα. Τέλος, θυμηθείτε ότι $h(t) \\in \\Re$, οπότε πρέπει να γνωρίζετε κάτι ακόμα για τους πόλους τέτοιων σημάτων, και τον τρόπο που εμφανίζονται στα διαγράμματα (ως ζεύγος...... ?).\n",
    "\n",
    "$\\rightarrow$ Ποιούς πόλους επιλέξατε από το $|H(s)|^2$ για το $H(s)$, και γιατί? Απαντήσετε παρακάτω:\n",
    "\n",
    "**(μην περιγράψετε αριθμητικά τους πόλους - περιγράψτε τους ως προς τη θέση τους στο μιγαδικό επίπεδο)**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Απάντηση:\n",
    "\n",
    "(γράψτε εδώ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(ε')** Επιστρέφοντας στη θεωρία του φίλτρου Butterworth... θα χρειαστούμε τους πόλους που βρήκατε στο ερώτημα **(γ')**. Παρατηρήστε ότι\n",
    "\n",
    "$$|H(s)|^2 \\Big|_{s=0} = H(s)H(-s)\\Big|_{s=0} = 1$$\n",
    "\n",
    "και ότι\n",
    "\n",
    "$$|H(s)|^2 = H(s)H(-s) = \\frac{(j2\\pi f_c)^{2N}}{\\prod_{k=0}^{N-1}(s+s_k)\\prod_{l=0}^{N-1}(s+s_l)}$$\n",
    "\n",
    "με τον παρονομαστή να αφορά τους πόλους του $H(s)$ και του $H(-s)$, αντίστοιχα.\n",
    "\n",
    "$\\rightarrow$ Υπολογίστε σε ένα πρόχειρο χαρτί το $H(s)$ για $N=1$ και $N=2$, βρίσκοντας το $|H(s)|^2$ από την παραπάνω σχέση και διαχωρίζοντας τα $H(s)$, $H(-s)$ με βάση την απαίτηση το $H(s)$ να είναι $\\textbf{ευσταθές και αιτιατό}$.\n",
    "\n",
    "$\\underline{Βοήθεια:}$\n",
    "Πρέπει να βρείτε: $$\\displaystyle H_{N=1}(s) = \\frac{2\\pi f_c}{s + 2\\pi f_c} , \\: \\: H_{N=2}(s) = \\frac{(2\\pi f_c)^2}{(s + 2\\pi f_c e^{j3\\pi/4})(s + 2\\pi f_c e^{j5\\pi/4})}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(στ)** Εν τέλει, πώς θα υλοποιούσαμε το φίλτρο Butterworth στην πραγματικότητα? Ένας τρόπος είναι μέσω μιας διαφορικής εξίσωσης που μοντελοποιεί ένα ηλεκτρικό κύκλωμα! \n",
    "\n",
    "$\\rightarrow$ Βρείτε τη διαφορική εξίσωση $\\textit{δευτέρας τάξης}$ (δηλ. για $N=2$) που αντιστοιχεί στο φίλτρο Butterworth που βρήκατε παραπάνω, χρησιμοποιώντας συχνότητα αποκοπής $f_c = \\frac{1}{2\\pi}$ Hz, και γράψτε τους συντελεστές κάθε όρου της παρακάτω:\n",
    "\n",
    "### Απάντηση:\n",
    "1. Συντελεστής του $y(t) = ? $\n",
    "\n",
    "2. Συντελεστής του $\\frac{d}{dt}y(t) = ? $\n",
    "\n",
    "3. Συντελεστής του $\\frac{d^2}{dt^2}y(t) = ?$\n",
    "\n",
    "4. Συντελεστής του $x(t) = ?$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(ζ)** Υλοποιήστε στην ``Python`` την απόκριση πλάτους $|H(f)|$ του φίλτρου σας για $f_c = 3000$ Hz, δειγματοληπτώντας έναν άξονα συχνοτήτων $f$ στο διάστημα $[-5000, 5000]$ ανά $Df = 0.1$ Hz, για \n",
    "- $N=6$, \n",
    "- $N=16$,\n",
    "- $N=46$. \n",
    "\n",
    "Για αυτό το ερώτημα, συμπληρώστε τον παρακάτω κώδικα:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = np.arange(start=-5000, stop=5000, step=0.1)           # Άξονας συχνοτήτων\n",
    "fc = 3000                                                 # Συχνότητα αποκοπής\n",
    "\n",
    "N = 6\n",
    "H1 = # INSERT CODE HERE - υλοποιήστε το φίλτρο |H(f)| για N=6\n",
    "N = 16\n",
    "H2 = # INSERT CODE HERE - υλοποιήστε το φίλτρο |H(f)| για N=16\n",
    "N = 46\n",
    "H3 = # INSERT CODE HERE - υλοποιήστε το φίλτρο |H(f)| για N=46\n",
    "\n",
    "plt.figure(figsize=(8,4))    # Γράφημα\n",
    "\n",
    "plt.plot(f, np.abs(H1))       # Το πρώτο φίλτρο (βάζουμε np.abs γιατί προκύπτουν πολύ πολύ μικρά φανταστικά μέρη)\n",
    "plt.plot(f, np.abs(H2))       # Το δεύτερο\n",
    "plt.plot(f, np.abs(H3))       # και το τρίτο\n",
    "\n",
    "plt.title('Απόκριση πλάτους')            # Ομορφαίνουμε\n",
    "plt.xlabel('Συχνότητα (Hz)')             # Ομορφαίνουμε\n",
    "plt.grid(True)                           # Ομορφαίνουμε\n",
    "plt.legend(['N=6', 'N=16', 'N=46'])      # Ομορφαίνουμε\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(η')** Τι παρατηρείτε στο φίλτρο σας όσο αυξάνεται η τάξη του, $N$? Γράψτε παρακάτω την απάντησή σας:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Απάντηση:\n",
    "\n",
    "(γράψτε εδώ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(θ)** Προτού παραδώσετε το φίλτρο σας στον προϊστάμενό σας ώστε να υλοποιηθεί σε κύκλωμα, θέλετε να βεβαιωθείτε ότι λειτουργεί όπως πρέπει, εξομοιώνοντάς το στην Python και βάζοντας ως είσοδο μια τυπική στρατιωτική διαταγή, δωρεά του Υπουργείου Άμυνας. Θα τη βρείτε στο αρχείο ``military.wav``, στο .RAR file του εργαστηρίου.\n",
    "\n",
    "Ακολουθήστε τα παρακάτω βήματα."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(θ.1')** Φορτώστε το αρχείο στην Python με τη - γνωστή πια - συνάρτηση $\\texttt{wav.read}$ της $\\texttt{SciPy}$ και δείτε το. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs, s = wav.read('./military.wav')                   # Ανάγνωση\n",
    "s = s/np.abs(np.max(s))                              # Κανονικοποίηση στο [-1, 1]\n",
    "t = np.arange(0, len(s)/fs, 1/fs)                    # Άξονας χρόνου\n",
    "\n",
    "plt.figure(figsize=(12,4))                           # Γράφημα\n",
    "plt.plot(t, s)                                       # Περιεχόμενο \n",
    "plt.xlabel('Χρόνος (s)')                             # Ομορφαίνουμε\n",
    "plt.ylabel('Πλάτος')                                 # Ομορφαίνουμε\n",
    "plt.grid()                                           # Ομορφαίνουμε\n",
    "plt.title('Σήμα στο χρόνο')                          # Ομορφαίνουμε"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ας το ακούσουμε:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(s, rate=fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(θ.2')** Η έτοιμη συνάρτηση ``butter`` υλοποιεί διάφορα φίλτρα Butterworth με τάξη $N$ την οποία παρέχετε εσείς ως όρισμα, όπως και τη συχνότητα αποκοπής $f_c$, και επιστρέφει τους συντελεστές του αριθμητή και του παρονομαστή για το φίλτρο σας στο χώρο του Laplace (ή και του Fourier)! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = 3000           # Συχνότητα αποκοπής\n",
    "N = 30              # Ας επιλέξουμε τάξη ίση με 30 - ένα \"αυστηρό\" φίλτρο Butterworth με απότομη μετάβαση\n",
    "b, a = butter(N=N, btype='low', Wn=fc, analog=False, output='ba', fs = fs)   # Δημιουργία φίλτρου"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Παρατηρήστε ότι ζητήσαμε από τη συνάρτηση ``butter`` να μας φτιάξει ένα lowpass (χαμηλοπερατό) φίλτρο, τάξης $N$, με συχνότητα απόκοπής $W_n = f_c$, το φίλτρο μας να είναι $\\textbf{ψηφιακό}$ κατευθείαν - θέτοντας $\\texttt{analog=False}$ - και να μας δώσει ως έξοδο τους συντελεστές $b, a$ της συνάρτησης μεταφοράς του φίλτρου μας (που όμως, ως ψηφιακό φίλτρο, δεν ορίζεται στο χώρο του Laplace)! Το διάνυσμα των συντελεστών του αριθμητή θα είναι το $b$ ενώ του παρονομαστή θα είναι το $a$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(θ.3')** Δείτε την απόκριση συχνότητας $H(f)$ (μόνο για $f > 0$) του (ψηφιακού πλέον) φίλτρου σας με χρήση των εντολών:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, H = freqz(b=b, a=a, fs=fs)          # Καλούμε τη συνάρτηση freqz για να μας βρει την απόκριση σε συχνότητα H(f)\n",
    "\n",
    "plt.figure(figsize=(10,6))             # Γράφημα\n",
    "\n",
    "plt.subplot(2, 1, 1)                   # 2 γραμμές, 1 στήλη στο γράφημα: 1η γραμμή\n",
    "plt.plot(f, np.abs(H))                 # Εμφάνιση απόκρισης πλάτους\n",
    "plt.title('Απόκριση πλάτους')          # Ομορφαίνουμε      \n",
    "plt.xlabel('Συχνότητα (Hz)')           # Ομορφαίνουμε\n",
    "plt.grid(True)                         # Ομορφαίνουμε\n",
    "\n",
    "plt.subplot(2, 1, 2)                   # 2 γραμμές, 1 στήλη στο γράφημα: 2η γραμμή\n",
    "plt.plot(f, np.angle(H))               # Εμφάνιση απόκρισης φάσης\n",
    "plt.title('Απόκριση Φάσης')            # Ομορφαίνουμε\n",
    "plt.xlabel('Συχνότητα (Hz)')           # Ομορφαίνουμε\n",
    "plt.grid(True)                         # Ομορφαίνουμε\n",
    "\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**(θ.4')** Ας χρησιμοποιήσουμε τη συνάρτηση ``lfilter`` της βιβλιοθήκης ``SciPy`` για την υλοποίηση του συστήματος - είναι παρόμοια με την ``np.convolve`` που έχετε δει σε άλλο Εργαστήριο, απλά παίρνει διαφορετικά ορίσματα."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = lfilter(b, a, x=s)       # Φιλτράρισμα (ουσιαστικά, συνέλιξη στο χρόνο του φίλτρου με το σήμα εισόδου)\n",
    "y = y/np.abs(np.max(y))      # Κανονικοποίηση (δεν είναι κι απαραίτητο :) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "για να φιλτράρουμε το σήμα μέσα από το σύστημα που φτιάξαμε! Ας το ακούσουμε!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(y, rate=fs) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Συγκρίνετε με το αρχικό:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(s, rate=fs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Το ίδιο και οπτικά!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,6))                      # Γράφημα\n",
    "\n",
    "plt.subplot(2, 1, 1)                            # 2 γραμμές, 1 στήλη στο γράφημα: 1η γραμμή\n",
    "plt.plot(t, s)                                  # Περιεχόμενο \n",
    "plt.grid()\n",
    "plt.title('Αρχικό σήμα $s(t)$')                 # Ομορφαίνουμε\n",
    "plt.xlabel('Χρόνος (s)')                        # Ομορφαίνουμε\n",
    "plt.ylabel('Πλάτος')                            # Ομορφαίνουμε\n",
    "\n",
    "plt.subplot(2, 1, 2)                            # 2 γραμμές, 1 στήλη στο γράφημα: 2η γραμμή\n",
    "plt.plot(t, y)                                  # Περιεχόμενο \n",
    "plt.grid()\n",
    "plt.title('Τελικό σήμα $y(t)$')                 # Ομορφαίνουμε\n",
    "plt.xlabel('Χρόνος (s)')                        # Ομορφαίνουμε\n",
    "plt.ylabel('Πλάτος')                            # Ομορφαίνουμε\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()\n",
    "plt.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Από τα σήματα στο χρόνο, ίσως να δυσκολεύεστε να καταλάβετε τη διαφορά (θα πρέπει να κάνουμε zoom σε κάποια μέρη για να δούμε περισσότερη λεπτομέρεια). \n",
    "\n",
    "Ας δούμε τα δυο σήματα στη συχνότητα, εκτελώντας το Μετασχηματισμό Fourier για καθένα σήμα. Επειδή όμως τα σήματα ομιλίας, μουσικής, και γενικότερα τα πραγματικά ηχητικά σήματα, αλλάζουν πολύ γρήγορα όσο περνάει ο χρόνος, προτιμάμε να εκτελούμε το μετασχ. Fourier σε **τμήματα** (segments ή frames) σήματος, όπως το έχετε ήδη δει σε προηγούμενες σειρές Εργαστηρίων. \n",
    "\n",
    "Αντί όμως να τυπώσουμε ένα σωρό μετασχ. Fourier (καθώς τα σήματά μας περιέχουν πάρα πολλά τμήματα), καλύτερα να:\n",
    "1. Χωρίσουμε το σήμα στα προαναφερθέντα τμήματα\n",
    "2. Εκτελέσουμε το μετασχ. Fourier στο καθένα\n",
    "3. Υπολογίσουμε το μέτρο του (δηλ. το φάσμα πλάτους)\n",
    "4. Μαζέψουμε το αποτέλεσμα σε μια μεγάλη δομή (διδιάστατο πίνακα, π.χ.)\n",
    "5. Τυπώσουμε όλα μαζί τα μέτρα των μετασχ. Fourier ως μια εικόνα, με άξονες $f$ και $t$, και το χρώμα κάθε pixel στις συντεταγμένες $(t_0, f_0)$ θα μας υποδηλώνει πόσο \"δυνατό\" είναι το μέτρο του μετασχηματισμού σε εκείνο το σημείο.\n",
    "\n",
    "Το τελικό αποτέλεσμα ονομάζεται **φασματογράφημα (spectrogram)**. \n",
    "\n",
    "**ΕΥΤΥΧΩΣ** όλα αυτά τα κάνει η συνάρτηση ``specgram`` της βιβλιοθήκης απεικονίσεων ``matplotlib`` (μεταξύ πολλών άλλων βιβλιοθηκών που έχουν παρόμοιες συναρτήσεις), και σας δίνουμε όλο το σχετικό κώδικα έτοιμο. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,8))\n",
    "\n",
    "plt.subplot(2, 2, 1)\n",
    "plt.title('Αρχικό σήμα')\n",
    "plt.plot(t, s)\n",
    "plt.grid()\n",
    "plt.ylabel('Πλάτος')\n",
    "\n",
    "plt.subplot(2, 2, 2)\n",
    "plt.specgram(s, Fs=fs)\n",
    "plt.title('Φασματογράφημα αρχικού σήματος')\n",
    "plt.ylabel('Συχνότητα (Hz)')\n",
    "\n",
    "plt.subplot(2, 2, 3)\n",
    "plt.title('Φιλτραρισμένο σήμα')\n",
    "plt.plot(t, y)\n",
    "plt.grid()\n",
    "plt.xlabel('Χρόνος (s)')\n",
    "plt.ylabel('Πλάτος')\n",
    "\n",
    "plt.subplot(2, 2, 4)\n",
    "plt.specgram(y, Fs=fs)\n",
    "plt.title('Φασματογράφημα φιλτραρισμένου σήματος')\n",
    "plt.xlabel('Χρόνος (s)')\n",
    "plt.ylabel('Συχνότητα (Hz)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Τι κερδίσαμε με το παραπάνω φιλτράρισμα? Το σήμα μας πλέον έχει **εύρος ζώνης (bandwidth)** περίπου ίσο με $3000$ Hz! (αντί του αρχικού, που ήταν $8000$ Hz). Άρα για να το μεταδώσουμε μέσα από ένα κανάλι χρειαζόμαστε το $1/4$ του αρχικού εύρους ζώνης. Αυτό σημαίνει ότι πληρώνουμε μικρότερο εύρος συχνοτήτων, άρα εξοικονομούμε χρήματα, χωρίς να αλλάζει το μήνυμα που φέρει το σήμα (πιστεύω ότι καταλαβαίνετε τις λέξεις στα σήματα και στις δυο περιπτώσεις)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$\\rightarrow$ Με βάση τόσο το ακουστικό αποτέλεσμα, όσο και τα παραπάνω σχήματα, **πώς θα χαρακτηρίζατε την ποιότητα του σήματος εξόδου σε σχέση με το αρχικό**; Σκεφτείτε πως επέδρασε το φίλτρο επάνω στις συχνότητες του σήματος εισόδου. Απαντήστε παρακάτω:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Απάντηση:\n",
    "\n",
    "(γράψτε εδώ)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3η Άσκηση\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Θυμηθείτε το 4ο εργαστήριο, την άσκηση με την αποθορυβοποίηση του μουσικού κομματιού \"Ρόζα\", των Αλκαίου/Μικρούτσικου, όπου ήταν μολυσμένο με ένα ισχυρό ημίτονο συγκεκριμένης συχνότητας που κληθήκατε να βρείτε κάνοντας συνεχείς μετασχηματισμούς Fourier πάνω στο σήμα - και με τα hints ότι το ημίτονο είχε $A = 0.01$, $\\phi = 0$, και η συχνότητα που ψάχνατε ήταν ακέραιο πολλαπλάσιο των $100$ Hz στο διάστημα $[2000, 3000]$ Hz 😃😃 !! \n",
    "\n",
    "Όπως είχαμε επισημάνει τότε, το παράδειγμα εκείνο ήταν πολύ εκπαιδευτικό, καθώς σας δινόταν το πλάτος, η φάση, και το \"εύρος\" που θα έπρεπε να ψάξετε την άγνωστη συχνότητα!\n",
    "\n",
    "Στην πράξη, δεν μπορεί κανείς να ξέρει το πλάτος και τη φάση του ημιτονοειδούς \"θορύβου\", ώστε να εκτελέσει τόσο απλά την αφαίρεση $signal(t) -noise(t)$ που κάναμε τότε στο 4ο εργαστήριο, αφού είχαμε κατασκευάσει ένα \"αντίγραφο\" του ημιτονοειδούς θορύβου. Έτσι, χρησιμοποιούμε τα λεγόμενα **notch** φίλτρα, που αναλαμβάνουν να αποκόψουν μια πολύ πολύ μικρή μπάντα συχνοτήτων, ανεξαρτήτως πλάτους και φάσης. Τα φίλτρα αυτά ουσιαστικά μηδενίζουν το πλάτος μιας ανεπιθύμητης συχνότητας $f_0$, με αποτέλεσμα η φάση να μη μας ενδιαφέρει - το $0\\cdot \\cos(2\\pi f_0 t + \\phi)$ ισούται με $0$, ανεξαρτήτως του $\\phi$. Η απόκριση πλάτους ενός notch φίλτρου με συχνότητα αποκοπής $f_c$ φαίνεται στο παρακάτω σχήμα.\n",
    "\n",
    "![notch](notch.png \"Notch Filer\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Η σχεδίαση τέτοιων συστημάτων γίνεται κατά κανόνα στο χώρο του μετασχ. Laplace. \n",
    "\n",
    "Η θεωρία της κατασκευής της συνάρτησης μεταφοράς $H(s)$ ενός συστήματος notch φίλτρου είναι απλή: για παράδειγμα, αν θέλουμε να μηδενίσουμε τη συχνότητα $f_c = 5000$ Hz του σήματος που θα δεχθεί ως είσοδο το σύστημά μας, τότε:\n",
    "- Μετατρέπουμε την $f_c$ σε $\\omega_c = 2\\pi f_c$ για ευκολία, δηλ. $\\omega_c = 10000\\pi$ rad/s\n",
    "- Τοποθετούμε δυο συζυγή μηδενικά **ακριβώς πάνω** στις συχνότητες $\\pm 10000\\pi$, δημιουργώντας έτσι τον αριθμητή $$P(s) = (s-j\\omega_c)(s+j\\omega_c) = s^2 + \\omega_c^2$$\n",
    "- Τοποθετούμε δυο συζυγείς πόλους στις ίδιες συχνότητες ΑΛΛΑ όχι επάνω στο φανταστικό άξονα: ΚΟΝΤΑ του όμως, από τη μεριά του αριστερού ημιεπιπέδου (γιατί?), δημιουργώντας έτσι τον παρονομαστή $$Q(s) = (s - s_1)(s - s_1^*) = s^2 + \\frac{\\omega_c}{Q}s + \\omega_c^2$$\n",
    "- Οι πόλοι ειναι της μορφής $$s_{1,1^*} = -\\frac{1}{2Q}\\omega_c \\pm j\\omega_c\\sqrt{1 - \\left(\\frac{1}{2Q}\\right)^2}$$\n",
    "Ο λόγος που οι πόλοι είναι γραμμένοι έτσι δεν είναι του παρόντος 😄 αλλά η μορφή τους μας επιτρέπει, \"παίζοντας\" με την τιμή του $Q$, να μπορούμε να ελέγχουμε την \"αυστηρότητα\" του φίλτρου στη συχνότητα αποκοπής.\n",
    "- Έτσι, η συνάρτηση μεταφοράς μας θα είναι η \n",
    "$$H(s) = \\frac{P(s)}{Q(s)} = \\frac{s^2 + \\omega_c^2}{s^2 + \\frac{\\omega_c}{Q}s + \\omega_c^2}$$\n",
    "\n",
    "Ας γράψουμε μια συνάρτηση που υλοποιεί ένα τέτοιο φίλτρο. Συμπληρώστε τον κώδικα."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_laplace_notch_filter(fc, Q):\n",
    "    \"\"\"\n",
    "    Δημιουργία ενός notch φίλτρου και υλοποίησή του στο χώρο του Laplace.\n",
    "    \n",
    "    Παράμετροι:\n",
    "    fc : Συχνότητα που πρέπει να κοπεί\n",
    "    Q : Συντελεστής ποιότητας. Όσο ψηλότερος, τόσο αυστηρότερο - και ακριβότερο - το φίλτρο στην αποκοπή\n",
    "        \n",
    "    Έξοδος:\n",
    "    numL, denL : Συντελεστές για αριθμητή και παρονομαστή του notch φίλτρου H(s)\n",
    "    \"\"\"\n",
    "    # Συνήθης μορφή του φίλτρου notch: H(s) = (s^2 + ωc^2) / (s^2 + (ωc/Q)s + ωc^2)\n",
    "    \n",
    "    # Γωνιακή συχνότητα ωc = 2πfc rad/s\n",
    "    omegac = 2 * np.pi * fc\n",
    "    \n",
    "    # H(s) : αριθμητής : συντελεστές του πολυωνύμου του s\n",
    "    num_laplace = # INSERT CODE HERE\n",
    "    \n",
    "    # H(s) : παρονομαστής : συντελεστές του πολυωνύμου του s\n",
    "    den_laplace = # INSERT CODE HERE\n",
    "    \n",
    "    return num_laplace, den_laplace"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Πριν συνεχίσουμε, ας δούμε γιατί η απλή τοποθέτηση δυο μηδενικών επάνω στο φανταστικό άξονα, στις συχνότητες $\\pm j\\omega_c$ ΔΕΝ είναι αρκετή, όπως ίσως θα αναμένατε... \n",
    "\n",
    "Ας καλέσουμε τη συνάρτηση που φτιάξαμε, αλλά ας αγνοήσουμε τον παρονομαστή (θεωρώντας τον σταθερό πολυώνυμο, μοναδιαίο), έτσι ώστε να έχει πόλους στο άπειρο και μηδενικά στις θέσεις που θέλουμε. \n",
    "\n",
    "Έστω η συχνότητα $f_c = 1560$ Hz που θέλουμε να κόψουμε."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = 1560 # Hz\n",
    "Q = 2\n",
    "numL, denL = create_laplace_notch_filter(fc, Q) # το κανονικό μας φίλτρο, αμέσως μετά θα αγνοήσουμε τον denL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Σας δίνεται έτοιμη συνάρτηση ``plot_laplace_frequency_response`` που αναλαμβάνει να απεικονίσει την απόκριση πλάτους μέσω της συνάρτησης μεταφοράς που φτιάξατε.\n",
    "\n",
    "Γράψτε ``help(plot_laplace_frequency_response)`` για να δείτε τη σύνταξή της και μερικές λεπτομέρειες. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(plot_laplace_frequency_response)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Ας την καλέσουμε."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Βάζουμε [1] ως παρονομαστή για να δηλώσουμε ότι το H(s) είναι πρακτικά \n",
    "# μόνο της μορφής H(s) = P(s)/1 = (s^2 + ωc^2)/1 = s^2 + ωc^2\n",
    "plot_laplace_frequency_response(numL, [1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Χμμ... πράγματι μηδενίζεται η ζητούμενη συχνότητα των $f_c = 1560$ Hz αλλά η συμπεριφορά του φίλτρου μας δεν είναι καλή, καθώς η απόκριση πλάτους φεύγει στο άπειρο όσο μεγαλώνουν οι συχνότητες, ενισχύοντάς τες πολύ. Ταυτόχρονα, οι συχνότητες στο $[0, f_c)$ έχουν πολύ χαμηλότερο πλάτος, συγκριτικά με τις υπόλοιπες. Άρα μόνο η τοποθέτηση δυο μηδενικών στις θέσεις $s = j\\omega_c = j2\\pi f_c$ ΔΕΝ αρκεί, πράγματι. \n",
    "\n",
    "Αν όμως κρατήσουμε και τους πόλους? Τι γίνεται? Ας το δούμε."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_laplace_frequency_response(numL, denL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ωωωωω! 😍😍 Τώρα το φίλτρο μας δουλεύει πολύ καλύτερα!! Η συχνότητα $f_c = 1560$ Hz έχει απόλυτο μηδενικό πλάτος, αλλά οι γειτονικές της έχουν επίσης χαμηλότερο από μοναδιαίο πλάτος: δείτε τη συχνότητα των $2000$ Hz, το πλάτος της είναι $\\approx 0.7$, αντί μοναδιαίου. Όμοια η συχνότητα των $1300$ Hz, με πλάτος $\\approx 0.6$, και το ίδιο συμβαίνει και με άλλες εκεί γύρω. Αυτό δεν είναι κάτι που μας αρέσει - αν μπορούμε να το φτιάξουμε - καθώς μειώνουμε το πλάτος \"χρήσιμων\" συχνοτήτων, δηλ. συχνοτήτων που δε θέλουμε να πειράξουμε. Στην ουσία, το φίλτρο μας δεν είναι πολύ \"αυστηρό\", \"κατεβαίνει απαλά και χαλαρά για να μηδενίσει τη συχνότητα $f_c$, και ανεβαίνει το ίδιο αργά και χαλαρά\" - θα έλεγε κάποιος κρατώντας ένα freddo espresso στο χέρι. 😄\n",
    "\n",
    "Πώς θα φτιάξουμε αυτό το πρόβλημα? Μα ελέγχοντας/αυξάνοντας την τιμή του $Q$ στο φίλτρο μας!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = 1560 # Hz\n",
    "Q = 20 # από 2 σε 20\n",
    "numL, denL = create_laplace_notch_filter(fc, Q)\n",
    "plot_laplace_frequency_response(numL, denL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wow!! 😍😍 Αυτό πράγματι είναι ένα ΠΟΛΥ καλό φίλτρο! Μπορούμε καλύτερα?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = 1560 # Hz\n",
    "Q = 200 # από 20 σε 200\n",
    "numL, denL = create_laplace_notch_filter(fc, Q)\n",
    "plot_laplace_frequency_response(numL, denL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Σκίσαμε! 😃😍 Αυτό είναι ένα εξαιρετικό - και πολύ ακριβό - φίλτρο που αποκόπτει τη συχνότητα των $1560$ Hz. \n",
    "\n",
    "Μη σας προβληματίζει ότι \"μοιάζει\" να μη μηδενίζει πλήρως τη συχνότητα $f_c$ (το σχήμα φαίνεται να σταματάει στο $0.2$ και να μην κατεβαίνει στο $0.0$), είναι ζήτημα απεικόνισης και όχι ουσίας. \n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ας δούμε τώρα ένα σήμα με ημιτονοειδή θόρυβο, αλλά πλέον δε θα έχετε ΚΑΜΙΑ πληροφορία για τις παραμέτρους του - πλην του ότι η συχνότητα θορύβου είναι ακέραιος αριθμός σε Hz!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fs, s = wav.read('BH-noise.wav')  # Φορτώστε το σήμα (BH-noise.wav, που σας δίνεται)\n",
    "s = s / np.max(np.abs(s))         # Κανονικοποιήστε τα πλάτη στο [-1, 1]\n",
    "\n",
    "Audio(s, rate=fs)  # Ας το ακούσουμε - ΧΑΜΗΛΑ ΤΗΝ ΕΝΤΑΣΗ!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Σε αυτήν την άσκηση, είναι προτιμότερο να μην ακολουθήσουμε το μοτίβο του 4ου εργαστηρίου, που διαλέξαμε 12 παράθυρα και είδαμε το μετασχ. Fourier και το φάσμα πλάτους τους ώστε να εντοπίσουμε την άγνωστη συχνότητα. Ας δούμε κάτι διαφορετικό. \n",
    "\n",
    "Ο παρακάτω κώδικας υπολογίζει το **φασματογράφημα** ενός σήματος. Όπως είπαμε και στην προηγούμενη άσκηση, αυτό δεν είναι τίποτε άλλο από μια συλλογή φασμάτων πλάτους, σε μια οπτική πιο βολική, που χρησιμοποιεί χρωματικό κώδικα για να δείξει ποιές συχνότητες είναι πιο ισχυρές από τις άλλες.\n",
    "\n",
    "Μπορείτε να αναγνωρίσετε τη συχνότητα (ή το εύρος στο οποίο βρίσκεται) που έχει \"μολύνει\" το ηχητικό σήμα?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16,6))\n",
    "plt.specgram(x=s, NFFT=2048, Fs=fs, noverlap=32)\n",
    "plt.title('Φασματογράφημα αρχικού σήματος')\n",
    "plt.ylabel('Συχνότητα (Hz)')\n",
    "plt.xlabel('Χρόνος (s)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Μάλλον είναι πιο δύσκολο απ' ό,τι νομίζατε... σίγουρα αναγνωρίζετε στο σχήμα **που** πρέπει να ψάξετε - σε ποιό **εύρος συχνοτήτων** δηλαδή. \n",
    "\n",
    "Στον παρακάτω κώδικα, ορίστε τα όρια $f_1, f_2$ που θέλετε να ψάξετε, και ο κώδικας θα σας επιστρέψει μια \"φέτα\" του παραπάνω φασματογραφήματος στις συχνότητες $[f_1, f_2]$, ώστε να αναγνωρίσετε τη συχνότητα που θέλετε. \n",
    "\n",
    "Ίσως χρειαστεί να πειραματιστείτε 3-4 φορές στις επιλογές σας των $f_1, f_2$ μέχρι να συγκλίνετε σε κάποιο αριθμό."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = # INSERT CODE HERE\n",
    "f2 = # INSERT CODE HERE\n",
    "\n",
    "NFFT = 8192 # Πλήθος των FFT bins (μεγάλος αριθμός για να κρατήσει η ανάλυση στο zoom)\n",
    "\n",
    "# Fast Fourier Transform ενός τυχαίου τμήματος του σήματος\n",
    "X = np.fft.fft(s[1600:2300], NFFT)\n",
    "\n",
    "# Ζουμάρουμε στο f = [f1, f2]\n",
    "Xzoom = X[int(2*f1/fs * NFFT/2):int(2*f2/fs * NFFT/2)]\n",
    "\n",
    "# Άξονας συχνοτήτων f = [f1, f2]\n",
    "f = np.linspace(f1, f2, len(Xzoom))\n",
    "\n",
    "# Απεικόνιση του φάσματος πλάτους με λεπτομέρεια\n",
    "plt.figure(figsize=(12,4))\n",
    "plt.plot(f, np.abs(Xzoom))\n",
    "plt.title('Φάσμα πλάτους |X(f)| στο [f1, f2] που επιλέξατε')\n",
    "plt.xlabel('Συχνότητα (Hz)')\n",
    "plt.ylabel('Πλάτος')\n",
    "plt.grid(which='major', linestyle='-', linewidth=0.5, alpha=0.7)\n",
    "plt.minorticks_on()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "vscode": {
     "languageId": "ini"
    }
   },
   "source": [
    "Έχετε μια καλή \"μαντεψιά\" για τη συχνότητα που θέλετε να κόψετε? Βάλτε την τιμή που νομίζετε παρακάτω."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fc = # INSERT CODE HERE - Ακέραιος αριθμός σε Hz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ώρα να εφαρμόσουμε το φίλτρο μας???!!! \n",
    "\n",
    "Μη βιάζεστε! 😄 Το φίλτρο μας είναι \"αναλογικό\": είναι χτισμένο στο χώρο του Laplace, που είναι μετασχηματισμός ενός σήματος συνεχούς χρόνου. Το \"μολυσμένο\" σήμα που σας δίνεται είναι \"ψηφιακό\". Πρέπει λοιπόν να μετατρέψουμε το αναλογικό φίλτρο σε ψηφιακό. \n",
    "\n",
    "**ΕΥΤΥΧΩΣ** η ``Python`` φροντίζει για όλα! "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = 200 # από 20 σε 200\n",
    "\n",
    "# Φίλτρο στο χώρο του Laplace\n",
    "# ---------------------------\n",
    "# Τροποποιημένη συχνότητα fc για σωστή λειτουργία της μετατροπής από αναλογικό σε ψηφιακό (don't bother, είναι advanced υλικό)\n",
    "fc_warped = (fs / np.pi) * np.tan(np.pi * fc / fs)\n",
    "\n",
    "# Αναλογικό φίλτρο\n",
    "numL, denL = create_laplace_notch_filter(fc_warped, Q)\n",
    "\n",
    "# Φίλτρο στον ψηφιακό κόσμο\n",
    "# -------------------------\n",
    "# Μετατροπή του αναλογικού notch φίλτρου στον \"ψηφιακό\" κόσμο μέσω bilinear transformation\n",
    "b, a = signal.bilinear(numL, denL, fs)\n",
    "\n",
    "# Φιλτράρισμα\n",
    "y = lfilter(b, a, s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ας ακούσουμε το αρχικό:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(s, rate=fs)  # Πριν τον καθαρισμό"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ας ακούσουμε το \"καθαρό\": "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Audio(y, rate=fs)  # Μετά τον καθαρισμό"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Τα καταφέρατε???? 😃 Αν ακούτε ακόμα το ενοχλητικό ημίτονο, δοκιμάστε μια διαφορετική συχνότητα $f_c$, και τρέξτε ξανά τα τελευταία τμήματα του κώδικα."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "***\n",
    "## Παραδοτέα:\n",
    "*** \n",
    "### Μετονομάστε το **συμπληρωμένο και με εκτελεσμένο κάθε κελί κώδικα** notebook σε \n",
    "\n",
    "## *HY215-Lab5-2025-CSDXXXX.ipynb*, \n",
    "\n",
    "### με *XXXX* ο αριθμός μητρώου σας. \n",
    "\n",
    "### Αν είστε με προσωρινό κωδικό λόγω μεταγραφής, τότε ονομάστε το αρχείο σας \n",
    "\n",
    "## *HY215-Lab5-2025-mtgXXXX.ipynb*, \n",
    "\n",
    "### με *XXXX* ο προσωρινός αριθμός μητρώου σας.\n",
    "\n",
    "### Αν είστε από άλλο τμήμα εκτός του CSD, τότε ονομάστε το αρχείο σας ως\n",
    "\n",
    "## *HY215-Lab5-2025-DPTXXXX.ipynb*, \n",
    "\n",
    "### με DTP ένα εκ των {MATH, PHYS} και *XXXX* ο αριθμός μητρώου σας. \n",
    "\n",
    "### Στείλτε το μετονομασμένο notebook ως επισυναπτόμενο αρχείο με e-mail στο sisamaki@csd.uoc.gr . \n",
    "\n",
    "### Στο e-mail θα αναγράφετε:\n",
    "\n",
    "### - Τίτλος: [ΗΥ215] 5ο εργαστήριο Python\n",
    "\n",
    "### και στο σώμα του e-mail τα παρακάτω: \n",
    "\n",
    "### - Ονοματεπώνυμο\n",
    "### - ΑΜ\n",
    "\n",
    "#### **Σημείωση**: μπορείτε να κάνετε επαναπαραδόσεις του εργαστηρίου - αν χρειαστεί. Χρησιμοποιήστε τίτλο\n",
    "\n",
    "#### [ΗΥ215] Επαναπαράδοση 5ου εργαστηρίου Python\n",
    "\n",
    "#### και κρατήστε τα ίδια στοιχεία στο σώμα του e-mail."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ΠΡΟΣΕΞΤΕ: Το παραδοτέο αρχείο σας πρέπει να έχει extension .ipynb, και ΟΧΙ κάτι άλλο!!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.8.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
