{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["
\n", "
\n", "
\n", "
Signal Processing
\n", "
Filtering Basics - Solution
\n", "
\n", "
\n", "
"]}, {"cell_type": "markdown", "metadata": {}, "source": ["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": "markdown", "metadata": {}, "source": ["

# Basics in filtering

\n", "
"]}, {"cell_type": "code", "execution_count": null, "metadata": {"code_folding": }, "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": {}, "source": ["## The filter\n", "\n", "A good definition of what a filter is can be found in the book 'Of poles and zeros' by Frank Scherbaum:\n", "\n", "> _Filters or systems are, in the most general sense, devices (in the physical world) or algorithms (in the mathematical world) which act on some input signal to produce a - possibly different - output signal_\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.\n", "\n", "A filter is characterised by its *frequency response function* which is the [Fourier transformation](fourier_transform.ipynb) of the output signal divided by the Fourier transformation of the input signal:\n", "\n", "$$T(j\\omega) = \\frac{Y(j\\omega)}{X(j\\omega)}$$\n", "\n", "For a simple lowpass filter it is given as:\n", "\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. \n", "\n", "In the time domain, filtering means to [convolve](convolution.ipynb) 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](fourier_transform.ipynb)."]}, {"cell_type": "markdown", "metadata": {}, "source": ["---\n", "### 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 Wettzell - Germany, is downloaded and [pre-processed](spectral_analysis+preprocessing.ipynb). 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", "2) Change the order of the filter (i.e the number of corners). What happens and why? "]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# Cell 1: prepare data from Tohoku earthquake. \n", "client = Client(\"BGR\")\n", "t1 = UTCDateTime(\"2011-03-11T05:00:00.000\")\n", "st = client.get_waveforms(\"GR\", \"WET\", \"\", \"BHZ\", 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.plot()"]}, {"cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": ["# Cell 2 - filter types\n", "npts = st.stats.npts # number of samples in the trace\n", "dt = st.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", "# 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)))\n", "\n", "# 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.data)\n", "LPspec = np.fft.rfft(stLP.data)\n", "HPspec = np.fft.rfft(stHP.data)\n", "BPspec = np.fft.rfft(stBP.data)\n", "BSspec = np.fft.rfft(stBS.data)\n", "\n", "# ---------------------------------------------------------------\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.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.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.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.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.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()\n"]}, {"cell_type": "markdown", "metadata": {}, "source": ["---\n", "\n", "### Causal versus 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", "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": {}, "outputs": [], "source": ["# Cell 3 - filter effects\n", "stf = read() # load example seismogram\n", "tr = stf # 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=True, corners=2)\n", "tr_filt2 = tr.copy() # causal filter / not zero phase. Order = set by ncorners\n", "tr_filt2.filter('lowpass', freq=f, zerophase=True, corners=ncorners)\n", "tr_filt3 = tr.copy() # acausal filter / zero phase. Order = set by ncorners\n", "tr_filt3.filter('lowpass', freq=f, zerophase=False, corners=ncorners)\n", "\n", "# 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": {}, "source": ["---\n", "### 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", "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": {}, "outputs": [], "source": ["# Cell 4 - get + preprocess data\n", "c = Client(\"IRIS\")\n", "tmp1 = UTCDateTime(\"1989-07-08T03:40:00.0\")\n", "tmp2 = UTCDateTime(\"1989-07-08T04:05: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": {}, "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=10.)\n", "\n", "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", "plt.subplots_adjust(hspace=0.3)\n", "plt.show()"]}, {"cell_type": "markdown", "metadata": {"tags": ["solution"]}, "source": []}], "metadata": {"kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"}}, "nbformat": 4, "nbformat_minor": 2}