Interpolazione trigonometrica e FFT

20 gennaio 2009

1 Approssimazione con funzioni trigonometriche

Consideriamo una funzione periodica ^{2}([0 2 ]), : � �, e cerchiamo la miglior approssimazione del tipo

( ) _{}( ) : _{}^{}_{}exp(  )  (1)

dove come al solito è la costante immaginaria. Osserviamo fin d’ora che in alcuni manuali, si studia il problema di cercare un’approssimante del tipo

( ) : _{}^{}_{}exp( )

il che corrisponde a un ordine diverso dei coefficienti _{}. I coefficienti _{}sono detti di Fourier. Di seguito poniamo

( )_{2} _{0}^{2} ( )\overline{( )}  (2)

dove \overline{( )} è il coniugio di ( ). L’insieme _{}^{2}([0 2 ]) delle funzioni di ^{2}([0 2 ]) periodiche una volta dotato del prodotto interno ( )_{2} è uno spazio euclideo, e

span(exp(  )  exp(  ))

è in effetti un suo sottospazio vettoriale di dimensione finita, in cui

exp(  )  exp(  )

è in particolare una sua base ortogonale. Vediamo perchè .

Che la funzione _{}sia periodica è conseguenza dell’uguaglianza di Eulero

exp( ) cos( ) sin( ) �  (3)

Infatti

exp( ( 2 ))  cos( 2 ) sin( 2 )
cos( ) sin( )  (4)

e la somma finita di funzioni periodiche è periodica. Quindi in effetti _{}è un sottoinsieme di _{}^{2}([0 2 ]). Senza molte difficoltà si vede che in effetti è un sottospazio vettoriale di _{}^{2}([0 2 ]). Passiamo a vedere la proprietà di ortogonalità tra elementi di _{}del tipo exp( ) con . Osserviamo che

(exp( )exp( ))_{2} _{0}^{2} exp( ( ) )

e che se

_{0}^{2} exp( ( ) ) _{0}^{2} 1 2

altrimenti,

_{0}^{2} exp( ( ) )  _{0}^{2} cos(( ) )  _{0}^{2} sin(( ) )
\frac{1}{( )} _{0}^{2 ( )} cos( ) \frac{}{( )} _{0}^{2 ( )} ( )
0  (5)

Poniamo _{1}( ) : exp( ) _{2 1} ( ) : exp( ), 2 1. Il teorema di miglior approssimazione in spazi euclidei, stabilisce che l’elemento di miglior approssimazione di _{}^{2}([0 2 ]) è _{}^{}_{1}^{}_{}_{}˜ ( ) dove

˜\frac{( _{})}{( _{}_{})} 1

Di conseguenza essendo

( _{}_{})_{2} 2 1

si ha

˜  1 \frac{( _{1} )}{( _{1}  _{1} )} \frac{1}{2} ( exp(  ))_{2}  (6)

Ora integriamo numericamente

( exp(  ))_{2} _{0}^{2} ( ) exp(  )

mediante la formula dei trapezi

_{0}^{2} ( )  _{}( ) : \frac{}{2} ( _{0}) _{1}^{1} ( _{}) \frac{}{2} ( _{})  \frac{2}{}_{}

Dalla periodicità di , essendo _{0} 0, _{}2 , abbiamo che

( ) : _{0}^{1} ( _{}) _{1}^{}( _{})

e quindi per

\frac{1}{2} ( exp( ))_{2} \frac{1}{2} _{}( exp(  ))
\frac{1}{2} \frac{2}{}_{0}^{1} ( _{})exp(  _{})

\frac{1}{}_{0}^{1} ( _{})exp  \frac{2}{}

\frac{1}{}_{1}^{}( _{1})exp  ( 1) \frac{2}{}(7)

In generale tutti i coefficienti _{}sono calcolati tramite un famoso algoritmo, la FFT (cioè la Fast Fourier Transform) (cf. [1, p.181], [4, p.277], [8, p.398]).

2 Interpolazione trigonometrica

Sia una funzione continua e periodica nell’intervallo [0 2 ] e si pongano

\frac{2}{2 1}  0  2

Il problema dell’interpolazione trigonometrica consiste nel calcolare i coefficienti ˜ , tali che

