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

add new Jupyter Notebook for exercise onf fourier transsform

Browse Source
This commit is contained in:
Kasper D. Fischer 2025-05-05 12:23:37 +02:00
parent 56d10e31b9
commit 35de751e04
1 changed files with 562 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

562
02-Fourier_Transformation_and_Series/Fourier_Transformation_and_Series.ipynb Normal file
View File

@ -0,0 +1,562 @@
{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"<div style='background-image: url(\"../images/header.svg\") ; padding: 0px ; background-size: cover ; border-radius: 5px ; height: 250px'>\n",
" <div style=\"float: right ; margin: 50px ; padding: 20px ; background: rgba(255 , 255 , 255 , 0.7) ; width: 50% ; height: 150px\">\n",
" <div style=\"position: relative ; top: 50% ; transform: translatey(-50%)\">\n",
" <div style=\"font-size: xx-large ; font-weight: 900 ; color: rgba(0 , 0 , 0 , 0.8) ; line-height: 100%\">Signal Processing</div>\n",
" <div style=\"font-size: large ; padding-top: 20px ; color: rgba(0 , 0 , 0 , 0.5)\">Fourier Transformation</div>\n",
" </div>\n",
" </div>\n",
"</div>\n",
"\n",
"Seismo-Live: http://seismo-live.org\n",
"\n",
"##### Authors:\n",
"* Stefanie Donner ([@stefdonner](https://github.com/stefdonner))\n",
"* Celine Hadziioannou ([@hadzii](https://github.com/hadzii))\n",
"* Ceri Nunn ([@cerinunn](https://github.com/cerinunn))\n",
"\n",
"Some code used in this tutorial is taken from [stackoverflow.com](http://stackoverflow.com/questions/4258106/how-to-calculate-a-fourier-series-in-numpy/27720302#27720302). We thank [Giulio Ghirardo](https://www.researchgate.net/profile/Giulio_Ghirardo) for his kind permission to use his code here.\n",
"\n",
"---"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"code_folding": [
0
],
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"# Cell 0 - Preparation: load packages, set some basic options \n",
"%matplotlib inline\n",
"from scipy import signal\n",
"from obspy.signal.invsim import cosine_taper \n",
"from matplotlib import rcParams\n",
"import numpy as np\n",
"import matplotlib.pylab as plt\n",
"plt.style.use('ggplot')\n",
"plt.rcParams['figure.figsize'] = 15, 3"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"<h1>Tutorial on Fourier transformation in 1D </h1>\n",
"<br>"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## The Fourier transformation\n",
"\n",
"In the world of seismology, we use the *Fourier transformation* to transform a signal from the time domain into the frequency domain. That means, we split up the signal and separate the content of each frequency from each other. Doing so, we can analyse our signal according to energy content per frequency. We can extract information on how much amplitude each frequency contributes to the final signal. In other words: we get a receipt of the ingredients we need to blend our measured signal. "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"### Fourier series\n",
"\n",
"The *Fourier transformation* is based on the *Fourier series*. With the *Fourier series* we can approximate an (unknown) function $f(x)$ by another function $g_n(x)$ which consists of a sum over $N$ basis functions weighted by some coefficients. The basis functions need to be orthogonal. $sin$ and $cos$ functions seem to be a pretty good choice because any signal can be filtered into several sinusoidal paths. In the period range of $[-T/2 ; T/2]$ the *Fourier series* is defined as:\n",
"\n",
"$$\n",
"f(t) \\approx g_n(t) = \\frac{1}{2} a_0 + \\sum_{k=1}^N \\left[ a_k \\cos \\left(\\frac{2\\pi k t}{T} \\right) + b_k \\sin\\left(\\frac{2\\pi k t}{T}\\right)\\right]\n",
"$$\n",
"\n",
"$$ \n",
"a_k = \\frac{2}{T} \\int_{-T/2}^{T/2} f(t) \\cos\\left(\\frac{2\\pi k t}{T}\\right)dt\n",
"$$\n",
"\n",
"$$\n",
"b_k = \\frac{2}{T} \\int_{-T/2}^{T/2} f(t) \\sin\\left(\\frac{2\\pi k t}{T}\\right)dt\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Fourier transform of periodic functions\n",
"\n",
"At this stage, we consider continuous, periodic and infinite functions. The more basis functions are used to approximate the unknown function, the better is the approximation, i.e. the more similar the unknown function is to its approximation. \n",
"\n",
"For a non-periodic function the interval of periodicity tends to infinity. That means, the steps between neighbouring frequencies become smaller and smaller and thus the infinite sum of the *Fourier series* turns into an integral and we end up with the integral form of the *Fourier transformation*:\n",
"\n",
"$$\n",
"F(\\omega) = \\frac{1}{2\\pi} \\int_{-\\infty}^{\\infty} f(t) e^{-i\\omega t} dt \\leftrightarrow f(t) = \\int_{-\\infty}^{\\infty} F(\\omega)e^{i\\omega t}dt\n",
"$$\n",
"\n",
"Attention: sign and factor conventions can be different in the literature!"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"In seismology, we do not have continuous but discrete time signals. Therefore, we work with the discrete form of the *Fourier transformation*:\n",
"\n",
"$$\n",
"F_k = \\frac{1}{N} \\sum_{j=0}^{N-1} f_j e^{-2\\pi i k j /N} \\leftrightarrow f_k = \\sum_{j=0}^{N-1} F_j e^{2\\pi i k j /N}\n",
"$$\n",
"\n",
"Some intuitive gif animations on what the *Fourier transform* is doing, can be found [here](https://en.wikipedia.org/wiki/File:Fourier_series_and_transform.gif), [here](https://en.wikipedia.org/wiki/File:Fourier_series_square_wave_circles_animation.gif), and [here](https://en.wikipedia.org/wiki/File:Fourier_series_sawtooth_wave_circles_animation.gif).\n",
"Further and more detailed explanations on *Fourier series* and *Fourier transformations* can be found [here](https://betterexplained.com/articles/an-interactive-guide-to-the-fourier-transform/) and [here](www.fourier-series.com).\n",
"\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### The Fourier series and its coefficients \n",
"\n",
"In the following two code cells, we first define a function which calculates the coefficients of the Fourier series for a given function. The function in the next cell does it the other way round: it is creating a function based on given coefficients and weighting factors."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"code_folding": [
1
],
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"# Cell 1: code by Giulio Ghirardo \n",
"def fourier_series_coeff(f, T, N):\n",
" \"\"\"Calculates the first 2*N+1 Fourier series coeff. of a periodic function.\n",
"\n",
" Given a periodic, function f(t) with period T, this function returns the\n",
" coefficients a0, {a1,a2,...},{b1,b2,...} such that:\n",
"\n",
" f(t) ~= a0/2+ sum_{k=1}^{N} ( a_k*cos(2*pi*k*t/T) + b_k*sin(2*pi*k*t/T) )\n",
" \n",
" Parameters\n",
" ----------\n",
" f : the periodic function, a callable like f(t)\n",
" T : the period of the function f, so that f(0)==f(T)\n",
" N_max : the function will return the first N_max + 1 Fourier coeff.\n",
"\n",
" Returns\n",
" -------\n",
" a0 : float\n",
" a,b : numpy float arrays describing respectively the cosine and sine coeff.\n",
" \"\"\"\n",
" # From Nyquist theorem we must use a sampling \n",
" # freq. larger than the maximum frequency you want to catch in the signal. \n",
" f_sample = 2 * N\n",
" \n",
" # We also need to use an integer sampling frequency, or the\n",
" # points will not be equispaced between 0 and 1. We then add +2 to f_sample.\n",
" t, dt = np.linspace(0, T, f_sample + 2, endpoint=False, retstep=True)\n",
" y = np.fft.rfft(f) / f.size\n",
" y *= 2\n",
" return y[0].real, y[1:-1].real[0:N], -y[1:-1].imag[0:N]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"code_folding": [
1
],
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"# Cell 2: code by Giulio Ghirardo \n",
"def series_real_coeff(a0, a, b, t, T):\n",
" \"\"\"calculates the Fourier series with period T at times t,\n",
" from the real coeff. a0,a,b\"\"\"\n",
" tmp = np.ones_like(t) * a0 / 2.\n",
" for k, (ak, bk) in enumerate(zip(a, b)):\n",
" tmp += ak * np.cos(2 * np.pi * (k + 1) * t / T) + bk * np.sin(\n",
" 2 * np.pi * (k + 1) * t / T)\n",
" return tmp"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Now, we can create an arbitrary function, which we use to experiment with in the following example. \n",
"1) When you re-run cell 3 several times, do you always see the same function? Why? What does it tell you about the Fourier series? "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"# Cell 3: create periodic, discrete, finite signal\n",
"\n",
"# number of samples (intial value: 3000)\n",
"samp = 3000\n",
"# sample rate (initial value: 1)\n",
"dt = 1\n",
"# period\n",
"T = 1.0 / dt\n",
"length = samp * dt\n",
"# number of coefficients (initial value: 100)\n",
"N = 100\n",
"# weighting factors for coefficients (selected randomly)\n",
"a0 = np.random.rand(1)\n",
"a = np.random.randint(low=1, high=11, size=N)\n",
"b = np.random.randint(low=1, high=11, size=N)\n",
"\n",
"# \n",
"a0= 0.0\n",
"a = np.exp(-np.linspace(0, 5, N))\n",
"b = np.zeros(N)\n",
"\n",
"t = np.linspace(0, length, samp) # time axis\n",
"sig = series_real_coeff(a0, a, b, t, T)\n",
"\n",
"# plotting\n",
"plt.plot(t, sig, 'r', label='arbitrary, periodic, discrete, finite signal')\n",
"plt.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
"plt.xlabel('time [sec]')\n",
"plt.ylabel('amplitude')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"Now, we can play with the signal and see what happens when we try to reconstruct it with a limited number of coefficients. \n",
"2) Run the cells 4 and 5. What do you observe? \n",
"3) Increase the number of coefficients $n$ step by step and re-run cells 4 and 5. What do you observe now? Can you explain? \n",
"4) In cell 5 uncomment the lines to make a plot which is not normalized (and comment the other two) and re-run the cell. What do you see now and can you explain it?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"# Cell 4: determine the first 'n' coefficients of the function using the code function of cell 1\n",
"T = 1 # period\n",
"n = 100 # number of coeffs to reconstruct\n",
"a0, a, b = fourier_series_coeff(sig, T, n)\n",
"a_ = a.astype(int)\n",
"b_ = b.astype(int)\n",
"print('coefficient a0 = ', int(a0))\n",
"print('array coefficients ak =', a_)\n",
"print('array coefficients bk =', b_)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"# Cell 5: reconstruct the function using the code in cell 2\n",
"g = series_real_coeff(a0, a, b, t, dt)\n",
"\n",
"# plotting\n",
"#plt.plot(t, sig, 'r', label='original signal') # NOT normalized \n",
"#plt.plot(t, g, 'g', label='reconstructed signal')\n",
"plt.plot(t, sig, 'r', label='original signal') # normalized \n",
"plt.plot(t, g, 'g', label='reconstructed signal')\n",
"\n",
"plt.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
"plt.xlabel('time [sec]')\n",
"plt.ylabel('amplitude')\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Fourier series, convergence and Gibb's phenomenon\n",
"\n",
"As seen above the convergence of the *Fourier series* can be tricky due to the fact that we work with signals of finite length. To analyse this effect in a bit more detail, we define a square wave in cell 6 and try to reconstruct it in cell 7. \n",
"5) First, we use only 5 coefficients to reconstruct the wave. Describe what you see. \n",
"6) Increase the number of coefficients $n$ in cell 7 step by step and re-run the cell. What do you see now? Can you explain it?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"# Cell 6: define a square wave of 5 Hz\n",
"freq = 5.\n",
"npts = 3000\n",
"dt_ = 0.002\n",
"length = npts * dt_\n",
"t_ = np.linspace(0, length, npts, endpoint=False)\n",
"square = signal.square(2 * np.pi * freq * t_)\n",
"\n",
"plt.plot(t_, square)\n",
"plt.xlabel('time [sec]')\n",
"plt.ylabel('amplitude')\n",
"plt.xlim(0, 1.05)\n",
"plt.ylim(-1.2, 1.2)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"source": [
"##### You may replace the square function by something else. What about a sawtooth function?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"signal.sawtooth?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"# Cell 7: reconstruct signal using convergence criterion\n",
"n = 5 # number of coefficients (initial: 5)\n",
"T_ = 1/freq # period of signal\n",
"\n",
"# determine coefficients\n",
"a0 = 0\n",
"a = []\n",
"b = []\n",
"for i in range(1,n):\n",
" if (i%2 != 0):\n",
" a_ = 4/(np.pi*i)\n",
" else:\n",
" a_ = 0\n",
" a.append(a_)\n",
" b_ = (2*np.pi*i)/T_\n",
" b.append(b_)\n",
"\n",
"# reconstruct signal\n",
"g = np.ones_like(t_) * a0\n",
"for k, (ak, bk) in enumerate(zip(a, b)):\n",
" g += ak * np.sin(bk*t_)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"# Cell 7b: plotting\n",
"plt.plot(t_, square, 'r', label='Analytisches Signal') \n",
"plt.plot(t_, g, 'g', label='Reihenentwicklung')\n",
"plt.ticklabel_format(axis='y', style='sci', scilimits=(-1,1))\n",
"plt.xlabel('time [sec]')\n",
"plt.ylabel('amplitude')\n",
"#plt.ylim(0.9,1.1)\n",
"plt.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"### Fourier transformation\n",
"\n",
"Let us now do the Fourier transformation of the signal created in cell 3 and have a look on the amplitude spectra. In computer science the transformation is performed as fast Fourier transformation (FFT). \n",
"\n",
"7) Why do we need to taper the signal before we perform the FFT? \n",
"8) How do you interpret the plot of the amplitude spectra? \n",
"9) Which frequency contributes most to the final signal? "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"# Cell 8: FFT of signal\n",
"# number of sample points need to be the same as in cell 3\n",
"print('samp =',samp,' Need to be the same as in cell 3.')\n",
"# number of sample points need to be the same as in cell 3\n",
"print('T =',T,' Need to be the same as in cell 3.')\n",
"# percentage of taper applied to signal (initial: 0.1)\n",
"taper_percentage = 0.1\n",
"taper = cosine_taper(samp,taper_percentage)\n",
"\n",
"sig_ = square * taper\n",
"Fsig = np.fft.rfft(sig_, n=samp)\n",
"\n",
"# prepare plotting\n",
"xf = np.linspace(0.0, 1.0/(2.0*T), (samp//2)+1)\n",
"rcParams[\"figure.subplot.hspace\"] = (0.8)\n",
"rcParams[\"figure.figsize\"] = (15, 9)\n",
"rcParams[\"axes.labelsize\"] = (15)\n",
"rcParams[\"axes.titlesize\"] = (20)\n",
"rcParams[\"font.size\"] = (12)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"#Cell 8b: plotting\n",
"plt.subplot(311)\n",
"plt.title('Time Domain')\n",
"plt.plot(t, square, linewidth=1)\n",
"plt.xlabel('Time [s]')\n",
"plt.ylabel('Amplitude')\n",
"\n",
"plt.subplot(312)\n",
"plt.title('Frequency Domain')\n",
"plt.plot(xf, 2.0/npts * np.abs(Fsig))\n",
"#plt.xlim(0, 0.04) \n",
"plt.xlabel('Frequency [Hz]')\n",
"plt.ylabel('Amplitude')\n",
"plt.show()"
]
}
],
"metadata": {
"anaconda-cloud": {},
"celltoolbar": "Slideshow",
"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.9.21"
}
},
"nbformat": 4,
"nbformat_minor": 4
}