Esercizio - laboratorio 4, coefficienti di Fourier con la FFT

Testo (laboratorio 4 del corso Teoria dei Segnali, UniPD, file lab4.m, parti 9-13). (9-10) Campionare un periodo ($T_p=1$ s) del segnale $s(t)=\cos(2\pi t)$ con $N=30$ campioni e calcolarne i coefficienti di Fourier con la FFT, con la scalatura $\frac{dt}{T_p}\,\text{fft}(s)$, l'asse delle frequenze $f=[-\frac N2,\dots,\frac N2-1]\,F$ e il riordino delle due metà del vettore. (11-13) Fare lo stesso per un'onda quadra (1 nella prima metà del periodo, 0 nella seconda) con $N=30$ e poi con $N=150$ campioni per periodo, confrontando con i coefficienti analitici $S_n=\frac12\operatorname{sinc}\frac n2$ (in modulo) e spiegando che cosa cambia con $N$.

Teoria usata: Serie di FourierUn segnale periodico di periodo $T_p$ si scrive come somma di esponenziali alle frequenze multiple della fondamentale $F=1/T_p$: $s(t)=\sum_nS_ne^{i2\pi nFt}$, con $S_n=\frac1{T_p}\int_{T_p}s(t)e^{-i2\pi nFt}dt$. Si basa sull'ortogonalità degli esponenziali su un periodo. $S_0$ è il valor medio; per segnali reali $S_{-n}=S_n^*$ e si passa alla forma con coseni e seni; vale il teorema di Parseval $P=\sum|S_n|^2$. La convoluzione ciclica diventa il prodotto $T_pX_nY_n$ dei coefficienti. Le somme troncate presentano il fenomeno di Gibbs vicino ai salti.Serie di Fourier →, Trasformata di Fourier discreta (DFT) e FFTUn segnale discreto periodico di periodo $NT$ ha una trasformata discreta e periodica, la DFT: $S(kF)=\sum_{n=0}^{N-1}T,s(nT)e^{-i2\pi kn/N}$, con $F=1/(NT)$, e $s(nT)=\sum_{k=0}^{N-1}F,S(kF)e^{i2\pi kn/N}$. Dipende da soli $N$ numeri, si calcola senza approssimazioni e con la FFT costa $N\log_2N$ invece di $N^2$. I coefficienti di Fourier del segnale periodico sono $S_k=F,S(kF)$. La DFT dà anche campioni della trasformata di un segnale discreto o continuo a durata limitata (con zero-padding), ma la scalatura $T$ e l'asse delle frequenze vanno gestiti con cura.Trasformata di Fourier discreta (DFT) e FFT →, 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 → (aliasing), 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 → (il $\operatorname{sinc}$).

Cosa si vuole calcolare e perché

