Interpolazione trigonometrica in forma esponenziale

[dissonance]

Ho una funzione \( \displaystyle f\colon \mathbb{R}\to \mathbb{R} \) periodica di periodo finito \( \displaystyle T \) . Ne prendo \( \displaystyle N \) campioni

\( \displaystyle y_k=f(\frac{T}{N}k) \) e ne faccio la DFT \( \displaystyle Y_j=\frac{1}{N}\sum_{k=0}^{N-1}y_ke^{\imath(\frac{2\pi}{T}j)k} \) ,

ottenendo il vettore \( \displaystyle \bold{Y} \) di \( \displaystyle \mathbb{C}^N \) . Ora come posso, tramite questo vettore, costruire il polinomio trigonometrico di interpolazione della \( \displaystyle f \) nei nodi \( \displaystyle x_k=k\dfrac{T}{N} \) in forma esponenziale (intendo una rappresentazione tipo \( \displaystyle P(t)=\sum_{k=-a}^bc_ke^{\imath\frac{2\pi}{T}bt} \) )? Sul libro che sto leggendo si costruisce solo la rappresentazione con seni e coseni nel caso particolare \( \displaystyle T=2\pi \) e mi sto incasinando con i conti.

[dissonance]

Ok, trovato. Che fatica, però! Riporto qui il risultato nel caso servisse a qualcuno. Per la definizione di \( \displaystyle f, \bold{y}, \bold{Y}, N, T \) vedere il post precedente.

L'idea è quella di definire un polinomio trigonometrico assegnando ad ogni entrata del vettore \( \displaystyle \bold{Y} \) una frequenza, secondo questo schema:

1) \( \displaystyle N \) è dispari, \( \displaystyle N=2M-1 \) , allora si suddivide \( \displaystyle \bold{Y} \) come segue:

\( \displaystyle Y_0, \left[ Y_1 \ldots Y_{M-1}\right],\ \left[Y_M \ldots Y_{N-1} \right] \)

assegnando ai coefficienti del primo gruppo le armoniche corrispondenti, e ai coefficienti del secondo gruppo le armoniche negative secondo la corrispondenza: a \( \displaystyle Y_{N-j} \) corrisponde la \( \displaystyle -j \) -esima armonica.

2) \( \displaystyle N \) è pari, \( \displaystyle N=2M \) , allora si suddivide \( \displaystyle \bold{Y} \) come segue:

\( \displaystyle Y_0, \left[ Y_1 \ldots Y_{M-1}\right],\ Y_M,\ \left[Y_{M+1} \ldots Y_{N-1} \right] \)

poi si procede come nel caso 1).



***Formalmente***:



Sia \( \displaystyle \omega_N=e^{\imath \frac{2\pi}{N}} \) . E' un fatto algebrico noto che \( \displaystyle \omega_N \) è radice primitiva N-esima dell'unità; in particolare l'aritmetica delle sue potenze è modulo \( \displaystyle N \) :

(*) \( \displaystyle \omega_N^k=\omega_N^h\quad \iff k \equiv h \mod N \) .

Per costruzione il polinomio \( \displaystyle P(w)=\sum_{j=0}^{N-1}Y_j w^j \) soddisfa le condizioni

\( \displaystyle P(\omega_N^k)=y_k,\quad k=0\ldots N-1 \) ;

da queste ricaviamo il polinomio trigonometrico di interpolazione della funzione \( \displaystyle f \) distinguendo due casi:

1) \( \displaystyle N=2M-1 \) :
in questo caso possiamo riscrivere la condizione \( \displaystyle P(\omega_N^k)=y_k \) come

\( \displaystyle y_k=\sum_{j=0}^{N-1}Y_j\omega_N^{kj}=\sum_{j=0}^{M-1}Y_j\omega_N^{kj}+\sum_{j=M}^{2M-2}Y_j\omega_N^{kj} \) ;

usando, nella seconda sommatoria, il cambiamento di indice \( \displaystyle h=N-j \) , e ricordando la (*), abbiamo allora

\( \displaystyle y_k=\sum_{j=0}^{M-1}Y_j\omega_N^{kj}+\sum_{h=1}^{M-1}Y_{N-h}\omega_N^{k(-h)} \) .

Quindi il polinomio trigonometrico

