Codifica numerica del segnale audio/Appendici

by Daniele Sereno/Pietro Valocchi

Appendici

../Capitolo 7 ../Bibliografia IncludiIntestazione 17 marzo 2021 75% Da definire

Questo testo è completo.

Daniele Sereno, Pietro Valocchi - Codifica numerica del segnale audio (1996)

Appendici

◄

Capitolo 7

Bibliografia

►

Appendice A

PROCESSI AUTOREGRESSIVI

---

A.1 SISTEMI LINEARI TEMPO INVARIANTI A TEMPO DISCRETO

Un sistema discreto è individuato dalla trasformazione T che lega l’ingresso x con l’uscita y

$y(n)=T\left\{x(n)\right\}$

(A.1)

Affinché un sistema discreto sia lineare deve risultare che la sua uscita ad una combinazione lineare di ingressi sia la stessa combinazione lineare alle singole uscite, cioè

$T\left\{a\ x_{1}(n)+b\ x_{2}(n)\right\}=T\left\{a\ x_{1}(n)\right\}+T\left\{b\ x_{2}(n)\right\}$

(A.2)

dove a e b sono delle costanti. Invece, un sistema è detto invariante alle traslazioni se, nota l’uscita y ad un certo ingresso x, la sua uscita ad un ingresso costituito dal segnale x traslato nel tempo coincida con la versione traslata della y

$y(n)=T\left\{x(n)\right\}$

$x(n-n_{0})=>y(n-n_{0})$

(A.3)

Rappresentando un segnale discreto x(n) come una serie di impulsi δ(n)

$x(n)=\sum _{k=-\infty }^{\infty }x(k)\delta (n-k)$

(A.3)

dove

