diff --git a/assignment_03/gibbs_incremental_reconstruction.ipynb b/assignment_03/gibbs_incremental_reconstruction.ipynb new file mode 100644 index 0000000..f050513 --- /dev/null +++ b/assignment_03/gibbs_incremental_reconstruction.ipynb @@ -0,0 +1,225 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "b95ff00b-7fff-4d33-87ef-0c03b1a2217d", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from setupFigure import SetupFigure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b17baa5-8598-4be8-ae73-67510d6cce6f", + "metadata": {}, + "outputs": [], + "source": [ + "def dft_coeff(x):\n", + " \"\"\"\n", + " Evaluate c_n = 1/N*sum_{k=0}^{N-1}x_k exp(-i*2*pi*nk/N)\n", + " :param x: samples of function with length N\n", + " \"\"\"\n", + " ns = np.size(x)\n", + " c = np.zeros(ns, dtype=np.complex64) # zero coefficients before summing\n", + " for n in range(ns):\n", + " for k in range(ns):\n", + " arg = -2*np.pi*1j*n*k/ns\n", + " c[n] = c[n]+x[k]*np.exp(arg)\n", + " c[n] = c[n]/ns\n", + " return c" + ] + }, + { + "cell_type": "markdown", + "id": "00aee50a-e3df-4511-bb56-b61ce5caf8f9", + "metadata": {}, + "source": [ + "In the next function, we use the property of the coefficients $c_{-n+N}=c_{-n}=c_{n}^*$. Thus, we have $c_0$ and the conjugated pairs\n", + "$c_n$, $c_{N-n}$. \n", + "\n", + "If $N$ is even, then $n$ can run from 1 to $N/2-1$. The last conjugated pair is $c_{N/2-1}, c_{N/2+1}$. Finally, we add $c_{N/2}$ which must be real. Thus, we find \n", + "\\begin{align}\n", + "x_k &= c_0+c_{N/2}\\exp(2\\pi ik\\frac{N/2}{N})+\\sum_{n=1}^{N/2-1}[c_n \\exp(2\\pi i k\\frac{n}{N}) + c_{N-n}\\exp(2\\pi i k\\frac{N-n}{N})] \\\\\n", + " &= c_0+c_{N/2}\\exp(\\pi ik)+\\sum_{n=1}^{N/2-1}[c_n \\exp(2\\pi i k\\frac{n}{N}) + c_{n}^*\\exp(-2\\pi i k\\frac{n}{N})] \\\\\n", + " &= c_0+c_{N/2}\\exp(\\pi ik)+2\\sum_{n=1}^{N/2-1}\\mathcal{Re}\\left[c_n \\exp(2\\pi i k\\frac{n}{N})\\right] \\,.\n", + "\\end{align}\n", + "\n", + "For odd $N$, we again have $c_0$ and the conjugated pairs $c_n$, $c_{N-n}$. Now, $n$ can run from 1 to $(N-1)/2$. The last conjugated pair is $c_{(N-1)/2}, c_{N-(N-1)/2}=c_{(N+1)/2}$. There is no coefficient for $n=N/2$ which should be added. Thus we get:\n", + "\\begin{align}\n", + "x_k = c_0+2\\sum_{n=1}^{(N-1)/2}\\mathcal{Re}\\left[c_n \\exp(2\\pi i k\\frac{n}{N})\\right] \\,.\n", + "\\end{align}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a64a1710-563d-4e3d-91c6-708dbc03318f", + "metadata": {}, + "outputs": [], + "source": [ + "def dft_synthesis(c, nc=0):\n", + " \"\"\"\n", + " Evaluate the Fourier series as described above\n", + " :param c: complex DFT coefficients\n", + " :param nc: number of conjugated pairs of coefficients to be used \n", + " (nc <= N/2-1 for even N or nc <= (N-1)/2 for odd N) \n", + " (default: 0 indicating all coefficients)\n", + " \"\"\"\n", + " ns = np.size(c)\n", + " \n", + " # initialize x_k with c_0\n", + " x = np.ones(ns)*np.real(c[0]) \n", + "\n", + " # distinguish even and odd N, we later sum up to nmax\n", + " if ns%2 == 0: \n", + " nmax = ns//2-1\n", + " else:\n", + " nmax = (ns-1)//2\n", + "\n", + " # if N is even and nc == 0, add coefficient c[N/2]*exp(i k pi)\n", + " # the exponential is 1 for even k and -1 for odd k\n", + " if ns%2 == 0 and nc == 0:\n", + " x = x+np.real(c[ns//2])*np.where(np.arange(0, ns)%2 == 0, 1, -1)\n", + " \n", + " # check input nc, reset to nmax if greater\n", + " if nc == 0: nc = nmax\n", + " if nc > nmax: nc = nmax\n", + "\n", + " # do synthesis\n", + " for n in range(1, nc+1):\n", + " for k in range(ns):\n", + " arg = +2*np.pi*1j*n*k/ns\n", + " x[k] = x[k]+2*np.real(c[n]*np.exp(arg))\n", + " return x" + ] + }, + { + "cell_type": "markdown", + "id": "1cd060d7-66bf-4b4d-8bc5-ff4a674ee7a3", + "metadata": {}, + "source": [ + "### The boxcar function\n", + "We now start from a boxcar function which we want to recover from its Fourier coefficients. We do this for an increasing number of coefficients to watch how the series converges." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6de62b07-ed02-4567-8f4f-7e73fb272313", + "metadata": {}, + "outputs": [], + "source": [ + "tmax = 1000.0\n", + "dt = 1.0\n", + "ns = int(tmax/dt)\n", + "period = ns*dt\n", + "t = dt*np.arange(0, ns)\n", + "boxcar = np.where(t < 0.3*tmax, 0, 1)*np.where(t > 0.5*tmax, 0, 1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66ffef33-dc8a-41ab-90ca-7749ead9947d", + "metadata": {}, + "outputs": [], + "source": [ + "fig1, ax1 = SetupFigure(10, 5, \"Time [s]\", \"boxcar\", \"Boxcar, DT=0.1, TMAX=1000\", 14)\n", + "ax1.plot(t, boxcar, color='blue', ls='-')" + ] + }, + { + "cell_type": "markdown", + "id": "d74a4190-29f5-41e9-9c02-567eabbb80c0", + "metadata": {}, + "source": [ + "### The Fourier coefficients of the boxcar function\n", + "We first compute the Fourier coefficients of the boxcar to later carry out an incremental reconstruction." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7fa0c640-b264-4039-bb69-2182e8bd39b8", + "metadata": {}, + "outputs": [], + "source": [ + "c = dft_coeff(boxcar)\n", + "f = np.linspace(0, (ns-1)/period, ns)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d1455fd2-71ed-4ba4-b74d-3e3a42f815ba", + "metadata": {}, + "outputs": [], + "source": [ + "nmax = 50\n", + "fig2, ax2 = SetupFigure(10, 5, \"Frequency [Hz]\", \"c_n\", \"Fourier coefficients of Boxcar, DT=1.0, TMAX=1000, Nmax={:d}\".format(nmax), 14)\n", + "#ax2.plot(f[0:nmax], np.absolute(c[0:nmax]), color='black', ls='-', label='abs')\n", + "ax2.plot(f[0:nmax], np.real(c[0:nmax]), color='red', ls='-', label='real')\n", + "ax2.plot(f[0:nmax], np.imag(c[0:nmax]), color='blue', ls='-', label='imag')\n", + "ax2.legend()" + ] + }, + { + "cell_type": "markdown", + "id": "18a7cd08-0e87-4d10-8dfc-edaa70dcb167", + "metadata": {}, + "source": [ + "### Incremental reconstruction of the boxcar function from the Fourier coefficients \n", + "Now we do the Fourier synthesis with an increasing number of coefficients. Watch how the reconstruction changes. What happens at the edges for high numbers of coefficients?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7f465dfe-3d22-45f5-8b15-c7bc5d59687c", + "metadata": {}, + "outputs": [], + "source": [ + "nc = 20\n", + "xr = dft_synthesis(c, nc)\n", + "fig3, ax3 = SetupFigure(10, 5, \"Time [s]\", \"xr(t)\", \"Recovered Boxcar, DT=1, TMAX=1000, Nc={:d}\".format(nc), 14)\n", + "ax3.plot(t, xr , color='blue', ls='-')\n", + "ax3.plot(t, boxcar, color='red', ls=':')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e670440-e383-4ae8-b13e-490e9e6dcde2", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}