La trasformata di fourier

La serie di Fourier consente di lavorare con funzioni periodiche,
rappresentandole come coefficienti di una serie, ovvero come ampiezze di
composizioni di seni e coseni di pari frequenza (notare che questo
coincide alla soma di coseni di date ampiezza e fase).

La possibilità di manipolare una più ampia famiglia di funzioni
chiaramente estende questo strumento a molte più applicazioni, la
trasformata di Fourier elimina l'assunzione di periodicità. Per ottenere
questo risultato è sufficiente far tendere il periodo T all'infinito.
(Per i dettagli si veda il notebook Mathematica allegato.)

Data una funzione f(t) nel dominio del tempo, la sua rappresentazione
nel dominio della frequenza diventa:

f(w) = \int_{-\infty}^{+\infty} f(t) e^{-i w t} dt

viceversa:

f(t) = 1/(2 pi) \int_{-\infty}^{+\infty} f(w) e^{i w t} dw

Molte proprietà sono ereditate dalla serie di Fourier, ad esempio la
trasformata di Fourier di funzioni reali è simmetrica (pari).

Trasformata di funzioni discretizzate
La convoluzione:
Filtri ideali.. ma non particolarmente!
Trasformata di Fourier discreta (DFT)
FFT
Separabilità in 2D

Trasformata di funzioni discretizzate

Intuitivamente, se una funzione viene discretizzata (ammettiamo ad esempio
un passo di discretizzazione = 1), l'integrale da valutare per la
trasformata di Fourier torna ad essere una serie, il risultato è una
trasformata periodica:

% Una funzione di esempio non periodica:
t = -10:0.01:10;
f_t = 1./(1+t.^2);

% Discretizzazione:
N = 50;
f_tq = f_t;
for ind=1:numel(f_t); f_tq(ind)=f_t(N*floor(ind/N)+1); end

% La trasformata di Fourier (discreta, vedremo):
f_w = fft(f_t);
f_wq = fft(f_tq);

%  Notiamo che la periodicità esiste, l'utilizzo della trasformata discreta
% invece che continua ha prodotto una sorta di "attenuazione":
figure;
subplot(221); plot(t,f_t); title('1/{1+t^2}');
subplot(222); plot(abs(f_w(100:300))); title('DFT(1/{1+t^2})');
subplot(223); plot(t,f_tq); title('Q_{50}(1/{1+t^2})');
subplot(224); plot(abs(f_wq(100:300))); title('DFT(Q_{50}(1/{1+t^2}))');

La convoluzione:

Immaginiamo di avere un segnale f(t) rappresentato in frequenza come
f(w), immaginiamo di voler "filtrare" determinate frequenze, ovvero
rimuoverle moltiplicando f(w) per una funzione h(w) = 0 per le frequenze
da eliminare e 1 per le altre. Indipendentemente dall'utilità del filtro
h ci interessa salere cosa avviene nel dominio spaziale, ovvero quale sia
l'antitrasformata di f(w)h(w):

1/(2 pi) \int_{-\infty}^{+\infty} f(w)h(w) e^{i w t} dw
 = 1/(2 pi) \int_{-\infty}^{+\infty} [\int_{-\infty}^{+\infty} f(tau) e^{-i w tau} h(w) e^{i w t} dtau] dw
 = \int_{-\infty}^{+\infty} f(tau) [1/(2 pi)\int_{-\infty}^{+\infty} e^{-i w (t-tau)} h(w) dw] dtau
 = \int_{-\infty}^{+\infty} f(tau) h(t-tau) dtau = f(t)*h(t)

Chiamiamo questo integrale "di convoluzione", esso permette di applicare
nel dominio del tempo un filtro lineare a una funzione.
Si osservi che tramite un cambio di variabile si dimostra:

f(t)*H(t) = \int_{-\infty}^{+\infty} f(tau) h(t-tau) dtau
  = \int_{-\infty}^{+\infty} h(tau) f(t-tau) dtau = h(t)*f(t)

