Pubblicato da: Dieke | 11/05/2011

DFT, FFT, power spectrum, filtraggio

FFT – Introduzione

La trasformata di Fourier veloce (FFT: Fast Fourier Transform) è un algoritmo efficiente  per il calcolo della DFT ed è diventato uno strumento molto potente per l’analisi dei segnali.

La sua (ri)scoperta avvenne nel 1965 da parte di Cooley e Tukey, ma era già stata utilizzata un secolo prima da C.F. Gauss per velocizzare il calcolo delle orbite di asteroidi, sebbene il suo uso passò praticamente inosservato per lungo tempo.

Oggi è utilizzata in moltissimi ambiti, dalla compressione audio/video al filtraggio ed eliminazione del rumore, dal pattern matching allo studio della struttura cristallina fino alla tomografia assiale e all’elaborazione di immagini, analisi di onde sismiche, di serie storiche e immagini radar.

Discrete Fourier Transform (DFT)

Come è noto, la DFT di una sequenza x[n] di N valori è la sequenza complessa X[k] di N valori:

La DFT consente di passare dalla rappresentazione di un segnale nel dominio del tempo al dominio delle frequenze, rivelando la composizione delle diverse frequenze che contribuiscono a formare il segnale originale, ovvero il suo spettro.

La sequenza X[k] risultante è una quantità complessa (quindi scomponibile in ampiezza e fase oppure, in modo equivalente, in parte reale e immaginaria) e periodica con periodo N.

Nel caso di segnali reali, la FFT presenta una simmetria coniugata pari, cioè X[k] = X*[N-k], quindi il suo spettro è completamente determinato dagli N/2+1 valori X[0]…X[N/2]

Il calcolo diretto richiede un numero di moltiplicazioni proporzionali a N*N.

Usando la FFT per il calcolo della DFT, il numero di moltiplicazioni si riduce a N*Log2(N).

DFT Inversa

Per ricostruire il segnale originale x[n] dalla sua sequenza X[k] si esegue l’operazione inversa (IDFT):

Anche in questo caso il segnale risultante è una quantità complessa, periodica con periodo N, e può essere ottenuto tramite una FFT, con la sola modifica di un segno ed un fattore di correzione finale.

Esempio di FFT di segnali campionati

Proviamo a generare un semplice segnale e a calcolarne lo spettro utilizzando uno dei software a disposizione: Scilab, Ocatve piuttosto che Matlab, software commerciale della Mathworks (le notazioni usate sono praticamente identiche a parte minime variazioni).

Consideriamo un segnale campionato a 8000 Hz costituito da tre sinusoidi di ampiezza 5, 3, 2 e rispettive frequenze di 100 Hz, 200 Hz, 400 Hz e numero di campioni N=1024.

% Numero di campioni

N=1024

% Frequenza di campionamento

fs=8000

% Tempi di campionamento

t=0: 1/fs: (N-1)/fs;

% Segnale campionato

x=5*sin(2*pi*t*100)+3*sin(2*pi*t*200)+2*sin(2*pi*t*400);

Grafico del segnale campionato (dominio del tempo)

% Trasformata di fourier del segnale (DFT)
X=fft(x)

% Grafico della sequenza  FFT (Valore assoluto dei valori complessi)

plot(abs(X))

% Segnale reale: lo spettro presenta una simmetria rispetto al valore N/2

% Tutta l’informazione è contenuta nei valori da 0 a N/2 (i vettori in Matlab/Scilab iniziano da 1)

plot(abs(X(1:N/2+1)))

Spettro del segnale campionato (frequenze non normalizzate)

% Creo il vettore delle frequenze

freq=(fs/2)*(0:N/2)/(N/2);

plot(freq(1:100),abs(X(1:100)));

Spettro del segnale campionato – frequnze riscalate

Si vedono bene 3 picchi centrati sulle frequenze del segnale originale a 100, 200 e 400 Hz con le ampiezze relative.

Tuttavia si può notare che i picchi hanno un’ampiezza a metà altezza (FWHM) diversa da zero. CI saremmo aspettati invece 3 ‘delta di Dirac’ centrate esattamente sulle frequenze delle sinusoidi del segnale originale.

Vedremo in un prossimo post la causa di questo effetto (chiamato leakage) e come attenuarlo tramite windowing (hannning, hamming, blackman e altri tipi).

Macchie solari e ciclo dell’attività del sole

Da Internet si puo ottenere il file del numero medio mensile di macchie solari sintetizzato in un indice chiamato Wolf Number storicizzato negli ultimi 300 anni http://solarscience.msfc.nasa.gov/greenwch/spot_num.txt

Dopo aver salvato in locale il file ed averlo caricato con il comando: load(‘filepath’), il grafico è il seguente

Come in precedenza, calcolando la FFT della serie storica si ottiene il grafico dello spettro (sono visualizzati solo i valori corrispondenti alle frequenze piu basse)

Si vede chiaramente il picco posizionato intorno ai 0.0076 Hz ovvero (la frequenza di campionamento è di 1 mese = 1/12 di anno) 1/0.0076/12 = 10.9 anni circa.

Filtraggio nel dominio delle frequenze

Questa volta generiamo un segnale campionato a 8000 Hz composto dalle frequenze di 100, 200, 400 Hz e corrotto da rumore bianco (prendiamo i soliti 1024 campioni)

N=1024;

fs=8000;

t=0: 1/fs: (N-1)/fs;

x=5*sin(2*pi*t*100)+5*sin(2*pi*t*200)+6*sin(2*pi*t*400)+40*(rand(1,1024)-1);

Il grafico di una realizzazione del segnale è questo (altre realizzazioni saranno per forza diverse !)

Prendiamo la trasformata del segnale:

X=fft(x);

Il segnale presenta tre picchi abbastanza visibili a 100, 200 e 400 Hz ma sono visibili altre frequenze che occupano tutta la banda fino alla frequenza di Nyquist di 4000 Hz = fs/2 (nel grafico sono visualizzate solo le frequenze da 0 Hz a 1000 Hz). Il rumore bianco si estende omogeneamente su tutta la banda disponibile.

Proviamo a filtrare il segnale con un filtro passa-basso cercando di ridurre la presenza del rumore.

(Continua…)


Lascia un commento

Inserisci i tuoi dati qui sotto o clicca su un'icona per effettuare l'accesso:

Logo WordPress.com

Stai commentando usando il tuo account WordPress.com. Chiudi sessione / Modifica )

Foto Twitter

Stai commentando usando il tuo account Twitter. Chiudi sessione / Modifica )

Foto di Facebook

Stai commentando usando il tuo account Facebook. Chiudi sessione / Modifica )

Google+ photo

Stai commentando usando il tuo account Google+. Chiudi sessione / Modifica )

Connessione a %s...

Categorie

%d blogger cliccano Mi Piace per questo: