\documentclass[10pt]{article}
\usepackage[T1]{fontenc}
\usepackage[latin2]{inputenc}
\usepackage{polski}
\usepackage{indentfirst}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{warn}
\usepackage[A4]{vmargin}
\usepackage{here}

\usepackage{listings}
\lstset{
  language=C,
  indent=0pt,
  extendedchars=true,
  commentstyle=\itshape,
  stringstyle=\ttfamily,
  stringspaces=false,
  mathescape,
  }

\pagestyle{plain}

\newcommand{\cc}[1]{{#1}^{\ast}}
\newcommand{\abs}[1]{\lvert#1\rvert}
\newcommand{\field}[1]{\mathbb{#1}}

\begin{document}
\title{
  Szybka transformata Fouriera (FFT)\\
  i~jej zastosowania\\
  \large
  Metody Numeryczne II A}
\date{19.04.2001 r.}

\maketitle

%%%
\lstdefinestyle{number}
  {labelstep=1, labelstyle=\tiny, labelsep=5pt}
\lstdefinestyle{nonumber}
  {labelstep=0}

%%%
\section{Szybka transformata Fouriera (FFT)}
\label{sec:wstep}


\subsection{Właściwości transformaty Fouriera}
\label{sec:wlasciwosci-FT}

Fizyczny proces (sygnał) może być opisany albo w~\emph{dziedzinie
  czasowej}, jako (być może zespolona) funkcja czasu $h(t)$, albo
w~\emph{dziedzinie częstości}, jako zespolona (by opisać fazę) funkcja
częstości $H(f)$.  Przejścia między tymi dwoma reprezentacjami można
dokonać za pomocą transformaty Fouriera i~odwrotnej transformaty
Fouriera:
\begin{subequations}
  \label{eq:transformaty-Fouriera}
  \begin{align}
    \label{eq:FT}
    H(f) &= \int_{-\infty}^{\infty} h(t) e^{2 \pi \imath f t} dt\\
    \label{eq:RFT}
    h(t) &= \int_{-\infty}^{\infty} H(f) e^{-2 \pi \imath f t} df\\
  \end{align}
\end{subequations}

Z~równań~\eqref{eq:FT} i~\eqref{eq:RFT} widać, że transformacja
Fouriera jest operacją liniową, t.j. transformata Fouriera kombinacji
liniowej \(\alpha f + \beta h = \alpha F + \beta H\), gdzie $F$ jest
transformatą Fouriera funcji $f$, a~$H$ -- funkcji $h$, zaś \(\alpha,
\beta \in \field{C}\).

Funkcja $h(t)$ może mieć pewne szczególne właściwości: może być czysto
rzeczywista lub czysto urojona, może być parzysta bądź nieparzysta
\emph{Uwaga:} w~przypadku dyskretnej skończonej transformaty Fouriera
(patrz rozdział~\ref{sec:DFT}) badana funkcja jest funkcją okresową
o~okresie $N \Delta$, tak więc $h(-x_k) = h(x_{N-k})$.  Właściwości
funkcji $h(t)$ przekładają się na odpowiednie właściwości jej
transformaty Fouriera $H(f)$ --- patrz tabela~\ref{tab:symetrie-FT}.

\begin{table}[htb!]
 \centering
 \begin{tabular}[t]{ll}
   \multicolumn{1}{c}{funkcja $h(t)$} &
   \multicolumn{1}{c}{transformata $H(f)$} \\
    \hline\raisebox{0 ex}[3 ex]{\null}
    $h(t)  = \cc{h(t)}$  rzeczywista &
    $H(f)  = \cc{H(-f)}$  \\
    $h(t)  = -\cc{h(t)}$ urojona &
    $H(f)  = -\cc{H(f)}$ \\
    $h(t)  = h(-t)$ parzysta &
    $H(f)  = H(-f)$ parzysta \\
    $h(t)  = -h(-t)$ nieparzysta &
    $H(f)  = -H(-f)$ nieparzysta \\
    $h(t)$ rzeczywista i~parzysta &
    $H(f)$ rzeczywista i~parzysta \\
    $h(t)$ rzeczywista i~nieparzysta &
    $H(f)$ urojona i~nieparzysta \\
    $h(t)$ urojona i~parzysta &
    $H(f)$ urojona i~parzysta \\
    $h(t)$ urojona i~nieparzysta &
    $H(f)$ rzeczywista i~nieparzysta 
  \end{tabular}
  \caption{Symetrie $h(t)$ a~symetrie $H(f)$}
  \label{tab:symetrie-FT}
\end{table}

Zachodzą ponadto następujące własności
\begin{subequations}
  \label{eq:wlasnosci-h-vs-H}
  \begin{align}
    h(at) &\Longleftrightarrow \frac{1}{\abs{a}} H(\frac{f}{a})
    &&\text{skalowanie czasu}
    \label{eq:FT-skalowanie-czasu}\\
    \frac{1}{\abs{b}} h(\frac{t}{b}) &\Longleftrightarrow H(bf)
    &&\text{skalowanie częstości}
    \label{eq:FT-skalowanie-czestosci}\\
    h(t - t_0) &\Longleftrightarrow H(f) e^{2 \pi \imath f t_0}
    &&\text{przesunięcie w~czasie}
    \label{eq:FT-przesuniecie-w-czasie}\\
    h(t) e^{-2 \pi \imath f_0 t} &\Longleftrightarrow H(f - f_0)
    &&\text{przesunięcie w~częstościach}
    \label{eq:FT-przesuniecie-w-czestosciach}\\
    \frac{d}{dt} h(t) &\Longleftrightarrow -2 \pi \imath f H(f)
    &&\text{różniczkowanie po czasie}
    \label{eq:FT-rozniczkowanie-po-czasie}
  \end{align}
\end{subequations}

\subsection{Dyskretna transformata Fouriera}
\label{sec:DFT}

W~najczęściej spotykanych sytuacjach (zwłaszcza przy obróbce sygnału)
funkcja $h(t)$ jest zadana przez wartości w~równoodległych chwilach
czasu:
\begin{equation}
  \label{eq:hk-infty}
  h_k = h(k \Delta) \quad n = \dots,-3,-2,-1,0,1,2,3,\dots
\end{equation}
Dla każdej wielkości przedziału $\Delta$ istnieje szczególna częstość
$f_c$ zwana \textbf{krytyczną częstością Nyquista}, dana wzorem
\begin{equation}
  \label{eq:Nyquist}
  f_c \equiv \frac{1}{2 \Delta}
\end{equation}
Funkcja o~częstości większej od częstości Nyquista ma te same wartości
w~punktach pomiaru co pewna funkcja o~częstości mniejszej od częstości
Nyquista.  To zjawisko nazywa się \textbf{aliasingiem}.  Aby zatem nie
dostać błędnych wyników częstość próbkowania musi być dwa razy większa
od najwyższej znaczącej częstości w~sygnale.  Oznacza to także, że dla
funkcji zadanej przez wartości w~dyskretnych punktach możemy określić
jednoznacznie transformatę Fouriera tylko dla $-f_c < f < f_c$.

W~obliczeniach komputerowych będziemy wyliczać transformatę Fouriera
na podstawie wielkości sygnału w~skończonej ilości równoodległych
punktów próbkowania.  Załóżmy, że mamy $N$ (dla uproszczenia niech $N$
będzie parzyste) takich punktów
\begin{equation}
  \label{eq:hk}
  h_k \equiv h(t_k), \quad t \equiv k \Delta,
  \quad k = 0,1,2,\dots,N-1
\end{equation}
Możemy wówczas określić (jednoznacznie) transformatę Fouriera tylko
w~$N$ punktach:
\begin{equation}
  \label{eq:fn}
  f_n \equiv \frac{n}{N \Delta}, 
  \quad n = -\frac{N}{2},\dots,\frac{N}{2}
\end{equation}
(uwaga: $f_{-N/2} = f_{N/2}$).

Transformata funkcji zadanej na skończonej liczbie punktów równej $N$
jest identyczna z~dyskretną transformatą Fouriera funkcji okresowej
o~okresie $N \Delta$ i~wartościach dla jednego okresu takich jak
zadana funkcja.  \emph{Uwaga:} zatem dla dyskretnej transformaty
Fouriera funkcja parzysta to taka, że \(h_k = h_{N-k} = h_{-k}\).

Przy obliczaniu transformaty Fouriera funkcji należy pamiętać, że
wartości ujemne argumentu odpowiadają odpowiednim indekson \(N-k\)
z~końca tablicy.  Jeśli zatem liczymy transformatę Fouriera funkcji
skoncentrowanej wokół zera (np. funkcja Gaussa), to aby uniknąć
mnożenia transformaty przez $e^{2 \pi \imath f \text{przesunięcie}}$ (patrz
wzór~\eqref{eq:FT-przesuniecie-w-czasie}) należy rozmieścić wartości
funkcji na ,,zawiniętej'' tablicy, t.j. \lstinline!h[k]! \(= h(k
\Delta)\) dla \(k < N/2\), zaś \lstinline!h[N-k]! \(= h(-k \Delta)\)
dla \(k < N/2\).

Dyskretna transformata Fouriera i~odwrotna dyskretna transformata
Fouriera zadane są wzorami
\begin{subequations}
  \label{eq:dyskretne-transformaty-Fouriera}
  \begin{align}
  \label{eq:DFT}
  H_n &= \sum_{k=0}^{N-1} h_k e^{2 \pi \imath k n/N} \\
  \label{eq:DRFT}
  h_k &= \frac{1}{N} \sum_{n=0}^{N-1} H_n e^{-2 \pi \imath k n/N} 
  \end{align}
\end{subequations}
Skorzystaliśmy tutaj z~tego, że (jak łatwo zauważyć) $H_{-n} =
H_{N-n}$, tak więc zamiast zadawać $H_n$ dla $n = -N/2,\dots,N/2$
możemy używać $n = 0,\dots,N-1$.  \emph{Uwaga:} różne procedury
(biblioteki) stosują różną konwencję znaków w~transformacie
i~odwrotnej transformacie Fouriera.  Podajemy tutaj konwencję za
,,Numerical Recipies''~\cite{NRC}.

Związek między dyskretną a~zwykłą transformatą Fouriera zadany jest
wzorem
\begin{equation}
  \label{eq:FT-vs-DFT}
  H(f_n) \approx \Delta H_n
\end{equation}



\subsection{Algorytm szybkiej transformaty Fouriera}
\label{sec:algorytm-FFT}

Najprostsza wersja algorytmu obliczania dyskretnej transformaty
Fouriera (wprost ze wzorów~\eqref{eq:dyskretne-transformaty-Fouriera})
wymaga ilości obliczeń proporcjonalnej do $N^2$ (gdzie $N$ jest liczbą
punktów, w~których zadana jest nasza funkcja).  Algorytmy szybkiej
transformaty Fouriera korzystają z~rozkładu $N$ na czynniki: $N = r_1
r_2 \dots r_p$ i~mają koszt (ilość obliczeń) proporcjonalną do $N(r_1
+ r_2 + \dots + r_p$.  Najprostszą wersję algorytmu uzyskujemy dla
ilości punktów będącej potęga dwójki: $N = 2^p$.

Algorytm obliczania szybkiej transformaty Fouriera dla $N = 2^k$ nosi
nazwę \textbf{metody Cooleya--Tukeya}.  Opera się ona na zastosowaniu
tzw. metody \emph{dziel i~rządź}: dzielimy problem na dwa \emph{takie
  same} problemy o~dwukrotnie mniejszym rozmiarze, a~następnie
składamy je do rozwiązania pełnego problemu.  

Dyskretna transformata Fouriera dla danych długości $N$ może być
zapisana jako kombinacja liniowa dwu transformat Fouriera, każdej
o~długości $N/2$; jednej złożonej z~punktów o~parzystych, a~drugiej
z~punktów o~nieparzystych indeksach w~oryginalnej transformacie:
\begin{equation}
  \label{eq:wyprowadzenie-FFT}
  \begin{split}
    H_n 
    &= \sum_{k=0}^{N-1} e^{2 \pi \imath k n/N} h_k = \\
    &= \sum_{\beta=0}^{N/2-1} e^{2 \pi \imath (2\beta) n/N} h_{2\beta} +
       \sum_{\beta=0}^{N/2-1} e^{2 \pi \imath (2\beta+1) n/N} h_{2\beta+1} = \\
    &= \sum_{\beta=0}^{N/2-1} e^{2 \pi \imath \beta n/(N/2)} h_{2\beta} +
       \left(e^{2 \pi \imath/N}\right)^n
       \sum_{\beta=0}^{N/2-1} e^{2 \pi \imath \beta n/(N/2)} h_{2\beta+1} = \\
    &=     H^{ \shortmid}_{n \bmod \frac{N}{2}} 
     + W^n H^{\nshortmid}_{n \bmod \frac{N}{2}} 
  \end{split}
\end{equation}
Wprowadziliśmy oznaczenie \(W = e^{2 \pi \imath/N}\).  Łatwo zauważyć, że
w~obu sumach wartość dla \(n = N/2 + n_1\) (gdzie $k \geq 0$) jest taka
sama jak wartość dla \(n_1\):
\begin{equation*}
  (W^2)^{n \beta} = (W^2)^{\bigl(\frac{N}{2} + k_1\bigr)\beta}
  = (W^2)^{\frac{N\beta}{2}} (W^2)^{k_1 \beta} 
  = (W^N)^{\beta} (W^2)^{k_1 \beta} = (W^2)^{k_1 \beta},
\end{equation*}
gdyż $W^N = 1$.

Dla każdej z~części możemy zastosować to samo rozumowanie, za każdym
razem zmniejszając rozmiar zagadnienia dwa razy.  Po $\log_2 N$
krokach dojdziemy do transformaty Fouriera o~długości 1.  Transformata
Fouriera dla $N = 1$ jest funkcją identycznościową.  Mamy zatem
\begin{equation*}
  H^{\shortmid\nshortmid\nshortmid\shortmid\dots\shortmid}_n 
  = h_l \quad \text{dla pewnego $l$}
\end{equation*}
(symbol $\mid$ odpowiada braniu częsci złożonej z~parzystych indeksów
w~kolejnym kroku ($2 \mid k$), symbol $\nmid$ z~nieparzystych ($2
\nmid k$)).

Pozostaje zatem znaleźć indeks $l$ odpowiadający danej sekwencji
wyborów indeksów parzystych/nieparzystych.  Widzimy, że przy pierwszym
kroku w~sumie oznaczonej symbolem $\mid$ pojawiają się tylko elementy
$h_k$ o~indeksach parzystych, czyli ostatni bit indeksu $l$ dla
sekwencji rozpoczynającej się od $\mid$ wynosi 0.  Następny podział
zachodzi ze względu na przedostatni bit indeksu tablicy $h$.  Zatem
indeks $l$ odpowiadający sekwencji symboli
$\mid\nmid\nmid\mid\dots\mid$ (określającej kolejne wybory indeksów
parzystych/nieparzystych w~kolejnej części) otrzymujemy, biorąc liczbę
binarną uzyskaną przez odwrócenie sekwencji i~zamianę $\mid$ na 0,
a~$\nmid$ na 1.

Aby uniknąć kłopotliwego wyliczania indeksów po których sumujemy
w~kolejnych podziałach (lub ich zapamiętywanie w~dodatkowej tablicy),
najprościej jest ułożyć (przestawić) elementy tablicy z~danymi (naszym
sygnałem) \lstinline!h[k]! (czyli $h_k$) tak, aby sumowanie rozciągało
się zawsze po kolejnych elementach tablicy.  Przykładowa procedura
rekurencyjna (używająca dadatkowego pomocniczego wektora) mogła by
wyglądać następująco:
\begin{lstlisting}{FFT-przyklad}
void fft_rec(fcomplex data[], float tmp[], int n, int isign)
{
  int i, imod; 
  int n2 = n/2; 
  fcomplex wn = { 0, 0 };

  if(n == 1) /* przypadek graniczny */
    return;

  /* przestawiamy elementy, by wywołania rekurencyjne */
  /* miały je ułożone po kolei */
  for (i = 0; i < n2; i++) {
    /* elementy o indeksach parzystych */
    tmp[i].re = data[2*i].re;
    tmp[i].im = data[2*i].im;

    /* elementy o indeksach nieparzystych */
    tmp[i + n2].re = data[2*i+1].re;
    tmp[i + n2].im = data[2*i+1].im;
  }
  
  /* wywołujemy rekurencyjnie funkcję dla */
  /* "indeksów parzystych" i "nieparzystych */
  fft_rec(tmp,    data, n2, isign);
  fft_rec(tmp+n2, data, n2, isign);

  /* wyliczamy kolejne elementy transformaty fouriera */
  for (i = 0; i < n; i++) {
    imod = i % n2;
    wn.re = cos(isign*2*M_PI*i/N);
    wn.im = sin(isign*2*M_PI*i/N);

    /* $H_n =  H^{o}_{n \bmod N/2} + W^n H^{e}_{n \bmod N/2}$ */ 
    data[i].re = tmp[imod].re + wn.re*tmp[n2 + imod].re - wn.re*tmp[n2 + imod].re; 
    data[i].re = tmp[imod].im + wn.re*tmp[n2 + imod].im + wn.im*tmp[n2 + imod].re;
  }
} 
\end{lstlisting}
(Jawnie zapisujemy dodatkowy wektor pomocniczy dla uproszczenia
[zapisu] procedury.)

Stosunkowo łatwo zauważyć, że metoda ta prowadzić będzie do
przestawiania niektórych elementów aż $\log N$ razy (tyle ile
głębokość rekurencji).  Gdyby udało się tego uniknąć i~ustawić
elementy od razu w~żądanej kolejności to rekurencję można by zastąpić
iteracją.  Składamy rozwiązanie dla transformaty Fouriera o~długości 2
z~transformat Fouriera o~długości 1:
\begin{equation*}
  H^{\dots}_k = H^{\dots\mid}_k + W^k H^{\dots\nmid},  
\end{equation*}
następnie dla długości 2 i~tak dalej; prosta iteracja.  Jeśli ułożymy
elementy tablicy $h$ w~ten sposób, by elementy odpowiadające
$H^{\dots\mid}$ i~$H^{\dots\nmid}$ leżały obok siebie, obliczanie
transformaty Fouriera będzie proste.  Należy zatem najpierw posortować
(ułożyć) elementy w~kolejności składania transformat, czyli
w~kolejności odwróconych bitów.  Ponieważ relacja odwrócenia bitowego
jest przemienna, zatem przy odpowiedniej organizacji programu
wystarczy zrobić co najwyżej jedno przemieszczenie na element.  Można
to zrobić ,,w miejscu''.

\medskip

Algorytm można zatem złożyć z~dwu części.  Pierwsza to takie
uporządkowanie elementów, by $h_i$ odpowiadające elementom z~dwu
najbardziej wewnętrznych podziałów leżały obok siebie itd.,
tzn. według kolejności $l$-ów.  Druga to kolejne składania cząstkowych
transformat w~jedną całość.

Bardziej szczegółowe informacje n.t. algorytmu Cooleya--Tukeya można
znaleźć w~,,Numerical Recipies in C''~\cite{NRC}.


\subsection{Opis procedury \texttt{four1} z NRC}
\label{sec:opis-four1}

Procedura \texttt{four1} (zamieszczona na stronie WWW) oblicza
dyskretną transformatę Fouriera tablicy \lstinline!data[1..2*nn]!,
zastępując nią oryginalne dane, dla \lstinline!isign! = 1.  Dla
\lstinline!isign! = -1 oblicza odwrotną transformatę Fouriera
(pomnożoną przez \lstinline!nn!).  Argument \lstinline!nn! oznacza
ilość liczb zespolonych.  Procedurę tę można zastosować bez zmian do
tablic indeksowanych od 0 (jak w~C), tzn. dla
\lstinline!data[0..2*nn-1]! wywołując ją w~następujący sposób
\begin{lstlisting}{}
  four1(data - 1, N, 1);
\end{lstlisting}
Dane rozmieszczone są jako kolejno część rzeczywista i~urojona
kolejnych liczb zespolonych.  \lstinline!data[1]! jest częścią
rzeczywistą pierwszej wartości $h$, zaś \lstinline!data[2]! jest jej
częścią urojoną.  Na macierz \lstinline!data! można patrzeć albo jako
na tablicę \lstinline!nn! liczb zespolonych \lstinline!fcomplex!, albo
jako na dwa razy wiekszą tablicę liczb rzeczywistych.

Uwaga: W~macierzy wynikowej w~elementach \lstinline!data[1]!
i~\lstinline!data[2]! znajduje się częstość zerowa, zaś w~elementach
\lstinline!data[2*nn-1]! i~\lstinline!data[2*nn]! częstość Nyquista.
Rozmiar tablicy (liczony w~liczbach zespolonych) \textbf{musi} być
całkowitą potęgą 2.


\subsection{Analiza procedury \texttt{four1} z NRC}
\label{sec:analiza-four1}

\begin{lstlisting}[style=number]{four1}
include <math.h>

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
\end{lstlisting}
Powyższa makrodefinicja zamienia dwie zmienne wartościami.  Uwaga: nie
da się jej użyć w~wyrażeniu, musi być zakończona średnikiem.
Podwójnie szacuje swoje argumenty (czyli nie należy jej np. wywoływać
z~argumentami \lstinline!SWAP(*src++, *dest++)!). 
\begin{lstlisting}[style=number]{four1}

void four1(float data[], int nn,int isign)
\end{lstlisting}
Zastępuje \lstinline!data[1..2*nn]! przez jej dyskretną transformatę
Fouriera, jeśli \lstinline!isign! jest równe $1$; lub zastępuje
\lstinline!data[1..2*nn]! przez \lstinline!nn! razy jej odwrotna
transformata Fouriera (nieznormalizowaną odwrotną transformatę
Fouriera), jeśli \lstinline!isign! jest równe $-1$. \lstinline!data!
jest tablicą liczb zespolonych o~długości \lstinline!nn!, lub patrząc
inaczej tablicą liczb rzeczywistych o~długości
\lstinline!2*nn!. Rozmiar tablicy \lstinline!nn!
\textbf{\uppercase{musi}} być całkowitą potęgą 2 (procedura tego nie
sprawdza!).
\begin{lstlisting}[style=number]{four1}
{
  unsigned long n, mmax, m, j, istep, i;
  double wtemp, wr, wpr, wpi, wi, theta; /* podwójna precyzja do
                                            rekurencji trygonometrycznej */ 
  float tempr, tempi;

  n = nn << 1; /* n = 2*nn */
  j = 1;

  /*** Odwracanie bitowe ***/
  for (i = 1; i < n; i += 2) {
    if (j > i) {
      /* zamienamy miejscami dwie liczby zespolone */
      SWAP(data[j], data[i]);
      SWAP(data[j+1], data[i+1]);
    }
    m = n >> 1;
    while (m >= 2 && j > m) {
      j -= m;
      m >>= 1;
    }
    j += m;
  }
\end{lstlisting}
Najpierw sortujemy elementy w~kolejności odwróconych bitów, tak by
kolejne kawałki, które składamy ze sobą zgodnie ze wzorem
Danielsona--Lanczosa~\eqref{eq:wyprowadzenie-FFT}, zajmowały ciągłe
kawałki tablicy i~leżały obok siebie.

Algorytm użyty w~procedurze \lstinline!four1! z~NRC jest trudny do
przeanalizowania.  Prostszy algorytm, chociaż nie tak szybki, mógłby
korzystać z~makrodefincji:
\begin{lstlisting}{makrodefinicje}
#define SET_BIT(Value, Bit)      \
  ((Value) |= ((unsigned long)1 << (Bit)))

#define CLEAR_BIT(Value, Bit)    \
  ((Value) &= ~((unsigned long)1 << (Bit)))

#define IS_BIT_SET(Value, Bit)   \
  (((Value) & ((unsigned long)1 << (Bit))) != 0)

/* poniższa definicja nie jest najlepsza */
#define NBITS(Value)   \
  ((int)(log((double)Value)/log(2.0)))
\end{lstlisting}
Należy pamiętać, że odwracamy bity do \lstinline!NBITS(nn)! wyłącznie
(czyli do \lstinline!NBITS(nn) - 1! $=$ \lstinline!NBITS(n/2)!
włącznie.
\begin{lstlisting}[style=number]{four1}

  /** Algorytm Danielsona-Lanczosa **/
  mmax = 2;
  while (n > mmax) { /* zewnętrzna pętla wykonywana $\log_2$ (*\lstinline!nn!) razy */
    istep = 2*mmax;

    /* inicjujemy rekurencję trygonometryczną */
    theta = 6.28318530717959/(isign*mmax);
    wtemp = sin(0.5*theta);
    wpr = -2.0*wtemp*wtemp;
    wpi = sin(theta);
    wr = 1.0;
    wi = 0.0;

\end{lstlisting}
Korzystamy tutaj ze wzoru na funkcje trygonometryczne kąta
podwojonego: \(\cos 2\theta = 1 - 2\sin^2 \theta\).  \lstinline!wp!
traktowane jako zmienna zespolona jest równe \(wp = (cos \theta - 1) +
i \sin \theta\), zaś \(w = e^{2 \pi \imath (n-1)}\).  Kąt $\theta$ jest tak
dobrany, że \(wp + 1 = e^{2 \pi \imath/N}\).
\begin{lstlisting}[style=number]{four1}
  
    for (m = 1; m < mmax; m += 2) { /* dwie wewnętrzne pętle */
      for (i = m;i <= n;i += istep) {
\end{lstlisting}
Pierwsza pętla, po $m$ przebiega po kolejnych elementach (liczbach
zespolonych) \(H^{\dots\shortmid}_k\) i~\(H^{\dots\nshortmid}\).
Druga pętla przebiega po kolejnych częściach (po kolejnych węzłach
tego samego poziomu w~drzewie rekurencji).  Indeks $i$ odpowiada
$\mid$, zaś $j$ -- $\nmid$.
\begin{lstlisting}[style=number]{four1}
        j = i + mmax;

        /* to jest wzór Danielsona-Lanczosa */
        tempr = wr*data[j]   - wi*data[j+1];
        tempi = wr*data[j+1] + wi*data[j];
        data[j]  = data[i]   - tempr;
        data[j+1]= data[i+1] - tempi;
        data[i]   += tempr;
        data[i+1] += tempi;
\end{lstlisting}
Obliczamy tutaj 
\begin{math}
  H^{\dots}_k = H^{\dots\shortmid}_k + W^k H^{\dots\nshortmid},
\end{math}
gdzie $W$ jest równe \lstinline!wr + $i$ wi!.  Korzystamy z tego, że
\(e^{\alpha i + 2 \pi \imath N/2} = -e^{\alpha i}\). 
\begin{lstlisting}[style=number]{four1}
      }

      /* rekurencja trygonometryczna */
      wr = (wtemp=wr)*wpr - wi*wpi + wr;
      wi = wi*wpr + wtemp*wpi + wi;
\end{lstlisting}
Tutaj wyliczamy nową potęgę $W$ za pomocą odpowiedniej rekurencji
trygonometrycznej \(w_{\text{nowe}} = w e^{2 \pi \imath/N} = w (1 + wp)\).
\begin{lstlisting}[style=number]{four1}
    }
    mmax = istep;
  }
}

#undef SWAP  
\end{lstlisting}

\section{Zastosowania dyskretnej transformaty Fouriera}
\label{sec:zastosowania-DFT}


\subsection{Znajdowanie dominującej częstości}
\label{sec:znajdowanie-omega}

Często jest tak, że w~doświadczeniach dostajemy sygnał w~reprezentacji
,,czasowej'', a~wynik odczytujemy z~dominującej częstości.  Jeśli nie
możemy bezpośrednio mierzyć częstości, jednym z~rozwiązań jest
wprowadzenie sygnału do komputera (za pomocą przetwornika
analogowo\dywiz{}cyfrowego) a~następnie potraktowanie go dyskretną
(zazwyczaj rzeczywistą) transformatą Fouriera.

Przykładem takiego doświadczenia może być pomiar prędkości cząstek
rozpraszających światło w~ośrodku za pomocą efektu Dopplera:
krzyżujemy w~punkcie pomiaru dwie wiązki laserowe z~tego samego lasera
(rozszczepione poprzednio) i~odczytujemy prędkość z~częstości
dudnień.  (Tak naprawdę to pomiar ten wyznacza składową prędkości;
przy użyciu dodatkowych wiązek możemy wyznaczyć cały wektor prędkości.)

Innym przykładem jest badanie zjawiska tzw. rezonansu stochastycznego
w~układach z~detektorami progowymi.  W~najprostszym sformułowaniu mamy
dwa detektory progowe o~progach $A$ i~$-A$.  Na początku sygnał
wyjściowy jest dla uproszczenia zerowy.  Jeśli na wejściu pojawi się
sygnał większy niż $A$ to na wyjściu z~układu pojawia się sygnał $1$,
aż do momentu gdy na wejściu pojawi się sygnał mniejszy niż $-A$ (próg
detektora ,,dolnego'').  Wtedy na wyjściu pojawia się sygnał $-1$, aż
do sygnału na wejściu większego niż próg $A$ itd.  Na wejściu mamy
sygnał okresowy (o~znanej częstości) o~amplitudzie [dużo] mniejszej od
wartości progu detektora, w~towarzystwie szumu.  Jeśli szum jest
dostatecznie duży wśród wartości wejściuwych pojawią się wartości
ponadprogowe.  Wykonujemy analizę fourierowską sygnału wyjściowego
i~badamy kształt widma wokół częstości podprogowego sygnału
okresowego.  Powinniśmy zaobserwować, że wraz ze wzrostem szumu
następuje coraz lepsze wykrywanie częstości sygnału, aż szum zacznie
dominować.  Ponieważ do badania potrzebujemy tylko fragment widma,
zadanie to można wygodnie rozwiązywać za pomocą zwykłego (powolnego)
algorytmu transformaty Fouriera.

\subsection{Filtry pasmowe}
\label{sec:filtry-pasmowe}

Najprostszym zastosowaniem transformaty Fouriera jest filtrowanie
przez wycinanie (bądź tłumienie) nieporządanych częstości.  Możemy
np. za pomocą filtru górnoprzepustowego próbować usunąć szum
z~sygnału, lub za pomocą filtru dolnoprzepustowego wydobyć sygnał
wolnozmienny (trend) z~wysokozmiennego sygnału.

Gdy przesyłamy dane za pomocą modyfikacji sygnału o~danej częstości,
czy to za pomocą modulacji amplitudy, czy to za pomocą modulacji fazy,
możemy poprawić stosunek sygnału do szumu traktując sygnał przed
,,zdekodowaniem'' filtrem pasmowym przenoszącym częstości wokół
częstości nośnej.  Oczywiście nie możemy dać filtru zbyt wąskiego,
inaczej możemy zakłócić przesyłane dane.

Jeśli badamy sygnał wolnozmienny (n.p. zależność przyspieszenia
grawitacyjnego od pływów Księżyca) możemy odfiltrować część szumu za
pomocą filtru dolnoprzepustowego.  Jeśli wiemy, że częstość sygnału
jest większa od jakiejś wartości możemy pozbyć się wolnozmiennego
dryfu aparaturowego za pomocą filtru górnoprzepustowego.

\subsection{Obliczanie splotu za pomocą transformaty Fouriera}
\label{sec:splot}

Splotem dwu funkcji $f(t)$ i~$g(t)$ nazywamy funkcję
\begin{equation}
  \label{eq:splot-definicja}
  (f \otimes g)(t) :=
  \int_{-\infty}^{\infty} d \tau f(\tau) g(t - \tau)
\end{equation}
Jeśli znamy transformaty Fouriera obu funkcji to splot można łatwo
wyliczyć korzystając z~tego, że transformata Fouriera splotu jest
iloczynem transformat czynników:
\begin{equation}
  \label{eq:splot-FT}
  \mathcal{FT}[f \otimes g](f) =
  \mathcal{FT}[f](t) \mathcal{FT}[g](t)
\end{equation}

Zazwyczaj w~zastosowaniach te dwie funkcje mają inne znaczenie.  Jedna
z~nich reprezentuje ,,prawdziwy'' sygnał i~jest funkcją przynajmniej
teoretycznie o~nieograniczonej domenie, a~druga jest funkcją
odpowiedzi aparatury i~jest wypikowana wokół zera.

Ponieważ dla dyskretnej (skończonej) transformaty Fouriera sygnały są
traktowane jak okresowe, więc gdy prawdziwy sygnał nie jest okresowy
to aby uniknąć efektu splatania końca danych z~początkiem należy dane
dopełnić na końcu zerami na długości około długości funkcji
splatanej.

Dla typowego przypadku splotu sygnału z~funkcją odpowiedzi, funkcję
odpowiedzi przedstawiamy w~następujący sposób: wartości dla dodatnich
czasów (opóźnień) na początku tablicy, wartości dla ujemnych opóźnień
na końcu tablicy (korzystając z~,,okresowości'': $h_{N-k} \equiv
h_{-k}$).  Dla uniknięcia efektów brzegowych sygnał dopełniamy zerami
na długości częsci ujemnej domeny funkcji odpowiedzi (t.j. na długości
$m$, gdzie \(m = \inf_k\{k\colon h_{-k} \neq 0\}\); plus ewentualne
dopełnienie zerami do najbliższej potęgi dwójki przy stosowani FFT
(algorytmu szybkiej transformaty Fouriera).

\subsection{Algorytm szybkiego mnożenia wielomianów}
\label{sec:mnozenie-wielomianow}

Niech
\begin{equation*}
  P(x) = \sum_{k=0}^{N-1} a_k x^k
\end{equation*}
będzie wielomianem stopnia $N - 1$.  Wielomian ten można opisać
jednoznacznie albo przez podanie wartości współczynników $a_k$, albo
przez podanie wartości $P(x_n)$ w~$N$ różnych punktach.  Gdy dwa
wielomiany $P(x)$ i~$Q(x)$ zadane są ciągami współczynników $(a_k)$
i~$(b_k)$, trudno jest policzyć ich iloczyn $P(x)Q(x)$ w~tej
reprezentacji.  Wiadomo natomiast, że $(P Q)(x_n) = P(x_n)Q(x_n)$.
Jeśli znajdziemy wartości iloczynu w~deg$(PQ) + 1$ punktach, będziemy
mogli znaleźć jego współczynniki.  Wystarczy policzyć wartości
wielomianów $P(x)$ i~$Q(x)$ w~tych samych deg$(PQ) + 1$ punktach.

Jeśli jako te punkty wybierzemy kolejne pierwiastki (zespolone)
odpowiedniego ($N$--tego, gdzie $N$ jest ilością wystarczającą do
wyznaczenia żądanego wielomianu) stopnia z~jedności:
\begin{equation}
  \label{eq:wielom-punkty}
  x_n = e^{2 \pi \imath n/N},
\end{equation}
to do obliczenia wartości wielomianu
przy danych jego wspólczynnikach można bedzie zastosować dyskretną
transformatę Fouriera
\begin{equation}
  \label{eq:wielom-FT}
  P(x_n) = P(e^{2 \pi \imath n/N}) 
  = \sum_{k=0}^{N-1} a_k e^{2 \pi \imath k n/N},
\end{equation}
a~do obliczenia wartości współczynników wielomianu mając zadane
wartości (do odtworzenia wielomianu) w~tych punktach odwrotną
transformatę Fouriera:

Zakładając, ze oba wielomany $P(x)$ i~$Q(x)$ są stopnia $N-1$,
wystarczy policzyć wartości w~$2(N - 1)$ punktach (gdyż takiego co
najwyżej stopnia może być ich iloczyn).  Rozszerzamy je w~tym celu do
wielomianów $2 N - 2$. stopnia (dodatkowe współczynniki przyjmując
równe zeru), przy pomocy szybkiej transformaty Fouriera wyliczamy
wartości w~zespolonych pierwiastkach z~jedności, mnożymy ,,punkt po
punkcie'' wynik i~odtransformujemy z~powrotem do reprezentacji za
pomocą współczynników.

To zadanie można także traktować jako przykład zastosowania splotu
dyskretnego dwu wektorów (w tym przypadku zakończonych zerami).
Mnożenie dwu wielomianów w~reprezentacji ,,współczynnikowej'' robione
jest za pomocą splotu.  Transformata Fouriera splotu przekłada się na
iloczyn transformat.  Jest to to samo zadanie przedstawione w~nieco
innym języku.

W~ten sam sposób można mnożyć przez siebie duże liczby, obliczając na
końcu przeniesienia.

\bigskip

\hfill{}Więcej informacji n.t. zastosowań transformaty Fouriera
w~następnej części

\bigskip
\bigskip

Jakub Narębski, Paweł Jankowski

\vfill

\begin{thebibliography}{9}
\bibitem{NRC} ,,Numerical Recipies in C: the Art of Scientific
  Computing'', \texttt{http://www.nr.com}

  
\bibitem{Macukow} Z. Fortuna, B. Macukow, J. Wąsowski.: ,,Metody
  numeryczne''
\end{thebibliography}

\end{document}
