Дискретна Фуријеова трансформација

Дискретна Фуријеова трансформација или ДФТ јесте Фуријеова трансформација дискретног и коначног (или периодичног) сигнала. Дискретна Фуријеова трансформација је тиме и специјални облик Z-трансформације код које се z налази на јединичном кругу. Често се користи при обради дигиталних сигнала, а најпознатији алгоритам за то је брза фуријеова трансформација (FFT, Fast Fourier Transformation, енгл.).

Дискретна Фуријеова трансформација може да послужи такође за апроксимацију (у одређеним случајевима и реконструкцију) функције која одговара сигналу или као имплементација дигиталних филтера.

Путем инверзне Фуријеове трансформације се из Фуријеових коефицијената склапа излазни сигнал, а повезивањем ДФТ и инверзне ДФТ можемо да манипулишемо фреквенцијама (налази примену при еквилајзерима и филтерима).

Дефиниција

Узмимо да је $R$ комутативан, унитаран прстен, у којем је број $N$ јединица. Даље, у $R$ је $w$ јединични корен.

За вектор $c=(c_0,\dots ,c_{N-1})\in R^N$ је дискретна фуријеова трансформација $\hat{c}$ на следећи начин дефинисана:

${\hat{c}}_k=\sum _{j=0}^{N-1}{w}^{j\cdot k}\cdot c_j$ за $k=0,\dots ,N-1$

А за ${\hat{c}}_k$, инверзна фуријеова трансформација је

$c_k=\frac{1}{N}\sum _{j=0}^{N-1}{w}^{-j\cdot k}\cdot {\hat{c}}_j$ за $k=0,\dots ,N-1$

ДФТ и ИДФТ у комплексном домену

У комплексном домену користимо $w=e^{\frac{2\ast \pi \ast k}{n}},k=0,1,\dots ,n-1$.

Онда је ДФТ за $c\in {ℂ}^N$: ${\hat{c}}_k=\sum _{j=0}^{N-1}{e}^{-2\pi \mathrm{i}\cdot \frac{jk}{N}}\cdot c_j$ за $k=0,\dots ,N-1$,

а ИДФТ за $\hat{c}\in {ℂ}^N$: $c_k=\frac{1}{N}\sum _{j=0}^{N-1}{e}^{2\pi \mathrm{i}\cdot \frac{jk}{N}}\cdot {\hat{c}}_j$ за $k=0,\dots ,N-1$

ДФТ и ИДФТ у реалном домену

Рачуница у реалном домену је:

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

Ојлеров идентитет гласи: $e^{2\pi \mathrm{i}}=1$. Такође важи $\overline{e^{\mathrm{i}\varphi }}=e^{-\mathrm{i}\varphi }$ и $\overline{z_1+z_2}=\overline{z_1}+\overline{z_2}$.

Стога можемо још упростити израз:

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

Што ће рећи, $\hat{c}\in {ℂ}^N$ није реалан, али само N независних вредности (уместо 2N).

За ИДФТ можемо закључити следеће: Уколико за $\hat{c}\in {ℂ}^N$ важи ${\hat{c}}_{N-k}=\overline{{\hat{c}}_k}$ за све $k=0,\dots ,N-1$, онда је ИДФТ реалан вектор $c\in {ℝ}^N$.

Померање и скалирање у времену и фреквенцији

Ако је сигнал периодичан, онда није битно да ли трансформишемо у опсегу $0,\dots ,N-1$ или $k,\dots ,N-1+k$. Индексна променљива j треба да обухвати N опсег, али није битно где он почиње односно где се завршава (ово важи само за случај да је сигнал периодичан, тј. да се вектор $k,\dots ,N-1+k$ периодично понавља). Присетимо се: за $w$ важи $w^N=1$. Онда $w^{N+k}=w^N\cdot w^k=1\cdot w^k=w^k$.

У пракси често желимо да разлика у индексима буде истовремено и разлика у времену или раздаљини два мерења