Inoltre, si vede che applicando il filtro h alla funzione delta di
dirac delta(t)=1 sse t=0, 0 altrimenti, si ottiene:

\int_{-\infty}^{+\infty} delta(tau) h(t-tau) dtau = h(t)

per questo motivo h(t) si dice "risposta all'impulso".

% Il filtro:
h_t = zeros(size(f_t));
from = round(numel(h_t)/2);
h_t(from-2:from+2) = [-1 0 3 0 -1];
h_w = fftshift(fft(fftshift(h_t)));

% Applichiamolo:
fh_t = conv(h_t,f_t,'same');
fh_w = h_w.*fftshift(f_w);
fh_t2 = ifft(fftshift(fh_w));

% I risultati sono identici:
figure;
subplot(211); plot(fh_t); title('fh_t');
subplot(212); plot(abs(fh_t2)); title('fh_t2');

Filtri ideali.. ma non particolarmente!

Proviamo ad utilizzare un filtro ideale in frequenza (passa-basso) per
attenuare le alte frequenze presenti nello spigolo della funzione rect:

% La funzione:
f_t = zeros(1,1000);
f_t(400:600) = 1;

% Il filtro:
h_w = zeros(size(f_t));
N = 50;
from = floor(numel(h_w)/2-N/2);
h_w(from:from+N-1) = 1;

% Applichiamolo:
f_w = fftshift(fft(f_t));
fh_w = f_w.*h_w;
fh_t = ifft(fftshift(fh_w));

% Vediamo:
figure;
subplot(211); plot(f_t);
subplot(212); plot(real(fh_t));

Le alte frequenze rimosse hanno lasciato traccia non solo in prossimità
delle discontinuità, l'influenza di un filtro lineare localizzato in
frequenza non è localizzato nel dominio temporale.

Vediamo lo stesso effetto nelle immagini:

% Carico una immagine:
img = imread('cameraman.tif');

% La trasformo:
imgf = fftshift(fft2(img));

% La filtro con un filtro passa-basso ideale:
H1 = false(size(img)); H2 = H1;
[X,Y] = meshgrid(-size(H1,1)/2:size(H1,1)/2-1,-size(H1,2)/2:size(H1,2)/2-1);
H1(sqrt(X.^2+Y.^2) < size(H1,1)/4) = true;
H2(sqrt(X.^2+Y.^2) < size(H2,1)/8) = true;
imgf1 = imgf.*H1;
imgf2 = imgf.*H2;

% Vediamo prima e dopo in frequenza:
figure;
subplot(131); imshowscale(log(abs(imgf))); title('Prima del filtro');
subplot(132); imshowscale(log(abs(imgf1)+eps)); title('Dopo il filtro 1');
subplot(133); imshowscale(log(abs(imgf2)+eps)); title('Dopo il filtro 2');

% Antitrasformata e visualizzazione:
img1 = ifft2(fftshift(imgf1));
img2 = ifft2(fftshift(imgf2));
figure;
subplot(131); imshowscale(img); title('Prima del filtro');
subplot(132); imshowscale(abs(img1)); title('Dopo il filtro 1');
subplot(133); imshowscale(abs(img2)); title('Dopo il filtro 2');

Si può cercare di creare filtri localizzati sia nel tempo che nelle
frequenza tassellando il piano tempo-frequenza, in quel caso si parla di
wavelet.

Trasformata di Fourier discreta (DFT)

In elaborazione dei segnali e delle immagini non abbiamo a che fare con
segnali continui ma a campionamento discreto, pertanto è importante
definire formalmente uno strumento simile alla trasformata di Fourier per
segnali continui che operi su segnali discreti. La definizione della DFT
proviene direttamente da quella della controparte continua considerando
la funzione analizzata campionata, sorvolando i dettagli algebrici,
considerando che l'integrale della trasformata di Fourier discretizzando
diviene una serie bilatera, otteniamo:

