\documentclass[10pt,a4paper]{article}

%% Polonizacja
\usepackage[T1]{fontenc}
\usepackage[latin2]{inputenc}
\usepackage{polski}
\usepackage{indentfirst}

%% Wzorki matematyczne
\usepackage{amsmath,amsfonts}

%% Odnośniki
\usepackage{xr}
\externaldocument[X:]{fft-splot}

\usepackage{varioref}
\renewcommand{\reftextbefore}{na poprzedniej stronie}
\renewcommand{\reftextafter}{na następnej stronie}
\renewcommand{\reftextcurrent}{na bieżącej stronie}
\renewcommand{\reftextfaraway}[1]{na \pageref{#1}. stronie}

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

%% Wykresy
\usepackage{pstricks,pst-plot,pst-node}

\newcommand{\abs}[1]{\lvert#1\rvert}

\begin{document}

\title{
  Zadania do Metod Numerycznych II A\\
  Grupy 2, 3. Seria 4.\\[1 ex]
  \large
  Splatanie i~odplatanie przy użyciu transformaty Fouriera
  }
\author{Jakub Narębski}
\date{09.05.2001~r.}
\maketitle



Do obliczania szybkiej transformaty Fouriera dwu funkcji rzeczywistych
na raz (co jest przydatne np. przy obliczaniu splotu za pomocą
transformaty Fouriera) mozna zastosować procedurę \lstinline!twofft!
z~NRC zadeklarowaną jako:
\begin{lstlisting}{}
void twofft(float data1[], float data2[], float fft1[], float fft2[],
  unsigned long n);
\end{lstlisting}
Oblicza ona transformatę Fouriera dwu rzeczywistych tablic wejściowych
\lstinline!data1[1..n]! i~\lstinline!data2[1..n]! o~rozmiarze $n$
i~zapisuje je w~dwu zespolonych tablicach wyjściowych
\lstinline!fft1[1..2n]!  i~\lstinline!fft2[1..2n]! o~rozmiarze
rzeczywistym $2n$.  Rozmiar tablicy $n$ musi byc całkowitą potęgą 2.
Uwaga: w~tej funkcji tablice indeksowane są od 1 do n!  (Można
stosować z~nią tablice indeksowane normalnie od 0, korzystając
z~faktu, że \lstinline!(arr-1)[1..n]! $\equiv$
\lstinline!arr[0..n-1]!.  Funkcja ta korzysta z~procedury
\lstinline!four1! z~NRC (która musi być dołączona do programu).




\section{Splatanie i~odplatanie}
\label{sec:splatanie-odplatanie}

Przetestować twierdzenie o~splocie dyskretnym.  Obliczyć splot sygnału
,,rozciągłego'' (np. funkcji okresowej piłokształtnej) z~funkcją
odpowiedzi impulsowej o~małej szerokości (np. funkcji Gaussa lub
prostokątnej bariery t.j. \(\Theta(a - \abs{x})\)) z~definicji:
\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}
oraz za pomocą transformaty Fouriera, korzystając z~twierdzenia
o~splocie dyskretnym:
\begin{equation}
  \label{eq:thm:splot-dyskretny}
  (R \otimes S)_j = R_j S_j.
\end{equation}
Funkcja odpowiedzi impulsowej warto przyjąć znormalizowaną,
tzn. przeskalować ją tak, aby $\sum r_k = 1$.


\psset{xunit=0.5cm,yunit=0.8cm}
\savedata{\signal}[{
  { 0,0} { 1,1} { 2,0} { 3,-1}
  { 4,0} { 5,1} { 6,0} { 7,-1}
  { 8,0} { 9,1} {10,0} {11,-1}
  {12,0} {13,1} {14,0} {15,-1}
  {16,0}
  }]
\savedata{\signalleft}[{
  {-1,-1} { 0,0}
  }]
\savedata{\signalright}[{
  {16, 0} {17,1}
  }]
\savedata{\signalother}[{
  { 0,0} {0.1,0.1} {0.5,0.5} {0.9,1.9} { 1,2} 
  {1.1,1.9} { 1.5,0.3} { 1.6,0} {2,-1} {2.8, -1.7}
  {3.5,-1} {4,-0.3} {4.5, -0.1} {5, 0.1} {6, 1.2} {8,1}
  {9,-0.3}
  {9.9,0.35} {10,0.4} {10.1,0.35}
  {11,-0.25}
  {11.25,0}
  {11.5,0.3}
  {12,0.5} {13,0.25} {14,0.125} {15,0.0625}
  {16,0}
  }]
\savedata{\response}[{
  {0,1} 
  {0.5,1}  {0.5,0} 
  {15.5,0} {15.5,1} 
  {16, 1}
  }]
\savedata{\responseext}[{
  {0,1} 
  {0.5,1} {0.5,0} 
  {16,0}  {16,1} 
  {16.5, 1}
  }]
\savedata{\responseleft}[{
  {-1,0} {-0.5,0} {-0.5,1} {0,1}
  }]
\savedata{\responseright}[{
  {16,1} {16.5,1} {16.5,0} {17,0}
  }]


\begin{figure}[htb]
  \begin{pspicture}(20,7)
    \rput(2.5,4){
      \dataplot[plotstyle=line,showpoints=false]{\signal}
      \dataplot[plotstyle=line,showpoints=false,linestyle=dotted]{\signalleft}
      \dataplot[plotstyle=line,showpoints=false,linestyle=dotted]{\signalright}
      \psline{<-}(0,1.5)(0,-1.5)
      \psline{->}(-1,0)(17,0)
      \psline[linestyle=dashed](16,1.4)(16,-1.4)
      \uput[0](0,1.5){$s(t)$}
      \pcline{<->}(0,-1.3)(16,-1.3)
      \lput*{:U}{$N$}}
    \rput(2.5,0.5){
      \dataplot[plotstyle=line,showpoints=false]{\response}
      \dataplot[plotstyle=line,showpoints=false,linestyle=dotted]{\responseleft}
      \dataplot[plotstyle=line,showpoints=false,linestyle=dotted]{\responseright}
      \psline{<-}(0,1.5)(0,0)
      \psline{->}(-1,0)(17,0)
      \psline[linestyle=dashed](16,1.4)(16,0)
      \uput[0](0,1.5){$r(t)$}
      \psline[linewidth=.4pt,linestyle=dashed](-0.5,0)(-0.5,-0.31)
      \psline[linewidth=.4pt,linestyle=dashed]( 0.5,0)( 0.5,-0.31)
      \pcline{<->}(-0.5,-0.3)( 0.5,-0.3)
      \bput*{:U}{$M$}}
  \end{pspicture}
  \caption{Sygnał okresowy}
  \label{fig:okresowa}
\end{figure}

\begin{figure}[htb]
  \begin{pspicture}(18,9)
    \rput(0.5,5){
      \dataplot[plotstyle=curve,showpoints=false]{\signalother}
      \psline{<-}(0,2.5)(0,-2)
      \psline{->}(0,0)(17.5,0)
      \psline[linestyle=dotted](16  ,0)(16  ,-1.9)
      \psline[linestyle=dotted](16.5,0)(16.5,-1.9)
      \uput[0](0,2.5){$s(t)$}
      \pcline{<->}(0,-1.85)(16,-1.85)
      \bput*{:U}{$N$}
      \pcline{|-|}(16,-1.85)(16.5,-1.85)
      \bput*{:U}{$\tfrac{M}{2}$}}
    \rput(0.5,0.5){
      \dataplot[plotstyle=line,showpoints=false]{\responseext}
      \psline{<-}(0,1.5)(0,0)
      \psline{->}(0,0)(17.5,0)
      \psline[linestyle=dotted](16.5,1.4)(16.5,0)
      \uput[0](0,1.5){$r(t)$}}
  \end{pspicture}
  \caption{Sygnał nieokresowy, zanikający}
  \label{fig:nieokresowa}
\end{figure}

Uwaga: sygnał możemy traktować albo jako funkcję okresową, albo jako
funkcję nieokresową (wycinek funkcji nieskończonej w~czasie, np. jeśli
sygnał jest funkcją zanikającą).  Należy postępować stosownie do
przypadku:
\begin{itemize}
\item Sygnał \textbf{okresowy}: \(s_{k + l N} = s_k\).  Przy liczeniu
  splotu z~definicji należy pamiętać o~tym, że \(s_{j-k} \equiv 
  s_{j-k \bmod N}\).  Funkcję odpowiedzi można traktować albo jako
  $r_k\colon -M/2 < k \leq M/2$ i~sumować po $k$ od $-M/2+1$ do $M/2$
  (sięgając do odpowiednich elementów tablicy $r[k]$), albo używać
  w~postaci ,,zawiniętej'' $r_k\colon 0 \leq k < N$ i~sumować po $k$
  od $0$ do $N-1$.  $M$ oznacza tutaj szerokość funkcji odpowiedzi;
  można przyjąć ją równą $N$.

  Przy liczeniu splotu za pomocą transformaty Fouriera nie dopełnia
  się sygnału zerami.  Funkcja odpowiedzi jednostkowej $r$ powinna być
  zapisana w~postaci ,,zawiniętej'', korzystając z~tożsamości \(r_{-k}
  = r_{N-k}\).  Niezerowe elementy powinny się znajdować na końcu i~na
  początku tablicy.

  Patrz także rysunek~\vref{fig:okresowa}.
  
\item Sygnał \textbf{nieokresowy}: \(s_k \neq 0,\quad k=0,\dots N\).
  Przy liczeniu splotu z~definicji należy przyjąć \(s_k = 0\) dla \(k
  < 0\) i~dla \(k \geq N\).  Funkcję odpowiedzi należy traktować jako
  $r_k\colon -M/2 < k \leq M/2$ i~sumować po $k$ od $-M/2+1$ do $M/2$
  (sięgając do odpowiednich elementów tablicy $r[k]$), lub odpowiednio
  przeliczyć indeksy.
  
  Przy liczeniu splotu za pomocą transformaty Fouriera należy dopełnić
  sygnał $s$ zerami (patrz rozdział~\ref{X:sec:zero-padding}
  w~\texttt{fft-splot.ps}) na długości połowy szerokości funkcji
  odpowiedzi, a~dokładniej na większej z~szerokości ,,dodatniej''
  i~,,ujemnej'' funkcji odpowiedzi jednostkowej $r$.  Nastepnie należy
  zapisać $r$ do tablicy o~długości równej długości tablicy
  zawierającej sygnał rozszerzony zerami, w~postaci ,,zawiniętej'',
  korzystając z~tożsamości \(r_{-k} = r_{N-k}\).  Uzupełniające zera
  powinny się znaleźć pośrodku tablicy.

  Patrz także rysunek~\vref{fig:nieokresowa}.

\end{itemize}

Należy pamiętać, że transformata Fouriera funkcji rzeczywistej jest
funkcją \emph{zespoloną} o~własności $F(-f) = [F(f)]^{\ast}$.

Funkcja odpowiedzi impulsowej powinna być skoncentrowana wokół zera
(zgodnie ze wzorem \(r_{-k} \equiv r_{N-k}\)), tak że niezerowe
elementy są na początku i~na końcu tablicy.  Można stosować dowolną
z~procedur do obliczania transformaty Fouriera: własną,
\lstinline!four1!, \lstinline!realft! czy \lstinline!twofft!.
Porównać wynik obliczony z~definicji z~wynikiem za pomocą transformaty
Fouriera.

\bigskip

Odpleść funkcję odpowiedzi impulsowej z~obliczonego splotu (najlepiej
z~definicji): $\widetilde{S}(f) = C(f)/R(f)$, gdzie $c(t) = r(t)
\otimes s(t)$.  Porównać wynik odplecenia $\widetilde{s}_k$ z~sygnałem
oryginalnym.  Uwaga na wartości zerowe w~transformacie Fouriera
funkcji odpowiedzi!  Wartośc zerowa $R_n$ oznacza, że przy splataniu
utraciliśmy pewne informacje: bezpieczne jest przyjąć $\widetilde{S}_k
= 0$ dla $R_k = 0$.  \emph{Dla chętnych:} sprawdzić wpływ dodania
szumu do wyniku splotu $c$ na błąd odzyskania oryginalnego sygnału za
pomocą odplecenia $\widetilde{s}_k - s_k$.

\bigskip

Odpleść funkcję odpowiedzi impulsowej za pomocą procedury iteracyjnej
zadanej wzorem
\begin{equation}
  \label{eq:odplatanie-iteracyjne}
  q_{\text{nowe}} = q + (c - r \otimes q),
\end{equation}
gdzie $q$ oznacza wynik w~kolejnej iteracji, zaś $c$ jest mierzonym
sygnałem (po splocie) od którego chcemy odpleść funkcję odpowiedzi
jednostkowej $r$.  Początkowo (w pierwszej iteracji) $q = c$.  Splot
$r \otimes q$ można wyliczyć dowolnym sposobem, nawet przy użyciu
procedury \lstinline!convlv! z~NRC splatającej dwie funkcje
rzeczywiste z~NRC.



\section{Filtr Wienera (nieobowiązkowe)}
\label{sec:filtr-Wienera}

\emph{Zadanie dodatkowe, nieobowiązkowe.}  Znaleźć oryginalny sygnał
na podstawie pomiaru ze znaną funkcją odpowiedzi jednostkowej,
obarczonego szumem.  Dany jest sygnał wynikowy \(c(t) = r(t) \otimes
s(t) + n(t)\) oraz funkcja odpowiedzi jednostkowej $r(t)$ (lub jej
dobre przybliżenie).  Zastosować w~tym celu filtr optymalny Wienera
wzorem:
\begin{equation}
  \label{eq:Wiener:stosowanie}
  \widetilde{S}(f) = \frac{C(f)\varPhi(f)}{R(f)}.
\end{equation}
Funkcja $\varPhi(f)$ filtru zadana jest wzorem
\begin{equation}
  \label{eq:Wiener:wzor-na-filtr}
  \varPhi(f) = \frac{\abs{H(f)}^2}{\abs{H(f)}^2 + \abs{N(f)}^2},
\end{equation}
gdzie \(H(f) = R(f)S(f)\) (czyli \(h(t) = r(t) \otimes s(t)\)).  Model
$\abs{H(f)}^2$ oraz oszacowanie widma szumu $\abs{N(f)}$ można znaleźć
na podstawie widma mocy $P_c(f) = \abs{C(f)}^2$ wynikowego sygnału
(patrz rozdział~\ref{X:sec:Wiener-przyblizenia}
w~\texttt{fft-splot.ps}), lub w~jakiś inny sposób (proszę go opisać).
Porównać filtrowanie Wienera ze stosowaniem kilku kroków metody
iteracyjnej odplatania, danej wzorem~\eqref{eq:odplatanie-iteracyjne}.

\bigskip

\emph{Dodatkowo} zamiast badania sygnału jednowymiarowego można badać
splatanie i~rozplatanie oraz filtr Wienera dla obrazów
jednokolorowych.  Dwuwymiarową transformatę Fouriera można obliczyć
wyliczając jednowymiarową transformatę Fouriera dla kolumn,
a~następnie dla wierszy (albo odwrotnie).  Wartości wykraczające poza
zakres mapy kolorów nalezy obciąć.  Podany tutaj iteracyjny sposób
odplatania stosowany był pierwotnie do analizy obrazów
astronomicznych.

\vfill

\begin{flushright}
  Jakub Narębski
\end{flushright}

\end{document}
