Diskretna furijeova transformacija

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

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

${\displaystyle {\hat {c}}_{k}=\sum _{j=0}^{N-1}w^{\,j\cdot k}\cdot c_{j}}$ za ${\displaystyle k=0,\dots ,N-1}$

A za ${\displaystyle {\hat {c}}_{k}}$ , inverzna furijeova transformacija je

${\displaystyle c_{k}={1 \over N}\sum _{j=0}^{N-1}w^{-j\cdot k}\cdot {\hat {c}}_{j}}$ za ${\displaystyle k=0,\dots ,N-1}$

DFT i IDFT u kompleksnom domenu

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

Onda je DFT za ${\displaystyle c\in \mathbb {C} ^{N}}$ :${\displaystyle {\hat {c}}_{k}=\sum _{j=0}^{N-1}e^{-2\pi \mathrm {i} \cdot {\frac {jk}{N}}}\cdot c_{j}}$ za ${\displaystyle k=0,\dots ,N-1}$ ,

a IDFT za ${\displaystyle {\hat {c}}\in \mathbb {C} ^{N}}$ :${\displaystyle c_{k}={\frac {1}{N}}\sum _{j=0}^{N-1}e^{2\pi \mathrm {i} \cdot {\frac {jk}{N}}}\cdot {\hat {c}}_{j}}$ za ${\displaystyle k=0,\dots ,N-1}$

DFT i IDFT u realnom domenu

Računica u realnom domenu je:

${\displaystyle {\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: ${\displaystyle e^{2\pi \mathrm {i} }=1}$ . Takođe važi ${\displaystyle {\overline {e^{\mathrm {i} \phi }}}=e^{-\mathrm {i} \phi }}$ i ${\displaystyle {\overline {z_{1}+z_{2}}}={\overline {z_{1}}}+{\overline {z_{2}}}}$ .

Stoga možemo još uprostiti izraz:

${\displaystyle {\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, ${\displaystyle {\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 ${\displaystyle {\hat {c}}\in \mathbb {C} ^{N}}$ važi ${\displaystyle {\hat {c}}_{N-k}={\overline {{\hat {c}}_{k}}}}$ za sve ${\displaystyle k=0,\dots ,N-1}$ , onda je IDFT realan vektor ${\displaystyle c\in \mathbb {R} ^{N}}$ .

Pomeranje i skaliranje u vremenu i frekvenciji

Ako je signal periodičan, onda nije bitno da li transformišemo u opsegu ${\displaystyle 0,\dots ,N-1}$ ili ${\displaystyle 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 ${\displaystyle k,\dots ,N-1+k}$ periodično ponavlja). Prisetimo se: za ${\displaystyle w}$ važi ${\displaystyle w^{N}=1}$ . Onda ${\displaystyle 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

${\displaystyle t_{n}=nT,n=1-M,\dots ,N-M}$ , ${\displaystyle T}$ je perioda našeg merenja.

Često želimo i da koeficijentima dodelimo frekvenciju tako da su centrirane oko ${\displaystyle 0}$

${\displaystyle \omega _{n}:=2\pi {\frac {n}{NT}},n=1-K,...,N-K}$ , ${\displaystyle K}$ je negde u blizini ${\displaystyle {\frac {N}{2}}}$ .

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

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

Iz toga sledi:

${\displaystyle 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

${\displaystyle 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})}$

Primer filtera

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

${\displaystyle t_{i}=\,}$	0.5000	0.6000	0.7000	0.8000	0.9000	1.0000	1.1000	1.2000	1.3000	1.4000
${\displaystyle 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 ${\displaystyle T=0.1\,}$ , a frekvencija ${\displaystyle 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 ${\displaystyle \omega _{i}\,}$ . Nas zanimaju samo vrednosti vezane za pozitivne indekse:

${\displaystyle n=1-K,\ldots ,N-K;K=N/2=5;\,}$

${\displaystyle n=-4,\ldots ,5\,}$

${\displaystyle N\cdot T=10\cdot .1=1}$

${\displaystyle \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:

${\displaystyle 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}}$

${\displaystyle 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}}$

${\displaystyle 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:

${\displaystyle F(\omega _{3}=6\pi )=2.5-3\mathrm {i} =y_{3}={\hat {c}}_{3}}$

${\displaystyle F(\omega _{4}=8\pi )=6.7586e-14+2.8240e-14i=y_{4}={\hat {c}}_{4}}$

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

${\displaystyle c={\hat {c}}/N\Rightarrow c_{i}={\hat {c}}_{i}/10}$

${\displaystyle c_{i}:\,}$ 10 -0.35i 1.5 - 0i 0.25 - 0.3i 0 + 0i

Znamo da važi: ${\displaystyle 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 ${\displaystyle a_{i}\,}$ i ${\displaystyle b_{i}\,}$ :

${\displaystyle a_{0}=10\,}$

${\displaystyle a_{1}-b_{1}\mathrm {i} =-0.7\mathrm {i} \Rightarrow a_{1}=0,b_{1}=0.7\,}$

Ostale amplitude:

${\displaystyle a_{2}=3\,}$

${\displaystyle a_{3}=0.5\,}$

${\displaystyle b_{3}=0.6\,}$

${\displaystyle a_{4}=b_{4}=0\,}$

Iz ${\displaystyle 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 ${\displaystyle A_{k}\,}$ za neku frekvenciju ${\displaystyle k}$ je ${\displaystyle A_{k}={\sqrt {(a_{k}^{2}+b_{k}^{2})}}}$ .

Sve ${\displaystyle a_{i}\,}$ i ${\displaystyle b_{i}\,}$ za koje važi ${\displaystyle {\sqrt {(a_{i}^{2}+b_{i}^{2})}}<1}$ izbacujemo i na kraju dobijamo rekonstruisanu i obrađenu funkciju:

${\displaystyle f(t)=10+3\cos(2\omega t)\,}$

Sada možemo da ponovo da izračunamo ${\displaystyle 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 25double 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;}

• Dualnost po Pontrjaginu