SeismologicalDataAnalysis2025/assignment_07/p7_seismometer_correlation.ipynb
2025-06-01 16:41:04 +02:00

309 lines
9.4 KiB
Plaintext

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