DataAnalysis2021/05-Spectrogram/time_frequency_analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Time Frequency Analysis\n",
    "## Moving window analysis\n",
    "One way to analyse the time-varying frequency content of a signal is to\n",
    "apply windows in the time domain to the signal and to calculate a Fourier spectrum\n",
    "of the windowed part. The window marches along the signal with defined overlap creating\n",
    "a series of Fourier spectra associated with the center times of the windows. The resulting amplitude\n",
    "spectra are the plotted versus window center time. In more detail:\n",
    "\n",
    "1. Choose windowing functions: $w(t,t_m)$ with $t_m$ the center of the window.\n",
    "2. Multiply windowing function with time series: $f_m(t) = f(t)w(t,t_m)$\n",
    "3. Detrend the windowed signal.\n",
    "4. Perform a DFT: $F_{km} = \\Delta t\\sum_{n=0}^N f_m(t)\\exp(-2\\pi i \\frac{kn}{N})$\n",
    "    and calculate the absolute value, $|F_{km}|$.\n",
    "5. Plot the resulting matrix: $|F_{km}|$ in the time-frequency plane."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Borehole Data and Test Site\n",
    "\n",
    "Within the scope of constructing new buildings for the International Geothermal Center in Bochum in 2013 wells were drilled for the installation of a geothermal heat exchanger. Bore holes were drilled next to the newly constructed building (Station GEO3). A downhole hammer was used with a diameter of 152 mm. The drill bit operates by water flushing through the drill rod. The water flow rate determines the working frequency of the hammer.\n",
    "\n",
    "To observe the drill bit noise a temporarily operating 2-D seismic network was installed around the drill site. Here, the noise of the used downhole hammer is investigated. An array of 16 seismological stations was installed in the test site. Four Mark L-4C-3D 1 Hz sensors, eight S-13 1 Hz sensors, one GS-13 1 Hz sensor and two Güralp CMG-3ESPC broad-band 120 sec – 50 Hz sensors were in use. Additionally an accelerometer was fixed to the drill rod (GEO11, blue triangle).\n",
    "\n",
    "Some of the stations were positioned within one of the infrastructure tunnels servicing the university containing water conduits, long-distance heat line and electric cables (e.g. Station GEO4 and GEO05). Thus, noise could be reduced that might disturb the recordings. Other stations are located within buildings (e.g. Station GEO2 and GEO03) ore outside (e.g. Station GEO6 and GEO07).\n",
    "\n",
    "![Map ot the stations at GZB](../images/karte3.jpg)\n",
    "\n",
    "Drill bit noise was recorded up to the maximum drilling depth of 200 m. A drilling cycle is characterised by switching on the water pump, followed by the drilling with higher amplitude signals, that lasts several minutes. The water pump\n",
    "is stopped about 5 to 15 minutes after the drilling finished depending on the drill depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preparation: load packages\n",
    "import os\n",
    "from obspy.core import read\n",
    "from obspy.core import UTCDateTime\n",
    "import numpy as np\n",
    "import matplotlib.pylab as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data from stations, Jan 23 2013\n",
    "\n",
    "record = '1400'                                     # either use record starting at 1300 or 1400\n",
    "station = 'GEO7'                                    # station code\n",
    "#  available stations are GEO2, GEO3 (only 14:00), GEO4, GEO5, GEO6, GEO7\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:17:00Z')      # use record starting at 1400 (14:14-14:23)\n",
    "    etime = UTCDateTime('2013-01-23 14:19:00Z')\n",
    "    ttitle = ', depth: 56.5 m'\n",
    "    \n",
    "datapath = os.path.join(os.path.expanduser('~'),'work', 'data','GE-stations', station, 'e*'+record+'*.HHZ')    # search string for database\n",
    "st = read(datapath)                                  # read file using obspy read\n",
    "print(st)\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\n",
    "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": [
    "# Cell 2: Preparation of window\n",
    "wlen = 4.                                    # window length in seconds\n",
    "npts = tr.stats.npts                         # number of samples\n",
    "dt = tr.stats.delta                          # sampling interval\n",
    "\n",
    "nwin = wlen/dt                                           # number of samples in window\n",
    "nwin = int(pow(2, np.ceil(np.log(nwin)/np.log(2))))      # find next higher power of 2 of nwin\n",
    "nwmax = int(pow(2, np.ceil(np.log(npts)/np.log(2))))/8   # maximum window length (about 1/8 of length)\n",
    "nwin = min(nwin,nwmax)                                   # limit nwin\n",
    "print(\"Samples in window: \",nwin)\n",
    "print(\"Total number of samples: \",npts)\n",
    "print(\"Sampling interval: \",dt)\n",
    "print(\"Length of trace [s]: \",npts*dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The Hann window function\n",
    "def hann(nw):\n",
    "    arg = 2.*np.pi*np.arange(0,nw)/nw     # argument of cosine\n",
    "    fwin = 0.5*(1.-np.cos(arg))           # Hann window\n",
    "    return fwin\n",
    "\n",
    "# The boxcar window function\n",
    "def boxcar(nw):\n",
    "    fwin = np.ones(nw)\n",
    "    return fwin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the window function and its effect on the time series\n",
    "#\n",
    "fwin = hann(nwin)           # Hann window\n",
    "plt.figure(figsize = [15,6])\n",
    "plt.subplot(221)\n",
    "plt.plot(fwin,'-k')\n",
    "plt.title(\"Hann window\")\n",
    "seg = tr.data[0:nwin]*fwin  # multiply data with window\n",
    "plt.subplot(222)\n",
    "plt.plot(seg,'-k')\n",
    "fwin = boxcar(nwin)         # Boxcar window\n",
    "plt.subplot(223)\n",
    "plt.plot(fwin,'-k')\n",
    "plt.title(\"Boxcar window\")\n",
    "seg = tr.data[0:nwin]*fwin  # multiply data with window\n",
    "plt.subplot(224)\n",
    "plt.plot(seg,'-k')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def movingWindowAnalysis(data,winfun,nwin,shift):\n",
    "    \"\"\"\n",
    "    Performs moving window analysis of a time series.\n",
    "    data:   data array\n",
    "    winfun: name of the window function to be called\n",
    "    nwin:   number of window samples (power of 2)\n",
    "    shift:  displacement of moving window in samples\n",
    "    \"\"\"\n",
    "    fwin = winfun(nwin)                    # compute window values\n",
    "    npts = len(data)                       # number of total samples\n",
    "    nseg = int((npts-nwin)/shift)+1        # total number of expected data segment\n",
    "    mwa = np.zeros((nwin//2+1,nseg))       # array for result (rfft returns N/2+1 samples)  \n",
    "    wa = 0                                 # start index of data segment\n",
    "    we = nwin                              # end index of data segment\n",
    "    jseg = 0                               # initialize data segment counter\n",
    "    while we < npts:                       # loop over segments\n",
    "        seg = data[wa:we]*fwin                     # multiply data segment with window\n",
    "        seg = seg-seg.mean()                       # subtract mean value of segment\n",
    "        ftseg = np.abs(np.fft.rfft(seg))           # abs value of Fourier transform\n",
    "        maxft = np.amax(ftseg)                     # max value of Fourier transform\n",
    "        ftseg = ftseg/maxft+1.e-10                 # normalize spectrum to its maximum, remove zeros\n",
    "        mwa[:,jseg] = np.power(ftseg,0.25)         # assign values to the matrix\n",
    "        wa = wa+shift                              # move window start by shift\n",
    "        we = we+shift                              # move window end by shift\n",
    "        jseg = jseg+1                              # increase segment counter\n",
    "    return nseg,mwa                         # return number of segments and moving window matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# carry out the moving window analysis\n",
    "tshift = 0.2*wlen\n",
    "shift = int(tshift/dt)                                     # displacement of moving window in samples\n",
    "nseg,mwa = movingWindowAnalysis(tr.data,hann,nwin,shift)   # compute spectrogram\n",
    "freq = np.fft.rfftfreq(nwin,dt)                            # Fourier frequencies\n",
    "print(\"Frequency range [Hz]: \",freq[0],freq[-1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the result\n",
    "fmaxplot = 250.                                            # maximum frequency to be plotted\n",
    "nf = len(freq)                                             # number of frequencies\n",
    "jfmax = int(fmaxplot/freq[-1]*nf)                          # max index for plotting   \n",
    "extent = (0.5*wlen,0.5*wlen+nseg*tshift,fmaxplot,0.)       # extent of matrix in true time and frequency\n",
    "plt.figure(figsize = [15,6])\n",
    "plt.imshow(mwa[0:jfmax,:],extent = extent,aspect=0.25)      # do plotting\n",
    "plt.xlabel('Window center time [s]')\n",
    "plt.ylabel('Frequency [Hz]')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise \n",
    "\n",
    "1. Write a little routine for the Hamming window and do the analysis using this window\n",
    "2. Use the sqrt of the Fourier amplitude spectrum as output and compare\n",
    "3. Use other powers ($< 1$) of the Fourier amplitude (using np.power)\n",
    "4. Use the logarithm of the Fourier amplitude as output and compare (using np.log)\n",
    "5. Extend the frequency range when using the logarithm or a power of the Fourier amplitude\n",
    "6. Do a moving window analysis for the other available stations. Which ones appear good and which ones are bad?\n",
    "7. Try also records starting at 14:00\n",
    "\n",
    "Hint: when you change something in the code, use \"Run all below\" to see the changes."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Hamming window\n",
    "$ w(t) = 0.54 - 0.46 \\cos \\left( \\frac{2 \\pi t}{L} \\right) ~ , t \\in [0, L]$"
   ]
  },
  {
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}