Salta al contenuto
Note per Studenti Esercizio - laboratorio 1, Python al posto di MATLAB

Esercizio - laboratorio 1, Python al posto di MATLAB

In questa pagina 14

Testo (laboratorio 1 del corso Teoria dei Segnali, UniPD, file lab1.m, parti 1-22). Il laboratorio d'esame è in MATLAB: qui si rifanno le stesse 22 parti in Python con numpy (vettori e matrici) e matplotlib (grafici), così da avere un'alternativa gratuita che fa gli stessi calcoli. Le parti coprono: variabili, numeri complessi con modulo, fase e coniugato, if, cicli for e while, resto della divisione, una funzione (distanza radar), vettori e matrici, prodotto riga per colonna e prodotto elemento per elemento, zeros/eye/ones/linspace, grafici con curve sovrapposte e subplot, derivata numerica con le differenze finite.

Teoria usata: Segnali a tempo continuo - definizioni e trasformazioniUn segnale a tempo continuo è una funzione complessa di variabile reale $s:\mathbb R\to\mathbb C$ (reale se $s=s^*$). Si descrive con parte reale e immaginaria oppure con modulo e fase (valore principale in $(-\pi,\pi]$). Sui segnali si fanno tre operazioni elementari: ribaltamento $s(-t)$, traslazione $s(t-t_0)$ e cambio di scala $s(at)$. Ogni segnale si scompone in una parte pari e una dispari, in una reale e una immaginaria, in una hermitiana e una antihermitiana. Estensione e durata dicono dove il segnale è diverso da zero; i segnali con estensione nel semiasse positivo sono causali.Segnali a tempo continuo - definizioni e trasformazioni → (un segnale s(t)s(t) si studia al calcolatore come l'elenco dei suoi valori s(tn)s(t_n) nei tempi tnt_n di una griglia, cioè come vettore: vedi Segnali a tempo discretoUn segnale a tempo discreto è una funzione complessa $s(nT)$ definita sui multipli interi del quanto temporale $T$ (insieme $\mathbb Z(T)$, velocità $F_p=1/T$). Le definizioni sono quelle dei segnali continui con la somma al posto dell'integrale e il quanto $T$ al posto di $dt$: area $\sum T,s(nT)$, energia $\sum T|s(nT)|^2$, convoluzione $\sum T,x(kT)y(nT-kT)$. L'impulso ideale discreto vale $1/T$ nell'origine. Esponenziali e sinusoidi discreti sono periodici solo se $f_0/F_p$ è razionale e hanno frequenza ambigua a meno di multipli di $F_p$. I segnali periodici con periodo $NT$ sono descritti da $N$ valori e si trattano al calcolatore.Segnali a tempo discreto →), Numeri complessiI numeri complessi estendono i reali introducendo l'unità immaginaria $i$ ($i^2 = -1$) e possono essere rappresentati in forma algebrica, trigonometrica o polare. Tramite la formula di Eulero e le proprietà del modulo e dell'argomento, è possibile calcolare agilmente prodotti, potenze e radici ennesime.Numeri complessi →, Derivata - definizione e significatoLa derivata è il limite del rapporto incrementale; geometricamente è la pendenza della retta tangente. f è derivabile in x0 se e solo se f(x) = f(x0) + f'(x0)(x − x0) + o(x − x0); derivabile implica continua, non viceversa. Derivata destra e sinistra, punti angolosi, flessi a tangente verticale, cuspidi.Derivata - definizione e significato →, Segnali notevoli - gradino, rect, tri, sinc ed esponenzialiI segnali di uso più frequente sono la costante, la sinusoide $A_0\cos(2\pi f_0t+\varphi_0)$ e l'esponenziale complesso $Ae^{i2\pi f_0t}$ (periodici, a potenza finita), il gradino $\mathbf 1(t)$ e il segno, e gli impulsi a energia finita: $\operatorname{rect}$ (area $D$), $\operatorname{tri}$, $\operatorname{sinc}$ (area $1$), la gaussiana $e^{-\pi t^2}$ e gli esponenziali smorzati. Per ciascuno si sanno a memoria forma, area ed energia; gli altri segnali si ottengono da questi con traslazioni, scalature, somme e differenze.Segnali notevoli - gradino, rect, tri, sinc ed esponenziali → (per il sinc⁡\operatorname{sinc} della parte 21).

Cosa si vuole fare e perché

In MATLAB (e in Python) un segnale a tempo continuo non esiste "in generale": esiste solo un vettore di tempi t e un vettore di valori s della stessa lunghezza, con sn=s(tn)s_n=s(t_n). Tutto il resto del corso al calcolatore (area, energia, convoluzione, trasformata) parte da qui. Questo primo laboratorio non calcola niente di fisico, serve ad imparare gli strumenti: come si scrivono i numeri, come si guardano i vettori, come si disegna un segnale e come si approssima una derivata. L'ultima parte (la derivata) è la prima che usa davvero un'idea di analisi: sostituire un limite con un rapporto incrementale a passo finito.

Prima di cominciare: le differenze MATLAB → Python che servono

