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 uredi

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 uredi

U kompleksnom domenu koristimo $w=e^{{2*\pi *k} \over n},\quad k=0,1,\ldots ,n-1$ .

Onda je DFT za $c\in {\mathbb {C}}^{N}$ : ${\hat {c}}_{k}=\sum _{j=0}^{N-1}e^{-2\pi \mathrm {i} \cdot {\frac {jk}{N}}}\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 uredi

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 \mathbb {R} ^{N}$ .

Pomeranje i skaliranje u vremenu i frekvenciji uredi

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}={\frac {1}{N}}\sum _{j=1-K}^{N-K}e^{-2\pi \mathrm {i} {\frac {nkT}{NT}}}y_{j}={\frac {1}{N}}\sum _{j=1-K}^{N-K}e^{\mathrm {i} \omega _{j}\cdot t_{n}}F(\omega _{j})$

Primeri uredi

Primer filtera uredi

Situacija: Zvuk koji želimo da snimimo ima sledeći oblik (kada bi ga snimao analogan mikrofon):

Pošto je naš mikrofon digitalan, mi možemo samo da snimimo pojedinačne vrednosti. Na našem kompjuteru dobijamo:

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

<math>t_i =\,</math>0.5000    0.6000    0.7000    0.8000    0.9000    1.0000    1.1000    1.2000    1.3000    1.4000

<math>f(t_i) =\,</math>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 $freq={\frac {1}{T}}=10Hz$ . 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 $\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 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})}}$ . U našem slučaju naš frekventni spektrum izgleda ovako:

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\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 uredi

#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;
}

Reference uredi

Brigham, E. Oran (1988). The fast Fourier transform and its applications. Englewood Cliffs, N.J.: Prentice Hall. ISBN 978-0-13-307505-2.
Oppenheim,, Alan V.; Schafer, Ronald W.; Buck, J. R. (1999). Discrete-time signal processing. Upper Saddle River, N.J.: Prentice Hall. ISBN 978-0-13-754920-7.
Smith, Steven W. (1999). „Chapter 8: The Discrete Fourier Transform”. The Scientist and Engineer's Guide to Digital Signal Processing (Second izd.). San Diego, Calif.: California Technical Publishing. ISBN 978-0-9660176-3-2.
Thomas H. Cormen; Leiserson, Charles E.; Rivest, Ronald L. & Stein, Clifford (2001). „Chapter 30: Polynomials and the FFT”. Introduction to Algorithms (Second izd.). MIT Press and McGraw-Hill. str. 822-848. ISBN 978-0-262-03293-3. esp. section 30.2: The DFT and FFT, pp. 830–838.
Duhamel, P.; Piron, B. & Etcheto, J. M. (1988). „On computing the inverse DFT”. IEEE Trans. Acoust., Speech and Sig. Processing. 36 (2): 285—286. doi:10.1109/29.1519.
McClellan, J. H. & Parks, T. W. (1972). „Eigenvalues and eigenvectors of the discrete Fourier transformation”. IEEE Trans. Audio Electroacoust. 20 (1): 66—74. doi:10.1109/TAU.1972.1162342.
Dickinson, Bradley W. & Steiglitz, Kenneth (1982). „Eigenvectors and functions of the discrete Fourier transform”. IEEE Trans. Acoust., Speech and Sig. Processing. 30 (1): 25—31. doi:10.1109/TASSP.1982.1163843.
Grünbaum, F. A. (1982). „The eigenvectors of the discrete Fourier transform”. J. Math. Anal. Appl. 88 (2): 355—363. doi:10.1016/0022-247X(82)90199-8.
Atakishiyev, Natig M. & Kurt Bernardo Wolf (1997). „Fractional Fourier-Kravchuk transform”. J. Opt. Soc. Am. A. 14 (7): 1467—1477. doi:10.1364/JOSAA.14.001467.
C. Candan; M. A. Kutay and H. M.Ozaktas (2000). „The discrete fractional Fourier transform”. IEEE Trans. on Signal Processing. 48 (5): 1329—1337. doi:10.1109/78.839980.
Magdy Tawfik Hanna, Nabila Philip Attalla Seif, and Waleed Abd El Maguid Ahmed (2004). „Hermite-Gaussian-like eigenvectors of the discrete Fourier transform matrix based on the singular-value decomposition of its orthogonal projection matrices”. IEEE Trans. Circ. Syst. I. 51 (11): 2245—2254. doi:10.1109/TCSI.2004.836850.
Gurevich, Shamgar & Hadani, Ronny (2009). „On the diagonalization of the discrete Fourier transform”. Applied and Computational Harmonic Analysis. 27 (1): 87—99. arXiv:0808.3281  . doi:10.1016/j.acha.2008.11.003. preprint at. Tekst „arXiv

” ignorisan (pomoć)

Gurevich, Shamgar; Hadani, Ronny & Sochen, Nir (2008). „The finite harmonic oscillator and its applications to sequences, communication and radar”. IEEE Transactions on Information Theory. 54 (9): 4239—4253. arXiv:0808.1495  . doi:10.1109/TIT.2008.926440. preprint at. Tekst „arXiv

” ignorisan (pomoć)

Vargas-Rubio, Juan G. & Santhanam, Balu (2005). „On the multiangle centered discrete fractional Fourier transform”. IEEE Sig. Proc. Lett. 12 (4): 273—276. doi:10.1109/LSP.2005.843762.
J. Cooley, P. Lewis, and P. Welch (1969). „The finite Fourier transform”. IEEE Trans. Audio Electroacoustics. 17 (2): 77—85. doi:10.1109/TAU.1969.1162036.
F.N. Kong (2008). „Analytic Expressions of Two Discrete Hermite-Gaussian Signals”. IEEE Trans. Circuits and Systems –II: Express Briefs. 55 (1): 56—60. doi:10.1109/TCSII.2007.909865.

Vidi još uredi

Dualnost po Pontrjaginu