Diskretna furijeova transformacija

**Diskretna Furijeova transformacija** ili DFT jeste Furijeova transformacija diskretnog i konačnog (ili periodičnog) signala. Diskretna Furijeova transformacija je time i specijalni oblik Z-transformacije kod koje se z nalazi na jediničnom krugu. Često se koristi pri obradi digitalnih signala, a najpoznatiji algoritam za to je brza furijeova transformacija (FFT, Fast Fourier Transformation, engl.).

Diskretna Furijeova transformacija može da posluži takođe za aproksimaciju (u određenim slučajevima i rekonstrukciju) funkcije koja odgovara signalu ili kao implementacija digitalnih filtera.

Putem inverzne Furijeove transformacije se iz Furijeovih koeficijenata sklapa izlazni signal, a povezivanjem DFT i inverzne DFT možemo da manipulišemo frekvencijama (nalazi primenu pri ekvilajzerima i filterima).

## Definicija

Uzmimo da je R komutativan, unitaran prsten, u kojem je broj N jedinica. Dalje, u R je w jedinični koren.

Za vektor c=(c_0,\dots,c_{N-1})\in R^N je diskretna furijeova transformacija \hat c na sledeći način definisana:

\hat c_k = \sum_{j=0}^{N-1} w^{\,j\cdot k}\cdot c_j za k=0,\dots,N-1

A za \hat c_k, inverzna furijeova transformacija je

c_k = {1\over N}\sum_{j=0}^{N-1} w^{-j\cdot k}\cdot \hat c_j za k=0,\dots,N-1

### DFT i IDFT u kompleksnom domenu

U kompleksnom domenu koristimo w= e^\cdot c_j za k=0,\dots,N-1,

a IDFT za \hat c\in\mathbb C^N:
c_k=\frac 1 N \sum_{j=0}^{N-1} e^{2\pi \mathrm{i}\cdot\frac{jk}{N}}\cdot \hat c_j za k=0,\dots,N-1

### DFT i IDFT u realnom domenu

Računica u realnom domenu je:


\hat c_{N-k} 
= \sum_{j=0}^{N-1}e^{-2\pi \mathrm{i} \cdot\frac{j(N-k)}{N}} \cdot c_j
= \sum_{j=0}^{N-1}e^{-2\pi \mathrm{i} \cdot(\frac{jN}{N}-\frac{k}{N})} \cdot c_j
= \sum_{j=0}^{N-1}e^{-2\pi \mathrm{i} \cdot(\frac{jN}{N})}e^{2\pi \mathrm{i} \cdot(\frac{k}{N})} \cdot c_j
= \sum_{j=0}^{N-1}e^{-2\pi \mathrm{i} \cdot j}e^{2\pi \mathrm{i} \cdot(\frac{k}{N})} \cdot c_j


Ojlerov identitet glasi: e^{2 \pi \mathrm{i}}=1. Takođe važi \overline{e^{\mathrm{i}\phi}}=e^{-\mathrm{i}\phi} i \overline{z_1+z_2}=\overline{z_1} + \overline{z_2}.

Stoga možemo još uprostiti izraz:


\hat c_{N-k}
= \sum_{j=0}^{N-1}1 \cdot e^{2\pi \mathrm{i} \cdot(\frac{k}{N})} \cdot c_j
= \sum_{j=0}^{N-1}e^{2\pi \mathrm{i} \cdot\frac{k}{N}} \cdot c_j
= \sum_{j=0}^{N-1}\overline{e^{-2\pi \mathrm{i} \cdot\frac{k}{N}}} \cdot c_j
= \overline{\hat c_{k}}



Što će reći, \hat c\in\mathbb C^N nije realan, ali samo N nezavisnih vrednosti (umesto 2N).

Za IDFT možemo zaključiti sledeće: Ukoliko za \hat c\in\mathbb C^N važi \hat c_{N-k}=\overline{\hat c_k} za sve k=0,\dots,N-1, onda je IDFT realan vektor c\in\R^N.

### Pomeranje i skaliranje u vremenu i frekvenciji

