Analisi di Fourier in Python

Il principio su cui si basa l'analisi di Fourier è la constatazione del fatto che spesso in natura un segnale oscillante si presenta non come una oscillazione ad una singola frequenza (o lunghezza d'onda) "pura", ma come una sovrapposizione di frequenze (o lunghezze d'onda) differenti. Esempi classici si trovano nell'ottica (la luce bianca è una sovrapposizione di differenti colori, cioè frequenze, che possono essere decomposte facendo passare la luce attraverso un prisma), in musica (i suoni "armonici" non sono altro che la sovrapposizione di una frequenza fondamentale del suono, con sovrapposte tutte le frequenze multiple della fondamentale), nelle comunicazioni (un segnale radio, o TV, è costituito da un segnale "portante" ad alta frequenza, che serve a trasportare il segnale a grandi distanze, con sovrapposta una modulazione a bassa frequenza, che è il segnale vero e proprio, da filtrare dal resto per poterlo estrarre), ecc.
Lo scopo dell'analisi di Fourier è di dare un metodo per separare le varie frequenze contenute in un segnale e analizzare qual'è il contributo delle singole frequenze alla ricostruzione del segnale di partenza.

Analisi di Fourier

Supponiamo di avere un segnale (che supponiamo variare nel tempo, ma le stesse considerazioni si possono fare se il segnale varia nello spazio!) rappresentato da una data funzione periodica \(f(t)\).

Teorema di Fourier:


Data una funzione \(f(t)\) a quadrato sommabile, periodica su un periodo \(T\), si può sempre scrivere che: \[ f(t) = a_0 + \sum_{n=1}^{+\infty} a_n \cos\left( \frac{2\pi}{T} nt \right) + b_n \sin\left( \frac{2\pi}{T} nt \right) \] Notiamo che le funzioni \(\cos\left( \frac{2\pi}{T} nt \right)\) e \(\sin\left( \frac{2\pi}{T} nt \right)\) sono funzioni che hanno presentano, nell'intervallo \(t \in [0,T]\), \(n\) oscillazioni!
Esempio (per \(T=2\pi\)):

In [3]:
import numpy as np
import matplotlib.pyplot as plt
x=np.linspace(0.0,2.0*np.pi,129)
%matplotlib inline
plt.subplot(221)
n=1
cx1=np.cos(n*x)
sx1=np.sin(n*x)
plt.title("n=1")
plt.plot(x,cx1)
plt.plot(x,sx1)

plt.subplot(222)
n=2
cx2=np.cos(n*x)
sx2=np.sin(n*x)
plt.title("n=2")
plt.plot(x,cx2)
plt.plot(x,sx2)

plt.subplot(223)
n=3
cx3=np.cos(n*x)
sx3=np.sin(n*x)
plt.title("n=3")
plt.plot(x,cx3)
plt.plot(x,sx3)

plt.subplot(224)
n=4
cx4=np.cos(n*x)
sx4=np.sin(n*x)
plt.title("n=4")
plt.plot(x,cx4)
plt.plot(x,sx4)

plt.show()
In sostanza, il teorema di Fourier stabilisce che ogni funzione periodica può scriversi come una sovrapposizione di una parte media (\(a_0\)), più onde con frequenza di oscillazione che sono multiple intere della frequenza fondamentale: \[ \omega = \frac{2\pi}{T} \] I coefficienti \(a_n\) e \(b_n\), di conseguenza, rappresentano il contributo che ciascuna delle frequenze (multiple di \(\omega\)) dà alla costruzione del segnale!

Sviluppo di Fourier complesso

Nella pratica, invece di dover trattare due serie di coefficienti reali \(a_n\) e \(b_n\), è più conveniente utilizzare una sola sequenza di coefficienti complessi.
Utilizzando le relazioni di Eulero:
\[\cos x = \frac{e^{ix} + e^{-ix}}{2}; \qquad \sin x = \frac{e^{ix} - e^{-ix}}{2i}\]
possiamo riscrivere la \(f(t)\) come:
\[f(t) = a_0 + \sum_{n=1}^{+\infty} a_n \left[\frac{e^{i\omega nt} + e^{-i\omega nt}}{2}\right] + b_n \left[\frac{e^{i\omega nt} - e^{-i\omega nt}}{2i}\right]\]
da cui otteniamo:
\[f(t) = a_0 + \sum_{n=1}^{+\infty} \left( \frac{a_n - i b_n}{2} \right) e^{i\omega n t} + \sum_{n=1}^{+\infty} \left( \frac{a_n + i b_n}{2} \right) e^{-i\omega n t}\]
Definendo le quantià:
\[ c_0 = a_0; \qquad c_n = \frac{a_n - i b_n}{2}; \qquad c_{-n} = \frac{a_n + i b_n}{2}, \qquad n=1,\ldots,+\infty \]
e notando che: \(e^{i\omega n t} = 1\) per \(n=0\), si ha che \(f(t)\) può quindi scriversi come:
\[ f(t) = \sum_{n=-\infty}^{+\infty} c_n e^{i\omega n t} \qquad \textrm{(1)} \]
Cioé, abbiamo ottenuto che lo sviluppo complesso di Fourier corrisponde ad una sovrapposizione di onde con frequenze sia positive che negative.

Condizione di realtà

Notiamo che, dalle relazioni che abbiamo scritto si ha:
\[ c_{-n} = c^*_{n} \]
cioé la parte negativa dei coefficienti dello sviluppo di Fourier coincide con i complessi coniugati della parte positiva.
Inoltre si può notare che \(c_0\) è una quantità reale!

Calcolo dei coefficienti di Fourier

Abbiamo quindi visto come i coefficienti di Fourier rappresentino il contributo che ciascuna frequenza dà alla costruzione di un segnale periodico. Come calcolarli?
Possiamo moltiplicare la (1) a sinistra e a destra per un fattore: \(e^{-i\omega mt}\) e integrare sul periodo di ripetizione \(T\) del segnale:
\[ \int_0^T f(t) e^{-i \omega mt} dt = \int_0^T \sum_{n=-\infty}^{+\infty} c_n e^{i\omega nt} e^{-i\omega mt} dt = \sum_{n=-\infty}^{+\infty} c_n \int_0^T e^{i(n-m)\omega t} dt \qquad \textrm{(2)}\]
Notiamo ora che, se \(n \neq m\):
\[ \int_0^T e^{i(n-m)\omega t} dt = \frac{1}{i(n-m)\omega} \int_0^{i(n-m)\omega T} e^z dz = \frac{1}{i(n-m)\omega} \left[ e^{i(n-m)\omega T}-1 \right] = 0 \]
avendo effettuato la sostituzione:
\[ z = i(n-m)\omega t; \qquad dz = i(n-m)\omega dt\]
e ricordando la definizione di \(\omega\), nonché il fatto che \(n\) ed \(m\) sono interi.
Invece, nel caso: \(n = m\), si ha:
\[ \int_0^T e^{i(n-m)\omega t} dt = \int_0^T dt = T \]
da cui concludiamo che:
\[ \int_0^T e^{i(n-m)\omega t} dt = T \delta_{nm} \]
La (2), quindi diventa:
\[ \int_0^T f(t) e^{-i \omega mt} dt = \sum_{n=-\infty}^{+\infty} c_n \int_0^T e^{i(n-m)\omega t} dt = \sum_{n=-\infty}^{+\infty} c_n T \delta_{nm} = T c_m \]
Cioè, ogni coefficiente di Fourier \(c_n\) può essere calcolato moltiplicando la \(f(t)\) per il fattore \(e^{-i\omega mt}\), integrando sull'intervallo di periodicità \(T\) e dividendo per lo stesso:
\[ c_n = \frac{1}{T} \int_0^T f(t) e^{-i\omega nt} dt \]

Teorema di Parceval e spettro di energia

In realtà, più che andare a calcolare il contributo in ampiezza del singolo coefficiente di Fourier \(c_n\) al segnale totale, è più conveniente andare a calcolarsi il suo contributo all'energia totale del sistema.
Il Teorema di Parceval serve proprio a questo scopo e afferma che:
\[ \int_0^T |f(t)|^2 dt = T \left( c_0^2 + 2 \sum_{n=1}^{+\infty} |c_n|^2 \right) \]
cioé l'energia totale del segnale può scriversi come la somma del coefficiente \(c_0\) al quadrato più il doppio della somma dei moduli quadri dei singoli coefficienti dello sviluppo di Fourier. In altre parole, ogni multiplo \(n\) della frequenza fondamentale (sia nella sua parte positiva che negativa) contribuisce tramite il modulo quadro del coefficiente \(c_n\) all'energia totale del sistema.
Dimostrazione
Sostituendo alla \(f(t)\) il suo sviluppo di Fourier:
\[ \int_0^T |f(t)|^2 dt = \int_0^T f(t)\,f^*(t) dt = \int_0^T \left( \sum_{n=-\infty}^{+\infty} c_n e^{-i\omega nt} \right) \left( \sum_{m=-\infty}^{+\infty} c^*_m e^{+i\omega mt} \right) dt = \int_0^T \left( \sum_{n=-\infty}^{+\infty} \sum_{m=-\infty}^{+\infty} c_n c^*_m e^{i(m-n)\omega t} \right) dt = \\ = \sum_{n=-\infty}^{+\infty} \sum_{m=-\infty}^{+\infty} c_n c^*_m \int_0^T e^{i(m-n)\omega t} dt\]
e ricordando che:
\[ \int_0^T e^{i(m-n)\omega t} dt = T \delta_{nm} \]
otteniamo:
\[ \int_0^T |f(t)|^2 dt = \sum_{n=-\infty}^{+\infty} \sum_{m=-\infty}^{+\infty} c_n c^*_m T \delta_{nm} = T \sum_{n=-\infty}^{+\infty} c_n c^*_n = T \sum_{n=-\infty}^{+\infty} |c_n|^2 \]
Possiamo ora scomporre la sommatoria in tre pezzi distinti: la parte con \(n<0\), \(n=0\) ed \(n>0\):
\[ \int_0^T |f(t)|^2 dt = T \left( \sum_{n=-\infty}^{-1} |c_n|^2 + |c_0|^2 + \sum_{n=1}^{+\infty} |c_n|^2\right) \]
e ricordando che: \(c_{-n} = c^*_n\), e che \(c_0\) è reale, si ha:
\[ \int_0^T |f(t)|^2 dt = T \left( c_0^2 + 2 \sum_{n=1}^{+\infty} |c_n|^2\right) \]
che è quello che si voleva dimostrare.
Generalmente si chiama spettro di energia del segnale, il grafico di una funzione \(E(n)\), definita come:
\[ E(n) = \begin{cases} c_0^2 & \textrm{per}\ n=0;\\ 2|c_n|^2 & \textrm{per}\ n=1, \ldots, +\infty \end{cases} \]
Quindi, la funzione \(E(n)\) ci dà proprio l'informazione che stavamo cercando, cioé il contributo che ogni singolo multiplo \(n\) della frequenza fondamentale \(\omega\) dà all'energia totale del sistema

Sviluppo di Fourier discreto

Quando ci si trova ad analizzare dati sperimentali, tuttavia, non si dispone dei valori del segnale \(\forall t\), ma in genere il segnale viene campionato su un insieme di punti discreto: \(t_j = j\Delta t\), con \(j=0,\ldots, N\) e \(\Delta t = T / N\). Il Teorema del campionamento (o di Nyquist-Shannon) afferma che, dato lo sviluppo di Fourier di una certa \(f(t)\) come:
\[ f(t) = \sum_{k=-\infty}^{+\infty} c_k e^{ik\omega t} \]
non tutti i valori di \(k\) (sia positivi che negativi) sono indipendenti tra di loro in un segnale campionato, ma solo quelli per cui \(-N_{max} \leq k \leq N_{max}\), dove \(N_{max} = N/2\), che viene detto Numero di Nyquist e la corrispondente frequenza: \(\omega_{max} = N_{max} \omega\), detta Frequenza di Nyquist del segnale campionato.
Qual'è il significato del teorema? Vuol dire che se io considero un segnale che contenga nel suo sviluppo di Fourier un valore di \(k\) pari a: \(k = N_{max} + m\), il contributo di quel \(k\) al segnale è identico a quello che darebbe un valore di \(k\) pari ad \(m-N_{max}\). Infatti, sui punti discreti \(t_j\):
\[ e^{ik\omega t_j} = e^{i(N_{max}+m)\omega t_j} = e^{iN_{max}\omega t_j} e^{im\omega t_j}\]
e ricordando che: \(\omega = \frac{2\pi}{T}\) e che: \(t_j = j \Delta t = \frac{jT}{N}\), si ha:
\[ e^{ik\omega t_j} = e^{i\frac{N}{2}\omega t_j} e^{im\omega t_j} = e^{i(N-\frac{N}{2})\frac{2\pi}{T}\frac{jT}{N}} e^{im\omega t_j} = e^{i2\pi j} e^{-i\frac{N}{2}\omega t_j} e^{im\omega t_j} = e^{i(m-N_{max})\omega t_j}\]
quindi, l'onda con frequenza \(k\omega\) passa per gli stessi punti di un'onda con frequenza (negativa!) \((m-N_{max})\omega\).
Esempio: nell'esempio seguente, si scelgono: \(\omega=1\) (\(T=2\pi\)), \(N=32\) punti di campionamento e quindi \(N_{max}=16\). Si noti come per \(m=12\) la funzione \(\cos(k\omega t)\) (linea blu), con \(k=m+N_{max}\) e la funzione \(\cos(k_a\omega t)\) (linea verde), con \(k_a=m-N_{max}\), quando vengono rappresentate su \(N=32\) punti discreti, si intersecano negli stessi punti (rappresentati dai pallini gialli per la prima funzione e dai triangoli blu per la seconda)...

In [4]:
T = 2.0 * np.pi
N = 32
NN = 1000
Nmax = N / 2
m = 12
k = m + Nmax
ka = m - Nmax
x = np.linspace( 0.0, T, N + 1 )
xx = np.linspace( 0.0, T, NN + 1 )
yk = np.cos(k*x)
ykk = np.cos(k*xx)
yka = np.zeros(N+1)
ykka = np.zeros(NN+1)
for j in range(0,N+1):
    yka[j] = np.cos(ka*x[j])
for j in range(0,NN+1):
    ykka[j] = np.cos(ka*xx[j])
plt.plot(x,yk,'yo',markersize=20)
plt.plot(xx,ykk,'-')
plt.plot(x,yka,'b^',markersize=10)
plt.plot(xx,ykka,'-')
plt.show()
Lo stesso accade naturalmente per una frequenza \(k = - N_{max} - m\), che risulta indistinguibile da una frequenza (positiva!) \(+ N_{max} - m\). Questo problema, noto sotto il nome di aliasing, rende impossibile distinguere le frequenze al di sopra della frequenza di Nyquist \(\omega_{max}\) e quelle al di sotto di \(-\omega_{max}\) da quelle che stanno all'interno dell'intervallo stesso. Di conseguenza, ha senso considerare nello sviluppo di Fourier discreto di una funzione, solo i valori di \(k \in [-N_{max}, +N_{max}]\):
\[ f(t_j) = \sum_{k=-N_{max}}^{+N_{max}} c_k e^{ik\omega t_j} \]
e, analogamente, il teorema di Parceval ci dirà che:
\[ \sum_{j=0}^{N} |f(t_j)|^2 = T \left( c_0^2 + 2 \sum_{k=1}^{N_{max}} |c_k|^2 \right) \]
Allo stesso modo, per calcolare i coefficienti di Fourier dovremo trasformare l'integrale in una somma discreta, ad esempio utilizzando il metodo dei trapezi, che ci dà:
\[ c_k = \frac{1}{T} \frac{h}{2} \left[ F_k(t_0)+F_k(t_N) + 2\sum_{j=1}^{N-1} F_k(t_j) \right] = \frac{h}{T} \sum_{j=0}^{N-1} F_k(t_j) \]
dove \(F_k(t_j) = f(t_j) e^{ik\omega t_j}\), e \(F_k(t_N) = F_k(t_0)\) per la periodicità della funzione.

Trasformata di Fourier discreta (dft) e Trasformata di Fourier veloce (fft)

Notiamo come il calcolo dei coefficienti \(c_k\) debba essere effettuato solo per i valori \(k>0\), essendo valida la proprietà: \(c_{-k} = c^*_k\). Questo equivale a calcolare \(c_k\) per \(k=0,\ldots,N/2\), cioé \(N/2+1\) valori complessi. In realtà, poiché \(c_0\) è reale risulta inutile calcolarne la parte immaginaria nulla. Si può inoltre far vedere che anche \(c_{N/2}\) è reale, poiché:
\[ c_{N/2} = \frac{h}{T} \sum_{j=0}^{N-1} F_k(t_j) = \frac{h}{T} \sum_{j=0}^{N-1} f(t_j) e^{i\frac{N}{2}\frac{2\pi}{T}\frac{jT}{N}} = \frac{h}{T} \sum_{j=0}^{N-1} f(t_j) e^{i\pi j} = \frac{h}{T} \sum_{j=0}^{N-1} (-1)^j f(t_j) \]
Di conseguenza, i valori effettivi da calcolare sono: \(2(N/2+1) - 2 = N\), in quanto le parti immaginarie (nulle!) di \(c_0\) e \(c_{N/2}\) non necessitano di essere calcolate. Siccome il calcolo richiede quindi la valutazione per \(N/2\) volte della funzione complessa: \(F_k(t_j),\ k=0,\ldots,N/2,\ j=0,\ldots,N-1\), in totale dovremo effettuare \(N^2\) prodotti per valutare i coefficienti indipendenti. Esistono algoritmi chiamati FFT (Fast Fourier Transform), che permettono di effettuare tale valutazione effettuando circa \(N\log_2 N\) prodotti, che per \(N\) grandi fa una enorme differenza!

Python fornisce svariate librerie per il calcolo delle FFT.
La libreria standard numpi.fft fornisce le seguenti funzioni: 
rfft(A), irfft(A), rfft2(A), irfft2(A), rfftn(A), irfftn(A)
che calcolano la trasformata di un vettore reale A 1, 2, ..., n - dimensionale, e le analoghe funzioni per vettori complessi: 
fft(A), ifft(A), fft2(A), ifft2(A), fftn(A), ifftn(A)
Esistono poi librerie specifiche (es. pyFFTW) che forniscono funzioni analoghe ma con performances di calcolo migliori.
Note:
  • la funzione ak=rfft(a) restituisce un vettore di numeri complessi di dimensione N/2+1, contenente solo le frequenze positive del vettore a, in quanto quelle negative sono le complesse coniugate di quelle positive;
  • i valori ak[0] e ak[N/2] sono reali per quanto visto sopra;
  • la trasformata diretta rfft è non-normalizzata, cioé gli ak sono moltiplicati per il numero N di punti discreti, mentre la trasformata inversa irfft è normalizzata (cioé i valori reali sono quelli giusti!)


Esempio: trasformata della funzione \(f(t) = \sin(x) + \cos(5x)\) su N=64 punti discreti...

In [5]:
N = 64                              # Numero di intervalli di discretizzazione
dpi = 2.0 * np.pi                   # 2 pi greco
dx = dpi / N                        # Spaziatura degli intervalli discretizzati
x = np.arange( 0.0, dpi, dx )       # Genera un vettore di 128 punti con valori tra 0 e 2 pi - dx
y = np.sin(x) + np.cos(5.0*x)       # Funzione di cui calcolare la trasformata
plt.xlabel(r"$t$")
plt.ylabel(r"$f(t)$")
plt.title(r"$f(t)=\sin(x)+\cos(5x)$")
plt.plot(x,y,'-')
plt.show()
In [6]:
ak = np.fft.rfft(y)
print(ak)
print(len(ak))

[ -2.16493490e-15 +0.00000000e+00j  -1.97355722e-14 -3.20000000e+01j
  -3.21460981e-15 -2.61516588e-15j  -3.40823918e-15 -6.22546217e-15j
  -3.34221411e-15 +1.47424471e-15j   3.20000000e+01 -4.78770100e-14j
   7.01669845e-15 +6.30487299e-16j   5.95648052e-15 -3.83253469e-15j
   5.20103039e-15 -1.41229022e-16j   6.15547725e-15 +2.12716763e-15j
   6.11863178e-15 +2.83590513e-16j  -9.43689571e-15 -1.52791457e-15j
   6.35420350e-16 +2.86618156e-15j   3.00579718e-15 +3.23151528e-16j
  -8.91588444e-16 -1.58261755e-16j  -1.30586997e-15 -1.99840144e-15j
  -2.05391260e-15 -8.88178420e-16j  -3.08222681e-15 -3.77475828e-15j
   3.61146158e-15 -6.09361102e-15j   1.40439217e-14 -1.47339306e-15j
  -1.53013515e-15 +7.41155958e-15j   5.88418203e-15 +3.30427141e-15j
   2.60640171e-15 -6.85996273e-15j   8.51481559e-15 +9.63295734e-15j
  -9.08681098e-15 +2.96739545e-15j   5.91298858e-15 +1.89689818e-15j
  -3.23565031e-17 +2.32870696e-15j   1.77635684e-15 +2.35786601e-15j
  -5.08894449e-15 +4.46531050e-15j  -1.00887660e-14 -3.47552402e-15j
   5.21346490e-15 -3.26145592e-15j   1.13662069e-15 +3.55271368e-15j
  -4.60742555e-15 +0.00000000e+00j]
33

Scriviamo ora uno script che calcoli lo spettro di una data funzione nota su N punti:

In [7]:
def spettro(f):
    N = len(f)                       # Lunghezza del vettore di partenza
    Nf = int( N / 2 )                # Frequenza di Nyquist
    om = np.linspace(0.0, Nf, Nf+1)  # Vettore di output contenente le frequenze
    sp = np.zeros(Nf+1)              # Vettore di output contenente lo spettro
    ak = np.fft.rfft(f) / N          # Trasformata del vettore di partenza (normalizzata!)
    for n in range(Nf+1):
        if n == 0:
            sp[n] = np.square(np.abs(ak[n]))
        else:
            sp[n] = 2.0 * np.square(np.abs(ak[n]))
    return om, sp

w, Ew = spettro(y)
plt.title("Spettro")
plt.xlabel(r"$\omega$")
plt.ylabel(r"$E(w)$")
plt.plot(w,Ew,'-')
plt.show()

Esempio:
generiamo una funzione che abbia uno spettro con tutte le armoniche e ricostruiamo il segnale di partenza usando solo le prime \(n \le N/2\) armoniche del segnale, per vedere come le varie armoniche contribuiscono al segnale stesso.
Scegliamo come funzione di cui calcolare lo spettro l'impulso triangolare:
\[ f(x) = \begin{cases} x & \textrm{per} \qquad 0 \le x < \pi; \\ 0 & \textrm{per} \qquad x = \pi; \\ x-2\pi & \textrm{per} \qquad \pi < x < 2\pi \\ \end{cases} \]

In [8]:
N2 = int( N / 2 )
y[0:N2] = x[0:N2]
y[N2] = 0.0
y[N2+1:] = x[N2+1:] - dpi
plt.title("Impulso triangolare")
plt.xlabel(r"$t$")
plt.ylabel(r"$f(t)$")
plt.plot(x,y,'-')
plt.show()
In [9]:
w, Ew = spettro(y)
plt.title("Spettro")
plt.xlabel(r"$\omega$")
plt.ylabel(r"$E(\omega)$")
plt.xscale("log")                  # Spettro in scala log-log
plt.yscale("log")
wteo = w[1:N2]**(-2.0)
plt.ylim(1.e-5,10)
plt.plot( w, Ew, '-', w[1:N2], wteo, 'g-' )
plt.show()

Definiamo quindi una funzione che effettua la ricostruzione del segnale (trasformata inversa) con un numero di armoniche \(n \le N/2\):

In [10]:
def reconstruct( ak, n ):
    Nyq = len( ak ) - 1         # Frequenza di Nyquist
    Ns = 2 * Nyq                # Numero di punti nel vettore di output
    sr = np.zeros( Ns )         # Vettore di output
    if ( n > Nyq ):
        print("Errore! Il numero di armoniche da utilizzare per la ricostruzione deve essere minore della frequenza di")
        print("Nyquist del segnale!")
        return sr
    else:
        akk = np.zeros(Nyq+1, dtype=complex)
        akk[:n+1] = ak[:n+1]
        sr = np.fft.irfft(akk)
    return sr

ak = np.fft.rfft(y)              # Calcolo dei coefficienti di Fourier
sr3 = reconstruct(ak,3)
plt.title("Segnale di partenza e ricostruito con 3 armoniche")
plt.xlabel(r"$t$")
plt.ylabel(r"$f(t), f_3(t)$")
plt.plot(x,y,'b-',x,sr3,'g-')
plt.show()
In [11]:
sr8 = reconstruct(ak,8)
plt.title("Segnale di partenza e ricostruito con 8 armoniche")
plt.xlabel(r"$t$")
plt.ylabel(r"$f(t), f_8(t)$")
plt.plot(x,y,'b-',x,sr8,'g-')
plt.show()
In [12]:
sr32 = reconstruct(ak,32)
plt.title("Segnale di partenza e ricostruito con 32 armoniche")
plt.xlabel(r"$t$")
plt.ylabel(r"$f(t), f_{32}(t)$")
plt.plot(x,y,'b-',x,sr32,'g-')
plt.show()