""" mathFunctions (:mod:`skrf.mathFunctions`) ============================================= Provides commonly used mathematical functions. Mathematical Constants ---------------------- Some convenient constants are defined in the :mod:`skrf.constants` module. Complex Component Conversion ---------------------------- .. autosummary:: :toctree: generated/ complex_2_reim complex_2_magnitude complex_2_db complex_2_db10 complex_2_radian complex_2_degree complex_2_magnitude complex_2_quadrature complex_components Magnitude and decibel --------------------- .. autosummary:: :toctree: generated/ magnitude_2_db mag_2_db mag_2_db10 db_2_magnitude db_2_mag db10_2_mag Phase Unwrapping ---------------- .. autosummary:: :toctree: generated/ unwrap_rad sqrt_phase_unwrap Unit Conversion --------------- .. autosummary:: :toctree: generated/ radian_2_degree degree_2_radian np_2_db db_2_np feet_2_meter meter_2_feet db_per_100feet_2_db_per_100meter Scalar-Complex Conversion ------------------------- These conversions are useful for wrapping other functions that don't support complex numbers. .. autosummary:: :toctree: generated/ complex2Scalar scalar2Complex Special Functions ----------------- .. autosummary:: :toctree: generated/ dirac_delta neuman null cross_ratio Random Number Generation ------------------------ .. autosummary:: :toctree: generated/ set_rand_rng rand_rng rand_c Various Utility Functions -------------------------- .. autosummary:: :toctree: generated/ psd2TimeDomain rational_interp ifft irfft is_square is_symmetric is_Hermitian is_positive_definite is_positive_semidefinite get_Hermitian_transpose sqrt_known_sign find_correct_sign find_closest inf_to_num rsolve nudge_eig """ from __future__ import annotations from collections.abc import Callable import numpy as np from numpy import imag, pi, real, unwrap from scipy import signal from .constants import ALMOST_ZERO, EIG_COND, EIG_MIN, INF, LOG_OF_NEG, NumberLike # simple conversions def complex_2_magnitude(z: NumberLike): """ Return the magnitude of the complex argument. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- mag : ndarray or scalar """ return np.abs(z) def complex_2_db(z: NumberLike): r""" Return the magnitude in dB of a complex number (as :math:`20\log_{10}(|z|)`).. The magnitude in dB is defined as :math:`20\log_{10}(|z|)` where :math:`z` is a complex number. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- mag20dB : ndarray or scalar """ return magnitude_2_db(np.abs(z)) def complex_2_db10(z: NumberLike): r""" Return the magnitude in dB of a complex number (as :math:`10\log_{10}(|z|)`). The magnitude in dB is defined as :math:`10\log_{10}(|z|)` where :math:`z` is a complex number. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- mag10dB : ndarray or scalar """ return mag_2_db10(np.abs(z)) def complex_2_radian(z: NumberLike): """ Return the angle complex argument in radian. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- ang_rad : ndarray or scalar The counterclockwise angle from the positive real axis on the complex plane in the range ``(-pi, pi]``, with dtype as numpy.float64. """ return np.angle(z) def complex_2_degree(z: NumberLike): """ Returns the angle complex argument in degree. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- ang_deg : ndarray or scalar """ return np.angle(z, deg=True) def complex_2_quadrature(z: NumberLike): r""" Take a complex number and returns quadrature, which is (length, arc-length from real axis) Arc-length is calculated as :math:`|z| \arg(z)`. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- mag : array like or scalar magnitude (length) arc_length : array like or scalar arc-length from real axis: angle*magnitude """ return (np.abs(z), np.angle(z)*np.abs(z)) def complex_2_reim(z: NumberLike): """ Return real and imaginary parts of a complex number. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- real : array like or scalar real part of input imag : array like or scalar imaginary part of input """ return (np.real(z), np.imag(z)) def complex_components(z: NumberLike): """ Break up a complex array into all possible scalar components. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- c_real : array like or scalar real part c_imag : array like or scalar imaginary part c_angle : array like or scalar angle in degrees c_mag : array like or scalar magnitude c_arc : array like or scalar arclength from real axis, angle*magnitude """ return (*complex_2_reim(z), np.angle(z,deg=True), *complex_2_quadrature(z)) def magnitude_2_db(z: NumberLike, zero_nan: bool = True): """ Convert linear magnitude to dB. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers zero_nan : bool, optional Replace NaN with zero. The default is True. Returns ------- z : number or array_like Magnitude in dB given by 20*log10(|z|) """ out = 20 * np.log10(z) if zero_nan: return np.nan_to_num(out, nan=LOG_OF_NEG, neginf=-np.inf) return out mag_2_db = magnitude_2_db def mag_2_db10(z: NumberLike, zero_nan:bool = True): """ Convert linear magnitude to dB. Parameters ---------- z : array_like A complex number or sequence of complex numbers zero_nan : bool, optional Replace NaN with zero. The default is True. Returns ------- z : array_like Magnitude in dB given by 10*log10(|z|) """ out = 10 * np.log10(z) if zero_nan: return np.nan_to_num(out, nan=LOG_OF_NEG, neginf=-np.inf) return out def db_2_magnitude(z: NumberLike): """ Convert dB to linear magnitude. Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- z : number or array_like 10**((z)/20) where z is a complex number """ return 10**((z)/20.) db_2_mag = db_2_magnitude def db10_2_mag(z: NumberLike): """ Convert dB to linear magnitude. Parameters ---------- z : array_like A complex number or sequence of complex numbers Returns ------- z : array_like 10**((z)/10) where z is a complex number """ return 10**((z)/10.) def magdeg_2_reim(mag: NumberLike, deg: NumberLike): """ Convert linear magnitude and phase (in deg) arrays into a complex array. Parameters ---------- mag : number or array_like A complex number or sequence of real numbers deg : number or array_like A complex number or sequence of real numbers Returns ------- z : array_like A complex number or sequence of complex numbers """ return mag*np.exp(1j*deg*pi/180.) def dbdeg_2_reim(db: NumberLike, deg: NumberLike): """ Converts dB magnitude and phase (in deg) arrays into a complex array. Parameters ---------- db : number or array_like A realnumber or sequence of real numbers deg : number or array_like A real number or sequence of real numbers Returns ------- z : array_like A complex number or sequence of complex numbers """ return magdeg_2_reim(db_2_magnitude(db), deg) def db_2_np(db: NumberLike): """ Converts a value in decibel (dB) to neper (Np). Parameters ---------- db : number or array_like A real number or sequence of real numbers Returns ------- np : number or array_like A real number of sequence of real numbers """ return (np.log(10)/20) * db def np_2_db(x: NumberLike): """ Converts a value in Nepers (Np) to decibel (dB). Parameters ---------- np : number or array_like A real number or sequence of real numbers Returns ------- db : number or array_like A real number of sequence of real numbers """ return 20/np.log(10) * x def radian_2_degree(rad: NumberLike): """ Convert angles from radians to degrees. Parameters ---------- rad : number or array_like Angle in radian Returns ------- deg : number or array_like Angle in degree """ return (rad)*180/pi def degree_2_radian(deg: NumberLike): """ Convert angles from degrees to radians. Parameters ---------- deg : number or array_like Angle in radian Returns ------- rad : number or array_like Angle in degree """ return (deg)*pi/180. def feet_2_meter(feet: NumberLike = 1): """ Convert length in feet to meter. 1 foot is equal to 0.3048 meters. Parameters ---------- feet : number or array-like, optional length in feet. Default is 1. Returns ------- meter: number or array-like length in meter See Also -------- meter_2_feet """ return 0.3048*feet def meter_2_feet(meter: NumberLike = 1): """ Convert length in meter to feet. 1 meter is equal to 0.3.28084 feet. Parameters ---------- meter : number or array-like, optional length in meter. Default is 1. Returns ------- feet : number or array-like length in feet See Also -------- feet_2_meter """ return 3.28084*meter def db_per_100feet_2_db_per_100meter(db_per_100feet: NumberLike = 1): """ Convert attenuation values given in dB/100ft to dB/100m. db_per_100meter = db_per_100feet * rf.meter_2_feet() Parameters ---------- db_per_100feet : number or array-like, optional Attenuation in dB/ 100 ft. Default is 1. Returns ------- db_per_100meter : number or array-like Attenuation in dB/ 100 m See Also -------- meter_2_feet feet_2_meter np_2_db db_2_np """ return db_per_100feet * 100 / feet_2_meter(100) def unwrap_rad(phi: NumberLike): """ Unwraps a phase given in radians. Parameters ---------- phi : number or array_like phase in radians Returns ------- phi : number of array_like unwrapped phase in radians """ return unwrap(phi, axis=0) def sqrt_known_sign(z_squared: NumberLike, z_approx: NumberLike): """ Return the square root of a complex number, with sign chosen to match `z_approx`. Parameters ---------- z_squared : number or array-like the complex to be square-rooted z_approx : number or array-like the approximate value of z. sign of z is chosen to match that of z_approx Returns ------- z : number, array-like (same type as z_squared) square root of z_squared. """ z = np.sqrt(z_squared) return np.where( np.sign(np.angle(z)) == np.sign(np.angle(z_approx)), z, z.conj()) def find_correct_sign(z1: NumberLike, z2: NumberLike, z_approx: NumberLike): r""" Create new vector from z1, z2 choosing elements with sign matching z_approx. This is used when you have to make a root choice on a complex number. and you know the approximate value of the root. .. math:: z1,z2 = \pm \sqrt(z^2) Parameters ---------- z1 : array-like root 1 z2 : array-like root 2 z_approx : array-like approximate answer of z Returns ------- z3 : np.array array built from z1 and z2 by z1 where sign(z1) == sign(z_approx), z2 else """ return np.where( np.sign(np.angle(z1)) == np.sign(np.angle(z_approx)),z1, z2) def find_closest(z1: NumberLike, z2: NumberLike, z_approx: NumberLike): """ Return z1 or z2 depending on which is closer to z_approx. Parameters ---------- z1 : array-like root 1 z2 : array-like root 2 z_approx : array-like approximate answer of z Returns ------- z3 : np.array array built from z1 and z2 """ z1_dist = abs(z1-z_approx) z2_dist = abs(z2-z_approx) return np.where(z1_dist>> def f(x): return x >>> f_c = rf.complexify(f) >>> z = 0.2 -1j*0.3 >>> f_c(z) """ def f_c(z, *args, **kw): if name is not None: kw_re = {name: real(z)} kw_im = {name: imag(z)} kw_re.update(kw) kw_im.update(kw) return f(*args, **kw_re) + 1j*f(*args, **kw_im) else: return f(real(z), *args,**kw) + 1j*f(imag(z), *args, **kw) return f_c # old functions just for reference def complex2Scalar(z: NumberLike): """ Serialize a list/array of complex numbers Parameters ---------- z : number or array_like A complex number or sequence of complex numbers Returns ------- re_im : array_like produce the following array for an input z: z[0].real, z[0].imag, z[1].real, z[1].imag, etc. See Also -------- scalar2Complex """ z = np.array(z) re_im = [] for k in z: re_im.append(np.real(k)) re_im.append(np.imag(k)) return np.array(re_im).flatten() def scalar2Complex(s: NumberLike): """ Unserialize a list/array of real and imag numbers into a complex array. Inverse of :func:`complex2Scalar`. Parameters ---------- s : array_like an array with real and imaginary parts ordered as: z[0].real, z[0].imag, z[1].real, z[1].imag, etc. Returns ------- z : Number or array_like A complex number or sequence of complex number See Also -------- complex2Scalar """ s = np.array(s) z = [] for k in range(0,len(s),2): z.append(s[k] + 1j*s[k+1]) return np.array(z).flatten() def flatten_c_mat(s: NumberLike, order: str = 'F'): """ Take a 2D (mxn) complex matrix and serialize and flatten it. by default (using order='F') this generates the following from a 2x2 [s11,s12;s21,s22]->[s11re,s11im,s21re,s12im, ...] Parameters ------------ s : ndarray input 2D array order : string, optional either 'F' or 'C', for the order of flattening """ return complex2Scalar(s.flatten(order='F')) _global_rng = np.random.default_rng() def set_rand_rng(rng: np.random.Generator) -> None: """ Set the global :mod:`numpy` random number generator instance for generating random numbers in scikit-rf. This is useful for fixing a random seed for reproducible Monte Carlo analysis. This function is expected to be called before using any scikit-rf features, since it's a global variable and thread-unsafe. To temporarily change the random number generator of a particular method (e.g. :meth:`skrf.media.Media.random`), use the ``rng`` argument instead. Parameters ----------- rng : :class:`numpy.random.Generator` Any random number generator accepted by :mod:`numpy`. Examples --------- >>> set_rand_rng(np.random.default_rng(seed=42)) """ global _global_rng _global_rng = rng def rand_rng() -> np.random.Generator: """ Obtain the global :mod:`numpy` random number generator instance. By default, it is :func:`numpy.random.default_rng`. """ return _global_rng def rand_c(*size, rng=None) -> np.ndarray: """ Creates a complex random array of shape s. The bounds on real and imaginary values are (-1,1) Parameters ----------- s : list-like shape of array rng : :class:`numpy.random.Generator` or None Any random number generator accepted by :mod:`numpy`. Examples --------- >>> x1 = rf.rand_c(2, 2) >>> x2 = rf.rand_c(2, 2, np.random.default_rng(seed=42)) """ if rng is None: rng = _global_rng return 1-2*rng.random(size) + \ 1j-2j*rng.random(size) def psd2TimeDomain(f: np.ndarray, y: np.ndarray, windowType: str = 'hamming'): """ Convert a one sided complex spectrum into a real time-signal. Parameters ---------- f : list or np.ndarray frequency array y : list of np.ndarray complex PSD array windowType: string windowing function, defaults to 'hamming'' Returns ------- timeVector : array_like inverse units of the input variable f, signalVector : array_like Note ---- If spectrum is not baseband then, `timeSignal` is modulated by `exp(t*2*pi*f[0])`. So keep in mind units. Also, due to this, `f` must be increasing left to right. """ # apply window function # make sure windowType exists in scipy.signal if callable(getattr(signal, windowType)) and (windowType != 'rect' ): window = getattr(signal, windowType)(len(f)) y = y * window #create other half of spectrum spectrum = (np.hstack([np.real(y[:0:-1]),np.real(y)])) + \ 1j*(np.hstack([-np.imag(y[:0:-1]),np.imag(y)])) # do the transform df = abs(f[1]-f[0]) T = 1./df timeVector = np.linspace(-T/2.,T/2,2*len(f)-1) signalVector = np.fft.ifftshift(np.fft.ifft(np.fft.ifftshift(spectrum))) #the imaginary part of this signal should be from fft errors only, signalVector= np.real(signalVector) # the response of frequency shifting is # exp(1j*2*pi*timeVector*f[0]) # but I would have to manually undo this for the inverse, which is just # another variable to require. The reason you need this is because # you can't transform to a bandpass signal, only a lowpass. # return timeVector, signalVector def rational_interp( x: np.ndarray, y: np.ndarray, d: int = 4, epsilon: float = 1e-9, axis: int = 0, assume_sorted: bool = False) -> Callable: """ Interpolates function using rational polynomials of degree `d`. Interpolating function is singular when xi is exactly one of the original x points. If xi is closer than epsilon to one of the original points, then the value at that points is returned instead. Implementation is based on [#]_. Parameters ---------- x : np.ndarray y : np.ndarray d : int, optional order of the polynomial, by default 4 epsilon : float, optional numerical tolerance, by default 1e-9 axis : int, optional axis to operate on, by default 0 assume_sorted : bool, optional If False, values of x can be in any order and they are sorted first. If True, x has to be an array of monotonically increasing values. Returns ------- fx : Callable Interpolate function Raises ------ NotImplementedError if axis != 0. References ------------ .. [#] M. S. Floater and K. Hormann, "Barycentric rational interpolation with no poles and high rates of approximation," Numer. Math., vol. 107, no. 2, pp. 315-331, Aug. 2007 """ if axis != 0: raise NotImplementedError("Axis other than 0 is not implemented") if not assume_sorted: sort_indices = np.argsort(x, axis=axis) x = x[sort_indices] y = y[sort_indices] n = len(x) if n <= d: raise ValueError('Not enough x-axis points') w = np.zeros(n) # Scaling to give close to 1 weights hd = (x[n//2] - x[n//2-1])**d for k in range(n): for i in range(max(0,k-d), min(k+1, n-d)): p, xk = hd, x[k] for j in range(i,min(n,i+d+1)): if j == k: continue p /= (xk - x[j]) if i % 2 == 1: w[k] -= p else: w[k] += p # Add dimensions to match y shape w_shape = [1]*len(y.shape) w_shape[0] = -1 w = w.reshape(w_shape) def fx(xi): # The method will divide by zero if new x value is exactly existing x value. # To avoid this we need to check for too close values and replace them with # y value at that position. idx = np.searchsorted(x, xi) idx[idx == len(x)] = len(x) - 1 nearest_idx = np.where(np.abs(x[idx] - xi) < epsilon)[0] nearest_value = y[idx[nearest_idx]] xi = xi.reshape(*w_shape) with np.errstate(divide='ignore', invalid='ignore'): assert axis == 0 tmp = [w[i] / (xi - x[i]) for i in range(n)] v = sum(y[i] * tmp[i] for i in range(n)) / sum(tmp) # Fix divide by zero errors for e, i in enumerate(nearest_idx): v[i] = nearest_value[e] return v return fx def ifft(x: np.ndarray) -> np.ndarray: """ Transforms S-parameters to time-domain bandpass. Parameters ---------- x : array_like S-parameters vs frequency array. Returns ------- X : array_like Fourier transformed array See Also -------- irfft """ return np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(x, axes=0), axis=0), axes=0) def irfft(x: np.ndarray, n:int = None) -> np.ndarray: """ Transforms S-parameters to time-domain, assuming complex conjugates for values corresponding to negative frequencies. Parameters ---------- x : array_like S-parameters vs frequency array. n : int, optional n parameter passed to :func:`numpy.fft.irfft`. Defaults to None. Returns ------- X : array_like Fourier transformed array See Also -------- ifft """ return np.fft.fftshift(np.fft.irfft(x, axis=0, n=n), axes=0) def is_square(mat: np.ndarray) -> bool: """ Tests whether mat is a square matrix. Parameters ---------- mat : np.ndarray Matrix to test for being square Returns ------- res : boolean See Also -------- is_unitary is_symmetric """ return mat.shape[0] == mat.shape[1] def is_unitary(mat: np.ndarray, tol: float = ALMOST_ZERO) -> bool: """ Tests mat for unitariness. Parameters ---------- mat : np.ndarray Matrix to test for unitariness tol : float Absolute tolerance. Defaults to :data:`ALMOST_ZERO` Returns ------- res : boolean or array of boolean See Also -------- is_square is_symmetric """ if not is_square(mat): return False return np.allclose(get_Hermitian_transpose(mat) @ mat, np.identity(mat.shape[0]), atol=tol) def is_symmetric(mat: np.ndarray, tol: int = ALMOST_ZERO) -> bool: """ Tests mat for symmetry. Parameters ---------- mat : np.ndarray Matrix to test for symmetry tol : float, optional Absolute tolerance. Defaults to :data:`ALMOST_ZERO` Returns ------- res : boolean or array of boolean See Also -------- is_square is_unitary """ if not is_square(mat): return False return np.allclose(mat, mat.transpose(), atol=tol) def get_Hermitian_transpose(mat: np.ndarray) -> np.ndarray: """ Returns the conjugate transpose of mat. Parameters ---------- mat : np.ndarray Matrix to compute the conjugate transpose of Returns ------- mat : np.ndarray """ return mat.transpose().conjugate() def is_Hermitian(mat: np.ndarray, tol: float = ALMOST_ZERO) -> bool: """ Tests whether mat is Hermitian. Parameters ---------- mat : np.ndarray Matrix to test for being Hermitian tol : float Absolute tolerance Returns ------- res : boolean See Also -------- is_positive_definite is_positive_semidefinite """ if not is_square(mat): return False return np.allclose(mat, get_Hermitian_transpose(mat), atol=tol) def is_positive_definite(mat: np.ndarray, tol: float = ALMOST_ZERO) -> bool: """ Tests mat for positive definiteness. Verifying that (1) mat is symmetric (2) it's possible to compute the Cholesky decomposition of mat. Parameters ---------- mat : np.ndarray Matrix to test for positive definiteness tol : float, optional Absolute tolerance. Defaults to :data:`ALMOST_ZERO` Returns ------- res : bool or array of bool See Also -------- is_Hermitian is_positive_semidefinite """ if not is_Hermitian(mat, tol=tol): return False try: np.linalg.cholesky(mat) return True except np.linalg.LinAlgError: return False def is_positive_semidefinite(mat: np.ndarray, tol: float = ALMOST_ZERO) -> bool: """ Tests mat for positive semidefiniteness. Checking whether all eigenvalues of mat are nonnegative within a certain tolerance Parameters ---------- mat : np.ndarray Matrix to test for positive semidefiniteness tol : float, optional Absolute tolerance in determining nonnegativity due to loss of precision when computing the eigenvalues of mat. Defaults to :data:`ALMOST_ZERO` Returns ------- res : bool or array of bool See Also -------- is_Hermitian is_positive_definite """ if not is_Hermitian(mat): return False try: v = np.linalg.eigvalsh(mat) except np.linalg.LinAlgError: return False return np.all(v > -tol) def rsolve(A: np.ndarray, B: np.ndarray) -> np.ndarray: r"""Solves x @ A = B. Calls numpy.linalg.solve with transposed matrices. Same as B @ np.linalg.inv(A) but avoids calculating the inverse and should be numerically slightly more accurate. Input should have dimension of similar to (nfreqs, nports, nports). Parameters ---------- A : np.ndarray B : np.ndarray Returns ------- x : np.ndarray """ return np.transpose(np.linalg.solve(np.transpose(A, (0, 2, 1)).conj(), np.transpose(B, (0, 2, 1)).conj()), (0, 2, 1)).conj() def nudge_eig(mat: np.ndarray, cond: float | None = None, min_eig: float | None = None) -> np.ndarray: r"""Nudge eigenvalues with absolute value smaller than max(cond * max(eigenvalue), min_eig) to that value. Can be used to avoid singularities in solving matrix equations. Input should have dimension of similar to (nfreqs, nports, nports). Parameters ---------- mat : np.ndarray Matrices to nudge cond : float, optional Minimum eigenvalue ratio compared to the maximum eigenvalue. Default value is set by `skrf.constants.EIG_COND`. min_eig : float, optional Minimum eigenvalue. Default value is set by `skrf.constants.EIG_MIN`. Returns ------- res : np.ndarray Nudged matrices """ # use current constants if not cond: cond = EIG_COND if not min_eig: min_eig = EIG_MIN # Eigenvalues and vectors eigw, eigv = np.linalg.eig(mat) # Max eigenvalue for each frequency max_eig = np.amax(np.abs(eigw), axis=1) # Calculate mask for positions where problematic eigenvalues are mask = np.logical_or(np.abs(eigw) < cond * max_eig[:, None], np.abs(eigw) < min_eig) if not mask.any(): # Nothing to do. Return the original array. return mat mask_cond = cond * np.repeat(max_eig[:, None], mat.shape[-1], axis=-1)[mask] mask_min = min_eig * np.ones(mask_cond.shape) # Correct the eigenvalues eigw[mask] = np.maximum(mask_cond, mask_min) # Now assemble the eigendecomposited matrices back e = np.zeros_like(mat) np.einsum('ijj->ij', e)[...] = eigw return rsolve(eigv, eigv @ e)