Ako je signal periodičan, onda nije bitno da li transformišemo u opsegu 0,\dots,N-1 ili k,\dots,N-1+k. Indeksna promenljiva j treba da obuhvati N opseg, ali nije bitno gde on počinje odnosno gde se završava (ovo važi samo za slučaj da je signal periodičan, tj. da se vektor k,\dots,N-1+k periodično ponavlja). 
Prisetimo se: za w važi w^N = 1. Onda w^{N+k} = w^N \cdot w^k = 1 \cdot w^k = w^k.

U praksi često želimo da razlika u indeksima bude istovremeno i razlika u vremenu ili razdaljini dva merenja
:t_n = nT, n=1-M,\dots,N-M, T je perioda našeg merenja.

Često želimo i da koeficijentima dodelimo frekvenciju tako da su centrirane oko 0
:\omega_n:=2\pi\frac{n}{NT}, n=1-K,..., N-K, K je negde u blizini \frac{N}{2}.

Uzmimo neku funkciju f kojoj dodeljujemo x\in\mathbb C^N tako da x_n=f(t_n).

DFT je onda y_n=\hat f(\omega_n)=F(\omega_n).

Iz toga sledi:
:
F(\omega_n)=\sum_{j=1-M}^{N-M}e^{-2\pi \mathrm{i}\frac{njT}{NT}}x_j
=\sum_{j=1-M}^{N-M}e^{-\mathrm{i}\,\omega_n\cdot t_j}f(t_j)


a IDFT je 
:
f(t_n)=x_n=\frac1N \sum_{j=1-K}^{N-K}e^{-2\pi \mathrm{i}\frac{nkT}{NT}}y_j
=\frac1N \sum_{j=1-K}^{N-K}e^{\mathrm{i}\omega_j\cdot t_n}F(\omega_j)


## Primeri
### Primer filtera
Naš cilj je da izbacimo sve frekvencije koje su „tiše“ (tj. koje imaju amplitudu) od 1 V. 
Prvo pravimo tabelu:

Imamo 10 vrednosti na 1 sekundu, što će reći perioda našeg merenja je T = 0.1\,, a frekvencija freq = \frac{1}{T} = 10 Hz. Stoga mi možemo da rekonstruišemo talas do 5&nbsp;Hz. Ukoliko u našem originalnom signalu ima frekvencija viših od 5&nbsp;Hz onda će naša rekonstrukcija imati grešku. Ali, kao i uvek u životu, čovek mora biti optimističan te ćemo mi pretpostaviti da nema viših frekvencija (to je uostalom i jedan od razloga zašto kompakt-disk ima frekvenciju od oko 41&nbsp;kHz; ljudsko uho može da registruje tek do 20&nbsp;kHz!).

Sledi izračunavanje \omega_i\,. Nas zanimaju samo vrednosti vezane za pozitivne indekse:
:
n = 1-K,\ldots,N-K; K = N / 2 = 5;\,


:
n = -4,\ldots, 5\,


:
N \cdot T = 10 \cdot .1 = 1


:
\omega_0 = 2\pi\frac{0}{NT} = 0, \omega_1 = 2\pi, \omega_2=4\pi, \ldots, \omega_5=10\pi 


Sada imamo sve vrednosti i možemo da počnemo sa računanjem:
:
F(\omega_0 = 0) = \sum_{j=0}^{N-1=9} e^{ -\mathrm{i} \cdot 0 \cdot t_j } \cdot f(t_j )
= \sum_{j=0}^{9} f(t_j) = 100 = y_0 = {\hat c}_0 


:
F(\omega_1 = 2\pi) = \sum_{j=0}^{N-1=9} e^{ -\mathrm{i} \cdot 2\pi \cdot t_j } \cdot f(t_j) = \ldots = 0 -3.5\mathrm{i} = y_1 = {\hat c}_1 


:
F(\omega_2 = 4\pi) = \sum_{j=0}^{N-1=9} e^{ -\mathrm{i} \cdot 4\pi \cdot t_j } \cdot f(t_j) = \ldots = 15 +0\mathrm{i} = y_2 = {\hat c}_2 


Izračunavanje ostalih koeficijenata ide analogno, te ćemo ih ovde samo navesti kao rezultate:
:
F(\omega_3 = 6\pi) = 2.5 - 3\mathrm{i} = y_3 = {\hat c}_3 

:
F(\omega_4 = 8\pi) =  6.7586e-14 + 2.8240e-14i = y_4 = {\hat c}_4 


