83075407-Trasformata-di-Fourier

5/16/2018 83075407-Trasformata-di-Fourier - slidepdf.com

http://slidepdf.com/reader/full/83075407-trasformata-di-fourier 1/7

DEV DEVeloping Software Solutions — n. 120 — giugno 2004

Trasformata di Fourierdi Cesare Bianchi

Una formula vecchia di due secoli e alla base di uno degli algoritmi piu usati nel processa-mento di segnali digitali

Cesare Bianchi

Si sta laureando in Psi-cologia e lavora come

programmatore ed in-segnante. Si occupa dimultimedia, intelligenzaartificiale e di varie al-tre cose riguardantil’informatica.



pubblicato su

WWW.INFOMEDIA.IT

stampa digitale da

Lulu Enterprises Inc.

stores.lulu.com / infomedia

Infomedia

Infomedia e l’impresa editoriale che da quasi venti an-

ni ha raccolto la voce dei programmatori, dei sistemi-

sti, dei professionisti, degli studenti, dei ricercatori e dei

professori d’informatica italiani.

Sono piu di 800 gli autori che hanno realizzato per le te-

state Computer Programming, Dev, Login, Visual Basic

Journal e Java Journal, molte migliaia di articoli tecnici,

presentazioni di prodotti, tecnologie, protocolli, strumen-ti di lavoro, tecniche di sviluppo e semplici trucchi e stra-

tagemmi. Oltre 6 milioni di copie distribuite, trentamila

pagine stampate, fanno di questa impresa la piu grande ed

influente realta dell’editoria specializzata nel campo della

programmazione e della sistemistica.

In tutti questi anni le riviste Infomedia hanno vissuto del-

la passione di quanti vedono nella programmazione non

solo la propria professione ma un’attivita vitale e un vero

divertimento.

Nel 2009, Infomedia e cambiata radicalmente adottando

un nuovo modello aziendale ed editoriale e si e organiz-

zata attorno ad una idea di Impresa Sociale di Comunita,

partecipata da programmatori e sistemisti, separando le

attivita di gestione dell’informazione gestite da un board

comunitario professionale e quelle di produzione gesti-

te da una impresa strumentale. Questo assetto e in linea

con le migliori esperienze internazionali e rende Infome-

dia ancora di piu parte della Comunita nazionale degli

sviluppatori di software.

Infomedia e media-partner di manifestazioni ed eventi in

ambito informatico, collabora con molti dei pi u impor-tanti editori informatici italiani come partner editoriale e

fornitore di servizi di localizzazione in italiano di testi in

lingua inglese.

L’impaginazione automatica di questa rivista e realizzata al

100% con strumenti Open Source usando OpenOffice,

Emacs, BHL, LaTeX, Gimp, Inkscape e i linguaggi Lisp,

Python e BASH

For copyright information about the contents of DEV, please see the section “Copyright” at the end of each ar-

ticle if exists, otherwise ask authors. Infomedia contents

is© 2004 Infomedia and released as Creative Commons

2.5 BY-NC-ND. Turing Club content is © 2004 Turing

Club released as Creative Commons 2.5 BY-ND.

Le informazioni di copyright sul contenuto di DEV so-

no riportate nella sezione “Copyright” alla fine di cia-

scun articolo o vanno richieste direttamente agli autori.

Il contenuto Infomedia e© 2004 Infomedia e rilasciato

con Licenza Creative Commons 2.5 BY-NC-ND. Il con-

tenuto Turing Club e © 2004 Turing Club e rilasciato

con Licenza Creative Commons 2.5 BY-ND. Si applicano

tutte le norme di tutela dei marchi e dei segni distintivi.

