new assignment 6

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