# 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.fftpack import ifft, ifftshift, fftshift
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 = ifftshift(ifft(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