Polinomi Trigonometrici Interpolazione Trigonometrica DFT & FFT

1. Polinomi TrigonometriciInterpolazione TrigonometricaDFT & FFT Seminario Metodi Matematici per l’Ottimizzazione Anno accademico 2010/2011 Francesco Cambria Fabio Catrini

2. Introduzione Nel campo dell'analisi numerica può risultare utile effettuare l'approssimazione di una funzione. Tale esigenza scaturisce da varie cause, l'assenza dell' espressione analitica della funzione che descriva un certo fenomeno fisico o pur essendo nota l'espressione analitica, essa risulta complicata da gestire. Nel primo caso si deve trovare una funzione approssimante a partire dai punti noti,mentre nel secondo caso occorre sostituire la funzione nota con una più semplice dal punto di vista operativo.Esistono due modi diversi di risolvere il problema dell'approssimazione di una funzione,questi prendono il nome di Interpolazione e di Approssimazione ai minimi quadrati. L'interpolazione è un metodo nel quale si richiede alla funzione approssimante di assumere lo stesso valore della funzione da approssimare.

3. Definizione polinomio trigonometrico Se P(x,y) è un polinomio nelle variabili x e y un polinomio trigonometrico è semplicemente P(cosx,sinx). Il tipico polinomio trigonometrico è la serie di Fourier una funzione periodica che può essere espressa come somma di armoniche di differenti frequenze e ampiezze : (1.1) I numeri complessi sono i coefficienti del polinomio trigonometrico e dipendono da f.

4. Funzioni Armoniche La funzione armonica ha forma: Dove: A -> Ampiezza w -> Pulsazione(2π/T) Φ -> Fase a -> AsinΦ b ->AcosΦ T -> periodo Frequenza -> 1/T numero di periodi compiuti in un secondo misurata in hertz (Hz)

5. Composizioni di armoniche

6. Calcolo dei coefficienti 1 Date le seguenti identità: (1.2) (1.3) (1.4) Supponiamo che valga la (1.1); dalla (1.1), moltiplichiamo ambo i membri per cos(lx), si ha:

7. Calcolo dei coefficienti 1.1 Integrando ambo i membri ed utilizzando la proprietà distribuita, si ha: (*) Distinguiamo i seguenti casi: k=l utilizzando la (1.2) e la (1,4), si ha: da cui

8. Calcolo dei coefficienti 1.2 2) k≠l questo non dà informazioni. Per calcolare consideriamo l’equazione: Riducendo abbiamo:

9. Calcolo dei coefficienti 1.3 Analogamente, dalla (1.1), moltiplichiamo ambo i membri per sen(lx), si ha: Integrando ambi i membri ed utilizzando la proprietà distributiva, si ha: (**)

10. Calcolo dei coefficienti 1.4 Distinguiamo i seguenti casi: k=l utilizzando la (1.3) e la (1.4), si ha: da cui: 2) k≠l questo non dà informazioni

11. Calcolo dei coefficienti 1.5 Quindi i coefficienti di Fourier in forma trigonometrica sono:

12. Polinomi trigonometrici di Fourier in forma complessa 2 Consideriamo l’intervallo (0,2π). Ricordiamo l’identità di Eulero: da cui si ottiene: I polinomi trigonometrici di Fourier in forma complessa sono dati da: dove i indica l’unità immaginaria. Si tratta pertanto di funzioni a valori complessi periodiche di periodo 2π

13. Polinomi trigonometrici di Fourier in forma complessa 2.1 La notazione per indicare l’insieme delle funzioni: Con prodotto scalare e norma definiti rispettivamente da: (f,g)= Se , la sua serie di Fourier è definita da :

14. Polinomi trigonometrici di Fourier in forma complessa 2.2 dove α(x) è la parte reale di f e β(x) è quella immaginaria.

15. Polinomi trigonometrici di Fourier in forma complessa 2.3 I coefficienti di Fourier di f possono essere scritti come : Infatti:

