Seminario Metodi Matematici per lOttimizzazione Anno accademico 2010/2011 Francesco Cambria Fabio...

Polinomi TrigonometriciInterpolazione Trigonometrica

DFT & FFTSeminario Metodi Matematici per

l’Ottimizzazione Anno accademico 2010/2011

Francesco CambriaFabio Catrini

IntroduzioneNel 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.

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.

f (x) 1

2a0 (ak coskx bk sinkx)

k1

a0,ak ,bk

x

Funzioni Armoniche

La funzione armonica ha forma:

Dove:A -> Ampiezzaw -> Pulsazione(2π/T)Φ -> Fasea -> AsinΦb ->AcosΦ

( ) sin( ) cos siny t A wt a wt b wt ( ) cos( ) cos siny t A wt b wt a wt

T -> periodo

Frequenza -> 1/T numero di periodi compiuti in un secondo misurata in hertz (Hz)

Composizioni di armoniche

Calcolo dei coefficienti 1

cos kx cos lx dx 0, se k l

, se k l

0

2

sen kx sen lx dx , se k l

0, altrimenti

0

2

cos kx sen lx dx 00

2

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:

k,l

k,l

k,l

f x cos lx a0

2cos lx ak cos kx bksen kx cos lx

k1

Integrando ambo i membri ed utilizzando la proprietà distribuita, si ha:

(*)

Distinguiamo i seguenti casi:

1) k=l utilizzando la (1.2) e la (1,4), si ha:

da cui

f x cos lx dx a0

2cos lx

0

2

dx ak cos kx cos lx dx0

2

k1

bk sen kx cos lx dx0

2

k1

0

2

0al 0

al 1

f x cos lx dx

0

2

Calcolo dei coefficienti 1.1

2) k≠l questo non dà informazioni.

Per calcolare consideriamo l’equazione:

Riducendo abbiamo:


f x dx a0

20

2

dx ak cos kx dx0

2

k1

bk sen kx dx0

2

k1

0

2

f x dx 2 a0

200

0

2

a0 1

f (x)dx

0

2

a0

Analogamente, dalla (1.1), moltiplichiamo ambo i membri per sen(lx), si ha:

Integrando ambi i membri ed utilizzando la proprietà distributiva, si ha:

(**)

f x sen lx a0

2sen lx ak cos kx bksen kx sen lx

k1

f x sen lx dx a0

2sen lx

0

2

dx ak cos kx sen lx dx0

2

k1

bk sen kx sen lx dx0

2

k1

0

2


Distinguiamo i seguenti casi:1) k=l utilizzando la (1.3) e la (1.4), si ha: da cui:

2) k≠l questo non dà informazioni


0 0bl

bl 1

f x sen lx dx

0

2

Quindi i coefficienti di Fourier in forma trigonometrica sono:

al 1

f x cos lx dx

0

2

bl 1

f x sen lx dx

0

2

l 0,1,2...

l 1,2...


a0 1

f (x)dx

0

2

Polinomi trigonometrici di Fourier in forma complessa 2Consideriamo 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π

e ix cos x isen x

cos x e ix e ix

2

sen x e ix e ix

2i

x

k x e ikx

k 0,1,2,...

e ix cos x isen x

Polinomi trigonometrici di Fourier in forma complessa 2.1La 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 :

L2 0,2

L2 0,2 f : 0,2 C | f x 2dx

0

2

f x g x dx0

2

,

fL2 0, 2 f , f

fL2

0,2

Ff ˆ f kk ,k

Polinomi trigonometrici di Fourier in forma complessa 2.2

dove α(x) è la parte reale di f e β(x) è quella immaginaria.

f x x i x

x 0,2

e ikx cos kx isen kx

ak 1

2 x cos kx x sen kx dx

0

2

bk 1

2 x sen kx x cos kx dx

0

2

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

Infatti:

ˆ f k ak ibk

k 0,1,2,...

ˆ f k 1

2f x e ikxdx

0

2

1

2 x i x cos kx isen kx dx

0

2

1

2 0

2

x cos kx x sen kx i x sen kx i x cos kx dx

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.

Interpolazione trigonometrica 3.1Un 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)

f x a0

2 ak cos kx bksen kx

k1

M

x j 2 /(n 1), j 0,...,n

f (x j ) f (x j ), j 0,....,n

f

f

f

f x a0

2 ak cos kx bksen kx

