notwa/dsp
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

update 33

Browse Source
This commit is contained in:
Connor Olding 2017-09-21 04:04:22 -07:00
parent c678fe4512
commit d41df7f132
18 changed files with 462 additions and 260 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

138
lib/__init__.py
View File

@ -12,141 +12,13 @@ from .bs import *
from .cepstrum import *
from .windowing import *
from .piir import *
import numpy as np
def analog(b, a):
import sympy as sym
w,s = sym.symbols('w s')
filt_expr = sym.Poly(b, s)/sym.Poly(a, s)
mag_expr = abs(filt_expr.subs({s: w*sym.I}))**2
return sym.lambdify(w, mag_expr, 'numpy')
def makemag(w0, ba, gain=0):
f = analog(*ba)
def magf(w):
a = f(w/w0)
a[0] = 1e-35
a = np.log10(a)*10 + gain
a[0] = a[1] # safety measure
return a
return magf
def test_filter_raw(ba, fc=1000, gain=0, precision=4096):
fig, ax = new_response(ymin=-24, ymax=24)
xs = xsp(precision)
ax.semilogx(xs, makemag(fc, ba, gain)(xs))
def test_filter(ff, A=toA(12), Q=toQ(1), **kwargs):
test_filter_raw(ff(A, Q), **kwargs)
def neonpink(xs):
lament("neonpink(): DEPRECATED; use tilter2(xs, 'raw') instead.")
return tilter2(xs, 'raw')
def c_render(cascade, precision=4096):
# TODO: deprecate in favor of tilter2
xs = xsp(precision)
return xs, tilter2(xs, cascade)
def c_render2(xs, cascade, phase=False):
"""c_render optimized and specifically for first/second-order filters"""
if phase:
return c_render3(xs, cascade, mode='phase')
else:
return c_render3(xs, cascade, mode='magnitude')
def c_render3(xs, cascade, mode='magnitude'):
"""c_render optimized and specifically for first/second-order filters"""
import numexpr as ne
j = np.complex(0, 1)
# obviously this could be extended to higher orders
eq2 = '(b0 + b1*s + b2*s**2)/(a0 + a1*s + a2*s**2)'
eq1 = '(b0 + b1*s)/(a0 + a1*s)'
if mode == 'magnitude':
fmt = 'real(log10(abs({})**2)*10 + gain)'
elif mode == 'phase' or mode == 'group delay':
fmt = '-arctan2(imag({0}), real({0}))' # gross
else:
raise Exception("c_render3(): unknown mode: {}".format(mode))
ys = np.zeros(len(xs))
for f in cascade:
freq, ba, gain = f
b, a = ba
if len(b) == 3 and len(a) == 3:
eq = fmt.format(eq2)
b2, b1, b0 = b
a2, a1, a0 = a
elif len(b) == 2 and len(a) == 2:
eq = fmt.format(eq1)
b1, b0 = b
a1, a0 = a
else:
raise Exception("incompatible cascade; consider using c_render instead")
if mode == 'group delay':
# approximate derivative of phase by slope of tangent line
step = 2**-8
fa = freq - step
fb = freq + step
s = xs/fa*j
ya = ne.evaluate(eq)
s = xs/fb*j
yb = ne.evaluate(eq)
slope = (yb - ya)/(2*step)
ys += -slope/(xs/freq*tau)
else:
s = xs/freq*j
ys += ne.evaluate(eq)
if mode == 'phase':
ys = degrees_clamped(ys)
return ys
def firize(xs, ys, n=4096, srate=44100, ax=None):
import scipy.signal as sig
if ax:
ax.semilogx(xs, ys, label='desired')
xf = xs/srate*2
yg = 10**(ys/20)
xf = np.r_[0, xf, 1]
yg = np.r_[0, yg, yg[-1]]
b = sig.firwin2(n, xf, yg, antisymmetric=True)
if ax:
_, ys = sig.freqz(b, worN=xs/srate*tau)
ys = 20*np.log10(np.abs(ys))
ax.semilogx(xs, ys, label='FIR ({} taps)'.format(n))
ax.legend(loc=8)
return b
def tilter(xs, ys, tilt):
"""tilts a magnitude plot by some decibels, or by equalizer curve."""
lament("tilter(): DEPRECATED; use ys -= tilter2(xs, tilt) instead.")
return xs, ys - tilter2(xs, tilt)
def tilter2(xs, tilt): # TODO: rename
noise = np.zeros(xs.shape)
if isinstance(tilt, str) and tilt in cascades:
tilt = cascades[tilt]
if isinstance(tilt, list):
c = [makemag(*f) for f in tilt]
for f in c:
noise += f(xs)
elif isinstance(tilt, int) or isinstance(tilt, float):
noise = tilt*(np.log2(1000) - np.log2(xs + 1e-35))
return noise
from .mag import *
from .plotwav import *
# this is similar to default behaviour of having no __all__ variable at all,
# but ours ignores modules as well. this allows for `import sys` and such
# without clobbering `from our_module import *`.
__all__ = [o for o in locals() if type(o) != 'module' and not o.startswith('_')]
__all__ = [
o for o in locals()
if type(o) != 'module' and not o.startswith('_')]

96
lib/bq.py
View File

