Stima dei parametri di un segnale dai campioni - algoritmi nel tempo e in frequenza

Il programma del corso comprende "metodi e algoritmi per le misure: trasformata di Fourier, algoritmi per la stima di parametri nel dominio del tempo e della frequenza". Un oscilloscopio o una scheda di acquisizione restituiscono una sequenza di campioni $x[n]=x(nT_s)$, $n=0,\dots,N-1$, e ogni parametro "misurato" (frequenza, ampiezza, tempo di salitatempo per passare dal 10% al 90% dell'ampiezza di un fronte) è il risultato di un algoritmo applicato a quella sequenza. Questa nota raccoglie gli algoritmi di base, con un'implementazione in Python e i risultati su segnali simulati (verificati), per collegare le formule delle altre note ai calcoli che fa uno strumento.

Nel dominio del tempo

Statistiche del segnale. Valore medio $\bar x=\frac1N\sum x[n]$, valore efficace $\sqrt{\frac1N\sum x^2[n]}$ (RMS-DC) e RMS-AC con $x-\bar x$ (Valore efficace - RMS-DC, RMS-AC e segnali periodiciIl valore efficace è $V_{rms}=\sqrt{\frac1T\int_0^Tv^2dt}$. Due versioni: RMS-DC (tutto il segnale, componente continua compresa) e RMS-AC (si toglie il valor medio), legate da $V_{rms,DC}^2=V_{DC}^2+V_{rms,AC}^2$. Sinusoide $V_0/\sqrt2$; onda rettangolare $\sqrt{V_H^2d+V_L^2(1-d)}$, AC $(V_H-V_L)\sqrt{d(1-d)}$; triangolare/dente di sega AC $V_{pp}/(2\sqrt3)$, DC $\sqrt{(V_{max}^2+V_{max}V_{min}+V_{min}^2)/3}$. Il multimetro true-RMS in AC dà solo l'RMS-AC; l'oscilloscopio ha le funzioni DC RMS e AC RMS, da usare su un numero intero di periodi (N cycles).Valore efficace - RMS-DC, RMS-AC e segnali periodici →); picco-picco $\max x-\min x$. Il picco-picco è molto sensibile al rumore (un solo campione anomalo lo cambia): per questo si usano i percentilivalori sotto i quali cade una data percentuale dei campioni, meno sensibili ai valori anomali del massimo e del minimo (per esempio il $2\%$ e il $98\%$) o la media di più acquisizioni.

Attraversamento di un livello. Per misurare periodi e durate si cercano i campioni consecutivi che sono uno sotto e uno sopra il livello $L$, e si interpola linearmentesi assume che tra due campioni il segnale vari lungo una retta, per stimare l'istante esatto di un attraversamento per migliorare la risoluzione temporale oltre $T_s$: se $x[n]<L\le x[n+1]$, $$t_c=nT_s+\frac{L-x[n]}{x[n+1]-x[n]}\,T_s.$$

Frequenza. Con $M$ attraversamenti con pendenza positiva allo stesso livello (per esempio il valor medio) agli istanti $t_{c,1},\dots,t_{c,M}$: $$\hat f=\frac{M-1}{t_{c,M}-t_{c,1}}\qquad(\text{periodo medio su }M-1\text{ periodi}),$$ con il vantaggio che l'incertezza sull'istante si divide per il numero di periodi (la regola dei cursori, Misure con i cursori e incertezza dell'oscilloscopioCon i cursori si misurano intervalli di tempo e di tensione; la risoluzione ($\Delta x=T_w/1000$, $\Delta y=FS/2^8$) vale come incertezza assoluta con $u=\Delta/\sqrt3$ (e il manuale dà formule per l'incertezza estesa, $k=2$). Periodo su $n$ periodi: $\Delta x/n$. Grandezze derivate si trattano con la propagazione: pendenza $S=\Delta y/\Delta x$ e guadagno $G=\Delta y_2/\Delta y_1$ con le incertezze relative in quadratura; duty cycle $D=\tau/T$ (la base dei tempi non conta, restano le risoluzioni); sfasamento $\varphi=360^\circ,\tau_d/T$; durata di un fronte $T_s=t_2-t_1$ con $u(T_s)=\sqrt2,u(t)$. Per $V_0$ e $V_{DC}$ da due cursori: $\frac{V_{C1}-V_{C2}}2$ e $\frac{V_{C1}+V_{C2}}2$.Misure con i cursori e incertezza dell'oscilloscopio →).

Tempo di salita e durata dell'impulso. Si stimano il livello basso $V_L$ e il livello alto $V_H$ (da istogrammaconteggio di quanti campioni cadono in ciascun intervallo di valori: mostra i due livelli di un'onda rettangolare o da medie dei tratti piatti) e l'ampiezza $V_H-V_L$; il tempo di salita è la differenza tra gli istanti in cui il segnale attraversa $V_L+0{,}1(V_H-V_L)$ e $V_L+0{,}9(V_H-V_L)$; durata dell'impulso e periodo si misurano al livello 50% (massima pendenza); il duty cyclerapporto tra durata dell'impulso e periodo è il rapporto.

