new assignment 11

- Python Crash Course

.gitattributes

@@ -0,0 +1,3 @@
+// This file is used to define attributes for files in the repository.
+*.mseed filter=lfs diff=lfs merge=lfs -text
+*.sac filter=lfs diff=lfs merge=lfs -text

.gitignore

@@ -10,7 +10,7 @@ profile_default/
 ipython_config.py
 
 # Remove previous ipynb_checkpoints
-#   git rm -r .ipynb_checkpoints/
+git rm -r .ipynb_checkpoints/
 
 # ---> Python
 # Byte-compiled / optimized / DLL files

01-Python_Introduction/Python_Crash_Course.ipynb

@@ -0,0 +1,771 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "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%\">Scientific Python</div>\n",
+    "            <div style=\"font-size: large ; padding-top: 20px ; color: rgba(0 , 0 , 0 , 0.5)\">A super quick crash course</div>\n",
+    "        </div>\n",
+    "    </div>\n",
+    "</div>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Seismo-Live: http://seismo-live.org\n",
+    "\n",
+    "##### Authors:\n",
+    "* Lion Krischer ([@krischer](https://github.com/krischer))\n",
+    "\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This notebook is a very quick introduction to Python and in particular its scientific ecosystem in case you have never seen it before. It furthermore grants a possibility to get to know the [IPython/Jupyter notebook](http://www.nature.com/news/interactive-notebooks-sharing-the-code-1.16261). [See here for the official documentation](http://nbviewer.jupyter.org/github/jupyter/notebook/blob/master/docs/source/examples/Notebook/Notebook%20Basics.ipynb) of the Jupyter notebook - a ton more information can be found online.\n",
+    "\n",
+    "\n",
+    "A lot of motivational writing on *Why Python?* is out there so we will not repeat it here and just condense it to a single sentence: **Python is a good and easy to learn, open-source, general purpose programming language that happens to be very good for many scientific tasks (due to its vast scientific ecosystem).**\n",
+    "\n",
+    "\n",
+    "#### Quick Reference on How to Use This Notebook\n",
+    "\n",
+    "\n",
+    "<img src=\"../images/notebook_toolbar.png\" style=\"width:70%\"></img>\n",
+    "\n",
+    "* `Shift + Enter`: Execute cell and jump to the next cell\n",
+    "* `Ctrl/Cmd + Enter`: Execute cell and don't jump to the next cell\n",
+    "\n",
+    "\n",
+    "#### Disclaimer\n",
+    "\n",
+    "The tutorials are employing Jupyter notebooks but these are only one way of using Python. Writing scripts to text files and executing them with the Python interpreter of course also works:\n",
+    "\n",
+    "```bash\n",
+    "$ python do_something.py\n",
+    "```\n",
+    "\n",
+    "Another alternative is interactive usage on the command line:\n",
+    "\n",
+    "```bash\n",
+    "$ ipython\n",
+    "```\n",
+    "\n",
+    "## Notebook Setup\n",
+    "\n",
+    "First things first: In many notebooks you will find a cell similar to the following one. **Always execute it!** They do a couple of things:\n",
+    "* Make plots appear in the browser (otherwise a window pops up)\n",
+    "* Printing things works like this: \n",
+    "\n",
+    "```python\n",
+    "print(\"Hello\")\n",
+    "```\n",
+    "\n",
+    "This essentially makes the notebooks work under Python 2 and Python 3.\n",
+    "\n",
+    "* Plots look quite a bit nicer (this is optional).\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Plots now appear in the notebook.\n",
+    "%matplotlib inline                       \n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "plt.style.use('ggplot')                            # Matplotlib style sheet - nicer plots!\n",
+    "plt.rcParams['figure.figsize'] = 12, 8             # Slightly bigger plots by default"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "---\n",
+    "\n",
+    "## Useful Links\n",
+    "\n",
+    "Here is collection of resources regarding the scientific Python ecosystem. They cover a number of different packages and topics; way more than we will manage today.\n",
+    "\n",
+    "If you have any question regarding some specific Python functionality you can consult the official [Python documenation](http://docs.python.org/).\n",
+    " \n",
+    "Furthermore a large number of Python tutorials, introductions, and books are available online. Here are some examples for those interested in learning more.\n",
+    " \n",
+    "* [Learn Python The Hard Way](http://learnpythonthehardway.org/book/)\n",
+    "* [Dive Into Python](http://www.diveintopython.net/)\n",
+    "* [The Official Python Tutorial](http://docs.python.org/2/tutorial/index.html)\n",
+    "* [Think Python Book](http://www.greenteapress.com/thinkpython/thinkpython.html)\n",
+    " \n",
+    "Some people might be used to Matlab - this helps:\n",
+    " \n",
+    "* [NumPy for Matlab Users Introdution](http://wiki.scipy.org/NumPy_for_Matlab_Users)\n",
+    "* [NumPy for Matlab Users Cheatsheet](http://mathesaurus.sourceforge.net/matlab-numpy.html)\n",
+    " \n",
+    " \n",
+    "Additionally there is an abundance of resources introducing and teaching parts of the scientific Python ecosystem.\n",
+    " \n",
+    "* [NumPy Tutorial](http://wiki.scipy.org/Tentative_NumPy_Tutorial)\n",
+    "* [Probabilistic Programming and Bayesian Methods for Hackers](http://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/): Great ebook introducing Bayesian methods from an understanding-first point of view with the examples done in Python.\n",
+    "* [Python Scientific Lecture Notes](http://scipy-lectures.github.io/): Introduces the basics of scientific Python with lots of examples.\n",
+    "* [Python for Signal Processing](http://python-for-signal-processing.blogspot.de/): Free blog which is the basis of a proper book written on the subject.\n",
+    "* [Another NumPy Tutorial](http://www.loria.fr/~rougier/teaching/numpy/numpy.html), [Matplotlib Tutorial](http://www.loria.fr/~rougier/teaching/matplotlib/matplotlib.html)\n",
+    " \n",
+    "You might eventually have a need to create some custom plots. The quickest way to success is usually to start from some example that is somewhat similar to what you want to achieve and just modify it. These websites are good starting points:\n",
+    " \n",
+    "* [Matplotlib Gallery](http://matplotlib.org/gallery.html)\n",
+    "* [ObsPy Gallery](http://docs.obspy.org/gallery.html)\n",
+    "* [Basemap Gallery](http://matplotlib.org/basemap/users/examples.html)\n",
+    "\n",
+    "\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Core Python Crash Course\n",
+    "\n",
+    "This course is fairly non-interactive and serves to get you up to speed with Python assuming you have practical programming experience with at least one other language. Nonetheless please change things and play around an your own - it is the only way to really learn it!\n",
+    "\n",
+    "The first part will introduce you to the core Python language. This tutorial uses Python 3 but almost all things can be transferred to Python 2. If possible choose Python 3 for your own work!\n",
+    "\n",
+    "\n",
+    "### 1. Numbers\n",
+    "\n",
+    "Python is dynamically typed and assigning something to a variable will give it that type."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Three basic types of numbers\n",
+    "a = 1             # Integers\n",
+    "b = 2.0           # Floating Point Numbers\n",
+    "c = 3.0 + 4j      # Complex Numbers, note the use of j for the complex part\n",
+    "\n",
+    "\n",
+    "# Arithmetics work as expected.\n",
+    "# Upcasting from int -> float -> complex\n",
+    "d = a + b         # (int + float = float)\n",
+    "print(d)\n",
+    "\n",
+    "e = c ** 2        # c to the second power, performs a complex multiplication\n",
+    "print(e)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 2. Strings"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Just enclose something in single or double quotes and it will become a string. On Python 3 it defaults to unicode strings, e.g. non Latin alphabets and other symbols."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# You can use single or double quotes to create strings.\n",
+    "location = \"New York\"\n",
+    "\n",
+    "# Concatenate strings with plus.\n",
+    "where_am_i = 'I am in ' + location\n",
+    "\n",
+    "# Print things with the print() function.\n",
+    "print(location, 1, 2)\n",
+    "print(where_am_i)\n",
+    "\n",
+    "# Strings have a lot of attached methods for common manipulations.\n",
+    "print(location.lower())\n",
+    "\n",
+    "# Access single items with square bracket. Negative indices are from the back.\n",
+    "print(location[0], location[-1])\n",
+    "\n",
+    "# Strings can also be sliced.\n",
+    "print(location[4:])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Exercise\n",
+    "\n",
+    "Save your name in all lower-case letters to a variable, and print a capitalized version of it. Protip: [Google for \"How to capitalize a string in python\"](http://www.google.com/search?q=how+to+capitalize+a+string+in+python). This works for almost any programming problem - someone will have had the same issue before!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 3. Lists"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Python has two main collection types: List and dictionaries. The former is just an ordered collection of objects and is introduced here."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# List use square brackets and are simple ordered collections of things.\n",
+    "everything = [a, b, c, 1, 2, 3, \"hello\"]\n",
+    "\n",
+    "# Access elements with the same slicing/indexing notation as strings.\n",
+    "# Note that Python indices are zero based!\n",
+    "print(everything[0])\n",
+    "print(everything[:3])\n",
+    "print(everything[2:-2])\n",
+    "\n",
+    "# Negative indices are counted from the back of the list.\n",
+    "print(everything[-3:])\n",
+    "\n",
+    "# Append things with the append method.\n",
+    "everything.append(\"you\")\n",
+    "print(everything)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 4. Dictionaries\n",
+    "\n",
+    "The other main collection type in Python are dictionaries. They are similiar to associative arrays or (hash) maps in other languages. Each entry is a key-value pair."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Dictionaries have named fields and no inherent order. As is\n",
+    "# the case with lists, they can contain anything.\n",
+    "information = {\n",
+    "    \"name\": \"Hans\",\n",
+    "    \"surname\": \"Mustermann\",\n",
+    "    \"age\": 78,\n",
+    "    \"kids\": [1, 2, 3]\n",
+    "}\n",
+    "\n",
+    "# Acccess items by using the key in square brackets.\n",
+    "print(information[\"kids\"])\n",
+    "\n",
+    "# Add new things by just assigning to a key.\n",
+    "print(information)\n",
+    "information[\"music\"] = \"jazz\"\n",
+    "print(information)\n",
+    "\n",
+    "# Delete things by using the del operator\n",
+    "del information[\"age\"]\n",
+    "print(information)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "collapsed": true
+   },
+   "source": [
+    "### 5. Functions\n",
+    "\n",
+    "The key to conquer a big problem is to divide it into many smaller ones and tackle them one by one. This is usually achieved by using functions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Functions are defined using the def keyword.\n",
+    "def do_stuff(a, b):\n",
+    "    return a * b\n",
+    "\n",
+    "# And called with the arguments in round brackets.\n",
+    "print(do_stuff(2, 3))\n",
+    "\n",
+    "# Python function also can have optional arguments.\n",
+    "def do_more_stuff(a, b, power=1):\n",
+    "    return (a * b) ** power\n",
+    "\n",
+    "print(do_more_stuff(2, 3))\n",
+    "print(do_more_stuff(2, 3, power=3))\n",
+    "\n",
+    "# For more complex function it is oftentimes a good idea to \n",
+    "#explicitly name the arguments. This is easier to read and less error-prone.\n",
+    "print(do_more_stuff(a=2, b=3, power=3))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 6. Imports\n",
+    "\n",
+    "To use functions and objects not part of the default namespace, you have import them. You will have to do this a lot so it is necessary to learn how to do it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Import anything, and use it with the dot accessor.\n",
+    "import math\n",
+    "\n",
+    "a = math.cos(4 * math.pi)\n",
+    "\n",
+    "# You can also selectively import things.\n",
+    "from math import pi\n",
+    "\n",
+    "b = 3 * pi\n",
+    "\n",
+    "# And even rename them if you don't like their name.\n",
+    "from math import cos as cosine\n",
+    "c = cosine(b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "How to know what is available?\n",
+    "\n",
+    "1. Read the [documentation](https://docs.python.org/3/library/math.html)\n",
+    "2. Interactively query the module"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "print(dir(math))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Typing the dot and the TAB will kick off tab-completion."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "math."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In the IPython framework you can also use a question mark to view the documentation of modules and functions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "math.cos?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 7. Control Flow\n",
+    "\n",
+    "Loops and conditionals are needed for any non-trivial task. Please note that **whitespace matters in Python**. Everything that is indented at the same level is part of the same block. By far the most common loops in Python are for-each loops as shown in the following. While loops also exist but are rarely used."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "temp = [\"a\", \"b\", \"c\"]\n",
+    "\n",
+    "# The typical Python loop is a for-each loop, e.g.\n",
+    "for item in temp:\n",
+    "    # Everything with the same indentation is part of the loop.\n",
+    "    new_item = item + \" \" + item\n",
+    "    print(new_item)\n",
+    "    \n",
+    "print(\"No more part of the loop.\")    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Useful to know is the range() function.\n",
+    "for i in range(5):\n",
+    "    print(i)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The second crucial control flow structure are if/else conditional and they work the same as in any other language."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# If/else works as expected.\n",
+    "age = 77\n",
+    "\n",
+    "if age >= 0 and age < 10:\n",
+    "    print(\"Younger ten.\")\n",
+    "elif age >= 10:\n",
+    "    print(\"Older than ten.\")\n",
+    "else:\n",
+    "    print(\"wait what?\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# List comprehensions are a nice way to write compact loops.\n",
+    "# Make sure you understand this as it is very common in Python.\n",
+    "\n",
+    "a = list(range(10))\n",
+    "print(a)\n",
+    "b = [i for i in a if not i % 2]\n",
+    "print(b)\n",
+    "\n",
+    "# Equivalant loop for b.\n",
+    "b = []\n",
+    "for i in a:\n",
+    "    if not i % 2:\n",
+    "        b.append(i)\n",
+    "print(b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 8. Error Messages\n",
+    "\n",
+    "You will eventually run into some error messages. Learn to read them! The last line is often the one that matters - reading upwards traces the error back in time and shows what calls led to it. If stuck: just google the error message!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "def do_something(a, b): \n",
+    "    print(a + b + something_else)\n",
+    "    \n",
+    "do_something(1, 2)    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The Scientific Python Ecosystem\n",
+    "\n",
+    "The [SciPy Stack](https://www.scipy.org/stackspec.html) forms the basis for essentially all applications of scientific Python. Here we will quickly introduce the three core libraries:\n",
+    "\n",
+    "* `NumPy`\n",
+    "* `SciPy`\n",
+    "* `Matplotlib`\n",
+    "\n",
+    "The SciPy stack furthermore contains `pandas` (library for data analysis on tabular and time series data) and `sympy` (package for symbolic math), both very powerful packages, but we will omit them in this tutorial."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 9. NumPy\n",
+    "\n",
+    "Large parts of the scientific Python ecosystem use NumPy, an array computation package offering N-dimensional, typed arrays and useful functions for linear algebra, Fourier transforms, random numbers, and other basic scientific tasks."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "# Create a large array with with 1 million samples.\n",
+    "x = np.linspace(start=0, stop=100, num=int(1E6), dtype=np.float64)\n",
+    "\n",
+    "# Most operations work per-element.\n",
+    "y = x ** 2\n",
+    "\n",
+    "# Uses C and Fortran under the hood for speed.\n",
+    "print(y.sum())\n",
+    "\n",
+    "# FFT and inverse\n",
+    "x = np.random.random(100)\n",
+    "large_X = np.fft.fft(x)\n",
+    "x = np.fft.ifft(large_X)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 10. SciPy\n",
+    "\n",
+    "`SciPy`, in contrast to `NumPy` which only offers basic numerical routines, contains a lot of additional functionality needed for scientific work. Examples are solvers for basic differential equations, numeric integration and optimization, spare matrices, interpolation routines, signal processing methods, and a lot of other things."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "from scipy.interpolate import interp1d\n",
+    "\n",
+    "x = np.linspace(0, 10, num=11, endpoint=True)\n",
+    "y = np.cos(-x ** 2 / 9.0)\n",
+    "\n",
+    "# Cubic spline interpolation to new points.\n",
+    "f2 = interp1d(x, y, kind='cubic')(np.linspace(0, 10, num=101, endpoint=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 11. Matplotlib\n",
+    "\n",
+    "Plotting is done using `Matplotlib`, a package for greating high-quality static plots. It has an interface that mimics Matlab which many people are familiar with."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "plt.plot(np.sin(np.linspace(0, 2 * np.pi, 2000)), color=\"green\",\n",
+    "         label=\"Some Curve\")\n",
+    "plt.legend()\n",
+    "plt.ylim(-1.1, 1.1)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Exercises\n",
+    "\n",
+    "#### Functions, NumPy, and Matplotlib\n",
+    "\n",
+    "A. Write a function that takes a NumPy array `x` and `a`, `b`, and `c` and returns\n",
+    "\n",
+    "$$\n",
+    "f(x) = a x^2 + b x + c\n",
+    "$$\n",
+    "\n",
+    "B. Plot the result of that function with matplotlib."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### 99 Bottles of Beer\n",
+    "\n",
+    "*(stolen from http://www.ling.gu.se/~lager/python_exercises.html)*\n",
+    "\n",
+    "\n",
+    "\"99 Bottles of Beer\" is a traditional song in the United States and Canada. It is popular to sing on long trips, as it has a very repetitive format which is easy to memorize, and can take a long time to sing. The song's simple lyrics are as follows:\n",
+    "\n",
+    "```\n",
+    "99 bottles of beer on the wall, 99 bottles of beer.\n",
+    "Take one down, pass it around, 98 bottles of beer on the wall.\n",
+    "```\n",
+    "\n",
+    "The same verse is repeated, each time with one fewer bottle. The song is completed when the singer or singers reach zero.\n",
+    "\n",
+    "Your task here is write a Python program capable of generating all the verses of the song.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Ceasar Cipher\n",
+    "\n",
+    "*(stolen from http://www.ling.gu.se/~lager/python_exercises.html)*\n",
+    "\n",
+    "In cryptography, a Caesar cipher is a very simple encryption techniques in which each letter in the plain text is replaced by a letter some fixed number of positions down the alphabet. For example, with a shift of 3, A would be replaced by D, B would become E, and so on. The method is named after Julius Caesar, who used it to communicate with his generals. ROT-13 (\"rotate by 13 places\") is a widely used example of a Caesar cipher where the shift is 13. In Python, the key for ROT-13 may be represented by means of the following dictionary:\n",
+    "\n",
+    "```python\n",
+    "key = {'a':'n', 'b':'o', 'c':'p', 'd':'q', 'e':'r', 'f':'s', 'g':'t', 'h':'u', \n",
+    "       'i':'v', 'j':'w', 'k':'x', 'l':'y', 'm':'z', 'n':'a', 'o':'b', 'p':'c', \n",
+    "       'q':'d', 'r':'e', 's':'f', 't':'g', 'u':'h', 'v':'i', 'w':'j', 'x':'k',\n",
+    "       'y':'l', 'z':'m', 'A':'N', 'B':'O', 'C':'P', 'D':'Q', 'E':'R', 'F':'S', \n",
+    "       'G':'T', 'H':'U', 'I':'V', 'J':'W', 'K':'X', 'L':'Y', 'M':'Z', 'N':'A', \n",
+    "       'O':'B', 'P':'C', 'Q':'D', 'R':'E', 'S':'F', 'T':'G', 'U':'H', 'V':'I', \n",
+    "       'W':'J', 'X':'K', 'Y':'L', 'Z':'M'}\n",
+    "```\n",
+    "\n",
+    "Your task in this exercise is to implement an decoder of ROT-13. Once you're done, you will be able to read the following secret message:\n",
+    "\n",
+    "```\n",
+    "Pnrfne pvcure? V zhpu cersre Pnrfne fnynq!\n",
+    "```\n",
+    "\n",
+    "**BONUS:** Write an encoder!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "obspy",
+   "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.8.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}

01-Python_Introduction/Python_Crash_Course_solution.ipynb

@@ -0,0 +1,916 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "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%\">Scientific Python</div>\n",
+    "            <div style=\"font-size: large ; padding-top: 20px ; color: rgba(0 , 0 , 0 , 0.5)\">A super quick crash course</div>\n",
+    "        </div>\n",
+    "    </div>\n",
+    "</div>"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Seismo-Live: http://seismo-live.org\n",
+    "\n",
+    "##### Authors:\n",
+    "* Lion Krischer ([@krischer](https://github.com/krischer))\n",
+    "\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This notebook is a very quick introduction to Python and in particular its scientific ecosystem in case you have never seen it before. It furthermore grants a possibility to get to know the [IPython/Jupyter notebook](http://www.nature.com/news/interactive-notebooks-sharing-the-code-1.16261). [See here for the official documentation](http://nbviewer.jupyter.org/github/jupyter/notebook/blob/master/docs/source/examples/Notebook/Notebook%20Basics.ipynb) of the Jupyter notebook - a ton more information can be found online.\n",
+    "\n",
+    "\n",
+    "A lot of motivational writing on *Why Python?* is out there so we will not repeat it here and just condense it to a single sentence: **Python is a good and easy to learn, open-source, general purpose programming language that happens to be very good for many scientific tasks (due to its vast scientific ecosystem).**\n",
+    "\n",
+    "\n",
+    "#### Quick Reference on How to Use This Notebook\n",
+    "\n",
+    "\n",
+    "<img src=\"../images/notebook_toolbar.png\" style=\"width:70%\"></img>\n",
+    "\n",
+    "* `Shift + Enter`: Execute cell and jump to the next cell\n",
+    "* `Ctrl/Cmd + Enter`: Execute cell and don't jump to the next cell\n",
+    "\n",
+    "\n",
+    "#### Disclaimer\n",
+    "\n",
+    "The tutorials are employing Jupyter notebooks but these are only one way of using Python. Writing scripts to text files and executing them with the Python interpreter of course also works:\n",
+    "\n",
+    "```bash\n",
+    "$ python do_something.py\n",
+    "```\n",
+    "\n",
+    "Another alternative is interactive usage on the command line:\n",
+    "\n",
+    "```bash\n",
+    "$ ipython\n",
+    "```\n",
+    "\n",
+    "## Notebook Setup\n",
+    "\n",
+    "First things first: In many notebooks you will find a cell similar to the following one. **Always execute it!** They do a couple of things:\n",
+    "* Make plots appear in the browser (otherwise a window pops up)\n",
+    "* Printing things works like this: \n",
+    "\n",
+    "```python\n",
+    "print(\"Hello\")\n",
+    "```\n",
+    "\n",
+    "This essentially makes the notebooks work under Python 2 and Python 3.\n",
+    "\n",
+    "* Plots look quite a bit nicer (this is optional).\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Plots now appear in the notebook.\n",
+    "%matplotlib inline                       \n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "plt.style.use('ggplot')                            # Matplotlib style sheet - nicer plots!\n",
+    "plt.rcParams['figure.figsize'] = 12, 8             # Slightly bigger plots by default"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "---\n",
+    "\n",
+    "## Useful Links\n",
+    "\n",
+    "Here is collection of resources regarding the scientific Python ecosystem. They cover a number of different packages and topics; way more than we will manage today.\n",
+    "\n",
+    "If you have any question regarding some specific Python functionality you can consult the official [Python documenation](http://docs.python.org/).\n",
+    " \n",
+    "Furthermore a large number of Python tutorials, introductions, and books are available online. Here are some examples for those interested in learning more.\n",
+    " \n",
+    "* [Learn Python The Hard Way](http://learnpythonthehardway.org/book/)\n",
+    "* [Dive Into Python](http://www.diveintopython.net/)\n",
+    "* [The Official Python Tutorial](http://docs.python.org/2/tutorial/index.html)\n",
+    "* [Think Python Book](http://www.greenteapress.com/thinkpython/thinkpython.html)\n",
+    " \n",
+    "Some people might be used to Matlab - this helps:\n",
+    " \n",
+    "* [NumPy for Matlab Users Introdution](http://wiki.scipy.org/NumPy_for_Matlab_Users)\n",
+    "* [NumPy for Matlab Users Cheatsheet](http://mathesaurus.sourceforge.net/matlab-numpy.html)\n",
+    " \n",
+    " \n",
+    "Additionally there is an abundance of resources introducing and teaching parts of the scientific Python ecosystem.\n",
+    " \n",
+    "* [NumPy Tutorial](http://wiki.scipy.org/Tentative_NumPy_Tutorial)\n",
+    "* [Probabilistic Programming and Bayesian Methods for Hackers](http://camdavidsonpilon.github.io/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/): Great ebook introducing Bayesian methods from an understanding-first point of view with the examples done in Python.\n",
+    "* [Python Scientific Lecture Notes](http://scipy-lectures.github.io/): Introduces the basics of scientific Python with lots of examples.\n",
+    "* [Python for Signal Processing](http://python-for-signal-processing.blogspot.de/): Free blog which is the basis of a proper book written on the subject.\n",
+    "* [Another NumPy Tutorial](http://www.loria.fr/~rougier/teaching/numpy/numpy.html), [Matplotlib Tutorial](http://www.loria.fr/~rougier/teaching/matplotlib/matplotlib.html)\n",
+    " \n",
+    "You might eventually have a need to create some custom plots. The quickest way to success is usually to start from some example that is somewhat similar to what you want to achieve and just modify it. These websites are good starting points:\n",
+    " \n",
+    "* [Matplotlib Gallery](http://matplotlib.org/gallery.html)\n",
+    "* [ObsPy Gallery](http://docs.obspy.org/gallery.html)\n",
+    "* [Basemap Gallery](http://matplotlib.org/basemap/users/examples.html)\n",
+    "\n",
+    "\n",
+    "---"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Core Python Crash Course\n",
+    "\n",
+    "This course is fairly non-interactive and serves to get you up to speed with Python assuming you have practical programming experience with at least one other language. Nonetheless please change things and play around an your own - it is the only way to really learn it!\n",
+    "\n",
+    "The first part will introduce you to the core Python language. This tutorial uses Python 3 but almost all things can be transferred to Python 2. If possible choose Python 3 for your own work!\n",
+    "\n",
+    "\n",
+    "### 1. Numbers\n",
+    "\n",
+    "Python is dynamically typed and assigning something to a variable will give it that type."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Three basic types of numbers\n",
+    "a = 1             # Integers\n",
+    "b = 2.0           # Floating Point Numbers\n",
+    "c = 3.0 + 4j      # Complex Numbers, note the use of j for the complex part\n",
+    "\n",
+    "\n",
+    "# Arithmetics work as expected.\n",
+    "# Upcasting from int -> float -> complex\n",
+    "d = a + b         # (int + float = float)\n",
+    "print(d)\n",
+    "\n",
+    "e = c ** 2        # c to the second power, performs a complex multiplication\n",
+    "print(e)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import math\n",
+    "math.sqrt(7*7+24*24)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 2. Strings"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Just enclose something in single or double quotes and it will become a string. On Python 3 it defaults to unicode strings, e.g. non Latin alphabets and other symbols."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# You can use single or double quotes to create strings.\n",
+    "location = \"New York\"\n",
+    "\n",
+    "# Concatenate strings with plus.\n",
+    "where_am_i = 'I am in ' + location\n",
+    "\n",
+    "# Print things with the print() function.\n",
+    "print(location, 1, 2)\n",
+    "print(where_am_i)\n",
+    "\n",
+    "# Strings have a lot of attached methods for common manipulations.\n",
+    "print(location.lower())\n",
+    "\n",
+    "# Access single items with square bracket. Negative indices are from the back.\n",
+    "print(location[0], location[-1])\n",
+    "\n",
+    "# Strings can also be sliced.\n",
+    "print(location[4:])"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Exercise\n",
+    "\n",
+    "Save your name in all lower-case letters to a variable, and print a capitalized version of it. Protip: [Google for \"How to capitalize a string in python\"](http://www.google.com/search?q=how+to+capitalize+a+string+in+python). This works for almost any programming problem - someone will have had the same issue before!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 3. Lists"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Python has two main collection types: List and dictionaries. The former is just an ordered collection of objects and is introduced here."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# List use square brackets and are simple ordered collections of things.\n",
+    "everything = [a, b, c, 1, 2, 3, \"hello\"]\n",
+    "\n",
+    "# Access elements with the same slicing/indexing notation as strings.\n",
+    "# Note that Python indices are zero based!\n",
+    "print(everything[0])\n",
+    "print(everything[:3])\n",
+    "print(everything[2:-2])\n",
+    "\n",
+    "# Negative indices are counted from the back of the list.\n",
+    "print(everything[-3:])\n",
+    "\n",
+    "# Append things with the append method.\n",
+    "everything.append(\"you\")\n",
+    "print(everything)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 4. Dictionaries\n",
+    "\n",
+    "The other main collection type in Python are dictionaries. They are similiar to associative arrays or (hash) maps in other languages. Each entry is a key-value pair."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Dictionaries have named fields and no inherent order. As is\n",
+    "# the case with lists, they can contain anything.\n",
+    "information = {\n",
+    "    \"name\": \"Hans\",\n",
+    "    \"surname\": \"Mustermann\",\n",
+    "    \"age\": 78,\n",
+    "    \"kids\": [1, 2, 3]\n",
+    "}\n",
+    "\n",
+    "# Acccess items by using the key in square brackets.\n",
+    "print(information[\"kids\"])\n",
+    "\n",
+    "# Add new things by just assigning to a key.\n",
+    "print(information)\n",
+    "information[\"music\"] = \"jazz\"\n",
+    "print(information)\n",
+    "\n",
+    "# Delete things by using the del operator\n",
+    "del information[\"age\"]\n",
+    "print(information)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 5. Functions\n",
+    "\n",
+    "The key to conquer a big problem is to divide it into many smaller ones and tackle them one by one. This is usually achieved by using functions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Functions are defined using the def keyword.\n",
+    "def do_stuff(a, b):\n",
+    "    return a * b\n",
+    "\n",
+    "# And called with the arguments in round brackets.\n",
+    "print(do_stuff(2, 3))\n",
+    "\n",
+    "# Python function also can have optional arguments.\n",
+    "def do_more_stuff(a, b, power=1):\n",
+    "    return (a * b) ** power\n",
+    "\n",
+    "print(do_more_stuff(2, 3))\n",
+    "print(do_more_stuff(2, 3, power=3))\n",
+    "\n",
+    "# For more complex function it is oftentimes a good idea to \n",
+    "#explicitly name the arguments. This is easier to read and less error-prone.\n",
+    "print(do_more_stuff(a=2, b=3, power=3))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 6. Imports\n",
+    "\n",
+    "To use functions and objects not part of the default namespace, you have import them. You will have to do this a lot so it is necessary to learn how to do it."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Import anything, and use it with the dot accessor.\n",
+    "import math\n",
+    "\n",
+    "a = math.cos(4 * math.pi)\n",
+    "\n",
+    "# You can also selectively import things.\n",
+    "from math import pi\n",
+    "\n",
+    "b = 3 * pi\n",
+    "\n",
+    "# And even rename them if you don't like their name.\n",
+    "from math import cos as cosine\n",
+    "c = cosine(b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "How to know what is available?\n",
+    "\n",
+    "1. Read the [documentation](https://docs.python.org/3/library/math.html)\n",
+    "2. Interactively query the module"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "print(dir(math))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Typing the dot and the TAB will kick off tab-completion."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "math."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In the IPython framework you can also use a question mark to view the documentation of modules and functions."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": true
+    }
+   },
+   "outputs": [],
+   "source": [
+    "math.cos?"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 7. Control Flow\n",
+    "\n",
+    "Loops and conditionals are needed for any non-trivial task. Please note that **whitespace matters in Python**. Everything that is indented at the same level is part of the same block. By far the most common loops in Python are for-each loops as shown in the following. While loops also exist but are rarely used."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "temp = [\"a\", \"b\", \"c\"]\n",
+    "\n",
+    "# The typical Python loop is a for-each loop, e.g.\n",
+    "for item in temp:\n",
+    "    # Everything with the same indentation is part of the loop.\n",
+    "    new_item = item + \" \" + item\n",
+    "    print(new_item)\n",
+    "    \n",
+    "print(\"No more part of the loop.\")    "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# Useful to know is the range() function.\n",
+    "for i in range(5):\n",
+    "    print(i)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The second crucial control flow structure are if/else conditional and they work the same as in any other language."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# If/else works as expected.\n",
+    "age = 77\n",
+    "\n",
+    "if age >= 0 and age < 10:\n",
+    "    print(\"Younger ten.\")\n",
+    "elif age >= 10:\n",
+    "    print(\"Older than ten.\")\n",
+    "else:\n",
+    "    print(\"wait what?\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# List comprehensions are a nice way to write compact loops.\n",
+    "# Make sure you understand this as it is very common in Python.\n",
+    "\n",
+    "a = list(range(10))\n",
+    "print(a)\n",
+    "b = [i for i in a if not i % 2]\n",
+    "print(b)\n",
+    "\n",
+    "# Equivalant loop for b.\n",
+    "b = []\n",
+    "for i in a:\n",
+    "    if not i % 2:\n",
+    "        b.append(i)\n",
+    "print(b)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 8. Error Messages\n",
+    "\n",
+    "You will eventually run into some error messages. Learn to read them! The last line is often the one that matters - reading upwards traces the error back in time and shows what calls led to it. If stuck: just google the error message!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def do_something(a, b): \n",
+    "    print(a + b + something_else)\n",
+    "    \n",
+    "do_something(1, 2)    "
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The Scientific Python Ecosystem\n",
+    "\n",
+    "The [SciPy Stack](https://www.scipy.org/stackspec.html) forms the basis for essentially all applications of scientific Python. Here we will quickly introduce the three core libraries:\n",
+    "\n",
+    "* `NumPy`\n",
+    "* `SciPy`\n",
+    "* `Matplotlib`\n",
+    "\n",
+    "The SciPy stack furthermore contains `pandas` (library for data analysis on tabular and time series data) and `sympy` (package for symbolic math), both very powerful packages, but we will omit them in this tutorial."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 9. NumPy\n",
+    "\n",
+    "Large parts of the scientific Python ecosystem use NumPy, an array computation package offering N-dimensional, typed arrays and useful functions for linear algebra, Fourier transforms, random numbers, and other basic scientific tasks."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "\n",
+    "# Create a large array with with 1 million samples.\n",
+    "x = np.linspace(start=0, stop=100, num=int(1E6), dtype=np.float64)\n",
+    "\n",
+    "# Most operations work per-element.\n",
+    "y = x ** 2\n",
+    "\n",
+    "# Uses C and Fortran under the hood for speed.\n",
+    "print(y.sum())\n",
+    "\n",
+    "# FFT and inverse\n",
+    "x = np.random.random(100)\n",
+    "large_X = np.fft.fft(x)\n",
+    "x = np.fft.ifft(large_X)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 10. SciPy\n",
+    "\n",
+    "`SciPy`, in contrast to `NumPy` which only offers basic numerical routines, contains a lot of additional functionality needed for scientific work. Examples are solvers for basic differential equations, numeric integration and optimization, spare matrices, interpolation routines, signal processing methods, and a lot of other things."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from scipy.interpolate import interp1d\n",
+    "\n",
+    "x = np.linspace(0, 10, num=11, endpoint=True)\n",
+    "y = np.cos(-x ** 2 / 9.0)\n",
+    "\n",
+    "# Cubic spline interpolation to new points.\n",
+    "f2 = interp1d(x, y, kind='cubic')(np.linspace(0, 10, num=101, endpoint=True))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### 11. Matplotlib\n",
+    "\n",
+    "Plotting is done using `Matplotlib`, a package for greating high-quality static plots. It has an interface that mimics Matlab which many people are familiar with."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "plt.plot(np.sin(np.linspace(0, 2 * np.pi, 2000)), color=\"green\",\n",
+    "         label=\"Some Curve\")\n",
+    "plt.legend()\n",
+    "plt.ylim(-1.1, 1.1)\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Exercises\n",
+    "\n",
+    "#### Functions, NumPy, and Matplotlib\n",
+    "\n",
+    "A. Write a function that takes a NumPy array `x` and `a`, `b`, and `c` and returns\n",
+    "\n",
+    "$$\n",
+    "f(x) = a x^2 + b x + c\n",
+    "$$\n",
+    "\n",
+    "B. Plot the result of that function with matplotlib."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "def my_calculation(x, a, b, c):\n",
+    "    f = a * x**2 + b * x + c\n",
+    "    return f\n",
+    "\n",
+    "data = np.array([1, 2, 3])\n",
+    "a = 1\n",
+    "b = 0\n",
+    "c = 0\n",
+    "my_calculation(data, a, b, c)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# define x values and calculate y values\n",
+    "x = np.linspace(-10, 50, 30)\n",
+    "y = my_calculation(x, a, b, c)\n",
+    "\n",
+    "# create a plot\n",
+    "plt.legend()\n",
+    "plt.plot(x, y, color=\"green\", label=\"My Calculation\")\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### 99 Bottles of Beer\n",
+    "\n",
+    "*(stolen from http://www.ling.gu.se/~lager/python_exercises.html)*\n",
+    "\n",
+    "\n",
+    "\"99 Bottles of Beer\" is a traditional song in the United States and Canada. It is popular to sing on long trips, as it has a very repetitive format which is easy to memorize, and can take a long time to sing. The song's simple lyrics are as follows:\n",
+    "\n",
+    "```\n",
+    "99 bottles of beer on the wall, 99 bottles of beer.\n",
+    "Take one down, pass it around, 98 bottles of beer on the wall.\n",
+    "```\n",
+    "\n",
+    "The same verse is repeated, each time with one fewer bottle. The song is completed when the singer or singers reach zero.\n",
+    "\n",
+    "Your task here is write a Python program capable of generating all the verses of the song.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "for no_bottles in range(99, 0, -1):\n",
+    "    print(f\"{no_bottles} bottles of beer on the wall, {no_bottles} bottles of beer.\")\n",
+    "    no_bottles -= 1\n",
+    "    print(f\"Take one down, pass it around, {no_bottles} bottles of beer on the wall.\")\n",
+    "    print()\n",
+    "print(\"All bottles gone.\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Ceasar Cipher\n",
+    "\n",
+    "*(stolen from http://www.ling.gu.se/~lager/python_exercises.html)*\n",
+    "\n",
+    "In cryptography, a Caesar cipher is a very simple encryption techniques in which each letter in the plain text is replaced by a letter some fixed number of positions down the alphabet. For example, with a shift of 3, A would be replaced by D, B would become E, and so on. The method is named after Julius Caesar, who used it to communicate with his generals. ROT-13 (\"rotate by 13 places\") is a widely used example of a Caesar cipher where the shift is 13. In Python, the key for ROT-13 may be represented by means of the following dictionary:\n",
+    "\n",
+    "```python\n",
+    "key = {'a':'n', 'b':'o', 'c':'p', 'd':'q', 'e':'r', 'f':'s', 'g':'t', 'h':'u', \n",
+    "       'i':'v', 'j':'w', 'k':'x', 'l':'y', 'm':'z', 'n':'a', 'o':'b', 'p':'c', \n",
+    "       'q':'d', 'r':'e', 's':'f', 't':'g', 'u':'h', 'v':'i', 'w':'j', 'x':'k',\n",
+    "       'y':'l', 'z':'m', 'A':'N', 'B':'O', 'C':'P', 'D':'Q', 'E':'R', 'F':'S', \n",
+    "       'G':'T', 'H':'U', 'I':'V', 'J':'W', 'K':'X', 'L':'Y', 'M':'Z', 'N':'A', \n",
+    "       'O':'B', 'P':'C', 'Q':'D', 'R':'E', 'S':'F', 'T':'G', 'U':'H', 'V':'I', \n",
+    "       'W':'J', 'X':'K', 'Y':'L', 'Z':'M'}\n",
+    "```\n",
+    "\n",
+    "Your task in this exercise is to implement an decoder of ROT-13. Once you're done, you will be able to read the following secret message:\n",
+    "\n",
+    "```\n",
+    "Pnrfne pvcure? V zhpu cersre Pnrfne fnynq!\n",
+    "```\n",
+    "\n",
+    "**BONUS:** Write an encoder!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# define decoding key\n",
+    "decode_key = {'a':'n', 'b':'o', 'c':'p', 'd':'q', 'e':'r', 'f':'s', 'g':'t', 'h':'u', \n",
+    "       'i':'v', 'j':'w', 'k':'x', 'l':'y', 'm':'z', 'n':'a', 'o':'b', 'p':'c', \n",
+    "       'q':'d', 'r':'e', 's':'f', 't':'g', 'u':'h', 'v':'i', 'w':'j', 'x':'k',\n",
+    "       'y':'l', 'z':'m', 'A':'N', 'B':'O', 'C':'P', 'D':'Q', 'E':'R', 'F':'S', \n",
+    "       'G':'T', 'H':'U', 'I':'V', 'J':'W', 'K':'X', 'L':'Y', 'M':'Z', 'N':'A', \n",
+    "       'O':'B', 'P':'C', 'Q':'D', 'R':'E', 'S':'F', 'T':'G', 'U':'H', 'V':'I', \n",
+    "       'W':'J', 'X':'K', 'Y':'L', 'Z':'M'}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "# decode a secret\n",
+    "secret = \"Pnrfne pvcure? V zhpu cersre Pnrfne fnynq!\"\n",
+    "text = \"\"\n",
+    "\n",
+    "for letter in secret:\n",
+    "    try:\n",
+    "        text += decode_key[letter]\n",
+    "    except KeyError:\n",
+    "        text += letter\n",
+    "print(text)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# revers decoding_key to get encoding_key\n",
+    "encode_key = {}\n",
+    "for k in decode_key:\n",
+    "    encode_key = {decode_key[k]: k}"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# encode a text\n",
+    "secret = \"\"\n",
+    "text = \"Python can be fun!\"\n",
+    "for letter in text:\n",
+    "    try:\n",
+    "        secret += decode_key[letter]\n",
+    "    except KeyError:\n",
+    "        secret += letter\n",
+    "print(secret)"
+   ]
+  }
+ ],
+ "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.9.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

02-Fourier_Transformation_and_Series/Fourier_Transformation_and_Series.ipynb

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

LICENSE

@@ -1,158 +1,437 @@
-Creative Commons Attribution-NonCommercial 4.0 International
+Attribution-NonCommercial-ShareAlike 4.0 International
 
- Creative Commons Corporation (“Creative Commons”) is not a law firm and does not provide legal services or legal advice. Distribution of Creative Commons public licenses does not create a lawyer-client or other relationship. Creative Commons makes its licenses and related information available on an “as-is” basis. Creative Commons gives no warranties regarding its licenses, any material licensed under their terms and conditions, or any related information. Creative Commons disclaims all liability for damages resulting from their use to the fullest extent possible.
+=======================================================================
+
+Creative Commons Corporation ("Creative Commons") is not a law firm and
+does not provide legal services or legal advice. Distribution of
+Creative Commons public licenses does not create a lawyer-client or
+other relationship. Creative Commons makes its licenses and related
+information available on an "as-is" basis. Creative Commons gives no
+warranties regarding its licenses, any material licensed under their
+terms and conditions, or any related information. Creative Commons
+disclaims all liability for damages resulting from their use to the
+fullest extent possible.
 
 Using Creative Commons Public Licenses
 
-Creative Commons public licenses provide a standard set of terms and conditions that creators and other rights holders may use to share original works of authorship and other material subject to copyright and certain other rights specified in the public license below. The following considerations are for informational purposes only, are not exhaustive, and do not form part of our licenses.
+Creative Commons public licenses provide a standard set of terms and
+conditions that creators and other rights holders may use to share
+original works of authorship and other material subject to copyright
+and certain other rights specified in the public license below. The
+following considerations are for informational purposes only, are not
+exhaustive, and do not form part of our licenses.
+
+     Considerations for licensors: Our public licenses are
+     intended for use by those authorized to give the public
+     permission to use material in ways otherwise restricted by
+     copyright and certain other rights. Our licenses are
+     irrevocable. Licensors should read and understand the terms
+     and conditions of the license they choose before applying it.
+     Licensors should also secure all rights necessary before
+     applying our licenses so that the public can reuse the
+     material as expected. Licensors should clearly mark any
+     material not subject to the license. This includes other CC-
+     licensed material, or material used under an exception or
+     limitation to copyright. More considerations for licensors:
+	wiki.creativecommons.org/Considerations_for_licensors
+
+     Considerations for the public: By using one of our public
+     licenses, a licensor grants the public permission to use the
+     licensed material under specified terms and conditions. If
+     the licensor's permission is not necessary for any reason--for
+     example, because of any applicable exception or limitation to
+     copyright--then that use is not regulated by the license. Our
+     licenses grant only permissions under copyright and certain
+     other rights that a licensor has authority to grant. Use of
+     the licensed material may still be restricted for other
+     reasons, including because others have copyright or other
+     rights in the material. A licensor may make special requests,
+     such as asking that all changes be marked or described.
+     Although not required by our licenses, you are encouraged to
+     respect those requests where reasonable. More_considerations
+     for the public:
+	wiki.creativecommons.org/Considerations_for_licensees
+
+=======================================================================
+
+Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
+Public License
+
+By exercising the Licensed Rights (defined below), You accept and agree
+to be bound by the terms and conditions of this Creative Commons
+Attribution-NonCommercial-ShareAlike 4.0 International Public License
+("Public License"). To the extent this Public License may be
+interpreted as a contract, You are granted the Licensed Rights in
+consideration of Your acceptance of these terms and conditions, and the
+Licensor grants You such rights in consideration of benefits the
+Licensor receives from making the Licensed Material available under
+these terms and conditions.
+
+
+Section 1 -- Definitions.
+
+  a. Adapted Material means material subject to Copyright and Similar
+     Rights that is derived from or based upon the Licensed Material
+     and in which the Licensed Material is translated, altered,
+     arranged, transformed, or otherwise modified in a manner requiring
+     permission under the Copyright and Similar Rights held by the
+     Licensor. For purposes of this Public License, where the Licensed
+     Material is a musical work, performance, or sound recording,
+     Adapted Material is always produced where the Licensed Material is
+     synched in timed relation with a moving image.
+
+  b. Adapter's License means the license You apply to Your Copyright
+     and Similar Rights in Your contributions to Adapted Material in
+     accordance with the terms and conditions of this Public License.
+
+  c. BY-NC-SA Compatible License means a license listed at
+     creativecommons.org/compatiblelicenses, approved by Creative
+     Commons as essentially the equivalent of this Public License.
+
+  d. Copyright and Similar Rights means copyright and/or similar rights
+     closely related to copyright including, without limitation,
+     performance, broadcast, sound recording, and Sui Generis Database
+     Rights, without regard to how the rights are labeled or
+     categorized. For purposes of this Public License, the rights
+     specified in Section 2(b)(1)-(2) are not Copyright and Similar
+     Rights.
+
+  e. Effective Technological Measures means those measures that, in the
+     absence of proper authority, may not be circumvented under laws
+     fulfilling obligations under Article 11 of the WIPO Copyright
+     Treaty adopted on December 20, 1996, and/or similar international
+     agreements.
+
+  f. Exceptions and Limitations means fair use, fair dealing, and/or
+     any other exception or limitation to Copyright and Similar Rights
+     that applies to Your use of the Licensed Material.
+
+  g. License Elements means the license attributes listed in the name
+     of a Creative Commons Public License. The License Elements of this
+     Public License are Attribution, NonCommercial, and ShareAlike.
+
+  h. Licensed Material means the artistic or literary work, database,
+     or other material to which the Licensor applied this Public
+     License.
+
+  i. Licensed Rights means the rights granted to You subject to the
+     terms and conditions of this Public License, which are limited to
+     all Copyright and Similar Rights that apply to Your use of the
+     Licensed Material and that the Licensor has authority to license.
+
+  j. Licensor means the individual(s) or entity(ies) granting rights
+     under this Public License.
+
+  k. NonCommercial means not primarily intended for or directed towards
+     commercial advantage or monetary compensation. For purposes of
+     this Public License, the exchange of the Licensed Material for
+     other material subject to Copyright and Similar Rights by digital
+     file-sharing or similar means is NonCommercial provided there is
+     no payment of monetary compensation in connection with the
+     exchange.
+
+  l. Share means to provide material to the public by any means or
+     process that requires permission under the Licensed Rights, such
+     as reproduction, public display, public performance, distribution,
+     dissemination, communication, or importation, and to make material
+     available to the public including in ways that members of the
+     public may access the material from a place and at a time
+     individually chosen by them.
+
+  m. Sui Generis Database Rights means rights other than copyright
+     resulting from Directive 96/9/EC of the European Parliament and of
+     the Council of 11 March 1996 on the legal protection of databases,
+     as amended and/or succeeded, as well as other essentially
+     equivalent rights anywhere in the world.
+
+  n. You means the individual or entity exercising the Licensed Rights
+     under this Public License. Your has a corresponding meaning.
+
+
+Section 2 -- Scope.
+
+  a. License grant.
+
+       1. Subject to the terms and conditions of this Public License,
+          the Licensor hereby grants You a worldwide, royalty-free,
+          non-sublicensable, non-exclusive, irrevocable license to
+          exercise the Licensed Rights in the Licensed Material to:
+
+            a. reproduce and Share the Licensed Material, in whole or
+               in part, for NonCommercial purposes only; and
+
+            b. produce, reproduce, and Share Adapted Material for
+               NonCommercial purposes only.
+
+       2. Exceptions and Limitations. For the avoidance of doubt, where
+          Exceptions and Limitations apply to Your use, this Public
+          License does not apply, and You do not need to comply with
+          its terms and conditions.
+
+       3. Term. The term of this Public License is specified in Section
+          6(a).
+
+       4. Media and formats; technical modifications allowed. The
+          Licensor authorizes You to exercise the Licensed Rights in
+          all media and formats whether now known or hereafter created,
+          and to make technical modifications necessary to do so. The
+          Licensor waives and/or agrees not to assert any right or
+          authority to forbid You from making technical modifications
+          necessary to exercise the Licensed Rights, including
+          technical modifications necessary to circumvent Effective
+          Technological Measures. For purposes of this Public License,
+          simply making modifications authorized by this Section 2(a)
+          (4) never produces Adapted Material.
+
+       5. Downstream recipients.
+
+            a. Offer from the Licensor -- Licensed Material. Every
+               recipient of the Licensed Material automatically
+               receives an offer from the Licensor to exercise the
+               Licensed Rights under the terms and conditions of this
+               Public License.
+
+            b. Additional offer from the Licensor -- Adapted Material.
+               Every recipient of Adapted Material from You
+               automatically receives an offer from the Licensor to
+               exercise the Licensed Rights in the Adapted Material
+               under the conditions of the Adapter's License You apply.
+
+            c. No downstream restrictions. You may not offer or impose
+               any additional or different terms or conditions on, or
+               apply any Effective Technological Measures to, the
+               Licensed Material if doing so restricts exercise of the
+               Licensed Rights by any recipient of the Licensed
+               Material.
+
+       6. No endorsement. Nothing in this Public License constitutes or
+          may be construed as permission to assert or imply that You
+          are, or that Your use of the Licensed Material is, connected
+          with, or sponsored, endorsed, or granted official status by,
+          the Licensor or others designated to receive attribution as
+          provided in Section 3(a)(1)(A)(i).
+
+  b. Other rights.
+
+       1. Moral rights, such as the right of integrity, are not
+          licensed under this Public License, nor are publicity,
+          privacy, and/or other similar personality rights; however, to
+          the extent possible, the Licensor waives and/or agrees not to
+          assert any such rights held by the Licensor to the limited
+          extent necessary to allow You to exercise the Licensed
+          Rights, but not otherwise.
+
+       2. Patent and trademark rights are not licensed under this
+          Public License.
+
+       3. To the extent possible, the Licensor waives any right to
+          collect royalties from You for the exercise of the Licensed
+          Rights, whether directly or through a collecting society
+          under any voluntary or waivable statutory or compulsory
+          licensing scheme. In all other cases the Licensor expressly
+          reserves any right to collect such royalties, including when
+          the Licensed Material is used other than for NonCommercial
+          purposes.
+
+
+Section 3 -- License Conditions.
+
+Your exercise of the Licensed Rights is expressly made subject to the
+following conditions.
+
+  a. Attribution.
+
+       1. If You Share the Licensed Material (including in modified
+          form), You must:
+
+            a. retain the following if it is supplied by the Licensor
+               with the Licensed Material:
+
+                 i. identification of the creator(s) of the Licensed
+                    Material and any others designated to receive
+                    attribution, in any reasonable manner requested by
+                    the Licensor (including by pseudonym if
+                    designated);
+
+                ii. a copyright notice;
+
+               iii. a notice that refers to this Public License;
+
+                iv. a notice that refers to the disclaimer of
+                    warranties;
+
+                 v. a URI or hyperlink to the Licensed Material to the
+                    extent reasonably practicable;
+
+            b. indicate if You modified the Licensed Material and
+               retain an indication of any previous modifications; and
+
+            c. indicate the Licensed Material is licensed under this
+               Public License, and include the text of, or the URI or
+               hyperlink to, this Public License.
+
+       2. You may satisfy the conditions in Section 3(a)(1) in any
+          reasonable manner based on the medium, means, and context in
+          which You Share the Licensed Material. For example, it may be
+          reasonable to satisfy the conditions by providing a URI or
+          hyperlink to a resource that includes the required
+          information.
+       3. If requested by the Licensor, You must remove any of the
+          information required by Section 3(a)(1)(A) to the extent
+          reasonably practicable.
+
+  b. ShareAlike.
 
-Considerations for licensors: Our public licenses are intended for use by those authorized to give the public permission to use material in ways otherwise restricted by copyright and certain other rights. Our licenses are irrevocable. Licensors should read and understand the terms and conditions of the license they choose before applying it. Licensors should also secure all rights necessary before applying our licenses so that the public can reuse the material as expected. Licensors should clearly mark any material not subject to the license. This includes other CC-licensed material, or material used under an exception or limitation to copyright. More considerations for licensors.
+     In addition to the conditions in Section 3(a), if You Share
+     Adapted Material You produce, the following conditions also apply.
 
-Considerations for the public: By using one of our public licenses, a licensor grants the public permission to use the licensed material under specified terms and conditions. If the licensor’s permission is not necessary for any reason–for example, because of any applicable exception or limitation to copyright–then that use is not regulated by the license. Our licenses grant only permissions under copyright and certain other rights that a licensor has authority to grant. Use of the licensed material may still be restricted for other reasons, including because others have copyright or other rights in the material. A licensor may make special requests, such as asking that all changes be marked or described. Although not required by our licenses, you are encouraged to respect those requests where reasonable. More considerations for the public.
+       1. The Adapter's License You apply must be a Creative Commons
+          license with the same License Elements, this version or
+          later, or a BY-NC-SA Compatible License.
 
