new assignment 11

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 inline\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=10)         # 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
+}

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
+

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
+}

assignment_03/p3_gibbs_incremental_reconstruction.ipynb

@@ -0,0 +1,221 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "2ec965c5-047b-48c9-9ffd-60f096e06c1e",
+   "metadata": {},
+   "source": [
+    "## Assignment 3"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e504320f-93bf-4344-9529-05c51538b4bd",
+   "metadata": {},
+   "source": [
+    "### Incremental reconstruction of a time series and the Gibbs phenomenon\n",
+    "The idea of this assignment is to calculate Fourier coefficients of a discontinuous boxcar function and perform Fourier synthesis with an increasing number of Fourier coefficients. You will find that the reconstruction oscillates at the jumps. This behaviour is called the Gibbs phenomenon.\n",
+    "\n",
+    "After the function definitions you find several tasks to be carried out by you."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "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": 3,
+   "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 do Fourier synthesis. We use the property of the coefficients $c_{-n+N}=c_{-n}=c_{n}^*$. Thus, we have $c_0$ and the conjugated pairs $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": 4,
+   "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": [
+    "### Task 1: The boxcar function\n",
+    "Write code to calculate samples of the boxcar function with $\\Delta t = 1$:\n",
+    "\\begin{align}\n",
+    "b(t) = \\left\\{\\begin{array}{l} \n",
+    "0 \\quad\\mathrm{for}\\quad 0 \\le t \\lt 300 \\\\ \n",
+    "1 \\quad\\mathrm{for}\\quad 300 \\le t \\lt 500 \\\\ \n",
+    "0 \\quad\\mathrm{for}\\quad 500 \\le t \\lt 1000 \\end{array}\\right.\n",
+    "\\end{align}\n",
+    "Assume that this function continues periodically to smaller and greater times. Take samples from $0$ up to $t=1000-\\Delta t = 999$ giving you 1000 samples and a period of $T=1000$.\n",
+    "\n",
+    "Produce a graph of the boxcar function."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6de62b07-ed02-4567-8f4f-7e73fb272313",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "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": [
+    "### Task 2: The Fourier coefficients of the boxcar function\n",
+    "Compute the Fourier coefficients of the boxcar to later carry out an incremental reconstruction. Call the function dft_coeff and setup an array with the frequencies associated with the coefficients. Produce a graph of the first 50 coefficients and plot real and imaginary part."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "7fa0c640-b264-4039-bb69-2182e8bd39b8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "18a7cd08-0e87-4d10-8dfc-edaa70dcb167",
+   "metadata": {},
+   "source": [
+    "### Task 3: Incremental reconstruction of the boxcar function from the Fourier coefficients \n",
+    "Do the Fourier synthesis with an increasing number of coefficients. \n",
+    "\n",
+    "Call the function dft_synthesis which allow to sum up a prescribed number of conjugated coefficient pairs. Vary the number of coefficients and watch how the reconstruction changes. What happens at the edges of the boxcar for high numbers of coefficients?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7f465dfe-3d22-45f5-8b15-c7bc5d59687c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  }
+ ],
+ "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
+}

assignment_03/setupFigure.py

@@ -0,0 +1,15 @@
+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

assignment_04/dftSlow.py

@@ -0,0 +1,52 @@
+import numpy as np
+
+
+def dft_coeff(x):
+    """
+    Evaluate c_n = 1/N*sum_{k=0}^{N-1}x_k exp(-i*2*pi*nk/N)
+    :param x: array of samples of function
+    """
+    ns = np.size(x)
+    c = np.zeros(ns, dtype=np.complex64)  # zero coefficients before summing
+    for n in range(ns):
+        for k in range(ns):
+            arg = -2*np.pi*1j*n*k/ns
+            c[n] = c[n]+x[k]*np.exp(arg)
+        c[n] = c[n]/ns
+    return c
+
+
+def dft_synthesis(c, nc=0):
+    """
+    Evaluate the Fourier series as described above
+    :param c: complex DFT coefficients
+    :param nc: number of conjugated pairs of coefficients to be used
+               (nc <= N/2-1 for even N or nc <= (N-1)/2 for odd N)
+               (default: 0 indicating all coefficients)
+    """
+    ns = np.size(c)
+
+    # initialize x_k with c_0
+    x = np.ones(ns) * np.real(c[0])
+
+    # distinguish even and odd N, we later sum up to nmax
+    if ns % 2 == 0:
+        nmax = ns // 2 - 1
+    else:
+        nmax = (ns - 1) // 2
+
+    # if N is even and nc == 0, add coefficient c[N/2]*exp(i k pi)
+    # the exponential is 1 for even k and -1 for odd k
+    if ns % 2 == 0 and nc == 0:
+        x = x + np.real(c[ns // 2]) * np.where(np.arange(0, ns) % 2 == 0, 1, -1)
+
+    # check input nc, reset to nmax if greater
+    if nc == 0: nc = nmax
+    if nc > nmax: nc = nmax
+
+    # do synthesis
+    for n in range(1, nc + 1):
+        for k in range(ns):
+            arg = +2 * np.pi * 1j * n * k / ns
+            x[k] = x[k] + 2 * np.real(c[n] * np.exp(arg))
+    return x

assignment_04/p4_fft_sampling_interpol.ipynb

@@ -0,0 +1,220 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1559cf77-eb5b-4e56-8bd9-ae61a108fe92",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "from setupFigure import SetupFigure\n",
+    "from dftSlow import dft_coeff, dft_synthesis"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b101d00c-d9a2-4222-879b-ac340f89f576",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def dft_fast_coeff(x):\n",
+    "    \"\"\"\n",
+    "    Evaluate Fourier coefficients using Numpy's fast Fourier transform.\n",
+    "    This routine only returns the coefficients for the positive frequencies.\n",
+    "    If N is even, it goes up to n=N/2.\n",
+    "    If N is odd, it goes up to n=(N-1)/2.\n",
+    "    :param x: array of function samples\n",
+    "    \"\"\"\n",
+    "    return np.fft.rfft(x, norm='forward')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "14467eda-e7f2-423d-9b78-b7cc9672f22c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def dft_fast_synthesis(fc, outnum='even'):\n",
+    "    \"\"\"\n",
+    "    Use numpy's fast Fourier synthesis taking only Fourier coefficients for positive frequencies\n",
+    "    :param fc: aray of coefficients for positive frequencies only.\n",
+    "    :param outnum: specifies if output time series has an even or odd number of samples (default: 'even')\n",
+    "    \"\"\"\n",
+    "    ns = 2*fc.size-2\n",
+    "    if outnum == 'odd': ns = 2*fc.size-1\n",
+    "    return np.fft.irfft(fc, ns, norm='forward')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "18ffbc73-5210-4b9a-b3e2-28d5aafd33d9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def boxcar(dt, period, tup, tdown):\n",
+    "    \"\"\"\n",
+    "    Calculate samples of a boxcar function\n",
+    "    :param dt: sampling interval\n",
+    "    :param period: time range is 0 <= t < period (multiple of dt)\n",
+    "    :param tup: time where boxcar goes from 0 to 1\n",
+    "    :param tdown: time where boxcar goes from 1 to 0\n",
+    "    \"\"\"\n",
+    "    ns = int(period/dt)\n",
+    "    t = dt*np.arange(0, ns)\n",
+    "    return t, np.where(t < tup, 0, 1)*np.where(t > tdown, 0, 1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3435e191-57a4-425a-9618-a7597b11916d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def gaussian(dt, period, tmax, hwidth):\n",
+    "    \"\"\"\n",
+    "    Calculate samples of a Gaussian function\n",
+    "    :param dt: sampling interval\n",
+    "    :param period: time range is 0 <= t < period (multiple of dt)\n",
+    "    :param tmax: time of maximum of Gaussian\n",
+    "    :param hwidth: half width of Gaussian\n",
+    "    \"\"\"\n",
+    "    ns = int(period/dt)\n",
+    "    t = dt*np.arange(0, ns)\n",
+    "    return t, np.exp(-(t-tmax)**2/hwidth**2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8a7e3fe9-484a-49cc-af35-03b532eb62cd",
+   "metadata": {},
+   "source": [
+    "## Task 1: Compare \"slow\" Fourier transform with the fast version of numpy"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "971cc8b7-7ade-4755-b5ad-9cc83e66fed5",
+   "metadata": {},
+   "source": [
+    "Again set up the boxcar function as in the previous assignment and use the provided functions to compute the Fourier coefficients by both methods. Verify that both methods yield the same results by printing the first 20 coefficients."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f5da468c-fded-40b4-a691-6adc93fbf110",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bc4a9050-a171-479e-90a5-df7b924143a9",
+   "metadata": {},
+   "source": [
+    "Also compare the slow and fast versions of Fourier synthesis by calling the provided functions. Print values at times around the discontinuities."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a8abee97-2424-445d-81c3-98353d5ac270",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1204f35f-7933-490e-94a7-af77edaca83c",
+   "metadata": {},
+   "source": [
+    "## Task 2: Interpolation by appending zeros to the Fourier coefficients"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "accff2d0-77ce-4048-bad0-8396c62ef675",
+   "metadata": {},
+   "source": [
+    "A band limited time series can be interpolated using a simple trick: First, the Fourier coefficients are computed up to the Nyquist frequency, $f_{Ny} = \\frac{N}{2}\\Delta f = \\frac{N}{2T}$. Then, $L$ zero coefficents are appended to increase the Nyquist frequency to $f'_{Ny} = (\\frac{N}{2}+L)\\Delta f$ and to decrease the sampling interval to $\\Delta t' = \\frac{1}{2f'_{Ny}}$. Subsequent Fourier synthesis produces an interpolated version of the original time series. These relations hold for even and odd number of samples.\n",
+    "When doing the Fourier synthesis, the routine should be called with outnum='odd' for odd N and with outnum='even' for even N, respectively.\n",
+    "\n",
+    "First set up a Gaussian using the provided function. Then compute the Fourier coefficients. Print out the number of samples, the Nyquist frequency and the number of Fourier coefficients. For example, choose dt=5, period=100, tmax=50, hwidth=20. \n",
+    "\n",
+    "Second, append some zeros (e.g. 20) to the array of coefficients and compute the new Nyquist frequency and the new sampling interval. Do the Fourier synthesis and print the new number of samples, the new Nyquist frequency and the new sampling interval.\n",
+    "\n",
+    "Third, plot the new and old time series into one graph."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7c5663c1-5605-45df-8342-2dcefb2ec4b0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "fc9fd425-811c-49d1-8372-1fb2963a8310",
+   "metadata": {},
+   "source": [
+    "## Task 3: Aliasing and the sampling theorem"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a1cf5112-c521-46d9-b642-4c3fbe3fa0ee",
+   "metadata": {},
+   "source": [
+    "First, calculate values for a Gaussian with dt=0.25, period=100, tmax=50 and hwidth=1. Plot the Gaussian.\n",
+    "\n",
+    "Second, in a loop, calculate the same Gaussian with dt = 0.5, 1.0, 1.5 and 2.0. Compute the fast Fourier coefficients and the frequencies associated with them. Plot the absolute value of the coefficients into one graph. Set the upper frequency axis limit to 1.0 and use different colors for the curves. Compare the spectra, what do you observe?"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9e8596d7-1eb5-4248-85c5-e90ed8beaeb8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  }
+ ],
+ "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
+}

assignment_04/setupFigure.py

@@ -0,0 +1,15 @@
+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

assignment_05/acausalFilters.py

@@ -0,0 +1,34 @@
+# -------------------------------------------------------------
+#  Functions for acausal filtering
+# -------------------------------------------------------------
+import numpy as np
+
+def boxcar_lowpass_filter(fc, df, nf):
+    """
+    Calculates samples of the boxcar lowpass filter transfer function (positive frequencies only)
+    :param fc: corner frequency in Hz
+    :param df: frequency stepping
+    :param nf: number of frequencies
+    """
+    return np.where(df*np.arange(0, nf) > fc, 0, 1)
+
+
+def boxcar_highpass_filter(fc, df, nf):
+    """
+    Calculates samples of the boxcar highpass filter transfer function (positive frequencies only)
+    :param fc: corner frequency in Hz
+    :param df: frequency stepping
+    :param nf: number of frequencies
+    """
+    return np.where(df*np.arange(0, nf) < fc, 0, 1)
+
+
+def hann_lowpass_filter(fc, df, nf):
+    """
+    Calculates samples of the Hann filter transfer function (positive frequencies only)
+    :param fc: corner frequency in Hz
+    :param df: frequency stepping
+    :param nf: number of frequencies
+    """
+    f = df*np.arange(0, nf)
+    return np.where(f < fc, np.cos(np.pi*f/fc)**2, 0)

assignment_05/dftFast.py

@@ -0,0 +1,27 @@
+# -------------------------------------------------------------
+#  Functions for "fast" Fourier transform and Fourier synthesis
+#  using numpy tools
+# -------------------------------------------------------------
+import numpy as np
+
+def dft_fast_coeff(x):
+    """
+    Evaluate Fourier coefficients using Numpy's fast Fourier transform.
+    This routine only returns the coefficients for the positive frequencies.
+    If N is even, it goes up to n=N/2.
+    If N is odd, it goes up to n=(N-1)/2.
+    :param x: array of function samples
+    """
+    return np.fft.rfft(x, norm='forward')
+
+
+def dft_fast_synthesis(fc, outnum='even'):
+    """
+    Use numpy's fast Fourier synthesis taking only Fourier coefficients for positive frequencies
+    :param fc: aray of coefficients for positive frequencies only.
+    :param outnum: specifies if output time series has an even or odd number of samples (default: 'even')
+    """
+    ns = 2*fc.size-2
+    if outnum == 'odd': ns = 2*fc.size-1
+    return np.fft.irfft(fc, ns, norm='forward')
+

assignment_05/p5-spectral-estimation.ipynb

@@ -0,0 +1,358 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0ff7abb1-c698-480b-a894-bdbb5d151fe6",
+   "metadata": {},
+   "source": [
+    "## Assignment 5\n",
+    "In this assignment, we look at spectral estimation and study the effects of acausal filters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "156c8712-11be-4721-9d39-ff2c47a66bcb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from setupFigure import SetupFigure\n",
+    "from dftFast import dft_fast_coeff, dft_fast_synthesis\n",
+    "from seismicTrace import SeismicTrace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f98951f5-96dc-4ef5-b365-85bffa21f700",
+   "metadata": {},
+   "source": [
+    "### Task 1: Functions for Hann window and filter transfer functions\n",
+    "#### 1.1 Write a function that calculates samples of a Hann window for given sampling interval and window length."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c36eb15e-26e7-4e6e-bef9-8d566dc07411",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def hann_window(dt, tlen):\n",
+    "    \"\"\"\n",
+    "    Calculate samples of the Hann window for given DT and TLEN.\n",
+    "    Number of samples is assumed to be int(tlen/dt)\n",
+    "    :param dt: sampling interval\n",
+    "    :param tlen: length of window\n",
+    "    \"\"\"\n",
+    "    # Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "60dffc7b-17a1-40aa-8a6f-918504e06af2",
+   "metadata": {},
+   "source": [
+    "#### 1.2 Write a function that calculates samples of the boxcar lowpass filter transfer function for positive frequencies only."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "be77830b-6b2f-4371-9711-e83a1598c363",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def boxcar_lowpass_filter(fc, df, nf):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the boxcar lowpass filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    \"\"\"\n",
+    "    # Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "33eb7b1d-4abd-42a2-99ec-cdd1765fbfde",
+   "metadata": {},
+   "source": [
+    "#### 1.3 Write a function that calculates samples of the boxcar highpass filter transfer function for positive frequencies only."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29c56269-025a-43e9-8748-c0514a72f70d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def boxcar_highpass_filter(fc, df, nf):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the boxcar highpass filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    \"\"\"\n",
+    "    # Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "937b6b14-4d8f-40ed-bc71-86888d6e726c",
+   "metadata": {},
+   "source": [
+    "#### 1.4 Write a function that calculates samples of the Hann filter transfer function for positive frequencies only."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0101ff10-d14d-4b96-83f8-b26f96e5ec8d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def hann_lowpass_filter(fc, df, nf):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the Hann filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    \"\"\"\n",
+    "    # Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d641182d-5526-45a9-a69a-0a6965fd942b",
+   "metadata": {},
+   "source": [
+    "### Task 2: Spectral estimation and windowing\n",
+    "Here, we first setup a time series that consists of a superposition of different exponentially damped sine functions simulating a series of free oscillations:\n",
+    "\\begin{align}\n",
+    "s(t) = \\sum_{k=0}^N A_k\\sin(2\\pi f_k t)\\,\\exp(-\\pi f_k t/Q_k) \\,.\n",
+    "\\end{align}\n",
+    "The aim is to find the frequencies by Fourier analysis."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eba91f73-1651-4fd7-ad90-70980b1c5d4d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "feig = np.array([0.3, 0.48, 0.64, 0.68, 0.82, 0.85, 0.93, 0.94])*1.e-3         # eigenfreqeuncies in Hz\n",
+    "amp =  np.array([1.0, 0.30, 2.00, 0.61, 4.00, 0.53, 3.00, 0.14])               # amplitudes\n",
+    "phi =  np.array([ 11,   46,  271,  123,  337,  173,   65,  221])*np.pi/180.    # phases\n",
+    "q =    np.array([600,  200, 2000, 1000,  600,  100,  800,  900])               # quality factors"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e455172a-50ad-4aa4-a1b3-d29d50a52b49",
+   "metadata": {},
+   "source": [
+    "#### 2.1 Write code that implements the above equation using the provided frequencies, amplitudes, phases and Q values. Choose dt=100 s and allow a varaible time series length of tens of hours. Plot the resulting time series."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c808e6c4-8e39-4e7c-a689-860230ff775a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7ae65e6e-0f73-42fc-b251-878f0ef7f36c",
+   "metadata": {},
+   "source": [
+    "#### 2.2 Produce a second time series by multiplying the original one with a Hann window and also plot it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b62c1a43-f57e-4397-a8f9-7ec813db1e9d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "26feab1b-41a6-464f-b578-95814286b915",
+   "metadata": {},
+   "source": [
+    "#### 2.3 Perform a fast Fourier analysis, calculate the associated frequencies and plot the resulting amplitude spectra of both the original and the windowed time series into one graph. Compare the results. Try different time series lengths (between 20 and 200 hours)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "164cd4c2-72ea-44e8-a236-95475118aad1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "689b3af4-1146-4041-9f4e-baed098d2a23",
+   "metadata": {},
+   "source": [
+    "### Task 3: Read a seismogram from file and compute spectrum\n",
+    "Here, we first read a recorded seismogram form a file, do a Fourier analysis and later apply diverse acausal filters. We make use of the class \"SeismicTrace\" that defines a data structure for a seismic trace and also some useful functions related to traces. You may want to look at the provided code in \"seismicTrace.py\"."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "912ae1f7-8cbc-4446-9ef0-b446b57d789d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "seisfile = \"CH.PANIX.HHZ\"\n",
+    "with open(seisfile, 'r') as fd:\n",
+    "    seistrace = SeismicTrace.readFromAscii(fd)\n",
+    "t = seistrace.tanf+seistrace.dt*np.arange(0, seistrace.nsamp)\n",
+    "seistrace.printInfo()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eaee9262-f577-49a1-acee-9b9c27f996be",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig4, ax4 = SetupFigure(12, 5, \"Time [s]\", \"s(t)\", \"Unfiltered seismic data\", 14)\n",
+    "ax4.plot(t, seistrace.vals, color='blue', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "934a4fa1-7f8d-4977-bf86-bfe094abc5e5",
+   "metadata": {},
+   "source": [
+    "#### 3.1 Perform a fast Fourier analysis, calculate the associated frequencies and plot the resulting amplitude spectrum."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c4a0c8c9-b51a-4116-9988-73d75b2f6704",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spec = dft_fast_coeff(seistrace.vals)\n",
+    "f = 1./seistrace.tlen*np.arange(0, spec.size)\n",
+    "print(\"Number of frequencies: \", spec.size, \" fmax = \", f[-1], \"Hz\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "82fb96d5-39ac-4c6a-a360-3544167ec26b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig5, ax5 = SetupFigure(12, 5, \"Frequency [Hz]\", \"S(f)\", \"Amplitude spectrum of seismogram\", 14)\n",
+    "ax5.set_xlim(0, 1)\n",
+    "ax5.plot(f, np.absolute(spec), color='blue', ls='-', label='boxcar')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8d36d0c3-6faf-4e0b-b934-8822100b2b2f",
+   "metadata": {},
+   "source": [
+    "### Task 4: Acausal low pass filtering"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "36b37ff1-3080-4651-a344-9cb4ede9109d",
+   "metadata": {},
+   "source": [
+    "Apply the boxcar lowpass filter to the data and plot the result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1e59a76-d251-4d66-b6ba-7c29a0626b3e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e1097e6d-c459-489d-aecd-75911419a319",
+   "metadata": {},
+   "source": [
+    "#### 4.1 Apply the Hann low pass filter to the data and plot the result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "be2e1d6d-0c4e-4bfa-8612-ce626393bdff",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Here comes your code"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3bbec473-2bf7-41d0-828d-9413dff147e8",
+   "metadata": {},
+   "source": [
+    "#### 4.2 Apply the boxcar highpass filter to the data and plot the result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c9ec91f3-cc1c-48e5-ac53-882002fc74f6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Here comes your code"
+   ]
+  }
+ ],
+ "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
+}

assignment_05/seismicTrace.py

@@ -0,0 +1,79 @@
+#----------------------------------------------------------------------------
+#  Python class for a seismic trace
+#
+import numpy as np
+
+class SeismicTrace(object):
+    """
+    Class defining a seismic trace by station name, network name, component,
+    number of samples, sampling interval and start time
+    """
+
+    def __init__(self, staname, comp, dt, tanf, vals):
+        """
+        Initialize seismic trace object from data:
+        :param staname: station name
+        :param comp: component
+        :param dt: sampling interval
+        :param tanf: start time
+        :param vals: values of trace (1D numpy array)
+        """
+        self.staname = staname
+        self.comp = comp
+        self.dt = dt
+        self.tanf = tanf
+        self.vals = vals
+        self.nsamp = vals.size
+        self.tlen = self.dt*self.nsamp
+
+
+    @classmethod
+    def readFromAscii(cls, fid, comp='None'):
+        """
+        Class method as second constructor to initialize trace from an Ascii file
+        :param fid: file descriptor of already open Ascii file
+        :param comp: optional component selection
+        :return: an instance of the class
+        """
+        staname = fid.readline().strip()
+        if not staname:
+            print("No station name found. Is file really an Ascii gather file")
+            return None
+        compf = fid.readline().strip()
+        if comp != 'None' and comp != compf:
+            print("Requested component %s not found.", comp)
+            return None
+        h = fid.readline().strip().split()
+        vals = np.array([float(v) for v in fid.readline().strip().split()])
+        dt = float(h[1]); tanf = float(h[2])
+        return cls(staname, compf, dt, tanf, vals)
+
+
+    def absoluteMaximum(self):
+        """
+        return absolute maximum of trace
+        """
+        return abs(max(self.vals))
+
+
+    def writeToOpenAsciiFile(self, fid):
+        """
+        write a seismic trace to an open Ascii file
+        """
+        fid.write(self.staname + '\n')
+        fid.write(self.comp + '\n')
+        fid.write("{:12d} {:15.9f} {:20.9f}\n".format(self.nsamp, self.dt, self.tanf))
+        for n in range(self.nsamp):
+            fid.write("{:17.8e}".format(self.vals[n]))
+        fid.write("\n")
+
+
+    def printInfo(self):
+        """
+        Print out some information about the trace
+        """
+        print("Identification: ", self.staname, self.comp)
+        print("Sampling interval: ", self.dt)
+        print("Number of samples: ", self.nsamp)
+        print("Start time after midnight: ", self.tanf)
+        print("Length of trace: (N*DT) ", self.tlen)

assignment_05/setupFigure.py

@@ -0,0 +1,15 @@
+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

assignment_05/spectral-estimation.ipynb

@@ -0,0 +1,466 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "0ff7abb1-c698-480b-a894-bdbb5d151fe6",
+   "metadata": {},
+   "source": [
+    "## Assignment 5\n",
+    "In this assignment, we look at spectral estimation and study the effects of acausal filters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "156c8712-11be-4721-9d39-ff2c47a66bcb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from setupFigure import SetupFigure\n",
+    "from dftFast import dft_fast_coeff, dft_fast_synthesis\n",
+    "from seismicTrace import SeismicTrace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f98951f5-96dc-4ef5-b365-85bffa21f700",
+   "metadata": {},
+   "source": [
+    "### Task 1: Functions for Hann window and filter transfer functions\n",
+    "Write a function that calculates samples of a Hann window for given sampling interval and window length."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c36eb15e-26e7-4e6e-bef9-8d566dc07411",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def hann_window(dt, tlen):\n",
+    "    \"\"\"\n",
+    "    Calculate samples of the Hann window for given DT and TLEN.\n",
+    "    Number of samples is assumed to be int(tlen/dt)\n",
+    "    :param dt: sampling interval\n",
+    "    :param tlen: length of window\n",
+    "    \"\"\"\n",
+    "    ns = int(tlen/dt)\n",
+    "    return 2*np.sin(np.pi*dt*np.arange(0, ns)/tlen)**2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "60dffc7b-17a1-40aa-8a6f-918504e06af2",
+   "metadata": {},
+   "source": [
+    "Write a function that calculates samples of the boxcar lowpass filter transfer function for positive frequencies only."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "be77830b-6b2f-4371-9711-e83a1598c363",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def boxcar_lowpass_filter(fc, df, nf):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the boxcar lowpass filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    \"\"\"\n",
+    "    return np.where(df*np.arange(0, nf) > fc, 0, 1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "33eb7b1d-4abd-42a2-99ec-cdd1765fbfde",
+   "metadata": {},
+   "source": [
+    "Write a function that calculates samples of the boxcar highpass filter transfer function for positive frequencies only."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "29c56269-025a-43e9-8748-c0514a72f70d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def boxcar_highpass_filter(fc, df, nf):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the boxcar highpass filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    \"\"\"\n",
+    "    return np.where(df*np.arange(0, nf) < fc, 0, 1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "937b6b14-4d8f-40ed-bc71-86888d6e726c",
+   "metadata": {},
+   "source": [
+    "Write a function that calculates samples of the Hann filter transfer function for positive frequencies only."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0101ff10-d14d-4b96-83f8-b26f96e5ec8d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def hann_lowpass_filter(fc, df, nf):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the Hann filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    \"\"\"\n",
+    "    f = df*np.arange(0, nf)\n",
+    "    return np.where(f < fc, np.cos(np.pi*f/fc)**2, 0)    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d641182d-5526-45a9-a69a-0a6965fd942b",
+   "metadata": {},
+   "source": [
+    "### Task 2: Spectral estimation and windowing\n",
+    "Here, we first setup a time series that consists of a superposition of different exponentially damped sine functions simulating a series of free oscillations:\n",
+    "\\begin{align}\n",
+    "s(t) = \\sum_{k=0}^N A_k\\sin(2\\pi f_k t)\\,\\exp(-\\pi f_k t/Q_k) \\,.\n",
+    "\\end{align}\n",
+    "The aim is to find the frequencies by Fourier analysis."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eba91f73-1651-4fd7-ad90-70980b1c5d4d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "feig = np.array([0.3, 0.48, 0.64, 0.68, 0.82, 0.85, 0.93, 0.94])*1.e-3         # eigenfreqeuncies in Hz\n",
+    "amp =  np.array([1.0, 0.30, 2.00, 0.61, 4.00, 0.53, 3.00, 0.24])               # amplitudes\n",
+    "phi =  np.array([ 11,   46,  271,  123,  337,  173,   65,  221])*np.pi/180.    # phases\n",
+    "q =    np.array([600,  200, 2000, 1000,  600,  100,  800,  900])               # quality factors"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e455172a-50ad-4aa4-a1b3-d29d50a52b49",
+   "metadata": {},
+   "source": [
+    "Write code that impements the above equation using the provided frequencies, amplitudes, phases and Q values. Choose dt=100 s and allow a varaible time series length of tens of hours. Plot the resulting time series."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c808e6c4-8e39-4e7c-a689-860230ff775a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nhour = 100\n",
+    "dt = 100.0\n",
+    "tlen = 3600*nhour\n",
+    "ns = int(tlen/dt)\n",
+    "print(\"Number of samples: \", ns)\n",
+    "t = dt*np.arange(0, ns)\n",
+    "seis = np.zeros(ns)\n",
+    "for j, f in enumerate(feig):\n",
+    "    seis = seis + amp[j]*np.sin(2*np.pi*f*t+phi[j])*np.exp(-np.pi*f*t/q[j])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3d0dbe36-34e8-486d-bc63-f5be4aff2d0c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig1, ax1 = SetupFigure(12, 5, \"Time [h|\", \"s(t)\", \"Ground motion\", 14)\n",
+    "ax1.plot(t/3600, seis, color='blue', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7ae65e6e-0f73-42fc-b251-878f0ef7f36c",
+   "metadata": {},
+   "source": [
+    "Produce a second time series by multiplying the original one with a Hann window and also plot it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b62c1a43-f57e-4397-a8f9-7ec813db1e9d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "win = hann_window(dt, tlen)\n",
+    "seisw = seis*win"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "41193b3d-110b-42a4-b2ee-9ee0dbd759af",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig2, ax2 = SetupFigure(12, 5, \"Time [h|\", \"swin(t)\", \"Windowed ground motion\", 14)\n",
+    "ax2.plot(t/3600, seisw, color='blue', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "26feab1b-41a6-464f-b578-95814286b915",
+   "metadata": {},
+   "source": [
+    "Perform a fast Fourier analysis, calculate the associated frequencies and plot the resulting amplitude spectra of both the original and the windowed time series. Compare the results."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "164cd4c2-72ea-44e8-a236-95475118aad1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spec = dft_fast_coeff(seis)\n",
+    "f = 1./tlen*np.arange(0, spec.size)\n",
+    "print(\"Number of frequencies: \", spec.size, \" fmax = \", f[-1]*1000, \"mHz\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0c88a619-8b9c-44ed-bc5f-f6ecce938af3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "specw = dft_fast_coeff(seisw)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "73b8fae2-e712-47d7-99a5-ebe5375084cf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig3, ax3 = SetupFigure(12, 5, \"Frequency [mHz]\", \"S(f)\", \"Amplitude spectrum of ground motion\", 14)\n",
+    "ax3.set_xlim(0,1.2)\n",
+    "ax3.plot(f*1000, np.absolute(spec), color='blue', ls='-', label='boxcar')\n",
+    "ax3.plot(f*1000, np.absolute(specw), color='red', ls='-', label='hann')\n",
+    "ax3.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "689b3af4-1146-4041-9f4e-baed098d2a23",
+   "metadata": {},
+   "source": [
+    "### Task 2: Read a seismogram from file and compute spectrum\n",
+    "Here, we first read a recorded seismogram form a file, do a Fourier analysis and later apply diverse acausal filters. We make use of the class \"SeismicTrace\" that defines a data structure for a seismic trace and also some useful functions related to traces. You may want to look at the provided code in \"seismicTrace.py\"."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "912ae1f7-8cbc-4446-9ef0-b446b57d789d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "seisfile = \"CH.PANIX.HHZ\"\n",
+    "with open(seisfile, 'r') as fd:\n",
+    "    seistrace = SeismicTrace.readFromAscii(fd)\n",
+    "t = seistrace.tanf+seistrace.dt*np.arange(0, seistrace.nsamp)\n",
+    "seistrace.printInfo()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "eaee9262-f577-49a1-acee-9b9c27f996be",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig4, ax4 = SetupFigure(12, 5, \"Time [s]\", \"s(t)\", \"Unfiltered seismic data\", 14)\n",
+    "ax4.plot(t, seistrace.vals, color='blue', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "934a4fa1-7f8d-4977-bf86-bfe094abc5e5",
+   "metadata": {},
+   "source": [
+    "Perform a fast Fourier analysis, calculate the associated frequencies and plot the resulting amplitude spectrum."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c4a0c8c9-b51a-4116-9988-73d75b2f6704",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "spec = dft_fast_coeff(seistrace.vals)\n",
+    "f = 1./seistrace.tlen*np.arange(0, spec.size)\n",
+    "print(\"Number of frequencies: \", spec.size, \" fmax = \", f[-1], \"Hz\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "82fb96d5-39ac-4c6a-a360-3544167ec26b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig5, ax5 = SetupFigure(12, 5, \"Frequency [Hz]\", \"S(f)\", \"Amplitude spectrum of seismogram\", 14)\n",
+    "ax5.set_xlim(0, 1)\n",
+    "ax5.plot(f, np.absolute(spec), color='blue', ls='-', label='boxcar')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "8d36d0c3-6faf-4e0b-b934-8822100b2b2f",
+   "metadata": {},
+   "source": [
+    "### Task 3: Acausal low pass filtering"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "36b37ff1-3080-4651-a344-9cb4ede9109d",
+   "metadata": {},
+   "source": [
+    "Apply the boxcar lowpass filter to the data and plot the result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a1e59a76-d251-4d66-b6ba-7c29a0626b3e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fc = 0.1\n",
+    "df = 1./seistrace.tlen\n",
+    "specfil = spec*boxcar_lowpass_filter(fc, df, spec.size)\n",
+    "seisfil = dft_fast_synthesis(specfil)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "532ad374-2728-435d-8dca-afd7fcf9c475",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig6, ax6 = SetupFigure(12, 5, \"Time [s]\", \"sf(t)\", \"Boxcar lowpass filtered seismic data, FC = {:5.2f} Hz\".format(fc), 14)\n",
+    "#ax6.set_xlim(64000, 64500)\n",
+    "ax6.plot(t, seisfil, color='blue', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e1097e6d-c459-489d-aecd-75911419a319",
+   "metadata": {},
+   "source": [
+    "Apply the Hann low pass filter to the data and plot the result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "be2e1d6d-0c4e-4bfa-8612-ce626393bdff",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "specfil2 = spec*hann_lowpass_filter(fc, df, spec.size)\n",
+    "seisfil2 = dft_fast_synthesis(specfil2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b5d4bf72-9bc2-4679-9178-4f7bc650e4ee",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig7, ax7 = SetupFigure(12, 5, \"Time [s]\", \"sf(t)\", \"Hann lowpass filtered seismic data, FC = {:5.2f}\".format(fc), 14)\n",
+    "#ax7.set_xlim(64000, 64500)\n",
+    "ax7.plot(t, seisfil2, color='blue', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3bbec473-2bf7-41d0-828d-9413dff147e8",
+   "metadata": {},
+   "source": [
+    "Apply the boxcar highpass filter to the data and plot the result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c9ec91f3-cc1c-48e5-ac53-882002fc74f6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "specfil3 = spec*boxcar_highpass_filter(fc, df, spec.size)\n",
+    "seisfil3 = dft_fast_synthesis(specfil3)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7a24950e-8438-4e3b-903e-75f2f4408105",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig8, ax8 = SetupFigure(12, 5, \"Time [s]\", \"sf(t)\", \"Boxcar highpass filtered seismic data, FC = {:5.2f}\".format(fc), 14)\n",
+    "#ax7.set_xlim(64000, 64500)\n",
+    "ax8.plot(t, seisfil3, color='blue', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "92e0ebc5-9f7d-4302-8ad1-9aed6c7f2d6a",
+   "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
+}

assignment_06/dftFast.py

@@ -0,0 +1,27 @@
+# -------------------------------------------------------------
+#  Functions for "fast" Fourier transform and Fourier synthesis
+#  using numpy tools
+# -------------------------------------------------------------
+import numpy as np
+
+def dft_fast_coeff(x):
+    """
+    Evaluate Fourier coefficients using Numpy's fast Fourier transform.
+    This routine only returns the coefficients for the positive frequencies.
+    If N is even, it goes up to n=N/2.
+    If N is odd, it goes up to n=(N-1)/2.
+    :param x: array of function samples
+    """
+    return np.fft.rfft(x, norm='forward')
+
+
+def dft_fast_synthesis(fc, outnum='even'):
+    """
+    Use numpy's fast Fourier synthesis taking only Fourier coefficients for positive frequencies
+    :param fc: aray of coefficients for positive frequencies only.
+    :param outnum: specifies if output time series has an even or odd number of samples (default: 'even')
+    """
+    ns = 2*fc.size-2
+    if outnum == 'odd': ns = 2*fc.size-1
+    return np.fft.irfft(fc, ns, norm='forward')
+

assignment_06/p6-causal_filters.ipynb

@@ -0,0 +1,284 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "d4e32124-6e9e-4467-8726-d2adbb3934dd",
+   "metadata": {},
+   "source": [
+    "## Assignment 6\n",
+    "Here we apply causal filters to a seismological data trace."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c5711338-01f0-4e88-996d-b6026eff5098",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from setupFigure import SetupFigure\n",
+    "from dftFast import dft_fast_coeff, dft_fast_synthesis\n",
+    "from seismicTrace import SeismicTrace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bff77f1c-1f8a-4908-8419-dfed94f9f177",
+   "metadata": {},
+   "source": [
+    "### Task 1: Function implementing causal filters"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0603d817-6c80-444f-a42d-8318f24fdbd7",
+   "metadata": {},
+   "source": [
+    "#### 1.1  Write a function that implements the first order causal low pass filter"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "88aba420-be35-4070-8b6d-386fd924d0e3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def lowpass_first_order(fc, df, nf, napp):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the lowpass first order causal filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    :param napp: number of applications of filter\n",
+    "    \"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b3e445bf-630a-44f0-ae04-7d4829531b58",
+   "metadata": {},
+   "source": [
+    "#### 1.2  Write a function for the causal 2nd order lowpass with filter angle"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "89fda9ef-eca1-484f-9ba2-8ac9046c85f0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def lowpass_second_order(fc, df, nf, phi):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the lowpass second order causal filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    :param phi: filter angle in degrees\n",
+    "    \"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9ee31665-36cb-4e5c-aa31-f8eb26cc2b30",
+   "metadata": {},
+   "source": [
+    "#### 1.3 Write a function fo the Butterworth lowpass of any order"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d21fcbb1-52a3-4fca-8fc8-921034122836",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def butterworth_lowpass(fc, df, nf, order):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the lowpass Butterworth filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    :param order: filter order\n",
+    "    \"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3786d68c-4348-4a24-b2f5-83db76c9b414",
+   "metadata": {},
+   "source": [
+    "### Task 2: Application of filters to seismological data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "259e787f-2f67-43c4-afda-0741087958ad",
+   "metadata": {},
+   "source": [
+    "#### 2.1  Read the trace of CH.PANIX as in the previous assignment and plot it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c9938cab-7ae5-48ed-93d5-83215a350c18",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b015a477-6372-4fc8-a91d-95cb8cbaa73f",
+   "metadata": {},
+   "source": [
+    "#### 2.2 Compute the fast Fourier coefficients of the data and print number of frequencies and max frequency"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "56eb0313-e1a2-404e-9bc1-0ded3791f345",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f4ec9eea-a114-4e33-ac3f-76b7acfbc0ab",
+   "metadata": {},
+   "source": [
+    "#### 2.3 Apply the first order lowpass, do Fourier synthesis and plot the result. Use fc=0.1 Hz and vary corner frequency later. Set axis limits to zoom into the P-wave train."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "afc7c27a-173d-4339-b4ac-abfbd6b5e974",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "c3ccc36a-e666-4cca-b894-5a4e7c9aa3b7",
+   "metadata": {},
+   "source": [
+    "#### 2.4 Apply the first order lowpass 2 times, do Fourier synthesis and plot the result."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "50e6cc4c-4c9c-47a5-ba4f-04808b64db31",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\" Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f7c19547-cd98-49a2-a6d1-7e14e5a058e5",
+   "metadata": {},
+   "source": [
+    "#### 2.5 Apply the first order low pass 3 times, do Fourier synthesis and plot the result"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a54f6ee3-480e-4178-b4d6-04a01d47a3c7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "6b8d3381-5bbe-4558-ae08-6f1c2292716d",
+   "metadata": {},
+   "source": [
+    "#### 2.6 Apply the second order lowpass for phi = 5 degrees, do Fourier synthesis and plot the result"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "45fd47e8-0400-4c3d-bb2a-1ba15bb27c84",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "251422f2-6af2-4aa6-83cb-aacd31ddde7b",
+   "metadata": {},
+   "source": [
+    "#### 2.7 Apply the second order lowpass for phi = 45 degrees, do Fourier synthesis and plot the result"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0b292eda-a26d-443c-8c58-b16011d63403",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "66f56f5c-6b64-4933-9004-ea2a450f6bf0",
+   "metadata": {},
+   "source": [
+    "#### 2.8 Setup a loop over filter orders from 1 to 8 and apply the low pass Butterworth filter to the data. Do Fourier synthesis and plot the resulting traces in one graph in such a way that the traces are vertically shifted relative to each other. Here, it is usefull to normalize the traces before plotting."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0de14535-8b80-4921-95b2-e4cb917f9bac",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Here comes you code\"\"\""
+   ]
+  }
+ ],
+ "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
+}

assignment_06/seismicTrace.py

@@ -0,0 +1,79 @@
+#----------------------------------------------------------------------------
+#  Python class for a seismic trace
+#
+import numpy as np
+
+class SeismicTrace(object):
+    """
+    Class defining a seismic trace by station name, network name, component,
+    number of samples, sampling interval and start time
+    """
+
+    def __init__(self, staname, comp, dt, tanf, vals):
+        """
+        Initialize seismic trace object from data:
+        :param staname: station name
+        :param comp: component
+        :param dt: sampling interval
+        :param tanf: start time
+        :param vals: values of trace (1D numpy array)
+        """
+        self.staname = staname
+        self.comp = comp
+        self.dt = dt
+        self.tanf = tanf
+        self.vals = vals
+        self.nsamp = vals.size
+        self.tlen = self.dt*self.nsamp
+
+
+    @classmethod
+    def readFromAscii(cls, fid, comp='None'):
+        """
+        Class method as second constructor to initialize trace from an Ascii file
+        :param fid: file descriptor of already open Ascii file
+        :param comp: optional component selection
+        :return: an instance of the class
+        """
+        staname = fid.readline().strip()
+        if not staname:
+            print("No station name found. Is file really an Ascii gather file")
+            return None
+        compf = fid.readline().strip()
+        if comp != 'None' and comp != compf:
+            print("Requested component %s not found.", comp)
+            return None
+        h = fid.readline().strip().split()
+        vals = np.array([float(v) for v in fid.readline().strip().split()])
+        dt = float(h[1]); tanf = float(h[2])
+        return cls(staname, compf, dt, tanf, vals)
+
+
+    def absoluteMaximum(self):
+        """
+        return absolute maximum of trace
+        """
+        return abs(max(self.vals))
+
+
+    def writeToOpenAsciiFile(self, fid):
+        """
+        write a seismic trace to an open Ascii file
+        """
+        fid.write(self.staname + '\n')
+        fid.write(self.comp + '\n')
+        fid.write("{:12d} {:15.9f} {:20.9f}\n".format(self.nsamp, self.dt, self.tanf))
+        for n in range(self.nsamp):
+            fid.write("{:17.8e}".format(self.vals[n]))
+        fid.write("\n")
+
+
+    def printInfo(self):
+        """
+        Print out some information about the trace
+        """
+        print("Identification: ", self.staname, self.comp)
+        print("Sampling interval: ", self.dt)
+        print("Number of samples: ", self.nsamp)
+        print("Start time after midnight: ", self.tanf)
+        print("Length of trace: (N*DT) ", self.tlen)

assignment_06/setupFigure.py

@@ -0,0 +1,15 @@
+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

assignment_07/dftFast.py

@@ -0,0 +1,27 @@
+# -------------------------------------------------------------
+#  Functions for "fast" Fourier transform and Fourier synthesis
+#  using numpy tools
+# -------------------------------------------------------------
+import numpy as np
+
+def dft_fast_coeff(x):
+    """
+    Evaluate Fourier coefficients using Numpy's fast Fourier transform.
+    This routine only returns the coefficients for the positive frequencies.
+    If N is even, it goes up to n=N/2.
+    If N is odd, it goes up to n=(N-1)/2.
+    :param x: array of function samples
+    """
+    return np.fft.rfft(x, norm='forward')
+
+
+def dft_fast_synthesis(fc, outnum='even'):
+    """
+    Use numpy's fast Fourier synthesis taking only Fourier coefficients for positive frequencies
+    :param fc: aray of coefficients for positive frequencies only.
+    :param outnum: specifies if output time series has an even or odd number of samples (default: 'even')
+    """
+    ns = 2*fc.size-2
+    if outnum == 'odd': ns = 2*fc.size-1
+    return np.fft.irfft(fc, ns, norm='forward')
+

assignment_07/p7_seismometer_correlation.ipynb

@@ -0,0 +1,308 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "17de1a34-2933-47ed-b1c5-ae9dd0466edd",
+   "metadata": {},
+   "source": [
+    "## Assignment 7"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "681b47c6-f9af-4877-b745-febd443e1877",
+   "metadata": {},
+   "source": [
+    "### Seismometer transfer function and correlation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "85cb9941-eb24-4eb2-9c46-8610844e4063",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from setupFigure import SetupFigure\n",
+    "from dftFast import dft_fast_coeff, dft_fast_synthesis\n",
+    "from seismicTrace import SeismicTrace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a627c30b-eb67-4319-930e-7f1acb8e3ee1",
+   "metadata": {},
+   "source": [
+    "#### Task 1: Write a function that calculates the transfer function of a seismometer for different damping $h=\\sin\\phi$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8cbc4a71-fc18-4c3e-a26a-0c02ae89314d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def transfer_seismometer_velocity(fc, df, nf, phi, sigma):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the ground velocity seismometer transfer function (positive frequencies only)\n",
+    "    :param fc: eigenfrequency of seismometer in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    :param phi: damping angle in degrees (h=sin(phi))\n",
+    "    :param sigma: generator constant\n",
+    "    :return: array of frequencies\n",
+    "    :return: array of complex values of transfer function\n",
+    "    \"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f7d857c-1bfd-4703-96a7-61c7af95bbdc",
+   "metadata": {},
+   "source": [
+    "#### You may use this function for the Butterworth lowpass filter (or your own from Assignment 6)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dd1b00e9-a953-4d64-b91f-7d7ad6949a3f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def butterworth_lowpass(fc, df, nf, order):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the lowpass Butterworth filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping in Hz\n",
+    "    :param nf: number of frequencies\n",
+    "    :param order: filter order\n",
+    "    \"\"\"\n",
+    "    f = df*np.arange(0, nf)\n",
+    "    h = np.ones(nf, dtype=np.complex64)\n",
+    "    for k in range(order):\n",
+    "        arg = 1j*(2*k+1)*np.pi/(2*order)\n",
+    "        h = h*(-1j)/(f/fc-np.exp(arg))\n",
+    "    return h"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "38d6be7b-6579-4726-abee-4666e6910f39",
+   "metadata": {},
+   "source": [
+    "#### This is a useful function to plot traces next to each other with a vertical offset"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "48a7e17e-eaf2-4162-89f2-0d81eb9b7631",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plotTraceGather(tracelist, ax):\n",
+    "    \"\"\"\n",
+    "    Plots a gather of traces with vertical offset\n",
+    "    :param tracelist: list of SeismicTrace objects\n",
+    "    :param ax: Matplotlib axis object\n",
+    "    \"\"\"\n",
+    "    for j, tr in enumerate(tracelist):\n",
+    "        t = tr.tanf+tr.dt*np.arange(0, tr.nsamp)\n",
+    "        norm = np.max(np.abs(tr.vals))\n",
+    "        ax.plot(t, 0.5*tr.vals/norm+j, color='b', ls='-')   # Normalization to 0.5*absmax and shift by multiples of 1\n",
+    "        ax.text(tr.tanf, j+0.1, tr.staname, fontsize=14)    # Add station name on top of trace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "68bf389c-4464-43dc-853d-ff28e5791f99",
+   "metadata": {},
+   "source": [
+    "#### Task 2.1: Write a function that computes the cross correlation function for a given range of time delays"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "47a92e2f-da8c-45ce-b58c-7b8d153e38f5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def crossCorrelation(x, y, dt, neglim, poslim):\n",
+    "    \"\"\"\n",
+    "    Correlate x with y according to c_n = Delta t*sum_(k=0)^(N-1) x_(n+k) y_k.\n",
+    "    Assumes that both traces have same DT!\n",
+    "    Consider x is a Gaussian with max at t2, y a Gaussian with max at t1 < t2.\n",
+    "    To bring x and y into a match, x neeeds to be shifted to the left implying a positive n=n_s.\n",
+    "    Thus, ns*dt is the delay of x relative to y or the advance of y relative to x.\n",
+    "    If n+k < 0 or n+k >= x.size, x_(n+k) is considered as zero.\n",
+    "    :param x: values of first trace\n",
+    "    :param y: values of second trace\n",
+    "    :param dt: sampling interval\n",
+    "    :param neglim: most negative value for n\n",
+    "    :param poslim: most positive value for n (poslim not included)\n",
+    "    :return: array of delay times associated with the array elements of the cross correlation\n",
+    "    :return: array of values of the cross correlation function\n",
+    "    \"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "79f44179-4943-4e87-9505-4c326c4237c2",
+   "metadata": {},
+   "source": [
+    "#### Task 2.2: Evaluate the amplitude spectrum of the seismometer transfer function for different angles\n",
+    "Use fc = 10 Hz, nf = 256, df = 0.2, sigma = 1. Choose angles 15, 25, 45, 60, 85 degrees. Plot the\n",
+    "amplitude spectra into one graph."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "56a730c7-ba3c-4404-a8e4-6d8c83d5fab3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b2576de7-c170-4829-a5a9-a107e6ad09bb",
+   "metadata": {
+    "jp-MarkdownHeadingCollapsed": true
+   },
+   "source": [
+    "#### Task 2.3: Evaluate the phase spectrum of the seismometer transfer function for different angles\n",
+    "Use fc = 10 Hz, nf = 256, df = 0.2, sigma = 1. Choose angles 15, 25, 45, 60, 85 degrees. Plot the\n",
+    "amplitude spectra into one graph."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a796b315-f5e4-48cc-ae2e-d98323cc4996",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d220be85-9530-418d-9069-6ad21d87872a",
+   "metadata": {},
+   "source": [
+    "#### Task 2.4: Setup two Gaussians $e^{-(t-\\mu)^2/\\sigma^2}$ of lengths 100 s and 80 s with dt=1 s, $\\mu$=55 and 45, and $\\sigma$ = 5. Plot the two into one graph. Then, compute the cross correlation between -30 and +30 s using your function from Task 2.1 and plot it into a separate graph. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f7537491-54c9-4f9c-ac22-e02d9cefae05",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1529c72b-6ae3-416a-b4b2-ed1fd0bffba3",
+   "metadata": {},
+   "source": [
+    "#### Task 2.5: For comparison, use numpy's routine correlate to calculate the correlation. Try the different modes offered. Change the lengths of x and y."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "da38f6a3-4d2e-4a69-a605-7de5a55f5cfd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5226e897-3874-45ba-86ab-883dd92fc13e",
+   "metadata": {},
+   "source": [
+    "#### Task 3.1: We want to prepare tome delay measurements between two seismograms by cross-correlation. First, read the data of stations CH.PANIX and Z3.A261A.HHZ. Put the SeismicTrace objects into a list."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0655529e-3212-4069-a60f-bef57644328c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9f38af03-65f9-49d8-a28d-5fc0e2c4d965",
+   "metadata": {},
+   "source": [
+    "#### Task 3.2: Plot the traces using the function plotTraceGather."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "56551bfa-e2da-4773-9247-a33cae96a3c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Your code\"\"\""
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cec5430e-c602-4a99-a7c3-f3371e990d82",
+   "metadata": {},
+   "source": [
+    "#### Task 3.3: Filter the data using a Butterworth low pass with fc=0.5 Hz and order 6. Use either the provided function for the low-pass transfer function or take your own. First, Fourier transform the data, then multiply with the filter transfer function, and finally do the back transform. Create new SeismicTrace objects for the filtered traces, put them into a list and plot them."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6ce5cc83-f5ff-4fb2-a9ac-d98595429c02",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "\"\"\"Your code\"\"\""
+   ]
+  }
+ ],
+ "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
+}

assignment_07/seismicTrace.py

@@ -0,0 +1,79 @@
+#----------------------------------------------------------------------------
+#  Python class for a seismic trace
+#
+import numpy as np
+
+class SeismicTrace(object):
+    """
+    Class defining a seismic trace by station name, network name, component,
+    number of samples, sampling interval and start time
+    """
+
+    def __init__(self, staname, comp, dt, tanf, vals):
+        """
+        Initialize seismic trace object from data:
+        :param staname: station name
+        :param comp: component
+        :param dt: sampling interval
+        :param tanf: start time
+        :param vals: values of trace (1D numpy array)
+        """
+        self.staname = staname
+        self.comp = comp
+        self.dt = dt
+        self.tanf = tanf
+        self.vals = vals
+        self.nsamp = vals.size
+        self.tlen = self.dt*self.nsamp
+
+
+    @classmethod
+    def readFromAscii(cls, fid, comp='None'):
+        """
+        Class method as second constructor to initialize trace from an Ascii file
+        :param fid: file descriptor of already open Ascii file
+        :param comp: optional component selection
+        :return: an instance of the class
+        """
+        staname = fid.readline().strip()
+        if not staname:
+            print("No station name found. Is file really an Ascii gather file")
+            return None
+        compf = fid.readline().strip()
+        if comp != 'None' and comp != compf:
+            print("Requested component %s not found.", comp)
+            return None
+        h = fid.readline().strip().split()
+        vals = np.array([float(v) for v in fid.readline().strip().split()])
+        dt = float(h[1]); tanf = float(h[2])
+        return cls(staname, compf, dt, tanf, vals)
+
+
+    def absoluteMaximum(self):
+        """
+        return absolute maximum of trace
+        """
+        return abs(max(self.vals))
+
+
+    def writeToOpenAsciiFile(self, fid):
+        """
+        write a seismic trace to an open Ascii file
+        """
+        fid.write(self.staname + '\n')
+        fid.write(self.comp + '\n')
+        fid.write("{:12d} {:15.9f} {:20.9f}\n".format(self.nsamp, self.dt, self.tanf))
+        for n in range(self.nsamp):
+            fid.write("{:17.8e}".format(self.vals[n]))
+        fid.write("\n")
+
+
+    def printInfo(self):
+        """
+        Print out some information about the trace
+        """
+        print("Identification: ", self.staname, self.comp)
+        print("Sampling interval: ", self.dt)
+        print("Number of samples: ", self.nsamp)
+        print("Start time after midnight: ", self.tanf)
+        print("Length of trace: (N*DT) ", self.tlen)

assignment_07/seismometer_correlation.ipynb

@@ -0,0 +1,448 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "17de1a34-2933-47ed-b1c5-ae9dd0466edd",
+   "metadata": {},
+   "source": [
+    "## Assignment 7"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "681b47c6-f9af-4877-b745-febd443e1877",
+   "metadata": {},
+   "source": [
+    "### Seismometer transfer function and correlation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "85cb9941-eb24-4eb2-9c46-8610844e4063",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from setupFigure import SetupFigure\n",
+    "from dftFast import dft_fast_coeff, dft_fast_synthesis\n",
+    "from seismicTrace import SeismicTrace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a627c30b-eb67-4319-930e-7f1acb8e3ee1",
+   "metadata": {},
+   "source": [
+    "#### Task 1: Write a function that calculates the transfer function of a seismometer for different damping $h=\\sin\\phi$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8cbc4a71-fc18-4c3e-a26a-0c02ae89314d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def transfer_seismometer_velocity(fc, df, nf, phi, sigma):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the ground velocity seismometer transfer function (positive frequencies only)\n",
+    "    :param fc: eigenfrequency of seismometer in Hz\n",
+    "    :param df: frequency stepping\n",
+    "    :param nf: number of frequencies\n",
+    "    :param phi: damping angle in degrees (h=sin(phi))\n",
+    "    :param sigma: generator constant\n",
+    "    :return array with frequencies, array of complex values of transfer function\n",
+    "    \"\"\"\n",
+    "    f = df*np.arange(0, nf)\n",
+    "    alfa = np.deg2rad(phi)\n",
+    "    return f, -sigma*(f/fc)/(f/fc-np.exp(1j*alfa)) * (f/fc)/(f/fc-np.exp(1j*(np.pi-alfa)))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f7d857c-1bfd-4703-96a7-61c7af95bbdc",
+   "metadata": {},
+   "source": [
+    "#### You may use this function for the Butterworth lowpass filter (or your own from Assignment 6)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dd1b00e9-a953-4d64-b91f-7d7ad6949a3f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def butterworth_lowpass(fc, df, nf, order):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the lowpass Butterworth filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping in Hz\n",
+    "    :param nf: number of frequencies\n",
+    "    :param order: filter order\n",
+    "    \"\"\"\n",
+    "    f = df*np.arange(0, nf)\n",
+    "    h = np.ones(nf, dtype=np.complex64)\n",
+    "    for k in range(order):\n",
+    "        arg = 1j*(2*k+1)*np.pi/(2*order)\n",
+    "        h = h*(-1j)/(f/fc-np.exp(arg))\n",
+    "    return h"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "38d6be7b-6579-4726-abee-4666e6910f39",
+   "metadata": {},
+   "source": [
+    "#### This is a useful function to plot traces next to each other with a vertical offset"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "48a7e17e-eaf2-4162-89f2-0d81eb9b7631",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plotTraceGather(tracelist, ax):\n",
+    "    \"\"\"\n",
+    "    Plots a gather of traces with vertical offset\n",
+    "    :param tracelist: list of SeismicTrace objects\n",
+    "    :param ax: Matplotlib axis object\n",
+    "    \"\"\"\n",
+    "    for j, tr in enumerate(tracelist):\n",
+    "        t = tr.tanf+tr.dt*np.arange(0, tr.nsamp)\n",
+    "        norm = np.max(np.abs(tr.vals))\n",
+    "        ax.plot(t, 0.5*tr.vals/norm+j, color='b', ls='-')   # Normalization to 0.5*absmax and shift by multiples of 1\n",
+    "        ax.text(tr.tanf, j+0.1, tr.staname, fontsize=14)    # Add station name on top of trace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "68bf389c-4464-43dc-853d-ff28e5791f99",
+   "metadata": {},
+   "source": [
+    "#### Task 2.1: Write a function that computes the cross correlation function for a given range of time delays"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "47a92e2f-da8c-45ce-b58c-7b8d153e38f5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def crossCorrelation(x, y, dt, neglim, poslim):\n",
+    "    \"\"\"\n",
+    "    Correlate x with y according to c_n = Delta t*sum_(k=0)^(N-1) x_(n+k) y_k.\n",
+    "    Assumes that both traces have same DT!\n",
+    "    Consider x is a Gaussian with max at t2, y a Gaussian with max at t1 < t2.\n",
+    "    To bring x and y into a match, x neeeds to be shifted to the left implying a positive n=n_s.\n",
+    "    Thus, ns*dt is the delay of x relative to y or the advance of y relative to x.\n",
+    "    If n+k < 0 or n+k >= x.size, x_(n+k) is considered as zero.\n",
+    "    :param x: values of first trace\n",
+    "    :param y: values of second trace\n",
+    "    :param dt: sampling interval\n",
+    "    :param neglim: most negative value for n\n",
+    "    :param poslim: most positive value for n (poslim not included)\n",
+    "    :return: array of delay times associated with the array elements of the cross correlation\n",
+    "    :return: array of values of the cross correlation function\n",
+    "    \"\"\"\n",
+    "    nsamp = poslim-neglim\n",
+    "    cc = np.zeros(nsamp)\n",
+    "    for n in range(neglim, poslim):\n",
+    "        nc = n-neglim                  # adjust index associating index 0 to neglim\n",
+    "        for k in range(y.size):\n",
+    "            if n+k >= 0 and n+k < x.size:\n",
+    "                cc[nc] = cc[nc]+x[n+k]*y[k]\n",
+    "    return dt*np.arange(neglim, poslim), cc*dt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "79f44179-4943-4e87-9505-4c326c4237c2",
+   "metadata": {},
+   "source": [
+    "#### Task 2.2: Evaluate the amplitude spectrum of the seismometer transfer function for different angles\n",
+    "Use fc = 10 Hz, nf = 256, df = 0.2, sigma = 1. Choose angles 15, 25, 45, 60, 85 degrees. Plot the\n",
+    "amplitude spectra into one graph."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "56a730c7-ba3c-4404-a8e4-6d8c83d5fab3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fc = 10\n",
+    "nf = 256\n",
+    "df = 0.2\n",
+    "sigma = 1\n",
+    "phi = [15, 25, 45, 60, 85]\n",
+    "col = ['black', 'red', 'blue', 'orange', 'magenta']"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "dcfdc609-996d-425a-9750-f10cc453ef89",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig1, ax1 = SetupFigure(12, 5, \"Frequency [Hz]\", \"$|H(f)|$\", \"Amplitude of seismometer transfer function\", 14)\n",
+    "for j, p in enumerate(phi):\n",
+    "    f, hseis = transfer_seismometer_velocity(fc, df, nf, p, sigma)\n",
+    "    ax1.plot(f, np.absolute(hseis), color=col[j], ls='-', label=\"Phi = {:4.0f}\".format(p))\n",
+    "ax1.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b2576de7-c170-4829-a5a9-a107e6ad09bb",
+   "metadata": {},
+   "source": [
+    "#### Task 2.3: Evaluate the phase spectrum of the seismometer transfer function for different angles\n",
+    "Use fc = 10 Hz, nf = 256, df = 0.2, sigma = 1. Choose angles 15, 25, 45, 60, 85 degrees. Plot the\n",
+    "amplitude spectra into one graph."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a796b315-f5e4-48cc-ae2e-d98323cc4996",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig2, ax2 = SetupFigure(12, 5, \"Frequency [Hz]\", \"Phase / PI\", \"Phase of seismometer transfer function\", 14)\n",
+    "for j, p in enumerate(phi):\n",
+    "    f, hseis = transfer_seismometer_velocity(fc, df, nf, p, sigma)\n",
+    "    ax2.plot(f, np.angle(hseis)/np.pi, color=col[j], ls='-', label=\"Phi = {:4.0f}\".format(p))\n",
+    "ax2.legend()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d220be85-9530-418d-9069-6ad21d87872a",
+   "metadata": {},
+   "source": [
+    "#### Task 2.4: Setup two Gaussians $e^{-(t-\\mu)^2/\\sigma^2}$ of lengths 100 s and 80 s with dt=1 s, $\\mu$=45 and 55 and $\\sigma$ = 5. Plot the two into one graph. Then, compute the cross correlation between -30 and +30 s using your function from Task 2 and plot it into a separate graph. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f7537491-54c9-4f9c-ac22-e02d9cefae05",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tx = np.arange(0,100)\n",
+    "ty = np.arange(0,80)\n",
+    "x = np.exp(-((tx-45)/5)**2)\n",
+    "y = np.exp(-((ty-55)/5)**2)\n",
+    "fig7, ax7 = SetupFigure(12, 5, \"Time [s]\", \"Amplitude\", \"Two Gaussians\", 14)\n",
+    "ax7.plot(tx, x, color='blue', ls='-', label='Time series x(t)')\n",
+    "ax7.plot(ty, y, color='red', ls='-', label='Time series y(t)')\n",
+    "ax7.legend()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "624805e5-2feb-4a42-bd91-5b48e798b01a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "tc, cc = crossCorrelation(x, y, 1.0, -30, 30)\n",
+    "fig8, ax8 = SetupFigure(12, 5, \"Time delay [s]\", \"Correlation\", \"Cross correlation\", 14)\n",
+    "ax8.plot(tc, cc, color='blue', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1529c72b-6ae3-416a-b4b2-ed1fd0bffba3",
+   "metadata": {},
+   "source": [
+    "#### Task 2.5: For comparison, use numpy's routine correlate to calculate the correlation. Try the different modes offered.\n",
+    "Change the lengths of x and y."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "da38f6a3-4d2e-4a69-a605-7de5a55f5cfd",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig9, ax9 = SetupFigure(12, 5, \"Time delay [s]\", \"Correlation\", \"Numpy cross correlation\", 14)\n",
+    "ccnp = np.correlate(x, y, mode='full')\n",
+    "tcnp = np.arange(-y.size+1, x.size)    # This seems to be the correct relation between lag and index\n",
+    "print(x.size, y.size)\n",
+    "ax9.set_xlim(-30,30)\n",
+    "ax9.plot(tcnp, ccnp, color='red', ls='-')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5226e897-3874-45ba-86ab-883dd92fc13e",
+   "metadata": {},
+   "source": [
+    "#### Task 3.1: We want to prepare time delay measurements between two seismograms by cross-correlation. First, read the data of stations CH.PANIX and Z3.A261A.HHZ. Put the SeismicTrace objects into a list."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0655529e-3212-4069-a60f-bef57644328c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "stalist = [\"CH.PANIX.HHZ\", \"Z3.A261A.HHZ\"]\n",
+    "rawlist = []\n",
+    "for sta in stalist:\n",
+    "    with open(sta, 'r') as fd:\n",
+    "        tr = SeismicTrace.readFromAscii(fd)\n",
+    "    tr.printInfo()\n",
+    "    rawlist.append(tr)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "9f38af03-65f9-49d8-a28d-5fc0e2c4d965",
+   "metadata": {},
+   "source": [
+    "#### Task 3.2: Plot the traces using the function plotTraceGather."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "56551bfa-e2da-4773-9247-a33cae96a3c1",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig3, ax3 = SetupFigure(12, 5, \"Time [s]\", \"Displacement\", \"Raw traces\", 14)\n",
+    "plotTraceGather(rawlist, ax3)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cec5430e-c602-4a99-a7c3-f3371e990d82",
+   "metadata": {},
+   "source": [
+    "#### Task 3.3: Filter the data using a Butterworth low pass with fc=0.5 Hz and order 6. Use either the provided function for the low-pass transfer function or take your own. First, Fourier transform the data, then multiply with the filter transfer function, and finally do the back transform. Create new SeismicTrace objects for the filtered traces, put them into a list and plot them."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "6ce5cc83-f5ff-4fb2-a9ac-d98595429c02",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fc = 0.5; order = 6\n",
+    "filtlist = []\n",
+    "for tr in rawlist:\n",
+    "    c = dft_fast_coeff(tr.vals)\n",
+    "    df = 1./tr.tlen\n",
+    "    hb = butterworth_lowpass(fc, df, c.size, order)\n",
+    "    fts = c*hb\n",
+    "    s = dft_fast_synthesis(fts, outnum='even')\n",
+    "    trfil = SeismicTrace(tr.staname, tr.comp, tr.dt, tr.tanf, s)\n",
+    "    filtlist.append(trfil)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2fd73161-ade3-47e9-9880-774459625017",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig4, ax4 = SetupFigure(12, 5, \"Time [s]\", \"Displacement\", \"Low pass filtered traces, fc = {:5.2f} Hz, Order = {:d}\".format(fc, order), 14)\n",
+    "plotTraceGather(filtlist, ax4)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ca21ecc8-0d60-46c4-8be8-c1e8fbe2c940",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "winlist = [[64285, 64300], [64280, 64295.1]]\n",
+    "cutlist = []\n",
+    "for j, tr in enumerate(filtlist):\n",
+    "    win = winlist[j]\n",
+    "    jfirst = int((win[0]-tr.tanf)/tr.dt)\n",
+    "    jlast = int((win[1]-tr.tanf)/tr.dt)\n",
+    "    nsamp = jlast-jfirst\n",
+    "    tanew = jfirst*tr.dt+tr.tanf\n",
+    "    print(tr.tanf, tanew, jfirst, jlast)\n",
+    "    cuttr = SeismicTrace(tr.staname, tr.comp, tr.dt, tanew, tr.vals[jfirst:jlast])\n",
+    "    cutlist.append(cuttr)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f006344c-4728-4595-988e-d018bc88f994",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig5, ax5 = SetupFigure(12, 5, \"Time [s]\", \"Displacement\", \"Cut low pass filtered traces\".format(fc, order), 14)\n",
+    "plotTraceGather(cutlist, ax5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "4993c95b-831d-490f-927c-cb38a2de1bfb",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "dt = cutlist[0].dt\n",
+    "neglim = -500; poslim = 500\n",
+    "t, cc = crossCorrelation(cutlist[0].vals, cutlist[1].vals, dt, neglim, poslim)\n",
+    "delay = (np.argmax(cc)+neglim)*dt+cutlist[0].tanf-cutlist[1].tanf\n",
+    "print(\"Tanf_0 = \", cutlist[0].tanf, \"Tanf_1 = \", cutlist[1].tanf, \"Index of max cc: \", np.argmax(cc)+neglim)\n",
+    "print(\"Delay of \", cutlist[0].staname,\" relative to \", cutlist[1].staname,\":\", delay)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "65cf3ace-f9ce-4da0-ac50-d348e7a129b3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig6, ax6 = SetupFigure(12, 5, \"Time delay [s]\", \"Correlation\", \"Cross correlation\", 14)\n",
+    "ax6.plot(t, cc, color='blue', ls='-')"
+   ]
+  }
+ ],
+ "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
+}

assignment_07/setupFigure.py

@@ -0,0 +1,15 @@
+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

assignment_08/dftFast.py

@@ -0,0 +1,27 @@
+# -------------------------------------------------------------
+#  Functions for "fast" Fourier transform and Fourier synthesis
+#  using numpy tools
+# -------------------------------------------------------------
+import numpy as np
+
+def dft_fast_coeff(x):
+    """
+    Evaluate Fourier coefficients using Numpy's fast Fourier transform.
+    This routine only returns the coefficients for the positive frequencies.
+    If N is even, it goes up to n=N/2.
+    If N is odd, it goes up to n=(N-1)/2.
+    :param x: array of function samples
+    """
+    return np.fft.rfft(x, norm='forward')
+
+
+def dft_fast_synthesis(fc, outnum='even'):
+    """
+    Use numpy's fast Fourier synthesis taking only Fourier coefficients for positive frequencies
+    :param fc: aray of coefficients for positive frequencies only.
+    :param outnum: specifies if output time series has an even or odd number of samples (default: 'even')
+    """
+    ns = 2*fc.size-2
+    if outnum == 'odd': ns = 2*fc.size-1
+    return np.fft.irfft(fc, ns, norm='forward')
+

assignment_08/p8_correlation_picking.ipynb

@@ -0,0 +1,287 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "17de1a34-2933-47ed-b1c5-ae9dd0466edd",
+   "metadata": {},
+   "source": [
+    "## Assignment 8"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "681b47c6-f9af-4877-b745-febd443e1877",
+   "metadata": {},
+   "source": [
+    "### Correlation picking"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 14,
+   "id": "85cb9941-eb24-4eb2-9c46-8610844e4063",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from setupFigure import SetupFigure\n",
+    "from dftFast import dft_fast_coeff, dft_fast_synthesis\n",
+    "from seismicTrace import SeismicTrace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7f7d857c-1bfd-4703-96a7-61c7af95bbdc",
+   "metadata": {},
+   "source": [
+    "#### You may use this function for the Butterworth lowpass filter (or your own from Assignment 6)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "id": "dd1b00e9-a953-4d64-b91f-7d7ad6949a3f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def butterworth_lowpass(fc, df, nf, order):\n",
+    "    \"\"\"\n",
+    "    Calculates samples of the lowpass Butterworth filter transfer function (positive frequencies only)\n",
+    "    :param fc: corner frequency in Hz\n",
+    "    :param df: frequency stepping in Hz\n",
+    "    :param nf: number of frequencies\n",
+    "    :param order: filter order\n",
+    "    \"\"\"\n",
+    "    f = df*np.arange(0, nf)\n",
+    "    h = np.ones(nf, dtype=np.complex64)\n",
+    "    for k in range(order):\n",
+    "        arg = 1j*(2*k+1)*np.pi/(2*order)\n",
+    "        h = h*(-1j)/(f/fc-np.exp(arg))\n",
+    "    return h"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "38d6be7b-6579-4726-abee-4666e6910f39",
+   "metadata": {},
+   "source": [
+    "#### This is a useful function to plot traces next to each other with a vertical offset"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "id": "48a7e17e-eaf2-4162-89f2-0d81eb9b7631",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def plotTraceGather(tracelist, ax):\n",
+    "    \"\"\"\n",
+    "    Plots a gather of traces with vertical offset\n",
+    "    :param tracelist: list of SeismicTrace objects\n",
+    "    :param ax: Matplotlib axis object\n",
+    "    \"\"\"\n",
+    "    for j, tr in enumerate(tracelist):\n",
+    "        t = tr.tanf+tr.dt*np.arange(0, tr.nsamp)\n",
+    "        norm = np.max(np.abs(tr.vals))\n",
+    "        ax.plot(t, 0.5*tr.vals/norm+j, color='b', ls='-')   # Normalization to 0.5*absmax and shift by multiples of 1\n",
+    "        ax.text(tr.tanf, j+0.1, tr.staname, fontsize=14)    # Add station name on top of trace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "68bf389c-4464-43dc-853d-ff28e5791f99",
+   "metadata": {},
+   "source": [
+    "#### You may use this function that computes the cross correlation function for a given range of time delays"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "id": "47a92e2f-da8c-45ce-b58c-7b8d153e38f5",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def crossCorrelation(x, y, dt, neglim, poslim):\n",
+    "    \"\"\"\n",
+    "    Correlate x with y according to c_n = Delta t*sum_(k=0)^(N-1) x_(n+k) y_k.\n",
+    "    Assumes that both traces have same DT!\n",
+    "    Consider x is a Gaussian with max at t2, y a Gaussian with max at t1 < t2.\n",
+    "    To bring x and y into a match, x neeeds to be shifted to the left implying a positive n=n_s.\n",
+    "    Thus, ns*dt is the delay of x relative to y or the advance of y relative to x.\n",
+    "    If n+k < 0 or n+k >= x.size, x_(n+k) is considered as zero.\n",
+    "    :param x: values of first trace\n",
+    "    :param y: values of second trace\n",
+    "    :param dt: sampling interval\n",
+    "    :param neglim: most negative value for n\n",
+    "    :param poslim: most positive value for n (poslim not included)\n",
+    "    :return: array of delay times associated with the array elements of the cross correlation\n",
+    "    :return: array of values of the cross correlation function\n",
+    "    \"\"\"\n",
+    "    nsamp = poslim-neglim\n",
+    "    cc = np.zeros(nsamp)\n",
+    "    for n in range(neglim, poslim):\n",
+    "        nc = n-neglim                  # adjust index associating index 0 to neglim\n",
+    "        for k in range(y.size):\n",
+    "            if n+k >= 0 and n+k < x.size:\n",
+    "                cc[nc] = cc[nc]+x[n+k]*y[k]\n",
+    "    return dt*np.arange(neglim, poslim), cc*dt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3404b2f6-d296-4b3e-8234-a01c2bd061ea",
+   "metadata": {},
+   "source": [
+    "### We repeat the first parts of task 3 from the previous assignment to make this assignment self-contained."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5226e897-3874-45ba-86ab-883dd92fc13e",
+   "metadata": {},
+   "source": [
+    "#### Task 1.1: We want to do time delay measurements between two seismograms by cross-correlation. First, read the data of stations CH.PANIX and Z3.A261A.HHZ. Put the SeismicTrace objects into a list."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b8182997-34f6-4091-8f71-caa942e9d1f2",
+   "metadata": {},
+   "source": [
+    "#### Task 1.2: Plot the traces using the function plotTraceGather. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "951456c6-a3d2-4737-a803-c88d813655e4",
+   "metadata": {},
+   "source": [
+    "This function plots the seismogram values versus time in seconds after midnight."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "db6ea78c-03ec-4e56-88c5-c7dc58efed9c",
+   "metadata": {},
+   "source": [
+    "#### Task 1.3: Filter the data using a Butterworth low pass with fc=0.5 Hz and order 6. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4f3c59ab-69b7-4ab0-bedc-e525134a687f",
+   "metadata": {},
+   "source": [
+    "Use either the provided function for the low-pass transfer function or take your own. First, Fourier transform the data, then multiply with the filter transfer function, and finally do the back transform. Create new SeismicTrace objects for the filtered traces, put them into a list and plot them."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3f1990b2-9684-4a1e-a678-448ebd0b4afe",
+   "metadata": {},
+   "source": [
+    "####  Task 1.4: Cut out the P-wave train from the recordings. Try different time windows for each trace. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "e619d197-db9b-48ef-840d-aa123b5ff884",
+   "metadata": {},
+   "source": [
+    "As a start, use the times [64285, 64300] for CH.PANIX and the times [64280, 64295] for Z3.A261A. For cutting out samples, proceed as follows: (1) Calculate the sample indices for the beginning and end of the window, (2) calculate the new starting time, (3) create a new SeismicTrace object with same station name, component and DT but adjusted start time (tanf) and the data values within the\n",
+    "precalculated index range. (4) Again, put the windowed seismograms into a list."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "058f0df7-d40a-4a57-bec0-052fa2eb1c86",
+   "metadata": {},
+   "source": [
+    "#### Task 1.5: Plot the windowed traces using plotTraceGather."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "59ff1b7f-40f9-42d4-9354-10b14a3516bc",
+   "metadata": {},
+   "source": [
+    "#### Task 1.6: Cross-correlate the windowed seismograms and compute the time delay"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "175aabb2-a1b0-4c85-a84c-7f1e84ed2546",
+   "metadata": {},
+   "source": [
+    "Compute the correlation for delays from -500 to +500 samples. To get the time delay, search the index jmax of the maximum of the cross-correlation. Since sample zero of the cross-correlation belongs to a time of -500*dt, the time delay is given by $\\tau=(jmax-500)dt$. This is correct if the start times of the traces are identical. If not, we must correct for the start time difference. For example, if you correlate PANIX with A261A in this order, you obtain the delay of PANIX relative to A261A which is given by $\\tau$ + ta(PANIX)-ta(A261A).\n",
+    "\n",
+    "Print the index of the CC maximum, print the start times of the traces and print the time delay. Does the result meet your expectations from the figure of Task 1.5?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "32263d3c-e959-4644-9fee-20c42c8b5fdb",
+   "metadata": {},
+   "source": [
+    "#### Task 1.7: Plot the cross-correlation function"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0a9469e7-7d61-4d5e-8188-9259d982d60f",
+   "metadata": {},
+   "source": [
+    "Use the true delay time on the x-axis of your plot. "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a87c81e1-5c21-4b27-83e7-73712d7b2be3",
+   "metadata": {},
+   "source": [
+    "#### Task 1.8: Compute the energy and rms of the windowed data of each trace"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "33e4046e-b168-4105-a861-3e1d9c7f903f",
+   "metadata": {},
+   "source": [
+    "The energy of the trace is given by E=$\\int s^2(t)dt$. The square root of the energy is the root-mean-square (rms)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "a19552cd-14d4-4d60-a3d3-6ee13fe94f2b",
+   "metadata": {},
+   "source": [
+    "#### Task 1.9: Normalize the cross correlation function by the rms values of the correlated traces and plot it."
+   ]
+  }
+ ],
+ "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
+}

assignment_08/seismicTrace.py

@@ -0,0 +1,79 @@
+#----------------------------------------------------------------------------
+#  Python class for a seismic trace
+#
+import numpy as np
+
+class SeismicTrace(object):
+    """
+    Class defining a seismic trace by station name, network name, component,
+    number of samples, sampling interval and start time
+    """
+
+    def __init__(self, staname, comp, dt, tanf, vals):
+        """
+        Initialize seismic trace object from data:
+        :param staname: station name
+        :param comp: component
+        :param dt: sampling interval
+        :param tanf: start time
+        :param vals: values of trace (1D numpy array)
+        """
+        self.staname = staname
+        self.comp = comp
+        self.dt = dt
+        self.tanf = tanf
+        self.vals = vals
+        self.nsamp = vals.size
+        self.tlen = self.dt*self.nsamp
+
+
+    @classmethod
+    def readFromAscii(cls, fid, comp='None'):
+        """
+        Class method as second constructor to initialize trace from an Ascii file
+        :param fid: file descriptor of already open Ascii file
+        :param comp: optional component selection
+        :return: an instance of the class
+        """
+        staname = fid.readline().strip()
+        if not staname:
+            print("No station name found. Is file really an Ascii gather file")
+            return None
+        compf = fid.readline().strip()
+        if comp != 'None' and comp != compf:
+            print("Requested component %s not found.", comp)
+            return None
+        h = fid.readline().strip().split()
+        vals = np.array([float(v) for v in fid.readline().strip().split()])
+        dt = float(h[1]); tanf = float(h[2])
+        return cls(staname, compf, dt, tanf, vals)
+
+
+    def absoluteMaximum(self):
+        """
+        return absolute maximum of trace
+        """
+        return abs(max(self.vals))
+
+
+    def writeToOpenAsciiFile(self, fid):
+        """
+        write a seismic trace to an open Ascii file
+        """
+        fid.write(self.staname + '\n')
+        fid.write(self.comp + '\n')
+        fid.write("{:12d} {:15.9f} {:20.9f}\n".format(self.nsamp, self.dt, self.tanf))
+        for n in range(self.nsamp):
+            fid.write("{:17.8e}".format(self.vals[n]))
+        fid.write("\n")
+
+
+    def printInfo(self):
+        """
+        Print out some information about the trace
+        """
+        print("Identification: ", self.staname, self.comp)
+        print("Sampling interval: ", self.dt)
+        print("Number of samples: ", self.nsamp)
+        print("Start time after midnight: ", self.tanf)
+        print("Length of trace: (N*DT) ", self.tlen)

assignment_08/setupFigure.py

@@ -0,0 +1,15 @@
+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

assignment_10/multipleFilterTechnique.ipynb

@@ -0,0 +1,320 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Multiple filter technique"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "An alternative to the moving window method is the multiple filter technique. Instead of moving a window in the time domain, a window is moved in the frequency domain. Of course, a window in the frequency domain is a filter, and hence the name of the method. Let $H_k(\\omega,\\omega_k)$ be the transfer function of the $k$-th filter centered on frequency $\\omega_k$. Applying the filter involves multiplication of the filter transfer function with the Fourier transform of the time signal and then performing an inverse Fourier transform:\n",
+    "\\begin{align}\n",
+    "f_k(t) = \\int_{-\\infty}^\\infty\\,H_k(\\omega,\\omega_k)F(\\omega)\\exp(i\\omega t)\\frac{d\\omega}{2\\pi} \\,.\n",
+    "\\end{align}\n",
+    "Basically, one could plot the resulting time series against the associated center frequency of the filter. However, $f_k(t)$ can also assume negative values. It would be better to plot the envelope of $f_k(t)$.\n",
+    "\n",
+    "### The analytic signal\n",
+    "This can be obtained by constructing the analytic signal $z_k(t)$ of $f_k(t)$. It is a complex\n",
+    "signal whose real part is the original signal and whose imaginary part is its Hilbert transform. Its absolute\n",
+    "value is the envelope or instantaneous amplitude and its phase is the instantaneous phase:\n",
+    "\\begin{align}\n",
+    "z_k(t) = a_k(t)\\exp(i\\Phi_k(t)) = f_k(t)+iq_k(t) \\ \\mathrm{with}\\ q_k(t) = HT\\{f_k(t)\\}\\,,\n",
+    "\\end{align}\n",
+    "and\n",
+    "\\begin{align}\n",
+    "a_k(t) = \\sqrt{f_k(t)^2+q_k(t)^2},\\ \\Phi_k(t) = \\arctan\\left(\\frac{q_k(t)}{f_k(t)}\\right) \\,.\n",
+    "\\end{align}\n",
+    "\n",
+    "The Hilbert transform can again easily be obtained because its Fourier transform is:\n",
+    "\\begin{align}\n",
+    "Q_k(\\omega) = i\\,\\mathrm{sign}(\\omega)F_k(\\omega)\\,.\n",
+    "\\end{align}\n",
+    "where $F_k(\\omega)$ is the Fourier transform of $f_k(t)$. Inverse Fourier transform of $Q_k(\\omega)$ gives the Hilbert transform of $f_k(t)$.\n",
+    "\n",
+    "### Work flow\n",
+    "So the work flow of the multiple filter technique is as follows:\n",
+    "1. Compute Fourier transform of $f(t)$: $F(\\omega)$.\n",
+    "2. Multiply $F(\\omega)$ with bandpass filters $H_k(\\omega,\\omega_k)$ to obtain $F_k(\\omega)$.\n",
+    "3. Compute Fourier transform of the Hilbert transform: $Q_k(\\omega)=iF_k(\\omega)$.\n",
+    "4. Transform both back to the time domain.\n",
+    "5. Calculate the instantaneous amplitude $a_k(t)$ and phase $\\Phi_k(t)$.\n",
+    "6. Plot $a_k(t)$ in the time-frequency plane.\n",
+    "\n",
+    "### Gaussian filters\n",
+    "One frequently made choice for the bandpass filters is a Gaussian function of the form:\n",
+    "\\begin{align}\n",
+    "H_k(\\omega,\\omega_k) = e^{-\\alpha\\left(\\frac{\\omega-\\omega_k}{\\omega_k}\\right)^2}\n",
+    "\\end{align}\n",
+    "with inverse Fourier transform (impulse response):\n",
+    "\\begin{align}\n",
+    "h_k(t) = \\frac{\\sqrt{\\pi}\\omega_k}{2\\alpha}e^{-\\frac{\\omega_k^2t^2}{4\\alpha}}\\cos\\omega_k t \\,.\n",
+    "\\end{align}\n",
+    "This Fourier transform pair demonstrates a fundamental property of spectral analysis:\n",
+    "Choosing $\\alpha$ large gives a very narrow-banded filter but a very\n",
+    "long impulse response. Choosing $\\alpha$ small gives a wide-banded filter but a very short \n",
+    "impulse response. This is an expression of the uncertainty principle. If the frequency of a signal\n",
+    "is known very well, the signal is distributed over time. On the other hand, an impulse in the time domain has\n",
+    "a constant spectrum, i.e. the frequency is entirely unknown. In quantum mechanics, there is the energy-time\n",
+    "uncertainty relation, where energy $E=\\hbar\\omega$.\n",
+    "\n",
+    "### Time and frequency resolution\n",
+    "The Gaussian Fourier transform pair has a special property. If we measure the duration of the impulse response by\n",
+    "\\begin{align}\n",
+    "D_t^2 = \\int_{-\\infty}^\\infty t^2 |h(t)|^2 dt\n",
+    "\\end{align}\n",
+    "and the spread of the filter transfer function by\n",
+    "\\begin{align}\n",
+    "D_\\omega^2 = \\int_{-\\infty}^\\infty \\omega^2 |H(\\omega)|^2 \\frac{d\\omega}{2\\pi}\\,,\n",
+    "\\end{align}\n",
+    "then the product\n",
+    "\\begin{align}\n",
+    "D_t D_\\omega \\ge \\frac{1}{2} \\,.\n",
+    "\\end{align}\n",
+    "For the Gaussian filter, equality holds. Hence, the Gaussian filter is the one with the\n",
+    "best compromise between bandwidth and duration."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "# Preparation: load packages\n",
+    "%matplotlib inline\n",
+    "from obspy.core import read\n",
+    "from obspy.core import UTCDateTime\n",
+    "import numpy as np\n",
+    "import matplotlib.pylab as plt\n",
+    "from setupFigure import SetupFigure"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "def multipleFilterAnalysis(data, alfa, cfreq, dt, ndec):\n",
+    "    \"\"\"\n",
+    "    Perform a multiple filter analysis of data.\n",
+    "    :param data:  Array of detrended and demeaned data \n",
+    "    :param alfa:  Width parameter of Gaussian bandpass filter\n",
+    "    :param cfreq: Array of center frequencies of Gaussian filter\n",
+    "    :param dt: sampling interval of data\n",
+    "    :param ndec: decimation factor for instantaneous amplitude output\n",
+    "    \"\"\"\n",
+    "    nd = len(data)\n",
+    "    ftd = np.fft.rfft(data,nd)                   # Fourier transform of entire data set (pos. freq.)\n",
+    "    freq = np.fft.rfftfreq(nd,dt)                # Fourier frequencies (positive frequencies)\n",
+    "    jf = 0                                       # center frequency counter\n",
+    "    mfa = np.zeros((len(cfreq),nd//ndec))        # numpy array for MFA result\n",
+    "    for cf in cfreq:\n",
+    "        hg = np.exp(-alfa*((freq-cf)/cf)**2)     # Gaussian filter (use f instead of omega here)\n",
+    "        fk = hg*ftd                              # multiply FT of data with filter\n",
+    "        qk = 1j*fk                               # FT of Hilbert transform of filtered data\n",
+    "        ftk = np.fft.irfft(fk)                   # filtered data\n",
+    "        qtk = np.fft.irfft(qk)                   # Hilbert transform of filtered data\n",
+    "        at = np.sqrt(ftk**2+qtk**2)              # instantaneous amplitude\n",
+    "        mfa[jf,:] = at[0:nd:ndec]                # store decimated original result\n",
+    "        jf = jf+1                                # increase center frequency count\n",
+    "    return mfa"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def normalizeMFT(mfa, mode = 'full', exp = 0.5):\n",
+    "    \"\"\"\n",
+    "    Normalize the result of the mutiple filtering operation.\n",
+    "    :param mfa: array with envelopes versus frequency (mfa(f,t))\n",
+    "    :param mode: normalization mode: \n",
+    "                 if 'time': normalize timewise\n",
+    "                 if 'freq': normalize frequencywise\n",
+    "                 else: normalize to absolute maximum of entire array\n",
+    "    :param exp:  exponent for modifying envelope using a power\n",
+    "    \"\"\"\n",
+    "    nf, nt = mfa.shape\n",
+    "    if mode == 'freq':\n",
+    "        for jf in range(0,nf):\n",
+    "            mfamax = np.amax(mfa[jf,:])\n",
+    "            mfa[jf,:] = np.power(mfa[jf,:]/mfamax+1.e-10, exp)          \n",
+    "    elif mode == 'time':\n",
+    "        for jt in range(0,nt):\n",
+    "            mfamax = np.amax(mfa[:,jt])\n",
+    "            mfa[:,jt] = np.power(mfa[:,jt]/mfamax+1.e-10, exp)\n",
+    "    else:        \n",
+    "            mfamax = np.amax(mfa)\n",
+    "            mfa = np.power(mfa/mfamax+1.e-10, exp)\n",
+    "    #        mfa = np.log10(mfa/mfamax+1.e-10)\n",
+    "    return mfa"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "record = '1400'                                     # either use record starting at 1300 or 1400\n",
+    "datapath = 'geo2.' + record+ '.HHZ'\n",
+    "\n",
+    "if record =='1300':\n",
+    "    stime = UTCDateTime('2013-01-23 13:16:00Z')      # use record starting at 1300\n",
+    "    etime = UTCDateTime('2013-01-23 13:25:00Z')\n",
+    "    ttitle = ', depth: 36.5 m'\n",
+    "else:\n",
+    "    stime = UTCDateTime('2013-01-23 14:14:00Z')      # use record starting at 1400 (14:14-14:23)\n",
+    "    etime = UTCDateTime('2013-01-23 14:23:00Z')\n",
+    "    ttitle = ', depth: 56.5 m'\n",
+    "    \n",
+    "st = read(datapath)                                  # read file using obspy read\n",
+    "\n",
+    "st.trim(stime,etime)      # trim data stream to desired start and end time\n",
+    "st.detrend('linear')      # do a linear detrending\n",
+    "st.detrend('demean')      # subtract mean value"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "tr = st[0]                # extract first trace from stream (there is only one)\n",
+    "tr.plot()                 # plot trace"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# prepare multiple filtering\n",
+    "alfa = 10.                                     # filter width parameter\n",
+    "nd = tr.stats.npts                                # number of samples in time series\n",
+    "dt = tr.stats.delta                               # sampling interval\n",
+    "fmax = 1./(2.*dt)                                 # Nyquist frequency\n",
+    "cfreq = np.linspace(1., fmax, 1000)               # array of filter center frequencies\n",
+    "freq = np.fft.rfftfreq(nd, dt)                    # Fourier frequencies"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "ndec = 100                                        # decimation factor for envelope\n",
+    "mfa = multipleFilterAnalysis(tr.data, alfa, cfreq, dt, ndec)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mfa = normalizeMFT(mfa, mode='full', exp = 0.5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot the result\n",
+    "extent = (0., nd*dt, fmax, 1.)                         # extent of matrix in true time and frequency\n",
+    "fig1, ax1 = SetupFigure(15, 6, \"Time [s]\", \"Center frequency [Hz]\", \"Multiple filter analysis\", 14)\n",
+    "ax1.imshow(mfa, extent = extent, aspect=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot Gaussian filter\n",
+    "cf = 30\n",
+    "fig2, ax2 = SetupFigure(15, 6, \"Center frequency [Hz]\", \"Amplitude\", \"Gaussian filter\", 14)\n",
+    "arg = (freq-cf)/cf\n",
+    "hg = np.exp(-alfa*((freq-cf)/cf)**2)     # Gaussian filter (use f instead of omega here)\n",
+    "ax2.plot(freq, hg)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot filter response\n",
+    "fig3, ax3 = SetupFigure(15, 6, \"Time [s]\", \"Amplitude\", \"Filter response\", 14)\n",
+    "wk = cf*2.*np.pi\n",
+    "ts = np.linspace(-2.5,2.5,5000)\n",
+    "hres = wk/np.sqrt(np.pi*alfa)*np.exp(-(0.25*(wk*ts)**2/alfa))*np.cos(wk*ts) \n",
+    "ax3.plot(ts, hres)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tasks"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "1. Study the effect of $\\alpha$ on the result of the multiple filtering and on the form of the filter impulse response.\n",
+    "2. Up to which value of $\\alpha$ do you need to go to achieve a frequency resolution comparable to the moving window method\n",
+    "3. Think about the numerical effort of MFT in comparison with the moving window technique.\n",
+    "4. Try different scalings of the MFT result using different powers and also the logarithm.\n",
+    "5. Try timewise and frequenc-wise normalization.\n",
+    "6. Do you have some ideas why the results look different? What role play scaling and normalization? \n",
+    "7. Plot the Gaussian filter transfer function for different values of $\\alpha$.\n",
+    "8. Plot the filter responses for different values of $\\alpha$."
+   ]
+  }
+ ],
+ "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": 4
+}

assignment_10/setupFigure.py

@@ -0,0 +1,15 @@
+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

assignment_11/p11_surface_wave_group_velocity.ipynb

@@ -0,0 +1,171 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "from obspy import *\n",
+    "from obspy.clients.fdsn import Client\n",
+    "from obspy.clients.fdsn import URL_MAPPINGS\n",
+    "import numpy as np\n",
+    "import matplotlib.pylab as plt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Task 1: Copy functions for the Hann window and the moving window analysis from previous assignments here"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Your function for a Hann window"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Your function for the moving window analysis"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Data from data archive using obspy webfdsn client\n",
+    "FDSN stands for the International Federal of Digital Seismic Networks (www.fdsn.org). Citation from their web site: \"The International Federation of Digital Seismograph Networks (FDSN) is a global organization. Its membership is comprised of groups responsible for the installation and maintenance of seismographs either within their geographic borders or globally. Membership in the FDSN is open to all organizations that operate more than one broadband station. Members agree to coordinate station siting and provide free and open access to their data. This cooperation helps scientists all over the world to further the advancement of earth science and particularly the study of global seismic activity.\"\n",
+    "BGR (Bundesanstalt für Geowissenschaften und Rohstoffe) operates broadband stations in Germany and operates a data archive. It is member of the FDSN. That is why we can freely download data from them.\n",
+    "There are other data archives around the world which are also member of the FDSN and allow opn access to their data (see OBSPY documentation "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "for key in sorted(URL_MAPPINGS.keys()):                       # eine Liste der Archive\n",
+    "    print(\"{0:<7} {1}\".format(key,  URL_MAPPINGS[key]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Get waveforms for the Tohoku earthquake for station TNS of the German Regional Seismic Network."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "client = Client(\"BGR\")                                          # data archive from which data are fetched\n",
+    "t1 = UTCDateTime(\"2018-01-14T09:19:00.000\")                     # desired start time of data\n",
+    "stc = client.get_waveforms(\"GR\", \"TNS\",\"\",\"LHZ\",t1,             # use fdsn web client to download data\n",
+    "                           t1 + 7200.,attach_response = True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "stc.detrend('linear')                       # take out linear trend\n",
+    "stc.plot()\n",
+    "print(stc)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "st = stc.copy()                                     # proceed with copy to avoid repeated downloading\n",
+    "ta = UTCDateTime(\"2018-01-14T10:00:00.000000Z\")     # cut to surface wave train\n",
+    "te = UTCDateTime(\"2018-01-14T10:29:59.000000Z\")\n",
+    "st.trim(ta,te)\n",
+    "st.taper(0.05,\"hann\")                               # apply a taper to bring signal to zero at ends\n",
+    "st.plot()                                           # dispersion is now well visible\n",
+    "print(st)\n",
+    "tr = st[0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Task 2: Apply the moving window analysis to the data. \n",
+    "Think about window length and shift time. Try different window lengths."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Task 3: Plot the moving window matrix and the time series below. \n",
+    "Use the same limits for the time axis. Choose a subrange for the frequencies. Does the time-frequency analysis reflect the properties of the time series?"
+   ]
+  }
+ ],
+ "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": 4
+}