Compare commits
9 Commits
06d2f0c31c
...
develop
| Author | SHA1 | Date | |
|---|---|---|---|
| 138d954abf | |||
| 8ce3df158c | |||
| 274ba54371 | |||
| 1a077cb628 | |||
| 3f22d1c923 | |||
| ea6168f6c2 | |||
| 29f0a4dc76 | |||
| 4310d0acf0 | |||
| 010335a394 |
911
03-Cross_Correlation/3-ambient_noise.ipynb
Normal file
911
03-Cross_Correlation/3-ambient_noise.ipynb
Normal file
@@ -0,0 +1,911 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"<div style='background-image: url(\"../images/header.svg\") ; padding: 0px ; background-size: cover ; border-radius: 5px ; height: 250px'>\n",
|
||||
" <div style=\"float: right ; margin: 50px ; padding: 20px ; background: rgba(255 , 255 , 255 , 0.7) ; width: 50% ; height: 150px\">\n",
|
||||
" <div style=\"position: relative ; top: 50% ; transform: translatey(-50%)\">\n",
|
||||
" <div style=\"font-size: xx-large ; font-weight: 900 ; color: rgba(0 , 0 , 0 , 0.8) ; line-height: 100%\">Ambient Seismic Noise Analysis</div>\n",
|
||||
" <div style=\"font-size: large ; padding-top: 20px ; color: rgba(0 , 0 , 0 , 0.5)\">Cross Correlation </div>\n",
|
||||
" </div>\n",
|
||||
" </div>\n",
|
||||
"</div>\n",
|
||||
"\n",
|
||||
"In this tutorial you will try to reproduce one of the figures in Shapiro _et al._. To see which one, execute the second code block below. \n",
|
||||
"\n",
|
||||
"Reference: *High-Resolution Surface-Wave Tomography from Ambient Seismic Noise*, Nikolai M. Shapiro, et al. **Science** 307, 1615 (2005);\n",
|
||||
"DOI: 10.1126/science.1108339\n",
|
||||
"\n",
|
||||
"##### Authors:\n",
|
||||
"* Celine Hadziioannou\n",
|
||||
"* Ashim Rijal\n",
|
||||
"---\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"### In this notebook\n",
|
||||
"We will reproduce figure B below. This figure compares: \n",
|
||||
"1) the seismogram from an event near station MLAC, recorded at station PHL (top)\n",
|
||||
"2) the \"Greens function\" obtained by correlating noise recorded at stations MLAC and PHL (center and bottom)\n",
|
||||
"\n",
|
||||
"All bandpassed for periods between 5 - 10 seconds. \n",
|
||||
"\n",
|
||||
"<img src=\"https://raw.github.com/ashimrijal/NoiseCorrelation/master/data/shapiro_figure.png\">\n",
|
||||
"\n",
|
||||
"---"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "slide"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Configuration step (Please run it before the code!)\n",
|
||||
"\n",
|
||||
"# MatplotLib.PyPlot\n",
|
||||
"%matplotlib inline\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"plt.style.use('ggplot') # Matplotlib style sheet - nicer plots!\n",
|
||||
"plt.rcParams['figure.figsize'] = 12, 8\n",
|
||||
"\n",
|
||||
"# NumPy\n",
|
||||
"import numpy as np\n",
|
||||
"\n",
|
||||
"# ObsPy\n",
|
||||
"from obspy.core import UTCDateTime, read\n",
|
||||
"from obspy.clients.fdsn import Client\n",
|
||||
"try: # depends on obspy version; this is for v1.1.0\n",
|
||||
" from obspy.geodetics import gps2dist_azimuth as gps2DistAzimuth\n",
|
||||
"except ImportError:\n",
|
||||
" from obspy.core.util import gps2DistAzimuth\n",
|
||||
"\n",
|
||||
"# ignore warnings from filters\n",
|
||||
"import warnings\n",
|
||||
"warnings.filterwarnings('ignore')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "slide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"### 1. Read in noise data \n",
|
||||
"Read the noise data for station MLAC into a stream. \n",
|
||||
"\n",
|
||||
"Then, **read in** noise data for station PHL.\n",
|
||||
"Add this to the stream created above.\n",
|
||||
"\n",
|
||||
"These two data files contain 90 days of vertical component noise for each station.\n",
|
||||
"#### If you need data for more than 90 days, it can be downloaded form IRIS database. ###"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Shapiro et al. use noise data from MLAC and PHL stations\n",
|
||||
"\n",
|
||||
"num_of_days = 90 # no of days of data: change if more than 90days of data is required\n",
|
||||
"if num_of_days <= 90:\n",
|
||||
" # get noise data for station MLAC\n",
|
||||
" stn = read('https://raw.github.com/ashimrijal/NoiseCorrelation/master/data/noise.CI.MLAC.LHZ.2004.294.2005.017.mseed')\n",
|
||||
" # get noise data for the station PHL and add it to the previous stream\n",
|
||||
" stn += read('https://raw.github.com/ashimrijal/NoiseCorrelation/master/data/noise.CI.PHL.LHZ.2004.294.2005.017.mseed')\n",
|
||||
" # if you have data stored locally, comment the stn = and stn += lines above\n",
|
||||
" # then uncomment the following 3 lines and adapt the path: \n",
|
||||
" # stn = obspy.read('./noise.CI.MLAC.LHZ.2004.294.2005.017.mseed')\n",
|
||||
" # stn += obspy.read('noise.CI.PHL.LHZ.2004.294.2005.017.mseed')\n",
|
||||
" # ste = obspy.read('event.CI.PHL.LHZ.1998.196.1998.196.mseed')\n",
|
||||
"else:\n",
|
||||
" # download data from IRIS database\n",
|
||||
" client = Client(\"IRIS\") # client specification\n",
|
||||
" t1 = UTCDateTime(\"2004-10-20T00:00:00.230799Z\") # start UTC date/time\n",
|
||||
" t2 = t1+(num_of_days*86400) # end UTC date/time\n",
|
||||
" stn = client.get_waveforms(network=\"CI\", station=\"MLAC\",location=\"*\", channel=\"*\",\n",
|
||||
" starttime=t1, endtime=t2) # get data for MLAC\n",
|
||||
" stn += client.get_waveforms(network=\"CI\", station=\"PHL\", location=\"*\", channel=\"*\",\n",
|
||||
" starttime=t1, endtime=t2) # get data for PHL and add it to the previous stream"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "slide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"### 2. Preprocess noise ###\n",
|
||||
"***Preprocessing 1***\n",
|
||||
"* Just to be sure to keep a 'clean' original stream, first **copy** the noise stream with [st.copy()](https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.copy.html)\n",
|
||||
"The copied stream is the stream you will use from now on. \n",
|
||||
"\n",
|
||||
"* In order to test the preprocessing without taking too long, it's also useful to first **trim** this copied noise data stream to just one or a few days. This can be done with [st.trim()](https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.trim.html), after defining your start- and endtime. \n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Many processing functionalities are included in Obspy. For example, you can remove any (linear) trends with [st.detrend()](https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.detrend.html), and taper the edges with [st.taper()](https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.taper.html). \n",
|
||||
"Different types of filter are also available in [st.filter()](https://docs.obspy.org/packages/autogen/obspy.core.stream.Stream.filter.html). \n",
|
||||
"\n",
|
||||
"* first **detrend** the data. \n",
|
||||
"* next, apply a **bandpass filter** to select the frequencies with most noise energy. The secondary microseismic peak is roughly between 0.1 and 0.2 Hz. The primary microseismic peak between 0.05 and 0.1 Hz. Make sure to use a zero phase filter! *(specify argument zerophase=True)*"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Preprocessing 1\n",
|
||||
"\n",
|
||||
"stp = stn.copy() # copy stream\n",
|
||||
"t = stp[0].stats.starttime\n",
|
||||
"stp.trim(t, t + 4 * 86400) # shorten stream for quicker processing\n",
|
||||
"\n",
|
||||
"stp.detrend('linear') # remove trends using detrend\n",
|
||||
"stp.taper(max_percentage=0.05, type='cosine') # taper the edges\n",
|
||||
"stp.filter('bandpass', freqmin=0.1, freqmax=0.2, zerophase=True) # filter data of all traces in the streams"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"***Preprocessing 2***\n",
|
||||
"\n",
|
||||
"Some additional useful processing functions are provided in the following **Functions**\n",
|
||||
"\n",
|
||||
"* For each trace in the stream, apply **spectral whitening** on the frequency range you chose before (either [0.1 0.2]Hz or [0.05 0.1]Hz), using function **``whiten``**. \n",
|
||||
"\n",
|
||||
"\n",
|
||||
"* For the **time normalization**, the simplest option is to use the one-bit normalization option provided in function **``normalize``**. \n",
|
||||
"\n",
|
||||
"* *Optional: play around with different normalization options, such as clipping to a certain number of standard deviations, or using the running absolute mean normalization.*"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"A brief *desription of individual* **functions** (see the next cell) are as follows:\n",
|
||||
"\n",
|
||||
"1) **whiten**:\n",
|
||||
"\n",
|
||||
" spectral whitening of trace `tr` using a cosine tapered boxcar between `freqmin` and `freqmax`\n",
|
||||
" (courtesy Gaia Soldati & Licia Faenza, INGV)\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0,
|
||||
3,
|
||||
46,
|
||||
57,
|
||||
93,
|
||||
130,
|
||||
150
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def whiten(tr, freqmin, freqmax):\n",
|
||||
" \n",
|
||||
" nsamp = tr.stats.sampling_rate\n",
|
||||
" \n",
|
||||
" n = len(tr.data)\n",
|
||||
" if n == 1:\n",
|
||||
" return tr\n",
|
||||
" else: \n",
|
||||
" frange = float(freqmax) - float(freqmin)\n",
|
||||
" nsmo = int(np.fix(min(0.01, 0.5 * (frange)) * float(n) / nsamp))\n",
|
||||
" f = np.arange(n) * nsamp / (n - 1.)\n",
|
||||
" JJ = ((f > float(freqmin)) & (f<float(freqmax))).nonzero()[0]\n",
|
||||
" \n",
|
||||
" # signal FFT\n",
|
||||
" FFTs = np.fft.fft(tr.data)\n",
|
||||
" FFTsW = np.zeros(n) + 1j * np.zeros(n)\n",
|
||||
"\n",
|
||||
" # Apodization to the left with cos^2 (to smooth the discontinuities)\n",
|
||||
" smo1 = (np.cos(np.linspace(np.pi / 2, np.pi, nsmo+1))**2)\n",
|
||||
" FFTsW[JJ[0]:JJ[0]+nsmo+1] = smo1 * np.exp(1j * np.angle(FFTs[JJ[0]:JJ[0]+nsmo+1]))\n",
|
||||
"\n",
|
||||
" # boxcar\n",
|
||||
" FFTsW[JJ[0]+nsmo+1:JJ[-1]-nsmo] = np.ones(len(JJ) - 2 * (nsmo+1))\\\n",
|
||||
" * np.exp(1j * np.angle(FFTs[JJ[0]+nsmo+1:JJ[-1]-nsmo]))\n",
|
||||
"\n",
|
||||
" # Apodization to the right with cos^2 (to smooth the discontinuities)\n",
|
||||
" smo2 = (np.cos(np.linspace(0., np.pi/2., nsmo+1))**2.)\n",
|
||||
" espo = np.exp(1j * np.angle(FFTs[JJ[-1]-nsmo:JJ[-1]+1]))\n",
|
||||
" FFTsW[JJ[-1]-nsmo:JJ[-1]+1] = smo2 * espo\n",
|
||||
"\n",
|
||||
" whitedata = 2. * np.fft.ifft(FFTsW).real\n",
|
||||
" \n",
|
||||
" tr.data = np.require(whitedata, dtype=\"float32\")\n",
|
||||
"\n",
|
||||
" return tr"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"2) **correlateNoise**:\n",
|
||||
" \n",
|
||||
" correlate two stations, using slices of 'corrwin' seconds at a time correlations are also stacked. \n",
|
||||
" NB hardcoded: correlates 1st with 2nd station in the stream only signals are merged - any data gaps are\n",
|
||||
" filled with zeros.\n",
|
||||
" st : stream containing data from the two stations to correlate\n",
|
||||
" stations : list of stations\n",
|
||||
" corrwin : correlation window length\n",
|
||||
" returns 'corr' (all correlations) and 'stack' (averaged correlations)\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0,
|
||||
3,
|
||||
46,
|
||||
57,
|
||||
93,
|
||||
130,
|
||||
150
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def correlateNoise(st, stations, corrwin):\n",
|
||||
"\n",
|
||||
" print ('correlating stations', (stations[0], stations[1]))\n",
|
||||
"\n",
|
||||
" # initialize sliding timewindow (length = corrwin) for correlation\n",
|
||||
" # start 1 corrwin after the start to account for different stream lengths\n",
|
||||
" timewin = st.select(station=stations[1])[0].stats.starttime + corrwin\n",
|
||||
"\n",
|
||||
" # loop over timewindows \n",
|
||||
" # stop 1 corrwin before the end to account for different stream lengths\n",
|
||||
" while timewin < st.select(station=stations[0])[-1].stats.endtime - 2*corrwin:\n",
|
||||
" sig1 = st.select(station=stations[0]).slice(timewin, timewin+corrwin)\n",
|
||||
" sig1.merge(method=0, fill_value=0)\n",
|
||||
" sig2 = st.select(station=stations[1]).slice(timewin, timewin+corrwin)\n",
|
||||
" sig2.merge(method=0, fill_value=0)\n",
|
||||
" xcorr = np.correlate(sig1[0].data, sig2[0].data, 'same')\n",
|
||||
"\n",
|
||||
" try: \n",
|
||||
" # build array with all correlations\n",
|
||||
" corr = np.vstack((corr, xcorr))\n",
|
||||
" except: \n",
|
||||
" # if corr doesn't exist yet\n",
|
||||
" corr = xcorr\n",
|
||||
" \n",
|
||||
" # shift timewindow by one correlation window length\n",
|
||||
" timewin += corrwin\n",
|
||||
"\n",
|
||||
" # stack the correlations; normalize\n",
|
||||
" stack = np.sum(corr, 0)\n",
|
||||
" stack = stack / float((np.abs(stack).max())) \n",
|
||||
" print (\"...done\")\n",
|
||||
"\n",
|
||||
" return corr, stack"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"3) **plotStack**:\n",
|
||||
"\n",
|
||||
" plots stack of correlations with correct time axis\n",
|
||||
" st: stream containing noise (and station information)\n",
|
||||
" stack: array containing stack \n",
|
||||
" maxlag: maximum length of correlation to plot (in seconds)\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0,
|
||||
3,
|
||||
46,
|
||||
57,
|
||||
93,
|
||||
130,
|
||||
150
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def plotStack(st, stack, maxlag, figurename=None):\n",
|
||||
"\n",
|
||||
" # define the time vector for the correlation (length of corr = corrwin + 1)\n",
|
||||
" limit = (len(stack) / 2.) * st[0].stats.delta\n",
|
||||
" timevec = np.arange(-limit, limit, st[0].stats.delta)\n",
|
||||
"\n",
|
||||
" plt.plot(timevec, stack, 'k')\n",
|
||||
" stations = list(set([_i.stats.station for _i in st]))\n",
|
||||
" plt.title(\"Stacked correlation between %s and %s\" % (stations[0], stations[1]))\n",
|
||||
" plt.xlim(-maxlag, maxlag)\n",
|
||||
" plt.xlabel('time [s]')\n",
|
||||
"\n",
|
||||
" if figurename is not None:\n",
|
||||
" fig.savefig(figurename, format=\"pdf\")\n",
|
||||
" else:\n",
|
||||
" plt.show() "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"4) **plotXcorrEvent**:\n",
|
||||
"\n",
|
||||
" plot the noise correlation (MLAC, PHL) alongside the 1998 event signal \n",
|
||||
" st : event stream\n",
|
||||
" stn : noise stream\n",
|
||||
" stack : noise correlation array\n",
|
||||
" maxlag : maximum length of correlation, in seconds\n",
|
||||
" acausal : set to True to use acausal part (=negative times) of the correlation\n",
|
||||
" figurename : if a filename is specified, figure is saved in pdf format\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0,
|
||||
3,
|
||||
46,
|
||||
57,
|
||||
93,
|
||||
130,
|
||||
150
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def plotXcorrEvent(st, stn, stack, maxlag, acausal=False, figurename=None):\n",
|
||||
"\n",
|
||||
" eventtime = UTCDateTime(1998,7,15,4,53,21,0) # event near MLAC\n",
|
||||
"\n",
|
||||
" # station locations\n",
|
||||
" latP, lonP = 35.41, -120.55 # station PHL\n",
|
||||
" latM, lonM = 37.63, -118.84 # station MLAC\n",
|
||||
" latE, lonE = 37.55, -118.809 # event 1998\n",
|
||||
" \n",
|
||||
" # calculate distance between stations\n",
|
||||
" dist = gps2DistAzimuth(latP, lonP, latM, lonM)[0] # between PHL and MLAC\n",
|
||||
" distE = gps2DistAzimuth(latP, lonP, latE, lonE)[0] # between event and PHL\n",
|
||||
" #\n",
|
||||
" # CROSSCORRELATION\n",
|
||||
" # reverse stack to plot acausal part (= negative times of correlation)\n",
|
||||
" if acausal:\n",
|
||||
" stack = stack[::-1]\n",
|
||||
" \n",
|
||||
" # find center of stack\n",
|
||||
" c = int(np.ceil(len(stack)/2.) + 1)\n",
|
||||
" \n",
|
||||
" #cut stack to maxlag\n",
|
||||
" stack = stack[c - maxlag * int(np.ceil(stn[0].stats.sampling_rate)) : c + maxlag * int(np.ceil(stn[0].stats.sampling_rate))]\n",
|
||||
" \n",
|
||||
" # find new center of stack\n",
|
||||
" c2 = int(np.ceil(len(stack)/2.) + 1)\n",
|
||||
"\n",
|
||||
" # define time vector for cross correlation\n",
|
||||
" limit = (len(stack) / 2.) * stn[0].stats.delta\n",
|
||||
" timevec = np.arange(-limit, limit, stn[0].stats.delta)\n",
|
||||
" # define timevector: dist / t\n",
|
||||
" timevecDist = dist / timevec\n",
|
||||
" \n",
|
||||
" # EVENT\n",
|
||||
" ste = st.copy()\n",
|
||||
" st_PHL_e = ste.select(station='PHL')\n",
|
||||
" \n",
|
||||
" # cut down event trace to 'maxlag' seconds\n",
|
||||
" dt = len(stack[c2:])/stn[0].stats.sampling_rate #xcorrlength\n",
|
||||
" st_PHL_e[0].trim(eventtime, eventtime + dt)\n",
|
||||
" \n",
|
||||
" # create time vector for event signal\n",
|
||||
" # extreme values:\n",
|
||||
" limit = st_PHL_e[0].stats.npts * st_PHL_e[0].stats.delta\n",
|
||||
" timevecSig = np.arange(0, limit, st_PHL_e[0].stats.delta)\n",
|
||||
"\n",
|
||||
" # PLOTTING\n",
|
||||
" fig = plt.figure(figsize=(12.0, 8.0))\n",
|
||||
" ax1 = fig.add_subplot(2,1,1)\n",
|
||||
" ax2 = fig.add_subplot(2,1,2)\n",
|
||||
"\n",
|
||||
" # plot noise correlation\n",
|
||||
" ax1.plot(timevecDist[c2:], stack[c2:], 'k')\n",
|
||||
" ax1.set_title('Noise correlation between MLAC and PHL')\n",
|
||||
"\n",
|
||||
" # plot event near MLAC measured at PHL\n",
|
||||
" ax2.plot(distE/timevecSig, st_PHL_e[0].data / np.max(np.abs(st_PHL_e[0].data)), 'r')\n",
|
||||
" ax2.set_title('Event near MLAC observed at PHL')\n",
|
||||
"\n",
|
||||
" ax2.set_xlim((0, 8000))\n",
|
||||
" ax1.set_xlim((0, 8000))\n",
|
||||
"\n",
|
||||
" ax2.set_xlabel(\"group velocity [m/s]\")\n",
|
||||
" \n",
|
||||
" if figurename is not None:\n",
|
||||
" fig.savefig(figurename, format=\"pdf\")\n",
|
||||
" else:\n",
|
||||
" plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"5) **Normalize**:\n",
|
||||
"\n",
|
||||
" Temporal normalization of the traces, most after Bensen 2007. NB. before this treatment, traces must be\n",
|
||||
" demeaned, detrended and filtered. Description of argument:\n",
|
||||
"\n",
|
||||
" norm_method=\"clipping\"\n",
|
||||
" signal is clipped to 'clip_factor' times the std\n",
|
||||
" clip_factor recommended: 1 (times std)\n",
|
||||
" \n",
|
||||
" norm_method=\"clipping_iter\"\n",
|
||||
" the signal is clipped iteratively: values above 'clip_factor * std' \n",
|
||||
" are divided by 'clip_weight'. until the whole signal is below \n",
|
||||
" 'clip_factor * std'\n",
|
||||
" clip_factor recommended: 6 (times std)\n",
|
||||
" \n",
|
||||
" \n",
|
||||
" norm_method=\"ramn\"\n",
|
||||
" running absolute mean normalization: a sliding window runs along the \n",
|
||||
" signal. The values within the window are used to calculate a \n",
|
||||
" weighting factor, and the center of the window is scaled by this \n",
|
||||
" factor. \n",
|
||||
" weight factor: w = np.mean(np.abs(tr.data[win]))/(2. * norm_win + 1) \n",
|
||||
" finally, the signal is tapered with a tukey window (alpha = 0.2).\n",
|
||||
"\n",
|
||||
" norm_win: running window length, in seconds.\n",
|
||||
" recommended: half the longest period\n",
|
||||
"\n",
|
||||
" norm_method=\"1bit\"\n",
|
||||
" only the sign of the signal is conserved\n",
|
||||
" "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0,
|
||||
3,
|
||||
46,
|
||||
57,
|
||||
93,
|
||||
130,
|
||||
150
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Functions\n",
|
||||
"# collection of functions used in noise correlation processing\n",
|
||||
"\n",
|
||||
"def normalize(tr, clip_factor=6, clip_weight=10, norm_win=None, norm_method=\"1bit\"): \n",
|
||||
" \n",
|
||||
" if norm_method == 'clipping':\n",
|
||||
" lim = clip_factor * np.std(tr.data)\n",
|
||||
" tr.data[tr.data > lim] = lim\n",
|
||||
" tr.data[tr.data < -lim] = -lim\n",
|
||||
"\n",
|
||||
" elif norm_method == \"clipping_iter\":\n",
|
||||
" lim = clip_factor * np.std(np.abs(tr.data))\n",
|
||||
" \n",
|
||||
" # as long as still values left above the waterlevel, clip_weight\n",
|
||||
" while tr.data[np.abs(tr.data) > lim] != []:\n",
|
||||
" tr.data[tr.data > lim] /= clip_weight\n",
|
||||
" tr.data[tr.data < -lim] /= clip_weight\n",
|
||||
"\n",
|
||||
" elif norm_method == 'ramn':\n",
|
||||
" lwin = tr.stats.sampling_rate * norm_win\n",
|
||||
" st = 0 # starting point\n",
|
||||
" N = lwin # ending point\n",
|
||||
"\n",
|
||||
" while N < tr.stats.npts:\n",
|
||||
" win = tr.data[st:N]\n",
|
||||
"\n",
|
||||
" w = np.mean(np.abs(win)) / (2. * lwin + 1)\n",
|
||||
" \n",
|
||||
" # weight center of window\n",
|
||||
" tr.data[st + lwin / 2] /= w\n",
|
||||
"\n",
|
||||
" # shift window\n",
|
||||
" st += 1\n",
|
||||
" N += 1\n",
|
||||
"\n",
|
||||
" # taper edges\n",
|
||||
" taper = get_window(tr.stats.npts)\n",
|
||||
" tr.data *= taper\n",
|
||||
"\n",
|
||||
" elif norm_method == \"1bit\":\n",
|
||||
" tr.data = np.sign(tr.data)\n",
|
||||
" tr.data = np.float32(tr.data)\n",
|
||||
"\n",
|
||||
" return tr"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"6) **get_window**:\n",
|
||||
"\n",
|
||||
" Return tukey window of length N\n",
|
||||
" N: length of window\n",
|
||||
" alpha: alpha parameter in case of tukey window.\n",
|
||||
" 0 -> rectangular window\n",
|
||||
" 1 -> cosine taper\n",
|
||||
" returns: window (np.array)\n",
|
||||
" \n",
|
||||
"Doc of [scipy.signal.get_window](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.get_window.html)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0,
|
||||
3,
|
||||
46,
|
||||
57,
|
||||
93,
|
||||
130,
|
||||
150
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def get_window(N, alpha=0.2):\n",
|
||||
"\n",
|
||||
" window = np.ones(N)\n",
|
||||
" x = np.linspace(-1., 1., N)\n",
|
||||
" ind1 = (abs(x) > 1 - alpha) * (x < 0)\n",
|
||||
" ind2 = (abs(x) > 1 - alpha) * (x > 0)\n",
|
||||
" window[ind1] = 0.5 * (1 - np.cos(np.pi * (x[ind1] + 1) / alpha))\n",
|
||||
" window[ind2] = 0.5 * (1 - np.cos(np.pi * (x[ind2] - 1) / alpha))\n",
|
||||
" return window"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "slide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"### Actual preprocessing happens here -- this can take a while!"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Preprocessing 2\n",
|
||||
"st = stp.copy() # copy stream\n",
|
||||
"\n",
|
||||
"for tr in st:\n",
|
||||
" tr = normalize(tr, norm_method=\"1bit\")\n",
|
||||
" tr = whiten(tr, 0.1, 0.2)\n",
|
||||
"print ('done!')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"#### Cross-correlation ####\n",
|
||||
"Once you're happy with the preprocessing, you can calculate the **cross-correlation** using **``correlateNoise``** function. The cross-correlation are computed by slices of a few hours each (specified in *corrwin*). \n",
|
||||
"\n",
|
||||
"**For correlateNoise function**\n",
|
||||
"* input: stream, list of stations (here: ['MLAC', 'PHL']), slice length in seconds\n",
|
||||
"* output: all individual correlations, stack\n"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Cross-correlate\n",
|
||||
"xcorr, stack = correlateNoise(st, ['MLAC','PHL'], 7200)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"The resulting stack can be **plotted** with **``plotStack``** function. Since it doesn't make much sense to look at a 2 hour long correlation signal, you can decide to plot only the central part by specifying a ``maxlag`` (in seconds). "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Plotting\n",
|
||||
"\n",
|
||||
"plotStack(st, stack, 400)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"If you're only working with a few days of noise (after trimming), this plot probably doesn't look very nice. You could go back to the code block named 'preprocessing 1', and keep a longer noise record (10 days works quite well already). "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "slide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"#### Compare to event trace ####\n",
|
||||
"In 1998, a M = 5.1 event occurred next to station MLAC. This event was recorded at PHL and we read this data.\n",
|
||||
"\n",
|
||||
"* **read** the event data to a separate stream"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"ste = read('https://raw.github.com/ashimrijal/NoiseCorrelation/master/data/event.CI.PHL.LHZ.1998.196.1998.196.mseed')\n",
|
||||
"# if data is stored locally, uncomment the following line and comment the line above:\n",
|
||||
"#ste = obspy.read('./event.CI.PHL.LHZ.1998.196.1998.196.mseed')"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "subslide"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"#### Preprocess event ####\n",
|
||||
"\n",
|
||||
"The event signal should be processed in a similar way to the noise. \n",
|
||||
"\n",
|
||||
"* **detrend** in the same way as before\n",
|
||||
"* **bandpass filter** for the same frequencies as chosen above\n",
|
||||
"* apply **spectral whitening** to each trace in the event stream to ensure the spectrum is comparable to the noise. "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Preprocessing\n",
|
||||
"\n",
|
||||
"ste.detrend('linear')\n",
|
||||
"ste.filter('bandpass', freqmin=0.1, freqmax=0.2, zerophase=True)\n",
|
||||
"for tr in ste:\n",
|
||||
" tr = whiten(tr, 0.1, 0.2)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"source": [
|
||||
"#### Plot! ####\n",
|
||||
"\n",
|
||||
"A plotting function is provided to plot both signals alongside: **``plotXcorrEvent``**. \n",
|
||||
"\n",
|
||||
"* input: event stream, noise stream, stack, maxlag"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"code_folding": [
|
||||
0
|
||||
],
|
||||
"slideshow": {
|
||||
"slide_type": "fragment"
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Plotting\n",
|
||||
"\n",
|
||||
"plotXcorrEvent(ste, stn, stack, 400)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"celltoolbar": "Slideshow",
|
||||
"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.8.2"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 2
|
||||
}
|
||||
230
04-FFT_DFT_and_Applications/resample.ipynb
Normal file
230
04-FFT_DFT_and_Applications/resample.ipynb
Normal file
@@ -0,0 +1,230 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Resampling of Time Series\n",
|
||||
"## Decimation\n",
|
||||
"The Discrete Fourier Transform (DFT) can be used for some very basic signal processing operations. We can either use is for decimation, which results in redcution of the number of samples of a given time series or we can use it for interpolation, which is the opposite of decimation.\n",
|
||||
"Decimation implies a reduction of time resolution and thus a decrease in the Nyquist frequency.\n",
|
||||
"If the original signal contains higher frequency than the new Nyquist frequency\n",
|
||||
"(after decimation) we get aliasing. Thus, a low pass filtering to the new Nyquist\n",
|
||||
"frequency is required before decimation. Decimation and low-pass filtering can be\n",
|
||||
"achieved in one step by using the DFT:\n",
|
||||
"1. Transform the time series to the frequency domain using the DFT.\n",
|
||||
"2. Reduce the Nyquist frequency by the desired factor (when using the Fast Fourier Transform, this should be a power of 2), and apply an appropriate taper. That’s the low-pass filter!\n",
|
||||
"3. Do an inverse DFT of the modified spectrum with reduced Nyquist frequency and obtain a time series with greater sampling interval (because lower Nyquist frequency implies lower time resolution)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Import all required packages\n",
|
||||
"%matplotlib inline\n",
|
||||
"import os\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from obspy import read\n",
|
||||
"from scipy import signal"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Read the data\n",
|
||||
"datapath = os.path.join(os.path.expanduser('~'),'work', 'data', 'DOR50', '*063')\n",
|
||||
"st = read(datapath)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def decimate(x, dt, factor):\n",
|
||||
" \"\"\"\n",
|
||||
" Easy decimation for real input data x. Note, since new_dt might not be old_dt * 2**n,\n",
|
||||
" the decimation is not correct.\n",
|
||||
" x: Data array\n",
|
||||
" dt: sampling rate of x\n",
|
||||
" factor: factor for decimation, thus new sampling rate is dt/factor\n",
|
||||
" \"\"\"\n",
|
||||
" # 1. Transform x to frequnecy domain\n",
|
||||
" ft = np.fft.rfft(x)\n",
|
||||
" freqs = np.fft.rfftfreq(n=len(x), d=dt)\n",
|
||||
" \n",
|
||||
" # 2. Get Nyquist Frequencies for both old dt an new dt\n",
|
||||
" \n",
|
||||
" \n",
|
||||
" # Find index of new Nyquist Frequnency in freqs\n",
|
||||
" index = np.where(freqs > new_f_ny)[0][0]\n",
|
||||
" \n",
|
||||
" # 3. Apply window function as lowpass filter and do inverse DFT \n",
|
||||
" new_ft = ft[:index]\n",
|
||||
" new_ft = new_ft * signal.tukey(len(new_ft), alpha=0.25)\n",
|
||||
" x_new = np.fft.irfft(ft[:index])\n",
|
||||
" \n",
|
||||
" return x_new, dt*factor"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Plot results for original and decimated data\n",
|
||||
"# Read data from st\n",
|
||||
"data = st[0].data[:100000]\n",
|
||||
"\n",
|
||||
"# Apply Decimation\n",
|
||||
"dec, dt_dec = decimate(data, dt=st[0].stats.delta, factor=2)\n",
|
||||
"\n",
|
||||
"# Create arrays for time to plot\n",
|
||||
"t_data = np.arange(0, len(data)) * st[0].stats.delta\n",
|
||||
"t_dec = np.arange(0, len(dec)) * dt_dec\n",
|
||||
"\n",
|
||||
"# Create plot widget\n",
|
||||
"plt.plot(t_data, data - np.mean(data), alpha=0.5)\n",
|
||||
"plt.plot(t_dec, dec - np.mean(dec), alpha=0.5)\n",
|
||||
"plt.xlabel(\"Time (s)\")\n",
|
||||
"plt.ylabel(\"Amplitude (a.u.)\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Interpolation\n",
|
||||
"The opposite of decimation is interpolation. Here, we want a finer sampling of the\n",
|
||||
"time series without changing the frequency content. We could do this in the time\n",
|
||||
"domain by some interpolation rule using neighbouring samples. We can also do it\n",
|
||||
"by using the DFT:\n",
|
||||
"1. Transform the time series to the frequency domain using the DFT.\n",
|
||||
"2. Append zeros to the spectrum thus increasing the Nyquist frequency.\n",
|
||||
"3. Do an inverse DFT of the extended spectrum and obtain a time series with smaller sampling interval (because higher Nyquist frequency implies higher time resolution). The frequency content is unchanged!"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def interpolation(x, dt, new_dt):\n",
|
||||
" \"\"\"\n",
|
||||
" Easy function for interpolation of a given data array to increase the the sampling rate.\n",
|
||||
" \n",
|
||||
" x: data array\n",
|
||||
" dt: sampling rate\n",
|
||||
" new_dt: new sampling rae\n",
|
||||
" \"\"\"\n",
|
||||
" # 1. Transform x to frequnecy domain\n",
|
||||
" ft = np.fft.rfft(x)\n",
|
||||
" \n",
|
||||
" # Determine the number of zeros to add\n",
|
||||
" \n",
|
||||
" \n",
|
||||
" # 2. Append zeros to spectrum\n",
|
||||
" ft = np.concatenate((ft, np.zeros(append_zeros)))\n",
|
||||
" \n",
|
||||
" # 3. Do inverse DFT and determine new dt\n",
|
||||
" x_new = np.fft.irfft(ft)\n",
|
||||
" new_dt = t_max / len(x_new)\n",
|
||||
" \n",
|
||||
" return x_new, new_dt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Create plot\n",
|
||||
"# Read data from st\n",
|
||||
"data = st[0].data[:100000]\n",
|
||||
"\n",
|
||||
"# Apply interpolation\n",
|
||||
"interp, dt_interp = interpolation(data, dt=st[0].stats.delta, new_dt=1/150)\n",
|
||||
"\n",
|
||||
"# Create time arrays\n",
|
||||
"t_data = np.arange(0, len(data)) * st[0].stats.delta\n",
|
||||
"t_interp = np.arange(0, len(interp)) * dt_interp\n",
|
||||
"\n",
|
||||
"# Create plot widget\n",
|
||||
"plt.plot(t_data, data - np.mean(data), alpha=0.5)\n",
|
||||
"plt.plot(t_interp, interp - np.mean(interp), alpha=0.5)\n",
|
||||
"plt.xlabel(\"Time (s)\")\n",
|
||||
"plt.ylabel(\"Amplitude (a.u.)\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Comparison with function scipy.signal.resample\n",
|
||||
"For documentation see https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Import funtion\n",
|
||||
"from scipy.signal import resample\n",
|
||||
"\n",
|
||||
"# Read data from st and create time array\n",
|
||||
"data = st[0].data[:100000]\n",
|
||||
"time = np.arange(0, len(data)) * st[0].stats.delta\n",
|
||||
"\n",
|
||||
"# Apply resample funtion\n",
|
||||
"resampled_x, resampled_t = resample(data, num=int(len(data)*2), t=time)\n",
|
||||
"\n",
|
||||
"# Create plot\n",
|
||||
"plt.plot(time, data - np.mean(data), alpha=0.5)\n",
|
||||
"plt.plot(resampled_t, resampled_x - np.mean(resampled_x), alpha=0.5)\n",
|
||||
"plt.xlabel(\"Time (s)\")\n",
|
||||
"plt.ylabel(\"Amplitude (a.u.)\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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
|
||||
}
|
||||
307
05-Spectrogram/multipleFilterTechnique.ipynb
Normal file
307
05-Spectrogram/multipleFilterTechnique.ipynb
Normal file
@@ -0,0 +1,307 @@
|
||||
{
|
||||
"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
|
||||
}
|
||||
165
05-Spectrogram/nonlinear_thresholding.ipynb
Normal file
165
05-Spectrogram/nonlinear_thresholding.ipynb
Normal file
@@ -0,0 +1,165 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Nonlinear Thresholding\n",
|
||||
"When noise and signal share a common frequency band, spectral filtering would lead to a loss in signal. For these cases, we need different techniques to reduce the disturibing noise. We start with a reocorde signal $x(t)$, contains some signal $s(t)$ and additive noise $n(t)$:\n",
|
||||
"$$\n",
|
||||
"x(t) = s(t) + n(t)\n",
|
||||
"$$\n",
|
||||
"We can rewrite the equation above in time-frequency domain as\n",
|
||||
"$$\n",
|
||||
"X(t, f) = S(t, f) + N(t, f) ~.\n",
|
||||
"$$\n",
|
||||
"Assuming that some time-frequnecy coefficients can be associated by noise, we set all coefficients, which a below this threshold, to zero:\n",
|
||||
"$$\n",
|
||||
"\\tilde{X}(t, f) = \\left\\{\n",
|
||||
" \\begin{array}{@{}ll@{}}\n",
|
||||
" X(t, f) & \\mathrm{if}~ |X(t, f)| \\geq \\beta(f) \\\\\n",
|
||||
" 0 & \\mathrm{otherwise}\n",
|
||||
" \\end{array}\\right.\n",
|
||||
" ~,\n",
|
||||
"$$\n",
|
||||
"where $\\beta(f)$ denotes the threshold function. The threshold function is defined as \n",
|
||||
"$$\n",
|
||||
"\\beta(f) = \\mathrm{ECDF}_f^{-1} (P = 0.99) ~,\n",
|
||||
"$$\n",
|
||||
"where $\\mathrm{ECDF}_f^{-1}$ denotes the inverse cumulative distribution function or quantile function, e.g. before the first arrival of earthquake waves."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Load all required packages\n",
|
||||
"import os\n",
|
||||
"import numpy as np\n",
|
||||
"from scipy.signal import stft, istft\n",
|
||||
"import matplotlib.pyplot as plt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Load our earthquake data and plot\n",
|
||||
"d = np.load(os.path.join(os.path.expanduser('~'),'work', 'data', 'events', 'bug2019mgoh_Z.npz'))\n",
|
||||
"#d = np.load(os.path.join(os.path.expanduser('~'),'work', 'data', 'events', 'bug2019gbbo_Z.npz'))\n",
|
||||
"#d = np.load(os.path.join(os.path.expanduser('~'),'work', 'data', 'events', 'bug2019ibsd_Z.npz'))\n",
|
||||
"plt.figure(figsize=(15, 8))\n",
|
||||
"plt.plot(d['data']);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Compute spectrogram of data\n",
|
||||
"f, t, X = stft(d['data'], fs=1/0.01, nfft=198, nperseg=99) # Be careful with choice of nfft and nperseg, \n",
|
||||
" # because istft does not work for all pairs\n",
|
||||
"print(X.shape)\n",
|
||||
"plt.figure(figsize=(15, 8))\n",
|
||||
"plt.pcolormesh(t, f, np.abs(X))\n",
|
||||
"plt.xlabel(\"Frequency (Hz)\")\n",
|
||||
"plt.ylabel(\"Time (s)\");"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Functions for thresholding\n",
|
||||
"def threshold(X, quantile=0.99):\n",
|
||||
" \"\"\"\n",
|
||||
" X: time-frequency coefficients\n",
|
||||
" q: quantile [0-1]\n",
|
||||
" \"\"\"\n",
|
||||
" # Loop over frequencies to build threshold function\n",
|
||||
" beta = np.zeros(X.shape[0])\n",
|
||||
" for i in range(X.shape[0]):\n",
|
||||
" beta[i] = np.quantile(np.abs(X[i, :]), q=quantile)\n",
|
||||
" \n",
|
||||
" return beta\n",
|
||||
"\n",
|
||||
"def modify_spectrogram(X, beta):\n",
|
||||
" # Loop over all items in X and apply threshold\n",
|
||||
" # Task: modify X!\n",
|
||||
" \n",
|
||||
" return X"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Compute STFT of data before first arrival of P\n",
|
||||
"_, _, X_thres = stft(d['data'][:2500], fs=1/0.01, nfft=198, nperseg=99)\n",
|
||||
"# Estimate threshold fucntion\n",
|
||||
"beta = threshold(X_thres)\n",
|
||||
"\n",
|
||||
"# Modifiy original spectrogram\n",
|
||||
"X_mod = modify_spectrogram(X, beta)\n",
|
||||
"\n",
|
||||
"# Plot modified spectrogram\n",
|
||||
"plt.figure(figsize=(15, 8))\n",
|
||||
"plt.pcolormesh(t, f, np.abs(X_mod))\n",
|
||||
"plt.xlabel(\"Frequency (Hz)\")\n",
|
||||
"plt.ylabel(\"Time (s)\");"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Inverse STFT of modified spectrogram\n",
|
||||
"t, x_mod = istft(X_mod, fs=1/0.01, nfft=198, nperseg=99)\n",
|
||||
"\n",
|
||||
"# Plot corrected seismogram\n",
|
||||
"plt.figure(figsize=(15, 8))\n",
|
||||
"plt.plot(t, x_mod);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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
|
||||
}
|
||||
272
05-Spectrogram/time_frequency_analysis.ipynb
Normal file
272
05-Spectrogram/time_frequency_analysis.ipynb
Normal file
@@ -0,0 +1,272 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Time Frequency Analysis\n",
|
||||
"## Moving window analysis\n",
|
||||
"One way to analyse the time-varying frequency content of a signal is to\n",
|
||||
"apply windows in the time domain to the signal and to calculate a Fourier spectrum\n",
|
||||
"of the windowed part. The window marches along the signal with defined overlap creating\n",
|
||||
"a series of Fourier spectra associated with the center times of the windows. The resulting amplitude\n",
|
||||
"spectra are the plotted versus window center time. In more detail:\n",
|
||||
"\n",
|
||||
"1. Choose windowing functions: $w(t,t_m)$ with $t_m$ the center of the window.\n",
|
||||
"2. Multiply windowing function with time series: $f_m(t) = f(t)w(t,t_m)$\n",
|
||||
"3. Detrend the windowed signal.\n",
|
||||
"4. Perform a DFT: $F_{km} = \\Delta t\\sum_{n=0}^N f_m(t)\\exp(-2\\pi i \\frac{kn}{N})$\n",
|
||||
" and calculate the absolute value, $|F_{km}|$.\n",
|
||||
"5. Plot the resulting matrix: $|F_{km}|$ in the time-frequency plane."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Borehole Data and Test Site\n",
|
||||
"\n",
|
||||
"Within the scope of constructing new buildings for the International Geothermal Center in Bochum in 2013 wells were drilled for the installation of a geothermal heat exchanger. Bore holes were drilled next to the newly constructed building (Station GEO3). A downhole hammer was used with a diameter of 152 mm. The drill bit operates by water flushing through the drill rod. The water flow rate determines the working frequency of the hammer.\n",
|
||||
"\n",
|
||||
"To observe the drill bit noise a temporarily operating 2-D seismic network was installed around the drill site. Here, the noise of the used downhole hammer is investigated. An array of 16 seismological stations was installed in the test site. Four Mark L-4C-3D 1 Hz sensors, eight S-13 1 Hz sensors, one GS-13 1 Hz sensor and two Güralp CMG-3ESPC broad-band 120 sec – 50 Hz sensors were in use. Additionally an accelerometer was fixed to the drill rod (GEO11, blue triangle).\n",
|
||||
"\n",
|
||||
"Some of the stations were positioned within one of the infrastructure tunnels servicing the university containing water conduits, long-distance heat line and electric cables (e.g. Station GEO4 and GEO05). Thus, noise could be reduced that might disturb the recordings. Other stations are located within buildings (e.g. Station GEO2 and GEO03) ore outside (e.g. Station GEO6 and GEO07).\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"Drill bit noise was recorded up to the maximum drilling depth of 200 m. A drilling cycle is characterised by switching on the water pump, followed by the drilling with higher amplitude signals, that lasts several minutes. The water pump\n",
|
||||
"is stopped about 5 to 15 minutes after the drilling finished depending on the drill depth."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Preparation: load packages\n",
|
||||
"import os\n",
|
||||
"from obspy.core import read\n",
|
||||
"from obspy.core import UTCDateTime\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pylab as plt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# data from stations, Jan 23 2013\n",
|
||||
"\n",
|
||||
"record = '1400' # either use record starting at 1300 or 1400\n",
|
||||
"station = 'GEO7' # station code\n",
|
||||
"# available stations are GEO2, GEO3 (only 14:00), GEO4, GEO5, GEO6, GEO7\n",
|
||||
"\n",
|
||||
"if record =='1300':\n",
|
||||
" stime = UTCDateTime('2013-01-23 13:16:00Z') # use record starting at 1300\n",
|
||||
" etime = UTCDateTime('2013-01-23 13:25:00Z')\n",
|
||||
" ttitle = ', depth: 36.5 m'\n",
|
||||
"else:\n",
|
||||
" stime = UTCDateTime('2013-01-23 14:17:00Z') # use record starting at 1400 (14:14-14:23)\n",
|
||||
" etime = UTCDateTime('2013-01-23 14:19:00Z')\n",
|
||||
" ttitle = ', depth: 56.5 m'\n",
|
||||
" \n",
|
||||
"datapath = os.path.join(os.path.expanduser('~'),'work', 'data','GE-stations', station, 'e*'+record+'*.HHZ') # search string for database\n",
|
||||
"st = read(datapath) # read file using obspy read\n",
|
||||
"print(st)\n",
|
||||
"\n",
|
||||
"st.trim(stime,etime) # trim data stream to desired start and end time\n",
|
||||
"st.detrend('linear') # do a linear detrending\n",
|
||||
"st.detrend('demean') # subtract mean value\n",
|
||||
"tr = st[0] # extract first trace from stream (there is only one)\n",
|
||||
"tr.plot(); # plot trace"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Cell 2: Preparation of window\n",
|
||||
"wlen = 4. # window length in seconds\n",
|
||||
"npts = tr.stats.npts # number of samples\n",
|
||||
"dt = tr.stats.delta # sampling interval\n",
|
||||
"\n",
|
||||
"nwin = wlen/dt # number of samples in window\n",
|
||||
"nwin = int(pow(2, np.ceil(np.log(nwin)/np.log(2)))) # find next higher power of 2 of nwin\n",
|
||||
"nwmax = int(pow(2, np.ceil(np.log(npts)/np.log(2))))/8 # maximum window length (about 1/8 of length)\n",
|
||||
"nwin = min(nwin,nwmax) # limit nwin\n",
|
||||
"print(\"Samples in window: \",nwin)\n",
|
||||
"print(\"Total number of samples: \",npts)\n",
|
||||
"print(\"Sampling interval: \",dt)\n",
|
||||
"print(\"Length of trace [s]: \",npts*dt)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# The Hann window function\n",
|
||||
"def hann(nw):\n",
|
||||
" arg = 2.*np.pi*np.arange(0,nw)/nw # argument of cosine\n",
|
||||
" fwin = 0.5*(1.-np.cos(arg)) # Hann window\n",
|
||||
" return fwin\n",
|
||||
"\n",
|
||||
"# The boxcar window function\n",
|
||||
"def boxcar(nw):\n",
|
||||
" fwin = np.ones(nw)\n",
|
||||
" return fwin"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Plot the window function and its effect on the time series\n",
|
||||
"#\n",
|
||||
"fwin = hann(nwin) # Hann window\n",
|
||||
"plt.figure(figsize = [15,6])\n",
|
||||
"plt.subplot(221)\n",
|
||||
"plt.plot(fwin,'-k')\n",
|
||||
"plt.title(\"Hann window\")\n",
|
||||
"seg = tr.data[0:nwin]*fwin # multiply data with window\n",
|
||||
"plt.subplot(222)\n",
|
||||
"plt.plot(seg,'-k')\n",
|
||||
"fwin = boxcar(nwin) # Boxcar window\n",
|
||||
"plt.subplot(223)\n",
|
||||
"plt.plot(fwin,'-k')\n",
|
||||
"plt.title(\"Boxcar window\")\n",
|
||||
"seg = tr.data[0:nwin]*fwin # multiply data with window\n",
|
||||
"plt.subplot(224)\n",
|
||||
"plt.plot(seg,'-k')\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def movingWindowAnalysis(data,winfun,nwin,shift):\n",
|
||||
" \"\"\"\n",
|
||||
" Performs moving window analysis of a time series.\n",
|
||||
" data: data array\n",
|
||||
" winfun: name of the window function to be called\n",
|
||||
" nwin: number of window samples (power of 2)\n",
|
||||
" shift: displacement of moving window in samples\n",
|
||||
" \"\"\"\n",
|
||||
" fwin = winfun(nwin) # compute window values\n",
|
||||
" npts = len(data) # number of total samples\n",
|
||||
" nseg = int((npts-nwin)/shift)+1 # total number of expected data segment\n",
|
||||
" mwa = np.zeros((nwin//2+1,nseg)) # array for result (rfft returns N/2+1 samples) \n",
|
||||
" wa = 0 # start index of data segment\n",
|
||||
" we = nwin # end index of data segment\n",
|
||||
" jseg = 0 # initialize data segment counter\n",
|
||||
" while we < npts: # loop over segments\n",
|
||||
" seg = data[wa:we]*fwin # multiply data segment with window\n",
|
||||
" seg = seg-seg.mean() # subtract mean value of segment\n",
|
||||
" ftseg = np.abs(np.fft.rfft(seg)) # abs value of Fourier transform\n",
|
||||
" maxft = np.amax(ftseg) # max value of Fourier transform\n",
|
||||
" ftseg = ftseg/maxft+1.e-10 # normalize spectrum to its maximum, remove zeros\n",
|
||||
" mwa[:,jseg] = np.power(ftseg,0.25) # assign values to the matrix\n",
|
||||
" wa = wa+shift # move window start by shift\n",
|
||||
" we = we+shift # move window end by shift\n",
|
||||
" jseg = jseg+1 # increase segment counter\n",
|
||||
" return nseg,mwa # return number of segments and moving window matrix"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# carry out the moving window analysis\n",
|
||||
"tshift = 0.2*wlen\n",
|
||||
"shift = int(tshift/dt) # displacement of moving window in samples\n",
|
||||
"nseg,mwa = movingWindowAnalysis(tr.data,hann,nwin,shift) # compute spectrogram\n",
|
||||
"freq = np.fft.rfftfreq(nwin,dt) # Fourier frequencies\n",
|
||||
"print(\"Frequency range [Hz]: \",freq[0],freq[-1])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# plot the result\n",
|
||||
"fmaxplot = 250. # maximum frequency to be plotted\n",
|
||||
"nf = len(freq) # number of frequencies\n",
|
||||
"jfmax = int(fmaxplot/freq[-1]*nf) # max index for plotting \n",
|
||||
"extent = (0.5*wlen,0.5*wlen+nseg*tshift,fmaxplot,0.) # extent of matrix in true time and frequency\n",
|
||||
"plt.figure(figsize = [15,6])\n",
|
||||
"plt.imshow(mwa[0:jfmax,:],extent = extent,aspect=0.25) # do plotting\n",
|
||||
"plt.xlabel('Window center time [s]')\n",
|
||||
"plt.ylabel('Frequency [Hz]')\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Exercise \n",
|
||||
"\n",
|
||||
"1. Write a little routine for the Hamming window and do the analysis using this window\n",
|
||||
"2. Use the sqrt of the Fourier amplitude spectrum as output and compare\n",
|
||||
"3. Use other powers ($< 1$) of the Fourier amplitude (using np.power)\n",
|
||||
"4. Use the logarithm of the Fourier amplitude as output and compare (using np.log)\n",
|
||||
"5. Extend the frequency range when using the logarithm or a power of the Fourier amplitude\n",
|
||||
"6. Do a moving window analysis for the other available stations. Which ones appear good and which ones are bad?\n",
|
||||
"7. Try also records starting at 14:00\n",
|
||||
"\n",
|
||||
"Hint: when you change something in the code, use \"Run all below\" to see the changes."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### The Hamming window\n",
|
||||
"$ w(t) = 0.54 - 0.46 \\cos \\left( \\frac{2 \\pi t}{L} \\right) ~ , t \\in [0, L]$"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": []
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.6.7"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
90
05-Spectrogram/time_meets_frequency.ipynb
Normal file
90
05-Spectrogram/time_meets_frequency.ipynb
Normal file
@@ -0,0 +1,90 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Dispersive Signals\n",
|
||||
"Up to now,we considered the spectral content of the entire time series. In many\n",
|
||||
"cases, the spectral content varies with time in a time series and one often wants to\n",
|
||||
"find out how this happens in detail. One important example are dispersed signals in\n",
|
||||
"which energy with different frequency content arrives at different times. Most\n",
|
||||
"prominent among these are seismic surface which exhibit significant dispersion if\n",
|
||||
"epicentral distances between seismic station and source are large. Quantifying the\n",
|
||||
"dispersion allows inferences on Earth structure. \n",
|
||||
"\n",
|
||||
"<img src=\"../images/seismogram.png\" alt=\"Drawing\" style=\"width: 700px;\"/>\n",
|
||||
"<img src=\"../images/stack-surface-waves.png\" alt=\"Drawing\" style=\"width: 700px;\"/>"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## The Chirp or Sweep"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib inline\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"import numpy as np\n",
|
||||
"plt.rcParams['figure.figsize'] = 12, 8 # Slightly bigger plots by default"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"from scipy.signal import chirp\n",
|
||||
"t = np.linspace(0, 30, 30001)\n",
|
||||
"w = chirp(t, f0=0.01, f1=3, t1=np.max(t), method='linear')\n",
|
||||
"plt.plot(t, w);"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Tasks\n",
|
||||
"1. What is special about the Chirp Signal?\n",
|
||||
"2. Plot the amplitude spectrum of the chirp signal, using numpys [FFT module](https://numpy.org/doc/stable/reference/routines.fft.html), including the correct frequncies on the x-axis!\n",
|
||||
"3. Think of a solution how to visualize the increasing frequency content using the Fourier tansform"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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
|
||||
}
|
||||
276
05-Spectrogram/uncertainty_instantaneous_frequency.ipynb
Normal file
276
05-Spectrogram/uncertainty_instantaneous_frequency.ipynb
Normal file
@@ -0,0 +1,276 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# The Uncertainty Problem\n",
|
||||
"Using the moving window analysis, you may have noticed that we are not able to get a good resoltion in time and frequency. This is called the uncertainty problem, which you might have heard from in physics.\n",
|
||||
"\n",
|
||||
"Here, we a closer look how the uncertainty problem effects our sweep signal, when we apply the moving window average on this signal.\n",
|
||||
"\n",
|
||||
"### Tasks\n",
|
||||
"1. Read the documentation about [scipy's STFT](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.chirp.html), especially the parameters nperseg, noverlap and nfft\n",
|
||||
"2. Play around with the parameters"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Preparation: load packages\n",
|
||||
"import os\n",
|
||||
"from obspy.core import read\n",
|
||||
"from obspy.core import UTCDateTime\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from scipy.signal import stft, chirp\n",
|
||||
"import pywt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Load Chirp signal\n",
|
||||
"dt = 0.01\n",
|
||||
"t = np.arange(0, 5000) * dt\n",
|
||||
"data = chirp(t=t, f0=0.5, f1=25, t1=t[-1])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Applying Short-Time Fourier Transform (STFT) on data\n",
|
||||
"freqs, time, S = stft(x=data, fs=1/dt, window='hann', nperseg=256, noverlap=None, nfft=None)\n",
|
||||
"\n",
|
||||
"plt.figure(figsize = [15,6])\n",
|
||||
"plt.pcolormesh(time, freqs, np.abs(S))\n",
|
||||
"plt.xlabel(\"Time (s)\")\n",
|
||||
"plt.ylabel(\"Frequency (Hz)\");"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# The instantaneous Frequency\n",
|
||||
"The instantaneous frequency is a measure of the dominating frequency at each time instant of a signal and is determined as the time derivative of the instantaneous phase $\\phi (t)$. The phase can be obtained by constructing the analytical signal of an arbitrary time series:\n",
|
||||
"\n",
|
||||
"$$ \n",
|
||||
"\\hat{s}(t) = s(t) + i \\cdot \\mathcal{H} \\left( (s(t) \\right) = |\\hat{s}(t)| \\cdot \\exp (i \\phi (t)),\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"where $\\mathcal{H}$ denotes the Hilbert-Transform. Note, the scipy function [hilbert()](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.hilbert.html) computes the analytical signal. Now, we can derive the instantaneous phase as\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\phi (t) = \\arg \\left( \\frac{Im(\\hat{s}(t))}{Re(\\hat{s}(t))} \\right).\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"Note, the instantaneous phase is called wrapped phase if it is in the interval $(−\\pi, \\pi]$ or $[0, 2 \\pi)$. To determine the instantaneous frequncy, we need to [unwrap](https://numpy.org/doc/1.18/reference/generated/numpy.unwrap.html) the phase to get a continous function. Finally, the instantaneous frequency can be determined by\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"f(t) = \\frac{1}{2 \\pi} \\frac{\\text{d} \\phi (t)}{\\text{d} t}\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"1. Create a time series with \n",
|
||||
" 1. increasing / decreasing frequncy content, e.g. a [Chirp-Signal](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.chirp.html)\n",
|
||||
" 2. increasing and decreasing frequency content, i.e. overlapping of two chirp signals (how does the spectrogram looks like?)\n",
|
||||
"3. Write a function that computes the instantaneous frequency (Hint: use np.diff for instantaneous phase derivative, np.unwrap to unwrap the phase.)\n",
|
||||
"4. What is the meaning of $|\\hat{s}(t)|$? (Hint: plot it!) \n",
|
||||
"5. Plot the results of the instantaneous frequnecy"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Define the signal\n",
|
||||
"dt = 0.01\n",
|
||||
"t = np.arange(0, 5000) * dt\n",
|
||||
"data = chirp(t=t, f0=0.5, f1=25, t1=t[-1])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Function for instantaneous frequency\n",
|
||||
"from scipy.signal import hilbert\n",
|
||||
"def istantaneous_frequency(x, dt):\n",
|
||||
" # To continue"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Plot results"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# The Continuous Wavelet Transform (CWT)\n",
|
||||
"\n",
|
||||
"We learned during the STFT about the uncertainty problem. Using the CWT results in a different time-frequency resolution. For high frequency we get poor frequency resolution by good time resoltution, whereas for low frequencies we get poor time resolution bur good frequency resolution. As you can see, there is no recipe when to use STFT, instantaneous frequnecy or CWT, but it is good to know where are the advantages and disadvantages.\n",
|
||||
"\n",
|
||||
"The CWT is defined as\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"\\text{CWT}_{a, b} \\left(s(t) \\right) = \\frac{1}{\\sqrt{|a|}} ~ \\int_{-\\infty}^{\\infty} s(t) ~ \\Psi^\\ast \\left( \\frac{t-b}{a} \\right) \\text{d} t ~,\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"where $a$ denotes the scaling factor of the mother wavelet $\\Psi(t)$ and $b$ denotes the translation. \n",
|
||||
"In other words, the CWT is a convolution of the signal $s(t)$ with a conjugate complex and scaled wavelet-function ([here](https://en.wikipedia.org/wiki/Continuous_wavelet_transform) you find a good visualisation). Due to the different scaled wavelets, the CWT reacts well on abruptly changes in the signal.\n",
|
||||
"Usally we get instead of the frequency on the y-axis the scaling factor $a$, but it depends on the center frequency $f_c$ of the scaled wavelet by\n",
|
||||
"\n",
|
||||
"$$\n",
|
||||
"f_a = \\frac{f_c}{a} .\n",
|
||||
"$$\n",
|
||||
"\n",
|
||||
"Luckily it exists a library [pywt](https://pywavelets.readthedocs.io/en/latest/) to compute the CWT with several wavelets.\n",
|
||||
"For the following example we use a [complex Morlet wavelet](https://en.wikipedia.org/wiki/Morlet_wavelet) per default.\n",
|
||||
"\n",
|
||||
"#### Tasks\n",
|
||||
"\n",
|
||||
"1. Try to plot a time-frequency representation of different signals and compare with results of STFT\n",
|
||||
"2. Use different [mother wavelets](https://pywavelets.readthedocs.io/en/latest/ref/wavelets.html)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Display a mother wavelet\n",
|
||||
"wav, a = wavelet = pywt.ContinuousWavelet(\"cmor2-0.95\").wavefun();\n",
|
||||
"plt.plot(a, wav.real, label=\"real part\");\n",
|
||||
"plt.plot(a, wav.imag, label=\"imag part\");\n",
|
||||
"plt.legend();"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def frequency2scale(fmin, fmax, waveletname, dt):\n",
|
||||
" \"\"\"\n",
|
||||
" Computes scales for given frequency values fmin and famx.\n",
|
||||
" :param fmin: frequency minimum\n",
|
||||
" :param fmax: frequency maximum\n",
|
||||
" :param waveletname: name of wavelet\n",
|
||||
" :param dt: sampling interval (s)\n",
|
||||
" \"\"\"\n",
|
||||
"\n",
|
||||
" # Determine lower bound of scales\n",
|
||||
" scale_start = 0.1\n",
|
||||
" scale_min = pywt.scale2frequency(waveletname, scale_start)\n",
|
||||
" while scale_min/dt > fmax:\n",
|
||||
" scale_start *= 1.1\n",
|
||||
" scale_min = pywt.scale2frequency(waveletname, scale_start)\n",
|
||||
" scale_min = round(scale_start, 2)\n",
|
||||
"\n",
|
||||
" # Determine upper bound of scales\n",
|
||||
" scale_start = 1\n",
|
||||
" scale_max = pywt.scale2frequency(waveletname, scale_start)\n",
|
||||
" while scale_max/dt > fmin:\n",
|
||||
" scale_start *= 1.1\n",
|
||||
" scale_max = pywt.scale2frequency(waveletname, scale_start)\n",
|
||||
" scale_max = round(scale_start, 2)\n",
|
||||
"\n",
|
||||
" return scale_min, scale_max\n",
|
||||
"\n",
|
||||
"def scaleogram(x, fs, waveletname=\"cmor2-0.95\", fmin=1, fmax=25, num=100):\n",
|
||||
" \"\"\"\n",
|
||||
" Computes CWT of an signal x with sampling rate fs in Hz \n",
|
||||
" :param x: time seris\n",
|
||||
" :param fs: sampling rate\n",
|
||||
" :param waveletname: name of wavelet\n",
|
||||
" :param fmin: minimum frequency\n",
|
||||
" :param fmax: maximum frequency\n",
|
||||
" :param num: length of scale array\n",
|
||||
" \"\"\"\n",
|
||||
"\n",
|
||||
" # Check all input parameters and change them is they do not fit\n",
|
||||
" dt = 1 / fs\n",
|
||||
" f_ny = 1 / (2 * dt) # Nyquist frequency\n",
|
||||
" if f_ny < fmax:\n",
|
||||
" fmax = f_ny\n",
|
||||
"\n",
|
||||
" if fmin >= fmax:\n",
|
||||
" msg = \"Min. frequency {} is greater than max. frequency {} to compute scaleograms.\".format(fmin, fmax)\n",
|
||||
" raise ValueError(msg)\n",
|
||||
"\n",
|
||||
" # Determine scales from fmin and fmax\n",
|
||||
" scale_min, scale_max = frequency2scale(fmin, fmax, waveletname=waveletname, dt=dt)\n",
|
||||
" scales = np.logspace(np.log10(scale_min), np.log10(scale_max), num=num)\n",
|
||||
"\n",
|
||||
" # Generate scaleogram\n",
|
||||
" coefficients, frequencies = pywt.cwt(x, scales, waveletname, sampling_period=dt)\n",
|
||||
" coefficients = np.abs(coefficients)**2\n",
|
||||
" time = np.linspace(0, dt * len(x), num=len(x))\n",
|
||||
" \n",
|
||||
" return frequencies, time, coefficients"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"f_cwt, t_cwt, coeff = scaleogram(data, fs=1/(t[1]-t[0]))\n",
|
||||
"plt.figure(figsize=(16, 8))\n",
|
||||
"plt.pcolormesh(t_cwt, f_cwt, np.abs(coeff))\n",
|
||||
"plt.xlabel(\"Time (s)\")\n",
|
||||
"plt.xlabel(\"Frequency (Hz)\");"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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
|
||||
}
|
||||
35
06-Surface_Waves/movingWindow.py
Normal file
35
06-Surface_Waves/movingWindow.py
Normal file
@@ -0,0 +1,35 @@
|
||||
import numpy as np
|
||||
|
||||
def movingWindowAnalysis(data,winfun,nwin,shift,exp):
|
||||
"""
|
||||
Performs moving window analysis of a time series.
|
||||
data: data array
|
||||
winfun: name of the window function to be called
|
||||
nwin: number of window samples (power of 2)
|
||||
shift: displacement of moving window in samples
|
||||
exp: exponent for taking power of spectrum
|
||||
"""
|
||||
fwin = winfun(nwin) # compute window values
|
||||
npts = len(data) # number of total samples
|
||||
nseg = int((npts-nwin)/shift)+1 # total number of expected data segment
|
||||
mwa = np.zeros((nwin//2+1,nseg)) # array for result (rfft returns N/2+1 samples)
|
||||
wa = 0 # start index of data segment
|
||||
we = nwin # end index of data segment
|
||||
jseg = 0 # initialize data segment counter
|
||||
while we < npts: # loop over segments
|
||||
seg = data[wa:we]*fwin # multiply data segment with window
|
||||
seg = seg-seg.mean() # subtract mean value of segment
|
||||
ftseg = np.abs(np.fft.rfft(seg)) # abs value of Fourier transform
|
||||
maxft = np.amax(ftseg) # max value of Fourier transform
|
||||
ftseg = ftseg/maxft+1.e-10 # normalize spectrum to its maximum, remove zeros
|
||||
mwa[:,jseg] = np.power(ftseg,exp) # assign values to the matrix
|
||||
wa = wa+shift # move window start by shift
|
||||
we = we+shift # move window end by shift
|
||||
jseg = jseg+1 # increase segment counter
|
||||
return nseg,mwa # return number of segments and moving window matrix
|
||||
#------------------------------------------------------------------------------------
|
||||
#
|
||||
def hann(nw):
|
||||
arg = 2.*np.pi*np.arange(0,nw)/nw # argument of cosine
|
||||
fwin = 0.5*(1.-np.cos(arg)) # Hann window
|
||||
return fwin
|
||||
47
06-Surface_Waves/multipleFilter.py
Normal file
47
06-Surface_Waves/multipleFilter.py
Normal file
@@ -0,0 +1,47 @@
|
||||
import numpy as np
|
||||
|
||||
def multipleFilterAnalysis(data,alfa,cfreq,dt,ndec):
|
||||
"""
|
||||
Perform a multiple filter analysis of data.
|
||||
data: Array of detrended and demeaned data whose length is power of 2
|
||||
alfa: Width parameter of Gaussian bandpass filter
|
||||
cfreq: Array of center frequencies of Gaussian filter
|
||||
dt: sampling interval of data
|
||||
ndec: decimation factor for instantaneous amplitude output
|
||||
"""
|
||||
npts = len(data)
|
||||
nd = int(pow(2, np.ceil(np.log(npts)/np.log(2)))) # find next higher power of 2 of npts
|
||||
ftd = np.fft.rfft(data,nd) # Fourier transform of entire data set (pos. freq.)
|
||||
# data are padded with zeros since npts <= nd
|
||||
freq = np.fft.rfftfreq(nd,dt) # Fourier frequencies (positive frequencies)
|
||||
nt = int(np.ceil(npts/ndec))
|
||||
mfa = np.zeros((len(cfreq),nt)) # numpy array for MFA result
|
||||
for jf,cf in enumerate(cfreq):
|
||||
hg = np.exp(-alfa*((freq-cf)/cf)**2) # Gaussian filter (use f instead of omega here)
|
||||
fk = hg*ftd # multiply FT of data with filter
|
||||
qk = np.complex(0,1)*fk # FT of Hilbert transform
|
||||
ftk = np.fft.irfft(fk) # filtered data
|
||||
qtk = np.fft.irfft(qk) # Hilbert transform of filtered data
|
||||
at = np.sqrt(ftk**2+qtk**2) # instantaneous amplitude
|
||||
mfa[jf,:] = at[0:npts:ndec] # store decimated original result
|
||||
return mfa
|
||||
#-----------------------------------------------------------------------
|
||||
# normalize multiple filter result either along time or frequency axis
|
||||
def normalizeMFT(mfa,mode,exp):
|
||||
"""
|
||||
Normalize the result of the mutiple filtering operation.
|
||||
mfa: array with instantaneous amplitudes versus frequency (mfa(f,t))
|
||||
mode: normalization mode: if 'time', normalize along time axis
|
||||
else normalize along frequency axis
|
||||
exp: exponent for modifying inst amp using a power less than 1
|
||||
"""
|
||||
nf,nt = mfa.shape
|
||||
if mode == 'time':
|
||||
for jf in range(0,nf):
|
||||
mfamax = np.amax(mfa[jf,:])
|
||||
mfa[jf,:] = np.power(mfa[jf,:]/mfamax+1.e-10,exp)
|
||||
else:
|
||||
for jt in range(0,nt):
|
||||
mfamax = np.amax(mfa[:,jt])
|
||||
mfa[:,jt] = np.power(mfa[:,jt]/mfamax+1.e-10,exp)
|
||||
return mfa
|
||||
268
06-Surface_Waves/phaseVelocityCrossCorr.ipynb
Normal file
268
06-Surface_Waves/phaseVelocityCrossCorr.ipynb
Normal file
@@ -0,0 +1,268 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Two-station method\n",
|
||||
"Here we compute the phase velocity form seismograms of two stations lying on a common great circle. We make use of the fact that the phase of the Fourier tranform of the cross-correlation of the two seismograms is equal to the phase difference between the surface wave trains from which we calclate the phase velocity."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Preparation: load packages, set some basic options \n",
|
||||
"%matplotlib inline\n",
|
||||
"from obspy import UTCDateTime\n",
|
||||
"from obspy.clients.fdsn import Client\n",
|
||||
"from obspy.geodetics import locations2degrees\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pylab as plt\n",
|
||||
"from multipleFilter import multipleFilterAnalysis, normalizeMFT \n",
|
||||
"plt.style.use('ggplot')\n",
|
||||
"plt.rcParams['figure.figsize'] = 15, 4\n",
|
||||
"plt.rcParams['lines.linewidth'] = 0.5"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Get data from earthquake (Gulf of Alaska, m=7.9, d = 24 km)\n",
|
||||
"client = Client(\"BGR\")\n",
|
||||
"t1 = UTCDateTime(\"2018-01-23T09:31:42.000\")\n",
|
||||
"st1raw = client.get_waveforms(\"GR\",\"TNS\",\"\",\"LHZ\",t1,t1 + 2 * 3600, # data of station TNS of GRSN\n",
|
||||
" attach_response = True)\n",
|
||||
"inv1 = client.get_stations(network=\"GR\",station=\"TNS\",\n",
|
||||
" channel='LHZ',level=\"channel\",starttime=t1) # inventory info about TNS\n",
|
||||
"client = Client(\"GFZ\")\n",
|
||||
"st2raw = client.get_waveforms(\"GE\",\"STU\",\"\",\"LHZ\",t1,t1 + 2 * 3600, # data of station STU of GEOFON\n",
|
||||
" attach_response = True)\n",
|
||||
"inv2 = client.get_stations(network=\"GE\",station=\"STU\",\n",
|
||||
" channel='LHZ',level=\"channel\",starttime=t1) # inventory info about STU"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"st1 = st1raw.copy() # copy data\n",
|
||||
"st2 = st2raw.copy()\n",
|
||||
"st1.detrend('linear') # detrend\n",
|
||||
"st2.detrend('linear')\n",
|
||||
"st1.remove_response(output='VEL',zero_mean = True, # remove response (pre-filter is important)\n",
|
||||
" pre_filt = (0.001, 0.005, 0.4, 0.5))\n",
|
||||
"st2.remove_response(output='VEL',zero_mean = True, # remove response(pre-filter is important)\n",
|
||||
" pre_filt = (0.001, 0.005, 0.4, 0.5))\n",
|
||||
"st1.plot() # plot data\n",
|
||||
"st2.plot()\n",
|
||||
"print(st1)\n",
|
||||
"print(st2)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"co = inv1.get_coordinates(\"GR.TNS..LHZ\") # extract coordinates from inventory\n",
|
||||
"tns = (co[\"longitude\"],co[\"latitude\"])\n",
|
||||
"co = inv2.get_coordinates(\"GE.STU..LHZ\")\n",
|
||||
"stu = (co[\"longitude\"],co[\"latitude\"])\n",
|
||||
"event = (-149.09,56.05) # coordinates of event\n",
|
||||
"dis = locations2degrees(tns[1],tns[0],stu[1],stu[0])*np.pi/180.*6371. # epicentral distance\n",
|
||||
"print(\"Distance in km: \",dis)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"<img src=\"../images/phase_velocity_map.png\">"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"tw1 = UTCDateTime(\"2018-01-23T10:05:00.000\") # cut seismograms to surface wave train\n",
|
||||
"tw2 = UTCDateTime(\"2018-01-23T10:16:00.000\")\n",
|
||||
"st1.trim(tw1,tw2)\n",
|
||||
"st2.trim(tw1,tw2)\n",
|
||||
"st1.plot()\n",
|
||||
"st2.plot()\n",
|
||||
"npts = st1[0].stats.npts\n",
|
||||
"dt = st1[0].stats.delta\n",
|
||||
"print(\"Number of samples: \",npts,\"Sampling interval: \",dt)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"cc = np.correlate(st2[0],st1[0],mode=\"full\") # Cross correlation: cc(k)=sum_n a(n+k)b(n)\n",
|
||||
"tm = np.arange(-(npts-1),npts)*dt # includes negative lags as well\n",
|
||||
"plt.plot(tm,cc,'-k')\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
},
|
||||
"scrolled": true
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"ns = 256 # limit length of cross correlation to first ns points\n",
|
||||
"ccp = cc[npts-1:npts-1+ns] # cross correlation for positive lag\n",
|
||||
"print(len(ccp))\n",
|
||||
"tm = np.arange(0,ns)*dt # associated time values\n",
|
||||
"plt.plot(tm,ccp,'-k')\n",
|
||||
"plt.xlabel(\"Time\")\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"nd = int(pow(2, np.ceil(np.log(ns)/np.log(2)))) # Fourier transfrom of the cross correlation\n",
|
||||
"fcp = np.fft.rfft(ccp,nd)\n",
|
||||
"freq = np.fft.rfftfreq(nd,dt)\n",
|
||||
"print(nd,dt,len(freq),freq[0:4])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"fmax = 0.1 # compute phase up to maximum frequency fmax\n",
|
||||
"fny = 1./(2.*dt)\n",
|
||||
"jmax = int(fmax/fny*(nd//2-1))\n",
|
||||
"print(jmax)\n",
|
||||
"spec = fcp[0:jmax]\n",
|
||||
"frs = freq[0:jmax]\n",
|
||||
"phase = np.unwrap(np.angle(spec)) # unwrap the phase owing to 2*pi multiples"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"plt.subplot(311)\n",
|
||||
"plt.plot(frs,np.abs(spec),\"-k\") # plot abs of spectrum\n",
|
||||
"plt.subplot(312)\n",
|
||||
"plt.plot(frs,np.angle(spec),\"-k\") # plot phase of spectrum\n",
|
||||
"plt.subplot(313)\n",
|
||||
"plt.plot(frs,phase,\"-k\") # plot unwrapped phase\n",
|
||||
"plt.xlabel(\"Frequency [Hz]\")\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"j1 = int(jmax*0.1)\n",
|
||||
"j2 = int(jmax*0.6)\n",
|
||||
"phavel = -2.*np.pi*frs[j1:j2]*dis/phase[j1:j2] # calculate phase velocity and plot\n",
|
||||
"plt.plot(frs[j1:j2],phavel,'-k')\n",
|
||||
"plt.xlabel(\"Frequency [Hz]\")\n",
|
||||
"plt.ylabel(\"Phase velocity [km/s]\")\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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": 4
|
||||
}
|
||||
256
06-Surface_Waves/surfaceWaveGroupVelocity.ipynb
Normal file
256
06-Surface_Waves/surfaceWaveGroupVelocity.ipynb
Normal file
@@ -0,0 +1,256 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Surface Wave Velocity (Single Station Approach)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Preparation: load packages, set some basic options \n",
|
||||
"%matplotlib inline\n",
|
||||
"from obspy import UTCDateTime\n",
|
||||
"from obspy.clients.fdsn import Client\n",
|
||||
"from obspy.clients.fdsn import URL_MAPPINGS\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pylab as plt\n",
|
||||
"from movingWindow import movingWindowAnalysis, hann # import the moving window analysis routines\n",
|
||||
"from multipleFilter import multipleFilterAnalysis, normalizeMFT # import the multiple filter analysis routines\n",
|
||||
"plt.style.use('ggplot')\n",
|
||||
"plt.rcParams['figure.figsize'] = 15, 4\n",
|
||||
"plt.rcParams['lines.linewidth'] = 0.5"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Data from data archive using obspy webfdsn client\n",
|
||||
"FDSN stands for the International Federal of Digital Seismic Networks (www.fdsn.org). Citation from their web site: \"The International Federation of Digital Seismograph Networks (FDSN) is a global organization. Its membership is comprised of groups responsible for the installation and maintenance of seismographs either within their geographic borders or globally. Membership in the FDSN is open to all organizations that operate more than one broadband station. Members agree to coordinate station siting and provide free and open access to their data. This cooperation helps scientists all over the world to further the advancement of earth science and particularly the study of global seismic activity.\"\n",
|
||||
"BGR (Bundesanstalt für Geowissenschaften und Rohstoffe) operates broadband stations in Germany and operates a data archive. It is member of the FDSN. That is why we can freely download data from them.\n",
|
||||
"There are other data archives around the world which are also member of the FDSN and allow opn access to their data (see OBSPY documentation "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"for key in sorted(URL_MAPPINGS.keys()): # eine Liste der Archive\n",
|
||||
" print(\"{0:<7} {1}\".format(key, URL_MAPPINGS[key]))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Get waveforms for the Tohoku earthquake for station TNS of the German Regional Seismic Network."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"client = Client(\"BGR\") # data archive from which data are fetched\n",
|
||||
"t1 = UTCDateTime(\"2018-01-14T09:19:00.000\") # desired start time of data\n",
|
||||
"stc = client.get_waveforms(\"GR\", \"TNS\", \"\", \"LHZ\", t1, # use fdsn web client to download data\n",
|
||||
" t1 + 7200., attach_response=True)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"stc.detrend('linear') # take out linear trend\n",
|
||||
"stc.plot()\n",
|
||||
"print(stc)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"st = stc.copy() # proceed with copy to avoid repeated downloading\n",
|
||||
"ta = UTCDateTime(\"2018-01-14T10:05:00.000000Z\") # cut to surface wave train\n",
|
||||
"te = UTCDateTime(\"2018-01-14T10:22:00.000000Z\")\n",
|
||||
"st.trim(ta,te)\n",
|
||||
"st.taper(0.05,\"hann\") # apply a taper to bring signal to zero at ends\n",
|
||||
"st.plot() # dispersion is now well visible\n",
|
||||
"print(st)\n",
|
||||
"tr = st[0]"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Moving window approach"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"wlen = 256.0 # window length in sec\n",
|
||||
"npts = tr.stats.npts # number of samples in the trace\n",
|
||||
"dt = tr.stats.delta # sample interval\n",
|
||||
"nwin = int(wlen/dt) # nuber of samples in window\n",
|
||||
"nwin = int(pow(2, np.ceil(np.log(nwin)/np.log(2)))) # find next higher power of 2 of nwin\n",
|
||||
"shift = nwin//8\n",
|
||||
"tshift = shift*dt\n",
|
||||
"nseg,mwa = movingWindowAnalysis(tr.data,hann,nwin,shift,2.0) # compute spectrogram, exponent=2.0\n",
|
||||
"freq = np.fft.rfftfreq(nwin,dt) # Fourier frequencies\n",
|
||||
"print(\"Frequency range [Hz]: \",freq[0],freq[-1])\n",
|
||||
"print(\"Number of frequency samples: \",len(freq))\n",
|
||||
"print(\"Window length in samples: \",nwin)\n",
|
||||
"print(\"Window shift in samples: \",shift)\n",
|
||||
"print(\"Number of segments: \",nseg)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# plot the result\n",
|
||||
"st.plot() # plot seismogram again\n",
|
||||
"fmax = 0.08 # maximum frequency to be plotted\n",
|
||||
"fmin = 0.02\n",
|
||||
"nf = len(freq)-1 # number of frequencies\n",
|
||||
"jfmax = int(fmax/freq[-1]*nf) # max index for plotting \n",
|
||||
"jfmin = int(fmin/freq[-1]*nf) # min index for plotting \n",
|
||||
"extent = (0.5*wlen,0.5*wlen+nseg*tshift,fmin,fmax) # extent of matrix in true time and frequency\n",
|
||||
"asp = 0.5*nseg*tshift/(fmax-fmin)\n",
|
||||
"fg = plt.figure(figsize = [12.5,12])\n",
|
||||
"plt.imshow(mwa[jfmin:jfmax,:],extent = extent, aspect = asp, \n",
|
||||
" origin = \"lower\", interpolation = \"bicubic\") # plotting with interpolation\n",
|
||||
"plt.xlabel('Window center time [s]')\n",
|
||||
"plt.ylabel('Frequency [Hz]')\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"Observe the change of frequency content with time from 0.03 Hz initially to 0.06 Hz towards the end of the record. Given a epicentral distance $d$, we could transform the time axis into a group velocity axis via $U=d/t$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Multiple filtering approach"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"alfa = 50.0 # filter width parameter\n",
|
||||
"cfreq = np.linspace(fmin,fmax,100) # array of filter center frequencies\n",
|
||||
"mfa = multipleFilterAnalysis(tr.data,alfa,cfreq,dt,1)\n",
|
||||
"mfa = normalizeMFT(mfa,'freq',2.0)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"jupyter": {
|
||||
"outputs_hidden": false
|
||||
}
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"st.plot()\n",
|
||||
"fg = plt.figure(figsize = [12.5,12])\n",
|
||||
"extent = (0.,npts*dt,fmin,fmax) # extent of matrix in true time and frequency\n",
|
||||
"asp = 0.5*npts*dt/(fmax-fmin)\n",
|
||||
"plt.imshow(mfa,extent = extent, aspect = asp, \n",
|
||||
" origin = \"lower\", interpolation = \"bicubic\") # plotting\n",
|
||||
"plt.xlabel('Time [s]')\n",
|
||||
"plt.ylabel('Frequency [Hz]')\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"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": 4
|
||||
}
|
||||
BIN
07-Receiver_Functions/data/PB01.BHE.20110225.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHE.20110225.sac
Normal file
Binary file not shown.
BIN
07-Receiver_Functions/data/PB01.BHE.20110306.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHE.20110306.sac
Normal file
Binary file not shown.
BIN
07-Receiver_Functions/data/PB01.BHE.20110513.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHE.20110513.sac
Normal file
Binary file not shown.
BIN
07-Receiver_Functions/data/PB01.BHN.20110225.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHN.20110225.sac
Normal file
Binary file not shown.
BIN
07-Receiver_Functions/data/PB01.BHN.20110306.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHN.20110306.sac
Normal file
Binary file not shown.
BIN
07-Receiver_Functions/data/PB01.BHN.20110513.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHN.20110513.sac
Normal file
Binary file not shown.
BIN
07-Receiver_Functions/data/PB01.BHZ.20110225.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHZ.20110225.sac
Normal file
Binary file not shown.
BIN
07-Receiver_Functions/data/PB01.BHZ.20110306.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHZ.20110306.sac
Normal file
Binary file not shown.
BIN
07-Receiver_Functions/data/PB01.BHZ.20110513.sac
Normal file
BIN
07-Receiver_Functions/data/PB01.BHZ.20110513.sac
Normal file
Binary file not shown.
473
07-Receiver_Functions/example_events.xml
Normal file
473
07-Receiver_Functions/example_events.xml
Normal file
@@ -0,0 +1,473 @@
|
||||
<?xml version='1.0' encoding='utf-8'?>
|
||||
<ns0:quakeml xmlns:ns0="http://quakeml.org/xmlns/quakeml/1.2" xmlns="http://quakeml.org/xmlns/bed/1.2">
|
||||
<eventParameters publicID="smi:service.iris.edu/fdsnws/event/1/query">
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3287729">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=10171447</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16917168</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>CENTRAL MID-ATLANTIC RIDGE</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=10171447">
|
||||
<time>
|
||||
<value>2011-05-15T13:08:15.420000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>0.4584</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-25.6088</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>18900.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16917168">
|
||||
<mag>
|
||||
<value>6.1</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=10171449</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3287620">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=10167818</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16912325</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>COSTA RICA</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=10167818">
|
||||
<time>
|
||||
<value>2011-05-13T22:47:55.340000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>10.1114</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-84.1889</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>76800.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16912325">
|
||||
<mag>
|
||||
<value>6.0</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=10167822</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3285786">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=10137185</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16870909</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>SOUTH OF PANAMA</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=10137185">
|
||||
<time>
|
||||
<value>2011-04-30T08:19:16.720000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>6.8511</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-82.3594</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>10000.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16870909">
|
||||
<mag>
|
||||
<value>6.2</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=10137190</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3284483">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=10109851</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16833821</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>SOUTH OF KERMADEC ISLANDS</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=10109851">
|
||||
<time>
|
||||
<value>2011-04-18T13:03:04.360000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>-34.286</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>179.9433</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>98100.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16833821">
|
||||
<mag>
|
||||
<value>6.5</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=10109854</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3282641">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=10082429</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16795144</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>CHIAPAS, MEXICO</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=10082429">
|
||||
<time>
|
||||
<value>2011-04-07T13:11:23.430000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>17.2651</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-94.1439</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>165100.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16795144">
|
||||
<mag>
|
||||
<value>6.7</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=10082436</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3281051">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=10053039</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16759916</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>FIJI ISLANDS REGION</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=10053039">
|
||||
<time>
|
||||
<value>2011-03-31T00:11:58.880000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>-16.5479</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-177.3915</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>19400.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16759916">
|
||||
<mag>
|
||||
<value>6.4</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=10053046</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3279149">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=9925460</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16631835</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>SOUTH SANDWICH ISLANDS REGION</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=9925460">
|
||||
<time>
|
||||
<value>2011-03-06T14:32:36.940000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>-56.3864</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-27.0253</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>92000.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16631835">
|
||||
<mag>
|
||||
<value>6.5</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=9925462</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3278515">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=9916282</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16620185</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>EASTER ISLAND REGION</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=9916282">
|
||||
<time>
|
||||
<value>2011-03-01T00:53:45.350000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>-29.6428</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-112.1246</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>3800.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16620185">
|
||||
<mag>
|
||||
<value>6.1</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=9916288</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3278477">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=9851774</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16541902</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>OAXACA, MEXICO</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=9851774">
|
||||
<time>
|
||||
<value>2011-02-25T13:07:26.980000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>17.8214</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-95.1708</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>130600.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16541902">
|
||||
<mag>
|
||||
<value>6.0</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=9851779</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3278416">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=9831995</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16530631</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>SOUTH ISLAND, NEW ZEALAND</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=9831995">
|
||||
<time>
|
||||
<value>2011-02-21T23:51:42.340000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>-43.4935</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>172.713</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>4800.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16530631">
|
||||
<mag>
|
||||
<value>6.1</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=9832004</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3278381">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=9831146</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16529492</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>SOUTH OF FIJI ISLANDS</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=9831146">
|
||||
<time>
|
||||
<value>2011-02-21T10:57:51.760000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>-26.0435</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>178.4765</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>551800.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16529492">
|
||||
<mag>
|
||||
<value>6.5</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=9831149</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3277925">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=9817429</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16511317</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>TONGA ISLANDS</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=9817429">
|
||||
<time>
|
||||
<value>2011-02-12T17:57:56.170000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>-20.8515</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-175.5845</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>85900.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16511317">
|
||||
<mag>
|
||||
<value>6.1</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=9817432</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
<event publicID="smi:service.iris.edu/fdsnws/event/1/query?eventid=3277104">
|
||||
<preferredOriginID>smi:service.iris.edu/fdsnws/event/1/query?originid=9794081</preferredOriginID>
|
||||
<preferredMagnitudeID>smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16481049</preferredMagnitudeID>
|
||||
<type>earthquake</type>
|
||||
<description>
|
||||
<text>TONGA ISLANDS</text>
|
||||
<type>Flinn-Engdahl region</type>
|
||||
</description>
|
||||
<origin publicID="smi:service.iris.edu/fdsnws/event/1/query?originid=9794081">
|
||||
<time>
|
||||
<value>2011-01-31T06:03:26.330000Z</value>
|
||||
</time>
|
||||
<latitude>
|
||||
<value>-21.9987</value>
|
||||
</latitude>
|
||||
<longitude>
|
||||
<value>-175.5367</value>
|
||||
</longitude>
|
||||
<depth>
|
||||
<value>69300.0</value>
|
||||
</depth>
|
||||
<creationInfo>
|
||||
<author>ISC</author>
|
||||
</creationInfo>
|
||||
</origin>
|
||||
<magnitude publicID="smi:service.iris.edu/fdsnws/event/1/query?magnitudeid=16481049">
|
||||
<mag>
|
||||
<value>6.0</value>
|
||||
</mag>
|
||||
<type>MW</type>
|
||||
<originID>smi:service.iris.edu/fdsnws/event/1/query?originid=9794089</originID>
|
||||
<creationInfo>
|
||||
<author>GCMT</author>
|
||||
</creationInfo>
|
||||
</magnitude>
|
||||
</event>
|
||||
</eventParameters>
|
||||
</ns0:quakeml>
|
||||
654
07-Receiver_Functions/receiverFuntionCalculation.ipynb
Normal file
654
07-Receiver_Functions/receiverFuntionCalculation.ipynb
Normal file
File diff suppressed because one or more lines are too long
18
README.md
18
README.md
@@ -18,3 +18,21 @@ This has been forked from [Seismo Live](http://seismo-live.org). The source code
|
||||
|
||||
1. Simple Cross Correlation: Use the cross-correlation to detect and determine the time shift of two test events with one template event.
|
||||
2. Enhanced Picker: This has been taken from the [ObsPy Tutorial](https://docs.obspy.org/master/tutorial/code_snippets/xcorr_pick_correction.html). ObsPy is licensed under the [LGPL v3.0](https://www.gnu.de/documents/lgpl-3.0.en.html)
|
||||
3. Ambient Seismic Noise Analysis from [Seismo Live](http://seismo-live.org). The source code is available at https://github.com/krischer/seismo_live (licensed under a ***CC BY-NC-SA 4.0 License***. © 2015-2019 Lion Krischer).
|
||||
|
||||
### 04 - FFT, DFT and Applications
|
||||
|
||||
Interpolation and Deciamtion of Time Series using the DFT
|
||||
|
||||
### 05 - Spectrogram
|
||||
|
||||
1. Analysing dispersive signals. How to visualize changes in frequency content over time (moving window analysis, mutiple filter technique, instantaneous frequency, continuous wavelet transform, the uncertainty problem).
|
||||
2. Filtering time series with non spectral filtering method.
|
||||
|
||||
### 06 - Surface Waves
|
||||
|
||||
Analysing phase and group velocity of recorded earthquake signals at one and two seismological stations.
|
||||
|
||||
### 07 - Receiver Functions
|
||||
|
||||
How to get receiver functions from seismograms.
|
||||
BIN
images/karte3.jpg
Normal file
BIN
images/karte3.jpg
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 749 KiB |
BIN
images/phase_velocity_map.png
Normal file
BIN
images/phase_velocity_map.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 22 KiB |
BIN
images/seismogram.png
Normal file
BIN
images/seismogram.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 127 KiB |
BIN
images/stack-surface-waves.png
Normal file
BIN
images/stack-surface-waves.png
Normal file
Binary file not shown.
|
After Width: | Height: | Size: 28 MiB |
Reference in New Issue
Block a user