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

update 20

Browse Source
This commit is contained in:
Connor Olding 2015-11-10 04:04:41 -08:00
parent 96f3250148
commit 77f38773d0
6 changed files with 83 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
lib/__init__.py
View File

@ -115,8 +115,8 @@ def firize(xs, ys, n=4096, srate=44100, ax=None):
xf = xs/srate*2
yg = 10**(ys/20)
xf = np.hstack((0, xf, 1))
yg = np.hstack((0, yg, yg[-1]))
xf = np.r_[0, xf, 1]
yg = np.r_[0, yg, yg[-1]]
b = sig.firwin2(n, xf, yg, antisymmetric=True)

8
lib/cepstrum.py
View File

@ -26,11 +26,11 @@ def fold(r):
elif n % 2 == 1:
nt = (n + 1)//2
rf = r[1:nt] + conj(r[-1:nt-1:-1])
rw = np.hstack((r[0], rf, np.zeros(n-nt)))
rw = np.r_[r[0], rf, np.zeros(n-nt)]
else:
nt = n//2
rf = np.hstack((r[1:nt], 0)) + np.conj(r[-1:nt-1:-1])
rw = np.hstack((r[0], rf, np.zeros(n-nt-1)))
rf = np.r_[r[1:nt], 0] + np.conj(r[-1:nt-1:-1])
rw = np.r_[r[0], rf, np.zeros(n-nt-1)]
return rw
def minphase(s, pad=True, os=False):
@ -39,7 +39,7 @@ def minphase(s, pad=True, os=False):
if pad:
s = pad2(s)
if os:
s = np.hstack((s, np.zeros(len(s))))
s = np.r_[s, np.zeros(len(s))]
cepstrum = np.fft.ifft(np.log(clipdb(np.fft.fft(s), -100)))
signal = np.real(np.fft.ifft(np.exp(np.fft.fft(fold(cepstrum)))))
if os:

2
lib/fft.py
View File

@ -19,7 +19,7 @@ def magnitudes(s, size=8192):
# blindly pad with zeros for friendlier ffts and overlapping
z = np.zeros(size)
s = np.hstack((s, z))
s = np.r_[s, z]
win_size = size

6
lib/sweeps.py
View File

@ -45,10 +45,10 @@ def tsp(N, m=0.5):
j = np.complex(0, 1)
H = np.exp(j*4*M*np.pi*nn2/NN**2)
H2 = np.hstack([H, np.conj(H[1:NN2][::-1])])
H2 = np.r_[H, np.conj(H[1:NN2][::-1])]
x = np.fft.ifft(H2)
x = np.hstack([x[NN2 - M:NN + 1], x[0:NN2 - M + 1]])
x = np.hstack([x.real, np.zeros(1, N - NN)])
x = np.r_[x[NN2 - M:NN + 1], x[0:NN2 - M + 1]]
x = np.r_[x.real, np.zeros(1, N - NN)]
return x

2
lib/util.py
View File

@ -15,7 +15,7 @@ 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.hstack((x, np.zeros(ceil2(len(x)) - len(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]

98
lib/windowing.py
View File

@ -1,41 +1,87 @@
from . import tau
import numpy as np
def gen_hamm(a0, a1=0, a2=0, a3=0, a4=0):
form = lambda x: a0 - a1*np.cos(x) + a2*np.cos(2*x) - a3*np.cos(3*x) + a4*np.cos(4*x)
dc = form(np.pi)
def f(N):
a = tau*np.arange(N)/(N - 1)
return form(a)/dc
return f
def _deco_win(f):
# gives scipy compatibility
def deco(N, *args, sym=True, **kwargs):
if N < 1:
return np.array([])
if N == 1:
return np.ones(1)
odd = N % 2
if not sym and not odd:
N = N + 1
# TODO: rename or something
sinc = lambda N: np.sinc((np.arange(N) - (N - 1)/2)/2)
w = f(N, *args, **kwargs)
triangular = lambda N: 1 - np.abs(2*np.arange(N)/(N - 1) - 1)
welch = lambda N: 1 - (2*np.arange(N)/(N - 1) - 1)**2
cosine = lambda N: np.sin(np.pi*np.arange(N)/(N - 1))
if not sym and not odd:
return w[:-1]
return w
return deco
hann = gen_hamm(0.5, 0.5)
hamming = gen_hamm(0.53836, 0.46164)
blackman = gen_hamm(0.42659, 0.49656, 0.076849)
blackman_harris = gen_hamm(0.35875, 0.48829, 0.14128, 0.01168)
blackman_nutall = gen_hamm(0.3635819, 0.4891775, 0.1365995, 0.0106411)
nutall = gen_hamm(0.355768, 0.487396, 0.144232, 0.012604)
# specifically the Stanford Research flat-top window (FTSRS)
flattop = gen_hamm(1, 1.93, 1.29, 0.388, 0.028)
def _gen_hamming(*a):
L = len(a)
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),
]
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))
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
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:
return 1 - np.abs((2*np.arange(N) + 1)/N - 1)
else:
return 1 - np.abs(2*(np.arange(N) + 1)/(N + 1) - 1)
@_deco_win
def parzen(N):
# TODO: verify this is accurate. the code here is kinda sketchy
odd = N % 2
n = np.arange(N/2)
if not odd:
n += 0.5
center = 1 - 6*(2*n/N)**2*(1 - 2*n/N)
outer = 2*(1 - 2*n/N)**3
center = center[:len(center)//2]
outer = outer[len(outer)//2:]
if N % 2 == 0:
return np.hstack((outer[::-1], center[::-1], center, outer))
if not odd:
return np.r_[outer[::-1], center[::-1], center, outer]
else:
return np.hstack((outer[::-1], center[::-1], center[1:], outer))
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)
# some windows reach 0 at either side, so this can avoid that
winmod = lambda f: lambda N: f(N + 2)[1:-1]