Applying time-variant filter in Python

Question

I'm attempting to apply a bandpass filter with time-varying cutoff frequencies to a signal, using Python. The routine I am currently using partitions my signal into equal-length time segments, then for each segment I apply a filter with time-specific parameters, before merging the signal back together. The parameters are based on pre-existing estimates.

The problem I seem to be having is that there are "ripples" at the edge of each time-segment that appear after the filter has been applied. This causes discontinuities in my signal, which interfere with my post-filtering data analysis.

I was hoping someone could inform me whether there are any existing routines for implementing filters with time-varying parameters in Python? Alternatively, advice on how I might get around this problem would be much appreciated.

EDIT -- example of what I want to do is added below.

Let's say I have a signal x(t). I want to filter the first half with a bandpass filter with parameters (100,200) Hz. The second half I want to filter with parameters (140, 240) Hz. I iterate over x(t), applying my filter to each half, then recombine the results. Some example code might look like:

outputArray = np.empty(len(x))
segmentSize = len(x) / 2
filtParams  = [(100, 200), (140, 240)]

for i in range(2):
    tempData     = x[i*segmentSize:(i+1)*segmentSize]
    tempFiltered = bandPassFilter(tempData, parameters=filtParams[i])
    outputArray[i*segmentSize:(i+1)*segmentSize] = tempFiltered

