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 trasformazioni (un segnale $s(t)$ si studia al calcolatore come l'elenco dei suoi valori $s(t_n)$ nei tempi $t_n$ di una griglia, cioè come vettore: vedi Segnali a tempo discreto), Numeri complessi, Derivata - definizione e significato, Segnali notevoli - gradino, rect, tri, sinc ed esponenziali (per il $\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 $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=7$ , $3-4=-1$ , $3\cdot4=12$ , $3/4=0{,}75$ , $3^4=81$ , $\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 complessi):

il modulo di $3+4i$ è $\sqrt{3^2+4^2}=5$ ;
la fase di $-3$ è $\pi\approx3{,}1416$ (la fase è sempre l'angolo principale, tra $-\pi$ e $\pi$ ; $-3$ sta sul semiasse reale negativo);
$(1+i)(3+4i)=3+4i+3i+4i^2=-1+7i$ ;
$\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);
$e^{i\pi/4}=\cos\frac\pi4+i\sin\frac\pi4=0{,}7071+0{,}7071\,i$ , di modulo 1;
$i^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ù $\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=57$ , $b=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>30$ , quindi vince il ramo else; $57>40$ rende vera la condizione "almeno uno è maggiore di 40"; $40\le57\le60$ è vera; $57=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,\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\cdot10=50$ . Il while ripete finché la condizione è vera, quindi il contatore va incrementato a mano. Altre note sul sito: Ciclo for e range, 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 $\Delta t$ torna l'eco. L'impulso viaggia alla velocità della luce $c=3\cdot10^8$ m/s e percorre il tragitto due volte (andata e ritorno): $2d=c\,\Delta t$ , quindi $d=\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 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è $12$ (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\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\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\times3$ , v2 è $3\times1$ . Il prodotto riga per colonna v1 @ v2 è un numero: $10\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\times3$ che ha in posizione $(r,c)$ il prodotto $v2_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\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\times3$ : le dimensioni "interne" devono coincidere, qui sono $3$ e $1$ . Invece a * b in Python funziona sempre.
Parte 17. z1 è il vettore trasposto in colonna. z2 mostra l'espansione (broadcasting): una riga $1\times3$ per una colonna $3\times1$ produce la matrice $3\times3$ i cui elementi sono $x_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/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=\frac{\text{fine}-\text{inizio}}{\text{passo}}+1$ (qui $\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\pi\cdot5\,t)$ , una sinusoide di frequenza $f=5$ Hz (periodo $0{,}2$ s) sull'intervallo $[0,1]$ , campionata con passo $t_s=1$ ms: ci vogliono $\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 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 $k$ (da 1, come in MATLAB, riempiendo per righe). np.sinc(x) è già il sinc normalizzato $\operatorname{sinc}(x)=\frac{\sin\pi x}{\pi x}$ , lo stesso del corso e del sinc di MATLAB; la riga stampata verifica $\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$ e $1$ .

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\to0}\frac{f(t+h)-f(t)}{h}$ (Derivata - definizione e significato). Su una griglia non si può far tendere $h$ a zero: si prende $h=\Delta t=t_{n+1}-t_n$ , il passo di campionamento, e si calcola

$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 $n$ 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\pi t)$ la derivata esatta è $f'(t)=2\pi\cos(2\pi t)$ , di ampiezza $2\pi\approx6{,}28$ . L'errore massimo con $\Delta t=1$ ms è $0{,}0197$ , cioè circa lo $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 $\Delta t$ si ottiene sempre $\approx19{,}7$ ). Si capisce con lo sviluppo di Taylor: $f(t+h)=f(t)+hf'(t)+\frac{h^2}{2}f''(t)+\dots$ , quindi il rapporto incrementale vale $f'(t)+\frac h2f''(t)+\dots$ e l'errore è $\frac{h}{2}f''$ ; per il seno $|f''|_{\max}=(2\pi)^2=39{,}48$ e l'errore massimo è $\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 $t_n+\frac{\Delta t}2$ : confrontandolo con $2\pi\cos\bigl(2\pi(t_n+\frac{\Delta t}2)\bigr)$ l'errore scende a $10^{-5}$ per $\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: $c_5=3^4=81$ ; distanza radar $\frac{3\cdot10^8\cdot5\cdot10^{-7}}2=75$ m; prodotto scalare $10+40+90=140$ ; $1001$ campioni per $[0,1]$ con passo $10^{-3}$ ; errore di derivata $\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 $i$ ; * 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$ , $\arg(-3)=\pi$ , $(1+i)(3+4i)=-1+7i$ , $\frac{1+i}{3+4i}=0{,}28-0{,}04i$ .
Radar. $d=\frac{c\,\Delta t}2$ (andata e ritorno): $\Delta t=0{,}5\ \mu$ s dà $75$ m.
Prodotti. $(1\times3)(3\times1)=140$ , $(3\times1)(1\times3)$ è $3\times3$ ; broadcasting riga per colonna.
Derivata. $f'(t_n)\approx\frac{f(t_{n+1})-f(t_n)}{\Delta t}$ (ciclo, np.diff, fette); errore $\approx\frac{\Delta t}2\lvert f''\rvert=19{,}74\,\Delta t$ per $\sin2\pi t$ : proporzionale al passo.
Parti 1-2, calcoli e stampa. $3+4=7$ , $3-4=-1$ , $3\cdot4=12$ , $3/4=0{,}75$ , $3^4=81$ , $\sqrt4=2$ . In Python 3/4 dà sempre un valore in virgola; {a:d} vuole un intero; print va a capo da solo.
Parte 3, complessi. Con $a=1+i$ e $b=3+4i$ : $\lvert b\rvert=5$ , $\arg(-3)=\pi$ (valore principale), $(1+i)(3+4i)=-1+7i$ , $\frac{1+i}{3+4i}=0{,}28-0{,}04i$ , $e^{i\pi/4}=0{,}7071+0{,}7071i$ . np.sqrt(-1) dà nan: serve -1+0j.
Parti 4-6, if e logica. Con $a=57$ e $b=30$ : vince else ( $a>b$ ); or è vera perché $57>40$ ; and è vera perché $40\le57\le60$ ; $57\%2=1$ , quindi dispari. Per confrontare vettori si usano & e | con le parentesi.
Parti 7-8, cicli. range(1, n+1) include $n$ ed esclude il finale. Cinque iterazioni di $+10$ danno $50$ . Il while richiede di incrementare il contatore a mano.
Parte 9, funzione radar. $d=\frac{c\,\Delta t}2$ : con $\Delta t=0{,}5\ \mu$ s la distanza è $75$ m. Un $\Delta t$ negativo solleva ValueError.
Parti 10-13, vettori e matrici. Indici da 0: v[1]=13, v[-1]=23. reshape riempie per righe; m[1,3]=12 e m.shape=(3,5). v.T non trasforma un array 1D; v.reshape(-1,1) dà la colonna.
Parti 14-17, prodotti. v1 @ v2 con $v_1=(10,20,30)$ e $v_2=(1,2,3)^T$ dà $140$ ; v2 @ v1 è $3\times3$ . * è elemento per elemento: $(300,200,600)$ . A @ B con due righe $1\times3$ dà errore di dimensioni. x.reshape(1,3)*y.reshape(3,1) usa il broadcasting e dà una matrice $3\times3$ .
Parte 18, matrici speciali. np.zeros((2,3)) vuole una tupla. np.linspace(1,3,30) dà $30$ punti con estremi inclusi. np.arange(15,17,0.2) ha $10$ elementi perché il $17$ è escluso; con il numero di campioni se ne contano $11$ .
Parti 19-20, grafici. Con $f=5$ Hz e $t_s=1$ ms su $[0,1]$ ci sono $1001$ campioni. plt.plot(t,y,t,w) disegna due curve nello stesso riquadro; '--' e 'r' hanno lo stesso significato che in MATLAB.
Parte 21, subplot. plt.subplot(2,3,k) numera da 1, per righe. np.sinc(0{,}5)=\frac2\pi=0{,}6366. Per la tangente si limita l'asse con plt.ylim(-10,10).
Parte 22, derivata numerica. $f'(t_n)\approx\frac{f_{n+1}-f_n}{\Delta t}$ , in tre scritture che coincidono (True True). Con $\Delta t=10^{-3}$ l'errore massimo è $0{,}01974$ , cioè $\frac{\Delta t}2\lvert f''\rvert=19{,}74\,\Delta t$ : proporzionale al passo. Confrontandolo con la derivata nel punto di mezzo, l'errore scende a $10^{-5}$ .
Operazioni elementari. Con $a=3$ e $b=4$ : $a+b=7$ , $a-b=-1$ , $ab=12$ , $a/b=0{,}75$ , $a^b=81$ , $\sqrt b=2$ .

