Janis Heuel
a47d87bc9b
Reviewed-on: #13 Co-authored-by: Janis Heuel <janis.heuel@ruhr-uni-bochum.de> Co-committed-by: Janis Heuel <janis.heuel@ruhr-uni-bochum.de>
308 lines
13 KiB
Plaintext
308 lines
13 KiB
Plaintext
{
|
|
"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",
|
|
"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": {
|
|
"scrolled": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# Multiple filter analysis\n",
|
|
"def multipleFilterAnalysis(data,alfa,cfreq,dt,ndec):\n",
|
|
" \"\"\"\n",
|
|
" Perform a multiple filter analysis of data.\n",
|
|
" data: Array of detrended and demeaned data whose length is power of 2 \n",
|
|
" alfa: Width parameter of Gaussian bandpass filter\n",
|
|
" cfreq: Array of center frequencies of Gaussian filter\n",
|
|
" dt: sampling interval of data\n",
|
|
" ndec: decimation factor for instantaneous amplitude output\n",
|
|
" \"\"\"\n",
|
|
" npts = len(data)\n",
|
|
" nd = int(pow(2, np.ceil(np.log(npts)/np.log(2)))) # find next higher power of 2 of npts\n",
|
|
" ftd = np.fft.rfft(data,nd) # Fourier transform of entire data set (pos. freq.)\n",
|
|
" # data are padded with zeros since npts <= nd\n",
|
|
" freq = np.fft.rfftfreq(nd,dt) # Fourier frequencies (positive frequencies)\n",
|
|
" jf = 0 # center frequency counter\n",
|
|
" mfa = np.zeros((len(cfreq),npts//ndec+1)) # 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 = np.complex(0,1)*fk # FT of Hilbert transform\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:npts: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": [
|
|
"# normalize multiple filter result either along time or frequency axis\n",
|
|
"def normalizeMFT(mfa,mode,exp):\n",
|
|
" \"\"\"\n",
|
|
" Normalize the result of the mutiple filtering operation.\n",
|
|
" mfa: array with instantaneous amplitudes versus frequency (mfa(f,t))\n",
|
|
" mode: normalization mode: if 'time', normalize along time axis\n",
|
|
" else normalize along frequency axis\n",
|
|
" exp: exponent for modifying inst amp using a power less than 1\n",
|
|
" \"\"\"\n",
|
|
" nf,nt = mfa.shape\n",
|
|
" if mode == 'time':\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",
|
|
" else:\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",
|
|
" return mfa"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"scrolled": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"# borehole data from stations, Jan 23 2013\n",
|
|
"station = 'GEO2' # station code\n",
|
|
"# available stations are GEO2, GEO3 (only 14:00), GEO4, GEO5, GEO6, GEO7\n",
|
|
"\n",
|
|
"record = '1400' # record starting at 1400\n",
|
|
"stime = UTCDateTime('2013-01-23 14:17:00Z') # use record starting at 1400\n",
|
|
"etime = UTCDateTime('2013-01-23 14:19:00Z')\n",
|
|
"ttitle = ', depth: 56.5 m'\n",
|
|
"\n",
|
|
"# search string for database\n",
|
|
"datapath = os.path.join(os.path.expanduser('~'),'work', 'data', 'GE-stations', station, 'e*'+record+'*.HHZ')\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": [
|
|
"# prepare multiple filtering\n",
|
|
"alfa = 1000.0 # filter width parameter\n",
|
|
"npts = tr.stats.npts # number of samples in time series\n",
|
|
"dt = tr.stats.delta # sampling interval\n",
|
|
"fmax = 1./(2.*dt) # Nyquist frequency\n",
|
|
"nd = int(pow(2, np.ceil(np.log(npts)/np.log(2)))) # find next higher power of 2 of npts\n",
|
|
"print(\"Number of samples: \",npts, # print some information about time series\n",
|
|
" \"\\nSampling interval: \",dt,\n",
|
|
" \"\\nNumber of samples with padding: \",nd,\n",
|
|
" \"\\nMax. frequency: \",fmax)\n",
|
|
"cfreq = np.linspace(1.,fmax,1000) # array of filter center frequencies\n",
|
|
"freq = np.fft.rfftfreq(nd,dt)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {
|
|
"scrolled": true
|
|
},
|
|
"outputs": [],
|
|
"source": [
|
|
"ndec = 100 # decimation factor for inst amplitude\n",
|
|
"mfa = multipleFilterAnalysis(tr.data,alfa,cfreq,dt,ndec)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"mfa = normalizeMFT(mfa,'freq',0.5)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# plot the result\n",
|
|
"extent = (0.,npts*dt,fmax,1.) # extent of matrix in true time and frequency\n",
|
|
"plt.figure(figsize = [15,6])\n",
|
|
"plt.imshow(mfa,extent = extent,aspect=0.25) # do plotting\n",
|
|
"plt.xlabel('Time [s]')\n",
|
|
"plt.ylabel('Center frequency [Hz]')\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# plot Gaussian filters\n",
|
|
"cf = 0.5*fmax\n",
|
|
"plt.figure(figsize = [15,4])\n",
|
|
"for alfa in [10,100,1000,10000,100000]:\n",
|
|
" hg = np.exp(-alfa*((freq-cf)/cf)**2) # Gaussian filter (use f instead of omega here)\n",
|
|
" plt.plot(freq,hg)\n",
|
|
"plt.show()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Exercise"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"1. Study the effect of $\\alpha$ on the result of the multiple filtering\n",
|
|
"2. Think about the numerical effort of MFT in comparison with the moving window technique\n",
|
|
"3. Perform the analysis for different stations\n",
|
|
"4. Do you have some ideas why the results look different? What role plays normalization? Try to change the code such that each spectrum is normalized and not each instantaneous amplitude series."
|
|
]
|
|
},
|
|
{
|
|
"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": 2
|
|
}
|