_{}^{}_{}˜ exp( _{}) ( _{})  0  2  (8)

Osserviamo che

_{0}^{2} exp( ( ) _{}) \begin{matrix}0 & \\ 2 1 & \end{matrix}

Il caso in cui  è ovvio e quindi poniamo attenzione al caso in cui  . Si vede

subito che posto exp( ( ) \frac{2}{2 1} ), essendo exp( 2 ) 1 per un intero ,

_{0}^{2} exp( ( ) _{})  _{0}^{2} exp( ( ) \frac{2}{2 1} )

_{0}^{2} ^{}( ^{2 1} 1)/( 1)

exp ( ) \frac{2}{2 1} ^{2 1} 1 /( 1)
exp( ( ) 2 ) 1 /( 1)
0/( 1) 0  (9)

Moltiplicando ambo i membri di (8) per exp( _{})e sommando in , otteniamo

_{0}^{2} ( _{})exp( _{})  _{0}^{2} exp( _{})_{}^{}_{}˜ exp( _{})

_{}^{}_{}˜_{0}^{2} exp( ( ) _{})
(2 1) _{} ˜  (10)

da cui

˜\frac{1}{2 1} _{0}^{2} ( _{})exp( _{})

Una prima lettura non banale consiste nel fatto che i coefficienti _{}sono noti esplicitamente, senza risolvere numericamente il sistema lineare del tipo dove

exp( _{}) _{}( _{})

Osserviamo ora che in (7), per avevamo approssimato numericamente con la formula dei trapezi fornendo

\frac{1}{}_{1}^{}( _{1})exp  ( 1) \frac{2}{}(11)

Nel caso speciale in cui 2 1, posto 1, abbiamo essendo _{}\frac{2}{2 1}

\frac{1}{2 1} _{0}^{2} ( _{})exp  \frac{2}{2 1} _{} ˜  (12)

cioè l’interpolante polinomiale ha quali coefficienti ˜ _{}^{}quelli ottenuti calcolando numericamente i coefficienti di Fourier _{}con una formula dei trapezi avente 2 1 punti.

3 FFT

Il problema consiste nel calcolare i coefficienti di Fourier _{}_{0 1} per una funzione _{}^{2}([0 2 ]) i cui valori sono noti nei punti \frac{2}{}con 0 1 1. Per quanto visto,

\frac{1}{}_{0}^{1}  \frac{2}{}exp \frac{2}{}0  1  (13)

Si ponga

exp \frac{2}{}_{}\frac{1}{}\frac{2}{}

Possiamo quindi riscrivere il problema come segue: per 0 1, calcolare

_{0}^{1} _{}^{}dove ^{}1  (14)

Se per operazione intendiamo una moltiplicazione e una somma in campo complesso, allora usando lo schema di Horner, si può risolvere questo problema in ^{2} operazioni. L’algoritmo noto come FFT è diventato popolare dopo uno studio di Cooley and Tukey del 1966 (cf. [2]), ma pare fosse già noto in forma limitata a Gauss (1805), permette di risolvere il problema con una complessità inferiore. Consideriamo il caso 2^{}. Sia 2 _{1} quando è pari, e 2 _{1} 1 quando è dispari. Da (14) abbiamo

_{_{1} 0}^{\frac{}{2} 1} _{2 _{1}} ( ^{2})^{_{1}} _{_{1} 0}^{\frac{}{2} 1} _{2 _{1} 1} ( ^{2})^{_{1}} ^{}

Sia \frac{}{2} _{1}. Allora, poichè ^{}1 (e quindi ( ^{})^{_{1}} 1)

( ^{2})^{_{1}} ( ^{2})^{\frac{_{1}}{2}} ( ^{2})^{_{1} _{1}} ( ^{})^{_{1}} ( ^{2})^{_{1} _{1}} ( ^{2})^{_{1} _{1}}

Così se poniamo

( _{1})  _{_{1} 0}^{\frac{}{2} 1} _{2 _{1}} ( ^{2})^{_{1}} ( _{1} 0 1  \frac{}{2} 1)  (15)

( _{1})  _{_{1} 0}^{\frac{}{2} 1} _{2 _{1} 1} ( ^{2})^{_{1}}  (16)

