dsp/lib/smoothfft.py

103 lines
3.2 KiB
Python
Raw Permalink Normal View History

2019-01-02 06:45:12 -08:00
from . import xsp, lament, ceil2
2015-10-18 23:06:39 -07:00
import numpy as np
2017-09-21 04:04:22 -07:00
2015-10-18 23:06:39 -07:00
def smoothfft(xs, ys, bw=1, precision=512):
"""performs log-lin smoothing on magnitude data,
generally from the output of averfft."""
2019-01-02 06:45:12 -08:00
lament("smoothfft(): DEPRECATED; use smoothfft4 instead.")
2015-10-18 23:06:39 -07:00
xs2 = xsp(precision)
ys2 = np.zeros(precision)
log_xs = np.log(xs)
for i, x in enumerate(xs2):
dist = np.exp(np.abs(log_xs - np.log(x + 1e-35)))
window = np.maximum(0, 1 - (dist - bw))
2015-10-18 23:33:46 -07:00
# at this point we could probably
2015-10-18 23:06:39 -07:00
# normalize our *triangular* window to 0-1
# and transform it into *another* windowing function
wsum = np.sum(window)
ys2[i] = np.sum(ys*window/wsum)
return xs2, ys2
2017-09-21 04:04:22 -07:00
2015-10-18 23:06:39 -07:00
def smoothfft2(xs, ys, bw=1, precision=512, compensate=True):
"""performs log-lin smoothing on magnitude data,
generally from the output of averfft."""
2019-01-02 06:45:12 -08:00
lament('smoothfft2: DEPRECATED; use smoothfft4 instead.')
2015-10-18 23:06:39 -07:00
xs2 = xsp(precision)
ys2 = np.zeros(precision)
log2_xs2 = np.log2(xs2)
for i, x in enumerate(xs):
2015-10-18 23:33:46 -07:00
# before optimizations: dist = np.abs(np.log2(xs2/(x + 1e-35)))/bw
2015-10-18 23:06:39 -07:00
dist = np.abs(log2_xs2 - np.log2(x + 1e-35))/bw
2018-02-20 04:04:25 -08:00
# window = np.maximum(0, 1 - dist) # triangular
window = np.exp(-dist**2/(0.5/2)) # gaussian (untruncated)
2015-10-18 23:06:39 -07:00
ys2 += ys[i]*window
if compensate:
2017-09-21 04:04:22 -07:00
_, temp = smoothfft2(xs, np.ones(len(xs)),
bw=bw, precision=precision, compensate=False)
2015-10-18 23:06:39 -07:00
ys2 /= temp
return xs2, ys2
2018-02-19 04:04:28 -08:00
2018-03-03 04:04:22 -08:00
def smoothfft_setup(size, precision=512, bw=1/6):
2019-01-02 06:45:12 -08:00
lament('smoothfft_setup(): DEPRECATED; use smoothfft_setup2 instead.')
2018-03-03 04:04:22 -08:00
dotme = np.zeros((size, precision))
2018-02-19 04:04:28 -08:00
2018-03-03 04:04:22 -08:00
xs = np.arange(0, 1, 1/size)
2018-02-20 04:04:25 -08:00
xs2 = np.logspace(-np.log2(size), 0, precision, base=2)
2018-02-19 04:04:28 -08:00
comp = np.zeros(precision)
log2_xs2 = np.log2(xs2)
for i, x in enumerate(xs):
dist = np.abs(log2_xs2 - np.log2(x + 1e-35)) / bw
2018-02-20 04:04:25 -08:00
window = np.exp(-dist**2 * 4) # gaussian (untruncated)
2018-02-19 04:04:28 -08:00
comp += window
2018-03-03 04:04:22 -08:00
dotme[i] = window
2018-02-19 04:04:28 -08:00
2018-03-03 04:04:22 -08:00
dotme /= comp
return xs2, dotme
def smoothfft3(ys, bw=1, precision=512, srate=None):
"""performs log-lin smoothing on magnitude data"""
2019-01-02 06:45:12 -08:00
lament('smoothfft3(): DEPRECATED; use smoothfft4 instead.')
2018-03-03 04:04:22 -08:00
xs2, dotme = smoothfft_setup(len(ys), precision, bw)
if srate is None:
return xs2, ys @ dotme
else:
return xs2 * (srate / 2), ys @ dotme
2019-01-02 06:45:12 -08:00
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