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).

reson.py

#!/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()

variable_butterworth.py

#!/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()