python

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)

Modulo e fase. $\lvert z\rvert=\sqrt{\mathrm{Re}(z)^2+\mathrm{Im}(z)^2}$ e $\arg(-3)=\pi$ , perché $-3$ sta sul semiasse reale negativo. Risultati: (4+5j), (-1+7j), (0.28-0.04j), modulo 5.0, fase 3.141592653589793, coniugato (1-1j).

python

g = np.sqrt(-1 + 0j)      # np.sqrt(-1) darebbe nan: l'argomento deve essere complesso
rapporto = a / b
modulo_b = np.abs(b)         # MATLAB: abs(b)
fase_c = np.angle(c)         # MATLAB: angle(c)

Confronto e resto. $57=2\cdot28+1$ , quindi $57\,\%\,2=1$ e il numero è dispari. Il confronto doppio 40 <= a <= 60 è vero perché $40\le57\le60$ .

python

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

python

if a % 2 == 0:
    print("a è pari")

Ciclo. Dopo $n$ passi $s=10n$ : con $n=5$ si ha $s=50$ . Il while non ha un contatore automatico, quindi i = i + 1 è obbligatorio.

python

for i in range(1, n_iter + 1):
    s = s + 10
    print(f"Iterazione n° {i}, s = {s:2.1f}")

python

while i <= n_iter:
    s = s + 10
    print(f"Iterazione n° {i}, s = {s:2.1f}")
    i = i + 1