$t_n=nT,n=1-M,\dots ,N-M$, $T$ је периода нашег мерења.

Често желимо и да коефицијентима доделимо фреквенцију тако да су центриране око $0$

${\omega }_n:=2\pi \frac{n}{NT},n=1-K,\mathrm{...},N-K$, $K$ је негде у близини $\frac{N}{2}$.

Узмимо неку функцију $f$ којој додељујемо $x\in {ℂ}^N$ тако да $x_n=f(t_n)$.

ДФТ је онда $y_n=\hat{f}({\omega }_n)=F({\omega }_n)$.

Из тога следи:

$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)$

а ИДФТ је

$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)$

Примери

Пример филтера

Ситуација: Звук који желимо да снимимо има следећи облик (када би га снимао аналоган микрофон): Датотека:Dft signal originalan.gif

Пошто је наш микрофон дигиталан, ми можемо само да снимимо појединачне вредности. На нашем компјутеру добијамо: Датотека:Dft signal diskretan.gif

Наш циљ је да избацимо све фреквенције које су „тише“ (тј. које имају амплитуду) од 1 V. Прво правимо табелу:

<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

Имамо 10 вредности на 1 секунду, што ће рећи периода нашег мерења је $T=0.1$, а фреквенција $freq=\frac{1}{T}=10Hz$. Стога ми можемо да реконструишемо талас до 5 -{Hz}-. Уколико у нашем оригиналном сигналу има фреквенција виших од 5 -{Hz}- онда ће наша реконструкција имати грешку. Али, као и увек у животу, човек мора бити оптимистичан те ћемо ми претпоставити да нема виших фреквенција (то је уосталом и један од разлога зашто компакт-диск има фреквенцију од око 41 kHz; људско ухо може да региструје тек до 20 kHz!).

Следи израчунавање ${\omega }_i$. Нас занимају само вредности везане за позитивне индексе:

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

$n=-4,\dots ,5$

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

${\omega }_0=2\pi \frac{0}{NT}=0,{\omega }_1=2\pi ,{\omega }_2=4\pi ,\dots ,{\omega }_5=10\pi$

Сада имамо све вредности и можемо да почнемо са рачунањем:

$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)=\dots =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)=\dots =15+0\mathrm{i}=y_2={\hat{c}}_2$

Израчунавање осталих коефицијената иде аналогно, те ћемо их овде само навести као резултате:

$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$

Имамо $\hat{c}$, сада желимо да избацимо све превише „тихе“ тонове. Требају нам $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

Знамо да важи: $c_0=a_0,c_{i>0}=\frac{1}{2}(a_i-b_i\mathrm{i})$. На тај начин можемо да израчунамо $a_i$ и $b_i$:

$a_0=10$

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

Остале амплитуде:

$a_2=3$

$a_3=0.5$

$b_3=0.6$

$a_4=b_4=0$

Из $a_4=b_4=0$ можемо да закључимо да фреквенција од 4 -{Hz}- нема у нашем сигналу. Често је врло згодно навести све амплитуде у графикону. Амплитуда $A_k$ за неку фреквенцију $k$ је $A_k=\sqrt{(a_k^2+b_k^2)}$. У нашем случају наш фреквентни спектрум изгледа овако:

Датотека:Dft signal amplitude.gif

Све $a_i$ и $b_i$ за које важи $\sqrt{(a_i^2+b_i^2)}<1$ избацујемо и на крају добијамо реконструисану и обрађену функцију:

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

Датотека:Dft signal obradjen.gif

Сада можемо да поново да израчунамо $f(t_i)$ или да се послужимо ИДФТ и тако прерађен сигнал снимимо у меморију.

Пример у -{C}--у

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

Референце

• Шаблон:Cite book
• Шаблон:Cite book
• Шаблон:Cite book esp. section 30.2: The DFT and FFT, pp. 830–838.
• Шаблон:Cite book
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal
• Шаблон:Cite journal

Види још

• Дуалност по Понтрјагину