` E in ogni caso ammessa la riproduzione parziale o tota-

le dei testi e delle immagini per scopo didattico purch´ e

vengano integralmente citati gli autori e la completa

identificazione della testata.

Manoscritti e foto originali, anche se non pubblicati, non

si restituiscono.

Contenuto pubblicitario inferiore al 45%.

La biografia dell’autore riportata nell’articolo e sul

sito www.infomedia.it e di norma quella disponibi-

le nella stampa dell’articolo o aggiornata a cu-

ra dell’autore stesso. Per aggiornarla scrivere a

[email protected] o farlo in autonomia all’indirizzo

http: // mags.programmers.net / moduli / biografia

http://www.infomedia.it/

http://stores.lulu.com/infomedia




http://www.infomedia.it/



DEV > n. 119 giugno 2004 55 <<

di Cesare Bianchi > [email protected]

programmazione Multimedia

Una formula vecchia di due secoli è alla base di uno degli algoritmi più

usati nel processamento di segnali

digitali

J

ean Baptiste Joseph Fourier

(Figura 1), l’autore della famosa

trasformata, nacque nel 1768 ad

Auxerre in Francia. La sua vita non

fu esattamente quella che ci si

potrebbe immaginare per un ricer-

catore di fisica e matematica, ma d’altron-

de il luogo e i tempi non potevano che coin-

volgere il giovane Fourier, già studioso di

matematica, nella rivoluzione, e successi-

vamente spingerlo a seguire Napoleone

nella campagna d’Egitto come consigliere

scientifico. Tornato in Francia fu insignito

del titolo di Prefetto di Isére, e fu lì che tra

il 1804 e il 1807 concepì la sua più grande

scoperta e la sua più grande invenzione. La

prima è la famosa teoria sulla propagazione

del calore, la seconda, necessaria per

dimostrare la prima, è l’uso delle serie tri-

gonometriche per la trasformazione delle funzioni.

Serie e Trasformata di FourierLa rivoluzione introdotta nel calcolo da Fourier si basa sul-

l’assunzione che qualsiasi funzione periodica x(t) sia scom-

ponibile in una somma infinita di coseni, ovvero:

1) x(t) = A0 + Ak cos(2kf 0 t + k )

dove f 0 è la frequenza propria di x(t) , k è la fase di ogni

cosinusoide che compone la funzione e Ak è l’ampiezza di

tale cosinusoide.

Tale serie può essere generalizzata per funzioni non perio-

diche tramite la trasformata di Fourier:

2) x(t) = A(f )e i 2 ft df

e quindi

3) A(f ) = x (t )e –i 2 ft dt

che permette di trovare la funzione che

descrive le ampiezze delle frequenze, ovve-

ro, in termini matematici, passare dal domi-

nio del tempo al dominio delle frequenze.

Prima che qualche lettore si spaventi è il

caso di fare un veloce ripasso di trigono-

metria.

Come tutti ricorderanno dai tempi della

scuola, seno e coseno sono funzioni perio-

diche di periodo 2, ovvero si ripetono ogni

360°. Se vogliamo un periodo più corto (frequenza maggio-

re) o più lungo (frequenza minore) basta moltiplicare il 2

per la frequenza desiderata.Se l’inizio di un’onda sinusoidale è spostato rispetto allo 0,

tale spostamento si chiama fase.

Per capirci, immaginiamo un oggetto appeso ad una molla:

se lo lasciamo cadere dandogli una spinta la sua posizione

oscillerà arrivando anche più in alto di dove l’abbiamo

lasciato. Tale differenza è esattamente la fase di tale oscil-

lazione.

Inoltre, poiché sia seno che coseno hanno come massimo

1 e come minimo –1, se vogliamo che l’oscillazione sia più

ampia dovremo moltiplicarli per qualche fattore.

Sforzando un po’ la memoria possiamo anche ritrovare le

famose formule di prostaferesi, ovvero:

a) cos( + ) = cos()cos() – sin()sin()

e

b) sin( + ) = sin()cos() + cos()sin()

Facendo qualche calcolo potremmo renderci conto che

c) Acos(t +) = a cos(t ) – b sin(t ); a = Acos( ); b =Asin( )

ovvero qualsiasi oscillazione, anche con fase, può essere

espressa come somma di un seno e un coseno. Per quan-

to riguarda la notazione esponenziale, basta ricordare che

d) e it = cos(t ) + i sin(t )

dove i è l’unità immaginaria (ovvero la radice quadrata di

–1).

Trasformata di Fourier

FIGURA 1 Jean BaptisteJoseph Fourier

FIGURA 2Rappresentazione nello spazio complessodel vettore risultante dalle sommatoriedi seni e coseni

∞

k =1

∞

–∞

∞

–∞




56 DEV > n. 119 giugno 2004>>


Se aggiungiamo che l’integrale di una funzione ( e) f (x )dx )è la somma dei valori di tale funzione per ogni possibile valo-

re di x (non me ne vogliano i matematici per la genericità ditali affermazioni), capiamo finalmente cosa dice la trasfor-mata di Fourier.A(f) è perciò una funzione che descrive lo spettro delle variefrequenze che compongono la funzione x(t) . Per capirci, è lafunzione che indica l’altezza delle varie colonnine diWinamp. Se la musica che stiamo ascoltando ha un bel rif f di basso, vedremo le colonnine a sinistra (bassa frequenza)schizzare in alto, mentre se è Aretha Franklin saranno lecolonnine delle frequenze più alte a dover salire.La grande potenza della trasformata di Fourier è che, perinfinite frequenze (ed estensione infinita della funzione dastudiare) è possibile non solo scomporre una funzione, maanche ricomporla senza perdita di informazione.

Giusto per completezza è giusto aggiungere che A(f) hacodominio complesso, ovvero per ogni frequenza viene datal’ampiezza sia nel campo dei numeri reali che immaginari,e che lo stesso può valere per x(t) . Vedremo in seguito cosaquesto significhi e che utilità abbia.

Discrete Fourier TransformChi era già pronto ad obiettare che quanto finora detto valesolo per funzioni continue e che perciò nell’ambito infor-matico non è di grande utilità sarà stato di certo fermatodal titolo di questo paragrafo. Fatte le dovute modifiche,

infatti, la trasformata di Fourier si adatta perfettamente alprocessamento digitale dei segnali, e anzi è soprattutto nel-l’ambito discreto che è di maggiore utilità.Anche qui tuttavia occorre rinfrescare qualche concetto.Come tutti sapranno, il suono è una vibrazione che si pro-paga attraverso la materia. Se la “intercettiamo” con unmicrofono, verrà trasformata in un segnale elettrico, chesarà possibile registrare in due modi dif ferenti. Quello ana-logico prevede la trasposizione di tale segnale continuo suun supporto continuo (ad esempio un nastro magnetico), el’intensità del campo magnetico inciso sarà proporzionalealla potenza del segnale in ingresso. Nulla viene trasfor-mato in numeri.La registrazione digitale invece avviene attraverso un pro-cesso chiamato campionamento: ad ogni intervallo ditempo (nell’ordine di millesimi di secondo) viene misuratala potenza del segnale e convertita in un numero. Talesegnale sarà perciò discreto sia nel tempo che nell’am-piezza. La quantità di campioni presi per unità di tempoviene chiamata frequenza di campionamento, e ovviamenteall’aumentare di tale frequenza aumenta anche la qualitàdel suono che sarà possibile ricostruire.Perciò per descrivere il segnale invece di una funzione con-tinua x(t) useremo un vettore x[t] , che ci dirà ad ogni istan-te t (che sarà sempre un numero intero, da 0 a N–1, doveN è il numero totale dei campioni) il valore del segnale.La trasformata di Fourier è ancora valida, ma in una formaleggermente modificata, visto che abbiamo un numero dis-creto di campioni:

4) A(f ) = x [t ]e –i 2 ft = x [t ]cos(2 ft ) – ix [t ]sin(2 ft )

Cerchiamo a questo punto di capire meglio in termini prati-ci cosa questo voglia dire, e cosa indichino la componentereale e immaginaria dell’ampiezza di una data frequenza.Poniamo di voler sapere quanto di una data frequenza (adesempio il La con cui si accorda la chitarra: 440Hz) è pre-sente in un suono. Se abbiamo un diapason basterà espor-lo al suono e inizierà a vibrare in risonanza. Tanto più vibre-rà, tanto più in quel suono è rappresentata quella frequen-za. Ciò di cui non ci possiamo accorgere è che il diapason

vibrerà in fase con quella frequenza. Ovvero la sua funzio-ne sarà p(t) =A cos(2 ft+ ).

Poniamo invece ora il caso di un diapason speciale che nonpuò cambiare la sua fase. È il caso dei nostri cos(2 ft) esin(2 ft) . Possono entrare in risonanza, ma solo nella faseche loro compete. La fortuna è che tra loro sono sfasati di90°, e perciò, presi insieme, possono comunque acquisiretutta l’informazione relativa al segnale.Non è un caso, tra l’altro, che vengano rappresentati comele due dimensioni di uno spazio complesso.Come illustrato nella Figura 2, le due somme di seni ecoseni sono rappresentabili esse stesse come seno ecoseno, e il vettore risultante indica l’ampiezza (lunghezza)e la fase (angolo) della frequenza cercata.

Giusto per completezza è il caso di far notare che tale fun-zione non è normalizzata, ovvero all’aumentare di N aumen-ta anche il numero di contributi alla frequenza che vengonosommati, e quindi il totale complessivo. Per normalizzarlabasterà comunque dividere il risultato per N , così da poter

FIGURA 3 Rappresentazione grafica dei vettori dellevarie frequenze da cercare nel segnale

N –1

t =0N –1

t =0



DEV > n. 119 giugno 2004 57 <<



rendere confrontabili i risultati trovati per porzioni di suono

più o meno lunghe.

Ora che ci siamo un po’ chiariti le idee e che sappiamo

come scoprire quanto una data frequenza è rappresentata

in un suono, cerchiamo di capire come sfruttare questa

conoscenza.

Abbiamo detto che con la trasformata di Fourier è possibi-

le non solo scomporre, ma anche ricomporre una funzione

senza perdita di informazione. Lo stesso, posti certi vinco-

li, è vero per la trasformata discreta.

Se infatti prendiamo una porzione di suono composta da N

campioni, è possibile dimostrare che calcolando tutte le fre-

quenze f/N con f intero, 0 < f < N-1 abbiamo ricostruito,

nello spazio delle frequenze (quello delle colonnine di

Winamp, per capirci), la stessa informazione contenuta nei

campioni di suono originali.Nella Figura 3 sono rappresentate le varie frequenze che

andremo a “cercare” in un ipotetico segnale di soli 8 cam-

pioni (a sinistra la componente reale, ovvero il coseno, a

destra quella immaginaria, cioè il seno). Notare che benché

l’onda sia teoricamente continua, in realtà noi andremo a

“campionare” lo stato dell’oscillazione dell’ipotetico diapa-

son agli stessi istanti in cui abbiamo campionato il suono

originale. In questo modo otterremo 8 vettori, che moltipli-

cheremo con il vettore del suono, trovando così il contribu-

to che ogni campione di suono dà alle varie frequenze cer-

cate.

Poiché gli 8 vettori, descritti dalla formula

5) K f [t ] = e –i 2 ft/N

con 0 < f,t < N–1

(nel nostro caso N = 8 ) sono tra loro ortogonali (il prodotto

scalare tra due qualsiasi vettori dà sempre 0), e formano

perciò una perfetta base ortonormale in N dimensioni su

cui abbiamo proiettato gli N valori del nostro segnale, se

applichiamo la trasformata inversa otterremo nuovamente

il segnale originale.

L’utilità di questo algoritmo è rappresentata dalla possibili-

tà non solo di analizzare lo spettro delle frequenze che com-

pongono un dato suono, ma anche di costruire facilmente

dei filtri che eliminano alcune frequenze (ad esempio al di

sopra o al di sotto di un certo livello), per poi ricostruire il

segnale dalle frequenze rimanenti.

Fast Fourier TransformUno dei grandi limiti della Discrete Fourier Transform è nella

potenza di calcolo richiesta per computarla. Sono infatti

necessarie, per trovare le N frequenze contenute in un

frammento di suono, N somme per ogni frequenza, ogni

somma composta da 2 prodotti, in ognuno dei quali c’è una

funzione trigonometrica. Perciò il numero di calcoli richiesti

è proporzionale al quadrato di N.

Per ovviare a questo problema è stata inventata nel 1965

da James W. Cooley e John W. Tukey la Fast FourierTransform, un algoritmo che sfrutta alcune proprietà dei

vettori usati nel calcolo per ridurre drasticamente il numero

di computazioni necessarie.

Se definiamo

6) W N = e –i 2 /N

possiamo infatti riscrivere la (5) come

7) K f [t ] = e –i 2 ft/N = W N tf

e quindi anche la (4) come

8) A[f ] = x [t ]e –i 2 ft/N = x [t ]W N tf

Fin qui ci siamo limitati a riscrivere qualche formula, ma il

divertimento deve ancora arrivare. Ora infatti consideriamo

cosa accade a calcolare separatamente le somme dei

campioni pari (t=0,2,4,...) e dispari (t=1,3,5 ,...).

A[f ] = x [2r ]W N f (2r ) + x [2r +1]W N

f (2r +1) =

= x [2r ]W N f (2r ) + x [2r +1]W N

f (2r )W N f =

9)

= x [2r ]W N f (2r ) + W N

f x [2r +1]W N f (2r )=

= x [2r ]W N/2 f r + W N

f x [2r +1]W N/2 f r

FIGURA 4Rappresentazione grafica di unatrasformata di Fourier a 2 valori:schema “a farfalla”

FIGURA 5 Rappresentazione grafica di unatrasformata di Fourier a 4 valori

N –1

t =0

N /2–1

r =0

N /2–1

r =0

N /2–1

r =0

N /2–1

r =0

N /2–1

r =0

N /2–1

r =0

N /2–1

r =0

N /2–1

r =0

N –1

t =0




58 DEV > n. 119 giugno 2004>>


dove abbiamo usato una identità fondamentale:

10) W N f (2r ) = e (–i 2 )(2fr/N ) = e W N /2

f r

A questo punto possiamo notare una proprietà interessan-

te: le due sommatorie non sono nient’altro che due tra-

sformate di Fourier, calcolate l’una sui numeri pari, l’altra

su quelli dispari. Per N potenza di 2 possiamo applicarericorsivamente questo metodo fino ad arrivare a trasforma-

te di due soli valori.

Proviamo invece ora a partire dalla trasformata di due valo-

ri e risalire:

A[0] = x [0]W 2 0 + x [1]W 2

0 = x [0] + x [1]W 2 0

11)

A[1] = x [0]W 2 0 + x [1]W 2

1 = x [0] + x [1]W 2 1

Possiamo rappresentare graficamente queste due equazio-

ni con la famosa “farfalla” di Figura 4. Volendo continuare,

una trasformata di quattro valori può essere ridotta a due

trasformate di due valori, una per i valori pari, una per quel-

li dispari. Gli elementi dispari, come ricordiamo, andrannoperò poi moltiplicati per W4

f .

Graficamente questo può essere rappresentato come due

livelli di “farfalle” (Figura 5). Tra l’altro, utilizzando l’identi-

tà WN/2n = WN

2n possiamo esprimere tutti i moltiplicatori

come potenze dello stesso WN (nell’ultimo caso N = 4).

Andando avanti diventa ovvio che tutte le “farfalle” hanno

forma simile (vedi Figura 6). Utilizzando l’identità WNs+N/2 =

= WNsWN

N/2 = -WNs (vera in quanto WN

N/2 = e-2i(N/2)/N = e-i =

= cos(–) + i sin(–) = –1 ) possiamo semplificare ulterior-

mente il grafico (Figura 7).

Ed eccoci finalmente al cuore della FFT. Il trucco principale

è che non occorre calcolare ogni componente della trasfor-

mata separatamente: far ciò significherebbe ripetere inutil-

mente un gran numero di calcoli. Invece dividiamo i nostricalcoli in stadi. Ad ogni stadio abbiamo N numeri che

appaiamo tramite il sistema della farfalla per ottenere un

nuovo insieme di N numeri (complessi), i quali diventano a

loro volta l’input dello stadio successivo. Poiché ogni stadio

necessita di N moltiplicazioni, N/2 inversioni di segno e 2N

addizioni, e che il numero di stadi è log2N, il numero di cal-

coli necessari è proporzionale a Nlog(N) (per N grande c’è

un’enorme differenza rispetto all’N2 della DFT).

Inoltre i calcoli possono essere fatti usando un unico buf-

fer di N numeri complessi: il trucco sta nell’inizializzarequesto buffer con campioni appropriatamente ordinati (per

N = 4 l’ordine sarà x[0], x[2], x[1], x[3]). In generale, baste-

rà dividere i campioni in due gruppi, pari e dispari, e appli-

care questa divisione in modo ricorsivo. Giusto per diverti-

mento vi faccio notare che se scriviamo in codice binario i

numeri di una sequenza così ordinata scopriremo che sono

in ordine inverso di bit: ad esempio una sequenza di otto

numeri sarà 000, 100, 010, 110, 001, 101, 011, 111.

Ricapitolando, ecco come implementare un algoritmo FFT

(per la precisione, questo è l’algoritmo FFT a decimazione

nel tempo, un’altra classe di FFT usa invece la decimazio-

ne nella frequenza):

1) Prendiamo un numero N di campioni potenza di 2

2) Li inseriamo in un buffer di numeri complessi di

grandezza N (lasciando le parti immaginarie a zero)

3) Ordiniamo il buffer in ordine inverso di bit

4) Calcoliamo il primo stadio di “farfalle” usando paia

di numeri adiacenti nel buffer

5) Calcoliamo il secondo stadio di “farfalle” stavolta

usando paia di numeri che sono separati di 26) Calcoliamo il terzo stadio su paia di numeri separa-

ti di 4

7) Continuiamo a calcolare fino a quando non abbiamo

una separazione di N/2

8) Il buffer contiene la nostra trasformata di Fourier.

Nel Listato 1 c’è un esempio di questo algoritmo imple-

mentato in Visual Basic. La variabile STEP è la “larghezza”

FIGURA 6 Forma standard delle farfallee semplificazione

FIGURA 7 Trasformata di Fourier dopo lasemplificazione: FFT a 4 valori

–i 2 fr

N /2

La Trasformata di Fourier permette

di trovare la funzione che descrive

le ampiezze delle frequenze



DEV > n. 119 giugno 2004 59 <<



della farfalla, ovvero la distanza tra i due ingressi. Inizia con

1 e raddoppia ad ogni stadio. Ad ogni stadio perciò dobbia-

mo calcolare STEP insiemi di farfalle, formati da farfalle

separate tra loro di una distanza pari ad INCREM (che è il

doppio di STEP ). Il calcolo è organizzato in vari insiemi poi-

ché ogni insieme condivide lo stesso moltiplicatore U (divi-

so nella sue parte reale Ur e immaginaria Ui ). Alla fine del-

l’algoritmo vengono inoltre calcolate le ampiezze e le fasi

dello spettro trovato, che non sono nient’altro che il modu-

lo del vettore complesso (ricordate la Figura 2?) e l’angolo

della sua “pendenza” rispetto all’asse reale.

Insomma, è stata dura ma abbiamo finalmente capito comefunziona uno degli algoritmi più utilizzati nel processamento

di segnali audio (e non solo). Nelle prossime puntate, molto

meno matematiche, scopriremo invece alcune semplici

applicazioni di questo algoritmo.

BibliografiaSu Joseph Fourier: http://www-history.mcs.st-andrews.ac.

uk/Mathematicians /Fourier.html

Sulla FFT: James W. Cooley and John W. Tukey, An algo-

rithm for the machine calculation of complex Fourier series,

Math. Compu t. 19, 297–301 (1965).

http://en.wikipedia. org/wiki/Fast_Fourier_transform

Si sta laureando in Psicologia e lavora come programmato-

re ed insegnante. Si occupa di multimedia, intelligenza arti-

ficiale e di varie altre cose riguardanti l’informatica.

LISTATO 1Implementazione in Visual Basic della FFT

Option Explicit

Public Const PI# = 3.14159265358979 ‘ 4# * Atn(1#)

Public Const PI2# = 2# * PIPublic Const Log2# = 0.693147180559945 ‘ Log(2#)Public Const Log10# = 2.30258509299405 ‘ Log(10#)Public Const Sqr2# = 1.4142135623731 ‘ Sqr(2#)

Sub FFT(D() As Double, R() As Double, X() As Double _, S() As Double, F() As Double, ByRef NVAL As Long, ByRef

NFRE As Long)

‘ Entra con:‘ D() = Vettore dei dati Reali del Segnale(t).‘ NVAL = Numero max. dei dati in D() da usare.‘‘ Esce con:‘ R() = Vettore dei valori Reali della Trasformata(f).‘ X() = Vettore dei valori Immaginari della Trasformata(f).‘ S() = Vettore dello Spettro di Amp. della Trasformata(f).‘ F() = Vettore delle Fasi della Trasformata(f)

[Radianti].‘ NVAL = Numero dei dati usati dalla Trasformata(f).‘ NFRE = Numero dei valori di Frequenza calcolati.‘‘ I valori della Trasformata in R(), X() NON vengono‘ normalizzati sulla lunghezza del Segnale 2^M.‘‘ I vettori R(), X(), S() e F() sono dimensionati in questa

routine.‘ Gli NN dati nei vettori sono organizzati come Vettore(0 To

NN - 1).

Dim M, I, NMN1, J, K, L, INCREM, STEP, IP, JF As LongDim Ur, Ui, Wr, Wi, Tr, Ti As Double

M = Int(Log(CDbl(NVAL) + 0.5) / Log2)NVAL = 2 ^ MNFRE = NVAL / 2NMN1 = NVAL - 1

ReDim R(0 To NVAL - 1), X(0 To NVAL - 1)

ReDim S(0 To NFRE), F(0 To NFRE)

For I = 0 To NVAL - 1R(I) = D(I)

Next I

‘ Riordino il buffer in ordine inverso di bitJ = 0For I = 0 To NMN1 - 1

If (I < J) ThenSwap R(J), R(I)

End IfK = NFREDo While (K - 1 < J)

J = J - KK = K / 2

LoopJ = J + K

Next I

‘ Calcolo la trasformataFor L = 1 To M

INCREM = 2 ^ L

STEP = INCREM / 2Ur = 1#Ui = 0#Wr = Cos(PI / STEP)Wi = Sin(PI / STEP)For J = 0 To STEP - 1

For I = J To NVAL - 1 Step INCREMIP = I + STEPTr = R(IP) * Ur - X(IP) * UiTi = R(IP) * Ui + X(IP) * UrR(IP) = R(I) - TrX(IP) = X(I) - TiR(I) = R(I) + TrX(I) = X(I) + Ti

Next ITr = UrTi = UiUr = Tr * Wr - Ti * WiUi = Tr * Wi + Ti * Wr

Next JNext L

‘ Calcolo dello Spettro di Ampiezza e delle FasiS(0) = Abs(R(0)) / CDbl(NVAL)F(0) = IIf(R(0) < 0#, PI, 0#)For JF = 1 To NFRE - 1

S(JF) = Sqr(R(JF) * R(JF) + X(JF) * X(JF)) / CDbl(NFRE)F(JF) = Atn4Q(-X(JF), R(JF))

Next JFS(NFRE) = Abs(R(NFRE)) / CDbl(NVAL)F(NFRE) = IIf(R(NFRE) < 0#, PI, 0#)

End Sub

Public Function Atn4Q(ByVal Y#, ByVal X#) As Double‘ Ritorna il Valore dell’ ArcoTangente di y/x su 4 Quadranti.‘ E’: -PI < Atn4Q <= PI.

Select Case XCase Is > 0#Atn4Q = Atn(Y / X)

Case Is < 0#If Y >= 0# Then

Atn4Q = Atn(Y / X) + PIElseIf Y < 0# Then

Atn4Q = Atn(Y / X) - PIEnd If

Case Is = 0#Atn4Q = Sgn(Y) * PI / 2#

End SelectEnd Function

Public Sub Swap(vN1 As Variant, vN2 As Variant)Dim vTn As VariantvTn = vN1vN1 = vN2vN2 = vTn

End Sub

Cesare Bianchi

83075407-Trasformata-di-Fourier

Documents

Transcript of 83075407-Trasformata-di-Fourier