There will be a webinar on pyGIMLi hosted by SEG on March 19, 2024 at 4 pm CET. Register for free here.

Source code for pygimli.utils.complex

#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""Pygimli base functions to handle complex arrays"""

from math import pi
from pygimli.core.base import isComplex

import numpy as np
import pygimli as pg


isComplex = pg.isComplex

[docs] def toComplex(amp, phi=None): """Convert real values into complex (z = a + ib) valued array. If no phases phi are given, assuming z = amp[0:N] + i amp[N:2N]. If phi is given in (rad) complex values are generated: z = amp*(cos(phi) + i sin(phi)) Parameters ---------- amp: iterable (float) Amplitudes or real unsqueezed real valued array. phi: iterable (float) Phases in neg radiant Returns ------- z: ndarray(dtype=np.complex) Complex values """ if phi is not None: return np.array(amp) * (np.cos(phi) + 1j *np.sin(phi)) N = len(amp) // 2 return np.array(amp[0:N]) + 1j * np.array(amp[N:])
#return np.array(pg.math.toComplex(amp[0:N], amps[N:]))
[docs] def toPolar(z): """Convert complex (z = a + ib) values array into amplitude and phase in radiant If z is real valued we assume its squeezed. Parameters ---------- z: iterable (floats, complex) If z contains of floats and squeezedComplex is assumed [real, imag] Returns ------- amp, phi: ndarray Amplitude amp and phase angle phi in radiant. """ if isComplex(z): return np.abs(z), np.angle(z) else: return toPolar(toComplex(z))
[docs] def squeezeComplex(z, polar=False, conj=False): """Squeeze complex valued array into [real, imag] or [amp, phase(rad)]""" if isinstance(z, (pg.matrix.CSparseMapMatrix, pg.matrix.CSparseMatrix, pg.matrix.CMatrix)): return toRealMatrix(z, conj=conj) if isComplex(z): vals = np.array(z) if conj: vals = np.conj(vals) if polar is True: vals = pg.cat(*toPolar(z)) else: vals = pg.cat(vals.real, vals.imag) return vals return z
[docs] def toRealMatrix(C, conj=False): """Convert complex valued matrix into a real valued Blockmatrix Parameters ---------- C: CMatrix Complex valued matrix conj: bool [False] Fill the matrix as complex conjugated matrix Returns ------- R : pg.matrix.BlockMatrix() """ R = pg.matrix.BlockMatrix() Cr = pg.math.real(A=C) Ci = pg.math.imag(A=C) rId = R.addMatrix(Cr) iId = R.addMatrix(Ci) # we store the mats in R to keep the GC happy after leaving the scope R.addMatrixEntry(rId, 0, 0, scale=1.0) R.addMatrixEntry(rId, Cr.rows(), Cr.cols(), scale=1.0) if conj == True: pg.warn('Squeeze conjugate complex matrix.') R.addMatrixEntry(iId, 0, Cr.cols(), scale=1.0) R.addMatrixEntry(iId, Cr.rows(), 0, scale=-1.0) else: R.addMatrixEntry(iId, 0, Cr.cols(), scale=-1.0) R.addMatrixEntry(iId, Cr.rows(), 0, scale=1.0) return R
[docs] def KramersKronig(f, re, im, usezero=False): """Return real/imaginary parts retrieved by Kramers-Kronig relations. formulas including singularity removal according to Boukamp (1993) """ from scipy.integrate import simps x = 2. * pi * f # omega im2 = np.zeros(im.shape) re2 = np.zeros(im.shape) drdx = np.diff(re) / np.diff(x) # d Re/d omega didx = np.diff(im) / np.diff(x) # d Im/d omega dRedx = np.hstack((drdx[0], (drdx[:-1] + drdx[1:]) / 2, drdx[-1])) dImdx = np.hstack((didx[0], (didx[:-1] + didx[1:]) / 2, didx[-1])) for num, w in enumerate(x): x2w2 = x**2 - w**2 x2w2[num] = 1e-12 fun1 = (re - re[num]) / x2w2 fun1[num] = dRedx[num] / 2 / w im2[num] = -2./pi * w * simps(fun1, x) if usezero: fun2 = (im * w / x - im[num]) / x2w2 re2[num] = 2./pi * w * simps(fun2, x) + re[0] else: fun2 = (im * x - im[num] * w) / x2w2 fun2[num] = (im[num] / w + dImdx[num]) / 2 re2[num] = 2./pi * simps(fun2, x) + re[-1] return re2, im2