k1

M

aM 1 cos((M 1)x)

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

f (x) ckeikx

k M

M

ak ck c k

e ikx cos(kx) isin(kx)

ckeikx ck (cos(kx) isin(kx))

k M

M

k M

M

k1

M

[ck (cos(kx) isin(kx))c k (cos(kx) isin(kx))]c0

Interpolazione trigonometrica 3.2

bk i(ck c k ),

k 0,1,......,M

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.

f (x) ckeikx

k (M 1)

M

cM 1 c (M 1) aM 1 /2

f (x) ckeikx

k (M )

M

Interpolazione trigonometrica 3.2

Comando interpft MATLABIl 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.

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 equispaziatie calcola i valori dell’interpolatore trigonometrico in 100 nodi equispaziati.La linea tratteggiata è la nostra funzione ed in corrispondenza della linea continual’interpolatore trigonometrico relativo a 10 nodi equispaziati.

f (x) x(x 2)e x

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 adapprossimazioni inaccurate. Solo aumentando il numero di nodi sarà possibile approssimare correttamente le funzioni di frequenza più elevata.

Trasformata di Fourier ed FFT

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.

Perché utilizzare l’analisi in frequenzeAd 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.

Spettro di un segnale armonico a frequenza equispaziateSe 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.

Serie di Fourier

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.

ikxeixe

2

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.

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).

Esempio

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

La Trasformata di Fourier

La trasformata di Fourier è l’operatore che permette di

trasformare il segnale originario espresso nel dominio del tempo

in una “somma” di segnali sinusoidali a tutte le frequenze.

Per l’approssimazione di segnali aperiodici si applica la trasformata di Fourier

La trasformata discreta di Fourier opera su sequenze di valori

discreti, equispaziati in un intervallo finito.

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 aumenta ogni volta di 2n+1.

12

2

njt j

,...2,1,0 j

DFT ed interpolazione

2...0 210 nttt jitj ez

jt

1.z

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:

Calcolo della DFT: cerchio unitario

2 2 1x y

Scriviamo il polinomio di interpolazione come:

e imponiamo le condizioni di interpolazione:

Si ha:

n

nk

iktkn ectp

jjn tftp ,...2,1,0 j

j

n

nk

iktk tfec j

nj 2,...,0


Per trovare i coefficienti , usiamo il seguente lemma.

LemmaPer ogni intero k :

2

0

2 1 1

0 1

j

j

j

ikt

nikt

iktj

n ee

e

kc


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:

j

n

nk

iktk tfec j

nj 2,...,0

jilte

nln

n

j

n

nkj

n

j

ilttlkik tfeec jj

2

0

2

0

lkn

lke

n

j

tlki j

12

02

0


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

n

j

n

nkj

n

j

ilttlkik tfeec jj

2

0

2

0

j

n

j

iltl tfe

nc j

2

012

1

nn cc ,...,

nnl ,...,

ntftf 20 ,...,

tpn


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

DFT E INTERPOLAZIONE

1

0

1 mkj

k j mj

d f wm

0,1,..., 1k m

mi

m ew/2

mw ( ) 1m

mw

kmk dd k

k

m

jj

jkm

m

jj

jmm

jkm

m

jj

mkjmmk dfw

mfww

mfw

md

1

0

1

0

1

0

1111

10 ,..., mdd

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 )

2 11 1 1

2 12 2 2

2 13 3 3

2 1

1

1

1

1

n

n

n

nm m m

V

1,

ji j ia 1

,i

i j ja

Consideriamo con

Possiamo scriverli in forma matriciale ottenendo:

Ponendo

con e ,

con

La DFT in forma matriciale sarà:

LA DFT IN FORMA MATRICIALE

1

0

1 mkj

k j mj

d a wm

0,1,..., 1k m

1

1

0

1

1

0

111

12

1

1

111

1

m

A

m

T

mmm

m

d

d

d

a

a

a

ww

ww

m

jkj

kjmkj w

mT

1 ,1,...,1,0, mjk

mmkj CT jj aA 1,...,1,0 mj

jkjk ATd

Poniamo per comodità. j jf a

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( )

Costo della DFT

0 1,...,

md d

2m

1

0

1 mkj

k j mj

d f wm

FFT

La FFT è un algoritmo sviluppato negli anni ’60 per il calcolo

veloce della DFT, grazie al quale si ha una riduzione rilevante del

