seismology-RUB/SeismologicalDataAnalysis2025
10
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: seismology-RUB:b3db046a3b929c60304e1ca5c4ceefa7bae35d32
Branches Tags
seismology-RUB:main
...
compare: seismology-RUB:0dd1763a09308afff410a78661bcba964f115abc
Branches Tags
seismology-RUB:main

3 Commits
b3db046a3b ... 0dd1763a09

Author SHA1 Message Date
friederichwo
0dd1763a09 Merge branch 'main' of git.geophysik.ruhr-uni-bochum.de:seismology-RUB/SeismologicalDataAnalysis2025 2025-05-08 11:05:27 +02:00
friederichwo
2fa3e673d4 new files in assignment 3 2025-05-08 11:03:36 +02:00
friederichwo
3e7a875991 added notebooks to assignment 2 2025-05-07 11:13:09 +02:00
4 changed files with 609 additions and 0 deletions
Show Stats Expand all files Collapse all files

357
assignment_02/dft.ipynb Normal file
View File

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

16
assignment_02/setupFigure.py Normal file
View File

@ -0,0 +1,16 @@
import matplotlib.pyplot as plt
def SetupFigure(wx, wy, xlabel, ylabel, title, tls):
"""
Create a new figure frame
"""
fig = plt.figure(figsize=(wx,wy))
ax = fig.add_subplot(1,1,1)
ax.set_title(title, fontsize=tls+2)
ax.grid(which='major', axis='both', ls=':')
ax.set_ylabel(ylabel, fontsize=tls)
ax.set_xlabel(xlabel, fontsize=tls)
ax.tick_params(labelsize=tls)
return fig, ax

221
assignment_03/p3_gibbs_incremental_reconstruction.ipynb Normal file
View File

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

15
assignment_03/setupFigure.py Normal file
View File

@ -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