f(n) = \sum_{t=-\infty}^{+\infty} f(t) e^{-i 2 pi n/T t}

La cui antitrasformata diventa:

f(t) = 1/(2 pi) \sum_{n=-\infty}^{+\infty} f(n) e^{i 2 pi n/T t}

Se il segnale a disposizione è nullo per n<0 e n>T (ovvero di supporto
finito nel dominio del tempo) la DFT diviene:

f(n) = \sum_{t=0}^(T-1) f(t) e^{-i 2 pi n/T t}

e diviene algoritmicamente calcolabile:

% Inizializzazione:
f_t = zeros(1,513);
f_t(253:261) = 1;

% Trasformata "vanilla" della funzione rettangolo:
f_w = zeros(size(f_t));
T = numel(f_t);
for n=1:T
    f_w(n) = 0;
    for t=1:T
        % Manteniamo l'origine nel centro:
        tc = (t-floor(T/2)-1);
        f_w(n) = f_w(n) + f_t(t)*exp(-1i*2*pi*n/T*tc);
    end
end

% Facciamo lo shift "a mano":
from = floor(T/2);
f_w = [f_w(from:end),f_w(1:from-1)];

% Vediamone il risultato... ovvero la funzione sinc:
figure;
subplot(211); plot(f_t);
title('Funzione rect');
subplot(212); plot(real(f_w));
title('Funzione sinc come trasformata (reale) della funzione rect');

% La parte immaginaria è nulla a meno della precisione macchina:
fprintf('Massimo valore immaginario (in valore assoluto): %g\n',max(abs(imag(f_w))));

Massimo valore immaginario (in valore assoluto): 6.66134e-016

FFT

Runge-Konig nel 1924 e Cooley-Tukey nel 1965) proposero degli algoritmi
di calcolo della DFT con complessità O(N log(N)) basati sula metodologia
del divide-et-impera, l'idea è di ricondurre la DFT di un segnale
discreto finito di T campioni nella composizione di due DFT calcolate su
due sottoinsiemi di dimensione T/2.

Si veda fftNostra per l'algoritmo, qui un esempio di utilizzo:

% Una funzione con una potenza di 2 campioni:
% f_t = [ones(1,128),zeros(1,128)];
f_t = [ones(1,32),zeros(1,256-32)];

% Sua trasformata:
f_w = fftNostra(f_t);

% Mostriamo abs e angle:
figure;
subplot(211); plot(abs(f_w)); title('Modulo');
subplot(212); plot(angle(f_w)); title('Fase');

Separabilità in 2D

In 2 dimensioni la DFT assume la forma:

f(n,m) = \sum_{y=0}^(T-1)\sum_{x=0}^(T-1) f(x,y) e^{-i 2 pi (n x + m y)/T}

che separando gli esponenziali complessi diventa:

f(n,m) = \sum_{y=0}^(T-1) [\sum_{x=0}^(T-1) f(x,y) e^{-i 2 pi n/T x}] e^{-i 2 pi m/T y}

ovvero f(n,m) = FFT_y(FFT_x(f(x,y)) con FFT_x e FFT_y trasformate sulle x e sulle y rispettivamente.

Proviamo utilizzando la nostra fft su una immagine:

% La "funzione" immagine:
f_xy = checkerboard(16);
f_nm = f_xy;

% Trasformo le colonne:
for x=1:size(f_xy,1)
    f_nm(:,x) = fftNostra(f_nm(:,x));
end

% Trasformo le righe:
for y=1:size(f_xy,2)
    f_nm(y,:) = fftNostra(f_nm(y,:));
end

% FFT shift "a mano":
f_nm = [f_nm(:,size(f_nm,2)/2+1:end),f_nm(:,1:size(f_nm,2)/2)];
f_nm = [f_nm(size(f_nm,1)/2+1:end,:);f_nm(1:size(f_nm,1)/2,:)];

% Vediamo:
figure;
subplot(121); imshow(f_xy);
subplot(122); imshowscale(abs(f_nm));