Imamo {\hat c}\,, sada želimo da izbacimo sve previše „tihe“ tonove. Trebaju nam c_k\,:

:
c = {\hat c} / N \Rightarrow c_i = {\hat c}_i / 10


c_i:\,
10 -0.35i 1.5 - 0i 0.25 - 0.3i 0 + 0i

Znamo da važi: c_0 = a_0, c_{i>0} = \frac{1}{2}(a_i - b_i \mathrm{i}). Na taj način možemo da izračunamo a_i\, i b_i\,:

:
a_0 = 10\,


:
a_1 - b_1 \mathrm{i} = -0.7 \mathrm{i} \Rightarrow a_1 = 0, b_1 = 0.7\,


Ostale amplitude:
:a_2 = 3\,
:a_3 = 0.5\,
:b_3 = 0.6\,
:a_4 = b_4 = 0\,

Iz a_4=b_4 = 0\, možemo da zaključimo da frekvencija od 4&nbsp;Hz nema u našem signalu. Često je vrlo zgodno navesti sve amplitude u grafikonu. Amplituda A_k\, za neku frekvenciju k je A_k = \sqrt{(a_k^2 + b_k^2 )}.

Sve a_i\, i b_i\, za koje važi \sqrt{(a_i^2 + b_i^2 )}  izbacujemo i na kraju dobijamo rekonstruisanu i obrađenu funkciju:

:f(t) = 10 + 3\cos(2 \omega t)\,

Sada možemo da ponovo da izračunamo f(t_i)\, ili da se poslužimo IDFT i tako prerađen signal snimimo u memoriju.

### Primer u C-u
 #include 
 #include 
 #include 
 
 #define pi 3.14159265
 #define N 1000
 #define T 0.001
 #define FREQ 25
 
 double my_function (double t )
 {
 	 /* violina svira ton od 25&nbsp;Hz */
 	 double ugaona_brzina = 2 * pi * FREQ;
 	 return 5 + 10 * cos(ugaona_brzina * t) + 15 * cos(2 * (ugaona_brzina * t) ) 
 			+ 20 * sin (3 * (ugaona_brzina * t) );
 
 }
 
 complex double get_fourier_coef (double omega_n, double* t, double* f  )
 {
 	 complex double coeff = 0;
 
 	int k = 0;
 
 	for (k = 0; k  nas signal je f_t !*/
 	t[0] = 0;
 	f[0] = my_function (t[0] );
 	omega[0] = 0;
 
 	for (n = 1; n  0.1 ){
 			printf ("# Koeficijent %d:  %e * e^i*%e\n", n, cabs(coeff[n]), carg(coeff[n]) ); 
 		}
 	}
 	
 	
 
 	/* krece inverzija: */		
 	a[0] = cabs(coeff[0]) / N;
 	phi[0] = 0;
 	for (n = 1; n  0.1 )
 		{
 			// c = 1/2 (a + ib ), zato a = 2 * |c|, b == 0 
 			a[n] = 2  * cabs(coeff[n]) / N; 
 
 			if (abs (carg(coeff[n])) > 0.001 )
 			{
 				phi[n] = carg(coeff[n] );
 			}
 			 else 
 			{
 				phi[n] = 0;
 			}
 		} 
 		else 
 		{
 			a[n] = 0;
 			phi[n] = 0;
 		}
 	}
 
 
 	/* predstavljanje rezultata: */
 	printf ("Nasa rekonstrukcija:\n f (t) = %e", a[0] );
 	for (n = 1; n < N/2+1; n++ )
 	{
 
 		if (a[n] )
 		{
 			if (phi[n] )
 			{
 				printf (" + %e * cos (%d * (2 * pi * t + %e) )", a[n], n, phi[n] );
 			}
 			else
 			{
 				printf ("+ %e * cos (%d * 2 * pi * t )", a[n], n );
 			}
 		}
 	}
 	printf ("\n" );
 	
 
 	return 0;
 }

## Literatura

* 
* 
* 
*  esp. section 30.2: The DFT and FFT, pp.&nbsp;830–838.
* 
* 
*  
* 
* 
* 
* 
* 
* 
* 
* 
* 

## Povezano
* Dualnost po Pontrjaginu