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} = \frac{1}{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^{\frac{2 * π * k}{n}}, k = 0, 1, \dots, n - 1$ .

Onda je DFT za $c \in ℂ^{N}$ : ${\hat{c}}_{k} = \sum_{j = 0}^{N - 1} e^{- 2 π i \cdot \frac{j k}{N}} \cdot c_{j}$ za $k = 0, \dots, N - 1$ ,

a IDFT za $\hat{c} \in ℂ^{N}$ : $c_{k} = \frac{1}{N} \sum_{j = 0}^{N - 1} e^{2 π i \cdot \frac{j k}{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 π i \cdot \frac{j (N - k)}{N}} \cdot c_{j} = \sum_{j = 0}^{N - 1} e^{- 2 π i \cdot (\frac{j N}{N} - \frac{k}{N})} \cdot c_{j} = \sum_{j = 0}^{N - 1} e^{- 2 π i \cdot (\frac{j N}{N})} e^{2 π i \cdot (\frac{k}{N})} \cdot c_{j} = \sum_{j = 0}^{N - 1} e^{- 2 π i \cdot j} e^{2 π i \cdot (\frac{k}{N})} \cdot c_{j}$

Ojlerov identitet glasi: $e^{2 π i} = 1$ . Takođe važi $\overline{e^{i ϕ}} = e^{- i ϕ}$ 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 π i \cdot (\frac{k}{N})} \cdot c_{j} = \sum_{j = 0}^{N - 1} e^{2 π i \cdot \frac{k}{N}} \cdot c_{j} = \sum_{j = 0}^{N - 1} \overline{e^{- 2 π i \cdot \frac{k}{N}}} \cdot c_{j} = \overline{{\hat{c}}_{k}}$

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

Za IDFT možemo zaključiti sledeće: Ukoliko za $\hat{c} \in ℂ^{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 ℝ^{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} = n T, 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$

ω_{n} : = 2 π \frac{n}{N T}, n = 1 - K, . . ., N - K

,

K

je negde u blizini

\frac{N}{2}

.

Uzmimo neku funkciju $f$ kojoj dodeljujemo $x \in ℂ^{N}$ tako da $x_{n} = f (t_{n})$ .

DFT je onda $y_{n} = \hat{f} (ω_{n}) = F (ω_{n})$ .

Iz toga sledi:

F (ω_{n}) = \sum_{j = 1 - M}^{N - M} e^{- 2 π i \frac{n j T}{N T}} x_{j} = \sum_{j = 1 - M}^{N - M} e^{- i ω_{n} \cdot t_{j}} f (t_{j})

a IDFT je

f (t_{n}) = x_{n} = \frac{1}{N} \sum_{j = 1 - K}^{N - K} e^{- 2 π i \frac{n k T}{N T}} y_{j} = \frac{1}{N} \sum_{j = 1 - K}^{N - K} e^{i ω_{j} \cdot t_{n}} F (ω_{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:

$t_{i} =$	0.5000	0.6000	0.7000	0.8000	0.9000	1.0000	1.1000	1.2000	1.3000	1.4000
$f (t_{i}) =$	12.5000	10.0995	7.6644	6.8554	9.7905	13.5000	11.7546	7.4815	8.2905	12.0636

Imamo 10 vrednosti na 1 sekundu, što će reći perioda našeg merenja je $T = 0.1$ , a frekvencija $f r e q = \frac{1}{T} = 10 H z$ . Stoga mi možemo da rekonstruišemo talas do 5 Hz. Ukoliko u našem originalnom signalu ima frekvencija viših od 5 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 kHz; ljudsko uho može da registruje tek do 20 kHz!).

Sledi izračunavanje $ω_{i}$ . Nas zanimaju samo vrednosti vezane za pozitivne indekse:

n = 1 - K, \dots, N - K; K = N / 2 = 5;

n = - 4, \dots, 5

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

ω_{0} = 2 π \frac{0}{N T} = 0, ω_{1} = 2 π, ω_{2} = 4 π, \dots, ω_{5} = 10 π

Sada imamo sve vrednosti i možemo da počnemo sa računanjem:

F (ω_{0} = 0) = \sum_{j = 0}^{N - 1 = 9} e^{- i \cdot 0 \cdot t_{j}} \cdot f (t_{j}) = \sum_{j = 0}^{9} f (t_{j}) = 100 = y_{0} = {\hat{c}}_{0}

F (ω_{1} = 2 π) = \sum_{j = 0}^{N - 1 = 9} e^{- i \cdot 2 π \cdot t_{j}} \cdot f (t_{j}) = \dots = 0 - 3.5 i = y_{1} = {\hat{c}}_{1}

F (ω_{2} = 4 π) = \sum_{j = 0}^{N - 1 = 9} e^{- i \cdot 4 π \cdot t_{j}} \cdot f (t_{j}) = \dots = 15 + 0 i = y_{2} = {\hat{c}}_{2}

Izračunavanje ostalih koeficijenata ide analogno, te ćemo ih ovde samo navesti kao rezultate:

F (ω_{3} = 6 π) = 2.5 - 3 i = y_{3} = {\hat{c}}_{3}

F (ω_{4} = 8 π) = 6.7586 e - 14 + 2.8240 e - 14 i = 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} i)$ . Na taj način možemo da izračunamo $a_{i}$ i $b_{i}$ :

a_{0} = 10

a_{1} - b_{1} i = - 0.7 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 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})} < 1$ izbacujemo i na kraju dobijamo rekonstruisanu i obrađenu funkciju:

f (t) = 10 + 3 \cos (2 ω 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 <stdio.h>
#include <math.h>
#include <complex.h>

#define pi 3.14159265
#define N 1000
#define T 0.001
#define FREQ 25

double my_function (double t )
{
	 /* violina svira ton od 25 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 < N; k ++ )
	{
		// f[k] == f(t[k] );
		coeff += cexp (- I * omega_n * t[k]) * f[ k] ;
	}
	return coeff;
}

int main()
{
	double t[N];
	double omega[N];
	double f[N];

	double a[N/2+1];
	double phi[N/2+1];
	int n = 0;

	complex double coeff[N];
	
	/* pripremi vektore t i f_t -> nas signal je f_t !*/
	t[0] = 0;
	f[0] = my_function (t[0] );
	omega[0] = 0;

	for (n = 1; n < N; n ++ )
	{
		omega[n] = 2 * pi * n / (N * T );
		t[n] = n * T;
		f[n] = my_function (t[n] );	
	}

 
	/* izracunavanje koeficijenata */
	for (n = 0; n < N/2+1; n ++ )
	{
		coeff[n] = get_fourier_coef (omega[n], t, f );
		if (cabs(coeff[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 < N/2+1; n++ )
	{
		if (cabs(coeff[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

Šablon:Refbegin

Šablon:Cite book
Šablon:Cite book
Šablon:Cite book
Šablon:Cite book esp. section 30.2: The DFT and FFT, pp. 830–838.
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal
Šablon:Cite journal

Šablon:Refend

Povezano

Dualnost po Pontrjaginu

Diskretna furijeova transformacija

Sadržaj

Definicija

DFT i IDFT u kompleksnom domenu

DFT i IDFT u realnom domenu

Pomeranje i skaliranje u vremenu i frekvenciji

Primeri

Primer filtera

Primer u C-u

Literatura

Povezano

Navigacija

Diskretna furijeova transformacija

Definicija

DFT i IDFT u kompleksnom domenu

DFT i IDFT u realnom domenu

Pomeranje i skaliranje u vremenu i frekvenciji

Primeri

Primer filtera

Primer u C-u

Literatura

Povezano

Navigacija

Pretraga