1. Strona główna
  2. Matematyka Obliczeniowa - zadania
  3. Algorytm FFT

Algorytm FFT

W tym rozdziale przedstawimy zadania laboratoryjne, w których przetestujemy
dyskretną transformatę Fouriera (DFT) oraz szybki algorytm jej obliczania, czyli algorytm szybkiej trasformaty Fouriera; FFT (ang. Fast Fourier Transform)

W opisie matematycznym przyjmujemy, że mamy do czynienia z wektorami okresowymi o okresie długości $ N $ indeksowanymi od zera i
$ x_k=x_{k+N} $ dla $ k $ ujemnych.

Przez $ \mathcal{F}_N x $ oznaczymy wynik działania operatora DFT na wektorze $ x $, a przez $ \mathcal{F}_N^{-1} x $ oznaczymy wynik działania operatora odwrotnego do DFT na wektorze $ x $.

  1. Funkcja octave'a z algorytmem FFT fft(), która oblicza wartość DFT i funkcja octave'a z algorytmem odwrotnym FFT ifft(), która oblicza wartość transformaty odwrotnej do DFT.

    Sprawdź, czy rzeczywiście operacje wykonywane przez fft() i ifft() są do siebie odwrotne, tzn. dla różnych losowych wektorów
    $ x $ o długości $ N=4,8,16,32 $ policz $ y= \mathcal{F}_N x $ za pomocą
    fft(), a następnie $ z= \mathcal{F}_N^{-1} y $ za pomocą
    ifft().

    Policz normy drugie różnicy $ z-x $.

  2. Skalowanie w funkcjach octave'a fft() i ifft().

    W literaturze zdarza się definicja DFT bez $ 1/N $ przed macierzą.
    Wtedy macierz odwrotną do DFT należy przemnożyć przez $ N $ albo obie przez $ 1/\sqrt{N} $.

    Sprawdź, jaka jest stała przed DFT obliczanym przez funkcję
    fft(), tzn. wyznacz $ C_N $ takie, że

    $$(\mathcal{F}_N x)_j=C_N \sum_{k=0}^{N-1}x_k\exp(-2\pi*j*k/N)$$

    dla $ \mathcal{F}_Nx $ obliczanego przez tę funkcję octave'a .

  3. Sprawdź, czy DFT obliczona przy pomocy funkcji fft() spełnia

    $$\|\mathcal{F}_N x\|_2= C_N\|x\|_2$$

    ze stałą $ C_N $ dla stu losowych wektorów $ x $, oraz czy analogicznie odwrotna DFT obliczona ifft() spełnia:

    $$\|\mathcal{F}_N^{-1} y\|_2= C_N^{-1}\|x\|_2.$$
  4. Algorytm szybkiego mnożenia wielomianów.

    Napisz funkcję octave'a

    function w=polymult(p,q)

    szybko mnożącą dwa wielomiany stopnia nie większego od $ N $, dla których znamy ich współczynniki w bazie potęgowej. Funkcja powinna korzystać z DFT i odwrotnej DFT, czyli dyskretnej transformaty Fouriera i odwrotnej dyskretnej transformaty Fouriera.

    Parametry funkcji powinny być:

    • wektor $ p=(p_k)_k $ ze współczynnikami wielomianu $ P(x)=\sum_kp_kx^k $,
    • wektor $ q=(q_k)_k $ ze współczynnikami wielomianu $ Q(x)=\sum_kq_kx^k $.

    Funkcja powinna zwrócić
    współczynniki wielomianu $ W(x)=P(x)*Q(x)=\sum_kw_kx^k $.

    Sprawdź działanie funkcji dla wielomianów $ P(x)=1+x $ i $ Q(x)=2+x $, czyli $ (P*Q)(x)=2+3x+x^2 $, a potem dla wielomianów dużych stopni, czyli np. dla $ Q(x)=x^{50} $ i $ P(X)=x^{50}+1 $ funkcja powinna zwrócić współczynniki wielomianu $ (P*Q)(x)=x^{100}+x^{50} $.

    Porównaj czas obliczeń tej funkcji z wynikiem funkcji octave'a conv() dla wielomianów różnych stopni z losowymi współczynnikami.

  5. DFT a dyskretny splot dwóch wektorów.

    Napisz dwie funkcje
    obliczające splot dwóch wektorów $ z=x\star y $

    • wprost z definicji, tzn. $ z_k=\sum_{j=0}^{N-1}x_ky_{k-j} $
    • z wykorzystaniem DFT - tu musimy skorzystać z tego, że
      $$<br /> (\mathcal{F}_N(x\star y))_k=\alpha_N(\mathcal{F}_N x)_k * (\mathcal{F}_N y)_k<br /> \qquad k=0,\ldots,N-1<br /> $$

      dla pewnej stałej $ \alpha_N $.
      Trzeba tę stałą wyznaczyć teoretycznie albo eksperymentalnie.

      Przetestuj obie funkcje dla losowych wektorów i różnych $ N=4,8,64,1024 $.
      Porównaj z wynikami funkcji octave'a conv(). Należy zapoznać się z pomocą do tej funkcji.

  6. Filtry a DFT.

    W niektórych zastosowaniach współczynniki wektora $ x $ mogą odpowiadać próbkom sygnału (np. dźwięku). Stosuje się wtedy filtry w postaci splotu $ x $ z zadanym wektorem $ F $ pełniącym rolę filtru:
    $ x\star F $.

    Rozważmy filtr postaci $ F=\frac{1}{4}*(1,1,1,1,0,\ldots,0)^T $ (pierwsze cztery współrzędne są równe $ \frac{1}{4} $, pozostałe - zero).

    Zastosuj ten filtr do wektora imitującego próbki sygnału z szumem: $ z=x+y $,
    gdzie $ x_k=\sin(k*h) $ dla $ h=2*\pi/N $ - imituje sygnał bez szumów, a
    $ y $ jest wektorem losowym o wartościach w $ [-0.1,0.1] $ imitującym losowe zaburzenie, czyli tzw. szumy. Przetestuj filtr dla dużych $ N $, np. $ N=1024 $ lub $ 2048 $, tzn.:

    • Narysuj wykres $ z $ przed zastosowaniem filtra i po zastosowaniu, tzn.
      wykresy $ z $ i $ z\star F $. Można narysować mały fragment wykresu, np dla $ 100\leq k \leq 120 $. Jaki efekt optyczny widać na wykresach?

    • Policz maksimum z modułu różnicy kolejnych elementów obu wektorów, tzn. $ \max_{k>0}|z_k-z_{k-1}| $ i
      $ \max_{k>0}|(z\star F)_k-(z\star F)_{k-1}| $. Im ta wartość dla danego wektora z próbkami sygnału jest mniejsza, tym sygnał jest gładszy, czyli możemy uznać go za sygnał z mniejszą ilością szumów.