655 lines
459 KiB
Plaintext
655 lines
459 KiB
Plaintext
|
{
|
||
|
|
"cells": [
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"# Calculation of receiver functions\n",
|
||
|
|
"In this notebook, we calculate a receiver function for a 3-component record of one event at a single station. It is a modification of the minimal example provided by Tom Eulenfeld (Richter) on Github."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 1,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"import numpy as np\n",
|
||
|
|
"import matplotlib.pylab as plt\n",
|
||
|
|
"from obspy.core import read, Stream, Trace, UTCDateTime\n",
|
||
|
|
"from obspy.taup import TauPyModel\n",
|
||
|
|
"from obspy import read_events\n",
|
||
|
|
"from obspy.signal.filter import lowpass\n",
|
||
|
|
"plt.style.use(\"ggplot\")\n",
|
||
|
|
"plt.rcParams['figure.figsize'] = 15, 4"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"We first read the QuakeML file which contains information about the avaiable seismic events using ObsPy. We get a catalogue object whih can be iterated through events and from which information can be extracted. Here, we are interested in the origin time."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 2,
|
||
|
|
"metadata": {
|
||
|
|
"scrolled": false
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"2011-05-15T13:08:15.420000Z\n",
|
||
|
|
"2011-05-13T22:47:55.340000Z\n",
|
||
|
|
"2011-04-30T08:19:16.720000Z\n",
|
||
|
|
"2011-04-18T13:03:04.360000Z\n",
|
||
|
|
"2011-04-07T13:11:23.430000Z\n",
|
||
|
|
"2011-03-31T00:11:58.880000Z\n",
|
||
|
|
"2011-03-06T14:32:36.940000Z\n",
|
||
|
|
"2011-03-01T00:53:45.350000Z\n",
|
||
|
|
"2011-02-25T13:07:26.980000Z\n",
|
||
|
|
"2011-02-21T23:51:42.340000Z\n",
|
||
|
|
"2011-02-21T10:57:51.760000Z\n",
|
||
|
|
"2011-02-12T17:57:56.170000Z\n",
|
||
|
|
"2011-01-31T06:03:26.330000Z\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"cat = read_events(\"example_events.xml\")\n",
|
||
|
|
"for ev in cat:\n",
|
||
|
|
" print(ev.origins[0].time)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"In the minimal example by Tom Eulenfeld, data of the following three events are provided in the SAC format (a fairly widespread format used with SAC (seismological analysis code) software In the following cell, you can choose one of the three events by setting the variable \"evnum\"."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 18,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"[Origin(resource_id=ResourceIdentifier(id=\"smi:service.iris.edu/fdsnws/event/1/query?originid=9925460\"), time=UTCDateTime(2011, 3, 6, 14, 32, 36, 940000), longitude=-27.0253, latitude=-56.3864, depth=92000.0, creation_info=CreationInfo(author='ISC'))]\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"evnum = 1 # coose an event (0,1,2)\n",
|
||
|
|
"avail_events = (cat[1],cat[6],cat[8]) # pick available events from catalogue\n",
|
||
|
|
"event_dates = (\"20110513\",\"20110306\",\"20110225\") # hack for reading associated SAC files\n",
|
||
|
|
"event = avail_events[evnum] # select an event\n",
|
||
|
|
"evdate = event_dates[evnum]\n",
|
||
|
|
"print(event.origins)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"Here, we read the SAC files, detrend and taper the data and then apply a bandpass filter. Afterwards, the traces are trimmed to a 90 seconds time window."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 19,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"st = read(\"data/PB01.BH?.\"+evdate+\".sac\", format = \"SAC\") # read data using ObsPy\n",
|
||
|
|
"st.detrend(\"linear\") # detrending\n",
|
||
|
|
"st.taper(0.05,type = \"hann\") # tapering\n",
|
||
|
|
"st.filter('bandpass', freqmin=0.4, freqmax=1.0) # bandpass filter\n",
|
||
|
|
"dt = st[0].stats.delta\n",
|
||
|
|
"for tr in st: # trim to 90 s window\n",
|
||
|
|
" ta = tr.stats.starttime\n",
|
||
|
|
" tr.trim(ta+5.,ta+95.)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"The \"stats\"-attribute of the data stream contains pretty much details about the traces and the event such as backazimuth, epicentral distance and other useful things."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 20,
|
||
|
|
"metadata": {
|
||
|
|
"scrolled": false
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
" network: CX\n",
|
||
|
|
" station: PB01\n",
|
||
|
|
" location: \n",
|
||
|
|
" channel: BHE\n",
|
||
|
|
" starttime: 2011-03-06T14:40:44.719539Z\n",
|
||
|
|
" endtime: 2011-03-06T14:42:14.719539Z\n",
|
||
|
|
" sampling_rate: 5.0\n",
|
||
|
|
" delta: 0.2\n",
|
||
|
|
" npts: 451\n",
|
||
|
|
" calib: 1.0\n",
|
||
|
|
" _format: SAC\n",
|
||
|
|
" processing: [\"ObsPy 1.1.0: detrend(options={}::type='linear')\", \"ObsPy 1.1.0: taper(max_length=None::max_percentage=0.05::side='both'::type='hann')\", \"ObsPy 1.1.0: filter(options={'freqmin': 0.4, 'freqmax': 1.0}::type='bandpass')\", 'ObsPy 1.1.0: trim(endtime=UTCDateTime(2011, 3, 6, 14, 42, 14, 719539)::fill_value=None::nearest_sample=True::pad=False::starttime=UTCDateTime(2011, 3, 6, 14, 40, 44, 719539))']\n",
|
||
|
|
" sac: AttribDict({'delta': 0.2, 'depmin': -4753.0, 'depmax': 7559.0, 'scale': 1.0, 'b': 0.00053899997, 'e': 120.00054, 'o': -482.77899, 'stla': -21.04323, 'stlo': -69.487396, 'stel': 900.0, 'evla': -56.386398, 'evlo': -27.025299, 'evdp': 92.0, 'mag': 6.5, 'dist': 5242.6309, 'az': 300.62534, 'baz': 149.24416, 'gcarc': 47.148113, 'depmen': 523.57239, 'nzyear': 2011, 'nzjday': 65, 'nzhour': 14, 'nzmin': 40, 'nzsec': 39, 'nzmsec': 719, 'nvhdr': 6, 'npts': 601, 'iftype': 1, 'iztype': 9, 'leven': 1, 'lpspol': 0, 'lovrok': 1, 'lcalda': 1, 'kstnm': 'PB01', 'kcmpnm': 'BHE', 'knetwk': 'CX'})\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwkAAALhCAYAAAAU++jXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcVGX/P/7XsA7DIooCyqaioIBLruVummmZpmZa5nJX\nLm12m7ZpX7U7Ne12yfKTlGWppWamlWYuuYSm5a4IgiAK4gKCDusMDDPX7w9+59wcmBlmQEXo9Xw8\nfBRzljln5gxcr3lf13VUQggBIiIiIiKi/59DTR8AERERERHdXxgSiIiIiIhIgSGBiIiIiIgUGBKI\niIiIiEiBIYGIiIiIiBQYEoiIiIiISIEhgYiIiIiIFBgSiIiIiIhIgSGBiIiIiIgUGBKIiIiIiEiB\nIYHIBkVFRXj++ecRHBwMLy8vPPjggzhy5IhinYULF6JRo0Zo0KAB3nrrLQgh5GVTpkxBixYtoFKp\ncODAAcV2mzdvxoMPPgi1Wo0JEyZUeizHjh1D27ZtodFo0Lt3b6SmpsrL5syZg6CgIHh6eqJly5ZY\nvXp1lfcFAN988w1atmwJDw8PtG7dGhcvXsSCBQvg4eEBDw8PuLq6wtnZWf55ypQpKCkpwVNPPYWg\noCCoVCpcvnzZ7HMfOXIEDg4OmDdvXqXnXFJSgjZt2qBFixZ2Hb85ly9fhpubG1588UWzywcNGgQn\nJyeL23/33Xfy+Xp4eMDNzQ0ODg7IysoCAKxcuRIdOnSAs7Mz5s6dq9hWr9fjtddeg7+/Pxo2bIiZ\nM2cqlqtUKri7u8v7XrBggbzM2nsrhMCsWbPQuHFj1K9fH0888QSuXbtW6WtRFdOnT0doaCg8PT3R\ntm1bbN++XbH8m2++QWBgILy8vPCvf/0LxcXF8rK5c+ciMjISDg4O+OabbxTbxcTEoHfv3vDw8ECf\nPn0qPY6LFy+ie/fu0Gg06NChA86cOSMvi46ORvPmzeHl5YXg4GB8+OGHVd4XAOzYsQNt2rSBu7s7\nQkNDcfjwYcV1oFar4ejoKP88aNAgANY/95LKrsfyyl+fmZmZGD16NBo3bgxvb2/0798f58+ft7i9\nra/zlClToFKpkJ6ebnb5wYMHFZ8Dd3d3qFQqnDhxAoD132tCCMydOxdBQUHw9vbGxIkTFddJ06ZN\nodFoFL9XJJW9t9HR0WjWrBk8PT3x1FNPQavVWjxHIqqEIKJK5efni/fff1+kpqYKo9EoNmzYIHx8\nfEReXp4QQohff/1VBAYGiuTkZHH9+nURFRUlvvzyS3n7lStXiv3794vmzZuL/fv3K/a9d+9e8cMP\nP4hp06aJ8ePHWz0OvV4vAgMDxapVq4ROpxMzZ84UPXr0kJdfuHBB5OTkCCGESExMFP7+/uLs2bNV\n2tf27dtF27ZtRVxcnDCZTCIpKUncunVLsY8PP/ywwjEbDAbx8ccfi8OHDwtXV1dx6dKlCs9tNBpF\n165dRZcuXcQHH3xg9ZyFEGLZsmWie/fuIjQ01Objt+TJJ58U3bp1Ey+88EKFZVu3bhXdu3cXjo6O\nle5HsnDhQtG7d2/FPn7++WcxatQoMWfOHMW6s2fPFr179xa3bt0SGRkZonPnzmLVqlXycgDiypUr\nZp/H2nu7efNmERQUJFJTU0VRUZEYP368GD16tM3nYI85c+aIxMREYTQaxb59+0S9evVESkqKEEKI\ns2fPCm9vb3H06FGh1WpFv379xHvvvSdvu27dOrFz507Rs2dP8fXXXyv2e+zYMfHtt9+KRYsWKV5P\nSzp37ixmz54tdDqd+Oyzz0SzZs2EwWAQQgiRkpIisrKyhBBCXLt2TURERIjt27dXaV+nT58WzZo1\nE0eOHBFGo1FcuXJFXL16VbH9hg0bzB6ztc+9xNr1WJ656/PixYti2bJl4saNG6KkpER89NFHomXL\nlhb3YcvrfOLECdGzZ0+r12N5GzduFE2bNhUmk0kIYf332urVq0VERIRIT08Xubm5YvDgwWLWrFny\n8pCQEHHw4EGzz2Ptvd23b5/w8/MT8fHxQq/Xi0mTJokxY8bYdPxEVBFDAlEVNW7cWBw/flwIIcTo\n0aMVjd2vv/5a9OrVq8I24eHhFhsL5hrc5e3cuVPRUC4oKBBubm5yI62sCxcuCH9/f/Hzzz9XaV9d\nunQRv//+u9XjqeyYLYWElStXiqlTp4rx48dXGhJu3LghWrduLbZv3644Xntei7LbDB06VMyZM6dC\no0yn04nIyEhx8OBBu0JCZGSkIhBKJk+eXCEkdOzYUfF+fPvtt6J79+7yz7Y2ysq/t4sXLxbPPvus\nvPzXX38Vbdu2tfkcquOhhx4SmzdvFkII8c477yhe1/3794vg4OAK2zz66KMVQoLEUoO7rISEBOHu\n7i70er38WEhIiNi3b1+Fda9duyaioqLE8uXLq7Svp59+2uz7a88xW/rcW7sey7P1+tTr9UKlUskN\naXuP2WQyie7du4vjx4/bFRIef/xxRSCUmPsdMWLECMX7cejQIREQECD/bC0klFX+vZ0+fbqYNm2a\nvDw9PV24uLiIgoICm86BiJTY3YioCpKSknDr1i25+0t8fDzatm0rL2/Tpg3i4uLu+POWfx6NRoPQ\n0FDFcy1cuBDu7u4ICwtDQEAA+vfvb/e+jEYjTp48iXPnziEoKAjNmzfHvHnzFF2oqio7Oxsff/wx\n3n///QrLDh06BG9vb8Vjb7/9NmbOnAl3d3ebjx8ofR0GDx4sLy8uLsabb76JJUuWmD2uhQsXYvTo\n0QgMDLT5XE6dOoWLFy9i5MiRNm9T9jUUQlS4Trp06YKAgABMmDAB2dnZFY7R3Hv71FNP4cKFC7h0\n6RJ0Oh02bNiAAQMG2HxMVXX79m2cO3cOERERAMx/DtLS0pCfn39Hnzc+Ph5hYWFwdXVVPFfZ13L9\n+vXw9PREkyZNUFhYaPE9qmxfR48exc2bN9GiRQsEBQXh3//+N4qKiqp9Dtaux7S0NHh7eyMtLU1+\nzNbr8+DBg/Dz84OPjw+A0teh7HtSmTVr1qBVq1bo2LGjzdtkZmZi165dGDt2rM3blP8cXL16FTk5\nOfJjTz31FPz8/DBs2LAK3Qitvbfl91tcXIykpCSbj4uI/ochgchOOp0Ozz33HN59913Uq1cPAJCf\nnw8vLy95HS8vrzveMDL3POae65133kF+fj7++usvjBgxAi4uLnbvKyMjAyUlJdi9ezdiY2Oxb98+\nrF27Ft9++221z2HWrFn497//XSEMAECPHj0UfYiPHDmCpKQkjBkzxq7jB0pfh7L95ZcuXYrHHnsM\noaGhFfZ1+fJlbNq0CTNmzLDrXNatW4ehQ4dWOA5LBg4ciCVLliArKwvXr1/H8uXLUVBQIC+PiYlB\namoqTp8+DZ1OV6Evt6X31t/fH126dEHz5s3h6emJuLg4zJ49265zsZfJZMK//vUvjBgxAq1btwZg\n/nMgPX4n2fI5ePbZZ5GXl4fY2FhMmDABnp6eVdrX1atXsXnzZhw8eBCnTp3CsWPH8NFHH1X7HKxd\nj8HBwdBqtQgODgZg+/WZlZWFyZMnY+HChfJjzz77LM6ePWvTMeXk5GDBggWKsTC22LhxIzp27Iiw\nsDCb1h84cCCio6ORmpoKrVYrjyuQPgvr16/H5cuXkZSUhODgYAwdOhQmk0lxTube24EDB2L9+vU4\nd+4cdDod5s6dC5VKpfiMEZHtGBKI7GAwGDBy5Ei0aNFC0Qjz8PBAbm6u/HNubi48PDyq/XyRkZHy\n4L20tLQKz2PpuVQqFbp27Ypr167hiy++sHtfbm5uAIC33noL3t7eaNq0KSZPnowdO3ZU63ykRtbE\niRMrXddkMmHq1Kn4+OOPoVKpKiy39bUASht6q1evxnvvvWf2uaZNm4YPPvgAarXaxjMBjEYjNmzY\ngHHjxtm8zaxZs9C2bVu0b98e3bt3x/DhwxXfDPfs2RPOzs5o1KgRPvnkE+zYsQN6vV6xD3Pv7fvv\nv4/4+HhkZmYiLy8PvXr
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.figure.Figure at 0x7f00ff600518>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"print(st[0].stats)\n",
|
||
|
|
"st.plot()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"Next, we rotate the ZNE-components into the LQT-system. To do this, we need the backazimuth and the incidence angle of the P-wave. Both are not perfectly known but can be estimated assuming the P-wave propagates in the plane of the source-receiver great circle and using a standard earth model to calculate travel time and incidence angle. The obspy.taup module helps us with this task. "
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 21,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"text/plain": [
|
||
|
|
"3 Trace(s) in Stream:\n",
|
||
|
|
"CX.PB01..BHT | 2011-03-06T14:40:44.719539Z - 2011-03-06T14:42:14.719539Z | 5.0 Hz, 451 samples\n",
|
||
|
|
"CX.PB01..BHQ | 2011-03-06T14:40:44.719538Z - 2011-03-06T14:42:14.719538Z | 5.0 Hz, 451 samples\n",
|
||
|
|
"CX.PB01..BHL | 2011-03-06T14:40:44.719539Z - 2011-03-06T14:42:14.719539Z | 5.0 Hz, 451 samples"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"execution_count": 21,
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "execute_result"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"model = TauPyModel(model=\"iasp91\") # create a model object for ray tracing\n",
|
||
|
|
"evdep = st[0].stats.sac[\"evdp\"] # extract event depth from trace stats\n",
|
||
|
|
"epidis = st[0].stats.sac[\"gcarc\"] # extract epicentral distance in degrees from trace stats\n",
|
||
|
|
"baz = tr.stats.sac[\"baz\"] # extract backazimuth from trace stats\n",
|
||
|
|
"arrivals = model.get_travel_times( # get travel time, incidence angle stored in arrivals-object\n",
|
||
|
|
" source_depth_in_km = evdep,distance_in_degree=epidis,phase_list=\"P\")\n",
|
||
|
|
"incl = arrivals[0].incident_angle # extract incidence angle\n",
|
||
|
|
"ttime = arrivals[0].time # extract arrival time\n",
|
||
|
|
"st.rotate('ZNE->LQT',back_azimuth=baz, inclination=incl) # rotate components"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 22,
|
||
|
|
"metadata": {
|
||
|
|
"scrolled": false
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwkAAADtCAYAAAD0pjmFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XlcVNX/P/DXsM4GosiibCoIypYKYqYsJrmUVkqpiaZm\nrqXlkmtf01xSS9Psk5S5L5UfxcpyFxXIDTQ3QGURAUUQlH2GGZjz+4Pf3A+XmYEZFkl6Px+PeSh3\nOffcbea871mugDHGQAghhBBCCCH/n1FzZ4AQQgghhBDyz0JBAiGEEEIIIYSHggRCCCGEEEIIDwUJ\nhBBCCCGEEB4KEgghhBBCCCE8FCQQQgghhBBCeChIIIQQQgghhPBQkEAIIYQQQgjhoSCBEEIIIYQQ\nwkNBAiGEEEIIIYSHggRC9FBeXo733nsPzs7OsLS0xIsvvogLFy7wllm9ejVsbGzQpk0bzJs3D4wx\nbt7UqVPh5uYGgUCAs2fP8tY7cOAAXnzxRQiFQowfP77OvMTFxcHX1xdisRjBwcG4f/8+N++zzz6D\nk5MTLCws0LlzZ2zbtq3eaQHAjh070LlzZ0ilUnTt2hWpqalYtWoVpFIppFIpzM3NYWpqyv09depU\nVFRU4K233oKTkxMEAgHS09O1bvvChQswMjLCihUr6tzniooK+Pj4wM3NzaD8a5Oeng6RSIT3339f\n6/zBgwfDxMRE5/p79+7l9lcqlUIkEsHIyAh5eXkAgM2bN6NHjx4wNTXF0qVLeevK5XLMmDED9vb2\naNu2LRYtWsSbLxAIIJFIuLRXrVrFzavt3DLGsHTpUjg5OcHKygqTJk2CQqGo81jUx5w5c+Dq6goL\nCwv4+vrijz/+4M3fsWMHHB0dYWlpiQkTJvDysXTpUnh5ecHIyAg7duzgrRcdHY3g4GBIpVKEhITU\nmY/U1FT06dMHYrEYPXr0wPXr17l5ERER6NSpEywtLeHs7Iwvvvii3mkBwJEjR+Dj4wOJRAJXV1ec\nP3+edx0IhUIYGxtzfw8ePBhA7fe9Wl3XY001r8/c3FyMGjUK7dq1g5WVFUJDQ5GUlKRzfX2P89Sp\nUyEQCJCVlaV1fkxMDO8+kEgkEAgEuHLlCoDav9fqul47dOgAsVjM+15Rq+vcRkREoGPHjrCwsMBb\nb72FgoICnftICKkDI4TUqaSkhC1btozdv3+fVVZWsp9++olZW1uz4uJixhhjf/75J3N0dGQpKSks\nOzubeXt7sx9//JFbf/PmzezMmTOsU6dO7MyZM7y0T58+zf773/+yWbNmsXHjxtWaD7lczhwdHdmW\nLVuYTCZjixYtYn379uXm3717lxUWFjLGGLtz5w6zt7dnN27cqFdaf/zxB/P19WUJCQlMpVKx5ORk\n9uTJE14aX3zxhUaelUol27BhAzt//jwzNzdn9+7d09h2ZWUl69WrFwsICGDLly+vdZ8ZY+zrr79m\nffr0Ya6urnrnX5c333yTvfTSS2zixIka8w4dOsT69OnDjI2N60xHbfXq1Sw4OJiXxm+//cZGjhzJ\nPvvsM96yS5YsYcHBwezJkycsJyeH9ezZk23ZsoWbD4BlZmZq3U5t53bbtm3M09OTZWVlsaKiIjZk\nyBC2ePFivffBEJ999hm7c+cOq6ysZFFRUaxVq1YsLS2NMcbYjRs3mJWVFbt8+TIrKChg/fv3Z59+\n+im37u7du9mxY8dYYGAg2759Oy/duLg4tmfPHrZmzRre8dSlZ8+ebMmSJUwmk7HvvvuOdezYkSmV\nSsYYY2lpaSwvL48xxtjDhw+Zp6cn++OPP+qV1rVr11jHjh3ZhQsXWGVlJcvMzGQPHjzgrf/TTz9p\nzXNt971abddjTdquz9TUVPb111+zR48esYqKCrZ27VrWuXNnnWnoc5yvXLnCAgMDa70ea/r5559Z\nhw4dmEqlYozV/r1W1/Xq4uLCYmJitG6ntnMbFRXF7OzsWGJiIpPL5Wzy5MksPDxcr/wTQjRRkEBI\nPbVr147Fx8czxhgbNWoUr7C7fft2FhQUpLGOh4eHzsKCtgJ3TceOHeMVlEtLS5lIJOIKadXdvXuX\n2dvbs99++61eaQUEBLBTp07Vmp+68qwrSNi8eTObOXMmGzduXJ1BwqNHj1jXrl3ZH3/8wcuvIcei\n+jpvvPEG++yzzzQKZTKZjHl5ebGYmBiDggQvLy9eQKg2ZcoUjSDBz8+Pdz727NnD+vTpw/2tb6Gs\n5rkNCwtjGzdu5ObHxsYyBwcHvfehIXr37s0OHDjAGGNswYIFvON65swZ5uzsrLHOwIEDNYIENV0F\n7upu377NJBIJk8vl3DQXFxcWFRWlsezDhw+Zt7c37/gYktaIESO0nl9D8qzrvq/teqxJ3+tTLpcz\ngUDAFaQNzbNKpWJ9+vRh8fHxBgUJr732Gi8gVNP2HVHX9VpbkFBdzXM7Z84cNmvWLG5+VlYWMzMz\nY6WlpXrtAyGEj5obEVIPycnJePLkCdf8JTExEb6+vtx8Hx8fJCQkNPp2a25HLBbD1dWVt63Vq1dD\nIpHA3d0dDg4OCA0NNTityspKXL16Fbdu3YKTkxM6deqEFStW8JpQ1Vd+fj42bNiAZcuWacyLjY2F\nlZUVb9r8+fOxaNEiSCQSvfMPVB2HIUOGcPMVCgU++eQTrFu3Tmu+Vq9ejVGjRsHR0VHvffn777+R\nmpqKt99+W+91qh9DxpjGdRIQEAAHBweMHz8e+fn5GnnUdW5rpvvgwQMUFhbqna/6ePr0KW7dugVP\nT08A2u+DjIwMlJSUNOp2ExMT4e7uDnNzc962qh/Lffv2wcLCAu3bt0dZWZnOc1RXWpcvX8bjx4/h\n5uYGJycnfPzxxygvL2/wPtR2PWZkZMDKygoZGRncNH2vz5iYGNjZ2cHa2hpA1XGofk7qsnPnTnTp\n0gV+fn56r5Obm4vjx49j7Nixeq9T1/X61ltvwc7ODsOGDdNoRljbua2ZrkKhQHJyst75IoT8DwUJ\nhBhIJpNhzJgxWLhwIVq1agUAKCkpgaWlJbeMpaVloxeMtG1H27YWLFiAkpISXLx4EWFhYTAzMzM4\nrZycHFRUVODEiRO4efMmoqKisGvXLuzZs6fB+7B48WJ8/PHHGsEAAPTt25fXhvjChQtITk5GeHi4\nQfkHqo5D9fby69evx6uvvgpXV1eNtNLT07F//37MnTvXoH3ZvXs33njjDY186DJo0CCsW7cOeXl5\nyM7OxsaNG1FaWsrNj46Oxv3793Ht2jXIZDKNtty6zu2gQYMQERGB+/fvo6CggGunXT3txqZSqTBh\nwgSEhYWha9euALTfB+rpjUmf+2D06NEoLi7GzZs3MX78eFhYWNQrrQcPHuDAgQOIiYnB33//jbi4\nOKxdu7bB+1Db9ejs7IyCggI4OzsD0P/6zMvLw5QpU7B69Wpu2ujRo3Hjxg298lRYWIhVq1bx+sLo\n4+eff4afnx/c3d31Wr6u63Xfvn1IT09HcnIynJ2d8cYbb0ClUvH2Sdu5HTRoEPbt24dbt25BJpNh\n6dKlEAgETXofENKSUZBAiAGUSiXefvttuLm5YcmSJdx0qVSKoqIi7u+ioiJIpdIGb8/Ly4vrvJeR\nkaGxHV3bEggE6NWrFx4+fIgffvjB4LREIhEAYN68ebCyskKHDh0wZcoUHDlypEH7oy5kTZo0qc5l\nVSoVZs6ciQ0bNkAgEGjM1/dYAFUFvW3btuHTTz/Vuq1Zs2Zh+fLlEAqFeu4JUFlZiZ9++gnvvvuu\n3ussXrwYvr6+6NatG/r06YPhw4fzngwHBgbC1NQUNjY2+Oabb3DkyBHI5XJeGtrO7XvvvYcRI0Yg\nODgY3t7eCA0NhampKezs7OrMU/XOp4aYPn06CgsLERERwUur5n2gnt4QgwcP5vKo7jCr77n39vaG\nWCzG559/Xq+0RCIRZsyYgXbt2qFt27aYPXt2g++Duq7HmvS5PouLizF48GCMHDkS48aNq1e+li5d\niilTpsDW1tag9Xbv3m3
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.figure.Figure at 0x7f00ff6a5828>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"st[2].plot() # plot new L-component"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"Next, we trim the traces to a 40 s window around the P-onset which we estimate by adding the travel time to the event origin time. We could also pick the arrival time. The resulting trace is again tapered to avoid artifacts in the discrete Fourier transform."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 23,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"2011-03-06T14:40:59.816249Z\n",
|
||
|
|
"3 Trace(s) in Stream:\n",
|
||
|
|
"CX.PB01..BHT | 2011-03-06T14:40:49.719539Z - 2011-03-06T14:41:29.719539Z | 5.0 Hz, 201 samples\n",
|
||
|
|
"CX.PB01..BHQ | 2011-03-06T14:40:49.719538Z - 2011-03-06T14:41:29.719538Z | 5.0 Hz, 201 samples\n",
|
||
|
|
"CX.PB01..BHL | 2011-03-06T14:40:49.719539Z - 2011-03-06T14:41:29.719539Z | 5.0 Hz, 201 samples\n"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAwIAAADtCAYAAAAMRsJyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3XdYFFf7P/73UneXBcEG0kQwNsTHgiUKluijYjSKJnZj\nwRajPrEnYCxRY/lEY4pKTGzRGGMQjF0TiaKJNWhQbCgqVVQMdXdh2b2/f/jb+bGwsAUQI/fruva6\n2Dlnzjkzs7PMvXPOGRERERhjjDHGGGM1ikV1N4AxxhhjjDH24nEgwBhjjDHGWA3EgQBjjDHGGGM1\nEAcCjDHGGGOM1UAcCDDGGGOMMVYDcSDAGGOMMcZYDcSBAGOMMcYYYzUQBwKMMcYYY4zVQBwIMMYY\nY4wxVgNxIMAYY4wxxlgNxIEAY0YoKCjAhAkT4OnpCQcHB3Tq1Annzp3TybNq1SrUq1cPtWvXxvz5\n80FEQtrUqVPRuHFjiEQinDp1Sme9iIgIdOrUCWKxGOPGjTPYlkuXLqFVq1aQSqXo1q0bHj58KKQt\nXrwYHh4esLe3x2uvvYatW7eaXRYAbN++Ha+99hpkMhmaN2+Oe/fu4dNPP4VMJoNMJoOtrS2sra2F\n91OnTkVRURHefvtteHh4QCQS4cGDB3rrPnfuHCwsLLB8+XKD21xUVAQ/Pz80btzYpPbr8+DBA0gk\nEkycOFFYRkRYsmQJPDw84OjoiEmTJqGwsFDv+j/88IOwvTKZDBKJBBYWFnj69CkAYNOmTWjbti2s\nra2xZMkSnXWVSiVmzJgBFxcX1K1bF6GhoTrpIpEIdnZ2QtmffvqpkFbesTWl/RU1Z84c+Pj4wN7e\nHq1atcKhQ4d00rdv3w53d3c4ODhg/PjxOu1YsmQJfH19YWFhge3bt+usFxMTg27dukEmk6F79+4G\n23Hv3j106dIFUqkUbdu2xd9//y2khYeHw9vbGw4ODvD09MTKlSvNLgsAjhw5Aj8/P9jZ2cHHxwd/\n/vmnzudALBbD0tJSeB8UFASg/PNeS9/nsTxBQUGwsrLSWVbefi1pzZo1aNGiBezt7dGkSRNs27ZN\nJz08PByNGjWCvb093n77bWRlZekt58yZMzrngZ2dHUQiEf766y8A5X+vGfq8enl5QSqV6nyvFG9f\necfW2PYzxv4/xBgzKC8vj5YuXUoPHz4ktVpNP/74I9WpU4dyc3OJiOjw4cPk7u5Od+/epfT0dGrZ\nsiV99913wvqbNm2i33//nby9ven333/XKfvkyZP0888/06xZs2js2LHltkOpVJK7uzt9++23pFAo\nKDQ0lAICAoT0O3fuUHZ2NhER3b59m1xcXCguLs6ssg4dOkStWrWi+Ph40mg0lJCQQM+ePdMpY+XK\nlaXarFKpaP369fTnn3+Sra0t3b9/v1TdarWaOnbsSB06dKBly5aVu81ERJ9//jl16dKFfHx8jG5/\nWQYNGkSdO3emkJAQYdnWrVupRYsWlJKSQjk5OdS/f38KCwszWBYR0apVq6hbt27C+6ioKPrll19o\n2LBhtHjxYp28ixYtom7dutGzZ88oIyOD2rdvT99++62QDoCSk5P11lPesa1I+021ePFiun37NqnV\naoqOjqZatWpRYmIiERHFxcWRo6MjXbx4kbKysqhnz560cOFCYd2dO3fSsWPHKDAwkLZt26ZT7qVL\nl2jXrl20evVqnf1Zlvbt29OiRYtIoVDQxo0bqVGjRqRSqYiIKDExkZ4+fUpERGlpadSiRQs6dOiQ\nWWVdvXqVGjVqROfOnSO1Wk3JycmUmpqqs/6PP/6ot83lnfda+j6PZYmKiqIuXbqQpaWlzvLy9mtJ\na9asoStXrlBRURH9/fffVL9+fTp79iwREUVHR5OzszPduHGDlEolTZ48mUaNGmWwXUREe/bsIS8v\nL9JoNERU/veaoc9rw4YN6cyZM3rrKe/YVqT9jNVUHAgwZqYGDRrQ5cuXiYho+PDhOhe027Zto65d\nu5Zap2nTpmVeEOi7qC7p2LFjOhfD+fn5JJFIhAux4u7cuUMuLi70yy+/mFVWhw4d6Lfffiu3PYba\nXFYgsGnTJpo5cyaNHTvWYCDw6NEjat68OR06dEinvabsi+LrDBw4kBYvXqxz4TVkyBD64osvhPdn\nz54lNze3ctul5evrqxP0aU2ZMqVUINCuXTud47Fr1y7q0qWL8L68QKC4kse2Iu2vqNdff50iIiKI\niOjDDz/U2a+///47eXp6llqnT58+ZV6wlnVRXdytW7fIzs6OlEqlsKxhw4YUHR1dKm9aWhq1bNlS\nZ/+YUtbQoUP1Hl9T2lzWeV/W51EfhUJBvr6+dObMmVKBgFZ5+7UsI0aMoM8++4yIiObMmUOzZs0S\n0lJSUsjGxoby8/MNlvPmm2/qBH1a+r4jDH1eywsEiit5bCvSfsZqKu4axJgZEhIS8OzZM6Gryo0b\nN9CqVSsh3c/PD/Hx8ZVeb8l6pFIpfHx8dOpatWoV7Ozs0KRJE7i5uaFXr14ml6VWqxEbG4vr16/D\nw8MD3t7eWL58uU53J3NlZmZi/fr1WLp0aam0s2fPwtHRUWfZggULEBoaCjs7O6PbDzzfD/379xfS\nCwsLMW/ePKxdu1Zvu4pvGxEhNTUV2dnZ5W7LlStXcO/ePbzzzjvl5iuvnpKfkw4dOsDNzQ3jxo1D\nZmamTlp5x9ac9lfUP//8g+vXr6NFixYA9J8HSUlJyMvLq9R6b9y4gSZNmsDW1lanruL7cvfu3bC3\nt4erqyvkcnmZx8hQWRcvXsSTJ0/QuHFjeHh44IMPPkBBQUGFt6G8z2NSUhIcHR2RlJQkLFu1ahWG\nDx8Od3d3k+rZvXu3zjEpTqVS4fz58/D19RWWlfwcFRYWIiEhodw6Hj9+jOPHj2PMmDFGt8vQ5/Xt\nt9+Gs7MzgoODS3X5K+/YmtN+xmoyDgQYM5FCocDo0aPx0UcfoVatWgCAvLw8ODg4CHkcHBwq/eJH\nXz366vrwww+Rl5eH8+fPY8iQIbCxsTG5rIyMDBQVFeHEiRO4du0aoqOj8f3332PXrl0V3oawsDB8\n8MEHpS74ASAgIECnT++5c+eQkJCAUaNGmdR+4Pl+KN5/fd26dejXrx98fHxKldW3b1+Eh4fj4cOH\nyMrKEvod5+fnl7stO3fuxMCBA0u1oyx9+/bF2rVr8fTpU6Snp+OLL77QqSMmJgYPHz7E1atXoVAo\nSvWtLuvYmtv+itBoNBg/fjyGDBmC5s2bA9B/HmiXVyZjzoORI0ciNzcX165dw7hx42Bvb29WWamp\nqYiIiMCZM2dw5coVXLp0CWvWrKnwNpT3efT09ERWVhY8PT0BPB9HsHfvXsydO9fkekaOHIm4uDi9\naXPmzIGXlxf69OkD4PnnaPfu3bh+/ToUCgWWLFkCkUhk8HO0Z88etGvXDk2aNDGqTYY+r7t378aD\nBw+QkJAAT09PDBw4EBqNRmeb9B1bc9vPWE3GgQBjJlCpVHjnnXfQuHFjLFq0SFguk8mQk5MjvM/J\nyYFMJqtwfb6+vsKAuaSkpFL1lFWXSCRCx44dkZaWhs2bN5tclkQiAQDMnz8fjo6O8PLywpQpU3Dk\nyJEKbY/2QmrSpEkG82o0GsycORPr16+HSCQqlW7svgCeX8xt3boVCxcu1FvXhAkTMHToUHTr1g0t\nW7ZEr169YG1tDWdn5zLbp1ar8eOPP+Ldd981uC1aYWFhaNWqFVq3bo0uXbpg8ODBOr/wBgYGwtra\nGvXq1cOXX36JI0eOQKlU6pSh79ia036t4gM+TTFt2jRkZ2cjPDxcp6yS54F2eUUEBQUJbdQOUjX2\n2Lds2RJSqRSffPKJWWVJJBLMmDEDDRo0QN26dTF79uwKnweGPo8lzZo1C8uWLYNYLK5QvcWtXLkS\n0dHRiIiIEM6vXr16YdGiRQgODoaXl5cwKNzQXYidO3eadB4Y+rx27twZYrEYDg4OWLduHe7evYt7\n9+6VKqfksTW3/YzVZBw
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.figure.Figure at 0x7f00ff7dbbe0>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"onset = event.origins[0].time+ttime # onset time = origin time + travel time\n",
|
||
|
|
"print(onset)\n",
|
||
|
|
"st.trim(onset-10.,onset+30.)\n",
|
||
|
|
"st.taper(0.125,type = \"hann\")\n",
|
||
|
|
"print(st)\n",
|
||
|
|
"st[2].plot()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 24,
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": true
|
||
|
|
},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"zcomp = st[2].data # rename traces\n",
|
||
|
|
"rcomp = st[1].data"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Deconvolution in the frequency domain using waterlevel method\n",
|
||
|
|
"Here comes the central routine for receiver function calculation which does the deconvolution, here using the water level method. Note that we add zeros to the traces to avoid the pitfalls when doing cross-correlation inthe frequency domain. We deconvolve both the Z- or L-component and the R-or Q-component."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 25,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"def deconvWaterLevel(rsp, src, dt, waterlevel=0.05, alfa=2.,\n",
|
||
|
|
" tshift=None,length=None,normalize=False):\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" Frequency-domain deconvolution using waterlevel method.\n",
|
||
|
|
" Deconvolve src from rsp.\n",
|
||
|
|
"\n",
|
||
|
|
" :param rsp: trace to be deconvolved (nominator)\n",
|
||
|
|
" :param src: trace to deconvolved with (denominator)\n",
|
||
|
|
" :param dt: sampling interval of data\n",
|
||
|
|
" :param waterlevel: waterlevel to stabilize the deconvolution\n",
|
||
|
|
" :param gauss: Gauss parameter of low-pass filter\n",
|
||
|
|
" :param tshift: shift zero lag of receiver function to the right by tshift\n",
|
||
|
|
" :param length: number of data points in results, optional\n",
|
||
|
|
" :param normalize: normalize to src (default = False)\n",
|
||
|
|
" :return: receiver function\n",
|
||
|
|
" \"\"\"\n",
|
||
|
|
" npts = len(rsp)\n",
|
||
|
|
" if length == None:\n",
|
||
|
|
" nout = npts\n",
|
||
|
|
" else:\n",
|
||
|
|
" nout = length # cut output to nout samples\n",
|
||
|
|
" nfft = int(pow(2, np.ceil(np.log(npts)/np.log(2)))) # next higher power of 2\n",
|
||
|
|
" nfft = 2*nfft # add zeros to avoid pitfalls with cc \n",
|
||
|
|
" freq = np.fft.rfftfreq(nfft,dt) # Fourier frequencies\n",
|
||
|
|
" spq = np.fft.rfft(rsp,nfft) # FT of R or Q component\n",
|
||
|
|
" spz = np.fft.rfft(src,nfft) # FT of Z or L component\n",
|
||
|
|
" spshift = np.exp(-1j*2.*np.pi*freq*tshift) # time shift (applied in the frequency domain)\n",
|
||
|
|
" cspz = np.conjugate(spz) # complex conjugate needed later\n",
|
||
|
|
" spzabs2 = np.abs(spz)**2 # absolute value squared\n",
|
||
|
|
" water = np.maximum(spzabs2,max(spzabs2)*waterlevel) # regularized denominator\n",
|
||
|
|
" gauss = np.exp((-np.pi*freq/alfa)**2) # Gaussian low-pass filter\n",
|
||
|
|
" sprfz = gauss*spshift*spzabs2/water # deconvolve Z or L component \n",
|
||
|
|
" rfz = np.fft.irfft(sprfz,nfft)[:npts] # Z or L receiver function intime domain\n",
|
||
|
|
" if normalize:\n",
|
||
|
|
" norm = 1.0/max(rfz)\n",
|
||
|
|
" else:\n",
|
||
|
|
" norm = 1.0\n",
|
||
|
|
" rfz = norm*rfz # normalize Z-receiver function to max of 1\n",
|
||
|
|
" sprfq = gauss*spshift*spq*cspz/water # spectrum of Q-receiver function\n",
|
||
|
|
" rfq = -np.fft.irfft(sprfq,nfft)*norm # Q-receiver function (normalized and sign reversed)\n",
|
||
|
|
" return rfq[0:nout],rfz[0:nout]"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 26,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"tshift = 5. # shift the zero lag to +5 seconds\n",
|
||
|
|
"rfq,rfl = deconvWaterLevel(rcomp, zcomp, dt, \n",
|
||
|
|
" waterlevel=0.05, alfa=20.,\n",
|
||
|
|
" length=150,tshift=tshift,normalize=True) # do the deconvolution\n",
|
||
|
|
"tr = np.arange(0,len(rfq))*dt-tshift # take shift into account when plotting time axis"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 27,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA4cAAAENCAYAAAC8WRABAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXlcVXX6x9/fAwgqogKioIiA4r7mbpqWaeu0TMu0/qZt\n2mZsapqlxaaaaZmpqWmZlmmbmvbFaTfNrEzNNHdREVRAEARBFGTnPL8/vpebKMtdzgXU7/v16pVy\nz7nfR+65957n+zzP56NERDAYDAaDwWAwGAwGw3GN1dYBGAwGg8FgMBgMBoOh7THJocFgMBgMBoPB\nYDAYTHJoMBgMBoPBYDAYDAaTHBoMBoPBYDAYDAaDAZMcGgwGg8FgMBgMBoMBkxwaDAaDwWAwGAwG\ngwGTHBoMBoPBYDAYDAaDAZMcGgwGg8FgMBgMBoMBkxwaDAaDwWAwGAwGgwGTHBoMBoPBYDAYDAaD\nAQhu6wBag927d7d1CEcQHR3N3r172zoMQxtgXvvjG/P6H9+Y1//4xbz2xzfm9T++aQ+vf1xcnEfH\nmcqhwWAwGAwGg8FgMBhMcmgwGAwGg8FgMBgMBpMcGgwGg+E4x57/PvZrT7d1GAaDwWAwtDkmOTQY\nDAbDcYuUHUA+fQdZ+iVyYF9bh2MwGAwGQ5tikkODwWAwHLfIN59DdRWIIOt+aOtwDAaDwWBoU0xy\naDAYDIbjEqmuQhZ/BsPHQo9eyNoVbR2SwWAwGAxtikkODQZDu8F+5QnsVd+1dRiG4wRZvhhK92PN\nPh81ehJs2YCUH2zrsAwGg8FgaDNMcmgwGNoFUlKELP8KeeUJJGdnW4djOMYRuw758kPoNwBShqJG\nT4C6WmTT6rYOzWAwGAyGNsMkhwaDoX2Qme7+o/3c35HK8jYMxnDMs/YHKMjDOu18lFKQNAgiuoFp\nLTUYDAbDcYxJDg0GQ7tAdqZDUBDWTXdAQR7y32cRkbYOy3AMIiLYC+ZBj14weiIAyrJQoyYgG1cj\nNdVtHKHBYDAYDG2DSQ4NBkO7QDLTIa4vatgJqJ9dgqz8Fln6ZVuHZTgWSd8MO7ehZp2LsoLcP1aj\nJkJVBWxZ34bBGQwGg8HQdpjk0GAwtDkiApnpqMQUANQZF8CQUchb/zbzhwbHsRfMg/AI1ORTGj4w\naASEdTSqpQaDwWA4bmk3yeG6deu45ZZb+M1vfsOHH354xOO5ubncddddXHrppXz88cdtEKHBcHwg\ntbWtP+9XkAflB7U4CKCsIKxrboNO4Wb+8BhDamuw33iOuj9ejZSXtf76u7NhwyrUjDNRHUIbPKZC\nQlDDxyLrVyJ2XavHZjAYDAZDW9MukkPbtnnppZe48847efzxx1m2bBk5OTkNjgkPD+eqq67i7LPP\nbqMoDYZjHxHBfu5h7L/+rnXXdYnRKFdyCKAiumFdd7ueP/zk7VaNxxAYZF8R9iN3auP54r0NRIha\nLYZFH0OHDqgZZzZ+wOhJULofMra2bmAGg8FgMLQD2kVymJGRQa9evejZsyfBwcFMnjyZVatWNTim\na9eu9O/fn6CgoCaexWAw+Iss/RLWr4Q9uUjx3tZbeOc26NAB4vo2+LEaOAxShiIZW1ovFkNAkG2p\n2H+9FXKzUFfcrH+Wtb3140jbBMPGorpENPq4GjYGgoNNa6nBYDAYjkvaRXJYXFxMVFSU++9RUVEU\nFxe3YUQGw/GHFOYj77ykFRxBJ2yttXZmOvRNRjWy+aN6J0ButlEuPUoREeyvPsF+7G7o2Bnrzkex\nps2GHr2QzIzWjaWmBvbmo+LimzxGdewEg0cha78315zBYDAYjjuC2zqAQLBo0SIWLVoEwMMPP0x0\ndHQbR3QkwcHB7TIuQ+Bpj6+91NWx7/G5SJBF5H1PUnTL5YTlZdMl+meBX7u2loJdO+k0+1y6NPJ7\nKU8ZQuniT4mUWoJ6xAY8nkDTHl//QFK+8ENK336B0PFTiZgzF6tzOAAlA4ZQk7GlVX8Xtbt2UmTb\ndBkwmI7NrFs+dSalzzxMt7ISQhIHNHmcLxxvr7/hJ8xrf3xjXv/jm6Pp9W8XyWFkZCRFRUXuvxcV\nFREZGenz882cOZOZM2e6/753byu2x3lIdHR0u4zLEHja42tvL/gfsnk96qrfUhISBn2TKE9dR1Ur\nxCnZO6C6ioqefRpdT7rqz4Li1PUoKyTg8QSa9vj6BxJ75TKIiqHmmt9RXFEJFZX6573ikeWLKczc\ngQpvvMXTaWTLRgDKOkdwsJnXQJIHg1Ls+3o+VpfujsZwvL3+hp8wr/3xjXn9j2/aw+sfFxfn0XHt\noq00OTmZvLw8CgoKqK2tZfny5YwdO7atwzIYWgXZuwf7zeeRyoq2WT83C/nwvzBqImrSDABU4kDI\nzkBqawO/fr0YTVMVGtccouRmBTwWg/NITib0TUJZDb9uVEKy/kN2680dSp5L6Kxn72aPUxHdIXmw\nmTs0GAwGw3FHu0gOg4KCuPrqq3nggQe49dZbmTRpEvHx8SxcuJCFCxcCUFJSwg033MBnn33GvHnz\nuOGGGygvN/L2hqMfmf8+8vVnbaLIKbU12C8/rmfBrrgJpZR+ICkFqquhNRKyzHToFA5NtIyqTuHQ\nPRpyswMfi8FRpKoSCnaj+iQe+WBCf31Ma4rS7MmF7tGosI4tHqqGjIKcTP1vMBgMBoPhOKFdtJUC\njBkzhjFjxjT42axZs9x/7tatG88991xrh2UwBBQpP4is+AY6hCKLPkImTkfFN3IjHaj1v/gAsndg\n3XwnKqKb++cqaSACyI60nyo8gYphZzr0G/BTYtoYcfHIblM5POrYnQ0ijV7TqnO4FqXJaj1RGsnP\nhV7NVw3rUXHxCEB+LgT4PWAwGAwGQ3uhXVQODYbjFfl+MVRXYd18lzZ8f/0ZxLZbb/0138PA4ahR\nExs+ENkDunaHHWmBXb+qCnZncai/YWOo3gmQl2OMyX1ERLCXfIH97kvY/36Eur//ibo7f0XdnEuQ\n1csDt+6unfoPffo1+rjqmwytVDkUEcjPQfXq49kJvbSiqeTtCmBUBoOzSH4uUrq/rcMwGAxHMSY5\nNBjaCBHRZuCJKagho1AXXg070pClC1tn/bo6yMtBudr7DkUpBYkDkQAnh+zaDraNSjwyhgbEJUBN\nNRTuCWw8xyqZGch/n0G+/UK3cSoLlZgCgKz7IXDr5uyEsI4Q3bPxxxP6w949yMHSwMVQz/59UFEO\nsR4mhz1jwbKgfk7RYGjnSG0t9kO3Y99zE7J+VcsnGAwGQyOY5NBgaCu2boD8XNT0MwC0GMzA4cgH\nryIH9gV+/cI8qK2B3n0bfVglpUDBbqTsQMBCqBejoV9Ks8cplyhNq8xAHoNIpvastO5/hqAHniPo\n9w9iXXc7pAxFAuhnKbsyoU+/JluG3S3LrVE9zNdJnmpBjKYeFRwCMbFInpl1NRwlZG6D8oOAwn76\nL9hvv6C9PQ0Gg8ELTHJoMLQR9jfzIbwLatyJgK7WWZfdCFVVyHuvBD4Al8CL6p3Q6MMqaaD+w870\nwMWwMwO6RaG6tWBd4zItN3OHPpKZAV26QmRDjyWVmAJ7cpGDZY4vKSKQm9m4GE09ruSwNeYOxZUc\n4mlbKejW0lauHIoI9guPYn/2bquuazj6ka0bQCmsPz+BOvks5KtPsB/+vZ61bYt47Dokda3uUjEY\nDEcNJjk0GNoAKd4L61agppyKCung/rmK7YM67XxkxTfIlvWBjSE3C5Ryz1YdQUJ/UFZAW0slcxu0\nMG8IoELDoEcvo1jqI5LZuOhPfWspgageFhXoNs74fk0eojp3geierSNKk58LoWHQPcrjU1RcPBTm\nIbWtWH3JTEdWLkE+fB174Yett67BMewfvqXusbmtq8QLyJYNEJ+E6haFdcmvsH59NxQXYv/11lZv\nM5XaWuTFx7D/+Wfkq49bdW2DweAfJjk0GNoA+W6hVnE86bQjHlNnXAg9emG/8VxAW4Jkdxb06IUK\nDW30cRXWEXonBCw5lIN
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.figure.Figure at 0x7f00ff509668>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"plt.subplot(211)\n",
|
||
|
|
"plt.plot(tr,rfq)\n",
|
||
|
|
"plt.ylabel(\"Q-RF\")\n",
|
||
|
|
"plt.subplot(212)\n",
|
||
|
|
"plt.plot(tr,rfl)\n",
|
||
|
|
"plt.xlabel(\"Time\")\n",
|
||
|
|
"plt.ylabel(\"L-RF\")\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Deconvolution in the time domain using Wiener filter\n",
|
||
|
|
"Now, we do the deconvolution in the time domain using the Wiener filter method. We need the cross-correlation and the autocorrelation to form the Toeplitz matrix and the right hand side of the equation system."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 28,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"length of z: 201\n",
|
||
|
|
"length of czr: 401\n",
|
||
|
|
"Center of correlation: 200\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"p = len(zcomp) # length of z and r\n",
|
||
|
|
"czr = np.correlate(rcomp,zcomp,mode = \"full\") # correlation of z with r according to our definition\n",
|
||
|
|
"azz = np.correlate(zcomp,zcomp,mode = \"full\") # auto-correlation of z\n",
|
||
|
|
"q = len(czr) # length of correlations = 2*p-1\n",
|
||
|
|
"jc = q//2 # center of autocorrelation sits at j=q//2 = p-1\n",
|
||
|
|
"print(\"length of z: \",p)\n",
|
||
|
|
"print(\"length of czr: \",q)\n",
|
||
|
|
"print(\"Center of correlation: \",jc)"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 29,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA4cAAAEUCAYAAACYgENSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8VFX6x/HPuekFUulVQJoKiCBNBIW1oLurK2LFrquu\nfVFRdy2LiqgIuj8UEUXXsnbW3hBBBRQEsYCiNOmkkd5zz++Pm0zaJBkQMgG+79eLF5k79945c2Yy\nmeee5zzHWGstIiIiIiIiclBzgt0AERERERERCT4FhyIiIiIiIqLgUERERERERBQcioiIiIiICAoO\nRUREREREBAWHIiIiIiIiAoQGuwGN4fHHH2fFihXExcUxderUevdNTU3liSeeIDs7m9jYWK699lqS\nkpIaqaUiIiIiIiLBcVCMHI4cOZLbb789oH2ff/55jj32WB5++GHGjh3LSy+9tI9bJyIiIiIiEnwH\nxchh7969SUlJqbZtx44dPP3002RnZxMREcFf//pX2rVrx5YtW7jgggsAOOyww3jooYeC0WQRERER\nEZFGdVCMHPoza9YsLrnkEqZMmcL48eOZPXs2AJ06dWLp0qUALF26lIKCAnJycoLZVBERERERkX3u\noBg5rKmwsJA1a9bwyCOP+LaVlpYCMH78eJ555hkWLFhAr169SExMxHEO2hhaREREREQOEgdlcOi6\nLjExMX5TRhMTE5kwYQLgBZFff/01MTExjd1EERERERGRRnVQDolFR0fTsmVLlixZAoC1lo0bNwKQ\nnZ2N67oAzJ07l+OOOy5YzRQREREREWk0xlprg92IfW369OmsXr2anJwc4uLiGDduHIcffjhPPfUU\nmZmZlJaWMmzYMMaOHctXX33FSy+9hDGGXr16cemllxIWFhbspyAiIiIiIrJPHRTBoYiIiIiIiNTv\noEwrFRERERERkeoUHIqIiIiIiMjBUa1027ZtwW5CLcnJyaSlpQW7GQcl9X3wqO+DR30fPOr74FHf\nB4/6PnjU98HVVPu/bdu2Ae2nkUMRERERERFRcCgiIiIiIiIHSVqpiIhIMLjz3iYnPwf+dF6wmyIi\nItKg/So4TEtLY8aMGWRmZmKMYfTo0YwZMybYzRIREfHLfvsVhWk7MAoORURkP7BfBYchISGMHz+e\nLl26UFBQwMSJE+nTpw/t27cPdtNERERqy8zAzczAKSvDhIQEuzUiIiL12q/mHCYkJNClSxcAoqKi\naNeuHRkZGUFulYiISG3WWshMB9eFnMxgN0dERKRB+9XIYVUpKSls2LCBbt261bpv3rx5zJs3D4AH\nHniA5OTkxm5eg0JDQ5tkuw4G6vvgUd8Hj/q+8bl5uaQWFwEQb8sIU/83Or3vg0d9Hzzq++Da3/t/\nvwwOCwsLmTp1KhdddBHR0dG17h89ejSjR4/23W6Ka4001TVQDgbq++BR3weP+r7x2W2bfD9n/rYe\nk9AyiK05OOl9Hzzq++BR3wdXU+3/A3adw9LSUqZOncrw4cMZNGhQsJsjIiLiX2bltAe7S1MgRESk\n6duvgkNrLTNnzqRdu3aceuqpwW6OiIhInWxmeuWNTAWHIiLS9O1XaaVr1qzh888/p2PHjtx8880A\nnHPOOfTv3z/ILRMREamhPCA0sc29wjQiIiJN3H4VHPbs2ZNXX3012M0QERFpWGY6RMcQ0rodpRo5\nFBGR/cB+lVYqIiKyv7CZGRCfREhSC40ciojIfkHBoYiIyL6QmQHxiTgJyZpzKCIi+wUFhyIiIvtC\nZgamYuQwPxdbvuahiIhIU6XgUEREZC+zrgtZ5SOHieWLIWv0UEREmjgFhyIiIntbTha4bnlw2MLb\npnmHIiLSxCk4FBER2dsqlrGoSCulvECNiIhIE6bgUEREZG+rGCWsllaqkUMREWnaFByKiIjsZb5R\nwvgkTHQshEdozqGIiDR5Cg5FRET2tsx0MA40j8cYA/GJCg5FRKTJU3AoIiKyt+VkQWwzTEiId7t5\nPDY7M7htEhERaYCCQxERkb2tqBAioypvR0aB1jkUEZEmTsGhiIhIAGxOFra0NLB9CwshIrJyQ0Qk\nFBYEdqy1WBWvERGRIFBwKCIi0gBbVob7j6uwn/wvsAOKCqoFhyYiyhtNDMSKJbgTL8NmpO5BS0VE\nRPacgkMREZGGpO+E/Fzsr6sD27+oECKqpJVGRAYcHNpfV0FZGWzdtAcNFRER2XMKDkVERBqyc7v3\n/+b1ge1fWACRVdJKI3cjrXTzBu//lG2700IREZHfTcGhiIhIA3yBWmZGYFVHiwox1eYcRkFZKba0\npP7HsRbKg0N2KjgUEZHGpeBQRESkIVUDtU0BjB4WFdROK4WGU0vTdkJBHqCRQxERaXwKDkVERBpg\nU7dDcivv50BSS/1VK63YXp+KwLNFa0jZvgctFRER2XMKDkVERBqycxum86GQ1LIy7bMOtrQUSktq\nzDksH0Usqn/eod28HhwHc+RgSEtpMA1VRERkb1JwKCIiUg9bWgrpKdCyLXTogm0orbS4fHSwSlqp\nCTCt1G5aD63bQ7vOYF0vzVRERKSRKDgUERGpT9pOcF1o1QbTsQukbMPWV3m0InW0ZkEaaLhi6eb1\nmI5dMK3aereVWioiIo1IwaGIiEh9ygvDmJZtveDQWthST2ppkZ/gMLLhkUObnQmZGdDhEGjZxtum\nojQiItKIFByKiIjUw1aM3rVqC+07e9vqW6C+fF6hiYyu3FYeKNr60kq3eec07Q+B2OYQFVO5vqKI\niEgjCA12A0RERJq0lG1eoBbb3Cs0A5CXU/f+FQFgpJ+00voK0lScs3k8xhho2UYjhyIi0qg0cigi\nIlIPu3M7tGyDMQYTFg6hYZCfV/cBfuccNryUha04Z3QMgDfvcKeCQxERaTwKDkVEROqTuh1TPgcQ\n8IK3grqDQ1sxOlilWimBVCutERzSsi1kpGo5CxERaTQKDkVEROpgrfWKxCQmV26MjoG83LoP8lOQ\nxoSEQFh4/Wml+XlgnMqgMiHJK36Tnfk7noGIiEjgGiU4vPbaa1m4cGGd91944YWN0QwREZHdU5AP\nJcXQPL5yW1QMtp6RQ18AWHXOIXjBYr0jh7kQHePNNwRMxWNmKTgUEZHG0SjBYVpaGi+++CKzZs2i\ntLS01v3W2sZohoiIyO6pGLVrnlC5LSY2sDmH4X6Cw3rmHJKfV5lSChCXUL0NIiIi+1ijBIfh4eFM\nmTKFzZs3c8cdd5Camlrt/oqrpCIiIo3FumUNVwPN3gVUGcUDTFRM/cFhUSGEhmFCaxQEj4yqnI/o\nrz0FeRAdW7mh/DFteRvqPC49FVtUVO8+IiIigWi0OYcJCQncc8899O7dm1tvvZUVK1Y01kOLiMgB\nzhbkY79bivvuy9hfVwd2zPx3cf95NTZtZ907VYzaxVUZOWygIA1FBbVTSiHgtFKfioC0npFDW1KM\n+6/rsK8/U/d5q+6/fTPuWy9hv/mywaBTREQOPo26zqHjOFx44YX06NGDf//735xwwgmcddZZSisV\nEZHfxZ31IPzoXXS0nZcRcsfUBo+xi+aD62K/X4Y5/lT/+1TM96s65zDaGzm01vrPfCksrF6ptEKD\nwWEexCX6bpqwcG99xax6grg1P3ptWfo5dtxlmLCwuvcF7PuvYb9agAVo0RrnvieVvSMiIj6NGhxW\nGDx4MB07duSRRx7hl19+wXXdYDRDREQOALaoEH7+HjP8BIiKxn7yNjY/D1N1FK7mMVs3wZYN3s/f\nL4M6gkOyd4HjQEyzym1RsVBWCsVF1dcyrNoeP9uJiKp//mCBnzbHxdc/cvj9Mu+H/DxYtRz6Da57\nX2uxP38P/YdgOnXDzn0edmyBNh3qbpOIiBxUGiWt1N/IYNu2bbn//vtJSkqipERrOImIyB765Uco\nLcUMGIY5/CiwLjSQWmqXfg7GwQweCWt+wBbm+98xOxOaxWOcKn8uY8oDuLrmHRYV+A0OTWQAI4c1\ng8Pm8XWmf1prveDw8KO
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.figure.Figure at 0x7f00ff3f0ba8>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"tm = np.arange(0,q)*dt-jc*dt # shift times for plotting\n",
|
||
|
|
"plt.subplot(211)\n",
|
||
|
|
"plt.xlabel(\"Time\")\n",
|
||
|
|
"plt.ylabel(\"AZZ\")\n",
|
||
|
|
"plt.plot(tm,azz)\n",
|
||
|
|
"plt.subplot(212)\n",
|
||
|
|
"plt.plot(tm,czr)\n",
|
||
|
|
"plt.xlabel(\"Time\")\n",
|
||
|
|
"plt.ylabel(\"CZR\")\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 30,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [],
|
||
|
|
"source": [
|
||
|
|
"m = 150 # length of receiver function with m <= p\n",
|
||
|
|
"toeplmat = np.zeros((m,m))\n",
|
||
|
|
"for i in range(0,m):\n",
|
||
|
|
" for j in range(0,m): \n",
|
||
|
|
" toeplmat[i,j] = azz[jc+i-j] # negative values of k-j are OK because we used \"full\"\n",
|
||
|
|
" if j == i:\n",
|
||
|
|
" toeplmat[i,j] = toeplmat[i,j]+0.05*azz[jc] # regularization of solution\n",
|
||
|
|
"if m <= 4: \n",
|
||
|
|
" print(toeplmat)\n",
|
||
|
|
"rfq = -np.linalg.solve(toeplmat,czr[jc:jc+m]) # solve equation system for receiver function\n",
|
||
|
|
"rfl = np.linalg.solve(toeplmat,azz[jc:jc+m]) # solve for deconvolution of z with itself\n",
|
||
|
|
"rfl = lowpass(rfl,0.75,3,zerophase = True) # low-pass filter result\n",
|
||
|
|
"rfq = lowpass(rfq,0.75,3,zerophase = True)\n",
|
||
|
|
"norm = np.max(rfl)\n",
|
||
|
|
"rfq = rfq/norm # normalize result\n",
|
||
|
|
"rfl = rfl/norm"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 31,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA40AAAENCAYAAACre4DIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXdYVVfWh999aIpdsKOCIopYscWWpjHNSZyUSe+mmjH9\nSy+m90zKZJKZGNPG9ExMj9EUY40dLDSxd8GCIAj3rO+PDcRCueWce1H3+zx5ItxzDj+45967115r\n/ZYSEcFgMBgMBoPBYDAYDIYqsEItwGAwGAwGg8FgMBgMdRcTNBoMBoPBYDAYDAaDoVpM0GgwGAwG\ng8FgMBgMhmoxQaPBYDAYDAaDwWAwGKrFBI0Gg8FgMBgMBoPBYKgWEzQaDAaDwWAwGAwGg6FaTNBo\nMBgMBoPBYDAYDIZqMUGjwWAwGAwGg8FgMBiqxQSNBoPBYDAYDAaDwWCoFhM0GgwGg8FgMBgMBoOh\nWsJDLSCUbNq0KdQSDiM2NpYdO3aEWobhKMbcY4ZgYO4zg9uYe8wQDMx9ZnCbUN9jbdu29eo4k2k0\nGAwGg8FgMBgMBkO1mKDRYDAYDAaDwWAwGAzVYoLGIwD752+w33gGEQm1FIPBYDAYDAaDwXCMYYLG\nOo7YNvLDF8jCWbBkXqjlGAwGg8FgMBgMhmMMEzTWdXJWwM4dEB6OPeW/iO0JtSKDwWAwGAwGg8Fw\nDGGCxjqOzJsBkVGoi2+AjWuR+TNDLclgMBgMBoPBYDAcQ5igsQ4jZaXIwlmoPsehho6EuHjkq8lI\nWVmopRkMBoPBYDAYDIZjBBM01mVWLIHCAtSg41GWhXX2JbBtMzLn51ArMxgMBoPBYDAYDMcIJmis\nw8i8GdCgEXTvo7/ReyAkJCHffISUloZWnMFgqNOIiOmBNhgMBoPB4AgmaKyjSEkxsmQuqt9QVHgE\nAEoprDGXQv4OZMYPIVZoMBjqMvLxW9gTbkGKi0ItxWAwGAwGwxGOCRrrKLJkHuwvQQ06/uAHkntD\n157It58gJcWhEWcwGOo0snUT8vO3sGkd8uk7oZZjMBgMBoPhCMcEjXUU+WMGNIuFxO4Hfb8y21iw\nWy8KDQaD4RDkm48gIhw1dAQy4wdk+eJQSzIYDAaDwXAEY4LGOojs3QPLF6EGDkdZhz9FKjEZevZH\nfvgcKSoMgUKDwVBXkc0bkHkzUCeeibrkRmjTHvvdV5GivaGWZjAYDAaD4QjFBI11EFk0Gzwe1MAT\nqj3GGnMJFO1FfpoSRGUGg6Guo7OMEahT/4qKiMS66lbYnY98PDHU0gwGg8FgMByhmKCxDiLzZkDr\nOGifUO0xqkNn6DcEmTYFKTNOqgaDAWTTOmT+76iTR6MaNwVAJXRBnXYuMns6snR+iBUaDAaDwWA4\nEjFBYx3Ds2MrZC/XsxmVqvFY1XcwFO+DLRuCpM5gODawf/0ezz+fwP5pCrJ+NWLboZbkFfL1RxBZ\nDzXqrwd9X42+ENp1xH7/NaSwIETqDAaDwWAwHKmEh1qA4WCKZ04HEdTA42s9VsUlIIBsWIuKqz4r\neaxj//A5REShTj6z1kDcYACQqf+DnXnIknkIQMNG0LUnqlsvVM8BqJgWoZZ4GLJhDbJwFur081CN\nGh/0mIqIwLr6Vuwn70Qm/xt17R0hUmkwGAwGg+FIxASNdYzi36dCQhKqZdvaD27VFsLDYcMa13Ud\nqcieXcj/3gfbhg2r4ZIbUeHmtjdUj+TvgO1bUBdcg0odimSmQ0YakrEUWTgbafQh1gvv1bkNCPvr\njyCqHmrUmCofVx06o874G/L1h0i/wajUIUFWaDAYDAaD4UilzqyelyxZwqRJk7BtmxEjRjBmzMEL\nHxFh0qRJLF68mKioKG666SY6deoEwLhx46hXrx6WZREWFsbTTz8dil8hYGTzBjy5WagLxnp1vAoP\nh9btkY1r3BV2BCMLZ4FtowafjMz8CcnfjnX93ajoBqGWZqijSPZyAFRSD1TzWNTgk2DwSYgI8vM3\nyEf/gW2b9aZNHUHWr4ZFs1GjL0A1aFTtceqM85Glf2B/8C+sXgPNBorBYDAYDAavqBM9jbZtM3Hi\nRO677z5eeuklZs2axYYNB/fpLV68mC1btvDKK69w3XXX8dZbbx30+MMPP8xzzz13xAaMUD6b0bJQ\n/Yd5fY6KizeZxhqQeb9Bu45YV9+KuuLvkJmO/ew9SN72UEsz1FWylkP9aIiLP+jbSilU1x4AyJrs\nEAirHvvrD6F+A9TIs2s8ToWHY51xHhTshrU5QVJnMBgMBoPhSKdOBI05OTm0bt2aVq1aER4ezpAh\nQ5g//2CXvwULFnD88docJikpicLCQnbu3Bkixe4g6QuI7JGKatrc+5Pi4mFXvp7tWAeQsjKktG64\nucqOrbAqo7I/1Bp2Ctb4hyF/O/ZTdyFrV4VYoaEmxONBduYF/+dmL4fE7igr7PAH23SAyChYnRV0\nXdUh61bB4rmokWehGjSs/YQuKfq8rOUuKzMYDAaDwXC0UCeCxvz8fGJiYiq/jomJIT8//7BjYmNj\nqz3mscce4+6772batGnuC3YJ6+6naXzzfT6do9p11P/YuNYFRd4j+Tuw//c+9l1XYj93b0i1VCDz\nfwdADRhe+T3VvQ/W3c9AWBj2c/ciyxaFSp6hBiRnBfbjt2Hfey2Sty14P3fPLti8HlUeWB2KCguD\njp3rVKZRvvsMohugRp7l1fGqcVNo0x7JWuayMoPBYDAYDEcLR0VDy2OPPUbz5s3ZvXs3jz/+OG3b\ntqV79+6HHTdt2rTKoPLpp58+KAitK4SHhxNbVub18Z5eqewAGuzcTnSQfx8RoXTlUoq+/YySub+B\n2IS1jsOzOotmeAiLbRVUPYeSt3AWqltPmif3OPiB2Fg8z01k5yO3wEf/JuZfn9Y5UxM3CQ8Pr5P3\nPoC9K5+C916n+JfvsJo2x/aU0WBNFtFdD389u0Fxdjq7gaYDhxJZzd+ooHtvir77nJimTUPeEyi2\nzfbMNOoNOZkmHTp6fd6eXv0pnvEjMc2aosLc+R3q8n1mODow95ghGJj7zOA2R8o9VieCxubNm5OX\n92cZWl5eHs2bNz/smB07dlR5TMX/mzRpwoABA8jJyakyaBw5ciQjR46s/PrA69UVYmNjfdIlItCw\nMXszl1N0XPB+H1mxBPuzSbB+dWWWQ514OrK/BB75O3m/T8caPipoeg7Tt2EN9tpVqIuuq/bvaR9/\nKvLfN9ixPA3Vul2QFYYOX++xYCC2B/ntR+TL96GkBHX6uXDmBfDQOAr++J2ift73+QaCvXAuREax\nu0kMqrr7plUclO5nx9KFqI6dg6KrOmTzBmRvASVtO/r0nNodOiP7itixeD4qvosr2urifWY4ugjV\nPWb//A0U7kWdcCqqcbOg/3xDcDHvZQa3CfU91ratd8Z+daI8tXPnzmzevJlt27ZRVlbG7Nmz6d+/\n/0HH9O/fnxkzZiAiZGVlER0dTbNmzSguLmbfvn0AFBcXk5aWRocOHULxa4QEpRTExSNBLE8VEexJ\nL0NRIeqym7CenYR1/lWoFq2hbQdoGoMsD23Z55+mQkOrPUalpOpjly8OlixDFcj2LdhP3oVMfgM6\ndMZ6+BWsc65ARdVDpfSFlWmIxxMcLVnLoHM3VHhEtcdUBFl1oURVcjMAUJ27+XSeSqroazQlqgaD\nL0hJMfLpJOSrydh3X4P9zsuIMaMzGAzHAHUi0xgWFsbVV1/NE088gW3bnHTSSbRv356pU6cCMGrU\nKPr27cuiRYsYP348kZGR3HTTTQDs3r2b559/HgCPx8OwYcPo06dPyH6XUKDadUR+n4rYnqrNO5xm\n6ybYlYe69Cas4087WItSqJS+yKI5iMeje8CCjIjooDG5d427wKpFa2jZRge4I0YHUaHhQOT7z3Qf\n4bV3ogYMP6hUWHXvg/w+FdZkg4+Bkc86CvfChjWov1xU84GxraBhY22Gc8JpNR/rNqsyILoBtI7z\n6TTVNEbf+1nLYdRfXRJ
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.figure.Figure at 0x7f00ff450240>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"tm = np.arange(0,m)*dt # shift times for plotting\n",
|
||
|
|
"plt.subplot(211)\n",
|
||
|
|
"plt.plot(tm,rfq)\n",
|
||
|
|
"plt.ylabel(\"Q-RF\")\n",
|
||
|
|
"plt.subplot(212)\n",
|
||
|
|
"plt.plot(tm,rfl)\n",
|
||
|
|
"plt.xlabel(\"Time\")\n",
|
||
|
|
"plt.ylabel(\"L-RF\")\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "markdown",
|
||
|
|
"metadata": {},
|
||
|
|
"source": [
|
||
|
|
"## Iterative deconvolution\n",
|
||
|
|
"Finally, we try iterative deconvolution where the time of the conversions is determined experimentally by searching for maxima of the cross-correlation between radial and vertical component."
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 50,
|
||
|
|
"metadata": {
|
||
|
|
"scrolled": false
|
||
|
|
},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"name": "stdout",
|
||
|
|
"output_type": "stream",
|
||
|
|
"text": [
|
||
|
|
"Iteration: 0 Location of max correlation: 45 Coefficient: -0.188064932852 RMS: 0.808067884685\n",
|
||
|
|
"Iteration: 10 Location of max correlation: 56 Coefficient: 0.0232878399425 RMS: 0.236482738022\n",
|
||
|
|
"Iteration: 20 Location of max correlation: 26 Coefficient: -0.0132327288213 RMS: 0.157940698442\n",
|
||
|
|
"Iteration: 30 Location of max correlation: 55 Coefficient: 0.00682842269855 RMS: 0.127478308995\n",
|
||
|
|
"Iteration: 40 Location of max correlation: 96 Coefficient: -0.00547868539688 RMS: 0.113284892005\n",
|
||
|
|
"Iteration: 50 Location of max correlation: 117 Coefficient: -0.00491286692813 RMS: 0.0989209837833\n",
|
||
|
|
"Iteration: 60 Location of max correlation: 45 Coefficient: 0.00363439542244 RMS: 0.0911542901353\n",
|
||
|
|
"Iteration: 70 Location of max correlation: 56 Coefficient: 0.00285406937047 RMS: 0.084300593319\n",
|
||
|
|
"Iteration: 80 Location of max correlation: 44 Coefficient: 0.00276115316265 RMS: 0.0791800623638\n",
|
||
|
|
"Iteration: 90 Location of max correlation: 45 Coefficient: 0.0022053005175 RMS: 0.0752592825827\n",
|
||
|
|
"Iteration: 100 Location of max correlation: 37 Coefficient: -0.00211030980769 RMS: 0.0713918083254\n",
|
||
|
|
"Iteration: 110 Location of max correlation: 92 Coefficient: 0.0022716970582 RMS: 0.0684080863292\n",
|
||
|
|
"Iteration: 120 Location of max correlation: 53 Coefficient: -0.00184308153988 RMS: 0.0657632835901\n",
|
||
|
|
"Iteration: 130 Location of max correlation: 124 Coefficient: 0.00138745498748 RMS: 0.0637757664872\n",
|
||
|
|
"Iteration: 140 Location of max correlation: 30 Coefficient: 0.00138958532611 RMS: 0.0622412721775\n",
|
||
|
|
"Iteration: 150 Location of max correlation: 45 Coefficient: 0.00145676616742 RMS: 0.0602417928777\n",
|
||
|
|
"Iteration: 160 Location of max correlation: 14 Coefficient: 0.0010268119389 RMS: 0.0586738835766\n",
|
||
|
|
"Iteration: 170 Location of max correlation: 12 Coefficient: -0.00139735070189 RMS: 0.0571199966087\n",
|
||
|
|
"Iteration: 180 Location of max correlation: 4 Coefficient: 0.00151602625314 RMS: 0.0553731496016\n",
|
||
|
|
"Iteration: 190 Location of max correlation: 45 Coefficient: 0.00130031181716 RMS: 0.0540243403414\n",
|
||
|
|
"Iteration: 200 Location of max correlation: 4 Coefficient: 0.00103009869364 RMS: 0.0528147877942\n",
|
||
|
|
"Iteration: 210 Location of max correlation: 37 Coefficient: -0.000981396173385 RMS: 0.0516961431046\n",
|
||
|
|
"Iteration: 220 Location of max correlation: 22 Coefficient: -0.000948845430963 RMS: 0.0504874410121\n"
|
||
|
|
]
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"itmax = 1000 # stop after itmax iterations\n",
|
||
|
|
"eps = 0.1 # or if coeff less than eps\n",
|
||
|
|
"zcomp = st[2] # rename traces\n",
|
||
|
|
"rcomp = st[1]\n",
|
||
|
|
"p = len(zcomp) # length of z and r\n",
|
||
|
|
"it = 0 # start with iteration 0\n",
|
||
|
|
"rnew = rcomp[:] # copy of radial component\n",
|
||
|
|
"norm = np.dot(rnew,rnew) # initial energy in radial component\n",
|
||
|
|
"zsh = np.zeros(p) # shifted zcomp allocated here\n",
|
||
|
|
"a = np.zeros(itmax) # space allocated for coefficients a\n",
|
||
|
|
"jmax = np.empty(itmax,dtype = int) # space allocated for conversion samples\n",
|
||
|
|
"while it < itmax: # iteration loop\n",
|
||
|
|
" czr = np.correlate(rnew,zcomp,mode = \"full\") # correlation of z with r according to our definition\n",
|
||
|
|
" jc = len(czr)//2 # positive lags start at j=q//2 = p-1\n",
|
||
|
|
" jmax[it] = np.argmax(np.abs(czr[jc:])) # location of absolute maximum of czr\n",
|
||
|
|
" zsh[:] = 0.0 # zero shifted z(t)\n",
|
||
|
|
" zsh[jmax[it]:] = zcomp[0:p-jmax[it]] # shifted z(t-tmax) cut at end\n",
|
||
|
|
" a[it] = np.dot(rnew,zsh)/np.dot(zsh,zsh) # compute coeficient\n",
|
||
|
|
" rnew = rnew-a[it]*zsh # subtract shifted z(t-tmax) from radial component\n",
|
||
|
|
" mf = np.sqrt(np.dot(rnew,rnew)/norm) # normalized residual amplitude in radial component\n",
|
||
|
|
" if it%10 == 0:\n",
|
||
|
|
" print(\"Iteration: \",it,\" Location of max correlation: \",jmax[it],\n",
|
||
|
|
" \" Coefficient: \",a[it],\" RMS: \",mf)\n",
|
||
|
|
" if abs(mf) < eps: break\n",
|
||
|
|
" it = it+1"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 51,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA5AAAAENCAYAAAB5K4qQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXl4VNXdxz/nZrKThCSTPYEAYZN9UQQFRGmte4v70qrV\nVmvdu1r3Wvu64S762lrr7mutorWuiCyCKKgoiwIhARKyzyQh+zL3vH+cmSEhM5k7CwTlfJ4nTyb3\nnjnnzL03d+7v/JavkFJKNBqNRqPRaDQajUajCYAx0BPQaDQajUaj0Wg0Gs13A21AajQajUaj0Wg0\nGo3GEtqA1Gg0Go1Go9FoNBqNJbQBqdFoNBqNRqPRaDQaS2gDUqPRaDQajUaj0Wg0ltAGpEaj0Wg0\nGo1Go9FoLKENSI1Go9FoNBqNRqPRWEIbkBqNRqPRaDQajUajsYQ2IDUajUaj0Wg0Go1GYwltQGo0\nGo1Go9FoNBqNxhK2gZ7AwUJFRcVAT6EPdrudurq6gZ7GIYk+9gOHPvYDhz72A4s+/gOHPvYDhz72\nA4c+9gPHwXrsc3NzLbXTHkiNRqPRaDQajUaj0VhCG5AajUaj0Wg0Go1Go7GENiA1Go1Go9EEjfx8\nNa77b0a2tw30VDQajUZzANEGpEaj0Wg0mqCRX30G33yFfPnJgZ6KRqPRaA4g2oDUaDQajUYTNLKm\nAgwDuepDzLUrB3o6Go1GozlAaANSo9FoNBpN8FTvRsycB8NHI59bhHTUDPSMNBqNRnMA0AakRqPR\naA4o0nQhXa6BnoYmDGRLEzQ3Qe4QjEt/A9LE/Pv9SFOfV41Go/m+ow1IjUaj0RxQ5Iv/i3nfnwZ6\nGppwqFbaySIrD5GRjTj/cijejHzn3wM8MY1Go9Hsb7QBqdFoNJoDhjRdyM9XwfZvkW2tAz0dTYhI\ntwFJphKdNo6cB2MmIj9bMYCz0mg0Gs2BQBuQGo1Gozlw7CxRoY9Swo5tAz0bTahU7wZhQEaWd5PI\nHQL1dQM4KY1Go9EcCGwDPQGNRqPRDDzS5aJ9zXLk8DEII2r/jbPpi72vS7Ygxk7ab2Np9iPVFWDP\nRNii925Ls0NbK7KtFRGfENHh5FefIct3QFw8xMYhEgbB+KmImNiIjqPRaDSawGgDUqPRaDSwfg2N\nT9yNuOhqxFHz99swctOXMLQI2tuQpVv32zia/YusqYCsvN4bU+3qt7MO8oZEbiwpMf/xILQ2790G\niHN/iTj25IiNo9FoNBpr6BBWjUaj0SB3Fqvfb/9rv1VIla0tUPItYtxUxPBRULIFKeV+GUuz/5BS\nQnUFIiu313aRlqFe1NdGdsCWJmhtRpx+IcYDz2Pc9XcYlATlOyI7jkaj0WgsoQ1IjUaj0SB3lUB0\nDNRUItd9vH8G+fYrME3EuCkwbDQ0NYLWDvzu0eiEjnbYx4AkTXkgpTPCeZA1lQCInALEoGREeibk\nFCAryyM7jkaj0WgsoQ1IjUajOcSRUsKuEuKOng+5Q5D/fQVpmpEfZ+MXEJ8Aw0crDyToMNbvIl4J\nj30MyJQ0VVgnwoV0pNuAJDPHu03kDIGKXdqDrdFoNAOANiA1Go3mUKfRCU2NRI8YjTjxTKgsgy8/\niegQUkqV/zhmIsJmg7xC5fEs0Qbkdw2vhMc+OZAiKgoGp4EjwiGsNZUgBNj3Vnwlt0DlRDY1RHYs\njUaj0QREG5AajUZzqFNWCoBt2EjE4UdDZi7mW69E1rtTtRuctYhxUwGUETl0BLJ0S+TG0BwYqivA\nFr23aE5P0uzISEt51FZCqh0RHePdJHIK1IuKssiOpdFoNJqAHBRVWOvq6njsscdoaGhACMH8+fM5\n8cQTaW5u5oEHHqC2tpaMjAyuu+46Bg0aBMDrr7/O0qVLMQyDiy++mMmTJwNQUlLCY489RmdnJ1Om\nTOHiiy9GCDGQH0+j0WgOauSuEgBshSMRrW2IE89E/vMh+HotTDoiMmNs+hxA5T+6EcNGIT96G9nd\n1VsOQnNQI6t3Q2YOwui7Bi1S7d7rKWLj1VT2Cl8FlAcSkJVliDETIzqeRqPRaPrnoPBARkVF8dOf\n/pQHHniAO++8k/fee4/y8nIWL17MhAkTePjhh5kwYQKLFy8GoLy8nNWrV3P//fdz44038tRTT2G6\n83X+9re/cdlll/Hwww9TVVXF+vXrB/KjaTQazUGP3FUCmTkYCYkAiBlzIT0T87+vRG6MTV9Cdh6i\nZxjisNHQ3aWrae4npOnC/O8ruO64DtnSHPgNVqmu6FtAx0OaHerrIuu9rqlE7GtApqRBfKL2QGo0\nGs0AcFAYkKmpqQwfPhyA+Ph48vLycDqdrF27lrlz5wIwd+5c1q5dC8DatWuZNWsW0dHRZGZmkp2d\nTXFxMfX19bS1tTFq1CiEEMyZM8f7Ho1Go9H4oawEUTDc+6ew2RBzjofSrciWprC7l50dsGWjN3zV\nO87w0Wq/LqQTcWRtFeY9NyAXPw+7tsOObZHp13RBbRViXw1ID6l26OqE5vCvGwDZ2gzNe/p4IIUQ\nkFuArNQGpEaj0RxoDgoDsic1NTWUlpZSVFREY2MjqampAAwePJjGxkYAnE4n6enp3vekpaXhdDr7\nbE9PT8fpdB7YD6DRaDTfIWRrC9RWQcGwXtu9BkIkZDa2bYauzj4GJGl2SEnVhXQijLn6Q8zbr4GK\nMsTZlwJEztBy1IKr268HUrilPCKmBVlbpfrNyOmzS+QUQMWuyIyj0Wg0GsscFDmQHtrb21m4cCEX\nXXQRCQkJvfYJISKay7hkyRKWLFkCwF133YXd7qMYwABjs9kOynkdCuhjP3DoY39g6dxYRj2QMn5K\nr2PfVTQKJ5DU2U5cmOejubqcFsA+cw4iLr7XvobR4+neWXzInfOOtR9jKywiKiPbuy0S1353RRmO\npx8i+rDJpFx7C4Y9i9r/vkJcfS3JETjGHbuKaQAGjxpLjI/+uoYVqeumuzPs6wag/dv1NAKDRx9G\n9D79tRSNofnjD0iLsWEkDw5rHH3fGTj0sR849LEfOL7rx/6gMSC7u7tZuHAhs2fPZsaMGQCkpKRQ\nX19Pamoq9fX1JCcnA8rj6HA4vO91Op2kpaX12e5wOEhLS/M53vz585k/f77377q6CFeNiwB2u/2g\nnNehgD72A4c+9gcWc+OXAOwZbCe2u9t77GWUqni5p7SY5qJx4Y1RvhOSB+NoboHmlt778gqRn62k\ndmcpIjEprHG+K8iWZsz/+SNkZGHccC9ikPpui8S1Lzeo8+lacCH1wgYOBzI7n7bSbXRG4P/K3PYN\nAI1xiQgf/ckoVQxpz47tNA8fG/5425V3usEW22c8maK+3x0bv0KMCu8a1fedgUMf+4FDH/uB42A9\n9rm5fvLb9+GgCGGVUvLEE0+Ql5fHySef7N0+ffp0li9fDsDy5cs5/PDDvdtXr15NV1cXNTU1VFZW\nUlRURGpqKvHx8WzduhUpJStWrGD69OkD8pk0Go3mO8GuEkhJRaSk9t6emASx8REJYZX1db4lH9ib\nB8mhlAdZ/A1IE2oqMZ+4G9ndFbGuZVW5etEjxFTk5ENFWWQK21RXQHwCJPnx+A1KAZsNnBF6MKqp\nhMFpiNi4vvtyhgARDM/VaDQajSUOCg/kli1bWLFiBUOGDOF3v/sdAOeeey4//vGPeeCBB1i6dKlX\nxgOgoKCAmTNncv3112MYBpdccgmGu5z4pZdeyqJFi+js7GTy5MlMmTLF77gajUZzsGO+vxjqqsGe\nhbBnqmIieYURC+mXZSV98h/BXaTEnomsqw5/kHoH+MhhA6CwCIRAlmxFjJ8W/ljfAeS2jWCzIc67\nHPnso8gX/xd++uvIdF5VDmn23qHCuQWw8n1oaoQwQz1ldQVk5vq9/oRhqMWCCGlBylofEh4e0uxq\nkUMbkIck5mcr4Ms1EBMLsbEQG4+YOQ+RO2Sgp6bRfO8J2oCUUtLU1ERSUlLEHmDGjBnDK6/4Lhd/\nyy23+Ny+YMECFixY0Gf
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.figure.Figure at 0x7f00ff4b94a8>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"tm = np.arange(0,p)*dt\n",
|
||
|
|
"plt.subplot(211)\n",
|
||
|
|
"plt.plot(tm,rcomp)\n",
|
||
|
|
"plt.ylabel(\"Q-Comp\")\n",
|
||
|
|
"plt.subplot(212)\n",
|
||
|
|
"plt.plot(tm,rnew)\n",
|
||
|
|
"plt.xlabel(\"Time\")\n",
|
||
|
|
"plt.ylabel(\"Residual\")\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": 52,
|
||
|
|
"metadata": {},
|
||
|
|
"outputs": [
|
||
|
|
{
|
||
|
|
"data": {
|
||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA40AAAENCAYAAACre4DIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl81NW9P/7X+cxkD9ljQtiEAAIKQowbLkCJ2NaldPO6\ndLlarfenVa9e+1WsaC3SS69aqrVWb0Vate21dqHVosVolVZUsIAoLiyCbEHISpbZP+f3x5nP7JPM\nJDP5fDLzej4ePDIz+czMyclkmPfn/T7vI6SUEkREREREREQxaGYPgIiIiIiIiKyLQSMRERERERHF\nxaCRiIiIiIiI4mLQSERERERERHExaCQiIiIiIqK4GDQSERERERFRXAwaiYiIiIiIKC4GjURERERE\nRBQXg0YiIiIiIiKKi0EjERERERERxWU3ewBmOnTokNlDiFJVVYXW1lazh5GVOPfm4dybh3NvHs69\nuTj/5uHcm4dzbx6rzn1dXV1CxzHTSERERERERHExaCQiIiIiIqK4GDQSERERERFRXAwaiYiIiIiI\nKC4GjURERERERBQXg0YiIiIiIiKKi0EjERERERERxcWgkYiSJtuOQL77L7OHQURERETDgEEjESVN\nvvI89P/9H7OHQURERETDgEEjESXP6QRcLkgpzR4JEREREaUZg0YiSp7XA0gd8HnNHgkRERERpRmD\nRiJKnsft/+oxdxxERERElHYMGokoadIIFo3gkYiIiIgyFoNGIkqe18g0MmgkIiIiynQMGokoecw0\nEhEREWUNBo1ElDyvETRyTSMRERFRpmPQSETJ87A8lYiIiChbMGgkouSxPJWIiIgoazBoJKLkMdNI\nRERElDUYNBJR8rimkYiIiChr2M0egGHr1q1YvXo1dF3HwoULsXjx4rDvHzx4EI888gj27NmDSy+9\nFBdffHHge9dffz3y8/OhaRpsNhtWrFgx3MMnyi7+YFG6XRAmD4WIiIiI0ssSQaOu61i1ahXuvPNO\nVFZWYsmSJWhsbMTYsWMDxxQXF+PKK6/Epk2bYj7G3XffjZKSkuEaMlF2M8pSvcw0EhEREWU6S5Sn\n7tq1C7W1taipqYHdbsfcuXOjgsPS0lJMnjwZNpvNpFESUQAb4RARERFlDUtkGtvb21FZWRm4XllZ\niZ07dyb1GMuWLYOmaTjvvPPQ1NQU85jm5mY0NzcDAFasWIGqqqrBDzpN7Ha7JceVDTj3iZE+L45I\nHQBQlJODohTMGefePJx783DuzcX5Nw/n3jyce/OM9Lm3RNA4VMuWLUNFRQW6urpw7733oq6uDjNm\nzIg6rqmpKSygbG1tHc5hJqSqqsqS48oGnPvESKcjcLm3swOOFMwZ5948nHvzcO7Nxfk3D+fePJx7\n81h17uvq6hI6zhLlqRUVFWhrawtcb2trQ0VFRVL3B1QJ66mnnopdu3alfIxE5BfaMZVrGomIiIgy\nniWCxvr6erS0tODIkSPwer3YsGEDGhsbE7qv0+mEw+EIXN62bRvGjx+fzuESZbfQdYxurmkkIiIi\nynSWKE+12Wy46qqrsHz5cui6jgULFmDcuHFYt24dAGDRokXo7OzE7bffDofDASEE1q5dix//+Mfo\n7u7G/fffDwDw+Xw4++yzMXv2bDN/HKLM5g0JFNkIh4iIiCjjWSJoBICGhgY0NDSE3bZo0aLA5bKy\nMjz66KNR9yssLMR9992X9vERkZ/HG7zM8lQiIiKijGeJ8lQiGkFCM41ul3njICIiIqJhwaCRiJIT\nUpIqPcw0EhEREWU6Bo1ElJzQQJFrGomIiIgyHoNGIkqOsY7RZueaRiIiIqIswKCRiJJjZBcLi7im\nkYiIiCgLMGgkoqQE1jEWFIWXqhIRERFRRmLQSETJCc00ck0jERERUcZj0EhEyTGyiwwaiYiIiLIC\ng0YiSo6/+Y1geSoRERFRVmDQSETJYXkqERERUVZh0EhEyWF5KhEREVFWYdBIRMnxuAF7DpCTB3jc\nkFKaPSIiIiIiSiMGjUSUHK8HyMlR/wDA6zV3PERERESUVgwaiSg5Ho8/05jrv+4ydzxERERElFYM\nGokoOR63ChgDQSM7qBIRERFlMgaNRJQcryciaGQzHCIiIqJMxqCRiJIiPW7Abg+uaWTQSERERJTR\nGDQSUXL8mUaRy0wjERERUTZg0EhEyfH4u6fauaaRiIiIKBswaCSi5EQ1wmGmkYiIiCiTMWgkouQY\nQSPLU4mIiIiyAoNGIkqO18tGOERERERZhEEjESXH44bIyQ2saZRc00hERESU0Rg0ElFyPBH7NLpd\n5o6HiIiIiNKKQSMRJcfjBuw5wTWNXmYaiYiIiDKZ3ewBGLZu3YrVq1dD13UsXLgQixcvDvv+wYMH\n8cgjj2DPnj249NJLcfHFFyd8XyJKIa9brWfkmkYiIiKirGCJTKOu61i1ahXuuOMOrFy5Eq+//joO\nHDgQdkxxcTGuvPJKXHTRRUnfl4hSQ0oJeLxqPaOxT6ObQSMRERFRJrNE0Lhr1y7U1taipqYGdrsd\nc+fOxaZNm8KOKS0txeTJk2Gz2ZK+LxGliM8HSB3IyYGw2wFNU2sciYiIiChjWaI8tb29HZWVlYHr\nlZWV2LlzZ8rv29zcjObmZgDAihUrUFVVNYRRp4fdbrfkuLIB535guqMXRwEUlZWjqKoKR3LzUJBj\nw6ghzhvn3jyce/Nw7s3F+TcP5948nHvzjPS5t0TQOFyamprQ1NQUuN7a2mriaGKrqqqy5LiyQbrm\nXr63GSgohKiflvLHHm6yuwsA0Ot2w9HaCmm3w9HVBdcQ542ve/Nw7s3DuTcX5988nHvzcO7NY9W5\nr6urS+g4S5SnVlRUoK2tLXC9ra0NFRUVab8v0XDQn30C+tpnzR5GahilqHZ/E5ycPDbCISIiIspw\nlgga6+vr0dLSgiNHjsDr9WLDhg1obGxM+32JhoWzD+jrNXsUqeH1B4jGHo05OVzTSERERJThLFGe\narPZcNVVV2H58uXQdR0LFizAuHHjsG7dOgDAokWL0NnZidtvvx0OhwNCCKxduxY//vGPUVhYGPO+\nRJbhdACODAka/QGiMLbbyMmF9LhMHBARERERpZslgkYAaGhoQENDQ9htixYtClwuKyvDo48+mvB9\niaxASgm4nCpwzASeiEyjnZlGIiIiokxnifJUoozl9ahtKjIs0wgj05ibyzWNRERERBmOQSNROhkZ\nRodDZR1HOm9kI5zc4G1ERERElJEYNBKlkxE0Sl2VqY50keWpObmAm2saiYiIiDIZg0aidHKFrGV0\n9Jk3jlSJKE8VOblc00hERESU4Rg0EqVTaAOcDFjXKI1Moz10yw2uaSQiIiLKZAwaLUwe3AfZ12P2\nMGgonBmWafRGNMLJYSMcIiIiokzHoNHC9PuWQL74R7OHQUMRuo4xE4LGyO6pLE8lIiIiyngMGi1K\netxAbzfQ2W72UGgIZEimUWZE0BhZnspMIxEREVGmY9BoVb2qLJXlqSNchq1pjO6emgN4PZC6bt6Y\niIiIiCitGDRalT9oDHylkSk0aHRmQKYxsE+jXX01gkfu1UhERESUsRg0WpWRYWSmcWRzOQCbDRAi\nc9Y02nMghFDXjaCR6xqJiIiIMpbd7AFQHL3d6iuDxpHN6QDyCgApMyNo9HqCgSIQEjS6ABSbMiQi\nIiIiSi8GjRYVWMvI8tSRzekE8gvU5UwIGj3uYOdUIHiZmUYiIiKijMWg0aqMoNHjhvS4IUKzOzRi\nSJdDBY1CQGZKI5ywTGNe8HYiIiIiykhc02hVoRlGZhtHLqcDyMsHCgozI9Po9YZlGkUg08igkYiI\niChTMWi0KgaNmcHlL08tKMqIoFF63IA9tDzVWNPIoJGIiIgoUzFotKrQQJHNcEYupypPFfkFGRE0\nRpensnsqERERUabjmkaLkn09aqsGn49B40jmdEDkFQC5eUBGrGn0RDTCYaaRiIiIKNMx02hVvd1A\nRTUAQBrbb9DIYzTCyZQ
|
||
|
|
"text/plain": [
|
||
|
|
"<matplotlib.figure.Figure at 0x7f00ff900630>"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
"metadata": {},
|
||
|
|
"output_type": "display_data"
|
||
|
|
}
|
||
|
|
],
|
||
|
|
"source": [
|
||
|
|
"rfq = np.zeros(p)\n",
|
||
|
|
"for j,v in zip(jmax[0:it+1],a[0:it+1]):\n",
|
||
|
|
" rfq[j] = rfq[j]-v\n",
|
||
|
|
"#rfq = lowpass(rfq,0.8,5,zerophase = True)\n",
|
||
|
|
"plt.plot(tm,rfq)\n",
|
||
|
|
"plt.xlabel(\"Time\")\n",
|
||
|
|
"plt.ylabel(\"Q-RF\")\n",
|
||
|
|
"plt.show()"
|
||
|
|
]
|
||
|
|
},
|
||
|
|
{
|
||
|
|
"cell_type": "code",
|
||
|
|
"execution_count": null,
|
||
|
|
"metadata": {
|
||
|
|
"collapsed": true
|
||
|
|
},
|
||
|
|
"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
|
||
|
|
}
|