Un segnale periodico $s(t)$ di periodo $T_p$ ha lo sviluppo in serie $s(t)=\sum_nS_n\,e^{i2\pi nFt}$, $F=\frac1{T_p}$, con coefficienti (Serie di FourierUn segnale periodico di periodo $T_p$ si scrive come somma di esponenziali alle frequenze multiple della fondamentale $F=1/T_p$: $s(t)=\sum_nS_ne^{i2\pi nFt}$, con $S_n=\frac1{T_p}\int_{T_p}s(t)e^{-i2\pi nFt}dt$. Si basa sull'ortogonalità degli esponenziali su un periodo. $S_0$ è il valor medio; per segnali reali $S_{-n}=S_n^*$ e si passa alla forma con coseni e seni; vale il teorema di Parseval $P=\sum|S_n|^2$. La convoluzione ciclica diventa il prodotto $T_pX_nY_n$ dei coefficienti. Le somme troncate presentano il fenomeno di Gibbs vicino ai salti.Serie di Fourier →)

$$S_n=\frac1{T_p}\int_{T_p}s(t)\,e^{-i2\pi nFt}\,dt.$$

Campionando un periodo con $N$ campioni equispaziati $t_m=m\,dt$, $dt=\frac{T_p}N$, l'integrale diventa una somma, $S_n\approx\frac1{T_p}\sum_{m=0}^{N-1}dt\;s(t_m)\,e^{-i2\pi nm/N}$, cioè

$$S_n\approx\frac{dt}{T_p}\,\text{fft}(s)[n]=\frac1N\,\text{fft}(s)[n]$$

(fft(s)[n] è proprio $\sum_ms_me^{-i2\pi nm/N}$). Il fattore $\frac{dt}{T_p}$ è lo stesso $\frac1N$ del laboratorio 3: la FFT dà direttamente i coefficienti di Fourier $S_n$, ognuno associato alla frequenza $nF$. Resta da capire l'asse delle frequenze (a quali $n$ corrispondono gli elementi del vettore) e quanto è fedele il risultato al coefficiente vero: la risposta sta nell'aliasing.

Parti 9-10 - un coseno

import numpy as np
import matplotlib.pyplot as plt

# PARTE 9 - un periodo di coseno campionato: N = 30 campioni, periodo Tp = 1 s
N = 30
Tp = 1

dt = Tp / N                        # passo di campionamento: 1/30 s
t = np.arange(N) * dt              # istanti 0, dt, ..., (N-1) dt  (l'istante Tp e' ESCLUSO: apparterrebbe al periodo dopo)

fcos = 1                           # frequenza del coseno: 1 Hz = 1/Tp
s = np.cos(2 * np.pi * fcos * t)

plt.figure()
plt.plot(t, s)
plt.xlabel('t [s]'); plt.grid(True)
plt.show()

print(f"dt = {dt:.4f} s, campioni per periodo = {N}, frequenza di campionamento Fp = {1 / dt:.0f} Hz")

Risultato:

dt = 0.0333 s, campioni per periodo = 30, frequenza di campionamento Fp = 30 Hz

Il coseno ha frequenza $1$ Hz $=\frac1{T_p}=F$: un periodo esatto nei $30$ campioni. Il vettore t va da $0$ a $(N-1)\,dt=0{,}9667$ s: l'istante $T_p$ è già il primo del periodo successivo e si lascia fuori (altrimenti il periodo avrebbe $N+1$ campioni). Ci sono $N=30$ campioni per periodo, cioè una frequenza di campionamento $F_p=\frac1{dt}=30$ Hz.

# PARTE 10 - coefficienti di Fourier con la fft
S = dt / Tp * np.fft.fft(s)             # dt/Tp = 1/N: coefficienti S_n = (1/N) somma di s e^(-i 2 pi n m / N)
F = 1 / Tp                              # passo in frequenza = 1/Tp = 1 Hz
f = np.arange(-N // 2, N // 2) * F      # -15, -14, ..., 14 Hz

S = np.concatenate((S[N // 2:], S[:N // 2]))     # riordino: [S(N/2+1:N), S(1:N/2)] in MATLAB; equivale a np.fft.fftshift(S)

plt.figure()
plt.stem(f, np.abs(S), basefmt=' ')
plt.xlabel('Frequenza [Hz]'); plt.ylabel('Modulo coeff. Fourier'); plt.ylim(-0.5, 1); plt.grid(True)
plt.show()

print("asse:", f[0], "...", f[-1], "Hz  (", len(f), "righe )")
print("righe con |S| > 1e-9:", [(float(ff), round(float(abs(c)), 6)) for ff, c in zip(f, S) if abs(c) > 1e-9])
print("valori complessi in -1 e +1 Hz:", np.round(S[f == -1][0], 6), np.round(S[f == 1][0], 6))
print("uguale a fftshift:", np.allclose(S, np.fft.fftshift(dt / Tp * np.fft.fft(s))))

Risultato:

asse: -15.0 ... 14.0 Hz  ( 30 righe )
righe con |S| > 1e-9: [(-1.0, 0.5), (1.0, 0.5)]
valori complessi in -1 e +1 Hz: (0.5+0j) (0.5-0j)
uguale a fftshift: True

I coefficienti. Si ha $\cos(2\pi t)=\frac12e^{i2\pi t}+\frac12e^{-i2\pi t}$ (formula di Eulero), quindi $S_{1}=S_{-1}=\frac12$ e tutti gli altri coefficienti sono nulli. È proprio ciò che stampa il codice: due sole righe, a $\pm1$ Hz, ciascuna di valore $0{,}5$ (e reale: la fase è nulla perché il coseno parte da un massimo). Una sinusoide reale ha sempre due righe, simmetriche (il segnale è reale: $S_{-n}=S_n^*$).

L'asse delle frequenze e il riordino. La fft restituisce $N$ numeri, indicati dall'indice $n=0,1,\dots,N-1$. Per la periodicità della DFT, l'indice $n$ e l'indice $n-N$ rappresentano la stessa riga: le righe da $n=\frac N2$ a $N-1$ sono in realtà quelle a frequenza negativa $-\frac N2,\dots,-1$. Per avere l'asse da $-15$ a $14$ Hz si mettono prima le posizioni da $N/2$ a $N-1$ e dopo quelle da $0$ a $N/2-1$: np.concatenate((S[N//2:], S[:N//2])), in MATLAB [S(N/2+1:N), S(1:N/2)]. È lo scambio delle due metà, cioè np.fft.fftshift (ultima riga stampata). L'asse np.arange(-N//2, N//2) * F ha $N$ punti, da $-15F$ a $14F$ (con $N$ pari c'è una riga in più a sinistra, $-15$ Hz $=-\frac{F_p}2$, che è il "doppione" di $+15$ Hz). Il massimo ordinabile è $\pm\frac{F_p}2=\pm15$ Hz: oltre non si vede niente, è il teorema 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 →).

Parti 11-13 - onda quadra, 30 campioni e 150 campioni

# PARTE 11-12 - onda quadra: 1 nella prima meta' del periodo, 0 nella seconda; N = 30
def onda_quadra(N, Tp=1):
    dt = Tp / N
    t = np.arange(N) * dt
    s = np.zeros(N)
    s[:N // 2] = 1                    # MATLAB s(1:N/2) = 1
    s[N // 2:] = 0                    # MATLAB s(N/2+1:N) = 0
    S = dt / Tp * np.fft.fft(s)
    S = np.concatenate((S[N // 2:], S[:N // 2]))         # riordino delle due meta'
    F = 1 / Tp
    f = np.arange(-N // 2, N // 2) * F
    return t, s, f, S

t30, s30, f30, S30 = onda_quadra(30)

plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(t30, s30); plt.ylim(-0.5, 1.5); plt.xlabel('t [s]'); plt.grid(True)
plt.subplot(1, 2, 2)
plt.stem(f30, np.abs(S30), basefmt=' '); plt.xlabel('Frequenza [Hz]'); plt.ylabel('Modulo coeff. Fourier'); plt.grid(True)
plt.show()

# analitico: S_n = (1/2) sinc(n/2) e^(-i pi n/2); in modulo (1/2) |sinc(n/2)|
def S_analitico(n):
    return 0.5 * np.sinc(n / 2) * np.exp(-1j * np.pi * n / 2)

print("   n   |S_n| con fft   (1/2)|sinc(n/2)|   differenza")
for n in [0, 1, 2, 3, 4, 5, 7, 9, 11, 14, -15]:
    num = abs(S30[f30 == n][0])
    an = abs(S_analitico(n))
    print(f"{n:>4}   {num:.5f}        {an:.5f}        {num - an:+.5f}")

Risultato per $N=30$:

   n   |S_n| con fft   (1/2)|sinc(n/2)|   differenza
   0   0.50000        0.50000        +0.00000
   1   0.31889        0.31831        +0.00058
   2   0.00000        0.00000        -0.00000
   3   0.10787        0.10610        +0.00177
   4   0.00000        0.00000        -0.00000
   5   0.06667        0.06366        +0.00300
   7   0.04982        0.04547        +0.00434
   9   0.04120        0.03537        +0.00583
  11   0.03649        0.02894        +0.00755
  14   0.00000        0.00000        -0.00000
 -15   0.03333        0.02122        +0.01211

L'onda quadra vale $1$ per $0\le t<\frac{T_p}2$ e $0$ per $\frac{T_p}2\le t<T_p$ (non è centrata nell'origine): i primi 15 campioni a 1, gli altri 15 a 0.

Il valore analitico. I coefficienti sono

$$S_n=\frac1{T_p}\int_0^{T_p/2}e^{-i2\pi nt/T_p}\,dt=\frac{1-e^{-i\pi n}}{i2\pi n}=\frac12\operatorname{sinc}\Bigl(\frac n2\Bigr)\,e^{-i\pi n/2},\qquad S_0=\frac12.$$

Il modulo è $\frac12\bigl|\operatorname{sinc}\frac n2\bigr|=\frac12\Bigl|\frac{\sin(\pi n/2)}{\pi n/2}\Bigr|$: vale $0$ per $n$ pari non nullo (lì $\sin\frac{\pi n}2=0$), $\frac1{\pi|n|}$ per $n$ dispari: $0{,}3183$ per $n=\pm1$, $0{,}1061$ per $\pm3$, $0{,}0637$ per $\pm5$, e così via come $\frac1n$. La fase $e^{-i\pi n/2}$ dipende solo da dove cade l'onda nel periodo (qui parte da $t=0$) e per questo si confrontano i moduli. Il valore medio $S_0=\frac12$ è la media dei campioni, $\frac{15}{30}$.

Cosa dà la FFT con $N=30$. La tabella mostra che i valori sono vicini a quelli veri ma sempre un po' più grandi: $0{,}3189$ contro $0{,}3183$ per $n=\pm1$ (errore $0{,}2\%$), $0{,}0412$ contro $0{,}0354$ per $n=\pm9$ (16%), e nell'ultima riga, $n=-15=-\frac N2$, $0{,}0333$ contro $0{,}0212$ (57%). Le righe pari sono zero (esattamente, come nel caso vero).

Perché c'è un errore: l'aliasing. L'onda quadra ha infinite armoniche, con $S_n\sim\frac1n$ che non si annullano mai: il suo spettro non è a banda limitata. Campionare con $N$ punti per periodo equivale a campionare nel tempo con frequenza $F_p=\frac N{T_p}$, e questo ripete lo spettro ogni $N$ righe e somma le code ripetute (formula del campionamento $S_c=\operatorname{rep}S$, 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 →):

$$S^{c}_n=\sum_{k=-\infty}^{\infty}S_{n+kN}.$$

Alla riga $n$ si aggiungono i contributi delle righe vere $n\pm N,\ n\pm2N,\dots$, tutte piccole ma non nulle. Più $n$ è vicino al bordo $\pm\frac N2$, più sono grandi: per $n\approx\frac N2$ la riga vera $n-N\approx-\frac N2$ ha indice quasi uguale in modulo a $n$, quindi ampiezza paragonabile. Il blocco finale della sezione lo verifica numericamente: se nei due punti di salto ($t=0$ e $t=\frac{T_p}2$) si mettono i valori di mezzo ($\frac12$, regola dell'emivalore) la fft coincide con la somma delle code fino a $10^{-6}$. Il laboratorio mette invece $1$ in $t=0$ e $0$ in $t=\frac{T_p}2$: questo aggiunge $\frac{1}{2N}\bigl(1-(-1)^n\bigr)$ ai coefficienti (cioè $\frac1N$ ai dispari), un altro errore dell'ordine di $\frac1N$.

# PARTE 12 (correzione) e 13 - stessa onda quadra con N = 150 campioni per periodo
t150, s150, f150, S150 = onda_quadra(150)

plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(t150, s150); plt.ylim(-0.5, 1.5); plt.xlabel('t [s]'); plt.grid(True)
plt.subplot(1, 2, 2)
plt.stem(f150, np.abs(S150), basefmt=' '); plt.xlabel('Frequenza [Hz]'); plt.ylabel('Modulo coeff. Fourier'); plt.grid(True)
plt.show()

print("asse con N = 150:", f150[0], "...", f150[-1], "Hz  (", len(f150), "righe )")
print("   n   N=30: fft  errore      N=150: fft  errore     (1/2)|sinc(n/2)|")
for n in [0, 1, 3, 5, 7, 9, 11]:
    a30 = abs(S30[f30 == n][0]); a150 = abs(S150[f150 == n][0]); an = abs(S_analitico(n))
    print(f"{n:>4}   {a30:.5f}  {a30 - an:+.5f}      {a150:.5f}  {a150 - an:+.5f}      {an:.5f}")
print("righe pari (n = 2, 4, 6) con N=30 e N=150:", np.abs(S30[[f30 == 2][0]]).round(12), np.abs(S150[[f150 == 2][0]]).round(12))
print(f"ultima riga (n = -N/2): N=30: {abs(S30[0]):.5f} contro {abs(S_analitico(-15)):.5f};  N=150: {abs(S150[0]):.5f} contro {abs(S_analitico(-75)):.5f}")

Risultato per $N=150$:

asse con N = 150: -75.0 ... 74.0 Hz  ( 150 righe )
   n   N=30: fft  errore      N=150: fft  errore     (1/2)|sinc(n/2)|
   0   0.50000  +0.00000      0.50000  +0.00000      0.50000
   1   0.31889  +0.00058      0.31833  +0.00002      0.31831
   3   0.10787  +0.00177      0.10617  +0.00007      0.10610
   5   0.06667  +0.00300      0.06378  +0.00012      0.06366
   7   0.04982  +0.00434      0.04564  +0.00016      0.04547
   9   0.04120  +0.00583      0.03558  +0.00021      0.03537
  11   0.03649  +0.00755      0.02919  +0.00026      0.02894
righe pari (n = 2, 4, 6) con N=30 e N=150: [0.] [0.]
ultima riga (n = -N/2): N=30: 0.03333 contro 0.02122;  N=150: 0.00667 contro 0.00424

Cosa cambia con $N=150$. Ci sono tre effetti dell'avere più campioni per periodo, a parità di $T_p=1$ s (stessa $F=1$ Hz):

1. L'asse delle frequenze si allarga. I punti sono $N=150$, sempre a distanza $F=1$ Hz, quindi l'asse va da $-75$ a $74$ Hz, cioè fino a $\pm\frac N{2T_p}=\pm\frac{F_p}2$: la frequenza di campionamento è cresciuta da $30$ a $150$ Hz e con essa il limite di Nyquist. Si vedono $150$ armoniche invece di $30$.
2. L'aliasing diminuisce. Le armoniche vere che si ripiegano sulla riga $n$ sono quelle a distanza $N$: più $N$ è grande, più sono lontane e più sono piccole ($\sim\frac1{\pi N}$). L'errore massimo sulle prime armoniche ($|n|\le9$) scende da $0{,}0058$ ($N=30$) a $0{,}0002$ ($N=150$): $0{,}31833$ contro $0{,}31831$ per $n=1$ (errore $2\cdot10^{-5}$), $0{,}03558$ contro $0{,}03537$ per $n=9$ (errore $2\cdot10^{-4}$).
3. Il difetto sul bordo resta, ma è più in alto e più piccolo. L'ultima riga ($n=-\frac N2=-75$) vale $0{,}0067$ contro $0{,}0042$ vero (sempre un errore del $57\%$ relativo, perché il confronto con le code è lo stesso) ma in valore assoluto è cinque volte più piccola: il modo giusto per fidarsi dei risultati è guardare le armoniche ben dentro all'asse, $|n|\ll\frac N2$.

# Errore massimo sulle prime armoniche, e perche' l'errore c'e'
for N, f_, S_ in [(30, f30, S30), (150, f150, S150)]:
    sel = np.abs(f_) <= 9
    err = np.max(np.abs(np.abs(S_[sel]) - np.abs(S_analitico(f_[sel]))))
    print(f"N = {N:>3}: errore massimo sulle righe con |n| <= 9: {err:.5f}")

# la fft di una forma d'onda campionata "ripete" lo spettro vero ogni N righe e somma le code (aliasing).
# Se nei due salti si mette il valore di meta' (0.5), il risultato e' esattamente la somma delle code:
N = 30
s_mezzo = np.zeros(N); s_mezzo[:N // 2] = 1; s_mezzo[0] = 0.5; s_mezzo[N // 2] = 0.5
S_mezzo = np.fft.fft(s_mezzo) / N
K = np.arange(-20000, 20001)
for n in [1, 3, 5]:
    somma_code = np.sum(S_analitico(n + K * N))        # somma su k dei coefficienti veri S_(n + kN)
    print(f"n = {n}: |fft| con valori di meta' = {abs(S_mezzo[n]):.6f}   |somma delle code| = {abs(somma_code):.6f}")

# il laboratorio mette 1 in t = 0 e 0 in t = Tp/2: aggiunge 0.5/N in t=0 e toglie 0.5/N in t=Tp/2
S_lab = np.fft.fft(s30) / N
corr = (0.5 / N) * (1 - np.exp(-1j * np.pi * np.arange(N)))
print("S_lab = S_mezzo + (0.5/N)(1 - (-1)^n):", np.allclose(S_lab, S_mezzo + corr))

N =  30: errore massimo sulle righe con |n| <= 9: 0.00583
N = 150: errore massimo sulle righe con |n| <= 9: 0.00021
n = 1: |fft| con valori di meta' = 0.317145   |somma delle code| = 0.317146
n = 3: |fft| con valori di meta' = 0.102589   |somma delle code| = 0.102590
n = 5: |fft| con valori di meta' = 0.057735   |somma delle code| = 0.057735
S_lab = S_mezzo + (0.5/N)(1 - (-1)^n): True

Le parti 12 (correzione) e 13 del file MATLAB contengono lo stesso codice (il blocco è duplicato). Nel file, nel grafico del tempo, ylabel dice 'Modulo coeff. Fourier' e nel grafico in frequenza xlabel dice 'Frequenza' senza unità: sono etichette scambiate (nella parte 12 correzione e nella 13), non influiscono sul calcolo.

Controllo

• Coseno: righe $\pm\frac12$ in $n=\pm1$, altre $<10^{-16}$.
• Onda quadra: $S_0=\frac12$; pari nulli; dispari $\frac1{\pi n}$ con errore da aliasing che passa da $0{,}0058$ ($N=30$) a $0{,}0002$ ($N=150$) per $|n|\le9$.
• La fft dei campioni con valori di mezzo ai salti coincide (a $10^{-6}$) con $\sum_kS_{n+kN}$.