16. Interpolazione Trigonometrica 3 In matematica, interpolazione trigonometrica è interpolazione con polinomi trigonometrici, cioè una somma di seni e coseni di dati periodi. . L'interpolazione è il processo di individuazione della funzione che passa attraverso alcuni dati punti di riferimenti. Questa forma è adatta particolarmente per interpolazione di funzioni periodiche.

17. Interpolazione trigonometrica 3.1 Un polinomio trigonometrico che la interpoli negli n+1 nodi ovvero tale che L’interpolatore trigonometrico si ottiene attraverso una combinazione lineare di seni e coseni. Assumerà le forme seguenti: se n pari (1) Dove M=n/2 mentre, se n è dispari, M=(n-1)/2 (2)

18. Interpolazione trigonometrica 3.2 Con n pari Possiamo scrivere la funzione (1) come dove i è l’unità immaginaria. I coefficienti ck sono legati ai coefficienti ak e bk nel modo seguente

19. Interpolazione trigonometrica 3.2 Con n dispari possiamo scrivere la funzione (2) come I coefficienti ck per k=0,….M sono determinati come prima, mentre scriviamo con μ=0 se n è pari, μ=1 se n è dispari.

20. Comando interpft MATLAB Il comando interpft calcola l’interpolatore trigonometrico di un insieme di dati. Richiede come parametri d’ingresso un numero intero N ed un vettore le cui componenti sono i valori d’ingresso assunti da una funzione periodica nei punti xi=ip/M, i=1,……,M-1. Il programma interpft restituisce gli N valori dell’interpolatore trigonometrico, ottenuto con la trasformata di Fourier, nei nodi ti=ip/N, i=0,…..,N-1.

21. Esempio x=pi/5*[0:9];y=x.*(x-2*pi).*exp(-x);z=interpft(y,100) Questo esempio prende una funzione in [0,2π] e valuta in 10 nodi equispaziati e calcola i valori dell’interpolatore trigonometrico in 100 nodi equispaziati. La linea tratteggiata è la nostra funzione ed in corrispondenza della linea continua l’interpolatore trigonometrico relativo a 10 nodi equispaziati.

22. Aliasing L’accuratezza dell’interpolazione trigonometrica può in certe situazione subire un forte degrado chiamato aliasing e si può manifestare ogni volta che in una stessa funzione coesistono componenti con frequenza diversa: finchè il numero di nodi non è sufficientemente alto per risolvere le frequenze più elevate, queste ultime potranno interferire con le frequenze più basse, dando origine ad approssimazioni inaccurate. Solo aumentando il numero di nodi sarà possibile approssimare correttamente le funzioni di frequenza più elevata.

23. Trasformata di Fourier ed FFT

24. Analisi spettrale • L’ analisi spettrale è la rappresentazione delle componenti in frequenza di un segnale. • Lo spettro è il vettore delle ampiezze delle componenti di un segnale, • disposte in funzione della loro frequenza. • Un segnale è in teoria rappresentato da una serie infinita di sinusoidi. • Lo spettro si stima tramite: • Metodo tradizionale (non parametrico): Trasformata di Fourier. • Metodo parametrico: basato su modelli (lineari) del segnale.

25. Perché utilizzare l’analisi in frequenze Ad esempio, in ottica, alcuni colori (rosso, giallo, blu), detti fondamentali, sono puri, cioè non ulteriormente scomponibili. A ciascuno di essi corrisponde una certa lunghezza d'onda (frequenza) del raggio luminoso, e il prisma (che scompone la luce bianca nei sette colori dello spettro luminoso) mostrerà solamente quella componente. La medesima cosa avviene per gli altri segnali.