-Creative Commons Attribution-NonCommercial 4.0 International Public License
+       2. You must include the text of, or the URI or hyperlink to, the
+          Adapter's License You apply. You may satisfy this condition
+          in any reasonable manner based on the medium, means, and
+          context in which You Share Adapted Material.
 
-By exercising the Licensed Rights (defined below), You accept and agree to be bound by the terms and conditions of this Creative Commons Attribution-NonCommercial 4.0 International Public License ("Public License"). To the extent this Public License may be interpreted as a contract, You are granted the Licensed Rights in consideration of Your acceptance of these terms and conditions, and the Licensor grants You such rights in consideration of benefits the Licensor receives from making the Licensed Material available under these terms and conditions.
+       3. You may not offer or impose any additional or different terms
+          or conditions on, or apply any Effective Technological
+          Measures to, Adapted Material that restrict exercise of the
+          rights granted under the Adapter's License You apply.
 
-Section 1 – Definitions.
 
-     a.	Adapted Material means material subject to Copyright and Similar Rights that is derived from or based upon the Licensed Material and in which the Licensed Material is translated, altered, arranged, transformed, or otherwise modified in a manner requiring permission under the Copyright and Similar Rights held by the Licensor. For purposes of this Public License, where the Licensed Material is a musical work, performance, or sound recording, Adapted Material is always produced where the Licensed Material is synched in timed relation with a moving image.
+Section 4 -- Sui Generis Database Rights.
 