Fit di una sinusoide a tre parametri

Se la frequenza $f_0$ è nota (è quella del generatore), il segnale $x[n]=A\sin(\omega_0nT_s+\varphi)+C$ si riscrive con $B_c=A\sin\varphi$, $B_s=A\cos\varphi$: $$x[n]=B_c\cos(\omega_0nT_s)+B_s\sin(\omega_0nT_s)+C,$$ lineare nei tre coefficienti. Si risolve ai minimi quadratimetodo che sceglie i parametri minimizzando la somma dei quadrati degli scarti tra modello e dati: $\mathbf x=\mathbf M\mathbf p$ con $\mathbf M=[\cos\ \ \sin\ \ 1]$ ($N\times3$) e $\hat{\mathbf p}=(\mathbf M^T\mathbf M)^{-1}\mathbf M^T\mathbf x$. Si ricavano ampiezza $\hat A=\sqrt{B_c^2+B_s^2}$, fasespostamento temporale della sinusoide espresso come angolo $\hat\varphi=\operatorname{atan2}(B_c,B_s)$ e offset $C$. È l'idea del metodo "a quattro parametri" dello standard per i convertitori A/D (che stima anche $f_0$). Con $N$ campioni e rumore $\sigma$ l'incertezza standard dell'ampiezza è circa $\sigma\sqrt{2/N}$: dimezzarla richiede quadruplicare i campioni (come la media di tipo A).

Nel dominio della frequenza

