Applying a time filter in Python

Question

I am trying to apply a band-pass filter with time-varying cutoff frequencies to a signal using Python. The routine that I am currently using divides my signal into equal time segments, then for each segment I apply a filter with time-specific parameters before combining the signal back together. Parameters are based on preexisting estimates.

The problem that I seem to encounter is that there are “ripples” on the edge of each time segment that appear after applying the filter. This causes gaps in my signal that interfere with the analysis of the data after filtering.

I was hoping someone could tell me if there are any existing procedures for implementing filters with time-varying parameters in Python? Alternatively, tips on how I can get around this issue would be greatly appreciated.

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

Say I have a signal x (t). I want to filter the first half with a band-pass filter with parameters (100,200) Hz. In the second half I want to filter with the parameters (140, 240) Hz. I iterate over x (t), applying my filter to each half, and then recombining the results. A sample code might look like this:

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, suppose I have a function that filters a band-pass filter).

As you can see, the data segments do not overlap and are simply “inserted” back into the new array.

EDIT 2 is some example code of my @HD problem

First of all, thanks for your significant contribution. The audio package looks like a great tool.

I thought it would be a little more useful if I described my goals in more detail. Since I posted elsewhere , I'm trying 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. However, the problem I came across is that IF is often unsteady. Therefore, using the “static” filter approach, I often need to use a wide bandwidth to capture all frequencies.

The following code demonstrates the effect of increasing the filter passband on the IF signal. It includes a signal generation function, an implementation of a bandpass filter using the scipy.signal package, and an IF extraction method for the resulting 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)) * fy = 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') 

The signal (at a frequency of 600 Hz) has one frequency component. The legend shows the filter parameters used in each case. Using a narrower “static” filter gives a “cleaner” output. But how narrow my filter is, it is limited by what the frequencies are. For example, consider the following signal with two frequency components (one at a frequency of 600 Hz, the other at a frequency of 650 Hz).

In this example, I had to use a wider bandpass filter, which led to the addition of additional noise to the IF data.

My idea is that by using a time-varying filter, I can “optimize” the filter for my signal with a specific time step. So, for the above signal, I can filter around 580-620 Hz for the first 5 seconds, then 630-670 Hz for the next 5 seconds. In essence, I want to minimize the noise in the final IF signal.

Based on the code you provided, I wrote a function that uses audio climbers to implement a Butterworth static filter by 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)) 

IF data obtained using this filter in combination with the Hilbert transform routine compares well with data 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') 

Now my question is: can I combine the Butterworth audio climber procedure with the time-varying properties that you described in your original posts?

Answer 1

AudioLazy works 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 to build (or analyze in general) with list(filt(sig)) or filt(sig).take(inf) . There are many other resources that can be useful, for example, applying time-varying coefficients directly in the Z-transform filter equation.

EDIT: Additional Information on AudioLazy Components

The following examples were performed using IPython.

Resonator is an instance of StrategyDict that links 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 the internal function resonator.poles_exp 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 the purpose of a detailed filter would be

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

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

Do this with freq = pi/4 and bw = pi/8 ( pi 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, not the one specified in the first example.

Another way to do this is to use the z object from the package. First we find the constants for this resonator of all poles:

 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 made directly using 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 time-varying coefficients

Filter coefficients can also be Stream instances (which can be inferred 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] 

Same thing, given entering a list instead of 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] 

The NumPy 1D array is also 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 time behavior.

EDIT 3: repetitive sequence behavior

The line function has a behavior similar to numpy.linspace , which gets a range of "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] 

Moreover, the filter equation has a different sample behavior per sample (1-sample "piece"). Anyway, you can use a repeater for large block 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 freq and bw will become:

 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: Emulation of time-varying filters / coefficients from LTI filters using a time-varying gain / envelope

Another way: use different “weights” for two different filtered versions of the signal and do 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 will 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) apply filt(sig1) and filt(sig2) , they should do the same.

EDIT 5: Temporary Butterworth Filter Option

I spent the last hours trying to let this Butterworth be personalized as your example, imposing order = 2 and providing half-power bandwidth (~ 3dB) directly. I made four examples, the code in this Gist , and I updated AudioLazy to include the gauss_noise Gaussian distributed noise stream. Please note that the code in gist has nothing optimized, it was made ony to work in this particular case, and the chirp example makes it very slow due to the behavior for determining the coefficient on the sample. The current frequency can be obtained from the [filtered] data in rad / sample using:

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

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

Here's what happens when Butterworth dramatically changes its coefficients, preserving its memory without noise:

There is noticeable behavior in this transition. You can use the moving median to “smooth” that at the lower frequency, while maintaining a sharp transition, but this will not work at a higher frequency. Well, this is what we expect from a perfect sine wave, but with noise (MUCH noise, Gaussian has a standard deviation equal to the sinusoidal amplitude), it becomes:

I tried to do the same with the chirp, exactly this:

This shows strange filtering behavior with a lower passband at the upper frequency. And with extra noise:

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

EDIT 6: Filter with a temporary version of the resonator

I added the example using resonators instead of Butterworth to the same Gist . They are in pure Python and not optimized, but they are faster than calling butter for each sample during the chirp, and it is much easier to implement, since all resonator strategies accept Stream instances as valid inputs. Here are the graphs for the cascade of two resonators (i.e. a 2nd order filter):

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

<T411>

These resonators have a gain of 1 (0 dB) at the center frequency, and that the oscillation diagram from the “Sharp Pure Sine Wave” plots occurs during the transition even without any filtering.

Answer 2

Try filtering the entire signal with each filter that you use, then merge the filtered signals accordingly. 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 that this will avoid the artifacts that you describe, which is the result of chopping the signal, then filtering.

Another possibility is to use the Short Four Four Transform (STFT) . Here's an implementation using numpy .

Basically, you can STFT your signal, filter your signal by working on a frequency-frequency array, and then invert the STFT array to get a filtered signal.

You can get even better results using the reversible wavelet transform . Here's a paywalled paper describing how to do pretty much what you are trying to do with a wavelet transform.

Answer 3

If you extract part of your signal using slices, then you effectively close your data with a rectangular window that “wobbles” around the edges due to sudden gaps. One way to fix this is by using a window that rings less, for example, shelter windows:

 import numpy as np signal = np.random.randn(222) length = 50 window = np.hanning(length) for i in range(0, len(signal)-length, length): do_filtering(signal[i:i+length] * window) 

Windows details: http://en.m.wikipedia.org/wiki/Window_function