numero di operazioni da compiere nel calcolo della DFT.

La FFT ci permette di ridurre il numero di operazioni da compiere

(in particolare il numero di moltiplicazioni). Questo è un aspetto

molto importante, soprattutto quando l’implementazione avviene su

un sistema che deve garantire un processing real-time.

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.

Cooley-Tukey

L'algoritmo FFT più diffuso è l'algoritmo di Cooley-Tukey.

Si basa sul principio di divide et impera: divide ricorsivamente

una DFT di qualsiasi dimensione N con N numero composto tale

che N=N1N2 in DFT più piccole di dimensioni N1 e N2, insieme a

O(n) moltiplicazioni per l'unità immaginaria, detti fattori twiddle.

Cooley-Tukey: Radix-2

L’algoritmo di cui parleremo è il Radix-2.

Caratteristiche dell’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

1. I campioni di indice pari nella prima sequenza

2. I campioni di indice dispari nella seconda sequenza moltiplicato per il

così detto fattore twiddle

2 k( i )K N

Nw e


Sia X[k] il coefficiente della DFT da calcolare. Suddiviamo la

sommatoria per il calcolo della DFT in due, sommando sugli indici

pari e sugli indici dispari:

Operiamo le sostituzioni n=2r per n pari e n=2r+1 per i dispari

ponendo


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.

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:

•

[ ] [ ] [ ]k

NX k G k W H k

2

NK K

N Nw w

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 è:

( ) [ ] [ ] 0,1,..., 1

2[ ] [ ] [ ]2

k

N

k

N

X k G k W H kN

kNX k G k W H k

Butterfly

• 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)

Decimazione del tempo per 8 punti

2( / 2)N

2( / 2) / 2N N

[ ]k

NW H k

La procedura di decimazione nel tempo risulta dunque vantaggiosa

e conviene riapplicarla al calcolo delle due DFT su N/2 punti G[k] e

H[k] della sottosequenza G[n] dei campioni di posto pari e H[n] di

campioni di posto dispari.

( /4 ) 1 ( /4 ) 1

/4 /2 /40 0

( /4 ) 1 ( /4 ) 1

/4 /2 /40 0

2 2 1

2 2 1

N Nnk k nk

N N Nn n

N Nnk k nk

N N Nn n

G k x r W W x r W

H k x r W W x r W

Decimazione del tempo: ricorsione

Analisi di complessità FFTSe operiamo ricorsivamente secondo la strategia appena

illustrata, dimezziamo la complessità ad ogni passo di ricorsione

dimezzando al contempo la lunghezza della sequenza numerica su

cui si opera. Così facendo, se :

quando la sequenza di ingresso si riduce ad un solo campione la

procedura ricorsiva è stata invocata un numero di volte pari a:

L = log2 N

Il calcolo della trasformata di un vettore a N componenti richiama

ricorsivamente il calcolo della trasformata di 2 vettori a N/2

componenti con O(N) operazioni di somma e prodotto aggiuntive.

2LN

Analisi di complessità FFT: equazione di ricorrenzaDetto T(N) il numero totale di operazioni per calcolare la trasformata

di un vettore a N componenti, vale allora la seguente equazione di

ricorrenza:

Poiché la soluzione alla precedente equazione di ricorrenza è T(N) =

O(N log N), l’algoritmo FFT permette una forte riduzione della

complessità in tempo per il calcolo della trasformata. In particolare,

l’analisi esposta sopra mostra che il calcolo della FFT richiede

esattamente (N/2) log2 N moltiplicazioni di numeri complessi.

Alcuni esempi in Matlab

Interpolatore trigonometrico

File: interp.m


Fenomeno di aliasing

File: aliasing.m


Approssimazione di un’onda quadra. Fenomeno di Gibbs

File: ondaquadra.m


Calcolo DFT ed FFT su f(x)= cos(2*pi*n/10).

File: esempio.m


Implementazione della DFT.

File: dft.m


Implementazione della FFT Radix2.

File: FFTRadix2.m


Confronto dei tempi di CPU tra DFT ed FFT

File: FFTflops.m

GRAZIE

Seminario Metodi Matematici per lOttimizzazione Anno accademico 2010/2011 Francesco Cambria Fabio...

Documents

Transcript of Seminario Metodi Matematici per lOttimizzazione Anno accademico 2010/2011 Francesco Cambria Fabio...