\( \displaystyle F(t)=\sum_{j=0}^{M-1}Y_je^{\imath(\frac{2\pi}{T}j)t}+\sum_{h=1}^{M-1}Y_{N-h}e^{\imath(\frac{2\pi}{T}(-h))t \)

verifica le identità

\( \displaystyle F(k\frac{T}{N})=y_k=f(k\frac{T}{N}),\quad k=0...N-1 \) .

2) \( \displaystyle N=2M \)
in questo caso considereremo a parte l'addendo corrispondente alla \( \displaystyle M \) -esima frequenza:

\( \displaystyle y_k=\sum_{j=0}^{M-1}Y_j\omega_N^{kj}+Y_M\omega_N^{kM}+\sum_{h=1}^{M-1}Y_{N-h}\omega_N^{k(-h)} \) .

Osserviamo poi che \( \displaystyle \omega_N^{kM}=\Re(\omega_N^{kM})=\cos(\pi k) \) .

Il polinomio trigonometrico cercato è

\( \displaystyle F(t)=\sum_{j=0}^{M-1}Y_je^{\imath(\frac{2\pi}{T}j)t}+Y_M\cos(\frac{2\pi}{T}Mt)+\sum_{h=1}^{M-1}Y_{N-h}e^{\imath(\frac{2\pi}{T}(-h))t \) .

_________________________________
NOTA BENE: Anche se nel corso di questi post ho usato l'articolo determinativo (il polinomio trigonometrico di interpolazione), non è vero che esso è unico. La scelta indicata è però la migliore per vari motivi. Vedi http://en.wikipedia.org/wiki/Discrete_F ... polynomial

[dissonance]

Per concludere aggiungo la traduzione dei polinomi trigonometrici trovati sopra nella rappresentazione reale in seni e coseni. Le ipotesi sono quelle precedenti, in particolare ora diventa necessario che \( \displaystyle f \) assuma valori reali.

Def.: Per ogni \( \displaystyle j=0\ldots N-1 \) , siano \( \displaystyle \frac{\alpha_j-\imath \beta_j}{2}=Y_j \) .
Risulta che

\( \displaystyle \alpha_j= \frac{2}{N}\sum_{k=0}^{N-1}y_k\cos(\frac{2\pi}{N}jk) \) ;

\( \displaystyle \beta_j= \frac{2}{N}\sum_{k=0}^{N-1}y_k\sin(\frac{2\pi}{N}jk) \) .

________________

1) \( \displaystyle N=2M-1 \)

In questo caso è

\( \displaystyle F(t)=Y_0+ \sum_{j=1}^{M-1}Y_je^{\imath (\frac{2\pi}{T}j)t}+\sum_{h=1}^{M-1}Y_{N-h}e^{\imath (\frac{2\pi}{T}(-h))t} \) ;

ma osserviamo che \( \displaystyle Y_{N-h}=\overline{Y_h} \) per ogni \( \displaystyle h=1\ldots M-1 \) perché il vettore \( \displaystyle \bold{y} \) è reale; dunque

\( \displaystyle F(t)=Y_0+\sum_{j=1}^{M-1}2\Re(Y_je^{\imath (\frac{2\pi}{T}j)t}) \)

e a conti fatti

\( \displaystyle F(t)=\frac{\alpha_0}{2}+\sum_{j=1}^{M-1}\left( \alpha_j\cos(\frac{2\pi}{T}jt) + \beta_j \sin(\frac{2\pi}{T}jt) \right) \) .

2) \( \displaystyle N=2M \)

In questo caso è

\( \displaystyle F(t)=Y_0+ \sum_{j=1}^{M-1}Y_je^{\imath (\frac{2\pi}{T}j)t}+Y_M\cos(\frac{2\pi}{T}Mt)+\sum_{h=1}^{M-1}Y_{N-h}e^{\imath (\frac{2\pi}{T}(-h))t} \) ;

si arriva, con ragionamento analogo a quello del caso 1), al risultato

\( \displaystyle F(t)=\frac{\alpha_0}{2}+\sum_{j=1}^{M-1}\left( \cos(\frac{2\pi}{T}j t)+\beta_j \sin(\frac{2\pi}{T}j t) \right) + \frac{\alpha_M}{2}\cos(\frac{2\pi}{T}M t) \) .

L'unico punto da osservare esplicitamente è l'identità

\( \displaystyle Y_M=\frac{\alpha_M}{2} \) ; questo è immediato se si osserva che, essendo

\( \displaystyle Y_M=\frac{1}{N}\sum_{k=0}^{N-1}y_k e^{-\imath (\frac{2\pi}{N}\frac{N}{2})k} \)

esso è un numero reale.