Last active
March 18, 2017 18:20
-
-
Save danilobellini/6419374 to your computer and use it in GitHub Desktop.
Trying to do a time-variant Butterworth filter using SciPy and AudioLazy. Also, a faster time-variant resonator (AudioLazy) for comparison. Both examples includes instantaneous frequency measurements (from Discrete Hilbert Transform).
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
# -*- coding: utf-8 -*- | |
# Created on Mon Sep 09 23:06:21 2013 | |
# @author: Danilo J. S. Bellini | |
""" | |
AudioLazy time-variant resonator helping instantaneous frequency measurements. | |
""" | |
from __future__ import print_function | |
from audiolazy import (sHz, sinusoid, gauss_noise, line, z, inf, resonator, | |
phase, pi, unwrap, Stream, repeat, thub, CascadeFilter) | |
from scipy.signal import hilbert | |
import pylab | |
# Basic configuration | |
rate = 1e4 | |
s, Hz = sHz(rate) | |
# Data/signal model | |
dur = 2 * s | |
time_axis = line(dur, 0, dur/s).take(inf) # For plots | |
freq1 = Stream(600 * Hz).limit(dur / 2).append(650 * Hz).take(dur) | |
freq2 = Stream( | |
repeat(200 * Hz).limit(dur / 2), | |
line(.1 * s, 200 * Hz, 900 * Hz), | |
repeat(900 * Hz), | |
).take(dur) | |
signals = { | |
"Abrupt noisy": sinusoid(freq1) + gauss_noise(), # Endless | |
"Abrupt pure sinusoid": sinusoid(freq1), | |
"Noisy": sinusoid(freq2) + gauss_noise(), | |
"Pure sinusoid": sinusoid(freq2), | |
} | |
limits = { | |
"Abrupt noisy": dict(xmin=.85, xmax=1.15, ymin=520, ymax=720), | |
"Abrupt pure sinusoid": dict(xmin=.95, xmax=1.05, ymin=590, ymax=680), | |
"Noisy": dict(xmin=.87, xmax=1.17, ymin=110, ymax=1080), | |
"Pure sinusoid": dict(xmin=.97, xmax=1.17, ymin=130, ymax=1080), | |
} | |
# Filter parameters | |
bandwidths = [200 * Hz, 100 * Hz, 40 * Hz] | |
order = 2 # Number of pole pairs (i.e., filter would have 2 * order roots) | |
# Helpers | |
diff = 1 - z ** -1 # Simplest difference filter (discrete derivative) | |
for name, sig in signals.items(): | |
data = sig.take(dur) # Don't try "inf" here! | |
freq = freq1 if name.startswith("Abrupt") else freq2 | |
# Starts plotting | |
title = "{} - Instantaneous frequency - Order {}".format(name, order) | |
pylab.figure(title) | |
for bw in bandwidths: | |
# Creates and apply the time-variant filters | |
f = thub(freq, order) | |
filt = CascadeFilter(resonator.z_exp(f, bw * order ** .5) | |
for idx in range(order)) | |
filtered_data = filt(data).take(inf) | |
# Finds the instant frequency | |
analytical_signal_phase = unwrap(phase(hilbert(filtered_data))) | |
inst_freq_stream = diff(analytical_signal_phase) | |
inst_freq = (inst_freq_stream / Hz).take(inf) # in Hz, as a list | |
# Some debug information | |
inst_freq_np = pylab.array(inst_freq) | |
print("=== {} - Bandwidth {} Hz ===".format(name, bw / Hz)) | |
print("Frequency - first half (in Hz):", inst_freq_np[:int(dur/2)].mean()) | |
print("... standard deviation (in Hz):", inst_freq_np[:int(dur/2)].std()) | |
print() | |
# Plotting | |
pylab.plot(time_axis, inst_freq, label="Bandwidth {} Hz".format(bw / Hz)) | |
# Finishes the plot | |
pylab.title(title) | |
pylab.legend(loc="best") | |
pylab.axis(**limits[name]) | |
fname = "reson_{}_{}.png".format(order, name.lower().replace(" ", "_")) | |
pylab.savefig(fname) | |
# Frequency response for the filters we're using | |
for f in [200 * Hz, 600 * Hz, 900 * Hz]: | |
for bw in bandwidths: | |
filt = CascadeFilter(resonator.z_exp(f, bw * order ** .5) | |
for idx in range(order)) | |
filt.plot(pylab.figure("f = {} Hz; bw = {} Hz".format(f/Hz, bw/Hz)), | |
rate = rate, | |
min_freq = max(f - bw, 0), | |
max_freq = min(f + bw, pi), | |
) | |
pylab.show() |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#!/usr/bin/env python | |
# -*- coding: utf-8 -*- | |
# Created on Mon Sep 2 18:40:22 2013 | |
# @author: Danilo J. S. Bellini | |
""" | |
Trying to do a time-variant Butterworth filter using SciPy and AudioLazy. | |
""" | |
from __future__ import print_function | |
from audiolazy import (sHz, sinusoid, gauss_noise, line, z, inf, ZFilter, | |
phase, pi, unwrap, Stream, repeat, AudioIO) | |
from scipy.signal import butter, hilbert | |
import pylab | |
# Basic configuration | |
rate = 1e4 | |
s, Hz = sHz(rate) | |
# Data/signal model | |
dur = 2 * s | |
time_axis = line(dur, 0, dur/s).take(inf) # For plots | |
freq1 = Stream(600 * Hz).limit(dur / 2).append(650 * Hz).take(dur) | |
freq2 = Stream( | |
repeat(200 * Hz).limit(dur / 2), | |
line(.1 * s, 200 * Hz, 900 * Hz), | |
repeat(900 * Hz), | |
).take(dur) | |
signals = { | |
"Abrupt noisy": sinusoid(freq1) + gauss_noise(), # Endless | |
"Abrupt pure sinusoid": sinusoid(freq1), | |
"Noisy": sinusoid(freq2) + gauss_noise(), | |
"Pure sinusoid": sinusoid(freq2), | |
} | |
limits = { | |
"Abrupt noisy": dict(xmin=.85, xmax=1.15, ymin=520, ymax=720), | |
"Abrupt pure sinusoid": dict(xmin=.95, xmax=1.05, ymin=590, ymax=680), | |
"Noisy": dict(xmin=.87, xmax=1.17, ymin=110, ymax=1080), | |
"Pure sinusoid": dict(xmin=.97, xmax=1.17, ymin=130, ymax=1080), | |
} | |
# Play the "Noisy" example | |
with AudioIO(True) as player: | |
player.play(signals["Noisy"].copy().append(0).limit(3 * s)) | |
print("Not playing anymore!") | |
print() | |
# Filter parameters | |
bandwidths = [200 * Hz, 100 * Hz, 40 * Hz] | |
order = 2 # Care should be done: this script doesn't design a butterworth from | |
# bandpass/bandstop range and gain restrictions, but uses the given | |
# bandwidth directly, centralized, and with this fixed order | |
# Helpers | |
diff = 1 - z ** -1 # Simplest difference filter (discrete derivative) | |
def bandlimits(f, bw): | |
half_bw = bw * .5 | |
return pylab.array([f - half_bw, f + half_bw]) | |
for name, sig in signals.items(): | |
data = sig.take(dur) # Don't try "inf" here! | |
freq = freq1 if name.startswith("Abrupt") else freq2 | |
used_freqs = set(freq) | |
# Starts plotting | |
title = "{} - Instantaneous frequency".format(name) | |
pylab.figure(title) | |
for bw in bandwidths: | |
# Creates a butterworth for all frequencies (slooow) | |
fdict = { | |
f: map(list, butter(order, bandlimits(f, bw) / pi, btype="bandpass")) | |
for f in used_freqs | |
} | |
# Gets the coefficients from the butterworth filters | |
nums, dens = zip(*[fdict[f] for f in freq]) | |
numerator = map(Stream, zip(*nums)) | |
denominator = map(Stream, zip(*dens)) | |
# Creates and apply the time-variant filters | |
filt_butter = ZFilter(numerator, denominator) | |
filtered_data = filt_butter(data).take(inf) | |
# Finds the instant frequency | |
analytical_signal_phase = unwrap(phase(hilbert(filtered_data))) | |
inst_freq_stream = diff(analytical_signal_phase) | |
inst_freq = (inst_freq_stream / Hz).take(inf) | |
# Some debug information | |
inst_freq_np = pylab.array(inst_freq) | |
print("=== {} - Bandwidth {} Hz ===".format(name, bw / Hz)) | |
print("Frequency - first half (in Hz):", inst_freq_np[:int(dur/2)].mean()) | |
print("... standard deviation (in Hz):", inst_freq_np[:int(dur/2)].std()) | |
print() | |
# Plotting | |
pylab.plot(time_axis, inst_freq, label="Bandwidth {} Hz".format(bw / Hz)) | |
# Finishes the plot | |
pylab.title(title) | |
pylab.legend(loc="best") | |
pylab.axis(**limits[name]) | |
fname = "variable_butterworth_{}.png".format(name.lower().replace(" ", "_")) | |
pylab.savefig(fname) | |
pylab.show() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment