SeismologicalDataAnalysis2025/assignment_03/p3_gibbs_incremental_reconstruction.ipynb

222 lines
7.7 KiB
Plaintext

{
"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
}