Esercizio - Conversione di frequenza da 44,1 kHz a 48 kHz

Testo (esempio della lezione 20 del corso Multimedia Signal Processing, UniPD). Un segnale audio campionato a $F_s=44{,}1$ kHz deve essere portato a $F_s'=48$ kHz.

1. Determinare i fattori $L$ e $M$ del sistema di conversione e la frequenza intermedia.
2. Disegnare lo schema a blocchi e dire quale taglio deve avere il filtro (in rad e in Hz) e quale guadagno.
3. Il filtro deve avere banda passante fino a $20$ kHz, banda oscura da $22{,}05$ kHz e attenuazione $80$ dB: stimare la lunghezza con la finestra di Kaiser e confrontare il numero di moltiplicazioni al secondo con la realizzazione ingenua e con quella polifase.
4. Cosa succederebbe nella conversione inversa ($48\to44{,}1$ kHz) senza filtro?

Teoria usata: Elaborazione multirate - decimazione e interpolazioneL'elaborazione multirate cambia la frequenza di campionamento di un segnale. Interpolazione per $L$ (da $F_s$ a $LF_s$): l'espansore inserisce $L-1$ zeri fra i campioni, $W(\hat\omega)=X(L\hat\omega)$, lo spettro si ripete $L$ volte nell'asse più ampio; un passabasso con taglio $\pi/L$ e guadagno $L$ elimina le immagini, e con il sinc ideale $h[n]=\operatorname{sinc}(n/L)$ i campioni originali sono conservati ($y[nL]=x[n]$); in pratica FIR a finestre o interpolatore lineare (triangolo). Decimazione per $M$ (da $F_s$ a $F_s/M$): $y[n]=x[nM]$, $Y(\hat\omega)=\frac1M\sum_{k=0}^{M-1}X\big(\frac{\hat\omega-2\pi k}M\big)$, repliche sovrapposte (aliasing) se la banda supera $\pi/M$: si fa precedere da un passabasso con taglio $\pi/M$. Conversione razionale $F_s'=\frac LMF_s$ ($L,M$ coprimi): espansore $L$, un solo passabasso di taglio $\min(\frac\pi L,\frac\pi M)$ e guadagno $L$, decimatore $M$; $y[n]=\sum_kx[k],g[nM-kL]$ non è una convoluzione e il sistema è periodicamente tempo-invariante. Realizzare il tutto in più stadi (fattori piccoli prima per l'interpolazione, grandi prima per la decimazione) e con strutture polifase riduce molto il costo.Elaborazione multirate - decimazione e interpolazione →, Progetto di filtri FIR con il metodo delle finestreUn filtro FIR a fase lineare di tipo I (ordine $N$ pari, $h[n]=h[N-n]$) ha risposta $H(e^{j\hat\omega})=e^{-j\hat\omega N/2}\bar H(\hat\omega)$ con ampiezza $\bar H(\hat\omega)=\sum_{n=0}^{N/2}p_n\cos(n\hat\omega)$. Per approssimare un'ampiezza desiderata $D(\hat\omega)$ (passa-basso: $1$ in banda passante, $0$ in banda oscura, con tolleranze $\delta_p,\delta_s$ e frequenze $\hat\omega_p,\hat\omega_s$) ci sono tre metodi. Finestre: si tronca la serie di Fourier di $D$, cioè $h[n]=h_d[n-N/2],w[n]$ con $h_d[n]=\frac{\hat\omega_0}{\pi}\operatorname{sinc}\frac{\hat\omega_0 n}{\pi}$; la rettangolare è ottima in errore quadratico ma dà il fenomeno di Gibbs (9%, $21$ dB), le finestre rastremate (Hann, Hamming, Blackman, Kaiser) abbassano i lobi laterali allargando la transizione ($\hat\omega_s-\hat\omega_p\simeq\alpha,2\pi/L$). Campionamento in frequenza: $h=\mathrm{IDFT}$ dei campioni di $D$, esatto solo sui campioni. Minimax (Parks-McClellan): errore pesato minimo nel caso peggiore, soluzione equiripple con almeno $r+2$ alternanze, ordine più basso a parità di specifiche.Progetto di filtri FIR con il metodo delle finestre →, 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 →.

1. Fattori di conversione

Il rapporto $\frac{F_s'}{F_s}=\frac{48000}{44100}$ si riduce ai minimi termini dividendo per il massimo comune divisore $300$: $$\frac{F_s'}{F_s}=\frac{160}{147}\ \Rightarrow\ L=160,\quad M=147,\quad\gcd(160,147)=1 \text{ (coprimi)} .$$ Si interpola per $L=160$ e si decima per $M=147$. La frequenza intermedia, la più alta del sistema, è $F_s''=L\,F_s=160\cdot44{,}1\text{ kHz}=7056\text{ kHz}$ ($=M\,F_s'=147\cdot48$ kHz, e coincide con $\operatorname{mcm}(44100,48000)=7\,056\,000$ Hz). L'ordine dei blocchi è: prima si interpola, poi si decima.

2. Schema e filtro

$$x[n]\ (44{,}1\text{ kHz})\ \to\ \boxed{\uparrow160}\ \to\ \boxed{g[m]}\ (7056\text{ kHz})\ \to\ \boxed{\downarrow147}\ \to\ y[n]\ (48\text{ kHz}).$$ Un solo filtro $g$ fa sia l'anti-immagine sia l'anti-aliasing, con taglio $\hat\omega_c=\min\big(\frac\pi L,\frac\pi M\big)=\min\big(\frac\pi{160},\frac\pi{147}\big)=\frac\pi{160}$ (rispetto a $F_s''$): in Hz $f_c=\hat\omega_c\frac{F_s''}{2\pi}=\frac{7056}{320}=22{,}05$ kHz, cioè $\frac{F_s}2$ (essendo $F_s<F_s'$ il taglio è fissato dall'ingresso: tutta la banda dell'ingresso passa, quella di uscita è più larga). Il guadagno in banda passante è $L=160$ (che compensa gli $L-1$ zeri inseriti).

3. Costo del filtro

Banda di transizione da $20$ a $22{,}05$ kHz, cioè $2{,}05$ kHz a $F_s''=7056$ kHz: $\Delta\hat\omega=2\pi\frac{2050}{7\,056\,000}=1{,}825\cdot10^{-3}$ rad. Con la formula di Kaiser: $$N\ge\frac{A_s-8}{14{,}36\,\frac{\Delta\hat\omega}{2\pi}}=\frac{72}{14{,}36\cdot2{,}905\cdot10^{-4}}=17\,258 ,$$ e scipy.signal.kaiserord dà $N=17\,275$, $\beta=7{,}857$: circa $17\,300$ coefficienti. Un filtro così lungo calcolato ingenuamente a $7056$ kHz richiede $$17\,276\ \text{moltiplicazioni}\times7\,056\,000\ \text{campioni/s}\simeq1{,}2\cdot10^{11}\ \text{moltiplicazioni/s}.$$ Impossibile in tempo reale. La struttura efficiente sfrutta due fatti: dopo l'espansore $L-1=159$ campioni su $160$ sono zero (non servono moltiplicazioni) e il decimatore butta via $146$ uscite su $147$ (non vanno calcolate). Restano solo $\frac{N+1}{L}\simeq108$ prodotti per campione d'uscita: $$108\times48\,000\simeq5{,}2\cdot10^{6}\ \text{moltiplicazioni/s},$$ cioè $1{,}2\cdot10^{11}/5{,}2\cdot10^6=23\,520=L\cdot M$ volte meno. I coefficienti usati variano periodicamente da un'uscita all'altra (sistema periodicamente tempo-invariante, con periodo $L$ uscite): è la struttura polifase.

Esempio. scipy.signal.resample_poly(x, 160, 147) fa esattamente questo, con un filtro di Kaiser ($\beta=5$) di $2\cdot10\cdot\max(L,M)+1=3201$ coefficienti per default. Convertendo $1$ s di un tono a $1$ kHz campionato a $44{,}1$ kHz si ottengono $48\,000$ campioni e, a parte i bordi, la differenza dal tono ideale a $48$ kHz è $6{,}7\cdot10^{-4}$ (il filtro scelto è più corto di quello da $80$ dB, quindi ha una transizione più larga).

4. Il caso inverso

Da $48$ kHz a $44{,}1$ kHz si ha $L=147$ e $M=160$, e il taglio è $\min\big(\frac\pi{147},\frac\pi{160}\big)=\frac\pi{160}$ rispetto a $7056$ kHz: $22{,}05$ kHz, fissato dall'uscita. Il filtro è necessario, perché le frequenze fra $22{,}05$ e $24$ kHz dell'ingresso non possono essere rappresentate a $44{,}1$ kHz. Senza filtro (per esempio prendendo semplicemente il campione più vicino all'istante $m/44100$) un tono a $23$ kHz verrebbe ripiegato (aliasing) a $44{,}1-23=21{,}1$ kHz con ampiezza piena: il picco dello spettro del segnale così ottenuto cade a $21\,100$ Hz (verificato numericamente).