Source code for skrf.mathFunctions

'''
.. currentmodule:: skrf.mathFunctions
=============================================
mathFunctions (:mod:`skrf.mathFunctions`)
=============================================


Provides commonly used mathematical functions.

Complex Component Conversion
---------------------------------
.. autosummary::
        :toctree: generated/

        complex_2_reim
        complex_2_magnitude
        complex_2_db
        complex_2_radian
        complex_2_degree
        complex_2_magnitude

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

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

'''
import numpy as npy
from numpy import pi,angle,unwrap, real, imag, array
from scipy import signal
from scipy.interpolate import interp1d

global LOG_OF_NEG
LOG_OF_NEG = -100

global INF
INF = 1e99

global ALMOST_ZERO
ALMOST_ZERO = 1e-12

## simple conversions
[docs]def complex_2_magnitude(input): ''' returns the magnitude of a complex number. ''' return abs(input)
[docs]def complex_2_db(input): ''' returns the magnitude in dB of a complex number. returns: 20*log10(|z|) where z is a complex number ''' return magnitude_2_db(npy.abs( input))
def complex_2_db10(input): ''' returns the magnitude in dB of a complex number. returns: 10*log10(|z|) where z is a complex number ''' return mag_2_db10(npy.abs( input))
[docs]def complex_2_radian(input): ''' returns the angle complex number in radians. ''' return npy.angle(input)
[docs]def complex_2_degree(input): ''' returns the angle complex number in radians. ''' return npy.angle(input, deg=True)
def complex_2_quadrature(z): ''' takes a complex number and returns quadrature, which is (length, arc-length from real axis) ''' return ( npy.abs(z), npy.angle(z)*npy.abs(z)) def complex_components(z): ''' break up a complex array into all possible scalar components takes: complex ndarray return: c_real: real part c_imag: imaginary part c_angle: angle in degrees c_mag: magnitude c_arc: arclength from real axis, angle*magnitude ''' return (npy.real(z), npy.imag(z), npy.angle(z,deg=True), complex_2_quadrature(z)[0], complex_2_quadrature(z)[1])
[docs]def complex_2_reim(z): ''' takes: input: complex number or array return: real: real part of input imag: imaginary part of input note: this just calls 'complex_components' ''' out = complex_components(z) return (out[0],out[1])
def magnitude_2_db(input,zero_nan=True): ''' converts linear magnitude to db db is given by 20*log10(|z|) where z is a complex number ''' if zero_nan: out = 20 * npy.log10(input) try: out[npy.isnan(out)] = LOG_OF_NEG except (TypeError): # input is a number not array-like if npy.isnan(out): return LOG_OF_NEG else: out = 20*npy.log10(input) return out mag_2_db = magnitude_2_db def mag_2_db10(input,zero_nan=True): ''' converts linear magnitude to db db is given by 10*log10(|z|) where z is a complex number ''' out = 10 * npy.log10(input) if zero_nan: try: out[npy.isnan(out)] = LOG_OF_NEG except (TypeError): # input is a number not array-like if npy.isnan(out): return LOG_OF_NEG return out def db_2_magnitude(input): ''' converts db to linear magnitude. returns: 10**((z)/20.) where z is a complex number ''' return 10**((input)/20.) db_2_mag = db_2_magnitude def db10_2_mag(input): ''' converts db to linear magnitude returns: 10**((z)/10.) where z is a complex number ''' return 10**((input)/10.) def magdeg_2_reim(mag,deg): ''' converts linear magnitude and phase (in deg) arrays into a complex array ''' return mag*npy.exp(1j*deg*pi/180.) def dbdeg_2_reim(db,deg): ''' converts db magnitude and phase (in deg) arrays into a complex array ''' return magdeg_2_reim(db_2_magnitude(db),deg)
[docs]def db_2_np(x): ''' converts a value in Nepers to dB ''' return (npy.log(10)/20) * x
[docs]def np_2_db(x): ''' converts a value in dB to Nepers ''' return 20/npy.log(10) * x
[docs]def radian_2_degree(rad): return (rad)*180/pi
[docs]def degree_2_radian(deg): return (deg)*pi/180.
[docs]def unwrap_rad(input): ''' unwraps a phase given in radians ''' return unwrap(input,axis=0)
def sqrt_known_sign(z_squared, z_approx): ''' Returns sqrt of complex number, with sign chosen to match `z_approx` Parameters ------------- z_squared : number, array-like the complex to be square-rooted z_approx : number, 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 = npy.sqrt(z_squared) return npy.where( npy.sign(npy.angle(z)) == npy.sign(npy.angle(z_approx)), z, z.conj()) def find_correct_sign(z1,z2,z_approx): ''' 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 : npy.array array built from z1 and z2 by z1 where sign(z1) == sign(z_approx), z2 else ''' return npy.where( npy.sign(npy.angle(z1)) == npy.sign(npy.angle(z_approx)),z1, z2) def find_closest(z1,z2,z_approx): ''' Returns 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 : npy.array array built from z1 and z2 ''' z1_dist = abs(z1-z_approx) z2_dist = abs(z2-z_approx) return npy.where(z1_dist<z2_dist,z1, z2)
[docs]def sqrt_phase_unwrap(input): ''' takes the square root of a complex number with unwrapped phase this idea came from Lihan Chen ''' return npy.sqrt(abs(input))*\ npy.exp(0.5*1j*unwrap_rad(complex_2_radian(input)))
# mathematical functions
[docs]def dirac_delta(x): ''' the Dirac function. can take numpy arrays or numbers returns 1 or 0 ''' return (x==0)*1.+(x!=0)*0.
[docs]def neuman(x): ''' neumans number 2-dirac_delta(x) ''' return 2. - dirac_delta(x)
[docs]def null(A, eps=1e-15): ''' calculates the null space of matrix A. i found this on stack overflow. ''' u, s, vh = npy.linalg.svd(A) null_space = npy.compress(s <= eps, vh, axis=0) return null_space.T
def inf_to_num(x): ''' converts inf and -inf's to large numbers Parameters ------------ x : array-like or number the input array or number Returns ------- ''' #TODO: make this valid for complex arrays try: x[npy.isposinf(x)] = INF x[npy.isneginf(x)] = -1*INF except(TypeError): x = npy.array(x) x[npy.isposinf(x)] = INF x[npy.isneginf(x)] = -1*INF def cross_ratio(a,b,c,d): ''' The cross ratio defined as .. math:: \frac{(a-b)(c-d)}{(a-d)*(c-b)} Parameters ------------- a,b,c,d : complex numbers, or arrays mm ''' return ((a-b)*(c-d))/((a-d)*(c-b)) def complexify(f, name=None): ''' make f(scalar) into f(complex) if the real/imag arguments are not first, then you may specify the name given to them as kwargs. ''' 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
[docs]def complex2Scalar(input): ''' Serializes a list/arary of complex numbers produces the following output for input list `x` x[0].real, x[0].imag, x[1].real, x[1].imag, etc ''' input= npy.array(input) output = [] for k in input: output.append(npy.real(k)) output.append(npy.imag(k)) return npy.array(output).flatten()
[docs]def scalar2Complex(input): ''' inverse of `complex2Scalar` ''' input= npy.array(input) output = [] for k in range(0,len(input),2): output.append(input[k] + 1j*input[k+1]) return npy.array(output).flatten()
def complex2dB(complx): dB = 20 * npy.log10(npy.abs( (npy.real(complx) + 1j*npy.imag(complx) ))) return dB def flatten_c_mat(s, order ='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 : ['F','C'] order of flattening ''' return complex2Scalar(s.flatten(order='F')) def complex2ReIm(complx): return npy.real(complx), npy.imag(complx) def complex2MagPhase(complx,deg=False): return npy.abs(complx), npy.angle(complx,deg=deg) def rand_c(*args): ''' 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 Examples --------- >>> x = rf.rand_c(2,2) ''' s = npy.array(args) return 1-2*npy.random.rand(npy.product(s)).reshape(s) + \ 1j-2j*npy.random.rand(npy.product(s)).reshape(s) def psd2TimeDomain(f,y, windowType='hamming'): '''convert a one sided complex spectrum into a real time-signal. takes f: frequency array, y: complex PSD array windowType: windowing function, defaults to rect returns in the form: [timeVector, signalVector] timeVector is in inverse units of the input variable f, 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 #TODO: make sure windowType exists in scipy.signal if (windowType != 'rect' ): exec("window = signal.%s(%i)" % (windowType,len(f))) y = y * window #create other half of spectrum spectrum = (npy.hstack([npy.real(y[:0:-1]),npy.real(y)])) + \ 1j*(npy.hstack([-npy.imag(y[:0:-1]),npy.imag(y)])) # do the transform df = abs(f[1]-f[0]) T = 1./df timeVector = npy.linspace(-T/2.,T/2,2*len(f)-1) signalVector = npy.fft.ifftshift(npy.fft.ifft(npy.fft.ifftshift(spectrum))) #the imaginary part of this signal should be from fft errors only, signalVector= npy.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, y, d=4, epsilon=1e-9, axis=0): """ 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 [0]_. References ------------ .. [0] 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 """ n = len(x) w = npy.zeros(n) for k in range(n): for i in range(max(0,k-d), min(k+1, n-d)): p = 1.0 for j in range(i,min(n,i+d+1)): if j == k: continue p *= 1/(x[k] - x[j]) w[k] += ((-1)**i)*p if axis != 0: raise NotImplementedError("Axis other than 0 is not implemented") def fx(xi): def find_nearest(a, values, epsilon): idx = npy.abs(npy.subtract.outer(a, values)).argmin(0) return npy.abs(a[idx] - values) < epsilon def find_nearest_value(a, values, y): idx = npy.abs(npy.subtract.outer(a, values)).argmin(0) return y[idx] nearest = find_nearest(x, xi, epsilon) nearest_value = find_nearest_value(x, xi, y) #There needs to be a cleaner way w_shape = [1]*len(y.shape) w_shape[0] = -1 wr = w.reshape(*w_shape) with npy.errstate(divide='ignore', invalid='ignore'): #nans will be fixed later v = npy.sum([y[i]*wr[i]/((xi - x[i]).reshape(*w_shape)) for i in range(n)], axis=0)\ /npy.sum([w[i]/((xi - x[i]).reshape(*w_shape)) for i in range(n)], axis=0) for e, i in enumerate(nearest): if i: v[e] = nearest_value[e] return v return fx def ifft(x): """ Transforms S-parameters to time-domain bandpass. """ return npy.fft.fftshift(npy.fft.ifft(x, axis=0), axes=0) def irfft(x, n=None): """ Transforms S-parameters to time-domain, assuming complex conjugates for values corresponding to negative frequencies. """ return npy.fft.fftshift(npy.fft.irfft(x, axis=0, n=n), axes=0)