dsp/lib/util.py

113 lines
2.1 KiB
Python
Raw Normal View History

2015-10-18 23:06:39 -07:00
import sys
import numpy as np
import scipy.signal as sig
isqrt2 = 1/np.sqrt(2)
tau = 2*np.pi
2015-10-27 04:04:29 -07:00
2017-09-21 04:04:22 -07:00
def dummy(*args, **kwargs):
return None
def lament(*args, **kwargs):
return print(*args, file=sys.stderr, **kwargs)
def toLK(x):
return -0.691 + 10*np.log10(x)
def toQ(bw):
return isqrt2/bw
def toA(db):
return 10**(db/40)
def unwarp(w):
return np.tan(w/2)
def warp(w):
return np.arctan(w)*2
def ceil2(x):
2018-02-16 04:04:24 -08:00
return np.power(2, np.ceil(np.log2(x)).astype(int))
2017-09-21 04:04:22 -07:00
def pad2(x):
return np.r_[x, np.zeros(ceil2(len(x)) - len(x))]
def rfft(src, size):
return np.fft.rfft(src, size*2)
def magnitude(src, size):
return 10*np.log10(np.abs(rfft(src, size))**2)[0:size]
2015-10-30 04:04:36 -07:00
# x axis for plotting above magnitude
2017-09-21 04:04:22 -07:00
def magnitude_x(srate, size):
return np.arange(0, srate/2, srate/2/size)
def degrees_clamped(x):
return ((x*180/np.pi + 180) % 360) - 180
2015-10-18 23:06:39 -07:00
2015-10-19 05:39:37 -07:00
2015-10-18 23:06:39 -07:00
def xsp(precision=4096):
2017-09-21 04:04:22 -07:00
"""
create #precision log-spaced points from
20 Hz (inclusive) to 20480 Hz (exclusive)
"""
xs = np.arange(0, precision)/precision
2015-10-18 23:06:39 -07:00
return 20*1024**xs
2017-09-21 04:04:22 -07:00
2015-10-18 23:06:39 -07:00
def blocks(a, step, size=None):
"""break an iterable into chunks"""
if size is None:
size = step
for start in range(0, len(a), step):
end = start + size
if end > len(a):
break
yield a[start:end]
2017-09-21 04:04:22 -07:00
2015-10-30 04:04:36 -07:00
def convolve_each(s, fir, mode='same', axis=0):
2017-09-21 04:04:22 -07:00
return np.apply_along_axis(
lambda s: sig.fftconvolve(s, fir, mode), axis, s)
2015-10-18 23:06:39 -07:00
def count_channels(s):
2017-03-22 04:04:21 -07:00
if s.ndim < 2:
2015-10-18 23:06:39 -07:00
return 1
return s.shape[1]
2017-09-21 04:04:22 -07:00
2015-10-18 23:06:39 -07:00
def monoize(s):
"""mixes an n-channel signal down to one channel.
technically, it averages a 2D array to be 1D.
existing mono signals are passed through unmodified."""
channels = count_channels(s)
if channels != 1:
2017-03-22 04:04:21 -07:00
s = np.average(s, axis=1)
2015-10-18 23:06:39 -07:00
return s
2016-11-03 04:04:22 -07:00
2017-09-21 04:04:22 -07:00
2016-11-03 04:04:22 -07:00
def div0(a, b):
"""division, whereby division by zero equals zero"""
# http://stackoverflow.com/a/35696047
a = np.asanyarray(a)
b = np.asanyarray(b)
with np.errstate(divide='ignore', invalid='ignore'):
c = np.true_divide(a, b)
2017-09-21 04:04:22 -07:00
c[~np.isfinite(c)] = 0 # -inf inf NaN
2016-11-03 04:04:22 -07:00
return c