From 2cec38a1c14cee1b95cbebc762a355b8193e752a Mon Sep 17 00:00:00 2001 From: Connor Olding Date: Wed, 2 Jan 2019 06:45:12 -0800 Subject: [PATCH] update 49 --- lib/__init__.py | 1 + lib/bs.py | 1 + lib/cepstrum.py | 2 +- lib/fft.py | 36 ++++++++++++++++++++++++------------ lib/plotwav.py | 28 ++++++++++++---------------- lib/smoothfft.py | 40 +++++++++++++++++++++++++++++++++++++--- lib/util.py | 16 ++++++++++++---- 7 files changed, 88 insertions(+), 36 deletions(-) diff --git a/lib/__init__.py b/lib/__init__.py index 5b774f5..33c6a69 100644 --- a/lib/__init__.py +++ b/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 * diff --git a/lib/bs.py b/lib/bs.py index 1a168ec..3a259fb 100644 --- a/lib/bs.py +++ b/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 diff --git a/lib/cepstrum.py b/lib/cepstrum.py index 69b9966..a423dd3 100644 --- a/lib/cepstrum.py +++ b/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] diff --git a/lib/fft.py b/lib/fft.py index 52d8cd8..3ae7f86 100644 --- a/lib/fft.py +++ b/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 - power = np.abs(np.fft.rfft(windowed, 2 * size))**2 - # this scraps the nyquist value to get exactly 'size' outputs - yield power[0:size] - count += 1 - - # assert(segs == count) # this is probably no good in a generator + if broken: + power = np.abs(np.fft.rfft(windowed, 2 * size))**2 + # this scraps the nyquist value to get exactly 'size' outputs + yield power[:size] + else: + power = np.abs(np.fft.rfft(windowed, size))**2 + # 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) - avg = np.zeros(size) - for power in magnitudes(s, size): - avg += power/segs + 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 diff --git a/lib/plotwav.py b/lib/plotwav.py index 48e8e5f..7312d1a 100644 --- a/lib/plotwav.py +++ b/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 = smoothfft(xs_raw, ys_raw, bw=bw) + xs, ys = smoothfft4(ys_raw, 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_s = smoother(xs_raw, ys_raw_s, bw=bw) + xs, ys_m = smoothfft4(ys_raw_m, bw=bw, srate=srate) + 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, - smoother=smoothfft2, **kwargs): +def plotwav2(fn, bw=1, size=8192, fix=False, **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 - with ewave.open(fno, 'w', sampling_rate=srate, - nchannels=count_channels(sf)) as f: - f.write(sf) + wav_write(fno, sf, srate, dtype='f') 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 diff --git a/lib/smoothfft.py b/lib/smoothfft.py index 95b9bdd..11d02ca 100644 --- a/lib/smoothfft.py +++ b/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 diff --git a/lib/util.py b/lib/util.py index 6d5d49d..a0b2acc 100644 --- a/lib/util.py +++ b/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): - return 10*np.log10(np.abs(np.fft.rfft(src, 2 * size))**2)[0:size] +def magnitude(src, size, broken=True): + 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): - return np.arange(0, srate/2, srate/2/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):