(17)

abbiamo

( _{1}) ^{}( _{1})  0  1

Si noti che il calcolo di ( _{1}) e ( _{1}) significa che vengono compiute due analisi di Fourier con \frac{}{2} 2^{1} termini invece di una con 2^{}. Ora applichiamo la stessa idea alle due analisi di Fourier. Si hanno così 4 analisi di Fourier, ognuna ha 2^{2} termini; queste possono essere ulteriormente divise in 8 analisi di Fourier con 2^{3} termini, eccetera. Il numero di operazioni richieste per calcolare _{}quando ( _{1}) e ( _{1}) sono stati calcolati, è al più 2 2^{}(un numero maggiore di calcoli può essere evitato se le potenze di sono immagazzinate in memoria). Sia _{}il numero totale di operazioni necessarie per calcolare i coefficienti quando 2^{}. Come prima, abbiamo

2 _{1} 2 2^{}1 2

Poichè _{0} 0, si ha per induzione, essendo 2^{}e quindi log_{2}( )

2 _{1} 2 2^{}2 2 _{2} 2 2^{1} 2 2^{}
4 _{2} 2 2 2^{}4 2 _{3} 2 2^{2} 4 2^{}
2^{3} _{3} 8 2^{2} 4 2^{}2^{3} _{3} 2 2^{}2^{2} 2^{}
2^{3} _{3} 2 3 2^{}2^{}_{0} 2 2^{}2 log_{2}  (18)

da paragonarsi con ^{2} operazioni del metodo di base.

Per capire i vantaggi, facciamo qualche conto con Matlab, mostrando per alcune potenze di 2 (prima colonna), il numero di operazioni per il metodo di base (seconda colonna) e dell’ FFT (terza colonna):

>> N=2.^(1:10);
>> N=N';
>> N2=N.^2;
>> Nlog=2*N.*log2(N);
>> [N N2 Nlog]
ans =
           2           4           4
           4          16          16
           8          64          48
          16         256         128
          32        1024         320
          64        4096         768
         128       16384        1792
         256       65536        4096
         512      262144        9216
        1024     1048576       20480
>>

4 FFT in Matlab

Vediamo come utilizzare in Matlab la FFT partendo dal relativo help

>> help fft
 FFT Discrete Fourier transform.
    FFT(X) is the discrete Fourier transform (DFT) of vector X.  For
    matrices, the FFT operation is applied to each column. For N-D
    arrays, the FFT operation operates on the first non-singleton
    dimension.
    FFT(X,N) is the N-point FFT, padded with zeros if X has less

    than N points and truncated if it has more.
    FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the
    dimension DIM.
    For length N input vector x, the DFT is a length N vector X,
    with elements
                     N
       X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
                    n=1
    The inverse DFT (computed by IFFT) is given by
                     N
       x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
                    k=1
    See also IFFT, FFT2, IFFT2, FFTSHIFT.
 Overloaded methods
    help qfft/fft.m

Si capisce che tale procedura calcola, a partire da un vettore ( ( ))_{1}, il vettore ( ( ))_{1} definito da

( ) _{1}^{}( ) exp \frac{2 ( 1)( 1)}{}1  (19)

(si noti che per errore nell’help di Matlab hanno indicato con la costante immaginaria ). Ci si domanda quale sia il legame di (19) con (13). Posto 1, 1, ( ) \frac{1}{}\frac{2 ( 1)}{}, abbiamo

1 _{1}^{}( )exp \frac{( 1)2 ( 1)}{}1

e quindi ( ) _{1} . Vediamo un esempio in Matlab, in cui analizziamo per 8 la funzione

( ) 5 exp(2 )

in cui evidentemente _{0} _{1} _{3} _{8} 0 e _{2} 5. Quindi ci aspettiamo un vettore _{}in cui poichè _{}_{1} si ha che _{3} 5 e tutte le altre componenti nulle. Salviamo nel file esempio_FFT

clear all;
format long;
n=8;

f=inline('5*exp(i*2*t)');
h=2*pi/n;
t=(0:n-1)'*h;
ft=feval(f,t);
c_fft=fft(ft,n)*h/(2*pi);

