From d03bc873687de1886e453143299749af705a57ac Mon Sep 17 00:00:00 2001 From: Connor Olding Date: Mon, 19 Oct 2015 05:39:37 -0700 Subject: [PATCH] update 5 --- autoupdate | 2 +- lib/__init__.py | 13 +++-- lib/plot.py | 33 ++++++++++++- lib/plotwav.py | 126 ++++++++++++++++++++++++++++++++++++++++++++++++ lib/util.py | 2 + 5 files changed, 170 insertions(+), 6 deletions(-) create mode 100644 lib/plotwav.py diff --git a/autoupdate b/autoupdate index 7a55142..21ba36f 100755 --- a/autoupdate +++ b/autoupdate @@ -10,7 +10,7 @@ if [ ! -d .git ]; then [ -e version ] && rm version git add .gitignore '*' '**/*' else - git add -u + git add -A changes="$(git status --porcelain | wc -l)" [ "$changes" -eq 0 ] && exit 0 fi diff --git a/lib/__init__.py b/lib/__init__.py index 89bfe43..164b3f4 100644 --- a/lib/__init__.py +++ b/lib/__init__.py @@ -60,13 +60,16 @@ def c_render(cascade, precision=4096): ys += f(xs) return xs, ys -def c_render2(xs, cascade): +def c_render2(xs, cascade, phase=False): """c_render optimized and specifically for first/second-order filters""" import numexpr as ne j = np.complex(0, 1) - eq2 = 'abs((b0 + j*b1*ws - b2*ws**2)/(a0 + j*a1*ws - a2*ws**2))**2' - eq1 = 'abs((b0 + j*b1*ws)/(a0 + j*a1*ws))**2' - fmt = 'real(log10({})*10 + gain)' + eq2 = '(b0 + j*b1*ws - b2*ws**2)/(a0 + j*a1*ws - a2*ws**2)' + eq1 = '(b0 + j*b1*ws)/(a0 + j*a1*ws)' + if not phase: + fmt = 'real(log10(abs({})**2)*10 + gain)' + else: + fmt = 'arctan2(imag({0}), real({0}))' # gross ys = np.zeros(len(xs)) for f in cascade: w0, ba, gain = f @@ -83,6 +86,8 @@ def c_render2(xs, cascade): raise Exception("incompatible cascade; consider using c_render instead") ws = xs/w0 ys += ne.evaluate(eq) + if phase: + ys = degrees_clamped(ys) return ys def firize(xs, ys, n=4096, srate=44100, plot=None): diff --git a/lib/plot.py b/lib/plot.py index 23efa67..c710415 100644 --- a/lib/plot.py +++ b/lib/plot.py @@ -1,6 +1,8 @@ import matplotlib.pyplot as plt from matplotlib import ticker +# TODO: remove set_size_inches calls, move them inline as necessary + def response_setup(ax, ymin=-24, ymax=24, yL=ticker.AutoMinorLocator(3)): ax.set_xlim(20, 20000) ax.set_ylim(ymin, ymax) @@ -10,6 +12,15 @@ 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) + ax.set_yticks(tuple(range(-180, 180 + 1, int(360/div)))) + ax.yaxis.set_minor_locator(yL) + ax.grid(True, 'both') + ax.set_xlabel('frequency (Hz)') + ax.set_ylabel('phase (degrees)') + def cleanplot(): fig, ax = plt.subplots() fig.set_size_inches((16, 16)) @@ -18,9 +29,29 @@ def cleanplot(): return fig, ax def new_response(*args, **kwargs): - #fig, ax = plt.subplots() fig = plt.figure() ax = fig.gca() response_setup(ax, *args, **kwargs) fig.set_size_inches(10, 6) return fig, ax + +def new_phase_response(*args, **kwargs): + fig = plt.figure() + ax = fig.gca() + phase_response_setup(ax, *args, **kwargs) + fig.set_size_inches(10, 6) + return fig, ax + +def new_bode(magnitude_offset=0): + fig = plt.figure() + ax1 = fig.gca() + ax2 = ax1.twinx() + ymin = -24 + magnitude_offset + ymax = 24 + magnitude_offset + response_setup(ax1, ymin, ymax) + phase_response_setup(ax2) + # ax1 and ax2 should have identical grids; + # disable ax2's so it doesn't overlap ax1's plot lines. + ax2.grid(False, which='both') + fig.set_size_inches(10, 6) + return fig, ax1, ax2 diff --git a/lib/plotwav.py b/lib/plotwav.py new file mode 100644 index 0000000..90df7e2 --- /dev/null +++ b/lib/plotwav.py @@ -0,0 +1,126 @@ +# 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 new_response, show, convolve_each, monoize, count_channels + +import numpy as np + +def plotwavsmooth(fn, ax, tilt=None, bw=1, size=8192, raw=False, fix=False, smoother=smoothfft2, **kwargs): + s, srate = wav_read(fn) + + s, rms = normalize(s, srate) + sm = monoize(s) + ss = monoize(s*np.array((1, -1))) + + xs_raw = np.arange(0, srate/2, srate/2/size) + ys_raw = averfft(sm, size=size) + + # tilting beforehand is negligible besides lowest frequencies, but eh + if tilt is not None: + ys_raw -= tilter2(xs_raw, tilt) + + xs, ys = smoother(xs_raw, ys_raw, bw=bw) + + if not 'label' in kwargs: + kwargs['label'] = fn + + if raw: + ax.semilogx(xs_raw, ys_raw, **kwargs) + ax.semilogx(xs, ys, **kwargs) + + if not fix: return + + fno = fn[:-4]+"-proc.wav" + fir = firize(xs, -ys, srate=srate) + sf = convolve_each(s/8, fir, mode='same') + + import ewave + with ewave.open(fno, 'w', sampling_rate=srate, nchannels=count_channels(sf)) as f: + f.write(sf) + print('wrote '+fno) + +def plotfftsmooth(s, srate, ax, bw=1, tilt=None, size=8192, window=0, raw=False, **kwargs): + sm = monoize(s) + + xs_raw = np.arange(0, srate/2, srate/2/size) + ys_raw = averfft(sm, size=size, mode=window) + + ys_raw -= tilter2(xs_raw, tilt) + + xs, ys = smoothfft(xs_raw, ys_raw, bw=bw) + + if raw: ax.semilogx(xs_raw, ys_raw, **kwargs) + ax.semilogx(xs, ys, **kwargs) + + return xs, ys + +def plotwav2(fn, ax, bw=1, size=8192, raw=False, fix=False, + smoother=smoothfft2, side_compensate=9, **kwargs): + s, srate = wav_read(fn) + + s, rms = normalize(s, srate) + sm = monoize(s) + ss = monoize(s*np.array((1, -1))) + + xs_raw = np.arange(0, srate/2, srate/2/size) + ys_raw = averfft(sm, size=size) + ys_raw_side = averfft(ss, size=size) + + # tilting beforehand is negligible besides lowest frequencies, but eh + ys_raw -= tilter2(xs_raw, 'np2') + ys_raw_side -= tilter2(xs_raw, 'np2s') + + xs, ys = smoother(xs_raw, ys_raw, bw=bw) + xs, ys_side = smoother(xs_raw, ys_raw_side, bw=bw) + + if not 'label' in kwargs: + kwargs['label'] = fn + + if raw: + ax.semilogx(xs_raw, ys_raw, **kwargs) + ax.semilogx(xs_raw, ys_raw_side + side_compensate, **kwargs) + ax.semilogx(xs, ys, **kwargs) + ax.semilogx(xs, ys_side + side_compensate, **kwargs) + + side_gain = np.average(ys_raw_side) - np.average(ys_raw) + print("side gain:", side_gain) + + if not fix: return + fno = fn[:-4]+"-proc.wav" + + fir = firize(xs, -ys, srate=srate) + smf = convolve_each(sm/8, fir, mode='same') + fir = firize(xs, -ys_side, srate=srate) + ssf = convolve_each(ss/8, fir, mode='same') + ssf *= 10**(side_gain/20) + sf = np.array((smf + ssf, smf - ssf)).T + + import ewave + with ewave.open(fno, 'w', sampling_rate=srate, nchannels=count_channels(sf)) as f: + f.write(sf) + print('wrote '+fno) + +def pw(fn, ax, **kwargs): + plotwavsmooth(fn, ax, tilt='np2', bw=1/6, **kwargs) + +def pwc(fn, **kwargs): + fig, ax = new_response(-18, 18) + ax.set_title('averaged magnitudes of normalized songs with tilt and smoothing') + + pw(fn, ax, fix=True, **kwargs) + fno = fn[:-4]+"-proc.wav" + pw(fno, ax, fix=False, **kwargs) + + ax.legend(loc=8) + show(fig) + +def pw2(fn, **kwargs): + fig, ax = new_response(-18, 18) + ax.set_title('averaged magnitudes of normalized songs with tilt and smoothing') + + plotwav2(fn, ax, fix=True, bw=1/6, **kwargs) + fno = fn[:-4]+"-proc.wav" + plotwav2(fno, ax, fix=False, bw=1/6, **kwargs) + + ax.legend(loc=8) + show(fig) diff --git a/lib/util.py b/lib/util.py index 6b40460..34d40a5 100644 --- a/lib/util.py +++ b/lib/util.py @@ -17,6 +17,8 @@ warp = lambda w: np.arctan(w)*2 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] +degrees_clamped = lambda x: ((x*180/np.pi + 180) % 360) - 180 + def xsp(precision=4096): """create #precision log-spaced points from 20 to 20480 Hz""" # i opt not to use steps or linspace here,