La FFT (Spettro dei segnali periodici e FFT dell'oscilloscopioUn segnale periodico di periodo $T$ ha uno spettro a righe alle frequenze $kf_0$, $f_0=1/T$. Per un'onda rettangolare $0\to A$ con duty cycle $d$ la riga $k$ ha ampiezza $2Ad,|\mathrm{sinc}(kd)|$ (inviluppo $|\sin x/x|$, primo zero a $f=1/\tau$); se $d\ne50%$ compaiono anche le armoniche pari. La FFT dell'oscilloscopio lavora sulla finestra visualizzata $T_w$: la risoluzione in frequenza è $\Delta f=1/T_w$ e si migliora solo aumentando $T_w$ (scala dei tempi più lenta), non con più campioni. Si ottengono righe pulite quando la finestra contiene un numero intero di periodi (campionamento coerente), altrimenti c'è dispersione spettrale.Spettro dei segnali periodici e FFT dell'oscilloscopio →, Trasformata di Fourier discreta (TFD)La TFD è la serie di Fourier dei segnali discreti periodici di periodo $N$: servono solo $N$ armoniche, $X(k)=\frac1N\sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N}$ e $x(n)=\sum_{k=0}^{N-1}X(k)e^{j2\pi kn/N}$. È un prodotto matrice-vettore in $\mathbb{C}^N$ (la matrice di Fourier); la convoluzione circolare diventa prodotto, $N X(k)Y(k)$. La FFT la calcola in $O(N\log N)$.Trasformata di Fourier discreta (TFD) →) dà un'ampiezza per ogni riga distanziata di $f_s/N$. Se $f_0$ non è un multiplo esatto di $f_s/N$ (campionamento non coerentein N campioni non entra un numero intero di periodi del segnale) l'energia si distribuisce sulle righe vicine e l'ampiezza sul picco è sottostimata fino a $36\%$ (caso peggiore con finestra rettangolare, $-3{,}9$ dB). Rimedi: finestra di Hann $w[n]=\frac12\big(1-\cos\frac{2\pi n}{N}\big)$, che attenua i bordi e riduce la dispersione, con le ampiezze normalizzate dividendo per $\frac12\sum w[n]$ (così una sinusoide di picco $A$ al centro di una riga vale $A$); e interpolazione parabolicastima della posizione del picco facendo passare una parabola per i tre valori attorno al massimo del picco sui logaritmi dei tre valori attorno al massimo ($k-1,k,k+1$): $\delta=\frac12\frac{a-c}{a-2b+c}$ con $a=\ln|X_{k-1}|$, $b=\ln|X_k|$, $c=\ln|X_{k+1}|$, e $\hat f=(k+\delta)\frac{f_s}N$.

Esempio completo

Sinusoide $v=1{,}5+2\sin(2\pi\,3517\,t+0{,}7)$ V, $f_s=100$ kS/s, $N=1000$ campioni (periodi non interi: $N f_0/f_s=35{,}17$), rumore gaussianorumore con distribuzione normale e valor medio nullo $\sigma=20$ mV.

import numpy as np
rng = np.random.default_rng(3)
fs, N, f0, A, C, ph = 100e3, 1000, 3517.0, 2.0, 1.5, 0.7
t = np.arange(N) / fs
x = C + A * np.sin(2 * np.pi * f0 * t + ph) + 0.02 * rng.standard_normal(N)

# nel tempo: media, picco-picco e frequenza dagli attraversamenti del valor medio
L = x.mean()
i = np.where((x[:-1] < L) & (x[1:] >= L))[0]
tc = t[i] + (L - x[i]) / (x[i + 1] - x[i]) / fs
print(x.mean(), x.max() - x.min(), (len(tc) - 1) / (tc[-1] - tc[0]))

# FFT con finestra rettangolare e di Hann, interpolazione parabolica
f = np.fft.rfftfreq(N, 1 / fs)
X = np.fft.rfft(x) / N
w = np.hanning(N)
Xw = np.fft.rfft((x - x.mean()) * w) / (w.sum() / 2)
k = np.argmax(abs(Xw[1:])) + 1
a, b, c = np.log(abs(Xw[k - 1:k + 2]))
delta = 0.5 * (a - c) / (a - 2 * b + c)
print(2 * abs(X[np.argmax(abs(X[1:])) + 1]), abs(Xw[k]), (k + delta) * fs / N)

# fit a tre parametri con f0 nota
M = np.column_stack([np.cos(2 * np.pi * f0 * t), np.sin(2 * np.pi * f0 * t), np.ones(N)])
Bc, Bs, Cc = np.linalg.lstsq(M, x, rcond=None)[0]
print(np.hypot(Bc, Bs), Cc)

Risultati ($A=2{,}000$ V, $C=1{,}500$ V, $f_0=3517{,}0$ Hz veri):

Il fit e il conteggio degli attraversamenti sono i più accurati perché sfruttano molti campioni e un modello del segnale; la FFT rettangolare non coerente è la peggiore per l'ampiezza e per la frequenza (risoluzione $f_s/N=100$ Hz).

Tempo di salita da un gradino del primo ordine. Per $y=1-e^{-t/\tau}$ con $\tau=1\ \mu$s, campionato a $100$ MS/s con rumore, gli attraversamenti interpolati del $10\%$ e del $90\%$ danno $T_r=2{,}19\ \mu$s, in accordo con $\tau\ln9=2{,}20\ \mu$s. Duty cycle di un'onda rettangolare con $T=192\ \mu$s e $\tau=58\ \mu$s a $50$ MS/s, dagli attraversamenti al $50\%$: periodo medio $192{,}0\ \mu$s e durata $58{,}0\ \mu$s ($D=30{,}2\%$), con risoluzione molto migliore dell'intervallo di campionamento di $20$ ns grazie all'interpolazione.

Errori comuni

• Leggere l'ampiezza dal massimo grezzo della FFT senza sapere se il campionamento è coerente.
• Calcolare il picco-picco con $\max-\min$ su un segnale rumoroso e crederlo il vero picco-picco.
• Usare gli attraversamenti senza interpolare: la risoluzione resta $T_s$.
• Interpolare il picco di una FFT con la finestra rettangolare: la formula parabolica vale con finestre a lobo principale largo (Hann).