Radar. $2d=c\,\Delta t$ , quindi $d=\frac{c\,\Delta t}2$ ; con $\Delta t=0{,}5\ \mu$ s si ottiene $75$ m.

python

c = 3e8                       # velocita' della luce [m/s]
d = (delta_t * c) / 2         # l'impulso fa andata e ritorno: meta' del percorso

Vettori e matrici. v[-1] è $23$ e v.reshape(-1, 1) dà la colonna $5\times1$ . Per la matrice $3\times5$ , m[1, 3] è la riga di indice $1$ e la colonna di indice $3$ , cioè $12$ .

python

print(v[-1])         # MATLAB v(end): l'ultimo
print(v.reshape(-1, 1))          # per ottenere la colonna (5x1) serve reshape (o v[:, None])

Prodotto riga per colonna. $v_1\,v_2=\sum_k(v_1)_k(v_2)_k=10\cdot1+20\cdot2+30\cdot3=140$ ; la matrice $v_2\,v_1$ ha in posizione $(r,c)$ il prodotto $(v_2)_r(v_1)_c$ .

python

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

Matrici speciali. np.linspace(1, 3, 30) ha passo $\frac{2}{29}\approx0{,}0690$ e il primo elemento $1$ ; np.zeros((2, 3)) è una matrice $2\times3$ di zeri.

python

matrice_zeri = np.zeros((2, 3))         # il formato e' UNA tupla: (righe, colonne)
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)

Campioni. $\frac{1-0}{0{,}001}+1=1001$ campioni; con passo $0{,}05$ s un periodo di $0{,}2$ s contiene solo $4$ campioni, e la forma d'onda si perde.

python

t = t_inizio + np.arange(round((t_fine - t_inizio) / ts) + 1) * ts
y = np.sin(2 * np.pi * f * t)

Derivata, tre scritture. Il ciclo, np.diff e le fette danno lo stesso risultato. Nelle fette f[1:] sono i successivi e f[:-1] i precedenti.

python

fd[n] = (f[n + 1] - f[n]) / (t[n + 1] - t[n])
fd[-1] = fd[-2]                           # all'ultimo istante manca il campione successivo

python

fd_vec = (f[1:] - f[:-1]) / (t[1:] - t[:-1])
fd_vec = np.append(fd_vec, fd_vec[-1])

Errore in funzione del passo. Con $\Delta t=0{,}1$ , $0{,}01$ , $0{,}001$ l'errore massimo è $1{,}941611$ , $0{,}197327$ , $0{,}019739$ : il rapporto errore/passo resta circa $19{,}7$ , come previsto da $\frac{\Delta t}2\,(2\pi)^2$ .

Errori tipici.

Parti 1-2: a ^ b in Python è l'XOR bit a bit ( $3\oplus4=7$ ), non la potenza; la potenza è a ** b.
Parte 3: np.sqrt(-1) dà nan; l'argomento va scritto -1 + 0j.
Parti 4-8: range(1, n) esclude $n$ ; while senza i = i + 1 non termina.
Parte 9: dimenticare la divisione per $2$ dell'andata e ritorno raddoppia la distanza ( $150$ m).
Parti 10-17: v[1] è il secondo elemento, non il primo; * al posto di @ non fa il prodotto riga per colonna.
Parte 18: np.arange(15, 17, 0.2) esclude il $17$ e può sbagliare di un elemento.
Parti 19-22: senza plt.show() la figura non si vede; fd[-1] lasciato a zero fa crollare il grafico della derivata all'ultimo istante.

Esercizio - laboratorio 1, Python al posto di MATLAB

Cosa si vuole fare e perché

Prima di cominciare: le differenze MATLAB → Python che servono

Parti 1-2 - variabili e stampa

Parte 3 - numeri complessi

Parti 4-6 - if, operatori logici, resto

Parti 7-8 - cicli for e while

Parte 9 - una funzione: il radar

Parti 10-13 - vettori e matrici

Parti 14-17 - prodotti tra vettori

Parte 18 - matrici e vettori speciali

Parti 19-20 - grafici e curve sovrapposte

Parte 21 - sei grafici in una figura

Parte 22 - la derivata numerica

Controllo

Teoria collegata

Parti 4-6 - `if`, operatori logici, resto

Parti 7-8 - cicli `for` e `while`