-     b.	Adapter's License means the license You apply to Your Copyright and Similar Rights in Your contributions to Adapted Material in accordance with the terms and conditions of this Public License.
+Where the Licensed Rights include Sui Generis Database Rights that
+apply to Your use of the Licensed Material:
 
-     c.	Copyright and Similar Rights means copyright and/or similar rights closely related to copyright including, without limitation, performance, broadcast, sound recording, and Sui Generis Database Rights, without regard to how the rights are labeled or categorized. For purposes of this Public License, the rights specified in Section 2(b)(1)-(2) are not Copyright and Similar Rights.
+  a. for the avoidance of doubt, Section 2(a)(1) grants You the right
+     to extract, reuse, reproduce, and Share all or a substantial
+     portion of the contents of the database for NonCommercial purposes
+     only;
 
-     d.	Effective Technological Measures means those measures that, in the absence of proper authority, may not be circumvented under laws fulfilling obligations under Article 11 of the WIPO Copyright Treaty adopted on December 20, 1996, and/or similar international agreements.
+  b. if You include all or a substantial portion of the database
+     contents in a database in which You have Sui Generis Database
+     Rights, then the database in which You have Sui Generis Database
+     Rights (but not its individual contents) is Adapted Material,
+     including for purposes of Section 3(b); and
 
