diff --git a/semester3/numcs/assets/01_interpolation/01_trigonometric/aliasing-in-music.png b/semester3/numcs/assets/01_interpolation/01_trigonometric/aliasing-in-music.png new file mode 100644 index 0000000..fcdeeef Binary files /dev/null and b/semester3/numcs/assets/01_interpolation/01_trigonometric/aliasing-in-music.png differ diff --git a/semester3/numcs/assets/01_interpolation/01_trigonometric/interpolation-error-convergence.png b/semester3/numcs/assets/01_interpolation/01_trigonometric/interpolation-error-convergence.png new file mode 100644 index 0000000..c216ee6 Binary files /dev/null and b/semester3/numcs/assets/01_interpolation/01_trigonometric/interpolation-error-convergence.png differ diff --git a/semester3/numcs/assets/01_interpolation/01_trigonometric/interpolation-error-examples.png b/semester3/numcs/assets/01_interpolation/01_trigonometric/interpolation-error-examples.png new file mode 100644 index 0000000..db0a648 Binary files /dev/null and b/semester3/numcs/assets/01_interpolation/01_trigonometric/interpolation-error-examples.png differ diff --git a/semester3/numcs/numcs-summary.pdf b/semester3/numcs/numcs-summary.pdf index f26699a..2f5f779 100644 Binary files a/semester3/numcs/numcs-summary.pdf and b/semester3/numcs/numcs-summary.pdf differ diff --git a/semester3/numcs/numcs-summary.tex b/semester3/numcs/numcs-summary.tex index d135cb4..acf18b6 100644 --- a/semester3/numcs/numcs-summary.tex +++ b/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} diff --git a/semester3/numcs/parts/01_interpolation/01_trigonometric/02_fft.tex b/semester3/numcs/parts/01_interpolation/01_trigonometric/02_fft.tex index e0976a0..bff51cd 100644 --- a/semester3/numcs/parts/01_interpolation/01_trigonometric/02_fft.tex +++ b/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} diff --git a/semester3/numcs/parts/01_interpolation/01_trigonometric/03_interpolation/00_intro.tex b/semester3/numcs/parts/01_interpolation/01_trigonometric/03_interpolation/00_intro.tex new file mode 100644 index 0000000..78f24e3 --- /dev/null +++ b/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} diff --git a/semester3/numcs/parts/01_interpolation/01_trigonometric/03_interpolation/01_zero-padding.tex b/semester3/numcs/parts/01_interpolation/01_trigonometric/03_interpolation/01_zero-padding.tex new file mode 100644 index 0000000..2466b2e --- /dev/null +++ b/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) diff --git a/semester3/numcs/parts/01_interpolation/01_trigonometric/04_error-estimation.tex b/semester3/numcs/parts/01_interpolation/01_trigonometric/04_error-estimation.tex new file mode 100644 index 0000000..63dccad --- /dev/null +++ b/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. +