diff --git a/lib/__init__.py b/lib/__init__.py index 1cfc20a..4405206 100644 --- a/lib/__init__.py +++ b/lib/__init__.py @@ -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) diff --git a/lib/cepstrum.py b/lib/cepstrum.py index 1e00dac..3fb1095 100644 --- a/lib/cepstrum.py +++ b/lib/cepstrum.py @@ -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: diff --git a/lib/fft.py b/lib/fft.py index a953953..90fd06b 100644 --- a/lib/fft.py +++ b/lib/fft.py @@ -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 diff --git a/lib/sweeps.py b/lib/sweeps.py index 1e6bbdf..dbd2ce4 100644 --- a/lib/sweeps.py +++ b/lib/sweeps.py @@ -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 diff --git a/lib/util.py b/lib/util.py index 7feb21c..c123b6f 100644 --- a/lib/util.py +++ b/lib/util.py @@ -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] diff --git a/lib/windowing.py b/lib/windowing.py index 4ee2926..a42c1e6 100644 --- a/lib/windowing.py +++ b/lib/windowing.py @@ -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]