-     e.	Exceptions and Limitations means fair use, fair dealing, and/or any other exception or limitation to Copyright and Similar Rights that applies to Your use of the Licensed Material.
+  c. You must comply with the conditions in Section 3(a) if You Share
+     all or a substantial portion of the contents of the database.
 
-     f.	Licensed Material means the artistic or literary work, database, or other material to which the Licensor applied this Public License.
+For the avoidance of doubt, this Section 4 supplements and does not
+replace Your obligations under this Public License where the Licensed
+Rights include other Copyright and Similar Rights.
 
-     g.	Licensed Rights means the rights granted to You subject to the terms and conditions of this Public License, which are limited to all Copyright and Similar Rights that apply to Your use of the Licensed Material and that the Licensor has authority to license.
 
-     h.	Licensor means the individual(s) or entity(ies) granting rights under this Public License.
+Section 5 -- Disclaimer of Warranties and Limitation of Liability.
 
-     i.	NonCommercial means not primarily intended for or directed towards commercial advantage or monetary compensation. For purposes of this Public License, the exchange of the Licensed Material for other material subject to Copyright and Similar Rights by digital file-sharing or similar means is NonCommercial provided there is no payment of monetary compensation in connection with the exchange.
+  a. UNLESS OTHERWISE SEPARATELY UNDERTAKEN BY THE LICENSOR, TO THE
+     EXTENT POSSIBLE, THE LICENSOR OFFERS THE LICENSED MATERIAL AS-IS
+     AND AS-AVAILABLE, AND MAKES NO REPRESENTATIONS OR WARRANTIES OF
+     ANY KIND CONCERNING THE LICENSED MATERIAL, WHETHER EXPRESS,
+     IMPLIED, STATUTORY, OR OTHER. THIS INCLUDES, WITHOUT LIMITATION,
+     WARRANTIES OF TITLE, MERCHANTABILITY, FITNESS FOR A PARTICULAR
+     PURPOSE, NON-INFRINGEMENT, ABSENCE OF LATENT OR OTHER DEFECTS,
+     ACCURACY, OR THE PRESENCE OR ABSENCE OF ERRORS, WHETHER OR NOT
+     KNOWN OR DISCOVERABLE. WHERE DISCLAIMERS OF WARRANTIES ARE NOT
+     ALLOWED IN FULL OR IN PART, THIS DISCLAIMER MAY NOT APPLY TO YOU.
 