$\delta (n)=\left\{{\begin{matrix}1&n=0\\0&n\neq 0\end{matrix}}\right.$

(A.5)

per sistemi lineari tempo-invarianti, la relazione

$y(n)=T\left\{x(n)\right\}=T\left\{\sum _{k=-\infty }^{\infty }x(k)\delta (n-k)\right\}$

(A.6)

può essere espressa, per la linearità, come

$y(n)=\sum _{k=-\infty }^{\infty }x(n)T\left\{x(k)\delta (n-k)\right\}$

(A.7)

Indicando con

$h(n)=T\left\{\delta (n)\right\}$

(A.8)

la risposta impulsiva del sistema, per l’invarianza alla traslazione si ottiene

$h(n-k)=T\left\{\delta (n-k)\right\}$

(A.9)

da cui

$y(n)=\sum _{k=-\infty }^{\infty }x(k)h(n-k)=x(n)\otimes h(n)$

(A.10)

Si può quindi concludere che un sistema discreto lineare tempo invariante è completamente caratterizzato dalla sua risposta impulsiva o, equivalentemente, dalla funzione di trasferimento ottenuta come trasformata z della stessa

$H(z)={\frac {Y(z)}{X(z)}}=\sum _{n=-\infty }^{\infty }h(n)z^{-n}$

(A.11)

Lo svantaggio fondamentale della sommatoria di convoluzione introdotta per l’analisi dei sistemi nel dominio del tempo risiede nel fatto che gli estremi della sommatoria risultano pari a ±∞. Per evitare tale inconveniente, il legame ingresso-uscita può essere approssimato tramite equazioni alle differenze. Una equazione lineare alle differenze a coefficienti costanti di ordine n è un legame tra tra gli ultimi M campioni dell’ingresso e gli ultimi N dell’uscita del tipo

$\sum _{k=0}^{N}a_{k}y(n-k)=\sum _{k=0}^{M}b_{k}x(n-k)$

(A.12)

A seconda della struttura dell’equazione alle differenze, si possono distinguere tre differenti tipi di sistemi. Se N = 0 (sistema Moving Average: MA), l’equazione alle differenze si riduce a

$y(n)=\sum _{k=0}^{M}{\frac {b_{k}}{a_{0}}}x(n-k)$

(A.13)

Tale equazione non è più ricorsiva, in quanto dipende esclusivamente dall'ingresso. Confrontando tale equazione con quella di definizione della risposta impulsiva di un sistema, si nota come i coefficienti b_k/a₀ coincidano con i campioni della risposta impulsiva h(k). Posto x(n) = δ(n) si ha che la risposta impulsiva è data da

$h(n)=\sum _{k=0}^{N}{\frac {b_{k}}{a_{0}}}\delta (n-k)$

(A.14)

ed è di durata finita: il sistema è definito FTR (Finite Impulse Response). Se M = 0 (sistema Auto Regressive: AR), l’equazione alle differenze si semplifica come

$y(n)=\sum _{k=1}^{N}{\frac {a_{k}}{a_{0}}}y(n-k)+{\frac {b_{0}}{a_{0}}}x(n)$

(A.15)

A causa della retrazione dell’uscita, la risposta impulsiva del filtro è di durata finita (sistema IIR: Infinite Impulse Response). Il caso più generale è quello in cui sia M che N non sono nulli. Il sistema è quindi governato dall'equazione alle differenze

$y(n)=\sum _{k=0}^{M}{\frac {b_{k}}{a_{0}}}x(n-k)-\sum _{k=1}^{N}{\frac {a_{k}}{a_{0}}}y(n-k)$

(A.16)

ed il sistema (di tipo IIR) è detto ARMA. Esprimendo l’equazione alle differenze

$\sum _{k=0}^{N}a_{k}y(n-k)=\sum _{k=0}^{M}b_{k}y(n-k)$

(A.17)

tramite la sua trasformata Z, si ottiene

$Y(z)\sum _{k=0}^{N}a_{k}z^{k}=X(z)\sum _{k=0}^{M}b_{k}z^{-k}$

(A.18)

con funzione di trasferimento

$H(z)={\frac {Y(z)}{X(z)}}={\frac {\sum _{k=0}^{M}b_{k}z^{-k}}{\sum _{k=0}^{N}a_{k}z^{-k}}}$

(A.19)

data da un rapporto di polinomi in z con a₀ = 1. Se b₀ ≠ 0, è possibile eliminare le potenze negative di z fattorizzando i termini b₀z^-M e a₀z^-N ed ottenendo

$H(z)=G\ z^{-M+N}{\frac {\prod _{k=1}^{M}(z-z_{k})}{\prod _{k=1}^{N}(z-p_{k})}}$

(A.20)

dove G = b₀/a₀ = b₀. Per sistemi FIR il denominatore non è presente e la funzione di trasferimento, è composta da soli zeri. Per sistemi AR la funzione di trasferimento è composta da soli poli.

Essendo la funzione di trasferimento dipendente dal contributo dei singoli poli e zeri, è possibile prevedere le caratteristiche della risposta impulsiva di un sistema a tempo discreto (e, quindi, il suo comportamento nel dominio del tempo) in funzione della loro localizzazione. Mentre il contributo del numeratore è poco intuitivo (ma sarà analizzato nel seguito), l’influenza della localizzazione dei poli è, invece, molto evidente. Imponendo il vincolo che l'uscita del sistema discreto sia una serie di valori reali, anche l’equazione alle differenze deve risultare a coefficienti reali. Nel caso di funzione di trasferimento con un solo polo, ciò comporta che il polo stesso si trovi sull'asse reale. Per quanto riguarda la risposta impulsiva si ha

$H(z)={\frac {1}{1-rz^{-1}}}\leftrightarrow ^{z}h(n)=r^{n}\ u(n)$

(A.21)

Nel caso in cui il polo si trovi sul semiasse positivo (r > 0), la risposta impulsiva risulta decadere esponenzialmente, essere costante o crescere esponenzialmente a secondo che il suo modulo risulti minore, uguale o maggiore di 1. Il valore di r determina la rapidità del decadimento. Nel caso in cui il polo si trovi sul semiasse negativo (r < 0), l’inviluppo della risposta impulsiva ha Io stesso andamento precedentemente descritto in funzione del modulo, solo che i campioni risultano a segni alternati. Il comportamento non cambia sostanzialmente nel caso di poli reali coincidenti

$H(z)={\frac {1}{\left(1-rz^{-1}\right)^{2}}}\leftrightarrow ^{z}h(n)=n\ r^{n}\ u(n)$

(A.22)

Per quanto riguarda poli complessi coniugati

$H(z)={\frac {r\ sin\omega _{0}z^{-1}}{\left(1-re^{j\omega _{0}}z^{-1}\right)\left(1-re^{-j\omega _{0}}z^{-1}\right)}}\leftrightarrow ^{z}h(n)=r^{n}\ sin(\omega _{0}n)\ u(n)$

(A.23)

la risposta impulsiva è oscillante con inviluppo decrescente, costante o crescente, secondo che i poli si trovino all'interno, sulla o all'esterno della circonferenza di raggio unitario.

Se un sistema è stabile, l’uscita prodotta a seguito di un ingresso limitato deve risultare limitata. Essendo l’uscita funzione della risposta impulsiva, dall'analisi precedente segue che un sistema è stabile se ha poli esclusivamente all'interno della circonferenza unitaria.

Per quanto riguarda la risposta in frequenza del sistema discreto LTI, definita come trasformata della sua risposta impulsiva, si nota come essa coincida con la H(z) calcolata sulla circonferenza di raggio unitario

$H(\omega )=\sum _{n=-\infty }{\infty }h(n)\ e^{-j\omega n}=H(z)|_{z=e^{j\omega }}$

$H(\omega )=b_{0}e^{j\omega (N-M)}{\frac {\prod _{k=1}^{M}\left(e^{j\omega }-z_{k}\right)}{\prod _{k=1}^{N}\left(e^{j\omega }-p_{k}\right)}}$

(A.24)

Interpretando i termini nelle produttorie, e^jω rappresenta il vettore uscente dall'origine con vertice sulla circonferenza unitaria nel punto di pulsazione CO, mentre z_k e p_k rappresentano i vettori che congiungono l’origine con gli zeri ed i poli (fig. A.l). I termini tra parentesi tonde, quindi, sono i vettori congiungenti i poli e gli zeri con la circonferenza, che in forma polare sono esprimibili come

${\begin{cases}e^{j\omega }-z_{k}=V_{k}(\omega )e^{j\Theta _{k}(\omega )}\\e^{j\omega }-p_{k}=U_{k}(\omega )e^{j\Phi _{k}(\omega )}\end{cases}}$

(A.25)

con

${\begin{cases}V_{k}(\omega )=|e^{j\omega }-z_{k}|;\Theta _{k}(\omega )=\angle \left(e^{j\omega }-z_{k}\right)\\U_{k}(\omega )=|e^{j\omega }-p_{k}|;\Phi _{k}(\omega )=\angle \left(e^{j\omega }-p_{k}\right)\end{cases}}$

(A.26)

Noti i moduli V_k e U_k dei vettori congiungenti i poli e gli zeri con la circonferenza unitaria e le loro rotazioni Θ_k e Φ_k, il modulo della funzione di trasferimento si ottiene come prodotto dei moduli

$|H(\omega )|=|b_{0}|{\frac {\prod _{k=1}^{M}V_{k}(\omega )}{\prod _{k=1}^{N}U_{k}(\omega )}}$

(A.27)

Analizzando il contributo dei singoli poli e zeri, si vede che il modulo della funzione di trasferimento si riduce per valori di © corrispondenti a punti della circonferenza unitaria prossimi ad uno zero, in quanto si riduce il relativo termine V_k. Per quanto riguarda i poli, invece, la riduzione del temine U_k si Fig. A.1 - Interpretazione geometrica del contributo di un polo al modulo ed alla fase della funzione di trasferimento.

traduce in un’esaltazione della risposta in frequenza, dato che tale termine si trova al denominatore della H(ω). È importante sottolineare come, a parità di modulo di due poli complessi coniugati, l’ampiezza del picco di risonanza dipenda dalla loro frequenza (fig. A.2). Viceversa, a parità di frequenza e interpretando il contributo di un polo come un filtraggio passa banda nell'intorno della sua frequenza centrale, la banda del filtro si riduca all'aumentare del modulo del polo. Tale miglioramento dal punto di vista della selettività in frequenza si paga, a causa della riduzione dello smorzamento, dal punto di vista della selettività temporale.

Per quanto riguarda la fase, essa dipende dalla differenza tra la somma delle fasi di ciascuno zero e quella dei poli

$\angle H(\omega )=\angle b_{0}+\omega \cdot (M-N)+\sum _{k=1}^{M}\Theta _{k}(\omega )-\sum _{k=1}^{N}\Phi _{k}(\omega )$

(A.28)

Nel caso di sistemi discreti LTI posti tra loro in cascata, si può verificare che la funzione di trasferimento complessiva è data dal prodotto delle funzioni di trasferimento dei singoli sistemi. Infatti

$Y_{2}(z)=X_{2}(z)\cdot H_{2}(z)=Y_{1}(z)\cdot H_{2}(z)=X_{1}(z)\cdot H_{1}(z)\cdot H_{2}(z)$

(A.29)

Ipotizzando di voler realizzare il sistema inverso ad uno dato. Ciò vuol dire che si vuole realizzare un sistema che messo in cascata ad un secondo produce una funzione di trasferimento complessiva unitaria. A tal fine è immediato verificare che è sufficiente realizzare un sistema che abbia come poli gli zeri del sistema da invertire e viceversa. Affinché anche il sistema inverso sia stabile, è necessario che i suoi poli siano all'interno della circonferenza unitaria, il che si traduce con il vincolo che il sistema originario, oltre ad avere i propri poli all'interno della circonferenza unitaria per la sua stabilità, vi abbia anche gli zeri. Sistemi con sia poli che zeri all'interno della circonferenza unitaria sono detti a fase minima.

A.2 FUNZIONE DI AUTOCORRELAZIONE DI UN PROCESSO AUTOREGRESSIVO

L’equazione caratteristica di un processo AR, soddisfatta sia dal processo u(n)

$u(n)=v(n)-\sum _{k=1}^{q}w_{k}\ u(n-k)$

(A.30)

che dall’autocorrelazione R(n)

$\sum _{k=0}^{q}w_{k}R_{uu}(n-k)=0;{\begin{cases}n=1,2,...,q\\w_{0}=1\end{cases}}$

(A.31)

è pari a

$1+\sum _{k=1}^{q}w_{k}z^{-k}=0;$

(A.32)

Fig. A.2a - Legame tra posizione dei poli, risposta impulsiva e risposta in frequenza. Fig. A.2b - Legame tra posizione dei poli, risposta impulsiva e risposta in frequenza. Per la stabilità, le radici p dell’equazione caratteristica debbo risiedere all’interno del circolo di raggio unitario. Si vuole mostrare come, variando i coefficienti del sistema AR, si riesca a variare le caratteristiche della funzione di autocorrelazione del processo generato.

Per un processo del secondo ordine, la funzione caratteristica è pari a

$1+w_{1}z^{-1}+w_{2}z^{-2}=0$

(A.33)

con poli

$p_{1}\cdot p_{2}={\frac {-w_{1}\pm {\sqrt {w_{1}^{2}-4w_{2}}}}{2}}{\begin{cases}|p_{1}|<1\\|p_{2}|<1\end{cases}}\to {\begin{cases}w_{2}+w_{1}\geq -1\\w_{2}-w_{1}\geq 1\\-1\leq w_{2}\leq 1\end{cases}}$

(A.34)

Imponendone la stabilità, si ricava la zona di esistenza dei coefficienti (fig. A.3)

${\begin{cases}|p_{1}|<1\\|p_{2}|<1\end{cases}}\to {\begin{cases}w_{2}+w_{1}\geq -1\\w_{2}-w_{1}\geq 1\\-1\leq w_{2}\leq 1\end{cases}}$

(A.35)

Fig. A.3 - Regione di esistenza per i coefficienti del processo AR. All'interno di tale zona, cambiando i coefficienti del processo, la funzione di autocorrelazione assume le seguenti caratteristiche

${\begin{cases}w_{1}^{2}--4w_{2}>0\to radici\ reali{\begin{cases}w_{1}<0\to {\begin{matrix}radice\ dominante\ positiva;\\decadimento\ esponenziale\end{matrix}}\\w_{1}>0\to {\begin{matrix}radice\ dominante\ negativa;\\decadimento\ esponenziale\ a\ segni\ alterni\end{matrix}}\end{cases}}\\w_{1}^{2}--4w_{2}<0\to radici\ complesse\ coniugate\to oscillazioni\end{cases}}$

(A.36)

L’andamento della funzione di autocorrelazione è riportata in figura A.4 per i tre casi. Fig. A.4a - Funzioni di autocorrelazione di processi autoregressivi. Fig. A.4b - Funzioni di autocorrelazione di processi autoregressivi.

Appendice B

METODO DEI MÌNIMI QUADRATI RICORSIVO

---

Il vantaggio principale dell’LMS è la sua semplicità. Per contro, il mantenere fisso il passo di aggiornamento è svantaggioso in quanto esso risulta essere solitamente troppo piccolo nelle fasi iniziali, ritardando la convergenza, ed eccessivo nelle fasi finali, deteriorando le prestazioni. Inoltre, dipendendo l’LMS solamente dall'errore istantaneo, non ha memoria dei precedenti passi di elaborazione. È possibile introdurre un differente algoritmo, derivato dall’LS, in cui l'aggiornamento dei coefficienti del predittore avviene indipendentemente e ricorsivamente per ciascuno di essi in funzione dell’andamento della stima (Recursive Least-Square: RLS). Ciò introduce memoria nell'algoritmo, come mostrato nel seguito, ma questo non rappresenta uno svantaggio nel caso di fenomeni stazionari. Nel caso di processi variabili, invece, è necessario attenuare la memoria sui campioni più remoti man mano che si procede nella stima. Ciò si ottiene con l’adozione di un nuovo criterio di ottimizzazione che introduce un fattore di pesatura w < 1 nella definizione dell’errore quadratico

${\hat {\epsilon }}(n)=\sum _{n=-0}^{t}w^{t-n}e^{(n)^{2}}$

(B.1)

È poi necessario introdurre la funzione di errore a priori, definita come

$\eta (n)=x(n)-\alpha ^{{n-1}^{T}}x(n-1)$

(B.2)

nella quale si stima il campione corrente in funzione dei pesi calcolati nel passo precedente, che risulta solitamente differente dalla funzione di errore a posteriori

$e^{(n)}=x(n)-\alpha ^{(n)^{T}}x(n-1)$

(B.3)

nella quale si utilizza il vettore dei coefficienti aggiornato. Ripetendo i passi già presentati nel capitolo 5, la m inim izzazione dell’errore in funzione dei coefficienti del predittore porta al passo n-esimo all’equazione normale

$\Gamma ^{(n)}=\alpha ^{(n)}\Theta ^{(n)}$

(B.4)

che risulta formalmente identica a quella utilizzata nella derivazione delle equazioni di Wiener-Hopf, ma nelle quale è implicitamente presente il fattore di pesatura w. La matrice di autocorrelazione dell’ingresso ed il vettore di cross-correlazione tra l’ingresso e l’uscita desiderata risultano, infatti, pari a

$\Gamma ^{(n)}=\sum _{i=1}^{n}w^{n-i}x(i-1)x^{T}(i-1)$

$\Theta ^{(n)}=\sum _{i=1}^{n}w^{n-i}x(i)x^{T}(i-1)$

(B.5)

La soluzione di tale sistema è ancora dato da

$\alpha ^{(n)}=\Gamma ^{(n)^{-1}}\Theta ^{(n)}$

(B.6)

anche se, in questo caso, ne ricerca una sua soluzione ricorsiva. Per quanto riguarda la matrice di autocorrelazione, è possibile dame una definizione ricorsiva isolandone dalla definizione il termine per i = n

$\Gamma ^{(n)}=w\left[\sum _{i=0}^{(n-1)}w^{n-1-i}x(i-1)x^{T}(i-1)\right]+x(n-1)x^{T}(n-1)$

$\Gamma ^{(n)}=w\Gamma ^{(n-1)}+x(n-1)x^{T}(n-1)$

(B.7)

Analogamente si ottiene

$\Theta (n)=w\Theta ^{(n-1)}+x(n)X^{T}(n-1)$

(B.8)

L'inversione della matrice di autocorrelazione può essere evitata, in quanto risulta [Hay86]

$\Gamma ^{(n)^{-1}}={\frac {\Gamma ^{(n-1)^{-1}}}{w}}-{\frac {\frac {\Gamma ^{(n-1)^{1}}x(n-1)x^{T}(n-1)\Gamma ^{(n-1)^{-1}}}{w^{2}}}{1+{\frac {x^{T}(n-1)\Gamma ^{(n-1)^{-1}}x(n-1)}{w}}}}$

(B.9)

Per semplificare la notazione si introduce la matrice di autocorrelazione inversa

$P(n)=\Gamma ^{(n)^{1}}$

(B.10)

ed il vettore dei guadagni

$k^{(n)}={\frac {\frac {P^{(n-1)}x(n-1)}{w}}{1+{\frac {x^{T}(n-1)P^{(n-1)}x(n-1)}{w}}}}$

(B.11)

che ha le stesse dimensioni dell'ordine del predittore. In tal modo, per l'equazione dell'inversa della matrice di autorrelazione si ottiene la relazione ricorsiva

$P^{(n)}={\frac {P^{(n-1)}-k^{(n)}x^{T}(n-1)P^{(n-1)}}{w}}$

(B.12)

Riorganizzando la definizione del vettore dei guadagni si ottiene, invece

$k^{(n)}={\frac {\left[P^{(n-1)}-k^{(n)}x^{T}(n-1)P^{(n-1)}\right]x(n-1)}{w}}=P^{(n)}x(n-1)$

(B.13)

con la quale si ottengono le seguenti relazioni ricorsive per i coefficienti del predittore

$\alpha ^{(n)}=P^{(n)}\Theta (n)=wP^{(n)}\Theta ^{(n-1)}+x(n)P^{(n)}x(n-1)$

$=P^{(n-1)}\Theta (n-1)=k^{(n)}x^{T}(n)P^{(n-1)}\Theta ^{(n-1)}+x(n)P^{(n)}x(n-1)$

$=\alpha ^{(n-1)}+x(n)P^{(n)}x(n-1)-k^{(n)}x^{T}(n-1)\alpha ^{(n-1)}$

$=\alpha ^{(n-1)}+k^{(n)}\left[x(n)-x^{T}(n-1)\alpha ^{(n-1)}\right]$

(B.14)

Sostituendo l’espressione dell’errore

$\eta ^{(n)}=x(n)-x^{T}(n-1)\alpha ^{(n-1)}$

(B.15)

si ottiene, infine

$\alpha ^{(n)}=\alpha ^{(n-1)}+k^{(n)}\eta ^{(n)}$

(B.16)

Da tale equazione si nota come il vettore dei guadagni k(n) ha le stesse funzioni del parametro p nell’LMS. A differenza di quest’ultimo, però, esso risulta variabile e, dipendendo dalla matrice di correlazione del segnale, ha memoria del processo di stima.

In tal modo sono state ottenute equazioni ricorsive per il calcolo di k, p, a e P. Al fine di ridurre la complessità computazionale, però, è opportuno definire la grandezza

$y^{(n)}={\frac {P^{(n-1)}x(n-1)}{w}}$

(B.17)

con la quale i passi per la predizione sono ridotti ai seguenti: - si calcola la grandezza

$y^{(n)}={\frac {P^{(n-1)}x(n-1)}{w}}$

(B.18)

- da questa si calcola il vettore dei guadagni

$k^{(n)}={\frac {y^{(n)}}{1+x^{T}(n-1)y^{(n)}}}$

(B.19)

errore quadratico medio Fig. B.1 - Curve di apprendimento per LMS e RLS.

• si calcola l’errore

$\eta (n)=x(n)=z^{T}(n-1)\alpha ^{(n-1)}$

(B.20)

• si aggiornano i parametri del predatore

$\alpha ^{(n)}=\alpha ^{(n-1)}+\eta ^{(n)}\ k^{(n)}$

(B.21)

• si calcola, infine, il valore aggiornato della matrice di autocorrelazione

$P^{(n)}={\frac {P^{(n-1)}}{w}}+k^{(n)}\ y^{(n)}$

(B.22)

Per quanto riguarda le condizioni iniziali, è necessario scegliere la matrice P(0) in modo tale che risulti non singolare. La scelta più immediata è quella di scegliere una costante δ positiva “piccola” (es.: 0.1) e fissare

$P^{(0)}={\frac {I}{\delta }}$

(B.23)

Per i coefficienti del predatore il loro valore iniziale è tipicamente nullo, mentre il coefficiente di smorzamento w deve risultare eventualmente minore, ma comunque molto prossimo ad 1 (es.: 1-10^-3).

Dal punto di vista delle prestazioni, la variabilità del vettore dei guadagni fa si che sia i tempi di convergenza che le fluttuazioni nell'intorno del vettore ottimo siano ridotti rispetto all’LMS (fig. B.l). In compenso, però, la complessità computazionale è aumentata, anche considerando versioni più sofisticate dell’algoritmo [Hay86].

Appendice C

ALGORITMI A BLOCCHI

PER IL FILTRAGGIO ADATTATIVO

---

Gli algoritmi per la soluzione dell’equazione normale della predizione a blocchi possono essere ricondotti, fondamentalmente, alle seguenti due classi (fig.C.1)

• metodo della covarianza;
• metodo dell’autocorrelazione;

La loro differenza risiede nel differente modo con il quale vengono stimati gli elementi della matrice di autocorrelazione Φ. La stima della funzione di auto correlazione viene calcolata come sommatoria dei prodotti ottenuti utilizzando due finestre di campioni di pari lunghezza progressivamente traslate. Nel secondo caso si utilizza un’unica finestra di campioni. Le traslazioni richieste nel calcolo dei prodotti avvengono estendendo tale finestra base con dei campioni nulli. Nel primo caso, invece, la finestra base di campioni viene estesa utilizzando effettivi campioni del segnale. Tale differenza porta poi ad una differente struttura della matrice Φ e quindi a differenti algoritmi per la soluzione del sistema lineare di equazioni associato.

Il metodo dell autocorrelazione è quello attualmente più utilizzato sia per l’esistenza di algoritmi più robusti nei riguardi di errori di quantizzazione (adottando strutture a traliccio), sia per resistenza di più efficienti algoritmi di calcolo (ricorsione di Schür). Fig. C.1 -Differenti finestre utilizzate nel metodo della covarianza e dell’autocorrelazione.

C.l METODO DELLA COVARIANZA

Nel metodo della covarianza gli elementi della matrice Φ vengono stimati utilizzando sequenze di “n” campioni del segnale, progressivamente sfasate. Ciascun elemento della matrice viene cioè calcolato come

$\Phi (i,j)={\frac {1}{n}}\sum _{k=0}^{n-1}x(k-i)x(k-j);{\begin{matrix}1\leq i\leq p\\1\leq j\leq p\end{matrix}}$

$={\frac {1}{n}}\sum _{k=-i}^{n-i-1}x(k)x(k+i-j)$

(C.1)

Considerare blocchi di “n” campioni del segnale vuol dire moltiplicare il segnale stesso per una finestra di pari lunghezza. La finestra utilizzata nel seguito è quella rettangolare, tralasciando considerazioni sull'opportunità di utilizzare finestre differenti.

Per il calcolo degli elementi Φ(i, k) da utilizzare nella predizione, è necessario considerare il prodotto di sequenze traslate nel tempo al massimo di “p” campioni. Se la finestra di interesse è quella relativa alla sequenza x(n), è necessario considerare fino a “p” campioni esterni a tale intervallo. Con tale approssimazione, la matrice Φ dell’equazione normale deterministica viene indicata come matrice di covarianza, da cui il nome dell’algoritmo. Si mette in evidenza come tale nome non indica che i termini della matrice rappresentino la covarianza del segnale, ma sono sempre termini di autocorrelazione. La soluzione dell’equazione si ottiene poi come

$\alpha =\Phi ^{-1}\Psi$

(C.2)

Es.: si consideri il vettore

$x(n)=[2\ -1\ 0\ -1\ 2]$

(C.3)

ottenuto considerando la funzione f = cos(t) + cos(2t) e prelevando 5 campioni nel suo periodo. Considerando un predittore del secondo ordine, l’equazione normale è

${\begin{bmatrix}\alpha _{1}\\\alpha _{2}\end{bmatrix}}={\begin{bmatrix}\Psi (1,1)&\Psi (1,2)\\\Psi (2,1)&\Psi (2,2)\end{bmatrix}}{\begin{bmatrix}\Psi (1,0)\\\Psi (2,0)\end{bmatrix}}$

(C.4)

Per il calcolo degli elementi della matrice di covarianza è necessario estendere l’intervallo di osservazione, includendo anche due campioni del precedente periodo. Nel nostro caso traslazioni verso destra o verso sinistra si equivalgono.

$x(n-1)=[-1\ 0\ -1\ 2\ -1]$

$x(n-2)=[0\ -1\ 2\ -1\ 0]$

(C.5)

È possibile, quindi, calcolare gli elementi della matrice di covarianza

$\Phi (1,0)={\frac {1}{5}}([-1\ 0\ -1\ 2\ -1][2\ -1\ 0\ -1\ 2]^{T})=-{\frac {6}{5}}$

$\Phi (2,0)={\frac {1}{5}}([0\ -1\ 2\ -1\ 0][2\ -1\ 0\ 1\ 2]^{T})={\frac {2}{5}}$

$\Phi (1,1)={\frac {1}{5}}([-1\ 0\ -1\ 2\ -1][-1\ 0\ -1\ 2\ -1]^{T})={\frac {7}{5}}$

$\Phi (2,2)={\frac {1}{5}}([0\ -1\ 2\ -1\ 0][0\ -1\ 2\ -1\ 0]^{T})={\frac {6}{5}}$

$\Phi (2,1)={\frac {1}{5}}([0\ -1\ 2\ -1\ 0][-1\ 0\ -1\ 2\ 1]^{T})=-{\frac {4}{5}}$

(C.6)

e calcolare il vettore dei coefficienti ottimi come

${\begin{bmatrix}\alpha _{1}\\\alpha _{2}\end{bmatrix}}={\begin{bmatrix}{\frac {7}{5}}&-{\frac {4}{5}}\\-{\frac {4}{5}}&{\frac {6}{5}}\end{bmatrix}}^{-1}{\begin{bmatrix}-{\frac {6}{5}}\\{\frac {2}{5}}\end{bmatrix}}={\begin{bmatrix}{\frac {15}{13}}&{\frac {10}{13}}\\{\frac {10}{13}}&{\frac {35}{26}}\end{bmatrix}}{\begin{bmatrix}-{\frac {6}{5}}\\{\frac {2}{5}}\end{bmatrix}}={\begin{bmatrix}-{\frac {14}{13}}\\-{\frac {5}{13}}\end{bmatrix}}$

(C.7)

Il calcolo della Φ può avvenire anche evitando di utilizzare campioni esterni al blocco corrente. A tal fine è necessario considerare sequenze di n-p campioni opportunamente sfasate ritagliate all'interno del blocco. Il calcolo della Φ può poi avvenire costruendo la matrice delle osservazioni

$H^{T}={\begin{bmatrix}x(p)&x(p+l)&\ldots &x(n)\\x(p-1)&x(p)&\ldots &x(n-1)\\\ldots &\ldots &\ldots &\ldots \\x(1)&x(2)&\ldots &x(n-p+1)\end{bmatrix}}$

(C.8)

dalla quale sì ricavano la matrice di auto-correlazione ed il vettore di cross-correlazione

$\Phi ={\frac {H^{T}H}{n-p}}$

$\Psi ={\frac {H^{T}x(n+1)}{n-p}}$

(C.9)

La matrice di covarianza risulta essere definita positiva e, per costruzione, simmetrica. La soluzione dell’equazione normale può essere quindi trovata evitando l’inversione della stessa, ma scomponendola tramite l’algoritmo di Cholesky nel prodotto

$\Phi =V\ D\ V^{T}$

(C.10)

dove V è una matrice triangolare inferiore con la diagonale principale unitaria, mentre D è una matrice diagonale

$V={\begin{bmatrix}1&0&0&\ldots &0\\V_{21}&1&0&\dots &0\\V_{31}&V_{32}&1&\ldots &0\\\dots &\dots &\dots &\dots &\dots \\V_{p1}&V_{p2}&V_{p3}&\ldots &1\end{bmatrix}};\ D={\begin{bmatrix}d_{1}&0&0&\ldots &0\\0&d_{2}&0&\dots &0\\0&0&d_{3}&\ldots &0\\\dots &\dots &\dots &\dots &\dots \\0&0&0&\ldots &d_{p}\end{bmatrix}}$

(C.11)

Per eseguire la fattorizzazione della matrice di covarianza si procede come segue: dato

$\Phi (i,j)=\sum _{k=1}^{j}V_{ik}\ d_{k}\ V_{jk};\ 1\leq j\leq i-1$

(C.12)

per gli elementi sulla diagonale segue

$d_{i}=\Phi (i,i)-\sum _{k=1}^{i-1}V_{ik}^{2}\ d_{k};\ i\geq 2$

(C.13)

a partire dalla condizione iniziale

$d_{1}=\varphi (1,1)$

(C.14)

mentre gli elementi della matrice V si ricavano riscrivendo l’espressione del generico elemento della matrice di covarianza isolando dalla sommatori il termine per k=j

$V_{ij}d_{j}=\Phi (i,j)\sum _{k=1}^{j-1}V_{ik}d_{k}V_{jk};\ 1\leq j\leq i-1$

$V_{ij}={\frac {\Phi (i,j)\sum _{k=1}^{j-1}V_{ik}d_{k}V_{jk}}{d_{j}}}$

(C.15)

In tal modo è possibile calcolare ricorsivamente gli elementi V_ij e d_j. Una volta fattorizzata la matrice di covarianza, anche la soluzione dell’equazione normale può essere ottenuta ricorsivamente. Infatti, definendo il vettore y come

$y=D\ V^{T}\alpha$

(C.16)

risulta

$\Phi \alpha =V\ D\ V^{T}\alpha =V\ y$

$V\ y=\Psi$

(C.17)

Essendo V una matrice triangolare, gli elementi della y possono essere calcolati ricorsivamente come

$y_{1}=\Psi _{1}$

$y_{i}=\Psi _{i}-\sum _{j=1}{i-1}V_{ij}y_{j};\ p\geq i\geq 2$

(C.18)

Una volta calcolati i componenti del vettore y è possibile ricavare i coefficienti ottimi del predittore in quanto, risultando

$V^{T}\alpha =D^{-1}y$

(C.19)

Si ottengono le relazioni ricorsive

$\alpha _{i}={\frac {y_{i}}{d_{i}}}-\sum _{j=i+1}^{p}V_{ji}\alpha _{j};\ 1\leq i\leq p-1$

(C.20)

con condizioni iniziali

$\alpha _{p}={\frac {y_{p}}{d_{p}}}$

(C.21)

Es. si consideri lo stesso vettore di campioni precedente. Applicando la

fattorizzazione di Cholesky risulta

$d_{1}=\Phi (1,1)={\frac {7}{5}}$

$V(2,1)={\frac {\Phi (2,1)}{d_{1}}}=-{\frac {4}{7}}$

$d_{2}=\Phi (2,2)-V(2,1)^{2}d_{1}={\frac {26}{35}}$

(C.22)

Fattorizzata la matrice di covarianza come

$\Phi =V\ D\ V^{T}$

$\Phi ={\begin{bmatrix}1&0\\{\frac {4}{7}}&1\end{bmatrix}}{\begin{bmatrix}{\frac {7}{5}}&0\\0&{\frac {26}{35}}\end{bmatrix}}{\begin{bmatrix}1&-{\frac {4}{7}}\\0&1\end{bmatrix}}$

(C.23)

si calcolano dapprima le variabili ausiliarie y

$V\ y=\Psi$

${\begin{bmatrix}1&0\\-{\frac {4}{7}}&1\end{bmatrix}}{\begin{bmatrix}y_{1}\\y_{2}\end{bmatrix}}={\begin{bmatrix}-{\frac {6}{5}}\\{\frac {2}{5}}\end{bmatrix}}$

(C.24)

tramite la relazione

$y_{i}=\Psi _{i}-\sum _{j=1}{i-1}V_{ij}y_{j};\ p\geq i\geq 2$

$y_{1}=\Phi (1,0)=-{\frac {6}{5}}$

$y_{2}=\Phi (2,0)-V(2,1)y_{1}=-{\frac {2}{7}}$

(C.25)

ed infine i coefficienti del predittore dalla relazione

$\alpha _{i}={\frac {y_{i}}{d_{i}}}-\sum _{j=i+1}^{p}V_{ji}\alpha _{j};\ 1\leq i\leq p-1$

$\alpha _{2}={\frac {y_{2}}{d_{2}}}={\frac {5}{13}}$

$\alpha _{1}={\frac {y_{1}}{d_{1}}}-V(2,1)\alpha _{2}={\frac {14}{13}}$

(C.26)

C.2 METODO DELL’AUTOCORRELAZIONE

II metodo dell'autocorrelazione stima gli elementi della matrice F utilizzando esclusivamente campioni all'interno di una finestra di lunghezza “n”. Ciò equivale ad assumere implicitamente che il segnale sia nullo al di fuori di tale intervallo e quindi la traslazione di sequenze di campioni si ottengono con il suo completamento tramite campioni nulli. In tal modo il calcolo dell’autocorrelazione è eseguito su di un numero di campioni non nulli via via inferiore man mano che si calcolano valori di autocorrelazione di indice crescente. Gli elementi della matrice di autocorrelazione andrebbero quindi calcolati come

$R(i)={\frac {1}{n-1}}\sum _{k=1}^{n-i}x(k)\ x(k+i);\ 1\leq i\leq p$

(C.27)

In tal modo, però, al diminuire del denominatore, verrebbe esaltato lo scostamento della stima dal valore effettivo dell'autocorrelazione, scostamento dovuto al ridursi del numero di campioni utili. Si preferisce, quindi, non eseguire il calcolo degli elementi della matrice della autocorrelazione (non polarizzata) secondo tale definizione, ma calcolare gli elementi della matrice della autocorrelazione (polarizzata) come

$R(i)={\frac {1}{n}}\sum _{k=1}^{n-i}x(k)\ x(k+i);\ 1\leq i\leq p$

(C.28)

Ciò permette di attenuare il contributo degli elementi esterni in quanto,

al ridursi delle dimensioni delle sequenze, il valore del denominatore rimane costante. L’equazione normale deterministica in forma matriciale diventa

${\begin{bmatrix}R(0)&R(1)&R(2)&\ldots &R(p-l)\\R(l)&R(0)&R(1)&\ldots &R(p-2)\\R(2)&R(1)&R(0)&\ldots &R(p-3)\\\ldots &\ldots &\ldots &\ldots &\ldots \\R(p-l)&R(p-2)&R(p-3)&\ldots &R(0)\end{bmatrix}}{\begin{bmatrix}\alpha _{1}\\\alpha _{2}\\\alpha _{3}\\\ldots \\\alpha _{p}\end{bmatrix}}={\begin{bmatrix}R(1)\\R(2)\\R(3)\\\ldots \\R(p)\end{bmatrix}}$

$R\alpha =r$

(C.29)

con soluzione

$\alpha =R^{-1}r$

(C.30)

Es.: ripetiamo il calcolo dei coefficienti del predittore del secondo ordine

${\begin{bmatrix}\alpha _{1}\\\alpha _{2}\end{bmatrix}}={\begin{bmatrix}R(0)&R(1)\\R(1)&R(0)\end{bmatrix}}^{-1}{\begin{bmatrix}R(1)\\R(2)\end{bmatrix}}$

(C.31)

per lo stesso segnale precedentemente considerato. Gli elementi della matrice di autocorrelazione sono dati da

$x(n)=[2\ -1\ 0\ -1\ 2]\ R(0)={\frac {1}{5}}([2\ -1\ 0\ -1\ 2][2\ -1\ 0\ -1\ 2]^{T})=2$

$x(n+1)=[-1\ 0\ -1\ 2\ 0]\ R(l)={\frac {1}{5}}([2\ -1\ 0\ -1\ 2][-1\ 0\ 1\ 2\ 0]^{T})={\frac {4}{5}}$

$x(n+2)=[0\ -1\ 2\ 0\ 0]\ R(2)={\frac {1}{5}}([2\ -1\ 0\ -1\ 2][0\ -1\ 2\ 0\ 0]^{T})={\frac {1}{5}}$

(C.32)

I coefficienti del predittore sono, quindi, calcolati come

${\begin{bmatrix}\alpha _{1}\\\alpha _{2}\end{bmatrix}}={\begin{bmatrix}2&-{\frac {4}{5}}\\-{\frac {4}{5}}&2\end{bmatrix}}^{-1}{\begin{bmatrix}-{\frac {4}{5}}\\{\frac {1}{5}}\end{bmatrix}}={\begin{bmatrix}{\frac {25}{42}}&{\frac {5}{21}}\\{\frac {5}{21}}&{\frac {25}{42}}\end{bmatrix}}{\begin{bmatrix}-{\frac {4}{5}}\\{\frac {1}{5}}=\end{bmatrix}}={\begin{bmatrix}-{\frac {9}{21}}\\-{\frac {3}{42}}\end{bmatrix}}$

(C.33)

Anche in questo caso è possibile utilizzare la matrice delle osservazioni,

avente dimensioni [n + m -1, m], che assume la forma

$H^{T}={\begin{bmatrix}x(1)&x(2)&\ldots &x(p)&x(p+l)&\ldots &x(n)&0&\ldots &0\\0&x(1)&\ldots &x(p-1)&x(p)&\ldots &x(n-1)&x(n)&\ldots &0\\\ldots &\ldots &\ldots &\ldots &\ldots &\ldots &\ldots &\ldots &\ldots &\ldots \\0&0&\ldots &x(1)&x(2)&\ldots &x(n-p+1)&x(n-p)&\ldots &x(n)\end{bmatrix}}$

(C.34)

Dalla H si ricavano la matrice di autocorrelazione ed il vettore di cross-correlazione

$R={\frac {H^{T}H}{n}}$

$r={\frac {H^{T}x(n+1)}{n}}$

(C.35)

La soluzione diretta dell’equazione di Yule-Walker tramite l’inversione della matrice di auto correlazione può essere anche in questo caso evitata, dato che essa risulta essere una matrice di Toeplitz, cioè simmetrica con tutti gli elementi di una diagonale identici. Ciò permette di ricavare ricorsivamente i coefficienti di un predittore di ordine p in funzione dei coefficienti di predittori di ordine inferiore tramite una relazione ricorsiva del tipo

$\alpha _{k}^{(m)}=\alpha _{k}^{(m-1)}+\delta _{k}^{(m-1)}$

(C.36)

dove α_k^(m) indica il coefficiente k-esimo del predittore di ordine “m” e quindi al passo m-esimo dell’algoritmo. Le condizioni iniziali sono quelle di un predittore del primo ordine, che possono essere ricavate direttamente dall'equazione normale

$R(0)\alpha _{1}^{(1)}=R(1)$

$\alpha _{1}^{(1)}={\frac {R(1)}{R(0)}}$

(C.37)

Dall’espressione del minimo dell’MSE, per un predatore del primo ordine si ricava

$e^{(1)}=R(0)-\alpha _{1}^{(1)}R(1)=R(0)\left\{1-\left[\alpha _{1}^{(1)}\right]^{2}\right\}$

(C.38)

Considerando poi il problema di calcolare i coefficienti di un predittore di ordine m in funzione dei coefficienti di un predittore di ordine (m-1) è necessario poter risolvere il sistema

$\alpha ^{(m)}={\begin{bmatrix}\alpha _{1}^{(m)}\\\alpha _{2}^{(m)}\\\ldots \\\alpha _{m}^{(m)}\end{bmatrix}}={\begin{bmatrix}\alpha ^{(m-1)}\\0\end{bmatrix}}+{\begin{bmatrix}\delta ^{(m-1)}\\\hline \\k^{(m)}\end{bmatrix}}$

(C.39)

con δ^(m-1) e k^(m) da determinare. In tal modo si ottiene l’equazione aumentata del predittore

${\begin{bmatrix}{\begin{array}{c|c}R^{(m-1)}&r_{b}^{(m-1)}\\\hline r_{b}^{(m-1)}&R(0)\end{array}}\end{bmatrix}}\left\{{\begin{bmatrix}\alpha ^{(m-1)}\\0\end{bmatrix}}+{\begin{bmatrix}\delta ^{(m-1)}\\\hline \\k^{(m)}\end{bmatrix}}\right\}={\begin{bmatrix}r^{(m-1)}\\R(m)\end{bmatrix}}$

(C.40)

dove r_b^(m-1) rappresenta il vettore di cross-correlazione al passo m-1 r^(m-1) con gli elementi in ordine invertito. Da questa relazione matriciale derivano le due relazioni

$R^{(m-1)}\alpha ^{(m-1)}+R^{(m-1)}\delta ^{(m-1)}+r^{(m)}r_{b}^{(m-1)}=r^{(m-1)}$

$r_{b}^{(m-1)^{T}}\alpha ^{(m-1)}+r_{b}^{(m-1)^{T}}\delta ^{(m-1)}+k^{(m)}R(0)=R(m)$

(C.41)

Per l’equazione del predittore di ordine (m-1)

$R^{(m-1)}\alpha ^{(m-1)}=r^{(m-1)}$

(C.42)

dalla prima equazione si ottiene

$\delta ^{(m-1)}=-k^{(m)}R^{(m-1)^{-1}}r_{b}^{(m-1)}$

$\delta ^{(m-1)}=-k^{(m)}\alpha _{b}^{(m-1)}$

(C.43)

dove α_b^(m-1) rappresenta il vettore dei coefficienti del predittore di ordine m-1 α^(m-1) con gli elementi in ordine invertito. Sostituendo l’espressione di δ^(m-1) nell’equazione iniziale è possibile ricavare ricorsivamente i coefficienti del predittore come

$\alpha ^{(m)}=\alpha ^{(m-1)}-k^{(m)}\alpha _{b}^{(m-1)}$

$\alpha _{j}^{(m)}=\alpha _{j}^{(m-1)}-k^{(m)}\alpha _{m-j}^{(m-1)},\ 1\leq j\leq m-1$

(C.44)

In tale espressione la costante k^(m) risulta ancora incognita. Per poterla ricavare si sostituisce l’espressione di nella seconda componente dell’equazione normale aumentata, ottenendo

$r_{b}^{(m-1)^{T}}\alpha ^{(m)}=k^{(m)}r_{b}^{(m-1)^{T}}\alpha _{b}^{(m-1)}-r^{(m)}R(0)=R(m)$

$k^{(m)}={\frac {R(m)-r_{b}^{(m-1)^{T}}\alpha ^{(m-1)}}{R(0)-r_{b}^{(m-1)^{T}}\alpha _{b}^{(m-1)}}}={\frac {R(m)-r_{b}^{(m-1)^{T}}\alpha ^{(m-1)}}{R(0)-r_{b}^{(m-1)^{T}}R^{(m-1)^{-1}}}}r_{b}^{(m-1)}={\frac {R(m)-r_{b}^{(m-1)^{T}}\alpha ^{(m-1)}}{R(0)-\alpha ^{(m-1)^{T}}r^{(m-1)}}}$

(C.45)

Riconoscendo nel denominatore l’espressione dell’errore minimo al passo m-1, si ha

$k^{(m)}={\frac {R(m)-r_{b}^{(m-1)^{T}}\alpha ^{(m-1)}}{e^{(m-1)}}}={\frac {R(m)-\sum _{j=1}^{m-1}\alpha _{j}^{(m-1)}R(m-j)}{e^{(m-1)}}}$

(C.46)

Anche l’errore può essere calcolato ricorsivamente, dato che

$e^{(m)}=R(0)-\alpha ^{(m-1)^{T}}R_{m-1}=R(0)\sum _{k=1}^{m}\alpha _{k}^{(m)}R(k)$

(C.47)

Sostituendo l’espressione ricorsiva degli a k si ottiene

$e^{(m)}=R(0)-\sum _{k=1}^{m-1}\alpha _{k}^{(m-1)}R(k)-\alpha _{m}^{(m)}\left[R(m)-\sum _{k=1}^{m-1}\alpha _{m-k}^{(m-1)}R(k)\right]$

(C.48)

Confrontando l’espressione in parentesi quadre con il numeratore dell’espressione che fornisce k( m ), si ottiene

$e^{(m)}=e^{(m-1)}-\alpha _{m}^{(m)}k^{(m)}e^{(m-1)}$

(C.49)

ed essendo $k^{(m)}=\alpha _{m}^{(m)}$, si ha, infine

$e^{(m)}=\left[1-k^{(m)^{2}}\right]e^{(m-1)}$

(C.50)

Riepilogando, è possibile calcolare ricorsivamente tramite l’algoritmo di Levinson-Durbin i coefficienti di un predatore di ordine p in funzione dei coefficienti di predittori di ordine inferiore. Ciò si ottiene tramite i seguenti passi:

• si determina il valore dell’errore dalla relazione

$e^{(m)}=R(0)-\sum _{k=1}^{m}\alpha _{k}^{(m)}R(k)$

$con\ valore\ iniziale\ e^{(0)}=R(0)$

(C.51)

• si ricava il valore del coefficiente k^(m) dalla relazione

$k^{(m)}={\frac {R(m)-\sum _{j=1}^{m-1}\alpha _{j}^{(m-1)}R(m-j)}{e^{(m-1)}}}$

$con\ valore\ iniziale\ k^{(1)}={\frac {R(1)}{e^{(0)}}}$

(C.52)

• si determinano i coefficienti del predittore

$\alpha _{m}^{(m)}=k^{(m)}$

$\alpha _{j}^{(m)}=\alpha _{j}^{(m-1)}-k^{(m)}\alpha _{m-j}^{(m-1)},\ 1\leq j\leq m-1$

(C.53)

• si determina il valore aggiornato dell’errore

$e^{(m)}=\left[1-k^{(m)^{2}}\right]e^{(m-1)}$

(C.54)

Tale procedura si ripete per un numero di passi pari all'ordine del predittore.

Es.: ripetiamo il calcolo dei coefficienti del predittore per lo stesso segnale precedentemente considerato. L’errore al passo iniziale è

$e^{(0)}=R(0)=2$

(C.55)

Il ciclo per m = 1 porta al calcolo delle grandezze di un predittore del primo ordine

$k^{(1)}={\frac {R(1)}{e^{(0)}}}={\frac {2}{5}}$

$\alpha _{1}^{(1)}=k^{(1)}={\frac {2}{5}}$

$e^{(1)}=\left(1-k^{(1)^{2}}\right)e^{(0)}={\frac {42}{25}}$

(C.56)

Per m = 2 si ottengono le grandezze del predittore del secondo ordine

$k^{(2)}={\frac {R(2)-\alpha _{1}^{(1)}R(1)}{e^{(1)}}}=-{\frac {3}{42}}$

$\alpha _{2}^{(2)}=k^{(2)}=-{\frac {3}{42}}$

$\alpha _{1}^{(2)}=\alpha _{1}^{(1)}-k^{(2)}\alpha _{1}^{(1)}=-{\frac {9}{21}}$

(C.57)

C.3 METODI UTILIZZANTI STRUTTURE A TRALICCIO

Sia il metodo della covarianza che il metodo dell’autocorrelazione precedentemente esposti sono concettualmente organizzati su due fasi: una stima preliminare della matrice di autocorrelazione del segnale ed il successivo calcolo dei coefficienti del predittore. Il metodo che si vuole presentare in questo paragrafo (lattice) può essere visto come un raffinamento del metodo dell’autocorrelazione, nel quale tali due fasi sono fuse in un’unica procedura ricorsiva. Per descrivere l’algoritmo è necessario riprendere l’equazione del predittore

${\hat {x}}(n)=\sum _{k=1}^{p}\alpha _{k}\ x(n-k)$

(C.58)

e della sua funzione di trasferimento

$P(z)=\sum _{k=1}^{p}\alpha _{k}\ z^{-k}$

(C.59)

Definita la funzione d’errore del predittore di ordine “p” come

$e^{(p)}(n)=x(n)-{\hat {x}}(n)=x(n)-\sum _{k=1}^{p}\alpha _{k}^{(p)}x(n-k)$

(C.60)

è possibile definire la funzione di trasferimento del sistema d’errore. Si considera, cioè, il sistema che, dato come ingresso il segnale, produce in uscita la funzione d’errore. Tale funzione di trasferimento è pari a

$A^{(p)}(z)=1-\sum _{k=1}^{p}\alpha _{k}z^{-k}$

(C.61)

In tal modo è possibile esprimere l’errore di predizione in z come

$E^{(p)}(z)=A^{(p)}(z)\ X(z)$

(C.62)

Ora ci si può porre il problema di esprimere la A(z) di un predittore di ordine “m” in funzione della A(z) di un predittore di ordine Per tale problema si può utilizzare la trattazione svolta per l’algoritmo di Levinson-Durbin che fornisce un’espressione ricorsiva per i coefficienti del predittore a k. Sostituendo tale relazione nella definizione di A(z) si ottiene

$A^{(m)}(z)=1-\sum _{j=1}^{m}\alpha _{j}^{(m)}z^{-j}=1-\sum _{j=1}^{m}\left[\alpha _{j}^{(m-1)}+k^{m}\alpha _{m-j}^{(m-1)}\right]z^{-j}$

$=1-\sum _{j=1}^{m-1}\left[\alpha _{j}^{(m-1)}-k^{m}\alpha _{m-j}^{(m-1)}\right]z^{j}-\alpha _{m}^{(m)}z^{m}$

$=\left[1-\sum _{j=1}^{m-1}\alpha _{j}^{(m-1)}z^{j}\right]-\left[\alpha _{m}^{(m)}z^{-m}-\sum _{j=1}^{m-1}k^{m}\alpha _{m-j}^{(m-1)}z^{-j}\right]$

$A^{(m)}(z)=A^{(m-1)}(z)-k^{(m)}z^{m}A^{(m-1)}\left(z^{-1}\right)$

(C.63)

Sostituendo nell'espressione dell’errore si ottiene

$E^{(m)}(z)=A^{(m-1)}(z)X(z)-k^{(m)}z^{-m}A^{(m-1)}\left(z^{-1}\right)X(z)$

(C.64)

Il primo termine del secondo membro rappresenta la trasformata z della funzione d’errore di un predittore di ordine “m - 1”. Tale predittore stima il campione x(n) utilizzando i campioni [ x(n - m)... x(n -1) ] e verrà nel seguito indicato come predittore in avanti. Per il secondo termine si può dare un’interpretazione simile se si definisce

$B^{(m)}(z)=z^{-m}A^{(m)}\left(z^{-1}\right)X(z)$

(C.65)

che ha un’antitrasformata del tipo

$b^{(m)}(n)=x(n-m)-\sum _{k=1}^{m}\alpha _{k}^{(m)}x(n+k-m)$

(C.66)

Questo è interpretabile come l’errore di un predittore che tenta di stimare x(n - m) in funzione degli “m-1” campioni [x(n m 1)... x(n) ]. Tale predittore verrà nel seguito indicato come predittore all'indietro. Sostituendo la B(z) nella E(z) si ottiene

$E^{(m)}(z)=E^{(m-1)}(z)-k^{(m)}B^{(m-1)}\left(z^{-1}\right)$

(C.67)

che, antitrasformata fornisce

$e^{(m)}(n)=e^{(m-1)}(n)-k^{(m)}b^{(m-1)}(n-1)$

(C.68)

Tale espressione rappresenta una relazione ricorsiva della funzione d’errore della predizione in avanti. Sostituendo l’espressione ricorsiva della A(z) nella definizione della B(z) si ottiene

$B^{(m)}(z)=z^{-m}A^{(m-1)}\left(z^{-1}\right)X(z)-k^{(m)}A^{(m-1)}(z)X(z)$

$=z^{-1}B^{(m-1)}(z)X(z)-k^{(m)}E^{(m-1)}(z)$

(C.69)

che, antitrasformata, fornisce la relazione ricorsiva della funzione d’errore della predizione all'indietro

$b^{(m)}(n)=b^{(m-1)}(n-1)-k^{(m)}e^{(m-1)}(n)$

(C.70)

Se si traccia il diagramma di flusso di tale algoritmo si ottiene una struttura a traliccio (lattice), da cui il nome dell’algoritmo (fig. C.2). È quindi possibile il calcolo ricorsivo delle due funzioni d’errore. Le condizioni iniziali sono quelle per un predittore di ordine zero, cioè senza alcuna predizione e quindi

$e^{(0)}(n)=b^{(0)}(n)=x(n)$

(C.71)

In forma matriciale le relazioni ricorsive in z sono esprimibili come

${\begin{bmatrix}A^{(m)}(z)\\B^{(m)}(z)\end{bmatrix}}={\begin{bmatrix}1&k^{(m)}\\k^{(m)}z^{-1}&z^{-1}\end{bmatrix}}{\begin{bmatrix}A^{(m-1)}(z)\\B^{(m-1)}(z)\end{bmatrix}}$

(C.72)

Per il calcolo dei coefficienti k( m ) si impone la simultanea minimizzazione degli errori quadratici di predizione in avanti e all'indietro

${\begin{cases}E\left\{\left[e^{(m-1)}(n)\right]^{2}\right\}=\sum _{i=0}^{n-1}\left[e^{(m-1)}(i)\right]^{2}\\E\left\{\left[b^{(m-1)}(n)\right]^{2}\right\}=\sum _{i=0}^{n-1}\left[b^{(m-1)}(i)\right]^{2}\end{cases}}$

(C.73)

La minimizzazione, ottenuta annullando le derivate rispetto ai k

${\begin{cases}{\frac {\partial E\left\{\left[e^{(m)}(n)\right]^{2}\right\}}{\partial k}}=0\\{\frac {\partial E\left\{\left[b^{(m)}(n)\right]^{2}\right\}}{\partial k}}=0\end{cases}}$

${\begin{cases}k^{(m)}E\left\{\left[b^{(m-1)}(n)\right]^{2}\right\}=E\left\{e^{(m-1)}(n)b^{(m-1)}(n-1)\right\}\\k^{(m)}E\left\{\left[e^{(m-1)}(n)\right]^{2}\right\}=E\left\{e^{(m-1)}(n)b^{(m-1)}(n-1)\right\}\end{cases}}$

(C.74)

porta alla seguente espressione ricorsiva per il calcolo dei coefficienti

$k^{(m)}={\frac {\sum _{i=0}^{n-1}e^{(m-1)}(i)b^{(m-1)}(i-1)}{\sqrt {\sum _{i=0}^{n-1}\left[e^{(m-1)}(i)\right]^{2}\sum _{i=0}^{n-1}\left[b^{(m-1)}(i)\right]^{2}}}}$

(C.75)

Dato che tale espressione ha la forma di una cross-correlazione tra le funzioni di errore di predizione in avanti e all’indietro, i coefficienti k( m ) vengono detti PARCOR (PARtial CORrelation). Noti i coefficienti k( m ) è possibile determinare i coefficienti del predittore come già visto nel Levinson-Durbin (dal quale il presente algoritmo deriva)

$\alpha _{m}^{(m)}=k^{(m)}$

$\alpha _{j}^{(m)}=\alpha _{j}^{(m-1)}-k^{(m)}\alpha _{m-j}^{(m-1)}\ 1\leq j\leq m-1$

(C.76)

Per chiarire ulteriormente la struttura degli errori di predizione è opportuno analizzare con un esempio la loro struttura. Per un predittore di ordine zero risulta

${\begin{cases}e_{0}(n)=x(n)\\b_{0}(n)=x(n)\end{cases}}$

(C.77)

Fig. C.2 - Predittore di ordine p implementato tramite una struttura a traliccio.

Per un predittore di ordine uno la predizione in avanti utilizza solamente il campione che precede il corrente, mentre la predizione all'indietro utilizza il campione corrente per predire il precedente. Infatti risulta

${\begin{cases}e^{(1)}(n)=e^{(0)}(n)+k^{(1)}b^{(0)}(n-1)=x(n)+\alpha _{1}^{(1)}x(n-1)\\b^{(1)}(n)=k^{(1)}e^{(0)}(n)+b^{(0)}(n-1)=\alpha _{1}^{(1)}x(n)+x(n-1)\end{cases}}$

(C.78)

Analogamente, per un predittore di ordine due risulta

${\begin{cases}e^{(2)}(n)=e^{(2)}(n)+k^{(2)}b^{(2)}(n-1)=x(n)+\alpha _{1}^{(2)}x(n-1)+\alpha _{2}^{(2)}x(n-2)\\b^{(2)}(n)=k^{(2)}e^{(1)}(n)+b^{(1)}(n-1)=\alpha _{2}^{(2)}x(n)+\alpha _{1}^{(2)}x(n-1)+x(n-2)\end{cases}}$

(C.79)

e così via. Riassumendo, i passi per il metodo lattice sono i seguenti:

• si impostano le condizioni iniziali

$e^{(0)}(n)=b^{(0)}(n)=x(n)$

(C.80)

• si calcolano i coefficienti PARCOR

$k^{(m)}={\frac {\sum _{m=0}^{n-1}e^{(m-1)}(n)b^{(m-1)}(n-1)}{\sqrt {\sum _{m=0}^{n-1}\left[e^{(m-1)}(n)\right]^{2}\sum _{m=0}^{n-1}\left[b^{(m-1)}(n)\right]^{2}}}}$

(C.81)

• si calcolano i coefficienti del predittore

$\alpha _{m}^{(m)}=k^{(m)}$

$\alpha _{j}^{(m)}=\alpha _{j}^{(m-1)}-k^{(m)}\alpha _{m-j}^{(m-1)}\ 1\leq j\leq m-1$

(C.82)

• si aggiornano le funzioni di errore

$e^{(m)}(n)=e^{(m-1)}(n)+k^{(m)}b^{(m-1)}(n-1)$

$b^{(m)}(n)=b^{(m-1)}(n-1)-k^{(m)}e^{(m-1)}(n)$

(C.83)

con vettori che, inizialmente di dimensioni pari ad “n”, assumono dimensioni crescenti fino a “n + p”.

Es.: si ripete il calcolo del predittore del secondo ordine con il segnale precedentemente utilizzato. Le condizioni iniziali sono

$x(n)=[2\ -1\ 0\ -1\ 2]$

$e^{(0)}(n)=[2\ -1\ 0\ -1\ 2]$

$b^{(0)}(n)=[2\ -1\ 0\ -1\ 2]$

$b^{(0)}(n-1)=[0\ 2\ -1\ 0\ -1]$

(C.84)

Il coefficiente PARCOR del primo ordine è

$k^{(1)}={\frac {\sum e^{(0)}(n)b^{(0)}(n-1)}{\sqrt {\sum e^{(0)}(n)^{2}b^{(0)}(n)^{2}}}}={\frac {\sum [2\ -1\ 0\ -1\ 2][0\ 2\ -1\ 0\ -1]}{\sqrt {\sum [2\ -1\ 0\ -1\ 2]^{2}[2\ -1\ 0\ -1\ 2]^{2}}}}$

$={\frac {\sum [0\ -2\ 0\ 0\ 2]}{\sqrt {\sum [4\ 1\ 0\ 1\ 4][4\ 1\ 0\ 1\ 4]}}}=-{\frac {2}{5}}$

(C.85)

che porta alle nuove funzioni di errore

$e^{(1)}(n)=e^{(0)}(n)+k^{(1)}b^{(0)}(n-1)=[2\ -1\ 0\ -1\ 2]+{\frac {2}{5}}[0\ 2\ -1\ 0\ -1\ 2]$

$=\left[2\ {\frac {-1}{5}}\ {\frac {2}{5}}\ -1\ {\frac {8}{5}}\ {\frac {4}{5}}\right]$

$b^{(1)}(n)=b^{(0)}(n-1)-k^{(1)}e^{(0)}(n)=[0\ 2\ -1\ 0\ -1]+{\frac {2}{5}}[2\ -1\ 0\ 1\ 2]$

$=\left[{\frac {4}{5}}\ {\frac {8}{5}}\ 1\ {\frac {-2}{5}}\ {\frac {-1}{5}}\ 2\right]$

(C.86)

ed al coefficiente

$\alpha _{1}^{(1)}=k^{(1)}=-{\frac {2}{5}}$

(C.87)

Per il secondo ordine si ha

$k^{(2)}={\frac {\sum e^{(1)}(n)b^{(1)}(n-1)}{\sqrt {\sum e^{(1)}(n)^{2}b^{(1)}(n)^{2}}}}={\frac {\sum \left[2\ {\frac {-1}{5}}\ {\frac {-2}{5}}\ -1\ {\frac {8}{5}}\ {\frac {4}{5}}\right]\left[0\ {\frac {4}{5}}\ {\frac {8}{5}}\ -1\ {\frac {-2}{5}}\ {\frac {-1}{5}}\right]}{\sqrt {\sum \left[2\ {\frac {-1}{5}}\ {\frac {-2}{5}}\ -1\ {\frac {8}{5}}\ {\frac {4}{5}}\right]^{2}\left[{\frac {4}{5}}\ {\frac {8}{5}}\ -1\ {\frac {-2}{5}}\ {\frac {-1}{5}}\ 2\right]^{2}}}}=-{\frac {3}{42}}$

(C.88)

Si possono quindi calcolare i coefficienti del precettore del secondo ordine

$\alpha _{2}^{(2)}=k^{(2)}=-{\frac {3}{42}}$

$\alpha _{1}^{(2)}=\alpha _{1}^{(1)}-k^{(2)}\alpha _{1}^{(1)}=-{\frac {9}{21}}$

(C.89)

Si noti come i valori dei coefficienti coincidono con quelli calcolati con

il metodo dell’auto-correlazione

C.4 METODO RICORSIVO DI SCHÜR

Si osserva che l’algoritmo di Levinson-Durbin fa uso di una relazione ricorsiva per Terrore al passo m-esimo in funzione dell’errore e dei coefficienti PARCOR al passo (m-1) del tipo

$e^{(m)}=\left[1-\left(k^{(m)}\right)^{2}\right]e^{(m-1)}$

(C.90)

È possibile esprimere tale errore il funzione dell’autocorrelazione esplicitandola come

$e^{(1)}=\left[1-\left(k^{(1)}\right)^{2}\right]R(0)$

$e^{(2)}=\left[1-\left(k^{(1)}\right)^{2}\right]\left[1-\left(k^{(m)}\right)^{2}\right]R(0)=\left[1-\left(k^{(1)}\right)^{2}\right]^{2}R(0)$

$\cdots$

$e^{(m)}=\left[1-\left(k^{(1)}\right)^{2}\right]^{m}R(0)$

(C.91)

Nell'algoritmo di Schür si esprimono tutte le relazioni al generico passo m-esimo in funzione degli elementi della matrice di autocorrelazione R(i) (i=0,l,...,m) e dei parametri k⁽ⁱ⁾ (i=l,2,...,m). Per l’algoritmo di Schür si definiscono innanzitutto gli elementi

$G_{0}={\begin{bmatrix}0&R(l)&R(2)&\cdots &R(p)\\R(0)&R(l)&R(2)&\cdots &R(p)\end{bmatrix}}$

$G_{1}={\begin{bmatrix}0&R(l)&R(2)&\cdots &R(p)\\0&R(0)&R(l)&\cdots &R(p-1)\end{bmatrix}}$

$k^{(1)}={\frac {G_{1}(1,2)}{G_{1}(2,2)}}$

$G_{0}={\begin{bmatrix}1&-k^{(1)}\\k^{(1)}&1\end{bmatrix}}$

C.92

dove le dimensioni di G₀ sono pari a 2 * (p+1), con p ordine del predittore. Si esegue poi il prodotto V₁G₁ e si definisce G₂ come il risultato della precedente operazione con la seconda riga traslata di un posto verso destra. Si ripete il procedimento per p passi, ottenendo al generico passo m-esimo le seguenti relazioni ricorsive

$k^{(m)}={\frac {G_{m}(1,m+1)}{G_{m}(2,m+1)}}$

$V_{m}={\begin{bmatrix}1&-k^{(m)}\\k^{(m)}&1\end{bmatrix}}$

C.93

Seguendo passo-passo Io svolgimento dell’algoritmo di Schür, si può facilmente verificare che al passo m-esimo l’elemento G_m(2,m+1), ovvero il denominatore della relazione che definisce K^(m), corrisponde al termine e^(m-1) nell'algoritmo di Levinson-Durbin. Ad esempio, se implementiamo la stima dei parametri di un predittore di ordine 3 la matrice G 0 è definita come

$V_{m}={\begin{bmatrix}0&R(l)&R(2)&R(3)\\R(0)&R(1)&R(2)&R(3)\end{bmatrix}}$

C.94

Si ottiene poi:

$G_{1}(2,2)=R(0)=e^{(0)}$

$G_{2}(2,3)=R(0)-k^{(1)}R(l)=[1-(k^{(1)})^{2}]R(0)=e^{(1)}$

$G_{3}(2,4)=R(0)-k^{(1)}R(l)-k^{(2)}[R(2)-k^{(1)}R(l)]=R(0)[1-(k^{(1)})^{2}][l-(k^{(2)})^{2}]=e^{(2)}$

C.95

dove, al generico passo m-esimo, si è ricavato R(m) invertendo la relazione che definisce k^(m) nell'algoritmo di Levinson Durbin.

Si nota anche che al generico passo m-esimo dell’algoritmo di Schür, la matrice G_m presenta le prime m colonne nulle, il che può essere sfruttato per ridurre sia la complessità di calcolo che l’occupazione di memoria.

Infatti, la versione dell’algoritmo di Schifi proposta dallo standard GSM organizza le informazioni contenute nelle matrici G m in due vettori (K e P) ed evita, ad ogni passo, di aggiornarne quegli elementi che non apportano contributo ai passi successivi. In questo modo non si utilizzano relazioni matriciali; tutte le operazioni sono ridefinite sugli elementi di K e di P, vettori rispettivamente di dimensione (p-1) e (p+1). Al generico passo m, la determinazione di r(m) è sempre effettuata come rapporto tra il secondo elemento di P ed il primo; vengono aggiornati tutti gli elementi di P e di K tranne gli ultimi m.

Appendice D

RICHIAMI SU FILTRI NUMERICI

---

D.1 STRUTTURE PER MULTI RATE DSP

Per Multirate Digital Signal Processing si intende la realizzazione di sistemi discreti tali che il periodo T’dei segnali in uscita sia differente da quello T dei segnali in ingresso. Nel caso in cui la frequenza di campionamento del segnale in uscita sia maggiore di quella del segnale d’ingresso, è necessaria un’operazione di interpolazione per ricostruire il segnale x(t) mancante. Viceversa, nel caso in cui la frequenza del segnale in uscita sia minore di quella d’ingresso, è necessario decimare il segnale, eliminando campioni. È, però, necessario far precedere alla decimazione un’operazione di filtraggio che limiti la banda del segnale da ricampionare, al fine di evitare aliasing. Indicando con h(n) la risposta impulsiva dei filtri coinvolti nel cambiamento di frequenza (anti aliasing o di interpolazione), l’uscita è prodotta come

$y(m)=x(t)|_{t=mT'}=\sum _{n=-\infty }^{\infty }x(n)h(t-nT)|_{t=mT'}=\sum _{n=-\infty }^{\infty }x(n)h(mT'-nT)$

(D.1)

In un sistema discreto è conveniente che il rapporto tra i periodi di campionamento in ingresso ed in uscita sia riconducibile ad un rapporto tra interi del tipo

${\frac {T'}{T}}={\frac {M}{L}}$

(D.2)

Inoltre, data la linearità del filtraggio, è possibile considerare separatamente i due casi di interpolazione ( M = 1, L > 1 ) e di decimazione (M>1,L=1). Il caso più generale { M > 1,L> 1 } è poi risolvibile come cascata di una interpolazione seguita da una decimazione.

Nel caso della decimazione, il ricampionamento si ottiene semplicemente eliminando M-l campioni del segnale originario ogni M. Al fine di evitare aliasing, il segnale dovrebbe essere preventivamente limitato in banda tramite un filtro passa basso ideale con funzione di trasferimento

$H(\omega )={\begin{cases}1&|\omega |\leq {\frac {\pi }{M}}\\0&|\omega |>{\frac {\pi }{M}}\end{cases}}$

(D.3)

Indicando con h(n) la corrispondente risposta impulsiva, l’uscita si ottiene come

${\begin{cases}w(n)=\sum _{k=-\infty }^{\infty }h(k)\ x(n-k)\\y(m)=w(Mm)=\sum _{n=-\infty }^{\infty }h(n)\ x(mM-n)\end{cases}}$

(D.4)

Con una trasformazione di variabili, l’equazione precedente può essere riscritta come

$y(m)=\sum _{n=-\infty }^{\infty }h(mM-n)\ x(n)$

(D.5)

dalla quale si nota come i campioni derivanti dalla decimazione si ottengano traslando progressivamente la risposta impulsiva del filtro passa basso con un passo pari ad M campioni. Dal punto di vista della rappresentazione in frequenza, l’operazione di decimazione riespande lo spettro del segnale, limitato dal filtraggio nell’intervallo [-π/M, π/M], ad occupare tutto l’intervallo [-π, π].

Per eseguire, invece, un’interpolazione, è necessario introdurre L-l campioni ogni due campioni adiacenti del segnale originale. Un modo per ottenere questo risultato è quello di considerare una versione sovracampionata di un fattore L del segnale PAM corrispondente alla x(n), ottenuto introducendo dei campioni nulli nel segnale

$w(m)={\begin{cases}x\left({\frac {M}{L}}\right)&m=0,\pm L,\pm 2L,\ldots \\0&altrimenti\end{cases}}$

(D.6)

Dal punto di vista della rappresentazione in frequenza, ciò non altera lo spettro del segnale, coincidente con quello del corrispondente segnale PAM, il quale presenta repliche a partire dalla pulsazione ω = π/L. Per passare allo spettro della x(t) sovracampionata con un fattore L è necessario eliminare tali repliche tramite un filtraggio passa basso, con frequenza di taglio corrispondente a ω = π/L. Nel tempo, tale filtraggio si traduce in un’interpolazione che genera i campioni mancanti nella x(n). Per mantenere nel segnale sovracampionato l’ampiezza del segnale originario, il filtro di interpolazione deve avere un guadagno non unitario, ma pari ad L. Dato che, nel filtraggio, la risposta impulsiva h(n) si trova ad essere sistematicamente combinata con campioni nulli della w(n), la sommatoria di convoluzione può essere semplificata come

$y(m)=\sum _{k=-\infty }^{\infty }w(k)\ h(m-k)=\sum _{n=-\infty }^{\infty }x(k)h(m-kL)$

(D.7)

Analizzando quali campioni della h(n) vengono interessati dalla sommatoria, si nota che vengono utilizzati periodicamente L suoi sottoinsiemi g_m, tra loro scalati di un campione (fig. D.l). Tale sfasamento corrisponde ad un offset temporale che che è una frazione del periodo T pari a

$\delta _{m}=\left[{\frac {m}{L}}\right]_{mod\ L}$

(D.8)

E quindi possibile riscrivere l’equazione che fornisce la y(m) come l’uscita di un filtro tempo variante con risposta impulsiva g m, tale che

${\begin{cases}g_{m}=h\left[(n+\delta _{m})T\right]\\y(m)=\sum _{n=-\infty }^{\infty }g_{m}(n)x\left({\frac {m}{L}}-n\right)\end{cases}}$

(D.9)

dove come indice della x deve considerarsi, ovviamente, un intero.

Analogamente, è possibile considerare i campioni della y(m) come prodotti da un banco di L filtri in parallelo, ciascuno con una differente risposta impulsiva g_m(n). L’uscita dell’interpolatore è poi ottenuto selezionando periodicamente l’uscita di uno degli L filtri che operano tra loro sfasati di un colpo di clock.

Fig. D.1 -Insieme di campioni della risposta impulsiva coinvolti in un’operazione di inteipolazione. Passando al problema dell’implementazione, è necessario considerare che i filtri utilizzati nella conversione di frequenza di campionamento sono tipicamente FIR a lunghezza finita N (fig. D.2), per cui l’uscita è ottenuta come

$y(n)=\sum _{k=0}^{N-1}h(k)\ x(n-k)$

(D.10)

Fig. D.2 - Strutture di filtri FIR per decimazione. Considerando per il momento il problema della decimazione, un inconveniente che si incontra implementando direttamente tale filtro è che esso lavora alla frequenza del segnale sovracampionato, mentre sarebbe auspicabile che lavorasse alla frequenza del segnale decimato. Per ottenere ciò è possibile anticipare l’operazione di decimazione prima del prodotto con i coefficienti del filtro come

$y(m)=\sum _{k=0}^{N-l}h(k)\ x(Mm-k)$

(D.11)

ottenendo un filtro con clock pari a quello dell’uscita. Imponendo poi che i filtri siano a fase lineare e quindi con risposta impulsiva risulta simmetrica, la lunghezza del filtro può essere dimezzata in quanto

$h(k)=h(N-l-k)\to y(n)=\sum _{k=0}^{N-l}h(k)\left\{x(n-k)+x[n-(N-1-k)]\right\}$

(D.12)

Strutture analoghe a quelle presentate nel caso della decimazione possono essere realizzate nel caso di interpolazione. In ogni caso, volendo avere filtri sufficientemente selettivi, le loro dimensioni risultano essere notevoli. L’efficienza computazionale può essere incrementata riducendo la lunghezza del filtro tramite filtri polifase (fig. D.3). È stato mostrato come i campioni prodotti da un interpolatore possano essere pensati come l’uscita di uno di L filtri in parallelo, ciascuno con risposta impulsiva g m (n). Dato che la risposta impulsiva di ciascun ramo del filtro coincide con una versione decimata della h(n), che è diversa da zero solo negli istanti L x T ed è nulla altrimenti, è possibile eliminare i contributi di tali campioni nulli considerando un insieme di L filtri di dimensioni N/L tali che

$p_{k}^{j}(n)=h(k+nL);{\begin{cases}k=0,1,\ldots ,L-1\\n=0,l,\ldots ,{\frac {N}{L}}-1\end{cases}}$

(D.13)

Ovviamente è opportuno che la lunghezza N del filtro prototipo sia multipla di L. Con una struttura polifase, per ogni campione dell’ingresso vengono prodotti L campioni d’uscita, uno per ogni ramo del filtro. L-l di tali campioni costituiscono l’interpolazione tra due campioni adiacenti dell’ingresso, necessari per il richiesto incremento della frequenza di campionamento.

Fig. D.3 - Strutture di filtri FIR polifase.

Per quanto riguarda la funzione di trasferimento di ciascuno degli L rami del filtro polifase, si ricorda che essi sono stati ottenuti come decimazione di un fattore L di un filtro passa basso ideale con frequenza di taglio corrispondente a ω = π/L. Dato che la decimazione riespande lo spettro di un fattore L, ciascun ramo del filtro polifase risulterà con funzione di trasferimento piatta nell'intervallo [-π, πJ. L’unica differenza tra i differenti rami è quindi data dalla fase, da cui il nome del filtro.

Tipicamente la struttura di un filtro polifase non può essere semplificata ulteriormente sfruttando eventuali simmetrie. Infatti, pur utilizzando una h(n) simmetrica, le L risposte impulsive dei rami del filtro polifase non lo sono. D’altra parte, l’indipendenza tra i vari rami del filtro, ne rende possibile un'implementazione hardware parallela. Per quanto riguarda la decimazione, la relazione

$y(m)=\sum _{k=-\infty }^{\infty }h(k)\ x(mM-k)$

(D.14)

può essere riscritta come

$y(m)=\sum _{k=0}^{M-1}\sum _{r=-\infty }^{\infty }h(rM+k)\ x[(m-r)M-k]$

(D.15)

Definendo i coefficienti del filtro polifase di decimazione come

$p_{k}^{d}(n)=h(nM+k);\ k=0,1,\ldots ,M-1$

(D.16)

si ottiene

$y(m)=\sum _{k=0}^{M-1}\sum _{r=-\infty }^{\infty }p_{k}^{d}(r)\ x_{k}(m-r)$

(D.17)

dove x_k(n) = x(nM - k) è l’ingresso del ramo k-esimo del filtro. Di tali filtri è possibile dare un’interpretazione tramite un commutatore che preleva campioni nel caso della interpolazione e distribuisce campioni nel caso della decimazione ruotando in senso antiorario. Ciò vuol dire, ad esempio, che nella decimazione il primo campione ingresso viene assegnato al ramo di ordine maggiore e così via. È possibile realizzare strutture con senso di rotazione del commutatore orario definendo i filtri polifase come

${\begin{cases}p_{k}^{j}(n)=h(nL-k);\ k=0,1,\ldots ,L-l\\p_{k}^{d}(n)=h(nM-k);\ k=0,1,\ldots ,M-l\end{cases}}$

(D.18)

Un’applicazione importante di multirate DSP si ha nell'elaborazione disegnali passa banda, cioè di segnali con una banda limitata ω_Δ, ma centrati in ω₀ ≠ 0 (fig. D.4). Tale esigenza si incontra, ad esempio, in un banco di filtri: il fine è quello di ricavare il corrispondente segnale in banda base X(m), cioè il segnale che, a parità di banda, risulta centrato nell'origine. Ciò è concettualmente ottenibile innanzitutto traslando il segnale passa banda nell'origine tramite una modulazione e quindi decimandolo dopo il necessario filtraggio antialiasing tramite un passa basso con frequenza di taglio corrispondente a ω_Δ/2. Indicando con h(n) la risposta impulsiva di tale filtro, si ottiene

$X(m)=\sum _{n=-\infty }^{\infty }h(mM-n)\left[x(n)\ e^{-jn\omega _{0}}\right]$

D.19

dove M è il fattore di decimazione che si ricava dal rapporto tra la frequenza centrale del segnale e la sua banda. Per ricostruire il segnale originale dalla sua versione decimata, invece, è innanzitutto necessario interpolarlo per poi modularlo, riportandolo su ω₀. Indicandolo con f(n) la risposta impulsiva del filtro di interpolazione, si ottiene

$x(n)=e^{jn\omega _{0}}\left[\sum _{m=-\infty }^{\infty }X_{k}(m)f(n-mM)\right]$

D.20

Fig. D.4 - Differenti modulazioni per segnali passabanda. Dato che il segnale d’uscita risulta essere complessa, a causa del prodotto con gli esponenziali, tale tecnica è detta di modulazione complessa (fig. D.5). La struttura del sistema di decimazione/interpolazione può essere anche differente da quella presentata, portando le modulazioni come ultima operazione nella decimazione e come prima della interpolazione. Con tale soluzione, però, è necessario l’uso di filtri passa banda centrati sul segnale stesso, al posto dei passa basso del caso precedente.

Fig. D.5 - Modulazione complessa ed a banda laterale unica.

Analiticamente l’espressione del decimatore può essere ottenuta dalla precedente moltiplicandola per $e^{jmM\omega _{0}}\cdot e^{-jmM\omega _{0}}$ e riarrangiando i termini come

$X(m)=e^{-jmM\omega _{0}}\sum _{n=-\infty }^{\infty }\left[h(mM-n)\ e^{j(mM-n)\omega _{0}}\right]x(n)$

(D.21)

dove il termine tra parentesi quadre corrisponde alla risposta impulsiva del filtro passa banda di ampiezza ω_Δ centrato su ω₀. Analogamente, per l’interpolazione

$x(n)=\sum _{m=-\infty }^{\infty }\left[X(m)e^{jmM\omega _{0}}\right]\left[f(n-mM)e^{j(n-mM)\omega _{0}}\right]$

(D.22)

Le relazioni presentate forniscono un mappaggio del segnale passa banda in un segnale complesso in banda base nell'intervallo [- ω_Δ/2, ω_Δ/2] e viceversa. Per applicazioni di codifica del segnale, invece, è più interessante considerare un mappaggio nell'intervallo [0, ω_Δ ], in grado di fornire un segnale reale. Per tale mappaggio, noto come modulazione SSB per l’analogia con le analoghe tecniche trasmissive, al segnale ottenuto dalla modulazione complessa viene applicata una traslazione in frequenza di metà larghezza di banda, tramite il prodotto con e e^jnω_Δ/2.

D.2 BANCHI DI FILTRI POLIFASE

La funzione di un banco di filtri può essere sia quella di separare un segnale nelle componenti che fanno capo a differenti bande spettrali (filtro di analisi), sia quella di ricombinare i contributi di bande differenti in un unico segnale (filtro di sintesi) (fig. D.6). Per tali operazioni è necessario considerare i differenti segnali passa banda relativi a ciascuna sottobanda, applicando ad essi i risultati ottenuti in Appendice D.l.

Per quanto ci riguarda, il caso di maggiore interesse è quello di un banco di filtri in grado di generare K bande uniformemente spaziate. In tal caso, la sottobanda relativa alle frequenze inferiori risulta modulata SSB. Tale disposizione delle sottobande è normalmente indicata come even-stacking, a differenza del caso in cui, pur avendo bande uniformi, la prima sottobanda è centrata nell'origine, nel qual caso si parla di struttura odd-stacking del banco di filtri.

Considerando modulazioni complesse, il segnale relativo alla banda k-esima X k (m) è quindi ottenibile come

${\begin{cases}X_{k}(m)=\sum _{n=-\infty }^{\infty }h(mM-n)\ x(n)\ e^{-jn\omega _{k}}\\\omega _{k}={\frac {2k\pi }{k}}\end{cases}}$

(D.23)

mentre per il filtro di sintesi, indicando con ${\hat {X}}_{k}$ il risultato dell’elaborazione eseguita sul segnale relativo alla banda k-esixna, si ricava

${\hat {x}}(n)={\frac {1}{K}}\sum _{k=0}^{K-1}{\hat {x}}_{k}(n)={\frac {1}{K}}\sum _{k=0}^{K-1}e^{jn\omega _{k}}\left[\sum _{m=-\infty }^{\infty }{\hat {X}}_{k}(m)\ f(n-nM)\right]$

(D.24)

dove h(m) è la risposta impulsiva del filtro passa basso di analisi e f(m) quella utilizzata nella sintesi. Riscrivendo il filtro di sintesi come

${\hat {x}}(n)=\sum _{m=-\infty }^{\infty }f(n-nM)\left[{\frac {1}{K}}\sum _{k=0}^{K-1}e^{jn\omega _{k}}\ {\hat {X}}_{k}(m)\right]$

(D.25)

si vede come le uscite dei banchi di analisi e sintesi siano ottenibili tramite trasformata discreta di Fourier, da cui il nome di DFT filter banks utilizzati nel caso di banchi con bande uniformemente spaziate. Il massimo valore ammissibile per il fattore di decimazione M coincide con il numero di sottobande K, nel qual caso si parla si bande critically-sampled.

Fig. D.6 - Struttura del sistema di analisi e sintesi. Utilizzando filtri passa banda complessi la struttura cambia, con le modulazioni che diventano l’ultima operazione nell’analisi e la prima della sintesi ed ottenendo

$X_{k}(m)=e^{-jmM\omega _{k}}\sum _{n=-\infty }^{\infty }\left[h(mM-n)e^{j(mM-n)\omega _{k}}\right]x(n)$

$x(n)={\frac {1}{K}}\sum _{k=0}^{K-1}\left\{\sum _{m=-\infty }^{\infty }\left[X_{k}(m)e^{jmM\omega _{k}}\right]\left[f(n-mM)e^{j(n-mM)\omega _{k}}\right]\right\}$

(D.26)

Nel caso di K bande uniformi, quindi, le funzioni di trasferimento di ciascun filtro passa banda sono derivate per traslazione di un unica funzione prototipo passa basso H(ffl), in modo tale che per la banda k-esima si abbia

${\begin{cases}H_{k}(\omega )=H_{k}(\omega -\omega _{k})\\h_{k}(n)=h(n)\cdot e^{jn\omega _{k}}\\\omega _{k}={\frac {2k\pi }{K}}\end{cases}}$

(D.27)

Un miglioramento dal punto di vista dell’efficienza computazionale si ha utilizzando filtri polifase (fig. D7). In tal caso, la risposta impulsiva del ramo di ordine p per la banda k-esima ed il suo ingresso sono pari a

${\begin{cases}P_{\rho ,k}(m)=h_{k}(mM-\rho )\\x_{\rho }(m)=x(mM-\rho )\end{cases}};\ \rho =,\ldots ,K-l$

(D.28)

Secondo tale impostazione sembrerebbe che i filtri polifase da utilizzare sarebbero in numero pari al prodotto del numero di sottobande per il fattore di decimazione. Fortunatamente, però, sfrattando la relazione

$h_{k}(n)=h(n)\cdot e^{jn{\frac {2k\pi }{K}}}=h(n)\cdot e^{jn{\frac {2k\pi }{M}}}$

(D.29)

si ricava

$p_{\rho ,k}(m)=h(mM-\rho )\cdot e^{j(mM-\rho ){\frac {2k\pi }{M}}}=h(mM-\rho )\cdot e^{-j{\frac {2k\pi }{M}}}$

(D.30)

Se si definisce il p-esimo ramo dell’equivalente filtro polifase passa basso come

$P_{\rho }(m)=h(mM-\rho )$

(D.31)

allora è possibile riscrivere il filtro polifase della banda k-esima come

$P_{\rho ,k}(m)=p_{\rho }(m)\cdot e^{-j{\frac {2k\rho \pi }{M}}}$

(D.32)

Sostituendo nelle equazioni precedenti ed indicando con q_ρ(m) il ρ-esimo ramo del polifase passa basso equivalente del filtro di sintesi, si ottiene

${\begin{cases}X_{k}(m)=\sum _{\rho =0}^{M-1}\left[\sum _{r=-\infty }^{\infty }p_{\rho }(r)\ x_{\rho }(m-r)\right]e^{-j{\frac {2k\rho \pi }{M}}}\\{\hat {x}}_{k}(n)=\sum _{m=-\infty }^{\infty }q_{\rho }(n-m)\left[{\frac {1}{M}}\sum _{k=0}^{M-1}{\hat {X}}_{k}(m)\ e^{j{\frac {2k\rho \pi }{M}}}\right]\end{cases}}$

(D.33)

Fig. D.7 - Banco di filtri polifase.

Dall'analisi di tali relazioni, si nota come l’efficienza computazionale nasca, innanzitutto, dal poter condividere il filtro pp(m) tra le differenti bande, Inoltre, si nota come l’uscita del filtro di analisi sia ottenibile come DFT delle uscite dei singoli filtri polifase e che il filtro di sintesi si ottiene applicando i filtri polifase alla DFT inversa degli X_k. Un ulteriore incremento di efficienza computazionale nasce, quindi, dall'uso di FFT. Un problema che si incontra nel filtraggio passa banda è la presenza di aliasing, introdotto dalla decimazione, a seguito dell’uso di filtri non ideali. In tal modo l’energia in ciascuna sottobanda è influenzata dal contributo delle bande adiacenti. Questo fenomeno è tanto più sensibile per quelle bande a bassa energia, per le quali il contributo dovuto alla non idealità dei filtri può essere predominante. La realizzazione ottimale del banco dei filtri è, quindi, un punto chiave e la scelta dei filtri h(n) ed f(n) deve essere tale da minimizzare l’aliasing. È possibile analizzare tale fenomeno sia nel dominio del tempo o in frequenza.

Analizzando l’aliasing nel dominio del tempo, si vuole vedere sotto quali condizioni ne è possibile 1*eliminazione. Imporre che x(n) = x(n) è possibile sostituendo l’espressione del filtro di analisi in quello di sintesi, ottenendo

${\hat {x}}(n)=\sum _{m=-\infty }^{\infty }f(n-nM)\left[{\frac {1}{K}}\sum _{k=0}^{K-1}e^{jn\omega _{k}}\sum _{r=-\infty }^{\infty }h(mM-r)x(r)e^{-jr\omega _{k}}\right]$

$=\sum _{m=-\infty }^{\infty }\sum _{r=-\infty }^{\infty }x(r)\ f(n-nM)\ h(mM-r)\left[{\frac {1}{K}}\sum _{k=0}^{K-1}e^{j(n-r)\omega _{k}}\right]$

(D.34)

Dato che l’esponenziale è pari ad 1 per r = n - sK, con s intero, l’espressione precedente può essere riscritta come

${\hat {x}}(n)=\sum _{s=-\infty }^{\infty }x(n-sK)\sum _{m=-\infty }^{\infty }f(n-nM)h(mM-n+sK)$

(D.35)

affinché sia valida l’uguaglianza x(n) = x(n), la seconda sommatoria deve annullarsi per s * 0 e deve valere uno per s = 0. Indicando le versioni traslate delle risposte impulsive come

${\begin{cases}h_{m}(n)=h(n-mM)\\f_{m}(n)=f(n-mM)\end{cases}}$

(D.36)

la relazione precedente si riscrive come

${\hat {x}}(n)=x(n)\sum _{m=-\infty }^{\infty }f_{m}(n)\ h_{m}(-n)+\sum _{s=-\infty ,s\neq 0}^{\infty }x(n-sK)\sum _{m=-\infty }^{\infty }f_{m}(n)\ h_{m}(sK-n)$

(D.37)

Questa equazione ci permette di definire quali siano le condizioni di perfetta

ricostruzione del segnale. Affinché sia valida l'uguaglianza x(n) = x(n),dalla prima sommatoria si ricava che

$\sum _{n=-\infty }^{\infty }f_{m}(n)\ h_{m}(-n)=1$

(D.38)

il che indica che la somma del prodotto delle risposte impulsive traslate deve essere unitaria per qualsiasi valore di n. Dalla seconda sommatoria, che è necessario annullare, si comprende come sia necessario annullare le infinite repliche poste ad intervalli di K campioni, introdotte dalla decimazione. Ciò è ottenibile imponendo

$\sum _{m=-\infty }^{\infty }f_{m}(n)\ h_{m}(sK-n)=0;\ s\neq 0$

(D.39)

il che si raggiunge solo con filtri di lunghezza inferiore a K. Dato che la lunghezza dei filtri risultante da tale vincolo sarebbe estremamente modesta, essi non risulterebbero soddisfacenti dal punto di vista dell’analisi in frequenza in quanto poco selettivi. Passando, infatti, all’analisi in frequenza, è conveniente considerare un modello a filtri passa banda complessi, semplificato eliminando le modulazioni terminali (che si elidono) e sostituendo la decimazione/interpolazione con un modulatore che moltiplica il segnale per la sequenza

$s(n)={\begin{cases}1&n=\pm mM,\ m\ intero\\0&altrimenti\end{cases}}$

(D.40)

La rappresentazione in frequenza del legame ingresso/uscita che si ottiene è pari a

${\hat {X}}(\omega )=X(\omega ){\frac {1}{K}}\sum _{k=0}^{K-1}F_{k}(\omega )\ H_{k}(\omega )+\sum _{l=0}^{K-1}X\left(\omega -{\frac {2l\pi }{M}}\right){\frac {1}{K}}\sum _{k=0}^{K-1}F_{k}(\omega )\ H_{k}\left(\omega -{\frac {2l\pi }{M}}\right)$

(D.41)

dove

${\begin{cases}H_{k}(\omega )=H(\omega -\omega _{k})\\Fk(\omega )=F(\omega -\omega _{k})\end{cases}}$

(D.42)

sono i filtri di analisi e sintesi della banda k-esima. Mentre il primo termine di tale equazione implica che la funzione di trasferimento complessiva sia unitaria, il secondo impone funzioni passa basso ideali, approssimabili solo con filtri di dimensione elevata. In conclusione, non è possibile raggiungere compietamente la cancellazione dell’aliasing né nel dominio del tempo né in frequenza.

Tutte le relazioni precedenti sono state ricavate per una struttura even-stacking del banco di filtri. Per un’organizzazione odd-stacking deve risultare

$\omega _{k}={\frac {2\pi }{K}}\left(k+{\frac {1}{2}}\right)$

(D.43)

Si è dimostrato come sia il filtro di analisi che quello di sintesi sono legati ad un mappaggio tramite DFT dal dominio del tempo alla frequenza e viceversa. Il mappaggio tramite DFT, però, è sempre tale da legate l’origine temporale n = 0 a quella in frequenza co = 0. Volendo permettere il mappaggio tra due qualsiasi origini nei due domini, è possibile definire una trasformata discreta di Fourier generalizzata come

${\begin{cases}X_{k}^{GDFT}=\sum _{n=0}^{K-1}x(n)\ e^{-j(k+k_{0})(n+n_{0}){\frac {2\pi }{K}}}\\x(n)={\frac {1}{K}}\sum _{n=0}^{K-1}X^{GDFT}\ e^{j(k+k_{0})(n+n_{0}){\frac {2\pi }{K}}}\end{cases}}$

(D.44)

dove no e ko sono le nuove origini nel tempo ed in frequenza. Applicando tale generalizzazione a banchi di filtri di analisi e sintési che, nel caso di struttura odd-stacking, sono caratterizzati ko = 1/2 , si ottiene

${\begin{cases}X_{k}^{GDFT}(m)=\sum _{n=-\infty }^{\infty }h(mM-n)\ x(n)\ e^{-j(n+n_{0})\left(k+{\frac {1}{2}}\right){\frac {2\pi }{K}}}\\x(n)=\sum _{m=-\infty }^{\infty }f(n-mM){\frac {1}{K}}\sum _{k=0}^{K-1}X_{k}^{GDFT}(m)\ e^{j(n+n_{0})\left(k+{\frac {1}{2}}\right){\frac {2\pi }{K}}}\end{cases}}$

(D.45)

Anche per tale banco, la struttura può evolvere verso strutture con filtri

passa banda complessi o polifase. Inoltre, al fine di ottenere segnali reali, è possibile trasformarlo applicando una modulazione SSB

$X_{k}^{SBB}(m)=Re\left[X_{k}^{GDFT}(m)e^{j{\frac {mM}{2}}\omega _{\Delta }}\right]$

(D.46)

da cui, posto M = K/2, deriva

$X_{k}(m)=cos\left({\frac {m\pi }{2}}\right)\sum _{n=-\infty }^{\infty }x(n)\cdot h\left(P-1+{\frac {mK}{2}}-n\right)\cdot cos\left[{\frac {2\pi }{K}}\cdot \left(k+{\frac {1}{2}}\right)\cdot (n+n_{0})\right]$

$+sin\left({\frac {m\pi }{2}}\right)\sum _{n=-\infty }^{\infty }x(n)\cdot h\left(P-1+{\frac {mK}{2}}-n\right)\cdot sin\left[{\frac {2\pi }{K}}\cdot \left(k+{\frac {1}{2}}\right)\cdot (n+n_{0})\right]$

(D.47)

Tale relazione può essere formalmente semplificata in

${\begin{cases}X_{k}(m)=cos(\pi mk)\sum _{n=-\infty }^{\infty }x(n)\cdot h(mM-n)\cdot cos\left[{\frac {2\pi }{K}}\left(k+{\frac {1}{2}}\right)(mM-n+n_{0})\right]\\{\hat {x}}_{k}(m)=\sum _{m=-\infty }^{\infty }\sum _{k=0}^{{\frac {K}{2}}-1}cos(\pi mk)\ {\hat {X}}_{k}(m)\cdot f(n-mM)\cdot cos\left[{\frac {2\pi }{K}}\left(k+{\frac {1}{2}}\right)(mM-n+n_{0})\right]\end{cases}}$

(D.48)

D.3 BANCHI DI FILTRI QMF

Una tecnica differente per la realizzazione del banco è quella dei Quadrature-Mirror Filters (QMF) (fig. D.8). Un banco QMF permette la divisione dello spettro del segnale in due sotto-bande. Ciò è ottenuto tramite un filtraggio passa-alto abbinato ad un filtraggio complementare passa-basso. Banchi con un numero maggiore di due sottobande si ottengono ripetendo il filtraggio su entrambi (per sottobande uniformi) o uno solo dei segnali d’uscita (fig. D.9). La struttura risultante è ad albero, bilanciato o meno. Al fine di minimizzare l’effetto dell’aliasing, la somma delle due funzioni di trasferimento deve essere piatta. Se si indica con H₁(e^jωT) e H₁(e^jωT) la Fig. D.8 - Diagramma di un banco QMF.

funzione di trasferimento del filtro passa alto e passa basso di analisi, i segnali da loro prodotti sono pari a

$X_{0}(\omega )={\frac {1}{2}}\left[X\left({\frac {\omega }{2}}\right)H_{0}\left({\frac {\omega }{2}}\right)+X\left({\frac {\omega +2\pi }{2}}\right)H_{0}\left({\frac {\omega +2\pi }{2}}\right)\right]$

$X_{1}(\omega )={\frac {1}{2}}\left[X\left({\frac {\omega }{2}}\right)H_{1}\left({\frac {\omega }{2}}\right)+X\left({\frac {\omega +2\pi }{2}}\right)H_{1}\left({\frac {\omega +2\pi }{2}}\right)\right]$

(D.49)

Indicando con F₀(ω) e F₁(ω) la funzione di trasferimento del filtro passa basso e passa alto di sintesi, in fase di ricostruzione si ottiene

${\hat {X}}(\omega )=X_{0}(2\omega )F_{0}(\omega )+X_{1}(2\omega )F_{1}(\omega )$

(D.50)

Sostituendo in tale espressione quella dei segnali prodotti da ciascuna sottobanda si ottiene

${\hat {X}}(\omega )={\frac {1}{2}}\left[H_{0}(\omega )F_{0}(\omega )+H_{1}(\omega )F_{1}(\omega )\right]X(\omega )+$

$+{\frac {1}{2}}\left[H_{0}(\omega +\pi )F_{0}(\omega )+H_{1}(\omega +\pi )F_{1}(\omega )\right]X(\omega +\pi )$

(D.51)

Imponendo che si annulli l’effetto dell’aliasing è necessario annullare il secondo termine della precedente equazione, ottenendo

$H_{0}(\omega +\pi )F_{0}(\omega )+H_{1}(\omega +\pi )F_{1}(\omega )=0$

(D.52)

Fig. D.9 - Filtraggio passa banda tramite cascata di filtri complementari.

Per trovare le condizioni tali da soddisfare tale vincolo è necessario ridurre le variabili del problema. A tal fine si imporre innanzitutto che il filtro passa altro sia una replica traslata del passa basso di analisi

${\begin{cases}H_{0}(\omega )=H(\omega )\\H_{1}(\omega )=H(\omega +\pi )\end{cases}}$

(D.53)

da cui deriva

$H(\omega +\pi )F_{0}(\omega )+H(\omega )F_{1}(\omega )=0\to {\begin{cases}F_{0}(\omega )=2H(\omega )\\F_{1}(\omega )=2H(\omega +\pi )\end{cases}}$

(D.54)

dove il fattore di scala 2 serve per rendere unitaria la funzione di trasferimento complessiva. Infatti

${\hat {X}}(\omega )=[H^{2}(\omega )-H^{2}(\omega +\pi )]X(\omega )$

(D.55)

L’esser riusciti ad eliminare il secondo termine dell’equazione precedente vuol dire che l’aliasing introdotto nel filtro di sintesi da parte della decimazione è cancellato dalle repliche introdotte dall’interpolazione nel filtro di sintesi. Nel dominio del tempo, le risposte impulsive dei filtri sono date da

${\begin{cases}h_{0}(n)=h(n)\\h_{1}(n)=(-1)^{n}h(n)\end{cases}};{\begin{cases}g_{0}(n)=2h(n)\\g_{1}(n)=-2(-1)^{n}h(n)\end{cases}}$

(D.56)

Rimane una sola incognita, data da H(co), per la determinazione della quale si possono seguire differenti criteri [Cro83]. Una classe di filtri particolarmente diffusa è quella dei filtri FIR a fase lineare. Il vantaggio dell’assenza di distorsioni di fase risulta particolarmente utile nel caso di connessioni in cascata, eliminando la necessità di equalizzazioni ad ogni stadio. La risposta impulsiva di un filtro FIR a fase lineare risulta essere simmetrica, per cui

$h(n)=h(N-1-n);\ n=0,1,\ldots ,N-l$

(D.57)

e la funzione di trasferimento è data da una funzione reale ed un termine di fase lineare

$H(\omega )=H_{r}(\omega )e^{-j\omega {\frac {N-1}{2}}}$

(D.58)

La funzione di trasferimento complessiva nel caso di N pari e dispari risulta

${\begin{cases}{\hat {X}}(\omega )=\left[|H(\omega )|^{2}+|H(\omega +\pi )|^{2}\right]e^{j\omega (N-1)}X(\omega );\ N\ pari\\{\hat {X}}(\omega )=\left[|H(\omega )|^{2}-|H(\omega +\pi )|^{2}\right]e^{j\omega (N-1)}X(\omega );\ N\ dispari\end{cases}}$

$A(\omega )=|H(\omega )|^{2}+|H(\omega +\pi )|^{2}=1$

(D.59)

Dato che con N dispari, la funzione di trasferimento si annulla per ω = π/2 in quanto H(π/2) = H(3π/2), si utilizzano solamente filtri con N pari. Fissato l’ordine del filtro, l’unico filtro che rende unitaria la funzione di trasferimento complessiva è il filtro

$|H(\omega )|^{2}=cos^{2}\alpha \omega$

(D.60)

le cui caratteristiche selettive, però, sono povere. Ogni altro filtro introduce delle distorsioni di ampiezza. Per la definizione di differenti funzioni di trasferimento si ricorre a procedure iterative per la contemporanea minimizzazione della distorsione introdotta dal filtro e dell’energia in banda soppressa, cioè dell’errore quadratico

$\epsilon =w\int _{\omega _{s}}^{\pi }|H(\omega )|^{2}d\omega +(1-w)\int _{0}^{\frac {\pi }{2}}[A(\omega )-1]^{2}d\omega ;\ 0\leq w\leq 1$

(D.61)

il cui primo termine rappresenta l’energia in banda soppressa, il secondo la distorsione in banda e w un coefficiente di pesatura.

Appendice E

RICHIAMI SU TRASFORMATE NUMERICHE

---

E.1 TRASFORMATA DI FOURIER DISCRETA

Data una funzione x(t) continua in -∞ < t < ∞, la sua trasformata ed antitrasformata di Fourier è data da

$X(\omega )=F[x(t)]\equiv {\sqrt {\frac {1}{2\pi }}}\int _{-\infty }^{\infty }x(t)\ e^{-j\omega t}dt$

$x(t)=F_{-1}[X(\omega )]\equiv {\sqrt {\frac {1}{2\pi }}}\int _{-\infty }^{\infty }X(\omega )\ e^{j\omega t}d\omega$

(E.1)

La trasformata di Fourier può essere estesa dall’analisi dei segnali continui ai segnali a tempo-discreto (Discrete Time Fourier Transform: DTFT) tramite la sostituzione della variabile continua t con l’indice n, ottenendo

$X(e^{j\omega })=\sum _{n=-\infty }^{\infty }x(n)\ e^{-j\omega n}$

(E.2)

La trasformata di Fourier rappresenta la sequenza x(n) come sovrapposizione di esponenziali complessi e^jωn di ampiezza

${\frac {X(e^{j\omega })d\omega }{2\pi }}$

(E.3)

ed è quindi una funzione complessa continua. Inoltre, data la periodicità dell’esponenziale complesso, la DTFT risulta essere una funzione periodica in frequenza. Ciò si ripercuote nell’antitrasformata, nella quale gli estremi di integrazione dell’integrale di definizione sono limitati in un intervallo di ampiezza pari a 2p

$x(n)={\frac {1}{2\pi }}\int _{-\pi }^{\pi }X(e^{j\omega })\ e^{j\omega n}d\omega$

(E.4)

Es.: si consideri la DTFT di una rect(n) discreta

$x(n)={\begin{cases}1&0\leq n\leq N-1\\0&altrove\end{cases}}$

(E.5)

la sua DTFT è data da

$X\left(e^{j\omega }\right)=\sum _{n=0}^{N-1}e^{-j\omega }={\frac {1-e^{-jN\omega n}}{1-e^{-j\omega }}}={\frac {sin(N\omega /2)}{sin(\omega /2)}}e^{-j(N-1){\frac {\omega }{2}}}$

(E.6)

II modulo della DTFT è dato dal solo rapporto tra sinusoidi, che è una funzione periodica, il cui primo zero è in ω = 2π/N. L’esponenziale è relativo alla sola traslazione della rect(n), che non risulta centrata nell'origine. È interessante confrontare questa DTFT con la trasformata di Fourier di un’analoga rect(t) continua, data da

$X\left(e^{j\omega }\right)=\int _{-{\frac {T}{2}}}^{\frac {T}{2}}A\ e^{-j\omega t}dt=A{\frac {sin(T\omega /2)}{\omega /2}}$

(E.7)

cioè una sinc(ω) il cui primo zero è in ω = 2π/T. Diagrammando le due funzioni in modo che gli zeri coincidano (che simula l’ottenimento della rect(n) dalla rect(t) campionata con N punti) si nota l’effetto dell’aliasing che fa discostare la DTFT dalla FT tanto più quanto ci si allontana dall'origine (fig. E.1).

Come la DTFT è l’analoga della trasformata, data una sequenza di N elementi periodica con periodo T, o una sequenza di durata finita periodicizzata, è possibile introdurre nel caso discreto l’equivalente della serie di Fourier (Discrete Fourier Series: DFS) ripetendo la sostituzione della variabile continua t con l’indice n. A causa della periodicità delle sequenze sinusoidali, Fig. E.1 - FT e DTFT della rect.

l’analisi in frequenza ha senso solamente nell'intervallo limitato [0 ≤ ω < 2π]. Per sequenze con pulsazione ω₀ rappresentate da N punti, inoltre, si possono considerare solamente le prime N armoniche, per cui gli indici delle sommatorie nelle definizioni della serie vengono limitati ad N.

$x(mT_{0})=\sum _{k=0}^{N-1}X(k)\ e^{jkm\left({\frac {2\pi }{N}}\right)}=\sum _{k=0}^{N-1}X(k)\ W_{N}^{-km}$

$X(k)={\frac {1}{N}}\sum _{m=0}^{N-1}x(mT_{0})\ e^{jkm\left({\frac {2\pi }{N}}\right)}={\frac {1}{N}}\sum _{m=0}^{N-1}x(mT_{0})\ W_{N}^{km}$

(E.8)

dove T₀=T/N e W^N=e^-j(2π/N) è radice n-esima dell’unità interpretabile come campionamento di funzioni trigonometriche. La limitazione in frequenza comporta un effetto di aliasing che, come per la DTFT, fa si che i coefficienti della serie siano differenti dagli analoghi dei sistemi continui. Il calcolo numerico della serie di Fourier è ricondotto alla soluzione del sistema di N equazioni espresse dalla definizione, una per ogni valore di m.

Esistendo una relazione diretta tra coefficienti della serie e valori dello spettro, le equazioni precedenti stabiliscono una corrispondenza tra N campioni di una sequenza x(n) e N campioni di una trasformata X(k), per cui la serie di Fourier per segnali tempo discreto è anche interpretabile come trasformata discreta di Fourier (DFT), definita come:

$x(n)={\frac {1}{N}}\sum _{k=0}^{N-1}X(k)\ W_{N}^{-kn}$

$X(k)=\sum _{k=0}^{N-1}x(n)\ W_{N}^{kn}$

(E.9)

Es.: data la serie x(n) = {0,2, 1, 1}, la sua DFT è data da

$W_{4}^{kn}=e^{-j{\frac {2\pi }{4}}kn}=cos\left({\frac {2\pi }{4}}kn\right)-j\ sin\left({\frac {2\pi }{4}}kn\right)$

$W_{4}^{0}=cos\left(0\right)-j\ sin\left(0\right)=1$

$W_{4}^{1}=cos\left({\frac {\pi }{2}}\right)-j\ sin\left({\frac {\pi }{2}}\right)=-j$

$W_{4}^{2}=cos\left(\pi \right)-j\ sin\left(\pi \right)=-1$

$W_{4}^{3}=cos\left({\frac {2}{2}}\pi \right)-j\ sin\left({\frac {3}{2}}\pi \right)=+j$

$X(0)=x(0)W_{4}^{0}+x(1)W_{4}^{0}+x(2)W_{4}^{0}+x(3)W_{4}^{0}=4$

$X(1)=x(0)W_{4}^{0}+x(1)W_{4}^{1}+x(2)W_{4}^{2}+x(3)W_{4}^{3}=1-j$

$X(2)=x(0)W_{4}^{0}+x(1)W_{4}^{2}+x(2)W_{4}^{0}+x(3)W_{4}^{2}=2$

$X(3)=x(0)W_{4}^{0}+x(1)W_{4}^{3}+x(2)W_{4}^{2}+x(3)W_{4}^{1}=1+j$

(E.10)

È possibile ora analizzare il legame esistente tra lo spettro di una funzione continua e la DFT di una sequenza finita, ottenuta valutando la Fig. E.2 - Legame tra DTFT e DFT su di una sequenza di durata finita. funzione stessa a multipli di un intervallo T (interpretabile come periodo di campionamento) (fig. E.2). Rappresentare una funzione tramite una serie (infinita) corrisponde a campionare la funzione stessa senza un preventivo filtraggio anti aliasing. Lo spettro della sequenza ottenuta dalla DTFT è quindi la ripetizione periodica dello spettro della funzione (eventualmente infinito) nell'intorno della frequenza F = 1/T, con un suo scalamento per la costante 1/T. A tale periodicizzazione è associato un fenomeno di aliasing. Passando alla DFT, considerare una serie finita di N elementi corrisponde a moltiplicare la serie precedentemente ottenuta dal campionamento con una rect(t) di ampiezza L = N*T. Dal punto di vista della rappresentazione in frequenza, tale troncamento corrisponde ad una convoluzione della DTFT con una sinc(f) (il cui primo zero è posto ad una frequenza pari a 1/L), che provoca una distorsione dello spettro tanto maggiore quanto più breve è la serie. La DFT di una serie di N elementi produce N campioni in frequenza ad intervalli F/N = 1/L. Tale campionamento in frequenza si riflette nella periodicizzazione nel tempo ad intervalli pari ad L. Fissato T, è possibile incrementare la definizione della rappresentazione in frequenza incrementando N tramite l’accodamento alla serie di campioni nulli (zero padding). Riducendo T, invece, si allontanano le repliche dello spettro, permettendo l’osservazione, di un maggiore range di frequenze (fig. E.3).

La diffusione della DFT è legata all'esistenza di una famiglia di algoritmi (Fast Fourier Transform: FFT) in grado di ridurre la complessità computazionale di una DFT da N 2 a N log 2 N: infatti, definita la DFT di un segnale come

$F(k)=\sum _{n=0}^{N-l}f(n)\ W_{N}^{nk}$

(E.11)

dove N è il numero di campioni, W N = eri ^ la radice N-esima dell’unità, f(n) l’n-esimo campione del segnale e F(k) il k-esimo campione dello spettro del segnale, l’applicazione diretta della (E.11) comporta l’effettuazione di N prodotti tra i campioni e le potenze di W N per ciascuno degli N campioni dello spettro, da cui la complessità di N 2 ricordata.

La FFT fornisce N campioni dello spettro del segnale in corrispondenza di frequenze, sia positive che negative, esprimibili come sottomultiple della frequenza di campionamento. Per segnali periodici, ciò vuol dire che campioFig. E.3 - DFT direct con numero di campioni e padding variabile.

nando N volte il segnale all'interno di un periodo, si avranno informazioni sulle prime N/2 armoniche.

Una prima classificazione di tali algoritmi può essere fatta in base alle radici dell’algoritmo stesso: si definiscono radici i fattori primi per mezzo dei quali è possibile scomporre il numero degli elementi della serie da trasformare. Ad esempio, gli algoritmi con radice due si applicano a blocchi di campioni in numero pari ad una potenza di due.

Una seconda classificazione degli algoritmi di FFT è tra algoritmi con decimazione nel tempo o in frequenza: tale distinzione si basa sul modo con il quale vengono organizzati i dati da elaborare. Con la decimazione nel tempo i campioni temporali del segnale da sottoporre a FFT vengono riordinati in modalità bit reversed nel dominio del tempo, mentre con la decimazione in frequenza la sequenza dei campioni non viene alterata, ma vanno riordinati i risultati della FFT. L’idea base di tutti gli algoritmi di FFT è quella di scomporre una sequenza di campioni da trasformare in serie più brevi, le cui DFT, ricombinate, forniscano la DFT del segnale originario. In una FFT con radice pari a due, le serie da trasformare vengono ripetutamente dimezzate fino ad ottenere serie di due soli elementi, dopo log₂N iterazioni. La DFT di una serie di due campioni risulta banalmente calcolabile, in quanto risulta essere dato da

${\begin{cases}F(0)=f(0)W_{2}^{0*0}+f(l)W_{2}^{1*0}=f(0)+f(l)\\F(1)=f(0)W_{2}^{0*1}+f(l)W_{2}^{1*1}=f(0)-f(l)\end{cases}}$

(E.12)

Si fa notare come per il calcolo di tale DFT non venga richiesto il calcolo di moltiplicazioni. Negli algoritmi di decimazione nel tempo la serie dei campioni è ripetutamente scomposta in due sottoserie, la prima delle quali contenente gli elementi di ordine pari della serie originaria, mentre la seconda quella di ordine dispari. L’ordinamento bit reserved deriva da tale progressiva scomposizione. È opportuno notare come le sottoserie generate possano essere viste come risultato di due campionamenti del segnale originario con una frequenza dimezzata ed istanti iniziali opportunamente sfasati. In tal modo la (E.11) può essere riscritta come

$F(k)=\sum _{n=0}^{N/2-1}f(2n)W_{N}^{2nk}+\sum _{n=0}^{N/2-1}f(2n+1)W_{N}^{(2n+1)k}$

$F(k)=\sum _{n=0}^{N/2-1}f_{1}(n)W_{N/2}^{nk}+W_{N}^{k}\sum _{n=0}^{N/2-1}f_{2}(n)W_{N}^{2nk}$

(E.13)

dove f₁ è la serie degli elementi pari, mentre f₂ è la serie degli elementi dispari. Il coefficiente W]^ tiene, appunto, conto dello sfasamento tra le due sottosequenze.

Dato che W_N²=W_N/2¹ (fig. E.4) la (E. 13) può essere scritta come

$F(k)=\sum _{n=0}^{N/2-1}f_{1}(n)W_{N/2}^{nk}+W_{N}^{k}\sum _{n=0}^{N/2-1}f_{2}(n)W_{N/2}^{nk}$

(E.14)

Fig. E.4 - Proprietà della radice n-esima dell’unità.

Applicando la (E.11), la (E. 14) diventa

$F(k)=F_{1}(k)+W_{N}^{k}\ F_{2}(k)$

(E.15)

cioè, la DFT di una serie può essere ottenuta ricombinando secondo la (E. 15) le trasformate di due sue opportune sottoserie ciò grazie alle proprietà di linearità e di traslazione nel tempo della trasformata di Fourier. Dato che la (E. 15) deve valere per 0 < k < N-l, essa va intesa come

$F(k)={\begin{cases}F_{1}(k)+W_{N}^{k}\ F_{2}(k)&0\leq k\leq {\frac {N}{2}}-1\\F_{1}\left(k-{\frac {N}{2}}\right)+W_{N}^{k}\ F_{2}\left(k-{\frac {N}{2}}\right)&{\frac {N}{2}}\leq k\leq N-1\end{cases}}$

(E.16)

Dato che $W_{N}^{k+N/2}=-W_{N}^{k}$, la (E.16) può essere riscritta come

$F(k)={\begin{cases}F_{1}(k)+W_{N}^{k}\ F_{2}(k)&0\leq k\leq {\frac {N}{2}}-1\\F_{1}\left(k-{\frac {N}{2}}\right)+W_{N}^{k-N/2}\ F_{2}\left(k-{\frac {N}{2}}\right)&{\frac {N}{2}}\leq k\leq N-1\end{cases}}$

(E.17)

In tal modo si vede, ad esempio, come la DFT di una serie di 8 campioni può essere ottenuta a partire dalle trasformate di due serie di quattro suoi elementi. Procedendo nella scomposizione di ciascuna sottoserie, si arriverà a poter applicare la (E. 12). Lo spettro completo del segnale sarà ottenibile ricombinando le trasformate tramite ripetute applicazioni della (E. 16).

Visto che la scomposizione del segnale in sottoserie si traduce nel riordinamento bit reversed dei campioni e che la valutazione della (E. 12) non richiede l’effettuazione di nessuna moltiplicazione, l’ordine di complessità N log₂N già ricordato deriva dall'applicazione della (E.16) a ciascuno dei log₂N livelli di scomposizione ottenuti per ognuno degli N campioni dello spettro.

Inoltre, ora si può comprendere che il riordinamento bit reversed ha il compito di riordinare la serie dei campioni in modo da rendere adiacenti i campioni ai quali applicare la (E. 12), che, in realtà, distano N/2 intervalli di campionamento; tali considerazioni valgono anche per le ulteriori sottosequenze in modo tale che la (E.16) possa essere applicata sequenzialmente, ottenendo in uscita la serie ordinata dei campioni dello spettro voluto (fig. E.5).

Come ultima considerazione si pone l’accento sul fatto che le (E.12) coincidono con le (E. 14) scritte per N=2; ciò vuol dire che, ai fini del calcolo, non si distingue tra calcolo delle trasformate ricombinazione delle sottoserie. Il programma di calcolo della FFT si ridurrà all’applicazione iterativa della (E.16) ai campioni opportunamente ordinati, con N crescente dal valore 2 ad N/2 secondo le potenze di 2. Es.: data la serie x(n) = {0, 2,1,1}, la sua FFT è data da

$W_{4}^{0}=1;\ W_{4}^{1}=-j$

$x(00_{b})\to x_{0}+x_{2}=1\to X(0)=4$

$x(10_{b})\to x_{0}+x_{2}=-1\to X(1)=-1-j$

$x(01_{b})\to x_{1}+x_{3}=3\to X(2)=2$

$x(11_{b})\to x_{1}+x_{3}=1\to X(3)=1+j$

(E.18)

E.2 DCT COME APPROSSIMAZIONE DELLA KLT

Si consideri un processo stazionario di Markov del primo ordine, per il quale la matrice di auto-covarianza è

$[R]_{jk}=\rho ^{|1-k|}\ i,k=0,1,....N-1$

(E.19)

per 0 ≤ ρ ≤ 1, dove p ρ il coefficiente di correlazione. La soluzione per la matrice Φ degli autovettori nei caso discreto è

$[\Phi ]_{mn}=\Phi _{m}(n)={\sqrt {\frac {2}{N+\mu _{m}}}}sin\left\{w_{a}\left[(n+1)-{\frac {N+1}{2}}\right]+(m+1){\frac {\pi }{2}}\right\}$

(E.20)

dove $\mu _{m}={\frac {1-\rho ^{2}}{1-2cos(w_{m})+\rho ^{2}}}$ sono gli autovalori e w_m sono le radici reali m positive dell’equazione:

$tan(Nw)={\frac {(1-\rho ^{2})sin(w)}{cos(w_{m})-2\rho +\rho ^{2}cos(w)}}$

(E.21)

Considerando il coefficiente di correlazione ρ → 1, si ha tan(Nw) = 0, da cui deriva

$w_{k}={\frac {k\pi }{N}},\ k=0,1,\cdots ,N-1$

(E.22)

Fig. E.5 - Flusso dei dati in una FFT. Gli autovalori μ_m sono uguali a zero quando i w_m sono diversi da zero. Per μ_m quando m = 0, l’espressione tende all’infinito. In conclusione, si ottiene che

$\sum _{m=0}^{N-1}[A]_{mn}=\sum _{m=0}^{N-1}\mu _{m}$

(E.23)

e, poiché gli elementi diagonali di [A] in un processo di Markov del primo ordine sono tutti unitari, abbiamo immediatamente che Po = N. Sostituendo si ottiene

${\begin{cases}[\varnothing ]_{n\ 0}={\frac {1}{\sqrt {N}}}\\\left[\varnothing \right]_{n\ m}={\sqrt {\frac {2}{N}}}sin\left[m\left(n+{\frac {1}{2}}\right){\frac {\pi }{N}}+{\frac {\pi }{2}}\right]\end{cases}}$

(E.24)

La seconda equazione può essere riscritta come

$[\varnothing ]_{n\ m}={\sqrt {\frac {2}{N}}}cos\left[m\left(n+{\frac {1}{2}}\right){\frac {\pi }{N}}\right],\ m\neq 0$

$={\sqrt {\frac {2}{N}}}k_{m}cos\left[m\left(n+{\frac {1}{2}}\right){\frac {\pi }{N}}\right],\ m,n=0,1,\ldots .N-1$

$k_{m}{\begin{cases}{\frac {1}{\sqrt {2}}}&m=0\\1&altrimenti\end{cases}}$

(E.25)

Questa equazione è la definizione della DCT. Quindi la DCT è asintoticamente equivalente alla KLT per processi Markoviani del 1° ordine quando ρ — > 1. Si noti che ciò è indipendente da N e che l’unico autovalore diverso da zero è μ₀, pari ad N. In termini di MSE, l’errore è compattato in un singolo autovalore.

La ragione è che R_-1 è una matrice simmetrica tri-diagonale, che per uno scalare β₂ = (1 - ρ₂ ) / (1 + ρ₂) ed α = ρ/(l + ρ²) soddisfa la relazione

$\beta ^{2}[A]{-1}={\begin{bmatrix}1-\rho \alpha &-\alpha &0&\ldots \\-\alpha &1&\ldots &\ldots \\\ldots &\ldots &1&-\alpha \\0&\ldots &-\alpha &1-\rho \alpha \end{bmatrix}}$

(E.26)

Si ricava l'approssimazione

$\rho \approx l\to \beta ^{2}R^{1}\approx Q_{y}$

(E.27)

Quindi gli autovettori di [A] e gli autovettori di Q_y, cioè, la trasformata coseno, saranno molto simili.