[NumCS] Finish Monday morning lecture summary

semester3/numcs/numcs-summary.tex

@@ -103,6 +103,9 @@ The numbering should match the script's numbering exactly (apart from the cases
 \input{parts/01_interpolation/01_trigonometric/01_dft/02_fftshift.tex}
 \input{parts/01_interpolation/01_trigonometric/01_dft/03_linalg.tex}
 \input{parts/01_interpolation/01_trigonometric/02_fft.tex}
+\input{parts/01_interpolation/01_trigonometric/03_interpolation/00_intro.tex}
+\input{parts/01_interpolation/01_trigonometric/03_interpolation/01_zero-padding.tex}
+\input{parts/01_interpolation/01_trigonometric/04_error-estimation.tex}
 
 
 \end{document}

semester3/numcs/parts/01_interpolation/01_trigonometric/02_fft.tex

@@ -21,7 +21,7 @@ dass man die Berechnung einer DFT der Länge $n$ auf die Berechnung vieler DFTs
 
 Für den Algorithmus müssen folgende vier Optionen betrachtet werden:
 \begin{enumerate}[label=\Roman*]
-    \item Vektoren der Länge $N = 2m \Longrightarrow$ Laufzeit ideal
+    \item Vektoren der Länge $N = 2m \Longrightarrow$ Laufzeit gut
     \item Vektoren der Länge $N = 2^L \Longrightarrow$ Laufzeit ideal
     \item Vektoren der Länge $N = pq$ mit $p, q \in \Z \Longrightarrow$ Etwas langsamer
     \item Vektoren der Länge $N$, mit $N$ prim $\Longrightarrow$ ca. $\tco{N^2}$, besonders für $N$ gross
@@ -36,3 +36,5 @@ Wir formen die Fourier-Transformation um für den ersten Fall ($N = 2m$):
 Der zweite Fall ist einfach eine rekursive Weiterführung des ersten Falls,
 bei welchem dann das $m$ kontinuierlich weiter dividiert wird bis zum Trivialfall mit einer $1 \times 1$-Matrix.
 
+\fhlc{Cyan}{In NumPy} gibt es die Funktionen \texttt{np.fft.fft} (Vorwärts FFT), \texttt{np.fft.ifft} (Rückwärts FFT). 
+\texttt{scipy.fft} liefert dieselben Funktionen und sie sind oft etwas schneller als die von \texttt{numpy}

semester3/numcs/parts/01_interpolation/01_trigonometric/03_interpolation/00_intro.tex

