From 98d91879187d558dd45c13fc83d4ed2fc6326742 Mon Sep 17 00:00:00 2001 From: "Kasper D. Fischer" Date: Sun, 18 Apr 2021 19:31:06 +0200 Subject: [PATCH 1/3] add solution for exercisse 01 (python crash course) --- .../Python_Crash_Course-with_solutions.ipynb | 796 ++++++++++++++++++ 1 file changed, 796 insertions(+) create mode 100644 01-Python_Introduction/Python_Crash_Course-with_solutions.ipynb diff --git a/01-Python_Introduction/Python_Crash_Course-with_solutions.ipynb b/01-Python_Introduction/Python_Crash_Course-with_solutions.ipynb new file mode 100644 index 0000000..7508793 --- /dev/null +++ b/01-Python_Introduction/Python_Crash_Course-with_solutions.ipynb @@ -0,0 +1,796 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "
\n", + "
\n", + "
\n", + "
Scientific Python
\n", + "
A super quick crash course
\n", + "
\n", + "
\n", + "
" + ] + }, + { + "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", + "\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": {}, + "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": {}, + "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": {}, + "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": {}, + "outputs": [], + "source": [ + "name = \"kasper\"\n", + "print(name.capitalize())" + ] + }, + { + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "def simple_poly(x, a, b, c):\n", + " return a * x ** 2 + b * x + c\n", + "\n", + "plt.plot(simple_poly(np.linspace(-5, 5), 10, 2, 2))\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": {}, + "outputs": [], + "source": [ + "print(\"99 bottles of beer on the wall, 99 bottles of beer.\")\n", + "for i in range(98, -1, -1):\n", + " print(\"Take one down, pass it around, %i bottles of beer on the wall.\" % i)" + ] + }, + { + "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": {}, + "outputs": [], + "source": [ + "sentence = \"Pnrfne pvcure? V zhpu cersre Pnrfne fnynq!\"\n", + "sentence=\"Trbculfvpf jvgu Clguba vf sha!\"\n", + "\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", + "result = \"\"\n", + "for letter in sentence:\n", + " if letter not in key:\n", + " result += letter\n", + " else:\n", + " result += key[letter]\n", + "print(result)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Encoder ###\n", + "Idea: Iterrate over all letters of a sentence and find the key (k) which contains the letter as value (v) in the code dictionary" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "sentence = \"Geophysics with Python is fun!\"\n", + "\n", + "result = \"\"\n", + "for letter in sentence:\n", + " encoded_letter = letter\n", + " for k, v in key.items():\n", + " if v == letter:\n", + " encoded_letter = k\n", + " result += encoded_letter\n", + "\n", + "print(result)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +} \ No newline at end of file From 0c0be93f924d34c4f0ef45a727c981864da1407e Mon Sep 17 00:00:00 2001 From: "Kasper D. Fischer" Date: Sun, 18 Apr 2021 20:58:41 +0200 Subject: [PATCH 2/3] add slideshow cell meta-data The notebook can now be viewed in slideshow mode --- .../Python_Crash_Course-with_solutions.ipynb | 384 ++++++++++++++---- 1 file changed, 295 insertions(+), 89 deletions(-) diff --git a/01-Python_Introduction/Python_Crash_Course-with_solutions.ipynb b/01-Python_Introduction/Python_Crash_Course-with_solutions.ipynb index 7508793..aea40d5 100644 --- a/01-Python_Introduction/Python_Crash_Course-with_solutions.ipynb +++ b/01-Python_Introduction/Python_Crash_Course-with_solutions.ipynb @@ -2,7 +2,11 @@ "cells": [ { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "
\n", "
\n", @@ -11,13 +15,8 @@ "
A super quick crash course
\n", "
\n", "
\n", - "" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ + "\n", + "\n", "Seismo-Live: http://seismo-live.org\n", "\n", "##### Authors:\n", @@ -28,14 +27,26 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "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", + "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).**" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ "#### Quick Reference on How to Use This Notebook\n", "\n", "\n", @@ -57,8 +68,17 @@ "\n", "```bash\n", "$ ipython\n", - "```\n", - "\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ "## 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", @@ -71,13 +91,17 @@ "\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" + "* Plots look quite a bit nicer (this is optional)." ] }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# Plots now appear in the notebook.\n", @@ -90,7 +114,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "---\n", "\n", @@ -100,6 +128,8 @@ "\n", "If you have any question regarding some specific Python functionality you can consult the official [Python documenation](http://docs.python.org/).\n", " \n", + "### Tutorials\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", @@ -107,12 +137,24 @@ "* [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", + "### Matlab\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", + "* [NumPy for Matlab Users Cheatsheet](http://mathesaurus.sourceforge.net/matlab-numpy.html)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### More\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", @@ -133,7 +175,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "## Core Python Crash Course\n", "\n", @@ -150,7 +196,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# Three basic types of numbers\n", @@ -170,22 +220,25 @@ }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 2. Strings" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ + "### 2. Strings\n", + "\n", "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": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# You can use single or double quotes to create strings.\n", @@ -210,7 +263,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "#### Exercise\n", "\n", @@ -220,7 +277,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "name = \"kasper\"\n", @@ -229,22 +290,25 @@ }, { "cell_type": "markdown", - "metadata": {}, - "source": [ - "### 3. Lists" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ + "### 3. Lists\n", + "\n", "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": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# List use square brackets and are simple ordered collections of things.\n", @@ -266,7 +330,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "### 4. Dictionaries\n", "\n", @@ -276,7 +344,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# Dictionaries have named fields and no inherent order. As is\n", @@ -304,7 +376,10 @@ { "cell_type": "markdown", "metadata": { - "collapsed": true + "collapsed": true, + "slideshow": { + "slide_type": "slide" + } }, "source": [ "### 5. Functions\n", @@ -315,7 +390,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# Functions are defined using the def keyword.\n", @@ -339,7 +418,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "### 6. Imports\n", "\n", @@ -349,7 +432,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# Import anything, and use it with the dot accessor.\n", @@ -369,7 +456,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "How to know what is available?\n", "\n", @@ -380,7 +471,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "print(dir(math))" @@ -388,7 +483,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "Typing the dot and the TAB will kick off tab-completion." ] @@ -396,7 +495,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "math." @@ -404,7 +507,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "In the IPython framework you can also use a question mark to view the documentation of modules and functions." ] @@ -412,7 +519,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "math.cos?" @@ -420,7 +531,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "### 7. Control Flow\n", "\n", @@ -430,7 +545,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "temp = [\"a\", \"b\", \"c\"]\n", @@ -447,7 +566,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# Useful to know is the range() function.\n", @@ -457,15 +580,23 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ - "The second crucial control flow structure are if/else conditional and they work the same as in any other language." + "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": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "# If/else works as expected.\n", @@ -479,15 +610,27 @@ " print(\"wait what?\")" ] }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "List comprehensions are a nice way to write compact loops. Make sure you understand this as it is very common in Python." + ] + }, { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "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", @@ -503,7 +646,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "### 8. Error Messages\n", "\n", @@ -513,7 +660,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "def do_something(a, b): \n", @@ -524,7 +675,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "## The Scientific Python Ecosystem\n", "\n", @@ -539,7 +694,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "### 9. NumPy\n", "\n", @@ -549,7 +708,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "import numpy as np\n", @@ -571,7 +734,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "### 10. SciPy\n", "\n", @@ -581,7 +748,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "from scipy.interpolate import interp1d\n", @@ -595,7 +766,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "### 11. Matplotlib\n", "\n", @@ -605,7 +780,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -619,7 +798,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, "source": [ "## Exercises\n", "\n", @@ -637,7 +820,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -652,7 +839,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "#### 99 Bottles of Beer\n", "\n", @@ -674,7 +865,12 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "scrolled": true, + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "print(\"99 bottles of beer on the wall, 99 bottles of beer.\")\n", @@ -684,7 +880,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "#### Ceasar Cipher\n", "\n", @@ -714,7 +914,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "outputs": [], "source": [ "sentence = \"Pnrfne pvcure? V zhpu cersre Pnrfne fnynq!\"\n", @@ -739,7 +943,11 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, "source": [ "### Encoder ###\n", "Idea: Iterrate over all letters of a sentence and find the key (k) which contains the letter as value (v) in the code dictionary" @@ -748,7 +956,11 @@ { "cell_type": "code", "execution_count": null, - "metadata": {}, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, "outputs": [], "source": [ "sentence = \"Geophysics with Python is fun!\"\n", @@ -763,16 +975,10 @@ "\n", "print(result)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { + "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3", "language": "python", @@ -793,4 +999,4 @@ }, "nbformat": 4, "nbformat_minor": 1 -} \ No newline at end of file +} From 39d38018790462cb3f099647479b562ef30c4b2d Mon Sep 17 00:00:00 2001 From: "Kasper D. Fischer" Date: Sun, 18 Apr 2021 21:39:19 +0200 Subject: [PATCH 3/3] add files for lecture 02 FFT and Basic Filtering --- .../1-fourier_transform.ipynb | 544 ++++++++++++++++ .../2-filter_basics.ipynb | 601 ++++++++++++++++++ README.md | 6 +- 3 files changed, 1150 insertions(+), 1 deletion(-) create mode 100644 02-FFT_and_Basic_Filtering/1-fourier_transform.ipynb create mode 100644 02-FFT_and_Basic_Filtering/2-filter_basics.ipynb diff --git a/02-FFT_and_Basic_Filtering/1-fourier_transform.ipynb b/02-FFT_and_Basic_Filtering/1-fourier_transform.ipynb new file mode 100644 index 0000000..8183d40 --- /dev/null +++ b/02-FFT_and_Basic_Filtering/1-fourier_transform.ipynb @@ -0,0 +1,544 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "
\n", + "
\n", + "
\n", + "
Signal Processing
\n", + "
Fourier Transformation
\n", + "
\n", + "
\n", + "
\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": [ + "

Tutorial on Fourier transformation in 1D

\n", + "
" + ] + }, + { + "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": [ + "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": [ + "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) / t.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": { + "scrolled": false, + "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(1, high=11, size=N)\n", + "b = np.random.randint(1, high=11, size=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/max(sig), 'r', label='original signal') # normalized \n", + "plt.plot(t, g/max(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 = 150 # 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_)\n", + "\n", + "# plotting\n", + "plt.plot(t_, square, 'r', label='original 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(-1.1,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": { + "scrolled": true, + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "#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", + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/02-FFT_and_Basic_Filtering/2-filter_basics.ipynb b/02-FFT_and_Basic_Filtering/2-filter_basics.ipynb new file mode 100644 index 0000000..1850435 --- /dev/null +++ b/02-FFT_and_Basic_Filtering/2-filter_basics.ipynb @@ -0,0 +1,601 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "
\n", + "
\n", + "
\n", + "
Signal Processing
\n", + "
Filtering Basics
\n", + "
\n", + "
\n", + "
\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", + "\n", + "---" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "# Cell 0 - Preparation: load packages, set some basic options \n", + "#%matplotlib inline\n", + "from obspy import *\n", + "from obspy.clients.fdsn import Client\n", + "import numpy as np\n", + "import matplotlib.pylab as plt\n", + "plt.style.use('ggplot')\n", + "plt.rcParams['figure.figsize'] = 15, 4\n", + "plt.rcParams['lines.linewidth'] = 0.5" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "# Basics in filtering" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## The Filter\n", + "\n", + "In seismology, filters are used to correct for the instrument response, avoid aliasing effects, separate 'wanted' from 'unwanted' frequencies, identify harmonic signals, model a specific recording instrument, and much more ...\n", + "\n", + "There is no clear classification of filters. Roughly speaking, we can distinguish linear vs. non-linear, analog (circuits, resistors, conductors) vs. digital (logical components), and continuous vs. discrete filters. In seismology, we generally avoid non-linear filters, because their output contains frequencies which are not in the input signal. Analog filters can be continuous or discrete, while digital filters are always discrete. Discrete filters can be subdivided into infinite impulse response (IIR) filters, which are recursive and causal, and finite impulse response (FIR) filters, which are non-recursive and causal or acausal. We will explain more about these types of filters below.\n", + "Some filters have a special name, such as Butterworth, Chebyshev or Bessel filters, but they can also be integrated in the classification described above." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Response Funtion\n", + "\n", + "A filter is characterised by its frequency response function $T(j\\omega)$ which relate the Fourier transformation of the output signal $Y(j\\omega)$ and the Fourier transformation of the input signal $X(j\\omega)$:\n", + "$$ Y(j\\omega) = T(j\\omega)X(j\\omega)$$\n", + "\n", + "For a simple lowpass filter it is given as:\n", + "$$ |T(j\\omega)| = \\sqrt{ \\frac{1}{1+(\\frac{\\omega}{\\omega_c})^{2n}} } $$\n", + "\n", + "with $\\omega$ indicating the frequency samples, $\\omega_c$ the corner frequency of the filter, and $n$ the order of the filter (also called the number of corners of the filter). For a lowpass filter, all frequencies lower than the corner frequency are allowed to pass the filter. This is the pass band of the filter. On the other hand, the range of frequencies above the corner frequency is called stop band. In between lies the transition band, a small band of frequencies in which the passed amplitudes are gradually decreased to zero. The steepness of the slope of this transition band is defined by the order of the filter: the higher the order, the steeper the slope, the more effectively 'unwanted' frequencies get removed." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "source": [ + "### Time Domain vs. Frequency Domain\n", + "\n", + "In the time domain, filtering means to convolve the data with the impulse response function of the filter. Doing this operation in the time domain is mathematically complex, computationally expensive and slow. Therefore, the digital application of filters is almost always done in the frequency domain, where it simplifies to a much faster multiplication between data and filter response. The procedure is as follows: transfer the signal into the frequency domain via FFT, multiply it with the filter's frequency response function (i.e. the FFT of the impulse response function), and transfer the result back to the time-domain. As a consequence, when filtering, we have to be aware of the characteristics and pit-falls of the Fourier transformation." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Filter types\n", + "\n", + "There are 4 main types of filters: a lowpass, a highpass, a bandpass, and a bandstop filter. Low- and highpass filters only have one corner frequency, allowing frequencies below and above this corner frequency to pass the filter, respectively. In contrast, bandpass and bandstop filters have two corner frequencies, defining a frequency band to pass and to stop, respectively.\n", + "Here, we want to see how exactly these filters act on the input signal. In Cell 1, the vertical component of the M$_w\\,$9.1 Tohoku earthquake, recorded at Bochum - Germany, is downloaded and pre-processed. In Cell 2, the four basic filters are applied to these data and plotted together with the filter functions and the resulting amplitude spectrum.\n", + "\n", + "1) Look at the figure and explain what the different filters do.\n", + "\n", + "2) Change the order of the filter (i.e the number of corners). What happens and why?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "outputs": [], + "source": [ + "# Cell 1: prepare data from Tohoku-Oki earthquake. \n", + "client = Client(\"BGR\")\n", + "t1 = UTCDateTime(\"2011-03-11T05:00:00.000\")\n", + "st = client.get_waveforms(\"GR\", \"BUG\", \"\", \"LHZ\", t1, t1 + 6 * 60 * 60, \n", + " attach_response = True)\n", + "st.remove_response(output=\"VEL\")\n", + "st.detrend('linear')\n", + "st.detrend('demean')\n", + "st[0].plot();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "print(st[0].stats)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "# Cell 2 - filter types\n", + "npts = st[0].stats.npts # number of samples in the trace\n", + "dt = st[0].stats.delta # sample interval\n", + "fNy = 1. / (2. * dt) # Nyquist frequency \n", + "time = np.arange(0, npts) * dt # time axis for plotting\n", + "freq = np.linspace(0, fNy, npts // 2 + 1) # frequency axis for plotting\n", + "corners = 4 # order of filter\n", + "\n", + "# several filter frequencies for the different filter types\n", + "f0 = 0.04\n", + "fmin1 = 0.04\n", + "fmax1 = 0.07\n", + "fmin2 = 0.03\n", + "fmax2 = 0.07\n", + "\n", + "# filter functions\n", + "LP = 1 / ( 1 + (freq / f0) ** (2 * corners))\n", + "HP = 1 - 1 / (1 + (freq / f0) ** (2 * corners))\n", + "wc = fmax1 - fmin1\n", + "wb = 0.5 * wc + fmin1\n", + "BP = 1/(1 + ((freq - wb) / wc) ** (2 * corners))\n", + "wc = fmax2 - fmin2\n", + "wb = 0.5 * wc + fmin2\n", + "BS = 1 - ( 1 / (1 + ((freq - wb) / wc) ** (2 * corners)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "fragment" + } + }, + "outputs": [], + "source": [ + "# filtered traces\n", + "stHP = st.copy()\n", + "stHP.filter('highpass', freq=f0, corners=corners, zerophase=True)\n", + "stLP = st.copy()\n", + "stLP.filter('lowpass', freq=f0, corners=corners, zerophase=True)\n", + "stBP = st.copy()\n", + "stBP.filter('bandpass', freqmin=fmin1, freqmax=fmax1, corners=corners, zerophase=True)\n", + "stBS = st.copy()\n", + "stBS.filter('bandstop', freqmin=fmin2, freqmax=fmax2, corners=corners, zerophase=True)\n", + "\n", + "# amplitude spectras\n", + "Ospec = np.fft.rfft(st[0].data)\n", + "LPspec = np.fft.rfft(stLP[0].data)\n", + "HPspec = np.fft.rfft(stHP[0].data)\n", + "BPspec = np.fft.rfft(stBP[0].data)\n", + "BSspec = np.fft.rfft(stBS[0].data)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "# ---------------------------------------------------------------\n", + "# plot\n", + "plt.rcParams['figure.figsize'] = 17, 17\n", + "tx1 = 3000\n", + "tx2 = 8000\n", + "fx2 = 0.12\n", + "\n", + "fig = plt.figure()\n", + "\n", + "ax1 = fig.add_subplot(5,3,1)\n", + "ax1.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(time, st[0].data, 'k')\n", + "plt.xlim(tx1, tx2)\n", + "plt.title('time-domain data')\n", + "plt.ylabel('original data \\n amplitude [ms$^-1$]')\n", + "\n", + "ax3 = fig.add_subplot(5,3,3)\n", + "plt.plot(freq, abs(Ospec), 'k')\n", + "plt.title('frequency-domain data \\n amplitude spectrum')\n", + "plt.ylabel('amplitude')\n", + "plt.xlim(0,fx2)\n", + "\n", + "ax4 = fig.add_subplot(5,3,4)\n", + "ax4.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(time, stLP[0].data, 'k')\n", + "plt.xlim(tx1, tx2)\n", + "plt.ylabel('LOWPASS \\n amplitude [ms$^-1$]')\n", + "\n", + "ax5 = fig.add_subplot(5,3,5)\n", + "plt.plot(freq, LP, 'k', linewidth=1.5)\n", + "plt.xlim(0,fx2)\n", + "plt.ylim(-0.1,1.1)\n", + "plt.title('filter function')\n", + "plt.ylabel('amplitude [%]')\n", + "\n", + "ax6 = fig.add_subplot(5,3,6)\n", + "plt.plot(freq, abs(LPspec), 'k')\n", + "plt.ylabel('amplitude ')\n", + "plt.xlim(0,fx2)\n", + "\n", + "ax7 = fig.add_subplot(5,3,7)\n", + "ax7.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(time, stHP[0].data, 'k')\n", + "plt.xlim(tx1, tx2)\n", + "plt.ylabel('HIGHPASS \\n amplitude [ms$^-1$]')\n", + "\n", + "ax8 = fig.add_subplot(5,3,8)\n", + "plt.plot(freq, HP, 'k', linewidth=1.5)\n", + "plt.xlim(0,fx2)\n", + "plt.ylim(-0.1,1.1)\n", + "plt.ylabel('amplitude [%]')\n", + "\n", + "ax9 = fig.add_subplot(5,3,9)\n", + "plt.plot(freq, abs(HPspec), 'k')\n", + "plt.ylabel('amplitude ')\n", + "plt.xlim(0,fx2)\n", + "\n", + "ax10 = fig.add_subplot(5,3,10)\n", + "ax10.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(time, stBP[0].data, 'k')\n", + "plt.xlim(tx1, tx2)\n", + "plt.ylabel('BANDPASS \\n amplitude [ms$^-1$]')\n", + "\n", + "ax11 = fig.add_subplot(5,3,11)\n", + "plt.plot(freq, BP, 'k', linewidth=1.5)\n", + "plt.xlim(0,fx2)\n", + "plt.ylim(-0.1,1.1)\n", + "plt.ylabel('amplitude [%]')\n", + "\n", + "ax12 = fig.add_subplot(5,3,12)\n", + "plt.plot(freq, abs(BPspec), 'k')\n", + "plt.ylabel('amplitude ')\n", + "plt.xlim(0,fx2)\n", + "\n", + "ax13 = fig.add_subplot(5,3,13)\n", + "ax13.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(time, stBS[0].data, 'k')\n", + "plt.xlim(tx1, tx2)\n", + "plt.xlabel('time [sec]')\n", + "plt.ylabel('BANDSTOPP \\n amplitude [ms$^-1$]');\n", + "\n", + "ax14 = fig.add_subplot(5,3,14)\n", + "plt.plot(freq, BS, 'k', linewidth=1.5)\n", + "plt.xlim(0,fx2)\n", + "plt.ylim(-0.1,1.1)\n", + "plt.ylabel('amplitude [%]')\n", + "plt.xlabel('frequency [Hz]')\n", + "\n", + "ax15 = fig.add_subplot(5,3,15)\n", + "plt.plot(freq, abs(BSspec), 'k')\n", + "plt.xlabel('frequency [Hz]')\n", + "plt.ylabel('amplitude ')\n", + "plt.xlim(0,fx2)\n", + "\n", + "plt.subplots_adjust(wspace=0.3, hspace=0.4)\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "## Causal vs. acausal\n", + "\n", + "Filters can be causal or acausal. The output of a causal filter depends only on past and present input, while the output also depends on future input. Thus, an acausal filter is always symmetric and a causal one not. In this exercise, we want to see the effects of such filters on the signal. In Cell 3, the example seismogram of ObsPy is loaded and lowpass filtered several times with different filter order $n$ and causality.\n", + "\n", + "3) Explain the effects of the different filters. You can also play with the order of the filter (variable $ncorners$).\n", + "\n", + "4) Zoom into a small window around the first onset (change the variables $start$, $end$ and $amp$). Which filter would you use for which purpose?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "# Cell 3 - filter effects\n", + "stf = read() # load example seismogram\n", + "tr = stf[0] # select the first trace in the Stream object\n", + "tr.detrend('demean') # preprocess data\n", + "tr.detrend('linear')\n", + "tr.filter(\"highpass\", freq=2) # removing long-period noise\n", + "print(tr)\n", + "t = tr.times() # time vector for x axis\n", + "\n", + "f = 15.0 # frequency for filters (intial: 10 Hz)\n", + "start = 4 # start time to plot in sec (initial: 4)\n", + "end = 8 # end time to plot in sec (initial: 8)\n", + "amp = 1500 # amplitude range for plotting (initial: 1500) \n", + "ncorners = 4 # number of corners/order of the filter (initial: 4)\n", + "\n", + "tr_filt = tr.copy() # causal filter / not zero phase. Order = 2\n", + "tr_filt.filter('lowpass', freq=f, zerophase=False, corners=2)\n", + "tr_filt2 = tr.copy() # causal filter / not zero phase. Order = set by ncorners\n", + "tr_filt2.filter('lowpass', freq=f, zerophase=False, corners=ncorners)\n", + "tr_filt3 = tr.copy() # acausal filter / zero phase. Order = set by ncorners\n", + "tr_filt3.filter('lowpass', freq=f, zerophase=True, corners=ncorners)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "outputs": [], + "source": [ + "# plot - comment single lines to better see the remaining ones\n", + "plt.rcParams['figure.figsize'] = 15, 4\n", + "plt.plot(t, tr.data, 'k', label='original', linewidth=1.)\n", + "plt.plot(t, tr_filt.data, 'b', label='causal, n=2', linewidth=1.2)\n", + "plt.plot(t, tr_filt2.data, 'r', label='causal, n=%s' % ncorners, linewidth=1.2)\n", + "plt.plot(t, tr_filt3.data, 'g', label='acausal, n=%s' % ncorners, linewidth=1.2)\n", + "\n", + "plt.xlabel('time [s]')\n", + "plt.xlim(start, end) \n", + "plt.ylim(-amp, amp)\n", + "plt.ylabel('amplitude [arbitrary]')\n", + "plt.legend(loc='lower right')\n", + "\n", + "plt.show();" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "slide" + } + }, + "source": [ + "### Frequency ranges - bandpass filter\n", + "\n", + "We will now look at an event which took place in Kazakhstan on July 8, 1989. We want to see how filtering helps to derive information from a signal. The signal has been recorded by a Chinese station and is bandpassed in several different frequency bands.\n", + "\n", + "5) What do you see in the different frequency bands?\n", + "\n", + "6) Play with the channel used for filtering in Cell 5. What do you not see? Can you guess what kind of event it is?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "# Cell 4 - get + preprocess data\n", + "c = Client(\"IRIS\")\n", + "tmp1 = UTCDateTime(\"1989-10-19T09:50:00.0\")\n", + "tmp2 = UTCDateTime(\"1989-10-19T10:20:00.0\")\n", + "dat = c.get_waveforms(\"CD\", \"WMQ\", \"\", \"BH*\", tmp1, tmp2, attach_response = True)\n", + "dat.detrend('linear')\n", + "dat.detrend('demean')\n", + "dat.remove_response(output=\"VEL\")\n", + "dat.detrend('linear')\n", + "dat.detrend('demean')\n", + "print(dat)\n", + "dat.plot();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "# Cell 5 - filter data in different frequency ranges\n", + "\n", + "chanel = 2\n", + "tm = dat[chanel].times()\n", + "xmin = 0\n", + "xmax = 700\n", + "\n", + "dat1 = dat[chanel].copy()\n", + "dat2 = dat[chanel].copy()\n", + "dat2.filter(type=\"bandpass\", freqmin=0.01, freqmax=0.05)\n", + "dat3 = dat[chanel].copy()\n", + "dat3.filter(type=\"bandpass\", freqmin=0.05, freqmax=0.1)\n", + "dat4 = dat[chanel].copy()\n", + "dat4.filter(type=\"bandpass\", freqmin=0.1, freqmax=0.5)\n", + "dat5 = dat[chanel].copy()\n", + "dat5.filter(type=\"bandpass\", freqmin=0.5, freqmax=1)\n", + "dat6 = dat[chanel].copy()\n", + "dat6.filter(type=\"bandpass\", freqmin=1., freqmax=5.)\n", + "dat7 = dat[chanel].copy()\n", + "dat7.filter(type=\"bandpass\", freqmin=5., freqmax=9.99)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "slideshow": { + "slide_type": "subslide" + } + }, + "outputs": [], + "source": [ + "plt.rcParams['figure.figsize'] = 17, 21\n", + "fig = plt.figure()\n", + "ax1 = fig.add_subplot(7,1,1)\n", + "ax1.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(tm, dat1.data, 'k')\n", + "plt.xlim(xmin, xmax)\n", + "plt.title('unfiltered')\n", + "plt.ylabel('amplitude \\n [m/s]')\n", + "ax2 = fig.add_subplot(7,1,2)\n", + "ax2.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(tm, dat2.data, 'k')\n", + "plt.xlim(xmin, xmax)\n", + "plt.title('0.01 - 0.05 Hz')\n", + "plt.ylabel('amplitude \\n [m/s]')\n", + "ax3 = fig.add_subplot(7,1,3)\n", + "ax3.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(tm, dat3.data, 'k')\n", + "plt.xlim(xmin, xmax)\n", + "plt.title('0.05 - 0.1 Hz')\n", + "plt.ylabel('amplitude \\n [m/s]')\n", + "ax4 = fig.add_subplot(7,1,4)\n", + "ax4.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(tm, dat4.data, 'k')\n", + "plt.xlim(xmin, xmax)\n", + "plt.title('0.1 - 0.5 Hz')\n", + "plt.ylabel('amplitude \\n [m/s]')\n", + "ax5 = fig.add_subplot(7,1,5)\n", + "ax5.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(tm, dat5.data, 'k')\n", + "plt.xlim(xmin, xmax)\n", + "plt.title('0.5 - 1.0 Hz')\n", + "plt.ylabel('amplitude \\n [m/s]')\n", + "ax6 = fig.add_subplot(7,1,6)\n", + "ax6.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(tm, dat6.data, 'k')\n", + "plt.xlim(xmin, xmax)\n", + "plt.title('1.0 - 5.0 Hz')\n", + "plt.ylabel('amplitude \\n [m/s]')\n", + "ax7 = fig.add_subplot(7,1,7)\n", + "ax7.ticklabel_format(style='sci', axis='y', scilimits=(-1,1))\n", + "plt.plot(tm, dat7.data, 'k')\n", + "plt.xlim(xmin, xmax)\n", + "plt.title('5.0 - 10.0 Hz')\n", + "plt.xlabel('time [sec]')\n", + "plt.ylabel('amplitude \\n [m/s]')\n", + "\n", + "plt.subplots_adjust(hspace=0.3)\n", + "plt.show();" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "# Answers" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "slideshow": { + "slide_type": "skip" + } + }, + "source": [ + "1) The lowpass removes frequencies below the corner frequency of the filter, while a highpass removes all frequencies above the corner frequency. The bandpass allows all frequencies within the two corner frequencies, while a bandstop allows all frequencies outside this range.\n", + "In the filtered signals not only the frequency content is changed (see plot of spectrum) but also the maximum amplitudes (see scaling of y-axis in time-domain plot.) In the frequency-domain plot, we see that each frequency contributes with a distinct amount of amplitude to the final signal. Depending on which frequencies remain after filtering, we have altered the amplitudes as well as the frequencies of the filtered signal.\n", + "Be aware: \"high-\" and \"low-\"pass refers to the unit 'frequency'! When you think in terms of 'period', you have to switch your mind.\n", + "\n", + "2) When setting the order of the filter very high, the slopes of the filter functions turn into steps and the amplitude spectra stop abruptly. Due to the harsh steps in the filter functions, we risk Gibbs effects during filtering. However, when setting the order of the filter lower, the slopes become too gentle and the unwanted frequencies are not removed effectively. Thus, we always have to balance between minimizing the Gibbs effects (low order) and maximizing the effective removal of frequencies (high order).\n", + "\n", + "3) The difference between the two causal filters (red and blue) are minor, even when the order of the filter is increased. But the acausal filter (green) shows a phase shift relative to the original data. This phase shift increases drastically when increasing the order of the filter.\n", + "\n", + "4) Zooming into a time window from 4.5 s to 5 s and setting the maximum amplitude to 400, we see that the phase shift is also shifting the first onset. We would pick it incorrectly. Therefore, whenever the scientific task involves the correct time picking of phases, we have to work with causal filters. When only the frequency content is important, we can also work with acausal filters.\n", + "\n", + "5) In the first two bands (0.01-0.05 Hz) we see Love waves with different amplitudes between 200 and 400 seconds. At 0.1-0.5 Hz we see Rayleigh waves after around 350 seconds and some scatter in the earlier times. At 0.5-1 Hz it is nearly all scattered energy. At 1-5 Hz we see a clear P-wave at around 25 seconds and some later scatter. The last band between 5 and 10 Hz removes the scatter completely with a very nice P-wave coda remaining. We can even distinguish the P-wave of a much smaller event at around 390 seconds.\n", + "\n", + "6) There is no clear S-wave visible. Also in the horizontal components (channel 0 and 1) almost no S-wave energy is visible. This is a clear hint for an explosive type of event. Indeed, this is the recording of a nuclear test explosion. The recording with exactly these bandpass filters was used as cover figure for the book 'Quantitative Seismology' by Keiiti Aki and Paul G. Richards, 2nd edition." + ] + } + ], + "metadata": { + "celltoolbar": "Slideshow", + "kernelspec": { + "display_name": "Python 3", + "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.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/README.md b/README.md index b0fb09c..22f5d3a 100644 --- a/README.md +++ b/README.md @@ -8,4 +8,8 @@ The content of this repository is licensed under Creative Commons Attribution 4. ### 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 (licensed under a ***CC BY-NC-SA 4.0 License***. © 2015-2019 Lion Krischer). \ No newline at end of file +This has been forked from [Seismo Live](http://seismo-live.org). The source code is available at https://github.com/krischer/seismo_live (licensed under a ***CC BY-NC-SA 4.0 License***. © 2015-2019 Lion Krischer). + +### 02 - Fourier-Transform and Basic Filtering + +This has been forked from [Seismo Live](http://seismo-live.org). The source code is available at https://github.com/krischer/seismo_live (licensed under a ***CC BY-NC-SA 4.0 License***. © 2015-2019 Lion Krischer).