26. Spettro di un segnale armonico a frequenza equispaziate Se le componenti sono in rapporto di frequenza intero con la componente di frequenza più bassa, si dicono armoniche. La componente a frequenza più bassa si chiama fondamentale o prima armonica e si indica con F0. La componente di frequenza doppia della fondamentale si chiama seconda armonica (y = sin(2x) ), la componente di frequenza tripla (y = sin(3x) ) della fondamentale si chiama terza armonica, e così via.

27. Serie di Fourier Qualunque segnale periodico può essere scomposto nella somma di un eventuale termine costante e di componenti sinusoidali, delle quali la prima si chiama prima armonica o fondamentale, e le altre, aventi periodi sottomultipli e quindi frequenze multiple, si chiamano armoniche superiori. La serie di Fourier è quindi una rappresentazione di una funzione periodica con periodo mediante una somma di funzioni periodiche della forma le quali sono le potenze di , cioè le sue armoniche.

28. Serie di Fourier: criteri di convergenza Affinchè la serie trigonometrica converga effettivamente a f(x) si deve rispettare il Criterio (o Condizioni) di Dirichlet che impone: f(x) deve essere definita nell'intervallo t0-t0+T; sono ammessi anche eventuali punti di discontinuità purchè in numero finito. f(x) e la sua derivata prima f'(x) devono essere continue a tratti nell'intervallo t0-t0+T.Se f(x) soddisfa alle ipotesi del Teorema di Dirichlet parleremo di sviluppo di f(x) in serie di Fourier anche se la serie di Fourier di f(x) può non coincidere nei punti di discontinuità con f(x).

29. Esempio Onda quadra di periodo 1 con ampiezza 4 approssimata con k=5 e k = 15

30. La Trasformata di Fourier Per l’approssimazione di segnali aperiodici si applica la trasformata di Fourier • La trasformata di Fourier è l’operatore che permette di trasformare il segnale originarioespresso nel dominio del tempo in una “somma” di segnali sinusoidali a tutte le frequenze. • La trasformata discreta di Fourier opera su sequenze di valori discreti, equispaziati in un intervallo finito.

31. DFT ed interpolazione Per il calcolo della DFT faremo uso dell’interpolazione in punti ugualmente spaziati. Più precisamente:  dove e i punti , sono ugualmente spaziati sul cerchio unitario I punti si ripetono quando j aumentaogni volta di 2n+1.

32. Calcolo della DFT: cerchio unitario Il cerchio unitario è il luogo dei punti del piano aventi una distanza minore o uguale all'unità da un punto detto centro del cerchio. In altri termini il cerchio unitario comprende la circonferenza unitaria e la parte di piano racchiusa dalla circonferenza stessa. Esso è indicato da:

33. DFT ed interpolazione Scriviamo il polinomio di interpolazione come: e imponiamo le condizioni di interpolazione: Si ha:

34. DFT ed interpolazione Per trovare i coefficienti , usiamo il seguente lemma. Lemma Per ogni intero k :

35. DFT ed interpolazione Dato il polinomio precedente: Moltiplicando ogni j-esima equazione per restringendo l’intervallo su e sommando su j otterremo: Invertiamo l’ ordine delle somme e usiamo il lemma per ottenere:

36. DFT ed interpolazione Sostituendo in: si avranno i coefficienti ottenuti come segue: Tali coefficienti sono detti trasformate discrete di Fourier (DTF) dei dati e danno una formula esplicita per il polinomio di interpolazione

37. DFT E INTERPOLAZIONE Effettuando un cambio di notazione i coefficienti possono essere scritti come: dove m=2n+1, k=l, è detta radice principale di ordine m dell’ unità, poiché Segue che per Infatti Quindi basta calcolare

38. Richiami di algebra lineare: • Matrice di Vandermonde E’ una matrice le cui righe (oppure le cui colonne) hanno elementi, a partire da 1, in progressione geometrica (oppure la trasposta )

39. LA DFT IN FORMA MATRICIALE Poniamo per comodità. Consideriamo con Possiamo scriverli in forma matriciale ottenendo: Ponendo con e , con La DFT in forma matriciale sarà:

