From 3e7a875991249e67e963305856522d1429566efd Mon Sep 17 00:00:00 2001 From: friederichwo Date: Wed, 7 May 2025 11:13:09 +0200 Subject: [PATCH] added notebooks to assignment 2 --- assignment_02/dft.ipynb | 357 +++++++++++++++++++++++++++++++++++ assignment_02/setupFigure.py | 16 ++ 2 files changed, 373 insertions(+) create mode 100644 assignment_02/dft.ipynb create mode 100644 assignment_02/setupFigure.py diff --git a/assignment_02/dft.ipynb b/assignment_02/dft.ipynb new file mode 100644 index 0000000..a5bb35e --- /dev/null +++ b/assignment_02/dft.ipynb @@ -0,0 +1,357 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "c271abd4-6464-4d0d-a598-18623c080570", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib widget\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from setupFigure import SetupFigure" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "33012d14-c18a-4e61-89c5-10ac079f4121", + "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": "code", + "execution_count": null, + "id": "80d7b404-00c7-49e8-9061-a262e7b84bb7", + "metadata": {}, + "outputs": [], + "source": [ + "def dft_inverse(c):\n", + " \"\"\"\n", + " Evaluate x_k = sum_{n=0}^{N-1}c_n exp(i*2*pi*nk/N)\n", + " :param c: complex DFT coefficients\n", + " \"\"\"\n", + " ns = np.size(c)\n", + " x = np.zeros(ns) # zero values before summing\n", + " for k in range(ns):\n", + " for n in range(ns):\n", + " arg = +2*np.pi*1j*n*k/ns\n", + " x[k] = x[k]+np.real(c[n]*np.exp(arg)) # since x[k] is real, we only sum the real parts\n", + " return x" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b24ed02f-2388-444a-943c-3d6341f9490a", + "metadata": {}, + "outputs": [], + "source": [ + "def dft_inverse_any_time(c, period, dt, nt, mc=1):\n", + " \"\"\"\n", + " Evaluate x(k*dt) = sum_{n=0}^{N-1}c_n exp(i*2*pi*n*(k*dt)/T)\n", + " :param c: complex DFT coefficients\n", + " :param period: period of time series\n", + " :param dt: new time sampling interval\n", + " :param nt: number of time samples\n", + " :param mc: take mc*N Fourier coefficients\n", + " \"\"\"\n", + " ns = np.size(c)\n", + " x = np.zeros(nt)\n", + " for k in range(nt):\n", + " for n in range(mc*ns):\n", + " arg = +2*np.pi*1j*n*(k*dt)/period\n", + " dum, n2 = divmod(n,ns) # use periodicity c[n+N] = c[n]\n", + " x[k] = x[k]+np.real(c[n2]*np.exp(arg)) # since x[k] result is real, we only sum the real parts\n", + " return x" + ] + }, + { + "cell_type": "markdown", + "id": "4112da9b-5ff0-4306-a4a4-d7cf3924b274", + "metadata": {}, + "source": [ + "### The sampled Gaussian function" + ] + }, + { + "cell_type": "markdown", + "id": "941b218a-4bbf-4588-962d-177a300a1051", + "metadata": {}, + "source": [ + "Setup of Gaussian function $x(t)=\\exp(-(t-3)^2)$ in range $0 < t < 6$. We take $\\Delta t=0.1$ producing 60 intervals and N=61 samples.\n", + "Note that the period of the function is then $T=N\\Delta t=6.1$ and not $6.0$; hence $x(T=N\\Delta t) = x(0)$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bd7afa4c-1091-4b0d-a211-38933215d782", + "metadata": {}, + "outputs": [], + "source": [ + "tmax = 6.0\n", + "dt = 0.1\n", + "ns = int(tmax/dt)+1\n", + "period = ns*dt\n", + "x = np.zeros(ns)\n", + "t = np.linspace(0, tmax, ns)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "03f1b9cb-be96-4e54-a479-1a3ed2da3d0a", + "metadata": {}, + "outputs": [], + "source": [ + "x = np.exp(-(t-3)**2) # compute ns samples of Gaussian" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe9d7345-81b0-4c49-afab-8d15acf62bbe", + "metadata": {}, + "outputs": [], + "source": [ + "fig1, ax1 = SetupFigure(10, 5, \"Time [s]\", \"x(t)\", \"Gaussian, DT=0.1, TMAX=6\", 14)\n", + "ax1.plot(t, x , color='blue', ls=':', marker='o', markersize=4.0)" + ] + }, + { + "cell_type": "markdown", + "id": "ab863cff-159d-4586-925e-580d5aedcaaf", + "metadata": {}, + "source": [ + "### The Fourier coefficients of the sampled Gaussian" + ] + }, + { + "cell_type": "markdown", + "id": "1e2933b2-adbf-4c20-96a0-34d5c0bb0077", + "metadata": {}, + "source": [ + "Compute the Fourier coefficients using $c_n = \\frac{1}{N}\\sum\\limits_{k=0}^{N-1} x_k \\exp(-2\\pi i nk/N)$. The frequencies to $c_n$ are $f_n = n/T$ or $\\omega_n = 2\\pi n/T$. The frequency to $c_{N-1}$ is $\\frac{N-1}{T}$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "015ba673-f453-46c8-b034-03bf72ffe965", + "metadata": {}, + "outputs": [], + "source": [ + "c = dft_coeff(x)\n", + "f = np.linspace(0, (ns-1)/period, ns)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ecdfb4c-28c0-447c-8110-52b3c6148992", + "metadata": {}, + "outputs": [], + "source": [ + "fig2, ax2 = SetupFigure(10, 5, \"Frequency [Hz]\", \"c_n\", \"Fourier coefficients of Gaussian, DT=0.1, TMAX=6\", 14)\n", + "ax2.plot(f, np.absolute(c), color='black', ls='--', marker='d', label='abs')\n", + "ax2.plot(f, np.real(c), color='red', ls=':', marker='*', markersize=4.0, label='real')\n", + "ax2.plot(f, np.imag(c), color='blue', ls=':', marker='o', markersize=4.0, label='imag')\n", + "ax2.legend()" + ] + }, + { + "cell_type": "markdown", + "id": "d8e7ff1d-c43b-418b-9d04-e104fd127a3f", + "metadata": {}, + "source": [ + "### The recovered Gaussian at the original samples" + ] + }, + { + "cell_type": "markdown", + "id": "d9fbd31d-af26-4e7c-8db9-eba38cd2d9b0", + "metadata": {}, + "source": [ + "Recover the samples by inverse discrete Fourier transform: $x_k = \\sum\\limits_{n=0}^{N-1} c_n \\exp(2\\pi i nk/N)$ belonging to times $t_k=k\\Delta t$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad695630-ed05-4508-abb7-802b1ef0556f", + "metadata": {}, + "outputs": [], + "source": [ + "xr = dft_inverse(c)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf73ed2b-5571-4d8e-bfb1-aa917dc351bc", + "metadata": {}, + "outputs": [], + "source": [ + "fig3, ax3 = SetupFigure(10, 5, \"Time [s]\", \"xr(t)\", \"Recovered Gaussian, DT=0.1, TMAX=6\", 14)\n", + "ax3.plot(t, xr , color='blue', ls='--', marker='o', markersize=4.0)" + ] + }, + { + "cell_type": "markdown", + "id": "d741cc78-d48e-4ffd-9029-811367b4b972", + "metadata": {}, + "source": [ + "### The recovered Gaussian at original sampling interval up to TMAX=12" + ] + }, + { + "cell_type": "markdown", + "id": "5008030a-6f9f-4613-9937-ebbb6e0d9298", + "metadata": {}, + "source": [ + "Evaluate the inverse transform at times beyond $T$ using $x_k = \\sum\\limits_{n=0}^{N-1} c_n \\exp(2\\pi in t/T)$. Here, we take the same sampling interval but extend the time range to $12$. We get a periodic continuation of the Gaussian." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fcef58bc-051d-49e9-b37f-40b6bc91a87e", + "metadata": {}, + "outputs": [], + "source": [ + "tmax = 12\n", + "nt = int(tmax/dt)+1\n", + "tlong = np.linspace(0, tmax, nt)\n", + "xrlong = dft_inverse_any_time(c, period, dt, nt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8e527baa-c2ea-43a0-b10c-8366048d6eac", + "metadata": {}, + "outputs": [], + "source": [ + "fig4, ax4 = SetupFigure(10, 5, \"Time [s]\", \"xrlong(t)\", \"Recovered Gaussian, DT=0.1, TMAX=12\", 14)\n", + "ax4.plot(tlong, xrlong , color='blue', ls='--', marker='o', markersize=4.0)" + ] + }, + { + "cell_type": "markdown", + "id": "aa22d27c-a959-4e3f-b7bf-5fc1561d9956", + "metadata": {}, + "source": [ + "### The recovered Gausian with dense sampling (DT=0.01)" + ] + }, + { + "cell_type": "markdown", + "id": "08123a90-df5b-4bf8-aade-4819b55d73db", + "metadata": {}, + "source": [ + "Here, we also evaluate the Fourier series for times in between the original samples by choosing times $t=0.1\\,k\\Delta t$. Evidently, this is not a suitable interpolation method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10458503-2f03-4a16-b20e-8fb4aa0c99b0", + "metadata": {}, + "outputs": [], + "source": [ + "tmax = 6\n", + "dtnew = 0.1*dt\n", + "nt = int(tmax/dtnew)+1\n", + "tdense = np.linspace(0, tmax, nt)\n", + "xdense = dft_inverse_any_time(c, period, dtnew, nt, mc=1) # dense sampling" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2cc0d1c9-f2a6-4bfd-84d9-bd8af31818d0", + "metadata": {}, + "outputs": [], + "source": [ + "fig5, ax5 = SetupFigure(10, 5, \"Time [s]\", \"xrdense(t)\", \"Recovered Gaussian, DT=0.01, TMAX=6\", 14)\n", + "ax5.plot(tdense, xdense , color='blue', ls=':', marker='o', markersize=4.0)\n", + "ax5.plot(t, x , color='red', ls='--', marker='d', markersize=4.0)" + ] + }, + { + "cell_type": "markdown", + "id": "061d0254-a311-4a70-b26f-f51f1ae523ac", + "metadata": {}, + "source": [ + "### The recovered Gausian with dense sampling and high order Fourier coefficients (DT=0.01, mc*N-1)" + ] + }, + { + "cell_type": "markdown", + "id": "20d7b156-6a75-43e7-b5d1-04d9663b15be", + "metadata": {}, + "source": [ + " Here, we again use dense sampling and, in addition, sum up the coefficients to $mc*N-1$ thus introducing higher frequencies. By doing this, we aproximate the time function $x(t)=\\sum_{k=0}^{N-1} x_k\\delta(t-k\\Delta t)$ with high peaks at the samples and zero values in between. Play with increasingly higher values of $mc$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "902ecafc-1f7e-4b16-a518-19f1969e2990", + "metadata": {}, + "outputs": [], + "source": [ + "xhigh = dft_inverse_any_time(c, period, dtnew, nt, mc=60) # dense sampling and higher freqeuncies" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d3e681e-49e5-4549-b864-bbf6df6f61de", + "metadata": {}, + "outputs": [], + "source": [ + "fig6, ax6 = SetupFigure(10, 5, \"Time [s]\", \"xrhigh(t)\", \"Recovered Gaussian, DT=0.01, TMAX=6, high frequencies\", 14)\n", + "ax6.plot(tdense, xhigh , color='blue', ls=':', marker='o', markersize=4.0)" + ] + } + ], + "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 +} diff --git a/assignment_02/setupFigure.py b/assignment_02/setupFigure.py new file mode 100644 index 0000000..e53468c --- /dev/null +++ b/assignment_02/setupFigure.py @@ -0,0 +1,16 @@ +import matplotlib.pyplot as plt + + +def SetupFigure(wx, wy, xlabel, ylabel, title, tls): + """ + Create a new figure frame + """ + fig = plt.figure(figsize=(wx,wy)) + ax = fig.add_subplot(1,1,1) + ax.set_title(title, fontsize=tls+2) + ax.grid(which='major', axis='both', ls=':') + ax.set_ylabel(ylabel, fontsize=tls) + ax.set_xlabel(xlabel, fontsize=tls) + ax.tick_params(labelsize=tls) + return fig, ax +