@ -4,39 +4,79 @@ import scipy.signal as sig
from .util import *
from .planes import s2z
bq_run = lambda bq, xs: sig.lfilter(*bq, x=xs, axis=0)
nfba = lambda b, a: (1/tau, (b, a), 0)
nf = lambda t, f, g, bw, mg: (f, t(toA(g), toQ(bw)), mg)
# PEP 8 fucking destroyed this file. I'm sorry.
def bq_run(bq, xs):
return sig.lfilter(*bq, x=xs, axis=0)
def nfba(b, a):
return (1/tau, (b, a), 0)
def nf(t, f, g, bw, mg):
return (f, t(toA(g), toQ(bw)), mg)
def LP1(A, Q):
return ((0, 1), (1, 1))
def HP1(A, Q):
return ((1, 0), (1, 1))
def LS1(A, Q):
return ((1, A), (1, 1/A))
def HS1(A, Q):
return ((A, 1), (1/A, 1))
LP1 = lambda A, Q: ((0,1),(1,1))
HP1 = lambda A, Q: ((1,0),(1,1))
LS1 = lambda A, Q: ((1,A),(1,1/A))
HS1 = lambda A, Q: ((A,1),(1/A,1))
# patterns observed, in case some simplification could be done:
# a always gets divided by A instead of multiplied
# b1 and a1 always /= Q
LP2 = lambda A, Q: ((0, 0, 1),
(1, 1/Q, 1))
HP2 = lambda A, Q: ((1, 0, 0),
(1, 1/Q, 1))
PE2 = lambda A, Q: ((1, A/Q, 1),
(1, 1/A/Q, 1))
AP2 = lambda A, Q: ((1, -1/Q, 1),
(1, 1/Q, 1))
BP2a= lambda A, Q: ((0, -A/Q, 0),
(1, 1/A/Q, 1))
BP2b= lambda A, Q: ((0,-A*A/Q, 0),
(1, 1/Q, 1))
NO2 = lambda A, Q: ((1, 0, 1),
(1, 1/Q, 1))
LS2 = lambda A, Q: ((1, np.sqrt(A)/Q, A),
(1, 1/np.sqrt(A)/Q, 1/A))
HS2 = lambda A, Q: ((A, np.sqrt(A)/Q, 1),
(1/A, 1/np.sqrt(A)/Q, 1))
def LP2(A, Q):
return ((0, 0, 1), (1, 1/Q, 1))
gen_filters = lambda cascade, srate: [
s2z(*f[1], fc=f[0], srate=srate, gain=10**(f[2]/20)) for f in cascade
]
def HP2(A, Q):
return ((1, 0, 0), (1, 1/Q, 1))
def PE2(A, Q):
return ((1, A/Q, 1), (1, 1/A/Q, 1))
def AP2(A, Q):
return ((1, -1/Q, 1), (1, 1/Q, 1))
def BP2a(A, Q):
return ((0, -A/Q, 0), (1, 1/A/Q, 1))
def BP2b(A, Q):
return ((0, -A*A/Q, 0), (1, 1/Q, 1))
def NO2(A, Q):
return ((1, 0, 1), (1, 1/Q, 1))
def LS2(A, Q):
return ((1, np.sqrt(A)/Q, A), (1, 1/np.sqrt(A)/Q, 1/A))
def HS2(A, Q):
return ((A, np.sqrt(A)/Q, 1), (1/A, 1/np.sqrt(A)/Q, 1))
def gen_filters(cascade, srate):
return [
s2z(*f[1], fc=f[0], srate=srate, gain=10**(f[2]/20)) for f in cascade
]

15
lib/bs.py
View File