Ogni file Python di queste note comincia con

python
import numpy as np
import matplotlib.pyplot as plt

np è la libreria dei vettori, plt quella dei grafici. In uno script, per vedere una figura serve plt.show() (MATLAB la apre da solo). Le differenze che contano, tutte incontrate nelle parti qui sotto:

cosa MATLAB Python
primo elemento di un vettore v(1) v[0] (gli indici partono da 0)
ultimo elemento v(end) v[-1]
fetta di vettore v(2:4) (estremi inclusi) v[1:4] (l'estremo finale è escluso)
unità immaginaria i, j 1j (un numero, non un nome)
potenza a^b a ** b
prodotto elemento per elemento a .* b a * b
prodotto riga per colonna A * B A @ B
trasposta v.' A.T (su un vettore 1D non fa niente)
resto della divisione mod(a,2) a % 2
"e" logico & and (sui vettori: &)
"o" logico barra verticale or (sui vettori: la barra verticale singola)
vettore con passo 1:0.5:3 np.arange(...) (estremo finale escluso)
matrice di zeri zeros(2,3) np.zeros((2,3)) (una tupla)
stampa fprintf("%d", a) print(f"{a:d}")
blocchi if ... end if ...: e indentazione

Parti 1-2 - variabili e stampa

python
# PARTE 1 - variabili e calcoli semplici
a = 3
b = 4

c1 = a + b
c2 = a - b
c3 = a * b
c4 = a / b          # in Python 3 la divisione / e' sempre in virgola mobile
c5 = a ** b         # potenza: ** al posto di ^
c6 = np.sqrt(b)     # (o anche b ** 0.5)

print(c1)
print(c2, c3, c4, c5, c6)

# PARTE 2 - stampare con segnaposto (al posto di fprintf)
print(f"Il prodotto tra a = {a:d} e b = {b:.4f} risulta:  {c3:.4f}")
# equivalente "alla C": print("... a = %d e b = %.4f ..." % (a, b, c3))

Risultato:

7
-1 12 0.75 81 2.0
Il prodotto tra a = 3 e b = 4.0000 risulta:  12.0000

Sono le cinque operazioni elementari: 3+4=73+4=7, 3−4=−13-4=-1, 3⋅4=123\cdot4=12, 3/4=0,753/4=0{,}75, 34=813^4=81, 4=2\sqrt4=2. In Python 3 / 4 dà sempre un numero con la virgola (in MATLAB tutto è già in virgola mobile). La stringa di fprintf con segnaposto %d (intero), %.4f (decimale con 4 cifre) diventa una f-string: {a:d}, {b:.4f}. Un dettaglio: :d pretende un intero (qui a=3 lo è); con un numero decimale si userebbe :.0f. Un'altra differenza: print va a capo da solo, fprintf no se manca \n.

Parte 3 - numeri complessi

python
# PARTE 3 - numeri complessi: l'unita' immaginaria e' 1j (non i, non j)
a = 1 + 1j
b = 3 + 4j
c = -3 + 0j
d = complex(10, 11)
e = np.exp(1j * np.pi / 4)
f = 1j ** 2
g = np.sqrt(-1 + 0j)      # np.sqrt(-1) darebbe nan: l'argomento deve essere complesso

somma = a + b
prodotto = a * b
rapporto = a / b

modulo_b = np.abs(b)         # MATLAB: abs(b)
fase_c = np.angle(c)         # MATLAB: angle(c)
coniugato_a = np.conj(a)     # MATLAB: conj(a)
reale_a = np.real(a)         # oppure a.real
immaginaria_b = np.imag(b)   # oppure b.imag

print(d)
print(f"{e:.4f}")
print(f)
print(g)
print(somma)
print(prodotto)
print(rapporto)
print(modulo_b)
print(fase_c)
print(coniugato_a)
print(reale_a)
print(immaginaria_b)

Risultato:

(10+11j)
0.7071+0.7071j
(-1+0j)
1j
(4+5j)
(-1+7j)
(0.28-0.04j)
5.0
3.141592653589793
(1-1j)
1.0
4.0

Spiegazione dei numeri (richiami in Numeri complessiI numeri complessi estendono i reali introducendo l'unità immaginaria $i$ ($i^2 = -1$) e possono essere rappresentati in forma algebrica, trigonometrica o polare. Tramite la formula di Eulero e le proprietà del modulo e dell'argomento, è possibile calcolare agilmente prodotti, potenze e radici ennesime.Numeri complessi →):

  • il modulo di 3+4i3+4i è 32+42=5\sqrt{3^2+4^2}=5;
  • la fase di −3-3 è π≈3,1416\pi\approx3{,}1416 (la fase è sempre l'angolo principale, tra −π-\pi e π\pi; −3-3 sta sul semiasse reale negativo);
  • (1+i)(3+4i)=3+4i+3i+4i2=−1+7i(1+i)(3+4i)=3+4i+3i+4i^2=-1+7i;
  • 1+i3+4i=(1+i)(3−4i)25=7−i25=0,28−0,04 i\dfrac{1+i}{3+4i}=\dfrac{(1+i)(3-4i)}{25}=\dfrac{7-i}{25}=0{,}28-0{,}04\,i (si moltiplica sopra e sotto per il coniugato del denominatore);
  • eiπ/4=cos⁡π4+isin⁡π4=0,7071+0,7071 ie^{i\pi/4}=\cos\frac\pi4+i\sin\frac\pi4=0{,}7071+0{,}7071\,i, di modulo 1;
  • i2=−1i^2=-1.

Due trappole. In Python np.sqrt(-1) dà nan (la radice di un reale negativo non è definita tra i reali): l'argomento va scritto complesso, -1 + 0j. Inoltre in MATLAB i e j sono nomi che si possono sovrascrivere (la parte 7 usa i come contatore del ciclo, e da quel momento i non è più −1\sqrt{-1}); in Python l'unità è il numero 1j e nessun nome la può rovinare.

Parti 4-6 - if, operatori logici, resto

python
# PARTE 4 - if / elif / else (i due punti e l'indentazione sostituiscono then ... end)
a = 57
b = 30

if a == b:
    print("a e b sono uguali")
elif a < b:
    print("a è minore di b")
else:
    print("a è maggiore di b")

# PARTE 5 - operatori logici: or / and (MATLAB: | e &)
if a > 40 or b > 40:
    print("almeno uno tra a e b è maggiore di 40")

if a >= 40 and a <= 60:      # equivale a 40 <= a <= 60
    print("a è compreso tra 40 e 60")
else:
    print("a non è compreso tra 40 e 60")

# PARTE 6 - resto della divisione: % (oppure np.mod)
if a % 2 == 0:
    print("a è pari")
if a % 2 == 1:
    print("a è dispari")

Risultato (con a=57a=57, b=30b=30):

a è maggiore di b
almeno uno tra a e b è maggiore di 40
a è compreso tra 40 e 60
a è dispari

Perché queste risposte: 57>3057>30, quindi vince il ramo else; 57>4057>40 rende vera la condizione "almeno uno è maggiore di 40"; 40≤57≤6040\le57\le60 è vera; 57=2⋅28+157=2\cdot28+1 ha resto 1, quindi è dispari. Per un singolo valore vero/falso si scrivono and/or; per confrontare vettori elemento per elemento (lo si fa spesso per costruire una finestra, vedi i laboratori successivi) si usano & e la barra verticale singola, con le parentesi: (t >= -1) & (t <= 1).

Parti 7-8 - cicli for e while

python
# PARTE 7 - ciclo for (range(1, n_iter+1) e' MATLAB 1:1:n_iter, estremo finale escluso)
n_iter = 5
s = 0
for i in range(1, n_iter + 1):
    s = s + 10
    print(f"Iterazione n° {i}, s = {s:2.1f}")

# PARTE 8 - ciclo while: stessa cosa di prima
s = 0
i = 1
while i <= n_iter:
    s = s + 10
    print(f"Iterazione n° {i}, s = {s:2.1f}")
    i = i + 1

Risultato (i due cicli stampano lo stesso):

Iterazione n° 1, s = 10.0
Iterazione n° 2, s = 20.0
Iterazione n° 3, s = 30.0
Iterazione n° 4, s = 40.0
Iterazione n° 5, s = 50.0
Iterazione n° 1, s = 10.0
Iterazione n° 2, s = 20.0
Iterazione n° 3, s = 30.0
Iterazione n° 4, s = 40.0
Iterazione n° 5, s = 50.0

range(1, n_iter + 1) produce 1,2,…,51,2,\dots,5: l'estremo finale è escluso, per questo si scrive n_iter + 1. È la differenza classica con 1:1:n_iter di MATLAB, che include l'ultimo valore. Dopo cinque passi s=5⋅10=50s=5\cdot10=50. Il while ripete finché la condizione è vera, quindi il contatore va incrementato a mano. Altre note sul sito: Ciclo for e rangefor scorre gli elementi di un iterabile; range per gli indici; enumerate, zip, reversed, sorted; cicli annidati; quando usare for e quando while.Ciclo for e range →, Ciclo whileCiclo while, cicli controllati da contatore, da sentinella e da condizione; break, continue, else; terminazione e invariante di ciclo.Ciclo while →.

Parte 9 - una funzione: il radar

python
# PARTE 9 - una funzione: distanza dal tempo di ritorno di un impulso radar
def calcolatore_distanza(delta_t):
    if delta_t < 0:
        raise ValueError("Il tempo di propagazione non può essere negativo!")
    c = 3e8                       # velocita' della luce [m/s]
    d = (delta_t * c) / 2         # l'impulso fa andata e ritorno: meta' del percorso
    return d

delta_t = 0.5e-6                  # 0,5 microsecondi
distanza = calcolatore_distanza(delta_t)
print(f"La distanza è {distanza:2.4f} metri")

Risultato:

La distanza è 75.0000 metri

Un radar manda un impulso e misura dopo quanto tempo Δt\Delta t torna l'eco. L'impulso viaggia alla velocità della luce c=3⋅108c=3\cdot10^8 m/s e percorre il tragitto due volte (andata e ritorno): 2d=c Δt2d=c\,\Delta t, quindi d=c Δt2=3⋅108⋅0,5⋅10−62=75d=\dfrac{c\,\Delta t}{2}=\dfrac{3\cdot10^8\cdot0{,}5\cdot10^{-6}}{2}=75 m. Il controllo delta_t < 0 evita un tempo senza senso fisico (in MATLAB si usa error(...), in Python si "solleva" un'eccezione con raise). Una funzione Python si definisce con def e restituisce con return; in MATLAB stava in un file a parte, qui basta definirla nello stesso script (Funzioni in Pythondef e return, parametri posizionali, con nome e con valore predefinito, *args e **kwargs; passaggio per riferimento a oggetto; ambito delle variabili LEGB; docstring, type hint, lambda.Funzioni in Python →).

Parti 10-13 - vettori e matrici

python
# PARTE 10 - vettori: gli indici partono da 0
v = np.array([11, 13, 17, 19, 23])

print(v)
print(v[1])          # MATLAB v(2): il secondo elemento
print(v[-1])         # MATLAB v(end): l'ultimo
print(f"La lunghezza del vettore è {len(v)}")

# PARTE 11 - operazioni tra vettori
w = np.array([1, 1, 2, 2, 0])
somma_v_w = v + w
sottrazione_v_w = v - w
print(w)
print(somma_v_w)
print(sottrazione_v_w)
print(v.T)                       # un array 1D non ha "riga" o "colonna": .T non fa nulla
print(v.reshape(-1, 1))          # per ottenere la colonna (5x1) serve reshape (o v[:, None])

# PARTE 12 - matrici (array 2D)
m = np.array([[2, 4, 6, 8, 10],
              [3, 6, 9, 12, 15],
              [4, 8, 12, 16, 20]])
print(m)
print(m.shape)                   # MATLAB size(m)

# da un vettore di 15 elementi: reshape riempie per RIGHE (MATLAB riempie per colonne)
vettore15 = np.array([2, 4, 6, 8, 10, 3, 6, 9, 12, 15, 4, 8, 12, 16, 20])
m2 = vettore15.reshape(3, 5)
print(np.array_equal(m, m2))

# PARTE 13 - accesso agli elementi di una matrice
print(m[1, 3])                   # MATLAB m(2,4)
print(m[0, :])                   # prima riga
print(m[:, 2])                   # terza colonna (esce come array 1D)

Risultato:

[11 13 17 19 23]
13
23
La lunghezza del vettore è 5
[1 1 2 2 0]
[12 14 19 21 23]
[10 12 15 17 23]
[11 13 17 19 23]
[[11]
 [13]
 [17]
 [19]
 [23]]
[[ 2  4  6  8 10]
 [ 3  6  9 12 15]
 [ 4  8 12 16 20]]
(3, 5)
True
12
[ 2  4  6  8 10]
[ 6  9 12]

Le regole importanti:

  • Indici da 0. v[1] è il secondo elemento (13), v[-1] l'ultimo (23). La matrice m ha 3 righe e 5 colonne: m[1, 3] è la riga di indice 1 (la seconda) e la colonna di indice 3 (la quarta), cioè 1212 (in MATLAB m(2,4)).
  • m.shape è (3, 5): è l'equivalente di size(m).
  • Un vettore numpy a una dimensione non è né riga né colonna: v.T lo lascia com'è. Per ottenere la colonna 5×15\times1 della parte 11 si usa v.reshape(-1, 1) (il -1 vuol dire "calcola tu quanti").
  • reshape riempie per righe, MATLAB per colonne: il vettore di 15 elementi dà la stessa matrice 3×53\times5 perché è stato scritto riga dopo riga.
  • m[:, 2] (terza colonna) è un array 1D e si stampa in orizzontale; in MATLAB si vedrebbe in colonna.

Parti 14-17 - prodotti tra vettori

python
# PARTE 14 - prodotto riga per colonna: operatore @ (non *)
v1 = np.array([[10, 20, 30]])        # 1x3  (riga)
v2 = np.array([[1], [2], [3]])       # 3x1  (colonna)

v3 = v1 @ v2          # (1x3)(3x1) -> 1x1 : 10*1 + 20*2 + 30*3
v4 = v2 @ v1          # (3x1)(1x3) -> 3x3
print(v3)
print(v4)

# PARTE 15 - prodotto elemento per elemento: in numpy e' il * semplice (MATLAB .*)
a = np.array([100, 200, 300])
b = np.array([3, 1, 2])
c = a * b
print(c)

# PARTE 16 - il prodotto tra matrici con dimensioni incompatibili da' errore
A = a.reshape(1, 3)
B = b.reshape(1, 3)
try:
    A @ B                       # (1x3)(1x3): le dimensioni interne non coincidono
except ValueError as errore:
    print("Errore:", errore)
print(a * b)                    # invece * funziona sempre (elemento per elemento)

# PARTE 17 - trasposta ed espansione: x (riga) per y colonna
x = np.array([1, 2, 3])
y = np.array([10, 10, 10])

z1 = (x * y).reshape(-1, 1)                  # MATLAB (x.*y).'  -> colonna 3x1
print(z1)
z2 = x.reshape(1, 3) * y.reshape(3, 1)       # MATLAB x.*y.'  -> matrice 3x3 (broadcasting)
print(z2)

Risultato:

[[140]]
[[10 20 30]
 [20 40 60]
 [30 60 90]]
[300 200 600]
Errore: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signature (n?,k),(k,m?)->(n?,m?) (size 1 is different from 3)
[300 200 600]
[[10]
 [20]
 [30]]
[[10 20 30]
 [10 20 30]
 [10 20 30]]
  • Parte 14. v1 è 1×31\times3, v2 è 3×13\times1. Il prodotto riga per colonna v1 @ v2 è un numero: 10⋅1+20⋅2+30⋅3=14010\cdot1+20\cdot2+30\cdot3=140 (ciascun elemento della riga per il suo corrispondente nella colonna, poi somma; è il prodotto scalare). Il prodotto v2 @ v1 è la matrice 3×33\times3 che ha in posizione (r,c)(r,c) il prodotto v2r v1cv2_r\,v1_c. Attenzione: in numpy l'operatore * non è il prodotto tra matrici ma quello elemento per elemento (con espansione), quindi per questi prodotti serve @.
  • Parte 15. Elemento per elemento: 100⋅3, 200⋅1, 300⋅2=300, 200, 600100\cdot3,\ 200\cdot1,\ 300\cdot2=300,\,200,\,600. In MATLAB è .*, in Python è semplicemente *.
  • Parte 16. In MATLAB a*b dà errore (due righe non si possono moltiplicare tra loro); in Python il corrispondente errore è quello di A @ B con due matrici 1×31\times3: le dimensioni "interne" devono coincidere, qui sono 33 e 11. Invece a * b in Python funziona sempre.
  • Parte 17. z1 è il vettore trasposto in colonna. z2 mostra l'espansione (broadcasting): una riga 1×31\times3 per una colonna 3×13\times1 produce la matrice 3×33\times3 i cui elementi sono xc yrx_c\,y_r. In MATLAB lo fa x.*y.'; in numpy funziona solo se i due array sono 2D come qui (con due array 1D si resterebbe in lunghezza 3).

Parte 18 - matrici e vettori speciali

python
# PARTE 18 - matrici e vettori speciali
matrice_zeri = np.zeros((2, 3))         # il formato e' UNA tupla: (righe, colonne)
matrice_identity = np.eye(4)
matrice_uno = np.ones((3, 2))
vec_linspace = np.linspace(1, 3, 30)    # 30 punti, estremi 1 e 3 inclusi (come MATLAB)
vec_step = 15 + 0.2 * np.arange(11)     # da 15 a 17 con passo 0,2 (11 valori, 17 incluso)

print(matrice_zeri)
print(matrice_identity)
print(matrice_uno)
print(vec_linspace[:5], "...", vec_linspace[-1])
print(vec_step)
print(len(np.arange(15, 17, 0.2)), "elementi con np.arange(15, 17, 0.2): il 17 e' escluso")

Risultato:

[[0. 0. 0.]
 [0. 0. 0.]]
[[1. 0. 0. 0.]
 [0. 1. 0. 0.]
 [0. 0. 1. 0.]
 [0. 0. 0. 1.]]
[[1. 1.]
 [1. 1.]
 [1. 1.]]
[1.         1.06896552 1.13793103 1.20689655 1.27586207] ... 3.0
[15.  15.2 15.4 15.6 15.8 16.  16.2 16.4 16.6 16.8 17. ]
10 elementi con np.arange(15, 17, 0.2): il 17 e' escluso

np.linspace(1, 3, 30) dà 30 punti equispaziati tra 1 e 3, estremi inclusi (passo 2/292/29), come in MATLAB. Con un passo decimale, invece, np.arange(15, 17, 0.2) non include il 17 e per gli arrotondamenti dei decimali può anche sbagliare di un elemento: l'asse tempi si costruisce meglio con il numero di campioni, inizio + np.arange(n_campioni) * passo, dove n=fine−iniziopasso+1n=\frac{\text{fine}-\text{inizio}}{\text{passo}}+1 (qui 17−150,2+1=11\frac{17-15}{0{,}2}+1=11). Questa scrittura è usata in tutte le note successive.

Parti 19-20 - grafici e curve sovrapposte

python
# PARTE 19 - grafico di una sinusoide
f = 5
t_inizio = 0
t_fine = 1
ts = 1 / 1000

# asse dei tempi: MATLAB [t_inizio:ts:t_fine]; in Python si conta il numero di campioni
# (np.arange con passo decimale puo' dare un elemento in piu' o in meno)
t = t_inizio + np.arange(round((t_fine - t_inizio) / ts) + 1) * ts
y = np.sin(2 * np.pi * f * t)
print(len(t), t[0], t[-1])

plt.figure()
plt.plot(t, y, '--', color='r')          # linea tratteggiata rossa
plt.xlabel("tempo [s]")
plt.ylabel("y")
plt.title("Grafico di una sinusoide")
plt.grid(True)
plt.show()

# solo pallini, passo grande (ts = 0.05): si vede che i campioni sono pochi
t_grande = np.arange(21) * 0.05
plt.figure()
plt.plot(t_grande, np.sin(2 * np.pi * f * t_grande), 'o')
plt.show()

# PARTE 20 - curve sovrapposte
w = np.cos(2 * np.pi * f * t)

plt.figure()
plt.plot(t, y)           # in matplotlib le curve successive restano sullo stesso asse
plt.plot(t, w)           # (non serve "hold on")
plt.xlabel("tempo [s]")
plt.title("Grafico con due curve")
plt.grid(True)
plt.show()

plt.figure()             # alternativa: due curve in una sola chiamata
plt.plot(t, y, t, w)
plt.show()

Risultato:

1001 0.0 1.0

Il segnale è y(t)=sin⁡(2π⋅5 t)y(t)=\sin(2\pi\cdot5\,t), una sinusoide di frequenza f=5f=5 Hz (periodo 0,20{,}2 s) sull'intervallo [0,1][0,1], campionata con passo ts=1t_s=1 ms: ci vogliono 1−00,001+1=1001\frac{1-0}{0{,}001}+1=1001 campioni, come stampato. Con soli pallini e un passo grande (50 ms, 4 campioni per periodo) si vede che la forma d'onda si perde: è l'anticipo del problema del campionamento (Campionamento e ricostruzioneIl campionamento $s_c(nT)=s(nT)$ trasforma un segnale continuo in uno discreto e, in frequenza, ripete lo spettro con periodo $F_c=1/T$: $S_c(f)=\sum_kS(f-kF_c)$. Se le repliche si sovrappongono si ha aliasing e il segnale non è recuperabile. Un interpolatore $\mathbb Z(T)\to\mathbb R$ con risposta impulsiva $g$ produce $\tilde s(t)=\sum_nT g(t-nT)s(nT)$ e $\tilde S=G,S_c$. Teorema del campionamento: se $s$ ha banda $B$ e $F_c\ge2B$, l'interpolatore ideale $G=\operatorname{rect}(f/F_c)$ ricostruisce esattamente $s(t)=\sum s(nT)\operatorname{sinc}(F_c(t-nT))$. Se le ipotesi non valgono c'è un errore, in banda e fuori banda, riducibile con un prefiltro anti-aliasing.Campionamento e ricostruzione →). In matplotlib le curve disegnate una dopo l'altra restano nello stesso grafico (non serve hold on) e una nuova figura si apre con plt.figure(). Le sigle 'b-', 'r--', 'g-.', 'm:', 'o' sono le stesse di MATLAB.

Grafico interattivo: y(t) = sin(2π·5·t) e w(t) = cos(2π·5·t) per t tra 0 e 1 secondo: 5 periodi in un secondo; il coseno è il seno anticipato di un quarto di periodo

Parte 21 - sei grafici in una figura

python
# PARTE 21 - sei grafici in una figura (subplot righe, colonne, posizione: da 1, come MATLAB)
t = np.linspace(-2 * np.pi, 2 * np.pi, 1000)

plt.figure(figsize=(12, 6))

plt.subplot(2, 3, 1)
plt.plot(t, np.sin(t), 'b-', linewidth=2)
plt.title('sen(t)'); plt.grid(True)

plt.subplot(2, 3, 2)
plt.plot(t, np.cos(t), 'r--', linewidth=2)
plt.title('cos(t)'); plt.grid(True)

plt.subplot(2, 3, 3)
plt.plot(t, np.tan(t), 'g-.', linewidth=2)
plt.ylim(-10, 10)
plt.title('tan(t)'); plt.grid(True)

plt.subplot(2, 3, 4)
plt.plot(t, np.exp(-t), 'm:', linewidth=2)
plt.title('exp(-t)'); plt.grid(True)

plt.subplot(2, 3, 5)
plt.plot(t, np.exp(t), 'c-', linewidth=2)
plt.title('exp(t)'); plt.grid(True)

plt.subplot(2, 3, 6)
plt.plot(t, np.sinc(t), 'k', linewidth=2)   # np.sinc(x) = sin(pi x)/(pi x), come sinc di MATLAB
plt.title('sinc(t)'); plt.grid(True)

for k in range(1, 7):
    plt.subplot(2, 3, k)
    plt.xlabel('Tempo (t)')
    plt.ylabel('Ampiezza')
plt.tight_layout()
plt.show()
print(np.sinc(0.5), 2 / np.pi)

Risultato:

0.6366197723675814 0.6366197723675814

plt.subplot(2, 3, k) divide la figura in 2 righe e 3 colonne e sceglie la casella kk (da 1, come in MATLAB, riempiendo per righe). np.sinc(x) è già il sinc normalizzato sinc⁡(x)=sin⁡πxπx\operatorname{sinc}(x)=\frac{\sin\pi x}{\pi x}, lo stesso del corso e del sinc di MATLAB; la riga stampata verifica sinc⁡(0,5)=sin⁡(π/2)π/2=2π≈0,6366\operatorname{sinc}(0{,}5)=\frac{\sin(\pi/2)}{\pi/2}=\frac2\pi\approx0{,}6366. La tangente ha asintoti verticali e i suoi valori esplodono, perciò si limita l'asse con plt.ylim(-10, 10). Il ylim [-10,10] nel grafico del sinc, nel file MATLAB, non serve: il sinc sta tra −0,22-0{,}22 e 11.

Parte 22 - la derivata numerica

python
# PARTE 22 - derivata numerica con le differenze finite
def calcolatore_derivata(f, t):
    """Derivata approssimata di f definita sui tempi t. Restituisce le tre versioni."""
    if len(f) != len(t):
        raise ValueError("f e t hanno lunghezze diverse!")

    # (1) con un ciclo: (f[n+1] - f[n]) / (t[n+1] - t[n])
    fd = np.zeros(len(t))
    for n in range(len(t) - 1):              # MATLAB: for n = 1:length(t)-1
        fd[n] = (f[n + 1] - f[n]) / (t[n + 1] - t[n])
    fd[-1] = fd[-2]                           # all'ultimo istante manca il campione successivo

    # (2) con np.diff (differenze tra elementi consecutivi)
    fd_diff = np.diff(f) / np.diff(t)
    fd_diff = np.append(fd_diff, fd_diff[-1])

    # (3) vettorizzata a mano con le fette: f[1:] sono i "successivi", f[:-1] i "precedenti"
    fd_vec = (f[1:] - f[:-1]) / (t[1:] - t[:-1])
    fd_vec = np.append(fd_vec, fd_vec[-1])

    return fd, fd_diff, fd_vec


t_start, t_end, t_step = 0, 6, 0.001
t = t_start + np.arange(round((t_end - t_start) / t_step) + 1) * t_step
f = np.sin(2 * np.pi * 1 * t)

fd, fd_diff, fd_vec = calcolatore_derivata(f, t)
teorica = 2 * np.pi * np.cos(2 * np.pi * 1 * t)

print(len(t), "campioni")
print("le tre versioni coincidono:", np.allclose(fd, fd_diff), np.allclose(fd, fd_vec))
print(f"errore massimo (derivata stimata - teorica): {np.max(np.abs(fd - teorica)):.5f}")

plt.figure()
plt.subplot(2, 1, 1)
plt.plot(t, f, 'b', linewidth=1)
plt.subplot(2, 1, 2)
plt.plot(t, fd, 'r--', linewidth=3, label="derivata stimata")
plt.plot(t, teorica, 'k', linewidth=1.5, label="derivata teorica")
plt.xlabel('Tempo (t)')
plt.ylabel('Ampiezza')
plt.legend()
plt.grid(True)
plt.show()

Risultato:

6001 campioni
le tre versioni coincidono: True True
errore massimo (derivata stimata - teorica): 0.01974

Idea. La derivata è il limite del rapporto incrementale, f′(t)=lim⁡h→0f(t+h)−f(t)hf'(t)=\lim_{h\to0}\frac{f(t+h)-f(t)}{h} (Derivata - definizione e significatoLa derivata è il limite del rapporto incrementale; geometricamente è la pendenza della retta tangente. f è derivabile in x0 se e solo se f(x) = f(x0) + f'(x0)(x − x0) + o(x − x0); derivabile implica continua, non viceversa. Derivata destra e sinistra, punti angolosi, flessi a tangente verticale, cuspidi.Derivata - definizione e significato →). Su una griglia non si può far tendere hh a zero: si prende h=Δt=tn+1−tnh=\Delta t=t_{n+1}-t_n, il passo di campionamento, e si calcola

f′(tn)≈f(tn+1)−f(tn)tn+1−tn.f'(t_n)\approx\frac{f(t_{n+1})-f(t_n)}{t_{n+1}-t_n}.

Il risultato ha un campione in meno (l'ultimo nn non ha un successivo): per mantenere la stessa lunghezza di t si ricopia l'ultimo valore, fd[-1] = fd[-2].

Tre scritture dello stesso calcolo. (1) un ciclo che scorre gli indici; (2) np.diff(f) calcola le differenze tra elementi consecutivi in un colpo solo, quindi np.diff(f) / np.diff(t) è il rapporto incrementale di tutti i campioni; (3) la stessa cosa con le "fette": f[1:] sono i campioni successivi e f[:-1] i precedenti. Le tre coincidono (True True). Le versioni (2) e (3) sono molto più veloci del ciclo con migliaia di campioni. Nel file MATLAB la funzione restituisce solo la versione con il ciclo e le altre due sono lasciate "per confronto"; qui restituiscono tutte e tre per poterle confrontare.

Confronto con la derivata teorica. Per f(t)=sin⁡(2πt)f(t)=\sin(2\pi t) la derivata esatta è f′(t)=2πcos⁡(2πt)f'(t)=2\pi\cos(2\pi t), di ampiezza 2π≈6,282\pi\approx6{,}28. L'errore massimo con Δt=1\Delta t=1 ms è 0,01970{,}0197, cioè circa lo 0,3%0{,}3\% dell'ampiezza: nel grafico le due curve sono sovrapposte.

Grafico interattivo: f(t) = sin(2πt) e la sua derivata teorica f'(t) = 2π cos(2πt) (dalla parte 22): la derivata massima, 2π, si ha dove f passa per zero in salita, e vale zero nei massimi e nei minimi di f

Ruolo del passo. Il secondo blocco di codice ripete il calcolo con quattro passi:

python
# Come dipende l'errore dal passo
print("passo    errore massimo   errore/passo")
for passo in [0.1, 0.01, 0.001, 0.0001]:
    tt = np.arange(round(6 / passo) + 1) * passo
    ff = np.sin(2 * np.pi * tt)
    d = np.diff(ff) / np.diff(tt)
    err = np.max(np.abs(d - 2 * np.pi * np.cos(2 * np.pi * tt[:-1])))
    # la differenza in avanti stima la derivata nel punto di mezzo t + passo/2
    err_mezzo = np.max(np.abs(d - 2 * np.pi * np.cos(2 * np.pi * (tt[:-1] + passo / 2))))
    print(f"{passo:<8} {err:<16.6f} {err / passo:<12.2f} errore a meta' passo: {err_mezzo:.2e}")
passo    errore massimo   errore/passo
0.1      1.941611         19.42        errore a meta' passo: 9.78e-02
0.01     0.197327         19.73        errore a meta' passo: 1.03e-03
0.001    0.019739         19.74        errore a meta' passo: 1.03e-05
0.0001   0.001974         19.74        errore a meta' passo: 1.03e-07

L'errore massimo è proporzionale al passo (dividendo per Δt\Delta t si ottiene sempre ≈19,7\approx19{,}7). Si capisce con lo sviluppo di Taylor: f(t+h)=f(t)+hf′(t)+h22f′′(t)+…f(t+h)=f(t)+hf'(t)+\frac{h^2}{2}f''(t)+\dots, quindi il rapporto incrementale vale f′(t)+h2f′′(t)+…f'(t)+\frac h2f''(t)+\dots e l'errore è h2f′′\frac{h}{2}f''; per il seno ∣f′′∣max⁡=(2π)2=39,48|f''|_{\max}=(2\pi)^2=39{,}48 e l'errore massimo è h2⋅39,48=19,74 h\frac h2\cdot39{,}48=19{,}74\,h, proprio i numeri in tabella. Il rapporto incrementale in avanti stima in realtà la derivata nel punto di mezzo tn+Δt2t_n+\frac{\Delta t}2: confrontandolo con 2πcos⁡(2π(tn+Δt2))2\pi\cos\bigl(2\pi(t_n+\frac{\Delta t}2)\bigr) l'errore scende a 10−510^{-5} per Δt=10−3\Delta t=10^{-3} (ultima colonna). Morale: un passo dieci volte più piccolo dà un errore dieci volte più piccolo, a prezzo di dieci volte più campioni.

Controllo

Tutti i numeri sopra sono stati ottenuti eseguendo il codice (con matplotlib.use("Agg"), che disegna senza aprire finestre). Verifiche a mano: c5=34=81c_5=3^4=81; distanza radar 3⋅108⋅5⋅10−72=75\frac{3\cdot10^8\cdot5\cdot10^{-7}}2=75 m; prodotto scalare 10+40+90=14010+40+90=140; 10011001 campioni per [0,1][0,1] con passo 10−310^{-3}; errore di derivata ≈Δt2(2π)2\approx\frac{\Delta t}2(2\pi)^2.

Versione ripasso

  • Differenze da ricordare. Indici da 0 e estremo finale escluso (v[1:4], range(1, n+1)); 1j per ii; * elemento per elemento, @ riga per colonna; np.zeros((2,3)) con una tupla; plt.show(); np.sinc è già normalizzato; asse tempi t0 + np.arange(n)*dt, mai arange con passo decimale.
  • Complessi. ∣3+4i∣=5|3+4i|=5, arg⁡(−3)=π\arg(-3)=\pi, (1+i)(3+4i)=−1+7i(1+i)(3+4i)=-1+7i, 1+i3+4i=0,28−0,04i\frac{1+i}{3+4i}=0{,}28-0{,}04i.
  • Radar. d=c Δt2d=\frac{c\,\Delta t}2 (andata e ritorno): Δt=0,5 μ\Delta t=0{,}5\ \mus dà 7575 m.
  • Prodotti. (1×3)(3×1)=140(1\times3)(3\times1)=140, (3×1)(1×3)(3\times1)(1\times3) è 3×33\times3; broadcasting riga per colonna.
  • Derivata. f′(tn)≈f(tn+1)−f(tn)Δtf'(t_n)\approx\frac{f(t_{n+1})-f(t_n)}{\Delta t} (ciclo, np.diff, fette); errore ≈Δt2∣f′′∣=19,74 Δt\approx\frac{\Delta t}2\lvert f''\rvert=19{,}74\,\Delta t per sin⁡2πt\sin2\pi t: proporzionale al passo.

Teoria collegata