\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{\norm}[1]{\lVert#1\rVert}
\newcommand{\field}[1]{\mathbb{#1}}

\newtheorem{tw}{Twierdzenie}

\begin{document}
\title{
  Zastosowania transformaty Fouriera\\
  część druga: splatanie, filtr Wienera\\[1 ex]
  \large
  Metody Numeryczne II A}
\author{Jakub Narębski}
\date{09.05.2001~r.}

\maketitle

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

%%%
\section{Splatanie i~odplatanie przy użyciu transformaty Fouriera}
\label{sec:splatanie-i-odplatanie}



\subsection{Własności splotu}
\label{sec:wprowadzenie}

Splotem dwu funkcji $f(t)$ i~$g(t)$ nazywamy funkcję
\begin{equation}
  \label{eq:def:splot}
  (f \otimes g)(t) :=
  \int_{-\infty}^{\infty} d \tau f(\tau) g(t - \tau)
\end{equation}
Zachodzi następujące twierdzenie:
\begin{tw}[o splocie]
  \label{thm:o-splocie}
  Transformata Fouriera splotu dwu funkcji $f$ i~$g$ jest iloczynem
  transformat Fouriera czynników:
  \begin{equation}
    \label{eq:splot-FT}
    \begin{gathered}
      \mathcal{FT}[f \otimes g](t) =
      \mathcal{FT}[f](t) \, \mathcal{FT}[g](t),\\
      (F \otimes G)(f) = F(f) G(f).
    \end{gathered}
  \end{equation}
\end{tw}

Splot jest operacją przemienną, tzn. \(r \otimes s = s \otimes r\).
Zazwyczaj jednak w~zastosowaniach te dwie splatane funkcje mają inne
znaczenie.  Jedna z~nich, powiedzmy $s$, reprezentuje ,,prawdziwy''
sygnał (niezaburzone wyniki) i~jest funkcją o~(przynajmniej
teoretycznie) nieograniczonej dziedzinie.  Funkcja $r$ (druga funkcja)
jest natomiast funkcją odpowiedzi aparatury, i~ma skończoną dziedzinę,
zazwyczaj odcinek zawierający zero.  Inaczej mówiąc, funkcja ta ma
skończoną szerokość.

Operacja splotu \(r \otimes s\) opisuje rozmycie sygnału (danych) $s$
zgodnie z~przepisem danym funkcją odpowiedzi $r$.  Np. dla funkcji
odpowiedzi będącej deltą Diraca dla czasu $t_0$ tzn. \(r(t) = \delta
(t - t_0)\), splot \(r \otimes s\) daje sygnał przesunięty o~$t_0$:
\((r \otimes s)(t) = s(t - t_0)\). 

\smallskip

W~przypadku dyskretnym sygnał $s(t)$ jest reprezentowany przez
wartości $s_k$ w~punktach próbkowania \(t_k = t_0 + k
\Delta\).  Funkcja odpowiedzi jest także zadana przez wartości
w~dyskretnych punktach $r_k$.  W~typowych zastosowaniach funkcja
odpowiedzi ma skończony czas trwania $M$, tzn. jest różna od zera
tylko dla \(-M/2 < k \leq M/2\).  Wówczas dyskretny splot zdefiniowany
jest wzorem:
\begin{equation}
  \label{eq:def:splot-dyskretny}
  (r \otimes s)_j := \negmedspace
  \sum_{k = -\frac{M}{2}+1}^{\frac{M}{2}} r_k s_{j-k} 
\end{equation}
Zachodzi następujące twierdzenie
\begin{tw}[o splocie dyskretnym]
  \label{thm:o-splocie-dyskretnym}
  Jeśli sygnał $s_j$ jest okresowy z~okresem $N$, tak że jest ściśle
  określony przez $N$ wartości \(s_0,\dots,s_{N-1}\), to transformata
  Fouriera jego splotu z~funkcją odpowiedzi $r_k$ o~czasie trwania
  (szerokości) $N$ (tzn. \(r_k \neq 0 \iff -N/2 < k \leq N/2\)) zadana
  jest wzorem:
  \begin{equation}
    \label{eq:splot-dyskretny-FT}
    \begin{gathered}
      \mathcal{FT}[r \otimes s]_j =
      \mathcal{FT}[r]_j \mathcal{FT}[s]_j,\\
      (R \otimes S)_j = R_j S_j,
    \end{gathered}
  \end{equation}
  gdzie $\mathcal{FT}[r]$ jest transformatą Fouriera wartości $r_k$
  \((k = 0,\dots,N-1)\) ułożonych zgodnie ze wzorem \(r_{N-k} \equiv
  r_{-k}\) dla \(k > N/2\).
\end{tw}

\smallskip

Zazwyczaj funkcja odpowiedzi $r$, czy to w~przypadku ciągłym, czy też
dyskretnym, jest znormalizowana, tzn. ma jednostkowe pole powierzchni.
Przykładem [nieznormalizowanych] funkcji odpowiedzi są funkcja Gaussa
\(\exp(-x^2/\sigma^2)\) i~funkcja schodkowa \(\Theta(a - \abs{x})\).

\subsection{Uwzględnianie efektów brzegowych za pomocą uzupełniania
  zerami}
\label{sec:zero-padding}

W~założeniach twierdzenia~\ref{thm:o-splocie-dyskretnym} o~splocie
dyskretnym znajduje się żądanie, by sygnał był okresowy oraz by
funkcja odpowiedzi miała szerokość równą długości (okresowi) sygnału.
Często chcemy stosować twierdzenie o~splocie (do splatania
i~rozplatania) dla funkcji (danych) nie spełniających tych warunków.

Jeśli chodzi o~warunek narzucany na funkcję odpowiedzi $r_k$, to nie
stanowi on kłopotu.  Szerokość funkcji odpowiedzi $M$ jest mniejsza,
zazwyczaj dużo mniejsza od długości sygnału $N$: $N \ll M$.  Wystarczy
zatem dopełnić funkcję $r$ zerami do szerokości $N$, tzn. postawić
\begin{math}
  r_k = 0 \text{ dla } -N/2 < k < M/2 \land M/2 \leq k \leq N/2
\end{math}
(pamiętając o~tym, że $r_k$ zawijamy, by mieć \(k = 0,\dots,N-1\)).

Nieco kłopotu jest natomiast z~warunkiem na $s$.  Zazwyczaj dane $s$
traktujemy jako wycinek nieskończonego w~czasie sygnału.  Dla funkcji
okresowej i~funkcji odpowiedzi o~skończonej szerokości wartości
z~początku okresu splatane są z~wartościami z~końca okresu.  Aby tego
uniknąć należy sygnał $s$ dopełnić zerami na odpowiedniej długości;
splot z~zerem jest równy zeru: \(f \otimes 0 = 0\).

Aby uniknąć fałszowania splotu na początku przedziału należy zatem
dopełnić sygnał tyloma zerami, ile wynosi największy ujemny indeks,
taki że $r_{-k}$ jest różne od zera.  Aby uniknąć fałszowania wyniku na
końcu przedziału (splatania ostatniej wartości z~wartościami
z~początku) należy natomiast dopełnić sygnał tyloma zerami, ile wynosi
największy dodatni indeks, taki że $r_k$ jest różne od zera.
Podsumowując: należy dopełnić większą z~szerokości ,,ujemnej'' $M_-$
i~,,dodatniej'' $M_+$ funkcji odpowiedzi $r$, gdzie
\begin{subequations}
  \begin{align}
    \label{eq:M-minus}
    M_- &= \min \{k \in \field{N}\colon r_{-k} \neq 0\},\\
    \label{eq:M-plus}
    M_+ &= \min \{k \in \field{N}\colon r_{k}{}_{\hphantom{-\mbox{}}} \neq 0\}.
  \end{align}
\end{subequations}
Dla symetrycznej funkcji odpowiedzi potrzeba dopełnić sygnał $M/2$
(gdzie $M$ jest szerokością funkcji $r$) zerami, by uniknąć
fałszowania splotu.

Przy stosowaniu algorytmu Cooleya--Tukeya szybkiej transformaty
Fouriera należy dodatkowo po dopełnieniu zerami w~celu uniknięcia
fałszowania splotu dopełnić zerami zarówno sygnał jak i~funkcję
odpowiedzi (pamiętając o~,,zawinięciu'' funkcji odpowiedzi) zerami do
najbliższej potęgi dwójki.

Wynik splotu dla indeksów, dla których $s$ było wypełnione zerami nie
jest znaczący (jest to wynik splatania końca z~danymi z~początku)
i~należy go pominąć.



\subsection{Odplatanie przy użyciu FFT}
\label{sec:odplatanie-FFT}

Mamy dany sygnał zadany za pomocą $N$ liczb rzeczywistych \(s_k, k =
0,\dots,N-1\) (odpowiednio dopełniony zerami) i~funkcję odpowiedzi
dopełnioną zerami do długości $N$, z~wartościami ,,zawiniętymi'', tak
że niezerowe wartości są skoncentrowane na początku i~na końcu
tablicy.  Możemy obliczyć dyskretny splot $r \otimes s$ używając FFT
do obliczenia dyskretnej transformaty Fouriera $r$ i~$s$, następnie
mnożąc te transformaty przez siebie (pamiętając, że składają się one
z~liczb \emph{zespolonych}) i~wyznaczając wynik wyliczając odwrotną
transformatę Fouriera (rzeczywistą) iloczynu.

\emph{Rozplataniem} (odplataniem) nazywamy procedurę odwrotną do
splatania.  Mamy dany wynik doświadczenia (z~aparatury) \(c = r
\otimes s\) oraz funkcję odpowiedzi układu pomiarowego $r$ (lub jej
przybliżenie) i~chcemy odzyskać sygnał $s$.  Dla splotu dyskretnego
rozplatanie w~dziedzinie czasowej wymaga rozwiązania układu $N$
równań; jest to zazwyczaj nierealistyczne.  Przy użyciu transformaty
Fouriera możemy uzyskać transformatę Fouriera oryginalnego sygnału
(niezaburzonych przez aparaturę pomiarową danych) przez podzielenie
transformatę Fouriera wyniku przez transformatę Fouriera funkcji
odpowiedzi.

Procedura ta może się nie powieść \emph{matematycznie} gdy
transformata Fouriera funkcji odpowiedzi $R_n$ jest równa dokładnie
zero dla pewnych wartości $n$, tak że nie możemy przez nią podzielić.
Oznacza to, że przy splataniu (stosując dany układ pomiarowy)
utraciliśmy informacje o~pewnych częstościach oryginalnego sygnału.
(Należy je zatem położyć równe zeru, jeśli napotkamy na taką sytuację
--- odzyskany sygnał nie będzie wtedy dokładny).  Poza tym metoda ta
wrażliwa jest na szumy w~danych wejściowych oraz na dokładność
znajomości funkcji odpowiedzi.  W~niektórych przypadkach odplatanie
może dać nonsensowne wyniki: należy wtedy zastosować \emph{optymalne
  filtrowanie}, które jest omówione w~rozdziale~\ref{sec:Wiener}.

\smallskip

Możemy przy użyciu algorytmu obliczającego zespoloną transformatę
Fouriera wyliczyć transformatę dwu funkcji (tablic) rzeczywistych.
Korzystamy z~faktu, że dla funkcji rzeczywistej $f$ ($f_k \in
\field{R}$) transformata Fouriera spełnia warunek
\begin{align}
  \label{eq:tworealft:real}
  F(-f_n) \equiv F_{N-n} &= (F_n)^{\ast} \equiv [F(f)]^{\ast},\\
\intertext{zaś dla funkcji czysto urojonej $g$ ($g_k \in \imath
  \field{R}$) mamy}
  \label{eq:tworealft:imag}
  G(-f_n) \equiv G_{N-n} &= -(G_n)^{\ast} \equiv -[G(f)]^{\ast}.
\end{align}
Zatem możemy obliczyć transformatę Fouriera dwu funkcji rzeczywistych
na raz (tak jak np. przy splataniu i~odplataniu) zapisując te funkcje
jako odpowiednio część rzeczywistą i~urojoną zespolonej funkcji
(tablicy) \(h_k = f_k + \imath g_k\), licząc zepoloną transformatę
Fouriera i~odtwarzając transformaty oryginalnych funkcji korzystając
z~własności~\eqref{eq:tworealft:real} i~\eqref{eq:tworealft:imag}:
\begin{align}
  \label{eq:tworealft:F}
  F_k &= 
  \frac{1}{2}(\Re H_k + \Re H_{N-k}) 
  + \imath \,
  \frac{1}{2}(\Im H_k - \Im H_{N-k}),\\
  \label{eq:tworealft:G}
  G_k &= 
  \frac{1}{2}(\Im H_k - \Im H_{N-k}) 
  - \imath \,
  \frac{1}{2}(\Re H_k + \Re H_{N-k}),
\end{align}
dla $0 < k < N/2$.  Po dalsze szczegóły odsyłam do~\cite{NRC}.



\section{Optymalny filtr Wienera przy użyciu FFT}
\label{sec:Wiener}

\subsection{Filtr Wienera}
\label{sec:filtr-Wienera}

Czasami potrzebujemy znaleźć oryginalny sygnał mając dany zaszumiony
sygnał z~aparatury pomiarowej (aparatura rozmywa sygnał zgodnie
z~funkcją odpowiedzi).  Oznaczając oryginalny sygnał przez $s(t)$,
funkcję odpowiedzi aparatury przez $r(t)$, a~szum przez $n(t)$, chcemy
znaleźć sygnał $s(t)$ (lub jego dobre przybliżenie) mając dane $r(t)$
oraz \(c(t) = r(t) \otimes s(t) + n(t)\).

Proste odplecenie zgodnie ze wzorem \(S(f) = C(f)/R(f) - N(f)/R(f)\)
nie daje oczekiwanych wyników, gdyż jeśli $R(f)$ jest małe to
wzmacniany jest zarówno sygnał, jak i~szum.

W~celu obliczenie $s(t)$ znajdujemy \emph{filtr optymalny} ($\phi(t)$
w~dziedzinie czasowej, $\varPhi(f)$ w~dziedzinie częstości) taki że
$\widetilde{s}(t)$ zadane wzorem (na transformatę Fouriera):
\begin{equation}
  \label{eq:Wiener-stosowanie}
  \widetilde{S}(f) = \frac{C(f) \varPhi(f)}{R(f)}
\end{equation}
(gdzie wielkie litery oznaczają odpowiednie transformaty Fouriera),
jest bliskie w~sensie odległości w~$\mathcal{L}^2$ oryginalnemu
sygnałowi $s(f)$, tzn. takie by
\begin{equation}
  \label{eq:odleglosc-L2}
  \norm{\widetilde{s} - s}^2_{L^2}
  = \int_{-\infty}^{\infty} \abs{\widetilde{s}(t) - s(t)}^2 \:dt
  = \int_{-\infty}^{\infty} \abs{\widetilde{S}(f) - S(f)}^2 \:df
  \quad \text{było minimalne}.
\end{equation}

Wprowadźmy oznaczenie \(h = r \otimes s\) na splot sygnału z~funkcją
odpowiedzi.  Zakładając, że $H$ jest nieskorelowane z~szumem $N$, oraz
że filtr $\varPhi$ jest funkcją rzeczywistą dostajemy następujący wzór
na filtr Wienera:
\begin{equation}
  \label{eq:Wiener-wzor}
  \varPhi(f) = \frac{\abs{H(f)}^2}{\abs{H(f)}^2 + \abs{N(f)}^2}.
\end{equation}



\subsection{Znajdowanie przybliżeń do obliczenia filtru Wienera}
\label{sec:Wiener-przyblizenia}

Aby znaleźć filtr Wienera musimy mieć jakiś sposób szacowania
niezależnie $\abs{H}^2$ i~$\abs{N}^2$.  Dany mamy tylko mierzony
sygnał $C$.  Możemy wyznaczyć te wielkości obliczając widmo mocy
sygnału $C$ (dla dodatnich częstości), które jest proporcjonalne do
sumy \(\abs{H}^2 + \abs{N}^2\):
\begin{equation}
  \label{eq:widmo-mocy}
  \abs{H(f)}^2 + \abs{N(f)}^2 \approx P_c(f) = \abs{C(f)}^2
  \quad 0 \leq f < f_c.
\end{equation}
Wykres widma mocy będzie zazwyczaj pokazywał widmo sygnału $H$
skoncentrowane wokół pewnej częstości (lub pewnych częstości)
wystające ponad ciągłe widmo szumu (najczęściej stałe lub
wolnozmienne).  Na podstawie tego widma możemy wyznaczyć jakiś model
na $H(f)$ oraz przybliżenie na widmo szumu $N$.  Filtr jest
symetryczny w~dziedzinie częstości; dopełniamy go ze wzoru
\(\varPhi(-f) = \varPhi(f)\).

Ponieważ wzór na filtr Wienera uzyskuje się za pomocą procedury
minimalizacji, jakość filtrowania zależy \emph{kwadratowo} od
dokładności filtru.  Zatem nawet zgrubnie wyznaczony filtr optymalny
daje dobre wyniki.  \textbf{Uwaga:} nie należy stosować filtrowania
optymalnego iteracyjnie: dalsze stosowanie nie przynosi poprawy,
a~procedura iteracyjna zbieżna jest do sygnału o~$H(f) \equiv 0$.



\subsection{Iteracyjne odplatanie bez potrzeby użycia FFT}
\label{sec:iteracyjne-odplatanie}

Można także próbować odtworzyć sygnał metodą iteracyjną.  W~tym celu
odejmujemy splot z~funkcją odpowiedzi.  Przyjrzyjmy się równaniu
iteracyjnemu\cite{Imag}:
\begin{equation}
  \label{eq:iteracyjne-odplatanie}
  q^{(k+1)} = q^{(k)} + (c - r \otimes q^{(k)}),
\end{equation}
gdzie $k$ jest indeksem iteracji (\(q^{(0)} = c\)), $q^{(k)}$
estymacją sygnału dla \mbox{$k$--tej} iteracji, $c$ mierzonym
sygnałem, a~$r$ funkcją odpowiedzi impulsowej.

Widzimy, że w~miarę wzrostu indeksu iteracji $k$ różnica \((c - r
\otimes q^{(k)})\) zbliża się do zera i~pozostaje $s$, czyli
oryginalny sygnał.  Stosując tę metodę nie musimy w~ogóle używać
transformaty Fouriera, gdyż splot możemy wyliczać bezpośrednio
z~definicji w~dziedzinie czasowej.  Zazwyczaj potrzeba niewielu
iteracji do uzyskania optymalnego wyniku; dalsze stosowanie wzoru
iteracyjnego powoduje ,,zniszczenie'' danych.

Oryginalnie metoda ta wymyślona została do odplatania funkcji
odpowiedzi (zazwyczaj funkcji Gaussa) z~obrazów CCD w~astronomii.

\bigskip
\bigskip

\begin{flushright}
Jakub Narębski, Paweł Jankowski  
\end{flushright}


\vfill

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

  
\bibitem{Imag} C. D. Watkins, A. Sadun, S. Marenka: ,,Nowoczesne
  metody przetwarzania obrazu''.

\end{thebibliography}

\end{document}