@ -3,6 +3,7 @@ from . import blocks, convolve_each, gen_filters, cascades, bq_run, toLK
import numpy as np
import matplotlib.pyplot as plt
def BS1770_3(s, srate, filters=None, window=0.4, overlap=0.75,
gate=10, absolute_gate=70, detail=False):
if filters is None:
@ -10,7 +11,7 @@ def BS1770_3(s, srate, filters=None, window=0.4, overlap=0.75,
sf = np.copy(s)
for f in filters:
if len(f) is 2: # dumb but effective
if len(f) is 2: # dumb way to tell what type we're given.
sf = bq_run(f, sf)
else:
sf = convolve_each(sf, f, 'same')
@ -35,6 +36,7 @@ def BS1770_3(s, srate, filters=None, window=0.4, overlap=0.75,
else:
return toLK(avg_g10), toLK(avg_g70), LKs, threshold
def BS_plot(ys, g10=None, g70=None, threshold=None, fig=None, ax=None):
if g10:
center = np.round(g10)
@ -47,25 +49,25 @@ def BS_plot(ys, g10=None, g70=None, threshold=None, fig=None, ax=None):
if ax is None:
ax = fig.gca()
if False: # histogram
if False: # histogram
ax.hist(ys, bins=bins, normed=True, facecolor='g', alpha=0.5)
ax.xlim(bins[0], bins[-1])
ax.ylim(0, 1)
ax.grid(True, 'both')
ax.xlabel('loudness (LKFS)')
ax.ylabel('probability')
fig.set_size_inches(10,4)
fig.set_size_inches(10, 4)
xs = np.arange(len(ys))
#ax.plot(xs, ys, color='#066ACF', linestyle=':', marker='d', markersize=2)
# ax.plot(xs, ys, color='#066ACF', linestyle=':', marker='d', markersize=2)
ax.plot(xs, ys, color='#1459E0')
ax.set_xlim(xs[0], xs[-1])
ax.set_ylim(-70, 0)
ax.grid(True, 'both', 'y')
ax.set_xlabel('bin')
ax.set_ylabel('loudness (LKFS)')
fig.set_size_inches(12,5)
#_, _, ymin, _ = ax.axis()
fig.set_size_inches(12, 5)
# _, _, ymin, _ = ax.axis()
if threshold:
ax.axhspan(-70, threshold, facecolor='r', alpha=1/5)
if g10:
@ -75,6 +77,7 @@ def BS_plot(ys, g10=None, g70=None, threshold=None, fig=None, ax=None):
return fig, ax
def normalize(s, srate):
"""performs BS.1770-3 normalization and returns inverted gain."""
db = BS1770_3(s, srate)

19
lib/cepstrum.py
View File

@ -1,12 +1,19 @@
import numpy as np
from .util import pad2
# fast cepstrum and inverted fast cepstrum
fcs = lambda s: np.fft.ifft(np.log(np.fft.fft(s)))
ifcs = lambda s: np.fft.fft(np.exp(np.fft.ifft(s)))
# magnitude
mcs = lambda s: (np.abs(np.fft.ifft(np.log(np.abs(np.fft.fft(s))**2)))**2)[:len(s)//2]
def fcs(s): # fast cepstrum
return np.fft.ifft(np.log(np.fft.fft(s)))
def ifcs(s): # inverted fast cepstrum
return np.fft.fft(np.exp(np.fft.ifft(s)))
def mcs(s): # magnitude
return (np.abs(np.fft.ifft(np.log(np.abs(np.fft.fft(s))**2)))**2
)[:len(s)//2]
def clipdb(s, cutoff=-100):
as_ = np.abs(s)
@ -16,6 +23,7 @@ def clipdb(s, cutoff=-100):
thresh = mas*10**(cutoff/20)
return np.where(as_ < thresh, thresh, s)
def fold(r):
# via https://ccrma.stanford.edu/~jos/fp/Matlab_listing_fold_m.html
# Fold left wing of vector in "FFT buffer format" onto right wing
@ -33,6 +41,7 @@ def fold(r):
rw = np.r_[r[0], rf, np.zeros(n-nt-1)]
return rw
def minphase(s, pad=True, os=False):
# via https://ccrma.stanford.edu/~jos/fp/Matlab_listing_mps_m.html
# TODO: actual oversampling

17
lib/data.py
View File

@ -14,18 +14,20 @@ cascades = {
(1501, HS2(toA(4), toQ(1)), 0),
(38.135457, HP2(0, 0.5003268), np.log10(1.004995)*20),
],
# "neon pink"
'raw': [
nf(LP1, 20, 0, 1, 29),
nf(HS1, 800, 12, 1, 0),
# i don't use the exact _bq2 coeffecients here for legacy reasons
( 45, HP2( 0, 1.32), 0.5), # roughly estimates
( 45, HP2( 0, 0.54), 0.5), # a 4-pole butterworth highpass
( 45, HP2( 0, 1.32), 0.5), # roughly estimates
( 45, HP2( 0, 0.54), 0.5), # a 4-pole butterworth highpass
nf(LP2, 14000, 0, 1.33, 0),
],
# like neon pink but for feeding into RMS
'raw2': [
(10000, HP1(0,0), 26),
(10000, HP1(0, 0), 26),
( 750, HS2(toA(-10), toQ(1.33)), 0),
( 45, HP2(0, 1.32), 0.5),
( 45, HP2(0, 0.54), 0.5),
@ -33,14 +35,16 @@ cascades = {
( 250, PE2(toA(3), toQ(1.33)), -1),
( 4000, PE2(toA(3), toQ(1.33)), -1),
],
# loosely based on the equal loudness contour at 60dB or something
'raw_ELC': [
( 40, HP2(0, toQ(1.33)), 0),
( 400, HP1(0,0), 6),
( 400, HP1(0, 0), 6),
( 1400, PE2(toA(-3), toQ(1.33)), 1),
( 4000, PE2(toA(5), toQ(1.00)),-1.5),
( 4000, LP2(0, toQ(1.33)), 1.5),
],
# here's the ideas written out:
# low (<40) freqs dont contribute much to ears (feel doesnt count.)
# high (>14000) freqs are mostly unheard.
@ -55,6 +59,7 @@ cascades = {
( 8000, PE2(toA(3), toQ(1.00)), 0.0),
(10000, LP2(0, toQ(0.50)),-0.5),
],
'tilt_test': [
(10000, HP1(0,0), 30),
( 1000, HS1(toA(-16), 0), 1.5),
@ -62,6 +67,7 @@ cascades = {
( 40, HP2(0, toQ(1.00)), 0.0),
(10000, LP1(0, 0), 0.0),
],
# average curve of my 227 favorite songs
'np2': [
nf(LP1, 20, 0, 1, 32),
@ -72,6 +78,7 @@ cascades = {
nf(LS2, 38, -9, 1.00, 0),
nf(PE2, 64, 4.5, 1.20, 0),
],
# same but for the side channel
'np2s': [
nf(LP1, 20, 0, 1, 32),
@ -79,7 +86,5 @@ cascades = {
nf(LP2, 14000, 0, 1.33, 0),
nf(HP2, 90, 0, 1.11, 0),
nf(PE2, 30, -9.5, 1.00, 0),
#(17500, LP2(0, _bq2a), 0),
#(17500, LP2(0, _bq2b), 0),
],
}

5
lib/fft.py
View File

@ -3,6 +3,7 @@ from . import rfft
import numpy as np
import scipy.signal as sig
def magnitudes_window_setup(s, size=8192):
L = s.shape[0]
overlap = 0.661
@ -10,6 +11,7 @@ def magnitudes_window_setup(s, size=8192):
segs = np.ceil(L/step)
return step, segs
def magnitudes(s, size=8192):
import scipy.linalg as linalg
@ -34,7 +36,8 @@ def magnitudes(s, size=8192):
yield power[0:size]
count += 1
#assert(segs == count) # this is probably no good in a generator
# assert(segs == count) # this is probably no good in a generator
def averfft(s, size=8192):
"""calculates frequency magnitudes by fft and averages them together."""

145
lib/mag.py Normal file
View File

@ -0,0 +1,145 @@
from . import toA, toQ, cascades
import numpy as np
def analog(b, a):
import sympy as sym
w, s = sym.symbols('w s')
filt_expr = sym.Poly(b, s)/sym.Poly(a, s)
mag_expr = abs(filt_expr.subs({s: w*sym.I}))**2
return sym.lambdify(w, mag_expr, 'numpy')
def makemag(w0, ba, gain=0):
f = analog(*ba)
def magf(w):
a = f(w/w0)
a[0] = 1e-35
a = np.log10(a)*10 + gain
a[0] = a[1] # safety measure
return a
return magf
def test_filter_raw(ba, fc=1000, gain=0, precision=4096):
fig, ax = new_response(ymin=-24, ymax=24)
xs = xsp(precision)
ax.semilogx(xs, makemag(fc, ba, gain)(xs))
def test_filter(ff, A=toA(12), Q=toQ(1), **kwargs):
test_filter_raw(ff(A, Q), **kwargs)
def neonpink(xs):
lament("neonpink(): DEPRECATED; use tilter2(xs, 'raw') instead.")
return tilter2(xs, 'raw')
def c_render(cascade, precision=4096):
# TODO: deprecate in favor of tilter2
xs = xsp(precision)
return xs, tilter2(xs, cascade)
def c_render2(xs, cascade, phase=False):
"""c_render optimized and specifically for first/second-order filters"""
if phase:
return c_render3(xs, cascade, mode='phase')
else:
return c_render3(xs, cascade, mode='magnitude')
def c_render3(xs, cascade, mode='magnitude'):
"""c_render optimized and specifically for first/second-order filters"""
import numexpr as ne
j = np.complex(0, 1)
# obviously this could be extended to higher orders
eq2 = '(b0 + b1*s + b2*s**2)/(a0 + a1*s + a2*s**2)'
eq1 = '(b0 + b1*s)/(a0 + a1*s)'
if mode == 'magnitude':
fmt = 'real(log10(abs({})**2)*10 + gain)'
elif mode == 'phase' or mode == 'group delay':
fmt = '-arctan2(imag({0}), real({0}))' # gross
else:
raise Exception("c_render3(): unknown mode: {}".format(mode))
ys = np.zeros(len(xs))
for f in cascade:
freq, ba, gain = f
b, a = ba
if len(b) == 3 and len(a) == 3:
eq = fmt.format(eq2)
b2, b1, b0 = b
a2, a1, a0 = a
elif len(b) == 2 and len(a) == 2:
eq = fmt.format(eq1)
b1, b0 = b
a1, a0 = a
else:
raise Exception(
"incompatible cascade; consider using c_render instead")
if mode == 'group delay':
# approximate derivative of phase by slope of tangent line
step = 2**-8
fa = freq - step
fb = freq + step
s = xs/fa*j
ya = ne.evaluate(eq)
s = xs/fb*j
yb = ne.evaluate(eq)
slope = (yb - ya)/(2*step)
ys += -slope/(xs/freq*tau)
else:
s = xs/freq*j
ys += ne.evaluate(eq)
if mode == 'phase':
ys = degrees_clamped(ys)
return ys
def firize(xs, ys, n=4096, srate=44100, ax=None):
import scipy.signal as sig
if ax:
ax.semilogx(xs, ys, label='desired')
xf = xs/srate*2
yg = 10**(ys/20)
xf = np.r_[0, xf, 1]
yg = np.r_[0, yg, yg[-1]]
b = sig.firwin2(n, xf, yg, antisymmetric=True)
if ax:
_, ys = sig.freqz(b, worN=xs/srate*tau)
ys = 20*np.log10(np.abs(ys))
ax.semilogx(xs, ys, label='FIR ({} taps)'.format(n))
ax.legend(loc=8)
return b
def tilter(xs, ys, tilt):
"""tilts a magnitude plot by some decibels, or by equalizer curve."""
lament("tilter(): DEPRECATED; use ys -= tilter2(xs, tilt) instead.")
return xs, ys - tilter2(xs, tilt)
def tilter2(xs, tilt): # TODO: rename
noise = np.zeros(xs.shape)
if isinstance(tilt, str) and tilt in cascades:
tilt = cascades[tilt]
if isinstance(tilt, list):
c = [makemag(*f) for f in tilt]
for f in c:
noise += f(xs)
elif isinstance(tilt, int) or isinstance(tilt, float):
noise = tilt*(np.log2(1000) - np.log2(xs + 1e-35))
return noise

12
lib/nsf.py
View File

@ -2,8 +2,10 @@
import numpy as np
def LPB(n):
# via https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source
# via:
# https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source
"""n-th order butterworth low-pass filter cascade
-3 dB at center frequency."""
@ -24,8 +26,10 @@ def LPB(n):
series += [(num, den)]
return series
def LPC(n, ripple, type=1):
# via https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source
# via:
# https://github.com/vinniefalco/DSPFilters/blob/master/shared/DSPFilters/source
# FIXME: type 2 has wrong center frequency?
"""n-th order chebyshev low-pass filter cascade
@ -46,9 +50,9 @@ def LPC(n, ripple, type=1):
v0 = np.arcsinh(1/eps)/n
else:
if type == 2:
v0 = 0 # allpass?
v0 = 0 # allpass?
else:
v0 = 1 # butterworth
v0 = 1 # butterworth
sinh_v0 = -np.sinh(v0)
cosh_v0 = np.cosh(v0)

5
lib/piir.py
View File

@ -57,6 +57,7 @@ halfband_c['olli'] = [
0.9987488452737**2,
]
class Halfband:
def __init__(self, c='olli'):
self.x = np.zeros(4)
@ -94,10 +95,10 @@ class Halfband:
sign = 1
if mode == 'hilbert':
#y[n] = c*(x[n] + y[n-2]) - x[n-2]
# y[n] = c*(x[n] + y[n-2]) - x[n-2]
pass
elif mode == 'filter':
#y[n] = c*(x[n] - y[n-2]) + x[n-2]
# y[n] = c*(x[n] - y[n-2]) + x[n-2]
sign = -1
in2 = self.x[2]

23
lib/planes.py
View File

@ -2,6 +2,7 @@ from . import tau
import numpy as np
# implements the modified bilinear transform:
# s <- 1/tan(w0/2)*(1 - z^-1)/(1 + z^-1)
# this requires the s-plane coefficients to be frequency-normalized,
@ -20,6 +21,7 @@ def zcgen_py(n, d):
zcs[i] += zcs[i - 1]
return zcs
def zcgen_sym(n, d):
import sympy as sym
z = sym.symbols('z')
@ -27,6 +29,7 @@ def zcgen_sym(n, d):
coeffs = expr.equals(1) and [1] or expr.as_poly().all_coeffs()
return coeffs[::-1]
def s2z_two(b, a, fc, srate, gain=1):
"""
converts s-plane coefficients to z-plane for digital usage.
@ -40,17 +43,18 @@ def s2z_two(b, a, fc, srate, gain=1):
cw = np.cos(w0)
sw = np.sin(w0)
zb = np.array((
b[2]*(1 - cw) + b[0]*(1 + cw) + b[1]*sw,
2*(b[2]*(1 - cw) - b[0]*(1 + cw)),
b[2]*(1 - cw) + b[0]*(1 + cw) - b[1]*sw,
(b[2]*(1 - cw) + b[0]*(1 + cw) + b[1]*sw),
(b[2]*(1 - cw) - b[0]*(1 + cw)) * 2,
(b[2]*(1 - cw) + b[0]*(1 + cw) - b[1]*sw),
))
za = np.array((
a[2]*(1 - cw) + a[0]*(1 + cw) + a[1]*sw,
2*(a[2]*(1 - cw) - a[0]*(1 + cw)),
a[2]*(1 - cw) + a[0]*(1 + cw) - a[1]*sw,
(a[2]*(1 - cw) + a[0]*(1 + cw) + a[1]*sw),
(a[2]*(1 - cw) - a[0]*(1 + cw)) * 2,
(a[2]*(1 - cw) + a[0]*(1 + cw) - a[1]*sw),
))
return zb*gain, za
def s2z1(w0, s, d):
"""
s: array of s-plane coefficients (num OR den, not both)
@ -67,6 +71,7 @@ def s2z1(w0, s, d):
y[i] += trig*zcs[i]*s[n]
return y
def s2z_any(b, a, fc, srate, gain=1, d=-1):
"""
converts s-plane coefficients to z-plane for digital usage.
@ -83,8 +88,10 @@ def s2z_any(b, a, fc, srate, gain=1, d=-1):
za = s2z1(w0, sa, cs - 1)
return zb*gain, za
# set our preference. zcgen_py is 1000+ times faster than zcgen_sym
# set our preference. zcgen_py is 1000+ times faster than zcgen_sym.
zcgen = zcgen_py
# s2z_any is only ~2.4 times slower than s2z_two and allows for filters of any degree
# s2z_any is only ~2.4 times slower than s2z_two
# and allows for filters of any degree.
s2z = s2z_any

12
lib/plot.py
View File

@ -1,6 +1,7 @@
import matplotlib.pyplot as plt
from matplotlib import ticker
def response_setup(ax, ymin=-24, ymax=24, yL=ticker.AutoMinorLocator(3)):
ax.set_xlim(20, 20000)
ax.set_ylim(ymin, ymax)
@ -10,6 +11,7 @@ def response_setup(ax, ymin=-24, ymax=24, yL=ticker.AutoMinorLocator(3)):
ax.set_xlabel('frequency (Hz)')
ax.set_ylabel('magnitude (dB)')
def phase_response_setup(ax, div=12, yL=ticker.AutoMinorLocator(2)):
ax.set_xlim(20, 20000)
ax.set_ylim(-180, 180)
@ -19,22 +21,26 @@ def phase_response_setup(ax, div=12, yL=ticker.AutoMinorLocator(2)):
ax.set_xlabel('frequency (Hz)')
ax.set_ylabel('phase (degrees)')
def cleanplot():
fig, ax = plt.subplots()
ax.set_axis_off()
ax.set_position([0,0,1,1])
ax.set_position([0, 0, 1, 1])
return fig, ax
def new_response(*args, **kwargs):
fig, ax = plt.subplots()
response_setup(ax, *args, **kwargs)
return fig, ax
def new_phase_response(*args, **kwargs):
fig, ax = plt.subplots()
phase_response_setup(ax, *args, **kwargs)
return fig, ax
def new_bode(magnitude_offset=0):
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
@ -51,8 +57,8 @@ def new_bode(magnitude_offset=0):
for tl in ax2.get_yticklabels():
tl.set_color(cc[1])
#ax1.hlines(0, 20, 40, linewidth=0.5, color=cc[0])
#ax2.hlines(0, 10000, 20000, linewidth=0.5, color=cc[1])
# ax1.hlines(0, 20, 40, linewidth=0.5, color=cc[0])
# ax2.hlines(0, 10000, 20000, linewidth=0.5, color=cc[1])
# share color cycles to prevent color re-use
ax2._get_lines.prop_cycler = ax1._get_lines.prop_cycler

16
lib/plotwav.py
View File

@ -5,7 +5,9 @@ from . import new_response, magnitude_x, convolve_each, monoize, count_channels
import numpy as np
def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192, window=0, raw=False, **kwargs):
def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192,
window=0, raw=False, **kwargs):
sm = monoize(s)
xs_raw = magnitude_x(srate, size)
@ -16,11 +18,13 @@ def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192, window=0, raw=F
xs, ys = smoothfft(xs_raw, ys_raw, bw=bw)
if ax:
if raw: ax.semilogx(xs_raw, ys_raw, **kwargs)
if raw:
ax.semilogx(xs_raw, ys_raw, **kwargs)
ax.semilogx(xs, ys, **kwargs)
return xs, ys
def plotwavinternal(sm, ss, srate, bw=1, size=8192, smoother=smoothfft2):
xs_raw = magnitude_x(srate, size)
ys_raw_m = averfft(sm, size=size)
@ -38,6 +42,7 @@ def plotwavinternal(sm, ss, srate, bw=1, size=8192, smoother=smoothfft2):
return xs, ys_m, ys_s
def plotwav2(fn, bw=1, size=8192, fix=False,
smoother=smoothfft2, **kwargs):
s, srate = wav_read(fn)
@ -64,19 +69,22 @@ def plotwav2(fn, bw=1, size=8192, fix=False,
sf = np.array((smf + ssf, smf - ssf)).T
import ewave
with ewave.open(fno, 'w', sampling_rate=srate, nchannels=count_channels(sf)) as f:
with ewave.open(fno, 'w', sampling_rate=srate,
nchannels=count_channels(sf)) as f:
f.write(sf)
print('wrote '+fno)
return xs, ys_m, ys_s
def pw2(fn, label=None, bw=1/6, **kwargs):
fno = fn[:-4]+"-proc.wav"
xs, ys_m, ys_s = plotwav2(fn, fix=True, bw=bw, **kwargs)
xs, ys_m, ys_s = plotwav2(fno, fix=False, bw=bw, **kwargs)
fig, ax = new_response(-18, 18)
ax.set_title('averaged magnitudes of normalized songs with tilt and smoothing')
ax.set_title(
'averaged magnitudes of normalized songs with tilt and smoothing')
label = label or fn
ax.semilogx(xs, ys_m + 0, label=label+' (mid)')
ax.semilogx(xs, ys_s + 9, label=label+' (side)')

9
lib/smoothfft.py
View File

@ -1,6 +1,7 @@
from . import xsp, lament
import numpy as np
def smoothfft(xs, ys, bw=1, precision=512):
"""performs log-lin smoothing on magnitude data,
generally from the output of averfft."""
@ -18,6 +19,7 @@ def smoothfft(xs, ys, bw=1, precision=512):
ys2[i] = np.sum(ys*window/wsum)
return xs2, ys2
def smoothfft2(xs, ys, bw=1, precision=512, compensate=True):
"""performs log-lin smoothing on magnitude data,
generally from the output of averfft."""
@ -28,10 +30,11 @@ def smoothfft2(xs, ys, bw=1, precision=512, compensate=True):
for i, x in enumerate(xs):
# before optimizations: dist = np.abs(np.log2(xs2/(x + 1e-35)))/bw
dist = np.abs(log2_xs2 - np.log2(x + 1e-35))/bw
#window = np.maximum(0, 1 - dist) # triangle window
window = np.exp(-dist**2/(0.5/2)) # gaussian function (non-truncated)
# window = np.maximum(0, 1 - dist) # triangle window
window = np.exp(-dist**2/(0.5/2)) # gaussian function (non-truncated)
ys2 += ys[i]*window
if compensate:
_, temp = smoothfft2(xs, np.ones(len(xs)), bw=bw, precision=precision, compensate=False)
_, temp = smoothfft2(xs, np.ones(len(xs)),
bw=bw, precision=precision, compensate=False)
ys2 /= temp
return xs2, ys2

77
lib/svf.py
View File

@ -2,8 +2,10 @@ from . import tau, unwarp
import numpy as np
# via http://nbviewer.ipython.org/urls/music-synthesizer-for-android.googlecode.com/git/lab/Second%20order%20sections%20in%20matrix%20form.ipynb
def svf_core(w0, Q, m, shelfA=1, gain=1):
# via:
# http://nbviewer.ipython.org/urls/music-synthesizer-for-android.googlecode.com/git/lab/Second%20order%20sections%20in%20matrix%20form.ipynb
# TODO: implement constant gain parameter
g = unwarp(w0)*shelfA
a1 = 1/(1 + g*(g + 1/Q))
@ -17,52 +19,79 @@ def svf_core(w0, Q, m, shelfA=1, gain=1):
C = v0*m[0] + v1*m[1] + v2*m[2]
return A, B, C
LP2S = lambda A, Q: (Q, [0, 0, 1], 1)
BP2S = lambda A, Q: (Q, [0, 1, 0], 1)
HP2S = lambda A, Q: (Q, [1, -1/Q, -1], 1)
#AP2S = lambda A, Q:
#BP2aS = lambda A, Q:
#BP2bS = lambda A, Q:
NO2S = lambda A, Q: (Q, [1, -1/Q, 0], 1)
PE2S = lambda A, Q: (Q*A, [1, (A*A - 1)/A/Q, 0], 1)
LS2S = lambda A, Q: (Q, [1, (A - 1)/Q, A*A - 1], 1/np.sqrt(A))
HS2S = lambda A, Q: (Q, [A*A, (1 - A)*A/Q, 1 - A*A], np.sqrt(A))
# actual peaking filter (not a bell?)
#PE2S = lambda A, Q: ([1, -1/Q, -2], 1)
# original uncompensated
#PE2S = lambda A, Q: (Q, [1, (A*A - 1)/Q, 0], 1)
#LS2S = lambda A, Q: (Q, [1, (A - 1)/Q, A*A - 1], 1/np.sqrt(A))
#HS2S = lambda A, Q: (Q, [A*A, (A - A*A)/Q, 1 - A*A], 1/np.sqrt(A))
def LP2S(A, Q):
return (Q, [0, 0, 1], 1)
def BP2S(A, Q):
return (Q, [0, 1, 0], 1)
def HP2S(A, Q):
return (Q, [1, -1/Q, -1], 1)
# TODO: AP2S
# TODO: BP2aS
# TODO: BP2bS
def NO2S(A, Q):
return (Q, [1, -1/Q, 0], 1)
def PE2S(A, Q):
return (Q*A, [1, (A*A - 1)/A/Q, 0], 1)
def LS2S(A, Q):
return (Q, [1, (A - 1)/Q, A*A - 1], 1/np.sqrt(A))
def HS2S(A, Q):
return (Q, [A*A, (1 - A)*A/Q, 1 - A*A], np.sqrt(A))
# actual peaking filter: (not a bell?)
# PE2S = lambda A, Q: ([1, -1/Q, -2], 1)
# original uncompensated:
# PE2S = lambda A, Q: (Q, [1, (A*A - 1)/Q, 0], 1)
# LS2S = lambda A, Q: (Q, [1, (A - 1)/Q, A*A - 1], 1/np.sqrt(A))
# HS2S = lambda A, Q: (Q, [A*A, (A - A*A)/Q, 1 - A*A], 1/np.sqrt(A))
def gen_filters_svf(cascade, srate):
return [
svf_core(tau*f[0]/srate, *f[1], gain=10**(f[2]/20)) for f in cascade
]
gen_filters_svf = lambda cascade, srate: [
svf_core(tau*f[0]/srate, *f[1], gain=10**(f[2]/20)) for f in cascade
]
def svf_run(svf, xs):
A, B, C = svf
result = []
y = np.zeros(2) # filter memory
y = np.zeros(2) # filter memory
for x in xs:
result.append(np.dot(C, np.concatenate([[x], y])))
y = B*x + np.dot(A, y)
return np.array(result)
def svf_mat(svf):
A, B, C = svf
AA = np.dot(A, A)
AB = np.dot(A, B)
CA = np.dot(C[1:], A)
cb = np.dot(C[1:], B)
return np.array([[ C[0], 0, C[1], C[2]],
[ cb, C[0], CA[0], CA[1]],
return np.array([[C[0], 0, C[1], C[2]],
[cb, C[0], CA[0], CA[1]],
[AB[0], B[0], AA[0][0], AA[0][1]],
[AB[1], B[1], AA[1][0], AA[1][1]]])
def svf_run4(mat, xs):
assert(len(xs) % 2 == 0)
out = np.zeros(len(xs))
y = np.zeros(4) # filter memory
y = np.zeros(4) # filter memory
for i in range(0, len(xs), 2):
y[0:2] = xs[i:i+2]
y = np.dot(mat, y)

9
lib/sweeps.py
View File

@ -2,18 +2,19 @@ from . import tau
import numpy as np
def sweep(amp, length, begin=20, end=20480, method='linear'):
method = method or 'linear'
xs = np.arange(length)/length
if method in ('linear', 'quadratic', 'logarithmic', 'hyperbolic'):
ys = amp*sig.chirp(xs, begin, 1, end, method=method)
elif method is 'sinesweep':
ang = lambda f: tau*f
# because xs ranges from 0:1, length is 1 and has been simplified out
domain = np.log(ang(end)/ang(begin))
ys = amp*np.sin(ang(begin)/domain*(np.exp(xs*domain) - 1))
domain = np.log((tau * end)/(tau * begin))
ys = amp*np.sin((tau * begin)/domain*(np.exp(xs*domain) - 1))
return ys
def tsp(N, m=0.5):
"""
OATSP(Optimized Aoshima's Time-Stretched Pulse) generator
@ -36,7 +37,7 @@ def tsp(N, m=0.5):
if N < 0:
raise Exception("The number of length must be the positive number")
NN = 2**np.floor(np.log2(N)) # nearest
NN = 2**np.floor(np.log2(N)) # nearest
NN2 = NN//2
M = np.round(NN2*m)

81
lib/util.py
View File

@ -2,33 +2,73 @@ import sys
import numpy as np
import scipy.signal as sig
dummy = lambda *args, **kwargs: None
lament = lambda *args, **kwargs: print(*args, file=sys.stderr, **kwargs)
toLK = lambda x: -0.691 + 10*np.log10(x)
isqrt2 = 1/np.sqrt(2)
toQ = lambda bw: isqrt2/bw
toA = lambda db: 10**(db/40)
tau = 2*np.pi
unwarp = lambda w: np.tan(w/2)
warp = lambda w: np.arctan(w)*2
ceil2 = lambda x: np.power(2, np.ceil(np.log2(x)))
pad2 = lambda x: np.r_[x, np.zeros(ceil2(len(x)) - len(x))]
rfft = lambda src, size: np.fft.rfft(src, size*2)
magnitude = lambda src, size: 10*np.log10(np.abs(rfft(src, size))**2)[0:size]
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):
return np.power(2, np.ceil(np.log2(x)))
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]
# x axis for plotting above magnitude
magnitude_x = lambda srate, size: np.arange(0, srate/2, srate/2/size)
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
degrees_clamped = lambda x: ((x*180/np.pi + 180) % 360) - 180
def xsp(precision=4096):
"""create #precision log-spaced points from 20 Hz (inclusive) to 20480 Hz (exclusive)"""
xs = np.arange(0,precision)/precision
"""
create #precision log-spaced points from
20 Hz (inclusive) to 20480 Hz (exclusive)
"""
xs = np.arange(0, precision)/precision
return 20*1024**xs
def blocks(a, step, size=None):
"""break an iterable into chunks"""
if size is None:
@ -39,14 +79,18 @@ def blocks(a, step, size=None):
break
yield a[start:end]
def convolve_each(s, fir, mode='same', axis=0):
return np.apply_along_axis(lambda s: sig.fftconvolve(s, fir, mode), axis, s)
return np.apply_along_axis(
lambda s: sig.fftconvolve(s, fir, mode), axis, s)
def count_channels(s):
if s.ndim < 2:
return 1
return s.shape[1]
def monoize(s):
"""mixes an n-channel signal down to one channel.
technically, it averages a 2D array to be 1D.
@ -56,6 +100,7 @@ def monoize(s):
s = np.average(s, axis=1)
return s
def div0(a, b):
"""division, whereby division by zero equals zero"""
# http://stackoverflow.com/a/35696047
@ -63,5 +108,5 @@ def div0(a, b):
b = np.asanyarray(b)
with np.errstate(divide='ignore', invalid='ignore'):
c = np.true_divide(a, b)
c[~np.isfinite(c)] = 0 # -inf inf NaN
c[~np.isfinite(c)] = 0 # -inf inf NaN
return c

7
lib/wav.py
View File

@ -1,15 +1,18 @@
import numpy as np
from .util import lament, count_channels
def wav_smart_read(fn):
lament('wav_smart_read(): DEPRECATED; use wav_read instead.')
import scipy.io.wavfile as wav # don't use this, it fails to load good files
# don't use this, it fails to load good files.
import scipy.io.wavfile as wav
srate, s = wav.read(fn)
if s.dtype != np.float64:
bits = s.dtype.itemsize*8
s = np.asfarray(s)/2**(bits - 1)
return srate, s
def wav_smart_write(fn, srate, s):
lament('wav_smart_write(): DEPRECATED; use wav_write instead.')
import scipy.io.wavfile as wav
@ -18,6 +21,7 @@ def wav_smart_write(fn, srate, s):
si += np.clip(s*2**(bits - 1), -32768, 32767)
wav.write(fn, srate, si)
def wav_read(fn):
import ewave
with ewave.open(fn) as f:
@ -30,6 +34,7 @@ def wav_read(fn):
s = np.asfarray(s)/2**(bits - 1)
return s, srate
def wav_write(fn, s, srate, dtype='h'):
import ewave
if dtype in ('b', 'h', 'i', 'l') and np.max(np.abs(s)) > 1:

36
lib/windowing.py
View File

@ -1,5 +1,6 @@
import numpy as np
def _deco_win(f):
# gives scipy compatibility
def deco(N, *args, sym=True, **kwargs):
@ -18,37 +19,46 @@ def _deco_win(f):
return w
return deco
def _gen_hamming(*a):
L = len(a)
a += (0, 0, 0, 0, 0) # pad so orders definition doesn't error
a += (0, 0, 0, 0, 0) # pad so orders definition doesn't error
orders = [
lambda fac: 0,
lambda fac: a[0],
lambda fac: a[0] - a[1]*np.cos(fac),
lambda fac: a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac),
lambda fac: a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac) - a[3]*np.cos(3*fac),
lambda fac: a[0] - a[1]*np.cos(fac) + a[2]*np.cos(2*fac) - a[3]*np.cos(3*fac) + a[4]*np.cos(4*fac),
lambda fac: a[0] - a[1]*np.cos(1*fac),
lambda fac: a[0] - a[1]*np.cos(1*fac) + a[2]*np.cos(2*fac),
lambda fac: a[0] - a[1]*np.cos(1*fac) + a[2]*np.cos(2*fac)
- a[3]*np.cos(3*fac),
lambda fac: a[0] - a[1]*np.cos(1*fac) + a[2]*np.cos(2*fac)
- a[3]*np.cos(3*fac) + a[4]*np.cos(4*fac),
]
f = orders[L]
return lambda N: f(np.arange(0, N)*2*np.pi/(N - 1))
def _normalize(*args):
a = np.asfarray(args)
return a/np.sum(a)
_h = lambda *args: _deco_win(_gen_hamming(*args))
def _h(*args):
return _deco_win(_gen_hamming(*args))
blackman_inexact = _h(0.42, 0.5, 0.08)
blackman = _h(0.42659, 0.49656, 0.076849)
blackman_nuttall = _h(0.3635819, 0.4891775, 0.1365995, 0.0106411)
blackman_harris = _h(0.35875, 0.48829, 0.14128, 0.01168)
nuttall = _h(0.355768, 0.487396, 0.144232, 0.012604)
flattop = _h(*_normalize(1, 1.93, 1.29, 0.388, 0.028)) # FTSRS
#flattop_weird = _h(*_normalize(1, 1.93, 1.29, 0.388, 0.032)) # ??? wtf
flattop_weird = _h(0.2156, 0.4160, 0.2781, 0.0836, 0.0069) # ??? scipy crap
flattop = _h(*_normalize(1, 1.93, 1.29, 0.388, 0.028)) # FTSRS
# flattop_weird = _h(*_normalize(1, 1.93, 1.29, 0.388, 0.032)) # ??? wtf
flattop_weird = _h(0.2156, 0.4160, 0.2781, 0.0836, 0.0069) # ??? scipy crap
hann = _h(0.5, 0.5)
hamming_inexact = _h(0.54, 0.46)
hamming = _h(0.53836, 0.46164)
@_deco_win
def triangular(N):
if N % 2 == 0:
@ -56,6 +66,7 @@ def triangular(N):
else:
return 1 - np.abs(2*(np.arange(N) + 1)/(N + 1) - 1)
@_deco_win
def parzen(N):
odd = N % 2
@ -71,17 +82,22 @@ def parzen(N):
else:
return np.r_[outer[::-1], center[::-1], center[1:], outer]
@_deco_win
def cosine(N):
return np.sin(np.pi*(np.arange(N) + 0.5)/N)
@_deco_win
def welch(N):
return 1 - (2*np.arange(N)/(N - 1) - 1)**2
# TODO: rename or something
@_deco_win
def sinc(N):
return np.sinc((np.arange(N) - (N - 1)/2)/2)
winmod = lambda f: lambda N: f(N + 2)[1:-1]
def winmod(f):
return lambda N: f(N + 2)[1:-1]