added notebooks to assignment 2
This commit is contained in:
parent
eaa9ef9d5b
commit
3e7a875991
357
assignment_02/dft.ipynb
Normal file
357
assignment_02/dft.ipynb
Normal file
@ -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
|
||||
}
|
||||
16
assignment_02/setupFigure.py
Normal file
16
assignment_02/setupFigure.py
Normal file
@ -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
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user