new assignment 11

assignment_10/multipleFilterTechnique.ipynb

@@ -0,0 +1,320 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Multiple filter technique"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "An alternative to the moving window method is the multiple filter technique. Instead of moving a window in the time domain, a window is moved in the frequency domain. Of course, a window in the frequency domain is a filter, and hence the name of the method. Let $H_k(\\omega,\\omega_k)$ be the transfer function of the $k$-th filter centered on frequency $\\omega_k$. Applying the filter involves multiplication of the filter transfer function with the Fourier transform of the time signal and then performing an inverse Fourier transform:\n",
+    "\\begin{align}\n",
+    "f_k(t) = \\int_{-\\infty}^\\infty\\,H_k(\\omega,\\omega_k)F(\\omega)\\exp(i\\omega t)\\frac{d\\omega}{2\\pi} \\,.\n",
+    "\\end{align}\n",
+    "Basically, one could plot the resulting time series against the associated center frequency of the filter. However, $f_k(t)$ can also assume negative values. It would be better to plot the envelope of $f_k(t)$.\n",
+    "\n",
+    "### The analytic signal\n",
+    "This can be obtained by constructing the analytic signal $z_k(t)$ of $f_k(t)$. It is a complex\n",
+    "signal whose real part is the original signal and whose imaginary part is its Hilbert transform. Its absolute\n",
+    "value is the envelope or instantaneous amplitude and its phase is the instantaneous phase:\n",
+    "\\begin{align}\n",
+    "z_k(t) = a_k(t)\\exp(i\\Phi_k(t)) = f_k(t)+iq_k(t) \\ \\mathrm{with}\\ q_k(t) = HT\\{f_k(t)\\}\\,,\n",
+    "\\end{align}\n",
+    "and\n",
+    "\\begin{align}\n",
+    "a_k(t) = \\sqrt{f_k(t)^2+q_k(t)^2},\\ \\Phi_k(t) = \\arctan\\left(\\frac{q_k(t)}{f_k(t)}\\right) \\,.\n",
+    "\\end{align}\n",
+    "\n",
+    "The Hilbert transform can again easily be obtained because its Fourier transform is:\n",
+    "\\begin{align}\n",
+    "Q_k(\\omega) = i\\,\\mathrm{sign}(\\omega)F_k(\\omega)\\,.\n",
+    "\\end{align}\n",
+    "where $F_k(\\omega)$ is the Fourier transform of $f_k(t)$. Inverse Fourier transform of $Q_k(\\omega)$ gives the Hilbert transform of $f_k(t)$.\n",
+    "\n",
+    "### Work flow\n",
+    "So the work flow of the multiple filter technique is as follows:\n",
+    "1. Compute Fourier transform of $f(t)$: $F(\\omega)$.\n",
+    "2. Multiply $F(\\omega)$ with bandpass filters $H_k(\\omega,\\omega_k)$ to obtain $F_k(\\omega)$.\n",
+    "3. Compute Fourier transform of the Hilbert transform: $Q_k(\\omega)=iF_k(\\omega)$.\n",
+    "4. Transform both back to the time domain.\n",
+    "5. Calculate the instantaneous amplitude $a_k(t)$ and phase $\\Phi_k(t)$.\n",
+    "6. Plot $a_k(t)$ in the time-frequency plane.\n",
+    "\n",
+    "### Gaussian filters\n",
+    "One frequently made choice for the bandpass filters is a Gaussian function of the form:\n",
+    "\\begin{align}\n",
+    "H_k(\\omega,\\omega_k) = e^{-\\alpha\\left(\\frac{\\omega-\\omega_k}{\\omega_k}\\right)^2}\n",
+    "\\end{align}\n",
+    "with inverse Fourier transform (impulse response):\n",
+    "\\begin{align}\n",
+    "h_k(t) = \\frac{\\sqrt{\\pi}\\omega_k}{2\\alpha}e^{-\\frac{\\omega_k^2t^2}{4\\alpha}}\\cos\\omega_k t \\,.\n",
+    "\\end{align}\n",
+    "This Fourier transform pair demonstrates a fundamental property of spectral analysis:\n",
+    "Choosing $\\alpha$ large gives a very narrow-banded filter but a very\n",
+    "long impulse response. Choosing $\\alpha$ small gives a wide-banded filter but a very short \n",
+    "impulse response. This is an expression of the uncertainty principle. If the frequency of a signal\n",
+    "is known very well, the signal is distributed over time. On the other hand, an impulse in the time domain has\n",
+    "a constant spectrum, i.e. the frequency is entirely unknown. In quantum mechanics, there is the energy-time\n",
+    "uncertainty relation, where energy $E=\\hbar\\omega$.\n",
+    "\n",
+    "### Time and frequency resolution\n",
+    "The Gaussian Fourier transform pair has a special property. If we measure the duration of the impulse response by\n",
+    "\\begin{align}\n",
+    "D_t^2 = \\int_{-\\infty}^\\infty t^2 |h(t)|^2 dt\n",
+    "\\end{align}\n",
+    "and the spread of the filter transfer function by\n",
+    "\\begin{align}\n",
+    "D_\\omega^2 = \\int_{-\\infty}^\\infty \\omega^2 |H(\\omega)|^2 \\frac{d\\omega}{2\\pi}\\,,\n",
+    "\\end{align}\n",
+    "then the product\n",
+    "\\begin{align}\n",
+    "D_t D_\\omega \\ge \\frac{1}{2} \\,.\n",
+    "\\end{align}\n",
+    "For the Gaussian filter, equality holds. Hence, the Gaussian filter is the one with the\n",
+    "best compromise between bandwidth and duration."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "# Preparation: load packages\n",
+    "%matplotlib inline\n",
+    "from obspy.core import read\n",
+    "from obspy.core import UTCDateTime\n",
+    "import numpy as np\n",
+    "import matplotlib.pylab as plt\n",
+    "from setupFigure import SetupFigure"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "def multipleFilterAnalysis(data, alfa, cfreq, dt, ndec):\n",
+    "    \"\"\"\n",
+    "    Perform a multiple filter analysis of data.\n",
+    "    :param data:  Array of detrended and demeaned data \n",
+    "    :param alfa:  Width parameter of Gaussian bandpass filter\n",
+    "    :param cfreq: Array of center frequencies of Gaussian filter\n",
+    "    :param dt: sampling interval of data\n",
+    "    :param ndec: decimation factor for instantaneous amplitude output\n",
+    "    \"\"\"\n",
+    "    nd = len(data)\n",
+    "    ftd = np.fft.rfft(data,nd)                   # Fourier transform of entire data set (pos. freq.)\n",
+    "    freq = np.fft.rfftfreq(nd,dt)                # Fourier frequencies (positive frequencies)\n",
+    "    jf = 0                                       # center frequency counter\n",
+    "    mfa = np.zeros((len(cfreq),nd//ndec))        # numpy array for MFA result\n",
+    "    for cf in cfreq:\n",
+    "        hg = np.exp(-alfa*((freq-cf)/cf)**2)     # Gaussian filter (use f instead of omega here)\n",
+    "        fk = hg*ftd                              # multiply FT of data with filter\n",
+    "        qk = 1j*fk                               # FT of Hilbert transform of filtered data\n",
+    "        ftk = np.fft.irfft(fk)                   # filtered data\n",
+    "        qtk = np.fft.irfft(qk)                   # Hilbert transform of filtered data\n",
+    "        at = np.sqrt(ftk**2+qtk**2)              # instantaneous amplitude\n",
+    "        mfa[jf,:] = at[0:nd:ndec]                # store decimated original result\n",
+    "        jf = jf+1                                # increase center frequency count\n",
+    "    return mfa"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def normalizeMFT(mfa, mode = 'full', exp = 0.5):\n",
+    "    \"\"\"\n",
+    "    Normalize the result of the mutiple filtering operation.\n",
+    "    :param mfa: array with envelopes versus frequency (mfa(f,t))\n",
+    "    :param mode: normalization mode: \n",
+    "                 if 'time': normalize timewise\n",
+    "                 if 'freq': normalize frequencywise\n",
+    "                 else: normalize to absolute maximum of entire array\n",
+    "    :param exp:  exponent for modifying envelope using a power\n",
+    "    \"\"\"\n",
+    "    nf, nt = mfa.shape\n",
+    "    if mode == 'freq':\n",
+    "        for jf in range(0,nf):\n",
+    "            mfamax = np.amax(mfa[jf,:])\n",
+    "            mfa[jf,:] = np.power(mfa[jf,:]/mfamax+1.e-10, exp)          \n",
+    "    elif mode == 'time':\n",
+    "        for jt in range(0,nt):\n",
+    "            mfamax = np.amax(mfa[:,jt])\n",
+    "            mfa[:,jt] = np.power(mfa[:,jt]/mfamax+1.e-10, exp)\n",
+    "    else:        \n",
+    "            mfamax = np.amax(mfa)\n",
+    "            mfa = np.power(mfa/mfamax+1.e-10, exp)\n",
+    "    #        mfa = np.log10(mfa/mfamax+1.e-10)\n",
+    "    return mfa"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "record = '1400'                                     # either use record starting at 1300 or 1400\n",
+    "datapath = 'geo2.' + record+ '.HHZ'\n",
+    "\n",
+    "if record =='1300':\n",
+    "    stime = UTCDateTime('2013-01-23 13:16:00Z')      # use record starting at 1300\n",
+    "    etime = UTCDateTime('2013-01-23 13:25:00Z')\n",
+    "    ttitle = ', depth: 36.5 m'\n",
+    "else:\n",
+    "    stime = UTCDateTime('2013-01-23 14:14:00Z')      # use record starting at 1400 (14:14-14:23)\n",
+    "    etime = UTCDateTime('2013-01-23 14:23:00Z')\n",
+    "    ttitle = ', depth: 56.5 m'\n",
+    "    \n",
+    "st = read(datapath)                                  # read file using obspy read\n",
+    "\n",
+    "st.trim(stime,etime)      # trim data stream to desired start and end time\n",
+    "st.detrend('linear')      # do a linear detrending\n",
+    "st.detrend('demean')      # subtract mean value"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "tr = st[0]                # extract first trace from stream (there is only one)\n",
+    "tr.plot()                 # plot trace"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# prepare multiple filtering\n",
+    "alfa = 10.                                     # filter width parameter\n",
+    "nd = tr.stats.npts                                # number of samples in time series\n",
+    "dt = tr.stats.delta                               # sampling interval\n",
+    "fmax = 1./(2.*dt)                                 # Nyquist frequency\n",
+    "cfreq = np.linspace(1., fmax, 1000)               # array of filter center frequencies\n",
+    "freq = np.fft.rfftfreq(nd, dt)                    # Fourier frequencies"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": true
+   },
+   "outputs": [],
+   "source": [
+    "ndec = 100                                        # decimation factor for envelope\n",
+    "mfa = multipleFilterAnalysis(tr.data, alfa, cfreq, dt, ndec)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "mfa = normalizeMFT(mfa, mode='full', exp = 0.5)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot the result\n",
+    "extent = (0., nd*dt, fmax, 1.)                         # extent of matrix in true time and frequency\n",
+    "fig1, ax1 = SetupFigure(15, 6, \"Time [s]\", \"Center frequency [Hz]\", \"Multiple filter analysis\", 14)\n",
+    "ax1.imshow(mfa, extent = extent, aspect=1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot Gaussian filter\n",
+    "cf = 30\n",
+    "fig2, ax2 = SetupFigure(15, 6, \"Center frequency [Hz]\", \"Amplitude\", \"Gaussian filter\", 14)\n",
+    "arg = (freq-cf)/cf\n",
+    "hg = np.exp(-alfa*((freq-cf)/cf)**2)     # Gaussian filter (use f instead of omega here)\n",
+    "ax2.plot(freq, hg)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot filter response\n",
+    "fig3, ax3 = SetupFigure(15, 6, \"Time [s]\", \"Amplitude\", \"Filter response\", 14)\n",
+    "wk = cf*2.*np.pi\n",
+    "ts = np.linspace(-2.5,2.5,5000)\n",
+    "hres = wk/np.sqrt(np.pi*alfa)*np.exp(-(0.25*(wk*ts)**2/alfa))*np.cos(wk*ts) \n",
+    "ax3.plot(ts, hres)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Tasks"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "1. Study the effect of $\\alpha$ on the result of the multiple filtering and on the form of the filter impulse response.\n",
+    "2. Up to which value of $\\alpha$ do you need to go to achieve a frequency resolution comparable to the moving window method\n",
+    "3. Think about the numerical effort of MFT in comparison with the moving window technique.\n",
+    "4. Try different scalings of the MFT result using different powers and also the logarithm.\n",
+    "5. Try timewise and frequenc-wise normalization.\n",
+    "6. Do you have some ideas why the results look different? What role play scaling and normalization? \n",
+    "7. Plot the Gaussian filter transfer function for different values of $\\alpha$.\n",
+    "8. Plot the filter responses for different values of $\\alpha$."
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

assignment_10/setupFigure.py

@@ -0,0 +1,15 @@
+import matplotlib.pyplot as plt
+
+
+def SetupFigure(wx, wy, xlabel, ylabel, title, tls):
+    """
+    Create a new figure frame
+    """
+    fig = plt.figure(figsize=(wx,wy))
+    ax = fig.add_subplot(1,1,1)
+    ax.set_title(title, fontsize=tls+2)
+    ax.grid(which='major', axis='both', ls=':')
+    ax.set_ylabel(ylabel, fontsize=tls)
+    ax.set_xlabel(xlabel, fontsize=tls)
+    ax.tick_params(labelsize=tls)
+    return fig, ax

assignment_11/p11_surface_wave_group_velocity.ipynb

@@ -0,0 +1,171 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "from obspy import *\n",
+    "from obspy.clients.fdsn import Client\n",
+    "from obspy.clients.fdsn import URL_MAPPINGS\n",
+    "import numpy as np\n",
+    "import matplotlib.pylab as plt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Task 1: Copy functions for the Hann window and the moving window analysis from previous assignments here"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Your function for a Hann window"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Your function for the moving window analysis"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Data from data archive using obspy webfdsn client\n",
+    "FDSN stands for the International Federal of Digital Seismic Networks (www.fdsn.org). Citation from their web site: \"The International Federation of Digital Seismograph Networks (FDSN) is a global organization. Its membership is comprised of groups responsible for the installation and maintenance of seismographs either within their geographic borders or globally. Membership in the FDSN is open to all organizations that operate more than one broadband station. Members agree to coordinate station siting and provide free and open access to their data. This cooperation helps scientists all over the world to further the advancement of earth science and particularly the study of global seismic activity.\"\n",
+    "BGR (Bundesanstalt für Geowissenschaften und Rohstoffe) operates broadband stations in Germany and operates a data archive. It is member of the FDSN. That is why we can freely download data from them.\n",
+    "There are other data archives around the world which are also member of the FDSN and allow opn access to their data (see OBSPY documentation "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "for key in sorted(URL_MAPPINGS.keys()):                       # eine Liste der Archive\n",
+    "    print(\"{0:<7} {1}\".format(key,  URL_MAPPINGS[key]))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Get waveforms for the Tohoku earthquake for station TNS of the German Regional Seismic Network."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "client = Client(\"BGR\")                                          # data archive from which data are fetched\n",
+    "t1 = UTCDateTime(\"2018-01-14T09:19:00.000\")                     # desired start time of data\n",
+    "stc = client.get_waveforms(\"GR\", \"TNS\",\"\",\"LHZ\",t1,             # use fdsn web client to download data\n",
+    "                           t1 + 7200.,attach_response = True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "stc.detrend('linear')                       # take out linear trend\n",
+    "stc.plot()\n",
+    "print(stc)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "jupyter": {
+     "outputs_hidden": false
+    }
+   },
+   "outputs": [],
+   "source": [
+    "st = stc.copy()                                     # proceed with copy to avoid repeated downloading\n",
+    "ta = UTCDateTime(\"2018-01-14T10:00:00.000000Z\")     # cut to surface wave train\n",
+    "te = UTCDateTime(\"2018-01-14T10:29:59.000000Z\")\n",
+    "st.trim(ta,te)\n",
+    "st.taper(0.05,\"hann\")                               # apply a taper to bring signal to zero at ends\n",
+    "st.plot()                                           # dispersion is now well visible\n",
+    "print(st)\n",
+    "tr = st[0]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Task 2: Apply the moving window analysis to the data. \n",
+    "Think about window length and shift time. Try different window lengths."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Task 3: Plot the moving window matrix and the time series below. \n",
+    "Use the same limits for the time axis. Choose a subrange for the frequencies. Does the time-frequency analysis reflect the properties of the time series?"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}