(To save space let's assume I have a function which performs bandpass filtering).

As you can see, the data segments do not overlap and are simply "pasted" back together in the new array.

EDIT 2 -- some sample code of my problem @H.D.

First of all, thanks for your significant input thus far. The audiolazy package looks like a great tool.

I thought it would be a bit more useful if I describe my goals in further detail. As I have posted elsewhere, I am attempting to extract the instantaneous frequency (IF) of a signal, using the Hilbert transform. My data contains significant noise but I have a good estimate of the bandwidth where my IF signal lies. A problem I have come up against, however, is that the IF is often nonstationary. Using a "static" filter approach I am often therefore required to use a broad bandpass region, to ensure all frequencies are captured.

The following code demonstrates the effect of increasing the filter bandwidth on an IF signal. It includes a signal generating function, an implementation of a bandpass filter using the scipy.signal package, and a method to extract the IF of the resultant filtered signal.

from audiolazy import *
import scipy.signal as sig
import numpy as np
from pylab import *

def sineGenerator( ts, f, rate, noiseLevel=None ):
    """generate a sine tone with time, frequency, sample rate and noise specified"""

    fs = np.ones(len(ts)) * f    
    y  = np.sin(2*np.pi*fs*ts)

    if noiseLevel: y = y + np.random.randn(len(y))/float(noiseLevel)
    return y

def bandPassFilter( y, passFreqs, rate, order ):
    """STATIC bandpass filter using scipy.signal Butterworth filter"""

    nyquist = rate / 2.0
    Wn      = np.array([passFreqs[0]/nyquist, passFreqs[1]/nyquist])    
    z, p, k = sig.butter(order, Wn, btype='bandpass', output='zpk')
    b, a    = sig.zpk2tf(z, p, k)

    return sig.lfilter(b, a, y)

if __name__ == '__main__':
     rate = 1e4
     ts   = np.arange(0, 10, 1/rate)

     # CHANGING THE FILTER AFFECTS THE LEVEL OF NOISE
     ys    = sineGenerator(ts, 600.0, 1e4, noiseLevel=1.0) # a 600Hz signal with noise
     filts = [[500, 700], [550, 650], [580, 620]]

    for f in filts:
        tempFilt = bandPassFilter( ys, f, rate, order=2 )
        tempFreq = instantaneousFrequency( tempFilt, rate )

        plot( ts[1:], tempFreq, alpha=.7, label=str(f).strip('[]') )

    ylim( 500, 750 )
    xlabel( 'time' )
    ylabel( 'instantaneous frequency (Hz)' )

    legend(frameon=False)
    title('changing filter passband and instantaneous frequency')
    savefig('changingPassBand.png')

There is a single frequency component in the signal (at 600Hz). The legend shows the filter parameters used in each case. Using a narrower "static" filter gives a "cleaner" output. But how narrow my filter can be is limited by what the frequencies are. For instance, consider the following signal with two frequency components (one at 600Hz, another at 650Hz).

In this example I have been forced to use a broader bandpass filter, which has resulted in extra noise creeping in to the IF data.

My idea is that by using a time varying filter, I can "optimise" the filter for my signal at certain time increments. So for the above signal I might want to filter around 580-620Hz for the first 5 seconds, then 630-670Hz for the next 5 seconds. Essentially I want to minimise noise in the final IF signal.

Based on the example code you sent I have written a function that uses audiolazy to implement a static Butterworth filter on a signal.

def audioLazyFilter( y, rate, Wp, Ws ):
    """implement a Butterworth filter using audiolazy"""

    s, Hz = sHz(rate)
    Wp    = Wp * Hz # Bandpass range in rad/sample
    Ws    = Ws * Hz # Bandstop range in rad/sample

    order, new_wp_divpi = sig.buttord(Wp/np.pi, Ws/np.pi, gpass=dB10(.6), gstop=dB10(.1))
    ssfilt = sig.butter(order, new_wp_divpi, btype='bandpass')
    filt_butter = ZFilter(ssfilt[0].tolist(), ssfilt[1].tolist())

    return list(filt_butter(y))

The IF data obtained using this filter in conjunction with the Hilbert transform routine compares well to those obtained using scipy.signal:

AL_filtered = audioLazyFilter( ys, rate, np.array([580, 620]), np.array([570, 630]) )
SP_filtered = bandPassFilter( ys, [580, 620], rate, order=2 )
plot(ts[1:], instantaneousFrequency( SP_filtered, 1e4 ), alpha=.75, label='scipy.signal Butterworth filter')
plot(ts[1:], instantaneousFrequency( AL_filtered, 1e4 ), 'r', alpha=.75, label='audiolazy Butterworth filter')

My question is now can I combine the audiolazy Butterworth routine with the time-varying properties you described in your original posts?

Answer 1

AudioLazy works natively with time varying filters

from audiolazy import sHz, white_noise, line, resonator, AudioIO

rate = 44100
s, Hz = sHz(rate)

sig = white_noise() # Endless white noise Stream

dur = 8 * s # Some few seconds of audio
freq = line(dur, 200, 800) # A lazy iterable range
bw = line(dur, 100, 240)

filt = resonator(freq * Hz, bw * Hz) # A simple bandpass filter

with AudioIO(True) as player:
  player.play(filt(sig), rate=rate)

You can also use it for plotting (or analysis, in general), by using list(filt(sig)) or filt(sig).take(inf). There are a lot of other resources that might be useful as well, such as applying time-varying coefficients directly in a Z-transform filter equation.

EDIT: More information about the AudioLazy components

The following examples were done using the IPython.

Resonator is a StrategyDict instance, which ties many implementations in one place.

In [1]: from audiolazy import *

In [2]: resonator
Out[2]: 
{('freq_poles_exp',): <function audiolazy.lazy_filters.freq_poles_exp>,
 ('freq_z_exp',): <function audiolazy.lazy_filters.freq_z_exp>,
 ('poles_exp',): <function audiolazy.lazy_filters.poles_exp>,
 ('z_exp',): <function audiolazy.lazy_filters.z_exp>}

In [3]: resonator.default
Out[3]: <function audiolazy.lazy_filters.poles_exp>

So resonator calls internally the resonator.poles_exp function, from which you can get some help

In [4]: resonator.poles_exp?
Type:       function
String Form:<function poles_exp at 0x2a55b18>
File:       /usr/lib/python2.7/site-packages/audiolazy/lazy_filters.py
Definition: resonator.poles_exp(freq, bandwidth)
Docstring:
Resonator filter with 2-poles (conjugated pair) and no zeros (constant
numerator), with exponential approximation for bandwidth calculation.

Parameters
----------
freq :
  Resonant frequency in rad/sample (max gain).
bandwidth :
  Bandwidth frequency range in rad/sample following the equation:

    ``R = exp(-bandwidth / 2)``

  where R is the pole amplitude (radius).

Returns
-------
A ZFilter object.
Gain is normalized to have peak with 0 dB (1.0 amplitude).

So a verbose filter assignment would be

filt = resonator.poles_exp(freq=freq * Hz, bandwidth=bw * Hz)

Where the Hz is just a number to change the unit from Hz to rad/sample, as used in most AudioLazy components.

Let's do so with freq = pi/4 and bw = pi/8 (pi is already in the audiolazy namespace):

In [5]: filt = resonator(freq=pi/4, bandwidth=pi/8)

In [6]: filt
Out[6]: 
              0.233921
------------------------------------
1 - 1.14005 * z^-1 + 0.675232 * z^-2


In [7]: type(filt)
Out[7]: audiolazy.lazy_filters.ZFilter

You can try using this filter instead of the one given in the first example.

Another way to do so would be using the z object from the package. First let's find the constants for that all-poles resonator:

In [8]: freq, bw = pi/4, pi/8

In [9]: R = e ** (-bw / 2)

In [10]: c = cos(freq) * 2 * R / (1 + R ** 2) # AudioLazy included the cosine

In [11]: gain = (1 - R ** 2) * sqrt(1 - c ** 2)

The denominator can be done directly by using the z in the equation:

In [12]: denominator = 1 - 2 * R * c * z ** -1 + R ** 2 * z ** -2

In [13]: gain / denominator
Out[14]: 
              0.233921
------------------------------------
1 - 1.14005 * z^-1 + 0.675232 * z^-2

In [15]: type(_) # The "_" is the last returned value in IPython
Out[15]: audiolazy.lazy_filters.ZFilter

EDIT 2: About the time varying coefficients

The filter coefficients can also be a Stream instance (which can be cast from any iterable).

In [16]: coeff = Stream([1, -1, 1, -1, 1, -1, 1, -1, 1, -1]) # Cast from a list

In [17]: (1 - coeff * z ** -2)(impulse()).take(inf)
Out[17]: [1.0, 0.0, -1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

The same, given a list input instead of the impulse() Stream:

In [18]: coeff = Stream((1, -1, 1, -1, 1, -1, 1, -1, 1, -1)) # Cast from a tuple 

In [19]: (1 - coeff * z ** -2)([1, 0, 0, 0, 0, 0, 0]).take(inf)
Out[19]: [1.0, 0.0, -1, 0, 0, 0, 0]

A NumPy 1D array is also an iterable:

In [20]: from numpy import array

In [21]: array_data = array([1, -1, 1, -1, 1, -1, 1, -1, 1, -1])

In [22]: coeff = Stream(array_data) # Cast from an array

In [23]: (1 - coeff * z ** -2)([0, 1, 0, 0, 0, 0, 0]).take(inf)
Out[23]: [0.0, 1.0, 0, 1, 0, 0, 0]

This last example shows the time-variant behaviour.

EDIT 3: Chunked-repeat sequences behaviour

The line function has a behaviour similar to the numpy.linspace, which gets the range "length" instead of "step".

In [24]: import numpy

In [25]: numpy.linspace(10, 20, 5) # Start, stop (included), length
Out[25]: array([ 10. ,  12.5,  15. ,  17.5,  20. ])

In [26]: numpy.linspace(10, 20, 5, endpoint=False) # Makes stop not included
Out[26]: array([ 10.,  12.,  14.,  16.,  18.])

In [27]: line(5, 10, 20).take(inf) # Length, start, stop (range-like)
Out[27]: [10.0, 12.0, 14.0, 16.0, 18.0]

In [28]: line(5, 10, 20, finish=True).take(inf) # Include the "stop"
Out[28]: [10.0, 12.5, 15.0, 17.5, 20.0]

With that, the filter equation has a different behaviour sample-per-sample (1-sample "chunk"). Anyhow, you can use a repeater for larger chunk sizes:

In [29]: five_items = _ # List from the last Out[] value

In [30]: @tostream
   ....: def repeater(sig, n):
   ....:     for el in sig:
   ....:         for _ in xrange(n):
   ....:             yield el
   ....:             

In [31]: repeater(five_items, 2).take(inf)
Out[31]: [10.0, 10.0, 12.5, 12.5, 15.0, 15.0, 17.5, 17.5, 20.0, 20.0]

And use it in the line from the first example, so that freq and bw becomes:

chunk_size = 100
freq = repeater(line(dur / chunk_size, 200, 800), chunk_size)
bw = repeater(line(dur / chunk_size, 100, 240), chunk_size)

EDIT 4: Emulating time-varying filters/coefficients from LTI filters using time-varying gain/envelope

Another way around would be using different "weights" for two different filtered versions of the signal, and making some "crossfade" math with the signal, something like:

signal = thub(sig, 2) # T-Hub is a T (tee) auto-copy
filt1(signal) * line(dur, 0, 1) + filt2(signal) * line(dur, 1, 0)

This would apply a linear envelope (from 0 to 1 and from 1 to 0) from different filtered versions of the same signal. If thub looks confusing, try sig1, sig2 = tee(sig, 2) applying filt(sig1) and filt(sig2) instead, these should do the same.

EDIT 5: Time-variant Butterworth filter

I spent the last hours trying to let that Butterworth be personalized as your example, imposing order = 2 and giving the half-power bandwidth (~3dB) directly. I've done four examples, the code is in this Gist, and I've updated AudioLazy to include a gauss_noise Gaussian-distributed noise stream. Please note that the code in gist has nothing optimized, it was done ony to work in this particular case, and the chirp example makes it really slow due to a "per sample" coefficient finding behaviour. The instant frequency can be get from the [filtered] data in rad/sample with:

diff(unwrap(phase(hilbert(filtered_data))))

where diff = 1 - z ** -1 or another approach to find derivatives in discrete time, hilbert is the function from scipy.signal that gives us the analytical signal (the Discrete Hilbert Transform is the imaginary part of its result) and the other two are helper functions from AudioLazy.

This is what happens when Butterworth changes its coefficients abruptly while keeping its memory, without noise:

It's noticeable a oscilatory behaviour in this transition. You can use a moving median to "smooth" that in the lower frequency side while keeping the abrupt transition, but that won't work with the higher frequency. Well, that was what we would expect from a perfect sinusoid, but with noise (a LOT of noise, the Gaussian has the standard deviation equals to the sinusoid amplitude), it becomes:

I tried then to do the same with a chirp, precisely this:

This shows a strange behaviour when filtering with the lower bandwidth, at the top frequency. And with the additive noise:

The code in gist also AudioIO().play this last noisy chirp.

EDIT 6: Time-variant resonator filter

I've added to the same Gist an example using resonators instead of Butterworth. They're in pure Python and aren't optimized, but performs faster than calling butter for each sample during a chirp, and is far easier to implement, as all the resonator strategies accepts Stream instances as valid inputs. Here are the plots for a cascade of two resonators (i.e., a 2nd order filter):

And the same for a cascade of three resonators (i.e., a 3nd order filter):

These resonators have gain equals to 1 (0 dB) at the center frequency, and that oscillation pattern from the "Abrupt pure sinusoid" plots in the transition happens even without any filtering at all.

Answer 2

Try filtering the whole signal with each filter you're using, then merge the filtered signals appropriately. Roughly, in pseudocode:

# each filter acts on the entire signal with each cutoff frequency
s1 = filter1(signal)
s2 = filter2(signal)
s3 = filter3(signal)

filtered_signal = s1[0:t1] + s2[t1:t2] + s3[t2:t3]

I think this will avoid the artifacts you're describing that result from chopping up the signal then filtering.

Another possibility is using a Short Time Fourier Transform (STFT). Here's an implementation using numpy.

Basically, you could STFT your signal, filter your signal by operating on the time-frequency array, then inverse STFT the array to get a filtered signal.

You might get even better results using an invertible wavelet transform. Here's a paywalled paper describing how to do pretty much what you're trying to do with a wavelet transform.