-     j.	Share means to provide material to the public by any means or process that requires permission under the Licensed Rights, such as reproduction, public display, public performance, distribution, dissemination, communication, or importation, and to make material available to the public including in ways that members of the public may access the material from a place and at a time individually chosen by them.
+  b. TO THE EXTENT POSSIBLE, IN NO EVENT WILL THE LICENSOR BE LIABLE
+     TO YOU ON ANY LEGAL THEORY (INCLUDING, WITHOUT LIMITATION,
+     NEGLIGENCE) OR OTHERWISE FOR ANY DIRECT, SPECIAL, INDIRECT,
+     INCIDENTAL, CONSEQUENTIAL, PUNITIVE, EXEMPLARY, OR OTHER LOSSES,
+     COSTS, EXPENSES, OR DAMAGES ARISING OUT OF THIS PUBLIC LICENSE OR
+     USE OF THE LICENSED MATERIAL, EVEN IF THE LICENSOR HAS BEEN
+     ADVISED OF THE POSSIBILITY OF SUCH LOSSES, COSTS, EXPENSES, OR
+     DAMAGES. WHERE A LIMITATION OF LIABILITY IS NOT ALLOWED IN FULL OR
+     IN PART, THIS LIMITATION MAY NOT APPLY TO YOU.
 
-     k.	Sui Generis Database Rights means rights other than copyright resulting from Directive 96/9/EC of the European Parliament and of the Council of 11 March 1996 on the legal protection of databases, as amended and/or succeeded, as well as other essentially equivalent rights anywhere in the world.
+  c. The disclaimer of warranties and limitation of liability provided
+     above shall be interpreted in a manner that, to the extent
+     possible, most closely approximates an absolute disclaimer and
+     waiver of all liability.
 
-     l.	You means the individual or entity exercising the Licensed Rights under this Public License. Your has a corresponding meaning.
 
-Section 2 – Scope.
+Section 6 -- Term and Termination.
 
-     a.	License grant.
+  a. This Public License applies for the term of the Copyright and
+     Similar Rights licensed here. However, if You fail to comply with
+     this Public License, then Your rights under this Public License
+     terminate automatically.
+
+  b. Where Your right to use the Licensed Material has terminated under
+     Section 6(a), it reinstates:
+
+       1. automatically as of the date the violation is cured, provided
+          it is cured within 30 days of Your discovery of the
+          violation; or
+
+       2. upon express reinstatement by the Licensor.
+
+     For the avoidance of doubt, this Section 6(b) does not affect any
+     right the Licensor may have to seek remedies for Your violations
+     of this Public License.
+
+  c. For the avoidance of doubt, the Licensor may also offer the
+     Licensed Material under separate terms or conditions or stop
+     distributing the Licensed Material at any time; however, doing so
+     will not terminate this Public License.
 
-          1. Subject to the terms and conditions of this Public License, the Licensor hereby grants You a worldwide, royalty-free, non-sublicensable, non-exclusive, irrevocable license to exercise the Licensed Rights in the Licensed Material to:
+  d. Sections 1, 5, 6, 7, and 8 survive termination of this Public
+     License.
 
-               A. reproduce and Share the Licensed Material, in whole or in part, for NonCommercial purposes only; and
 
-               B. produce, reproduce, and Share Adapted Material for NonCommercial purposes only.
+Section 7 -- Other Terms and Conditions.
+
+  a. The Licensor shall not be bound by any additional or different
+     terms or conditions communicated by You unless expressly agreed.
+
+  b. Any arrangements, understandings, or agreements regarding the
+     Licensed Material not stated herein are separate from and
+     independent of the terms and conditions of this Public License.
+
+
+Section 8 -- Interpretation.
 
-          2. Exceptions and Limitations. For the avoidance of doubt, where Exceptions and Limitations apply to Your use, this Public License does not apply, and You do not need to comply with its terms and conditions.
+  a. For the avoidance of doubt, this Public License does not, and
+     shall not be interpreted to, reduce, limit, restrict, or impose
+     conditions on any use of the Licensed Material that could lawfully
+     be made without permission under this Public License.
 
