update 49

lib/__init__.py

@@ -1,5 +1,6 @@
 from .util import *
 from .bq import *
+from .svf import *
 from .data import *
 from .nsf import *
 from .sweeps import *

lib/bs.py

@@ -22,6 +22,7 @@ def BS1770_3(s, srate, filters=None, window=0.4, overlap=0.75,
     means = np.array([
         np.sum(np.mean(b**2, axis=0)) for b in blocks(sf, stepsize, blocksize)
     ])
+    means[means < 1e-10] = 1e-10  # clip at -100 dB to avoid negative infinity
     LKs = toLK(means)
 
     gated = LKs > -absolute_gate

lib/cepstrum.py

@@ -10,7 +10,7 @@ def ifcs(s):  # inverted fast cepstrum
     return np.fft.fft(np.exp(np.fft.ifft(s)))
 
 
-def mcs(s):  # magnitude
+def mcs(s):  # magnitude (actually power)
     return (np.abs(np.fft.ifft(np.log(np.abs(np.fft.fft(s))**2)))**2
             )[:len(s)//2]
 

lib/fft.py

@@ -1,6 +1,8 @@
 import numpy as np
 import scipy.signal as sig
 
+from .util import lament
+
 
 def magnitudes_window_setup(s, size=8192, overlap=0.661):
     # note: the default overlap value is only
@@ -11,7 +13,10 @@ def magnitudes_window_setup(s, size=8192, overlap=0.661):
     return step, segs
 
 
-def magnitudes(s, size=8192):
+def magnitudes(s, size=8192, broken=True):
+    if broken:
+        lament("magnitudes(broken=True): DEPRECATED")
+
     step, segs = magnitudes_window_setup(s, size)
 
     L = s.shape[0]
@@ -25,24 +30,31 @@ def magnitudes(s, size=8192):
     win = sig.blackmanharris(win_size)
     win /= np.sqrt(np.sum(np.square(win)))
 
-    count = 0
     for i in range(0, L - 1, int(step)):
         windowed = s[i:i+win_size]*win
+        if broken:
             power = np.abs(np.fft.rfft(windowed, 2 * size))**2
             # this scraps the nyquist value to get exactly 'size' outputs
-        yield power[0:size]
+            yield power[:size]
-        count += 1
+        else:
-
+            power = np.abs(np.fft.rfft(windowed, size))**2
-    # assert(segs == count)  # this is probably no good in a generator
+            # this scraps the 0 Hz value to get exactly size//2 outputs
+            yield power[1:]
 
 
-def averfft(s, size=8192):
+def averfft(s, size=8192, broken=True):
     """calculates frequency magnitudes by fft and averages them together."""
     step, segs = magnitudes_window_setup(s, size)
 
+    if broken:
+        lament("averfft(broken=True): DEPRECATED")
         avg = np.zeros(size)
         for power in magnitudes(s, size):
             avg += power/segs
+    else:
+        avg = np.zeros(size//2)
+        for power in magnitudes(s, size):
+            avg += power/segs
 
     avg_db = 10*np.log10(avg)
     return avg_db

lib/plotwav.py

@@ -1,6 +1,7 @@
 # this is a bunch of crap that should really be reduced to one or two functions
 
-from . import wav_read, normalize, averfft, tilter2, smoothfft2, firize
+from . import wav_read, wav_write
+from . import normalize, averfft, tilter2, smoothfft4, firize
 from . import new_response, magnitude_x, convolve_each, monoize, count_channels
 
 import numpy as np
@@ -12,10 +13,8 @@ def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192,
 
     xs_raw = magnitude_x(srate, size)
     ys_raw = averfft(sm, size=size, mode=window)
-
     ys_raw -= tilter2(xs_raw, tilt)
-
+    xs, ys = smoothfft4(ys_raw, bw)
-    xs, ys = smoothfft(xs_raw, ys_raw, bw=bw)
 
     if ax:
         if raw:
@@ -25,7 +24,7 @@ def plotfftsmooth(s, srate, ax=None, bw=1, tilt=None, size=8192,
     return xs, ys
 
 
-def plotwavinternal(sm, ss, srate, bw=1, size=8192, smoother=smoothfft2):
+def plotwavinternal(sm, ss, srate, bw=1, size=8192):
     xs_raw = magnitude_x(srate, size)
     ys_raw_m = averfft(sm, size=size)
     ys_raw_s = averfft(ss, size=size)
@@ -37,24 +36,23 @@ def plotwavinternal(sm, ss, srate, bw=1, size=8192, smoother=smoothfft2):
     if bw <= 0:
         return xs_raw, xs_raw_m, xs_raw_s
 
-    xs, ys_m = smoother(xs_raw, ys_raw_m, bw=bw)
+    xs, ys_m = smoothfft4(ys_raw_m, bw=bw, srate=srate)
-    xs, ys_s = smoother(xs_raw, ys_raw_s, bw=bw)
+    xs, ys_s = smoothfft4(ys_raw_s, bw=bw, srate=srate)
 
     return xs, ys_m, ys_s
 
 
-def plotwav2(fn, bw=1, size=8192, fix=False,
+def plotwav2(fn, bw=1, size=8192, fix=False, **kwargs):
-             smoother=smoothfft2, **kwargs):
     s, srate = wav_read(fn)
 
     s, rms = normalize(s, srate)
     sm = monoize(s)
-    if s.shape[1] == 2:
+    if s.ndim > 1 and s.shape[1] == 2:
         ss = monoize(s*np.array((1, -1)))
     else:
         ss = np.zeros(len(s))
 
-    xs, ys_m, ys_s = plotwavinternal(sm, ss, srate, bw, size, smoother)
+    xs, ys_m, ys_s = plotwavinternal(sm, ss, srate, bw, size)
 
     side_gain = np.average(ys_s) - np.average(ys_m)
 
@@ -66,12 +64,9 @@ def plotwav2(fn, bw=1, size=8192, fix=False,
         smf = convolve_each(sm/8, fir_m, mode='same')
         ssf = convolve_each(ss/8, fir_s, mode='same')
         ssf *= 10**(side_gain/20)
-        sf = np.array((smf + ssf, smf - ssf)).T
+        sf = np.c_[smf + ssf, smf - ssf]
 
-        import ewave
+        wav_write(fno, sf, srate, dtype='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
@@ -89,3 +84,4 @@ def pw2(fn, label=None, bw=1/6, **kwargs):
     ax.semilogx(xs, ys_m + 0, label=label+' (mid)')
     ax.semilogx(xs, ys_s + 9, label=label+' (side)')
     ax.legend(loc=8)
+    return fig, ax

lib/smoothfft.py

@@ -1,11 +1,11 @@
-from . import xsp, lament
+from . import xsp, lament, ceil2
 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."""
-    lament("smoothfft(): DEPRECATED; use smoothfft2 instead.")
+    lament("smoothfft(): DEPRECATED; use smoothfft4 instead.")
     xs2 = xsp(precision)
     ys2 = np.zeros(precision)
     log_xs = np.log(xs)
@@ -23,7 +23,7 @@ def smoothfft(xs, ys, bw=1, precision=512):
 def smoothfft2(xs, ys, bw=1, precision=512, compensate=True):
     """performs log-lin smoothing on magnitude data,
     generally from the output of averfft."""
-    # this is probably implementable with FFTs now that i think about it
+    lament('smoothfft2: DEPRECATED; use smoothfft4 instead.')
     xs2 = xsp(precision)
     ys2 = np.zeros(precision)
     log2_xs2 = np.log2(xs2)
@@ -41,6 +41,7 @@ def smoothfft2(xs, ys, bw=1, precision=512, compensate=True):
 
 
 def smoothfft_setup(size, precision=512, bw=1/6):
+    lament('smoothfft_setup(): DEPRECATED; use smoothfft_setup2 instead.')
     dotme = np.zeros((size, precision))
 
     xs = np.arange(0, 1, 1/size)
@@ -61,8 +62,41 @@ def smoothfft_setup(size, precision=512, bw=1/6):
 
 def smoothfft3(ys, bw=1, precision=512, srate=None):
     """performs log-lin smoothing on magnitude data"""
+    lament('smoothfft3(): DEPRECATED; use smoothfft4 instead.')
     xs2, dotme = smoothfft_setup(len(ys), precision, bw)
     if srate is None:
         return xs2, ys @ dotme
     else:
         return xs2 * (srate / 2), ys @ dotme
+
+
+def smoothfft_setup2(size, precision=512, bw=1/6):
+    # tweaked/fixed to drop 0 Hz
+    size -= size % 2
+    assert size == ceil2(size), size
+
+    dotme = np.zeros((size, precision))
+    xs = np.arange(1, size + 1) / size
+    xs2 = np.logspace(-np.log2(size), 0, precision, base=2)
+    comp = np.zeros(precision)
+
+    log2_xs2 = np.log2(xs2)
+    for i, x in enumerate(xs):
+        dist = np.abs(log2_xs2 - np.log2(x)) / bw
+        window = np.exp(-dist**2 * 4)  # gaussian (untruncated)
+        comp += window
+        dotme[i] = window
+
+    dotme /= comp
+    return xs2, dotme
+
+
+def smoothfft4(ys, bw=1, precision=512, srate=None):
+    # tweaked/fixed to drop 0 Hz
+    if len(ys) % 2 == 1:
+        ys = ys[1:]
+    xs2, dotme = smoothfft_setup2(len(ys), precision, bw)
+    if srate is None:
+        return xs2, ys @ dotme
+    else:
+        return xs2 * (srate / 2), ys @ dotme

lib/util.py

@@ -45,13 +45,21 @@ def pad2(x):
     return np.r_[x, np.zeros(ceil2(len(x)) - len(x), x.dtype)]
 
 
-def magnitude(src, size):
+def magnitude(src, size, broken=True):
-    return 10*np.log10(np.abs(np.fft.rfft(src, 2 * size))**2)[0:size]
+    if broken:
+        lament("magnitude(broken=True): DEPRECATED")
+        return 10*np.log10(np.abs(np.fft.rfft(src, 2 * size))**2)[:size]
+    else:
+        return 10*np.log10(np.abs(np.fft.rfft(src, size))**2)[1:]
 
 
 # x axis for plotting above magnitude
-def magnitude_x(srate, size):
+def magnitude_x(srate, size, broken=True):
+    if broken:
+        lament("magnitude_x(broken=True): DEPRECATED")
         return np.arange(0, srate/2, srate/2/size)
+    else:
+        return np.arange(1, size // 2 + 1) / (size // 2) * (srate / 2)
 
 
 def degrees_clamped(x):