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

Из пројекта Википедија

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

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

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

Садржај

1 Дефиниција
2 Примери
- 2.1 Пример филтера
- 2.2 Пример у C-u

[уреди] Дефиниција

Узмимо да је $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 = {1\over N}\sum_{j=0}^{N-1} w^{-j\cdot k}\cdot \hat c_j$ за $k=0,\dots,N-1$

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

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

Онда је ДФТ за $c\in\Bbb C^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\Bbb 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$ за $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πi = 1$ . Такође важи $\overline{e^{\mathrm{i}\phi}}=e^{-\mathrm{i}\phi}$ и $\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\Bbb C^N$ није реалан, али само N независних вредности (уместо 2N).

За ИДФТ можемо закључити следеће: Уколико за $\hat c\in\Bbb C^N$ важи $\hat c_{N-k}=\overline{\hat c_k}$ за све $k=0,\dots,N-1$ , онда је ИДФТ реалан вектор $c\in\R^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,..., N-K$ , $K$ је негде у близини $\frac{N}{2}$ .

Узмимо неку функцију $f$ којој додељујемо $x\in\Bbb C^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=\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)$

[уреди] Примери

[уреди] Пример филтера

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

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

Наш циљ је да избацимо све фреквенције које су "тише" (тј. које имају амплитуду) од 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} = 10 Hz$ . Стога ми можемо да реконструишемо талас до 5 Hz. Уколико у нашем оригиналном сигналу има фреквенција виших од 5 Hz онда ће наша реконструкција имати грешку. Али, као и увек у животу, човек мора бити оптимистичан те ћемо ми претпоставити да нема виших фреквенција (то је уосталом и један од разлога зашто компакт диск има фреквенцију од око 41 kHz; људско ухо може да региструје тек до 20 kHz!).

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

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

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

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

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

$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 )}$ . У нашем случају наш фреквентни спектрум изгледа овако:

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

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

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

[уреди] Пример у 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;
}

Добављено из "http://sr.wikipedia.org../../../%D0%B4/%D0%B8/%D1%81/%D0%94%D0%B8%D1%81%D0%BA%D1%80%D0%B5%D1%82%D0%BD%D0%B0_%D0%A4%D1%83%D1%80%D0%B8%D1%98%D0%B5%D0%BE%D0%B2%D0%B0_%D1%82%D1%80%D0%B0%D0%BD%D1%81%D1%84%D0%BE%D1%80%D0%BC%D0%B0%D1%86%D0%B8%D1%98%D0%B0_2e0e.html"

Категорије страница: Дискретна математика · Математичка анализа