-          3. Term. The term of this Public License is specified in Section 6(a).
-
-          4. Media and formats; technical modifications allowed. The Licensor authorizes You to exercise the Licensed Rights in all media and formats whether now known or hereafter created, and to make technical modifications necessary to do so. The Licensor waives and/or agrees not to assert any right or authority to forbid You from making technical modifications necessary to exercise the Licensed Rights, including technical modifications necessary to circumvent Effective Technological Measures. For purposes of this Public License, simply making modifications authorized by this Section 2(a)(4) never produces Adapted Material.
-
-          5. Downstream recipients.
-
-               A. Offer from the Licensor – Licensed Material. Every recipient of the Licensed Material automatically receives an offer from the Licensor to exercise the Licensed Rights under the terms and conditions of this Public License.
-
-               B. No downstream restrictions. You may not offer or impose any additional or different terms or conditions on, or apply any Effective Technological Measures to, the Licensed Material if doing so restricts exercise of the Licensed Rights by any recipient of the Licensed Material.
-
-          6. No endorsement. Nothing in this Public License constitutes or may be construed as permission to assert or imply that You are, or that Your use of the Licensed Material is, connected with, or sponsored, endorsed, or granted official status by, the Licensor or others designated to receive attribution as provided in Section 3(a)(1)(A)(i).
-
-     b.	Other rights.
-
-          1. Moral rights, such as the right of integrity, are not licensed under this Public License, nor are publicity, privacy, and/or other similar personality rights; however, to the extent possible, the Licensor waives and/or agrees not to assert any such rights held by the Licensor to the limited extent necessary to allow You to exercise the Licensed Rights, but not otherwise.
-
-          2. Patent and trademark rights are not licensed under this Public License.
-
-          3. To the extent possible, the Licensor waives any right to collect royalties from You for the exercise of the Licensed Rights, whether directly or through a collecting society under any voluntary or waivable statutory or compulsory licensing scheme. In all other cases the Licensor expressly reserves any right to collect such royalties, including when the Licensed Material is used other than for NonCommercial purposes.
-
-Section 3 – License Conditions.
-
-Your exercise of the Licensed Rights is expressly made subject to the following conditions.
-
-     a.	Attribution.
-
-          1. If You Share the Licensed Material (including in modified form), You must:
-
-               A. retain the following if it is supplied by the Licensor with the Licensed Material:
-
-                    i.	identification of the creator(s) of the Licensed Material and any others designated to receive attribution, in any reasonable manner requested by the Licensor (including by pseudonym if designated);
-
-                    ii.	a copyright notice;
-
-                    iii. a notice that refers to this Public License;
-
-                    iv.	a notice that refers to the disclaimer of warranties;
-
-                    v.	a URI or hyperlink to the Licensed Material to the extent reasonably practicable;
-
-               B. indicate if You modified the Licensed Material and retain an indication of any previous modifications; and
-
-               C. indicate the Licensed Material is licensed under this Public License, and include the text of, or the URI or hyperlink to, this Public License.
-
-          2. You may satisfy the conditions in Section 3(a)(1) in any reasonable manner based on the medium, means, and context in which You Share the Licensed Material. For example, it may be reasonable to satisfy the conditions by providing a URI or hyperlink to a resource that includes the required information.
-
-          3. If requested by the Licensor, You must remove any of the information required by Section 3(a)(1)(A) to the extent reasonably practicable.
-
-          4. If You Share Adapted Material You produce, the Adapter's License You apply must not prevent recipients of the Adapted Material from complying with this Public License.
-
-Section 4 – Sui Generis Database Rights.
-
-Where the Licensed Rights include Sui Generis Database Rights that apply to Your use of the Licensed Material:
-
-     a.	for the avoidance of doubt, Section 2(a)(1) grants You the right to extract, reuse, reproduce, and Share all or a substantial portion of the contents of the database for NonCommercial purposes only;
-
-     b.	if You include all or a substantial portion of the database contents in a database in which You have Sui Generis Database Rights, then the database in which You have Sui Generis Database Rights (but not its individual contents) is Adapted Material; and
-
-     c.	You must comply with the conditions in Section 3(a) if You Share all or a substantial portion of the contents of the database.
-For the avoidance of doubt, this Section 4 supplements and does not replace Your obligations under this Public License where the Licensed Rights include other Copyright and Similar Rights.
-
-Section 5 – Disclaimer of Warranties and Limitation of Liability.
-
-     a.	Unless otherwise separately undertaken by the Licensor, to the extent possible, the Licensor offers the Licensed Material as-is and as-available, and makes no representations or warranties of any kind concerning the Licensed Material, whether express, implied, statutory, or other. This includes, without limitation, warranties of title, merchantability, fitness for a particular purpose, non-infringement, absence of latent or other defects, accuracy, or the presence or absence of errors, whether or not known or discoverable. Where disclaimers of warranties are not allowed in full or in part, this disclaimer may not apply to You.
-
-     b.	To the extent possible, in no event will the Licensor be liable to You on any legal theory (including, without limitation, negligence) or otherwise for any direct, special, indirect, incidental, consequential, punitive, exemplary, or other losses, costs, expenses, or damages arising out of this Public License or use of the Licensed Material, even if the Licensor has been advised of the possibility of such losses, costs, expenses, or damages. Where a limitation of liability is not allowed in full or in part, this limitation may not apply to You.
-
-     c.	The disclaimer of warranties and limitation of liability provided above shall be interpreted in a manner that, to the extent possible, most closely approximates an absolute disclaimer and waiver of all liability.
-
-Section 6 – Term and Termination.
-
-     a.	This Public License applies for the term of the Copyright and Similar Rights licensed here. However, if You fail to comply with this Public License, then Your rights under this Public License terminate automatically.
-
-     b.	Where Your right to use the Licensed Material has terminated under Section 6(a), it reinstates:
-
-          1. automatically as of the date the violation is cured, provided it is cured within 30 days of Your discovery of the violation; or
-
-          2. upon express reinstatement by the Licensor.
-
-     For the avoidance of doubt, this Section 6(b) does not affect any right the Licensor may have to seek remedies for Your violations of this Public License.
-
-     c.	For the avoidance of doubt, the Licensor may also offer the Licensed Material under separate terms or conditions or stop distributing the Licensed Material at any time; however, doing so will not terminate this Public License.
-
-     d.	Sections 1, 5, 6, 7, and 8 survive termination of this Public License.
-
-Section 7 – Other Terms and Conditions.
-
-     a.	The Licensor shall not be bound by any additional or different terms or conditions communicated by You unless expressly agreed.
-
-     b.	Any arrangements, understandings, or agreements regarding the Licensed Material not stated herein are separate from and independent of the terms and conditions of this Public License.
-
-Section 8 – Interpretation.
-
-     a.	For the avoidance of doubt, this Public License does not, and shall not be interpreted to, reduce, limit, restrict, or impose conditions on any use of the Licensed Material that could lawfully be made without permission under this Public License.
-
-     b.	To the extent possible, if any provision of this Public License is deemed unenforceable, it shall be automatically reformed to the minimum extent necessary to make it enforceable. If the provision cannot be reformed, it shall be severed from this Public License without affecting the enforceability of the remaining terms and conditions.
-
-     c.	No term or condition of this Public License will be waived and no failure to comply consented to unless expressly agreed to by the Licensor.
-
-     d.	Nothing in this Public License constitutes or may be interpreted as a limitation upon, or waiver of, any privileges and immunities that apply to the Licensor or You, including from the legal processes of any jurisdiction or authority.
-
-Creative Commons is not a party to its public licenses. Notwithstanding, Creative Commons may elect to apply one of its public licenses to material it publishes and in those instances will be considered the “Licensor.” Except for the limited purpose of indicating that material is shared under a Creative Commons public license or as otherwise permitted by the Creative Commons policies published at creativecommons.org/policies, Creative Commons does not authorize the use of the trademark “Creative Commons” or any other trademark or logo of Creative Commons without its prior written consent including, without limitation, in connection with any unauthorized modifications to any of its public licenses or any other arrangements, understandings, or agreements concerning use of licensed material. For the avoidance of doubt, this paragraph does not form part of the public licenses.
-
-Creative Commons may be contacted at creativecommons.org.
+  b. To the extent possible, if any provision of this Public License is
+     deemed unenforceable, it shall be automatically reformed to the
+     minimum extent necessary to make it enforceable. If the provision
+     cannot be reformed, it shall be severed from this Public License
+     without affecting the enforceability of the remaining terms and
+     conditions.
+
+  c. No term or condition of this Public License will be waived and no
+     failure to comply consented to unless expressly agreed to by the
+     Licensor.
+
+  d. Nothing in this Public License constitutes or may be interpreted
+     as a limitation upon, or waiver of, any privileges and immunities
+     that apply to the Licensor or You, including from the legal
+     processes of any jurisdiction or authority.
+
+=======================================================================
+
+Creative Commons is not a party to its public
+licenses. Notwithstanding, Creative Commons may elect to apply one of
+its public licenses to material it publishes and in those instances
+will be considered the “Licensor.” The text of the Creative Commons
+public licenses is dedicated to the public domain under the CC0 Public
+Domain Dedication. Except for the limited purpose of indicating that
+material is shared under a Creative Commons public license or as
+otherwise permitted by the Creative Commons policies published at
+creativecommons.org/policies, Creative Commons does not authorize the
+use of the trademark "Creative Commons" or any other trademark or logo
+of Creative Commons without its prior written consent including,
+without limitation, in connection with any unauthorized modifications
+to any of its public licenses or any other arrangements,
+understandings, or agreements concerning use of licensed material. For
+the avoidance of doubt, this paragraph does not form part of the
+public licenses.
+
+Creative Commons may be contacted at creativecommons.org.

README.md

@@ -1,2 +1,23 @@
-# SeismologicalDataAnalysis2025
+<!--
+ SeismologicalDataAnalysis2025 (c) by kasper
+ 
+ SeismologicalDataAnalysis2025 is licensed under a
+ Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
+ 
+ You should have received a copy of the license along with this
+ work. If not, see <https://creativecommons.org/licenses/by-nc-sa/4.0/>.
+-->
 
+# Seismological Data Analysis 2025
+
+Jupyter Notebooks for class "Seismological Data Analysis" (summer term 2025)
+
+## License
+
+The content of this repository is licensed under Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License (CC-BY-NC-SA 4.0) ***unless otherwise stated*** (see below)
+
+## Lectures
+
+### 01 - Python Introduction
+
+This has been forked from [Seismo Live](http://seismo-live.org). The source code is available at [https://github.com/krischer/seismo_live](https://github.com/krischer/seismo_live). It is licensed under a ***CC BY-NC-SA 4.0 License***. &copy; 2015-2019 Lion Krischer.

SeismologicalDataAnalysis2025.code-workspace

@@ -0,0 +1,14 @@
+{
+	"folders": [
+		{
+			"path": "."
+		},
+		{
+			"path": "../DataAnalysis2022"
+		}
+	],
+	"settings": {
+		"licenser.author": "Kasper D. Fischer <kasper.fischer@rub.de>",
+		"licenser.license": "CC-BY-NC-SA-4"
+	}
+}

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

images/header.svg

@@ -0,0 +1,28 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<svg width="850.39pt" height="428.03pt" version="1.1" viewBox="158.17 116.17 850.39 428.03" xmlns="http://www.w3.org/2000/svg">
+ <g transform="translate(-3.4095 -2.8898)" fill="none">
+  <title>Canvas 1</title>
+  <g>
+   <title>Layer 1</title>
+   <rect x="161.57" y="119.06" width="850.39" height="428.03" fill="#d2314b"/>
+   <text transform="translate(201.57 168.57)" fill="#ffffff" font-family="'Helvetica Neue'" style="line-height:0%"><tspan x="0" y="33" font-size="35px" font-weight="500">Seismo-Live</tspan><tspan x="0" y="84" fill-opacity=".52" font-size="20px" font-weight="300">Live Jupyter Notebooks</tspan><tspan x="0" y="108" fill-opacity=".52" font-size="20px" font-weight="300">for Seismology</tspan></text>
+   <rect x="731.34" y="119.06" width="280.63" height="428.03" fill="#6ec2ab"/>
+   <text transform="translate(771.34 309.07)" fill="#d3d3d3"/>
+   <path d="m507.32 119.06h106.54l-26.634 428.03h-106.54z" fill="#56be2a"/>
+   <path d="m590.24 119.06h106.54l-26.634 428.03h-106.54z" fill="#277d00"/>
+   <path d="m602.8 119.06h106.54l-26.634 428.03h-106.54z" fill="#3baa04"/>
+   <path d="m615.37 119.06h106.54l-26.634 428.03h-106.54z" fill="#78dc4e"/>
+   <path d="m630.44 119.06h106.54l-26.634 428.03h-106.54z" fill="#9feb7e"/>
+   <path d="m643.66 119.06h106.54l-26.634 428.03h-106.54z" fill="#de7036"/>
+   <path d="m726.57 119.06h106.54l-26.634 428.03h-106.54z" fill="#913500"/>
+   <path d="m739.14 119.06h106.54l-26.634 428.03h-106.54z" fill="#c64f12"/>
+   <path d="m751.7 119.06h106.54l-26.634 428.03h-106.54z" fill="#ff985f"/>
+   <path d="m766.78 119.06h106.54l-26.634 428.03h-106.54z" fill="#ffb78d"/>
+   <path d="m781.85 119.06h106.54l-26.634 428.03h-106.54z" fill="#258468"/>
+   <path d="m864.77 119.06h106.54l-26.634 428.03h-106.54z" fill="#0b5840"/>
+   <path d="m877.33 119.06h106.54l-26.634 428.03h-106.54z" fill="#137757"/>
+   <path d="m889.89 119.06h106.54l-26.634 428.03h-106.54z" fill="#3e9d80"/>
+   <path d="m904.97 119.06h106.54l-26.634 428.03h-106.54z" fill="#6ec2ab"/>
+  </g>
+ </g>
+</svg>