new assignment 7
This commit is contained in:
parent
9327a231a9
commit
d0ec9f7e40
4
assignment_07/CH.PANIX.HHZ
Normal file
4
assignment_07/CH.PANIX.HHZ
Normal file
File diff suppressed because one or more lines are too long
4
assignment_07/Z3.A261A.HHZ
Normal file
4
assignment_07/Z3.A261A.HHZ
Normal file
File diff suppressed because one or more lines are too long
27
assignment_07/dftFast.py
Normal file
27
assignment_07/dftFast.py
Normal file
@ -0,0 +1,27 @@
|
||||
# -------------------------------------------------------------
|
||||
# Functions for "fast" Fourier transform and Fourier synthesis
|
||||
# using numpy tools
|
||||
# -------------------------------------------------------------
|
||||
import numpy as np
|
||||
|
||||
def dft_fast_coeff(x):
|
||||
"""
|
||||
Evaluate Fourier coefficients using Numpy's fast Fourier transform.
|
||||
This routine only returns the coefficients for the positive frequencies.
|
||||
If N is even, it goes up to n=N/2.
|
||||
If N is odd, it goes up to n=(N-1)/2.
|
||||
:param x: array of function samples
|
||||
"""
|
||||
return np.fft.rfft(x, norm='forward')
|
||||
|
||||
|
||||
def dft_fast_synthesis(fc, outnum='even'):
|
||||
"""
|
||||
Use numpy's fast Fourier synthesis taking only Fourier coefficients for positive frequencies
|
||||
:param fc: aray of coefficients for positive frequencies only.
|
||||
:param outnum: specifies if output time series has an even or odd number of samples (default: 'even')
|
||||
"""
|
||||
ns = 2*fc.size-2
|
||||
if outnum == 'odd': ns = 2*fc.size-1
|
||||
return np.fft.irfft(fc, ns, norm='forward')
|
||||
|
||||
308
assignment_07/p7_seismometer_correlation.ipynb
Normal file
308
assignment_07/p7_seismometer_correlation.ipynb
Normal file
@ -0,0 +1,308 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "17de1a34-2933-47ed-b1c5-ae9dd0466edd",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"## Assignment 7"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "681b47c6-f9af-4877-b745-febd443e1877",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"### Seismometer transfer function and correlation"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "85cb9941-eb24-4eb2-9c46-8610844e4063",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%matplotlib inline\n",
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"from setupFigure import SetupFigure\n",
|
||||
"from dftFast import dft_fast_coeff, dft_fast_synthesis\n",
|
||||
"from seismicTrace import SeismicTrace"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "a627c30b-eb67-4319-930e-7f1acb8e3ee1",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Task 1: Write a function that calculates the transfer function of a seismometer for different damping $h=\\sin\\phi$."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "8cbc4a71-fc18-4c3e-a26a-0c02ae89314d",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def transfer_seismometer_velocity(fc, df, nf, phi, sigma):\n",
|
||||
" \"\"\"\n",
|
||||
" Calculates samples of the ground velocity seismometer transfer function (positive frequencies only)\n",
|
||||
" :param fc: eigenfrequency of seismometer in Hz\n",
|
||||
" :param df: frequency stepping\n",
|
||||
" :param nf: number of frequencies\n",
|
||||
" :param phi: damping angle in degrees (h=sin(phi))\n",
|
||||
" :param sigma: generator constant\n",
|
||||
" :return: array of frequencies\n",
|
||||
" :return: array of complex values of transfer function\n",
|
||||
" \"\"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "7f7d857c-1bfd-4703-96a7-61c7af95bbdc",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### You may use this function for the Butterworth lowpass filter (or your own from Assignment 6)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "dd1b00e9-a953-4d64-b91f-7d7ad6949a3f",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def butterworth_lowpass(fc, df, nf, order):\n",
|
||||
" \"\"\"\n",
|
||||
" Calculates samples of the lowpass Butterworth filter transfer function (positive frequencies only)\n",
|
||||
" :param fc: corner frequency in Hz\n",
|
||||
" :param df: frequency stepping in Hz\n",
|
||||
" :param nf: number of frequencies\n",
|
||||
" :param order: filter order\n",
|
||||
" \"\"\"\n",
|
||||
" f = df*np.arange(0, nf)\n",
|
||||
" h = np.ones(nf, dtype=np.complex64)\n",
|
||||
" for k in range(order):\n",
|
||||
" arg = 1j*(2*k+1)*np.pi/(2*order)\n",
|
||||
" h = h*(-1j)/(f/fc-np.exp(arg))\n",
|
||||
" return h"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "38d6be7b-6579-4726-abee-4666e6910f39",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### This is a useful function to plot traces next to each other with a vertical offset"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "48a7e17e-eaf2-4162-89f2-0d81eb9b7631",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def plotTraceGather(tracelist, ax):\n",
|
||||
" \"\"\"\n",
|
||||
" Plots a gather of traces with vertical offset\n",
|
||||
" :param tracelist: list of SeismicTrace objects\n",
|
||||
" :param ax: Matplotlib axis object\n",
|
||||
" \"\"\"\n",
|
||||
" for j, tr in enumerate(tracelist):\n",
|
||||
" t = tr.tanf+tr.dt*np.arange(0, tr.nsamp)\n",
|
||||
" norm = np.max(np.abs(tr.vals))\n",
|
||||
" ax.plot(t, 0.5*tr.vals/norm+j, color='b', ls='-') # Normalization to 0.5*absmax and shift by multiples of 1\n",
|
||||
" ax.text(tr.tanf, j+0.1, tr.staname, fontsize=14) # Add station name on top of trace"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "68bf389c-4464-43dc-853d-ff28e5791f99",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Task 2.1: Write a function that computes the cross correlation function for a given range of time delays"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "47a92e2f-da8c-45ce-b58c-7b8d153e38f5",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def crossCorrelation(x, y, dt, neglim, poslim):\n",
|
||||
" \"\"\"\n",
|
||||
" Correlate x with y according to c_n = Delta t*sum_(k=0)^(N-1) x_(n+k) y_k.\n",
|
||||
" Assumes that both traces have same DT!\n",
|
||||
" Consider x is a Gaussian with max at t2, y a Gaussian with max at t1 < t2.\n",
|
||||
" To bring x and y into a match, x neeeds to be shifted to the left implying a positive n=n_s.\n",
|
||||
" Thus, ns*dt is the delay of x relative to y or the advance of y relative to x.\n",
|
||||
" If n+k < 0 or n+k >= x.size, x_(n+k) is considered as zero.\n",
|
||||
" :param x: values of first trace\n",
|
||||
" :param y: values of second trace\n",
|
||||
" :param dt: sampling interval\n",
|
||||
" :param neglim: most negative value for n\n",
|
||||
" :param poslim: most positive value for n (poslim not included)\n",
|
||||
" :return: array of delay times associated with the array elements of the cross correlation\n",
|
||||
" :return: array of values of the cross correlation function\n",
|
||||
" \"\"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "79f44179-4943-4e87-9505-4c326c4237c2",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Task 2.2: Evaluate the amplitude spectrum of the seismometer transfer function for different angles\n",
|
||||
"Use fc = 10 Hz, nf = 256, df = 0.2, sigma = 1. Choose angles 15, 25, 45, 60, 85 degrees. Plot the\n",
|
||||
"amplitude spectra into one graph."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "56a730c7-ba3c-4404-a8e4-6d8c83d5fab3",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\"\"\"Your code\"\"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "b2576de7-c170-4829-a5a9-a107e6ad09bb",
|
||||
"metadata": {
|
||||
"jp-MarkdownHeadingCollapsed": true
|
||||
},
|
||||
"source": [
|
||||
"#### Task 2.3: Evaluate the phase spectrum of the seismometer transfer function for different angles\n",
|
||||
"Use fc = 10 Hz, nf = 256, df = 0.2, sigma = 1. Choose angles 15, 25, 45, 60, 85 degrees. Plot the\n",
|
||||
"amplitude spectra into one graph."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "a796b315-f5e4-48cc-ae2e-d98323cc4996",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\"\"\"Your code\"\"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "d220be85-9530-418d-9069-6ad21d87872a",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Task 2.4: Setup two Gaussians $e^{-(t-\\mu)^2/\\sigma^2}$ of lengths 100 s and 80 s with dt=1 s, $\\mu$=55 and 45, and $\\sigma$ = 5. Plot the two into one graph. Then, compute the cross correlation between -30 and +30 s using your function from Task 2.1 and plot it into a separate graph. "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "f7537491-54c9-4f9c-ac22-e02d9cefae05",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\"\"\"Your code\"\"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "1529c72b-6ae3-416a-b4b2-ed1fd0bffba3",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Task 2.5: For comparison, use numpy's routine correlate to calculate the correlation. Try the different modes offered. Change the lengths of x and y."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "da38f6a3-4d2e-4a69-a605-7de5a55f5cfd",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\"\"\"Your code\"\"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "5226e897-3874-45ba-86ab-883dd92fc13e",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Task 3.1: We want to prepare tome delay measurements between two seismograms by cross-correlation. First, read the data of stations CH.PANIX and Z3.A261A.HHZ. Put the SeismicTrace objects into a list."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "0655529e-3212-4069-a60f-bef57644328c",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\"\"\"Your code\"\"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "9f38af03-65f9-49d8-a28d-5fc0e2c4d965",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Task 3.2: Plot the traces using the function plotTraceGather."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "56551bfa-e2da-4773-9247-a33cae96a3c1",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\"\"\"Your code\"\"\""
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"id": "cec5430e-c602-4a99-a7c3-f3371e990d82",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"#### Task 3.3: Filter the data using a Butterworth low pass with fc=0.5 Hz and order 6. Use either the provided function for the low-pass transfer function or take your own. First, Fourier transform the data, then multiply with the filter transfer function, and finally do the back transform. Create new SeismicTrace objects for the filtered traces, put them into a list and plot them."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"id": "6ce5cc83-f5ff-4fb2-a9ac-d98595429c02",
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"\"\"\"Your code\"\"\""
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"kernelspec": {
|
||||
"display_name": "Python 3 (ipykernel)",
|
||||
"language": "python",
|
||||
"name": "python3"
|
||||
},
|
||||
"language_info": {
|
||||
"codemirror_mode": {
|
||||
"name": "ipython",
|
||||
"version": 3
|
||||
},
|
||||
"file_extension": ".py",
|
||||
"mimetype": "text/x-python",
|
||||
"name": "python",
|
||||
"nbconvert_exporter": "python",
|
||||
"pygments_lexer": "ipython3",
|
||||
"version": "3.12.4"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 5
|
||||
}
|
||||
79
assignment_07/seismicTrace.py
Normal file
79
assignment_07/seismicTrace.py
Normal file
@ -0,0 +1,79 @@
|
||||
#----------------------------------------------------------------------------
|
||||
# Python class for a seismic trace
|
||||
#
|
||||
import numpy as np
|
||||
|
||||
class SeismicTrace(object):
|
||||
"""
|
||||
Class defining a seismic trace by station name, network name, component,
|
||||
number of samples, sampling interval and start time
|
||||
"""
|
||||
|
||||
def __init__(self, staname, comp, dt, tanf, vals):
|
||||
"""
|
||||
Initialize seismic trace object from data:
|
||||
:param staname: station name
|
||||
:param comp: component
|
||||
:param dt: sampling interval
|
||||
:param tanf: start time
|
||||
:param vals: values of trace (1D numpy array)
|
||||
"""
|
||||
self.staname = staname
|
||||
self.comp = comp
|
||||
self.dt = dt
|
||||
self.tanf = tanf
|
||||
self.vals = vals
|
||||
self.nsamp = vals.size
|
||||
self.tlen = self.dt*self.nsamp
|
||||
|
||||
|
||||
@classmethod
|
||||
def readFromAscii(cls, fid, comp='None'):
|
||||
"""
|
||||
Class method as second constructor to initialize trace from an Ascii file
|
||||
:param fid: file descriptor of already open Ascii file
|
||||
:param comp: optional component selection
|
||||
:return: an instance of the class
|
||||
"""
|
||||
staname = fid.readline().strip()
|
||||
if not staname:
|
||||
print("No station name found. Is file really an Ascii gather file")
|
||||
return None
|
||||
compf = fid.readline().strip()
|
||||
if comp != 'None' and comp != compf:
|
||||
print("Requested component %s not found.", comp)
|
||||
return None
|
||||
h = fid.readline().strip().split()
|
||||
vals = np.array([float(v) for v in fid.readline().strip().split()])
|
||||
dt = float(h[1]); tanf = float(h[2])
|
||||
return cls(staname, compf, dt, tanf, vals)
|
||||
|
||||
|
||||
def absoluteMaximum(self):
|
||||
"""
|
||||
return absolute maximum of trace
|
||||
"""
|
||||
return abs(max(self.vals))
|
||||
|
||||
|
||||
def writeToOpenAsciiFile(self, fid):
|
||||
"""
|
||||
write a seismic trace to an open Ascii file
|
||||
"""
|
||||
fid.write(self.staname + '\n')
|
||||
fid.write(self.comp + '\n')
|
||||
fid.write("{:12d} {:15.9f} {:20.9f}\n".format(self.nsamp, self.dt, self.tanf))
|
||||
for n in range(self.nsamp):
|
||||
fid.write("{:17.8e}".format(self.vals[n]))
|
||||
fid.write("\n")
|
||||
|
||||
|
||||
def printInfo(self):
|
||||
"""
|
||||
Print out some information about the trace
|
||||
"""
|
||||
print("Identification: ", self.staname, self.comp)
|
||||
print("Sampling interval: ", self.dt)
|
||||
print("Number of samples: ", self.nsamp)
|
||||
print("Start time after midnight: ", self.tanf)
|
||||
print("Length of trace: (N*DT) ", self.tlen)
|
||||
15
assignment_07/setupFigure.py
Normal file
15
assignment_07/setupFigure.py
Normal file
@ -0,0 +1,15 @@
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
def SetupFigure(wx, wy, xlabel, ylabel, title, tls):
|
||||
"""
|
||||
Create a new figure frame
|
||||
"""
|
||||
fig = plt.figure(figsize=(wx,wy))
|
||||
ax = fig.add_subplot(1,1,1)
|
||||
ax.set_title(title, fontsize=tls+2)
|
||||
ax.grid(which='major', axis='both', ls=':')
|
||||
ax.set_ylabel(ylabel, fontsize=tls)
|
||||
ax.set_xlabel(xlabel, fontsize=tls)
|
||||
ax.tick_params(labelsize=tls)
|
||||
return fig, ax
|
||||
Loading…
x
Reference in New Issue
Block a user