Differences between revisions 8 and 9

Deletions are marked like this. Additions are marked like this.
Line 1: Line 1:
This sample code implements a zero phase delay filter that processes the signal in the forward and backward direction removing the phase delay. The order of the filter is the double of the original filter order.
The function also computes the initial filter parameters in order to provide a more stable response (via lfilter_zi).
The following code has been tested with Python 2.4.4 and Scipy 0.5.1.
This sample code demonstrates the use of the function scipy.signal.filtfilt, a linear filter that achieves zero phase delay by applying an [http://en.wikipedia.org/wiki/Infinite_impulse_response IIR filter] to a signal twice, once forwards and once backwards. The order of the filter is twice the original filter order. The function also computes the initial filter parameters in order to provide a more stable response (via lfilter_zi).

For comparison, this script also applies the same IIR filter to the signal using scipy.signal.lfilter; for these calculations, lfilter_zi is used to choose appropriate initial conditions for the filter. Without this, these plots would have long transients near 0. As it is, they have long transients near the initial value of the signal.
Line 8: Line 9:
from numpy import vstack, hstack, eye, ones, zeros, linalg, \
newaxis, r_, flipud, convolve, matrix, array
from scipy.signal import lfilter
   
def lfilter_zi(b,a):
    #compute the zi state from the filter parameters. see [Gust96].
    
    #Based on:
    # [Gust96] Fredrik Gustafsson, Determining the initial states in forward-backward
    # filtering, IEEE Transactions on Signal Processing, pp. 988--992, April 1996,
    # Volume 44, Issue 4
    
    n=max(len(a),len(b))
    
    zin = ( eye(n-1) - hstack( (-a[1:n,newaxis],
                                 vstack((eye(n-2), zeros(n-2))))))
    
    zid= b[1:n] - a[1:n]*b[0]
    
    zi_matrix=linalg.inv(zin)*(matrix(zid).transpose())
    zi_return=[]
Line 30: Line 10:
    #convert the result into a regular array (not a matrix)
    for i in range(len(zi_matrix)):
      zi_return.append(float(zi_matrix[i][0]))
    
    return array(zi_return)
from numpy import sin, cos, pi, linspace
from numpy.random import randn
from scipy.signal import lfilter, lfilter_zi, filtfilt, butter

from matplotlib.pyplot import plot, legend, show, hold, grid, figure, savefig
Line 37: Line 17:
# Generate a noisy signal to be filtered.
t = linspace(-1, 1, 201)
x = (sin(2 * pi * 0.75 * t*(1-t) + 2.1) + 0.1*sin(2 * pi * 1.25 * t + 1) +
    0.18*cos(2 * pi * 3.85 * t))
xn = x + randn(len(t)) * 0.08
Line 38: Line 23:
     def filtfilt(b,a,x):
    #For now only accepting 1d arrays
    ntaps=max(len(a),len(b))
    edge=ntaps*3
        
    if x.ndim != 1:
        raise ValueError, "Filiflit is only accepting 1 dimension arrays."
# Create an order 3 lowpass butterworth filter.
b, a = butter(3, 0.05)
Line 47: Line 26:
    #x must be bigger than edge
    if x.size < edge:
        raise ValueError, "Input vector needs to be bigger than 3 * max(len(a),len(b)."
        
    if len(a) < ntaps:
 a=r_[a,zeros(len(b)-len(a))]
# Apply the filter to xn. Use lfilter_zi to choose the initial condition
# of the filter.
zi = lfilter_zi(b, a)
z, _ = lfilter(b, a, xn, zi=zi*xn[0])
Line 54: Line 31:
    if len(b) < ntaps:
 b=r_[b,zeros(len(a)-len(b))]
    
    zi=lfilter_zi(b,a)
    
    #Grow the signal to have edges for stabilizing
    #the filter with inverted replicas of the signal
    s=r_[2*x[0]-x[edge-1:0:-1],x,2*x[-1]-x[-2:-edge-1:-1]]
    #in the case of one go we only need one of the extrems
    # both are needed for filtfilt
    
    (y,zf)=lfilter(b,a,s,-1,zi*s[0])
# Apply the filter again, to have a result filtered at an order
# the same as filtfilt.
z2, _ = lfilter(b, a, z, zi=zi*z[0])
Line 67: Line 35:
    (y,zf)=lfilter(b,a,flipud(y),-1,zi*y[-1])
    
    return flipud(y[edge-1:-edge+1])
    
    
# Use filtfilt to apply the filter.
y = filtfilt(b, a, xn)
Line 73: Line 38:
if __name__=='__main__':
Line 75: Line 39:
    from scipy.signal import butter
    from scipy import sin, arange, pi, randn
    
    from pylab import plot, legend, show, hold
        
    t=arange(-1,1,.01)
    x=sin(2*pi*t*.5+2)
    #xn=x + sin(2*pi*t*10)*.1
    xn=x+randn(len(t))*0.05
    
    [b,a]=butter(3,0.05)
    
    z=lfilter(b,a,xn)
    y=filtfilt(b,a,xn)
    
    
    
    plot(x,'c')
    hold(True)
    plot(xn,'k')
    plot(z,'r')
    plot(y,'g')
 
    legend(('original','noisy signal','lfilter - butter 3 order','filtfilt - butter 3 order'))
    hold(False)
    show()
# Make the plot.
figure(figsize=(10,5))
hold(True)
plot(t, xn, 'b', linewidth=1.75, alpha=0.75)
plot(t, z, 'r--', linewidth=1.75)
plot(t, z2, 'r', linewidth=1.75)
plot(t, y, 'k', linewidth=1.75)
legend(('noisy signal',
        'lfilter, once',
        'lfilter, twice',
        'filtfilt'),
        loc='best')
hold(False)
grid(True)
show()
#savefig('plot.png', dpi=65)
Line 105: Line 59:
inline:filfilt.jpg inline:filfilt2.jpg

This sample code demonstrates the use of the function scipy.signal.filtfilt, a linear filter that achieves zero phase delay by applying an IIR filter to a signal twice, once forwards and once backwards. The order of the filter is twice the original filter order. The function also computes the initial filter parameters in order to provide a more stable response (via lfilter_zi).

For comparison, this script also applies the same IIR filter to the signal using scipy.signal.lfilter; for these calculations, lfilter_zi is used to choose appropriate initial conditions for the filter. Without this, these plots would have long transients near 0. As it is, they have long transients near the initial value of the signal.

Code

   1 
   2 from numpy import sin, cos, pi, linspace
   3 from numpy.random import randn
   4 from scipy.signal import lfilter, lfilter_zi, filtfilt, butter
   5 
   6 from matplotlib.pyplot import plot, legend, show, hold, grid, figure, savefig
   7 
   8 
   9 # Generate a noisy signal to be filtered.
  10 t = linspace(-1, 1, 201)
  11 x = (sin(2 * pi * 0.75 * t*(1-t) + 2.1) + 0.1*sin(2 * pi * 1.25 * t + 1) +
  12     0.18*cos(2 * pi * 3.85 * t))
  13 xn = x + randn(len(t)) * 0.08
  14 
  15 # Create an order 3 lowpass butterworth filter.
  16 b, a = butter(3, 0.05)
  17 
  18 # Apply the filter to xn.  Use lfilter_zi to choose the initial condition
  19 # of the filter.
  20 zi = lfilter_zi(b, a)
  21 z, _ = lfilter(b, a, xn, zi=zi*xn[0])
  22 
  23 # Apply the filter again, to have a result filtered at an order
  24 # the same as filtfilt.
  25 z2, _ = lfilter(b, a, z, zi=zi*z[0])
  26 
  27 # Use filtfilt to apply the filter.
  28 y = filtfilt(b, a, xn)
  29 
  30 
  31 # Make the plot.
  32 figure(figsize=(10,5))
  33 hold(True)
  34 plot(t, xn, 'b', linewidth=1.75, alpha=0.75)
  35 plot(t, z, 'r--', linewidth=1.75)
  36 plot(t, z2, 'r', linewidth=1.75)
  37 plot(t, y, 'k', linewidth=1.75)
  38 legend(('noisy signal',
  39         'lfilter, once',
  40         'lfilter, twice',
  41         'filtfilt'),
  42         loc='best')
  43 hold(False)
  44 grid(True)
  45 show()
  46 #savefig('plot.png', dpi=65)

Figure

Cookbook/FiltFilt (last edited 2011-11-22 04:54:53 by WarrenWeckesser)