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

removed solutions of 03 and 04 from repository

Browse Source
This commit is contained in:
Wolfgang Friederich 2025-05-09 12:37:28 +02:00
parent b15de68602
commit 2feb219487
2 changed files with 0 additions and 528 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

225
assignment_03/gibbs_incremental_reconstruction.ipynb
View File

@ -1,225 +0,0 @@
{
"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
}

303
assignment_04/fft_sampling_interpol.ipynb
View File

@ -1,303 +0,0 @@
{
"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": [
"t, bx = boxcar(1, 1000, 300, 500)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fe9a16bd-1f5c-4a28-9b4b-ef7274c5f529",
"metadata": {},
"outputs": [],
"source": [
"c = dft_coeff(bx) # slow DFT\n",
"fc = dft_fast_coeff(bx) # fast DFT\n",
"print(t.size, c.size, fc.size)\n",
"for n in range(t.size//50):\n",
" print(\"{:6d} {:15.6e} {:15.6e}\".format(n, c[n], fc[n]))"
]
},
{
"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": [
"xr = dft_synthesis(c)\n",
"if t.size%2 == 0:\n",
" xrf = dft_fast_synthesis(fc, outnum='even')\n",
"else:\n",
" xrf = dft_fast_synthesis(fc, outnum='odd')\n",
"print(xr.size, xrf.size)\n",
"for k in range(290,311):\n",
" print(\"{:6d} {:15.6e} {:15.6e}\".format(k, xr[k], xrf[k]))"
]
},
{
"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, hw=20. \n",
"\n",
"Second, append some zeros (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": [
"dt = 5; period = 100; tmax = 50; hw = 20; nz = 20\n",
"t, g = gaussian(dt, period, tmax, hw)\n",
"c = dft_fast_coeff(g)\n",
"print(\"N = \", t.size,\" Nyquist frequency = \", t.size/(2*period),\" Ncoeff = \", c.size)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2ad7ff56-8713-4386-9b4a-9228fb63fde0",
"metadata": {},
"outputs": [],
"source": [
"cz = np.append(c, np.zeros(nz))\n",
"fny = (t.size/2+nz)/period\n",
"dtnew = 1./(2*fny)\n",
"if t.size%2 != 0:\n",
" xr = dft_fast_synthesis(cz, outnum='odd')\n",
"else:\n",
" xr = dft_fast_synthesis(cz, outnum='even')\n",
"tz = dtnew*np.arange(0, xr.size)\n",
"print(\"N = \", tz.size,\" Nyquist frequency = \", fny)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2755badf-b7d1-47f0-ab3b-1a057faaf921",
"metadata": {},
"outputs": [],
"source": [
"fig1, ax1 = SetupFigure(10, 5, \"Time [s]\", \"xr(t)\", \"Interpolated Gaussian\", 14)\n",
"ax1.plot(tz, xr, color='blue', ls='-', lw=0.25, marker='d', markersize=4, zorder=10)\n",
"ax1.plot(t, g, color='red', ls='', marker='o', markersize=6, zorder=5)"
]
},
{
"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 sme 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": [
"t, g = gaussian(0.25, 100, 50, 1)\n",
"fig3, ax3 = SetupFigure(10, 3, \"Time [s]\", \"g(t)\", \"Sharp Gaussian\", 14)\n",
"ax3.plot(t, g, color='blue', ls='-', lw=1.0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "83511cfe-45e9-4216-96c7-f932d6199034",
"metadata": {},
"outputs": [],
"source": [
"fig2, ax2 = SetupFigure(10, 5, \"Frequency [s]\", \"G(f)\", \"Fourier coefficients\", 14)\n",
"ax2.set_xlim(0, 1.0)\n",
"col = ['grey', 'orange', 'red', 'magenta', 'blue']\n",
"lw = [4, 1, 1, 1, 1]\n",
"for i, dt in enumerate([0.25, 0.5, 1.0, 1.5, 2.0]):\n",
" t, g = gaussian(dt, 100, 50, 1)\n",
" c = dft_fast_coeff(g)\n",
" f = 1./period*np.arange(0, c.size)\n",
" ax2.plot(f, np.absolute(c), color=col[i], ls='-', lw=lw[i], label=\"DT = {:5.2f}\".format(dt))\n",
"ax2.legend() "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b8b80cd5-d9b7-40e3-858c-ed2086167c37",
"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
}