40. Costo della DFT Il calcolo diretto di dk richiede (m-1) moltiplicazioni complesse. Per calcolare l’intera sequenza saranno necessarie m(m-1) operazioni per una complessità finale di O( )

41. FFT La FFT è un algoritmosviluppatoneglianni ’60 per ilcalcolovelocedella DFT, grazie al qualesi ha unariduzionerilevante del numerodioperazionidacompierenelcalcolodella DFT. La FFT cipermettediridurreilnumerodioperazionidacompiere (in particolareilnumerodimoltiplicazioni). Questo è un aspettomoltoimportante, soprattuttoquandol’implementazioneavvienesu un sistemachedevegarantire un processing real-time.

42. FFT: approccio “divide et impera” Il cuore dell’algoritmo FFT sta nella scomposizione in due sommatorie della sommatoria per il calcolo dei coefficienti della DFT selezionando in modo “opportuno” gli indici. L’algoritmo più famoso è quello di Cooley-Tukey.

43. Cooley-Tukey L'algoritmo FFT piùdiffuso è l'algoritmodi Cooley-Tukey. Si basasul principio didivide et impera: divide ricorsivamenteuna DFT diqualsiasidimensione N con N numerocomposto tale che N=N1N2 in DFT piùpiccoledidimensioni N1 e N2, insieme a O(n) moltiplicazioni per l'unitàimmaginaria, dettifattori twiddle.

44. Cooley-Tukey: Radix-2 L’algoritmodi cui parleremo è il Radix-2. Caratteristichedell’algoritmo: Sfrutta sia la simmetria che la periodicità dei fattori twiddle Considera il caso particolare di N campioni di dati, con potenza intera di 2 Si separa l’input x[n] in due sequenze di lunghezza N/2 I campioni di indice pari nella prima sequenza I campioni di indice dispari nella seconda sequenza moltiplicato per il così detto fattore twiddle

45. Cooley-Tukey: Radix-2 SiaX[k]ilcoefficientedella DFT dacalcolare. Suddiviamo la sommatoria per ilcalcolodella DFT in due, sommandosugliindicipari e sugliindicidispari: Operiamo le sostituzioni n=2r per n pari e n=2r+1 per i dispari ponendo

46. Cooley-Tukey: Radix-2 Un segnale con N punti viene decomposto in N segnali ognuno contenente un solo segnale. Ad ogni passo il segnale viene suddiviso tra le sue componenti pari e dispari.

47. Cooley-Tukey: Radix-2 • Le due sommatorie rappresentano rispettivamente la DFT G[k] su N/2 punti della sequenza x[2n] cioè della sequenza di campioni di posto pari e la DFT H[k] su N/2 punti della sequenza x[2r+1] cioè della sequenza dei campioni di posto dispari: • G[k] e H[k] sono sequenze periodiche di periodo N/2, cioè : • G[k + N/2] = G[k] • H[k + N/2] = H[k] • Inoltre vale che:

48. Radix-2 : decimazione del tempo • Avremo dunque che il calcolo della DFT X[k] su N punti a partire da quello delle due DFT su N/2 punti G[k] e H[k] è dato da: • Il grafo corrispondente a tale coppia di equazioni è: Butterfly

49. Decimazione del tempo per 8 punti • Due DFT a N/2 punti: Osserviamo che G[k], così come H[k] richiede moltiplicazioni complesse ; il numero di moltiplicazioni richieste per è invece N/2. Complessivamente per calcolare X[k] occorrono Quindi più efficiente della DFT diretta di un fattore 2 circa (per N grande)

50. Decimazione del tempo: ricorsione La proceduradidecimazionenel tempo risultadunquevantaggiosa e convieneriapplicarla al calcolodelle due DFT su N/2 punti G[k] e H[k] dellasottosequenzaG[n] deicampionidipostopari e H[n] dicampionidipostodispari.