@@ -0,0 +1,64 @@
+\newsection
+\subsection{Trigonometrische Interpolation}
+\subsubsection{Von Approximation zur Interpolation}
+Wir erinnern uns daran, dass wir die Fourier-Approximation durch den Abbruch der unendlichen Fourier-Reihe erhalten, oder in anderen Worten, wir verkleinern die Limiten der Summe.
+
+\fancyremark{DFT mit $N = 2n$ Koeffizienten an Punkten $\frac{l}{N}$ für $l = 0, 1, \ldots, N - 1$}
+
+Der Shift ist hier gegeben durch (für $k \geq 0$ ist $\gamma_k = \hat{f}_N(k)$ und für $k < 0$ ist $\gamma_k = \hat{f}_N(N + k)$)
+\begin{align*}
+    f_{N - 1}(x)                 & = \sum_{k = -n}^{n - 1} \gamma_k e^{2 \pi ikx} = \sum_{k = 0}^{n - 1} \gamma_k e^{2\pi ikx} + \sum_{k = -n}^{-1} \gamma_k e^{2\pi ikx} \\
+    \Leftrightarrow f_{N - 1}(x) & = \frac{1}{N} \left( \sum_{j = 0}^{N - 1} \left( f\left( \frac{j}{n} \right)
+        \sum_{k = -n}^{n - 1} e^{2\pi ik \left( x - \frac{j}{N} \right)} \right) \right)
+\end{align*}
+
+\vspace{-1pc}
+
+Wenn wir die Funktion nun an der Stelle $\frac{l}{N}$ auswerten so erhalten wir:
+\rmvspace
+\begin{align*}
+    f_{N - 1}\left( \frac{l}{N} \right) = \ldots = f\left( \frac{l}{N} \right)
+\end{align*}
+
+\vspace{-1.8pc}
+was aufgrund der Orthogonalität der diskreten Fourier-Vektoren funktioniert, welche besagt, dass $\displaystyle \sum_{k = -n}^{n - 1} \omega_N^{k(j - l)} = 0$, für alle $j \neq l$.
+Für $j = l$ ergibt die Summe $N$.
+
+Dies heisst also, dass die Fourier-Approximation die Interpolationsbedingungen an den Punkten $\frac{l}{N}$ erfüllt,
+also können wir die Lösung der Interpolationsaufgabe $p_{N - 1} \left( \frac{l}{N} \right) = f\left( \frac{l}{N} \right)$ f $l = 0, 1, \ldots, N - 1$ im Raum
+\rmvspace
+\begin{align*}
+    \mathcal{T}_N = \text{span}\{ e^{2\pi ijt} \divides j = - \floor{\frac{N - 1}{2}}, \ldots, \floor{\frac{N}{2}} \}
+\end{align*}
+
+\rmvspace\rmvspace
+folgendermassen finden können: 
+\begin{enumerate}[label=(\arabic*)]
+    \item Mittels Gleichungssystem $\sum_{j} \gamma_j e^{2\pi ijt_l} = f(t_l)$ für $l = 0, \ldots, N - 1$. Operationen: $\tco{N^3}$
+    \item Mittels FFT in $\tco{N \log(N)}$ Operationen, aber nur falls die Punkte äquidistant sind, also $t_l = \frac{l}{N}$.
+        Dann ist die Matrix des obigen Gleichungssystems $F^{-1}_N$
+\end{enumerate}
+
+\vspace{0.2cm}
+
+Unten findet sich Python code der mit den unterschiedlichen Methoden die Koeffizienten des Trigonometrischen Polynoms bestimmt.
+\rmvspace
+\begin{code}{python}
+    def get_coeff_trig_poly(t: np.ndarray, y: np.ndarray):
+    N = y.shape[0]
+    if N % 2 == 1:
+    n = (N - 1.0) / 2.0
+    M = np.exp(2 * np.pi * 1j * np.outer(t, np.arange(-n, n + 1)))
+    else:
+    n = N / 2.0
+    M = np.exp(2 * np.pi * 1j * np.outer(t, np.arange(-n, n)))
+    c = np.linalg.solve(M, y)
+    return c
+
+    N = 2**12
+    t = np.linspace(0, 1, N, endpoint=False)
+    y = np.random.rand(N)
+    direct = get_coeff_trig_poly(t, y)
+    using_fft = np.fft.fftshift(np.fft.fft(y) / N)
+    using_ifft = np.conj(np.fft.fftshift(np.fft.ifft(y)))
+\end{code}

semester3/numcs/parts/01_interpolation/01_trigonometric/03_interpolation/01_zero-padding.tex

@@ -0,0 +1,7 @@
+\newpage
+\subsubsection{Zero-Padding-Auswertung}
+Ein trigonometrisches Polynom $p_{N - 1}(t)$ kann effizient an den äquidistanten Punkten $\frac{k}{M}$ mit $M > N$ ausgewertet werden, für $k = 0, \ldots, M - 1$.
+Dazu muss das Polynom $p_{N - 1} \in \mathcal{T}_N \subseteq \mathcal{T}_M$ in der trigonometrischen Basis $\mathcal{T}_M$ neugeschrieben werden,
+in dem man \bi{Zero-Padding} verwendet, also Nullen im Koeffizientenvektor an den Stellen höheren Frequenzen einfügt.
+
+\TODO Insert cleaned up code from Page 95 (part of exercises)

semester3/numcs/parts/01_interpolation/01_trigonometric/04_error-estimation.tex