e otteniamo come prevedibile

>> c_fft
c_fft =
 -0.00000000000000 + 0.00000000000000i
 -0.00000000000000 + 0.00000000000000i
  5.00000000000000 - 0.00000000000000i
  0.00000000000000 + 0.00000000000000i
  0.00000000000000 + 0.00000000000000i
  0.00000000000000 + 0.00000000000000i
                 0 + 0.00000000000000i
 -0.00000000000000 + 0.00000000000000i

Ora consideriamo ( ) 10 exp( 2 )

e osserviamo che il risultato è

>> esempio_FFT
>> c_fft
c_fft =
 -0.00000000000000 - 0.00000000000000i
 -0.00000000000000 - 0.00000000000000i
                 0 - 0.00000000000000i
  0.00000000000000 - 0.00000000000000i
  0.00000000000000 - 0.00000000000000i
  0.00000000000000 - 0.00000000000000i
 10.00000000000000 + 0.00000000000000i
 -0.00000000000000 - 0.00000000000000i

A prima vista la cosa è sorprendente, visto che il coefficiente uguale a 10 appare alla settima componente, quella che dovrebbe essere di _{6} e non _{2} . Vediamo perchè . Si ha in generale per ( ) exp( )

( )  _{1}^{}( ) exp \frac{2 ( 1)( 1)}{}

\frac{1}{}_{1}^{}exp  \frac{2 ( 1)}{}exp \frac{2 ( 1)( 1)}{}

\frac{1}{}_{1}^{}exp 2 \frac{2 ( 1 )( 1)}{}(20)

sommatoria che vale se 1 è multiplo di (con avente valori tra 1 e N) e 0 altrimenti. Quindi nell’esempio precedente in cui 10, 2, 8, la sommatoria vale 10 se 1 2 1 è multiplo di 8 (cioè 7) e 0 altrimenti.

Per convincersi di questo, posto 8 applichiamo la FFT a

( ) 10 exp( 6 ) 20 exp(2 )

Ci aspettiamo che il primo termine contribuisca alla sommatoria, cosicche’ la terza componente valga 10 e il secondo termine contribuisca alla sommatoria, cosicche’ la terza componente valga 20. Per linearità avremo quindi una terza componente uguale a 30, in quanto in generale

_{1}^{}( ( ) ( )) exp( ) _{1}^{}( ) exp( ) _{1}^{}( ) exp( )

>> esempio_FFT
>> c_fft
c_fft =
 -0.00000000000000 - 0.00000000000000i
  0.00000000000000 - 0.00000000000000i
 30.00000000000000 + 0.00000000000000i
 -0.00000000000000 - 0.00000000000000i
  0.00000000000000 - 0.00000000000000i
 -0.00000000000000 + 0.00000000000000i
                 0 - 0.00000000000000i
  0.00000000000000 + 0.00000000000000i
>>

5 Online

Citiamo alcuni siti web che espongono parte dell’argomento:

http://en.wikipedia.org/wiki/Cooley-Tukey_FFT_algorithm
http://en.wikipedia.org/wiki/Fast_Fourier_transform
http://it.wikipedia.org/wiki/Rappresentazione_spettrale_dei_segnali

References

[1] K. Atkinson, An Introduction to Numerical Analysis, Wiley, (1989).

[2] J.W. Cooley e J.W. Tukey, An algorithm for the machine calculation of complex Fourier series, Math. Comput. 19, p. 297-301, (1965).

[3] G. Dahlquist e A. Bjorck , Numerical methods, Dover, (2003).

[4] S.D. Conte e C. de Boor, Elementary Numerical Analysis, An Algorithmic Approach, McGraw-Hill, (1981).

[5] G. Gilardi Analisi Due, seconda edizione, McGraw-Hill, (1996).

[6] D.H. Griffel, Applied functional analysis, Dover publications, 2002.

[7] A.N. Kolmogorov e S.V. Fomin, Introductory Real Analysis, Dover publications, 1970.

[8] A. Quarteroni, R. Sacco e F. Saleri Matematica Numerica, Springer, (1998).

[9] Wikipedia, (Serie di Fourier), http://it.wikipedia.org/wiki/Serie_di_Fourier .