{ "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 }