@@ -0,0 +1,99 @@
+\newsection
+\subsection{Fehlerabschätzungen}
+
+\begin{definition}[]{Konvergenz}
+    \begin{multicols}{2}
+        \fhl{Algebraische Konvergenz}\\
+        Wenn der Fehler $E(n) = \tco{\frac{1}{n^p}}$ mit $p > 0$ ist
+
+        \fhl{Exponentielle Konvergenz}\\
+        Wenn der Fehler $E(n) = \tco{q^n}$ mit $0 \leq q < 1$
+    \end{multicols}
+\end{definition}
+
+\inlineex Zur Fehlerbetrachtung verwenden wir drei Funktionen $f : [0, 1] \rightarrow \R$, welche wir mit trigonometrischer Interpolation an den Punkten $\frac{k}{N}$ approximieren:
+\begin{enumerate}[label=(\Roman*)]
+    % FIXME: Possibly wrong function definition in script
+    \item Stufenfunktion (periodische Fortsetzung von $f$) $f : [0, 1] \rightarrow \R$ mit $f(t) = \begin{cases}
+                  0 & \text{für } \left| t - \frac{1}{2} \right| > \frac{1}{4}    \\
+                  1 & \text{für } \left| t - \frac{1}{2} \right| \leq \frac{1}{4}
+              \end{cases}$
+    \item Periodische, glatte Funktion $h: \R \rightarrow \R$ mit $h(t) = \displaystyle \frac{1}{\sqrt{1 + \frac{1}{2} \sin(2\pi t)}}$
+          % TODO: Is it $h$ or $g$ here? $g$ makes little sense imho
+    \item Hutfunktion (periodische Fortsetzung von $h$) $g : [0, 1] \rightarrow \R$ mit $g(t) = \left| t - \frac{1}{2} \right|$
+\end{enumerate}
+Die untenstehende Abbildung \ref{fig:interpolation-error-examples} beinhaltet einen Plot, auf dem die Konvergenz in Abhängigkeit des Grades des Interpolationspolynoms aufgetragen ist.
+
+\begin{figure}[h!]
+    \begin{center}
+        \includegraphics[width=0.6\textwidth]{assets/01_interpolation/01_trigonometric/interpolation-error-examples.png}
+    \end{center}
+    \caption{Interpolierungsfehler der Beispiele. Algebraische Konvergenz für (I) und (III), exponentielle für (II)}
+    \label{fig:interpolation-error-examples}
+\end{figure}
+Auch hier tritt das Gibbs-Phänomen wieder an den Sprungstellen von $f(t)$ auf.
+Dies verursacht die Verlangsamung der Konvergenz in den Stellen, in welchen die Funktion nicht glatt ist.
+
+\newpage
+\stepcounter{all}
+\inlineex Sei für $\alpha \in [0, 1)$ $\displaystyle f(t) = \frac{1}{\sqrt{1 - \alpha \sin(2\pi t)}}$.
+Die Konvergenz ist exponentiell in $n$ und je kleiner $\alpha$, desto schneller ist sie.
+In der untenstehenden Abbildung \ref{fig:interpolation-error-convergence} sind einige Beispiele aufgetragen:
+\begin{figure}[h!]
+    \begin{center}
+        \includegraphics[width=0.6\textwidth]{assets/01_interpolation/01_trigonometric/interpolation-error-convergence.png}
+    \end{center}
+    \caption{Fehler bei der trigonometrischen Interpolation}\label{fig:interpolation-error-convergence}
+\end{figure}
+
+
+\setcounter{all}{6}
+\begin{theorem}[]{Aliasing}
+    Der k-te Fourier-Koeffizient des $N$-ten trigonometrischen Interpolationspolynoms unterscheidet sich vom $k$-ten Fourier-Koeffizienten von $f$ 
+    gerade um die Summe aller Fourier-Koeffizienten, die um ganze Vielfache von $N$ vom $k$-ten Fourier-Koeffizienten verschoben sind:
+    \begin{align*}
+        \hat{p}_N(k) - \hat{f}(k) = \sum_{j \neq 0 \in \Z} \hat{f}(k + jN)
+    \end{align*}
+\end{theorem}
+% FIXME: On page 98, just below the above theorem, there is a text I have no idea what he meant to say... in all honesty, I don't think he was sober when he wrote that
+
+\inlinecorollary Für $f \in \C^p([0, 1])$ mit $p \geq 1$ und $f$ $1$-periodisch, dann gilt: $|\hat{p}_N(k) - \hat{f}(k)| = \tco{(N^{-p})}$ für $|k| \leq \frac{N}{2}$
+
+Das heisst also, dass die Fourier-Koeffizienten von $f$ bei kleinen Frequenzen $\left( \text{hier } |k| < \frac{N}{2} \right)$
+sehr gut durch die Fourier-Koeffizienten des trigonometrischen Interpolationspolynoms approximiert werden.
+
+\begin{theorem}[]{Fehler der trigonometrischen Interpolation}
+    Falls $f$ $1$-periodisch ist und die Reihe $\sum_{k \in \Z} |\hat{f}(k)|$ absolut konvergiert, dann ist der Approximationsfehler definiert als:
+    \begin{align*}
+        |p_N(x) - f(x)| \leq 2 \sum_{|k| \geq \frac{N}{2}} |\hat{f}(k)| \mediumhspace \forall x \in \R
+    \end{align*}
+\end{theorem}
+Da durch diesen Satz die obere Schranke für den Approximationsfehler durch die schwer approximierbaren Fourier-Koeffizienten $\hat{f}(k)$ gegeben ist,
+heisst das folgendes für die Approximation von Polynomen von Grad $\deg(P(x)) < n$ für unser Approximationspolynom von Grad $\deg(P_N(x)) = n$:
+
+\fancycorollary{Abtasttheorem} Sei $f$ $1$-periodisch mit maximaler Frequenz $m$, also $\hat{f}(k) = 0 \smallhspace \forall |k| > m$. Falls $N > 2m$, dann gilt $p_N(x) = f(x) \smallhspace \forall x$
+
+\inlineex Ein Beispiel aus der Musik: Wir haben ein analoges Signal und wollen es digitalisieren. 
+Wir messen die Spannungswerte in äquidistanten Punkten. 
+Falls wir jedoch die Frequenz der Messung zu niedrig wählen, so kann ein total falsches Interpolationspolynom entstehen,
+wie in der untenstehenden Abbildung \ref{fig:aliasing} zu sehen:
+\begin{figure}[h!]
+    \begin{center}
+        \includegraphics[width=0.95\textwidth]{assets/01_interpolation/01_trigonometric/aliasing-in-music.png}
+    \end{center}
+    \caption{Aliasing für $f(t) = \cos(2\pi \cdot 19t)$}\label{fig:aliasing}
+\end{figure}
+Für unser Signal bedeutet das also, dass wir eine Art Verzerrung auf der Aufnahme haben, oder für Autoräder, dass es so scheint, als würden sich die Räder rückwärts drehen.
+
+\begin{theorem}[]{Fehlerabschätzung}
+    Sei $f^{(k)} \in L^2(0, 1) \smallhspace \forall k \in \N$, dann gilt:
+    \rmvspace
+    \begin{align*}
+        ||f - p_N(f)||_{L^2(0, 1)} \leq \sqrt{1 + c_k} N^{-k} ||f^{(k)}||_{L^2(0, 1)} \text{ wobei } c_k = 2 \sum_{l = 1}^{\infty} (2l - 1)^{-2k}
+    \end{align*}
+\end{theorem}
+Also, je mehr Ableitungen in $L^2(0, 1)$ liegen, desto kleiner ist der Fehler.
+
+
+Im Skript auf Seiten 101 und 102 gibt es einige Abbildungen die eine gewisse Intuition hinter der Approximation und den entstandenen Fehlern gibt.
+