filtering - Applying time-variant filter in Python -

i'm attempting apply bandpass filter time-varying cutoff frequencies signal, using python. routine using partitions signal equal-length time segments, each segment apply filter time-specific parameters, before merging signal together. parameters based on pre-existing estimates.

the problem seem having there "ripples" @ edge of each time-segment appear after filter has been applied. causes discontinuities in signal, interfere post-filtering data analysis.

i hoping inform me whether there existing routines implementing filters time-varying parameters in python? alternatively, advice on how might around problem appreciated.

edit -- example of want added below.

let's have signal x(t). want filter first half bandpass filter parameters (100,200) hz. second half want filter parameters (140, 240) hz. iterate on x(t), applying filter each half, recombine results. example code might like:

outputarray = np.empty(len(x)) segmentsize = len(x) / 2 filtparams  = [(100, 200), (140, 240)]  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 have function performs bandpass filtering).

as can see, data segments not overlap , "pasted" in new array.

edit 2 -- sample code of problem @h.d.

first of all, significant input far. audiolazy package looks great tool.

i thought bit more useful if describe goals in further detail. have posted elsewhere, attempting extract instantaneous frequency (if) of signal, using hilbert transform. data contains significant noise have estimate of bandwidth if signal lies. problem have come against, however, if nonstationary. using "static" filter approach therefore required use broad bandpass region, ensure frequencies captured.

the following code demonstrates effect of increasing filter bandwidth on if signal. includes signal generating function, implementation of bandpass filter using scipy.signal package, , method extract if of resultant filtered signal.

from audiolazy import * import scipy.signal sig import numpy np pylab import *  def sinegenerator( ts, f, rate, noiselevel=none ):     """generate sine tone time, frequency, sample rate , 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,    = sig.zpk2tf(z, p, k)      return sig.lfilter(b, a, y)  if __name__ == '__main__':      rate = 1e4      ts   = np.arange(0, 10, 1/rate)       # changing filter affects level of noise      ys    = sinegenerator(ts, 600.0, 1e4, noiselevel=1.0) # 600hz signal noise      filts = [[500, 700], [550, 650], [580, 620]]      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 , instantaneous frequency')     savefig('changingpassband.png') 

there single frequency component in signal (at 600hz). legend shows filter parameters used in each case. using narrower "static" filter gives "cleaner" output. how narrow filter can limited frequencies are. instance, consider following signal 2 frequency components (one @ 600hz, @ 650hz).

in example have been forced use broader bandpass filter, has resulted in noise creeping in if data.

my idea using time varying filter, can "optimise" filter signal @ time increments. above signal might want filter around 580-620hz first 5 seconds, 630-670hz next 5 seconds. want minimise noise in final if signal.

based on example code sent have written function uses audiolazy implement static butterworth filter on signal.

def audiolazyfilter( y, rate, wp, ws ):     """implement 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 filter in conjunction hilbert transform routine compares 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 can combine audiolazy butterworth routine time-varying properties described in original posts?

audiolazy works natively 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 # few seconds of audio freq = line(dur, 200, 800) # lazy iterable range bw = line(dur, 100, 240)  filt = resonator(freq * hz, bw * hz) # simple bandpass filter  audioio(true) player:   player.play(filt(sig), rate=rate) 

you can use plotting (or analysis, in general), using list(filt(sig)) or filt(sig).take(inf). there lot of other resources might useful well, such applying time-varying coefficients directly in z-transform filter equation.

edit: more information audiolazy components

the following examples done using ipython.

resonator strategydict instance, ties many implementations in 1 place.

in [1]: 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 resonator.poles_exp function, can help

