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_degree
complex_2_magnitude

Phase Unwrapping
--------------------------------
.. autosummary::
:toctree: generated/

sqrt_phase_unwrap

Unit Conversion
--------------------------------
.. autosummary::
:toctree: generated/

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))

'''
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)

'''
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
'''

[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

return (deg)*pi/180.

'''
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

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

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))*\

# 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)