in [4]: resonator.poles_exp? type:       function string form:<function poles_exp @ 0x2a55b18> file:       /usr/lib/python2.7/site-packages/audiolazy/lazy_filters.py definition: resonator.poles_exp(freq, bandwidth) docstring: resonator filter 2-poles (conjugated pair) , no zeros (constant numerator), exponential approximation bandwidth calculation.  parameters ---------- freq :   resonant frequency in rad/sample (max gain). bandwidth :   bandwidth frequency range in rad/sample following equation:      ``r = exp(-bandwidth / 2)``    r pole amplitude (radius).  returns ------- zfilter object. gain normalized have peak 0 db (1.0 amplitude). 

so verbose filter assignment be

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

where hz number change unit hz rad/sample, used in audiolazy components.

let's freq = pi/4 , bw = pi/8 (pi in 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 filter instead of 1 given in first example.

another way using z object package. first let's find constants 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 cosine  in [11]: gain = (1 - r ** 2) * sqrt(1 - c ** 2) 

the denominator can done directly using z in 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(_) # "_" last returned value in ipython out[15]: audiolazy.lazy_filters.zfilter 

edit 2: time varying coefficients

the filter coefficients can stream instance (which can cast iterable).

in [16]: coeff = stream([1, -1, 1, -1, 1, -1, 1, -1, 1, -1]) # cast 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 list input instead of impulse() stream:

in [18]: coeff = stream((1, -1, 1, -1, 1, -1, 1, -1, 1, -1)) # cast 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 iterable:

in [20]: 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 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-variant behaviour.

edit 3: chunked-repeat sequences behaviour

the line function has behaviour similar numpy.linspace, gets 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 "stop" out[28]: [10.0, 12.5, 15.0, 17.5, 20.0] 

with that, filter equation has different behaviour sample-per-sample (1-sample "chunk"). anyhow, can use repeater larger chunk sizes:

in [29]: five_items = _ # list last out[] value  in [30]: @tostream    ....: def repeater(sig, n):    ....:     el in sig:    ....:         _ 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 in line first example, freq , 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 lti filters using time-varying gain/envelope

another way around using different "weights" 2 different filtered versions of signal, , making "crossfade" math signal, like:

signal = thub(sig, 2) # t-hub t (tee) auto-copy filt1(signal) * line(dur, 0, 1) + filt2(signal) * line(dur, 1, 0) 

this apply linear envelope (from 0 1 , 1 0) different filtered versions of same signal. if thub looks confusing, try sig1, sig2 = tee(sig, 2) applying filt(sig1) , filt(sig2) instead, these should same.

edit 5: time-variant butterworth filter

i spent last hours trying let butterworth personalized example, imposing order = 2 , giving half-power bandwidth (~3db) directly. i've done 4 examples, code in gist, , i've updated audiolazy include gauss_noise gaussian-distributed noise stream. please note code in gist has nothing optimized, done ony work in particular case, , chirp example makes slow due "per sample" coefficient finding behaviour. instant frequency can [filtered] data in rad/sample with:

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

where diff = 1 - z ** -1 or approach find derivatives in discrete time, hilbert function scipy.signal gives analytical signal (the discrete hilbert transform imaginary part of result) , other 2 helper functions audiolazy.

this happens when butterworth changes coefficients abruptly while keeping memory, without noise:

it's noticeable oscilatory behaviour in transition. can use moving median "smooth" in lower frequency side while keeping abrupt transition, won't work higher frequency. well, expect perfect sinusoid, noise (a lot of noise, gaussian has standard deviation equals sinusoid amplitude), becomes:

i tried same chirp, precisely this:

this shows strange behaviour when filtering lower bandwidth, @ top frequency. , additive noise:

the code in gist audioio().play last noisy chirp.

edit 6: time-variant resonator filter

i've added the same gist example using resonators instead of butterworth. they're in pure python , aren't optimized, performs faster calling butter each sample during chirp, , far easier implement, resonator strategies accepts stream instances valid inputs. here plots cascade of 2 resonators (i.e., 2nd order filter):

and same cascade of 3 resonators (i.e., 3nd order filter):

these resonators have gain equals 1 (0 db) @ center frequency, , oscillation pattern "abrupt pure sinusoid" plots in transition happens without filtering @ all.