SF.net SVN: matplotlib: [4674] trunk/py4science/examples

[<jd...@us...>]

Revision: 4674
          http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4674&view=rev
Author:   jdh2358
Date:     2007-12-08 08:23:40 -0800 (Sat, 08 Dec 2007)

Log Message:
-----------
added butter and filtfilt

Modified Paths:
--------------
    trunk/py4science/examples/pyrex/movavg/movavg_ringbuf.py
    trunk/py4science/examples/pyrex/sums/numpyx.pyx

Added Paths:
-----------
    trunk/py4science/examples/butter_filter.py
    trunk/py4science/examples/filtilt_demo.py

Added: trunk/py4science/examples/butter_filter.py
===================================================================
--- trunk/py4science/examples/butter_filter.py	                        (rev 0)
+++ trunk/py4science/examples/butter_filter.py	2007-12-08 16:23:40 UTC (rev 4674)
@@ -0,0 +1,41 @@
+import numpy as n
+import scipy.signal as signal
+from pylab import figure, show
+
+dt = 0.01
+t = n.arange(0, 2, dt)
+s = n.sin(2*n.pi*t)
+
+# sine corrupted wih gaussian white noise
+sn = s + 0.1*n.random.randn(len(s))  # noisy sine
+
+# the nyquist frequency 1/(2dt)
+Nyq = 0.5/dt
+
+cornerfreq = 2.  # the corner frequency
+stopfreq = 5.    # the stop frequency
+
+# the corner and stop freqs as fractions of the nyquist
+ws = cornerfreq/Nyq
+wp = stopfreq/Nyq
+
+
+# the order and butterworth natural frequency for use with butter
+N, wn = signal.buttord(wp, ws, 3, 16)
+
+# return the butterworth filter coeffs for the given order and frequency
+b, a = signal.butter(N, wn)
+
+# filter the data
+sf = signal.lfilter(b, a, sn)
+
+fig = figure()
+ax = fig.add_subplot(111)
+ax.plot(t, s, label='original signal')
+ax.plot(t, sn, label='noisy signal')
+ax.plot(t, sf, label='filtered signal')
+ax.legend()
+ax.set_title('low pass butterworth filter of sine')
+ax.set_xlabel('time (s)')
+ax.grid()
+show()

Added: trunk/py4science/examples/filtilt_demo.py
===================================================================
--- trunk/py4science/examples/filtilt_demo.py	                        (rev 0)
+++ trunk/py4science/examples/filtilt_demo.py	2007-12-08 16:23:40 UTC (rev 4674)
@@ -0,0 +1,102 @@
+"""
+Cookbook / FiltFilt : http://www.scipy.org/Cookbook/FiltFilt
+
+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.
+
+"""
+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=[]
+
+    #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)
+
+
+
+
+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."
+
+    #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))]
+
+    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:-1],x,2*x[-1]-x[-1:-edge:-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])
+
+    (y,zf)=lfilter(b,a,flipud(y),-1,zi*y[-1])
+
+    return flipud(y[edge-1:-edge+1])
+
+
+
+if __name__=='__main__':
+
+    import scipy.signal as signal
+    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)
+
+    # add some noise to the signa
+    xn = x+randn(len(t))*0.05
+
+    # parameters for a butterworth lowpass filter
+    [b,a] = signal.butter(3,0.05)
+
+    z = lfilter(b, a, xn)
+    y = filtfilt(b, a, xn)
+
+    plot(x, 'c', label='original')
+    plot(xn, 'k', label='noisy signal')
+    plot(z, 'r', label='lfilter - butter 3 order')
+    plot(y, 'g', label='filtfilt - butter 3 order')
+    legend(loc='best')
+    show()
+

Modified: trunk/py4science/examples/pyrex/movavg/movavg_ringbuf.py
===================================================================
--- trunk/py4science/examples/pyrex/movavg/movavg_ringbuf.py	2007-12-08 14:00:28 UTC (rev 4673)
+++ trunk/py4science/examples/pyrex/movavg/movavg_ringbuf.py	2007-12-08 16:23:40 UTC (rev 4674)
@@ -4,7 +4,7 @@
 import numpy
 import ringbuf
 
-r = ringbuf.Ringbuf(30)
+r = ringbuf.Ringbuf(31)
 x = numpy.random.rand(10000)
 
 data = []

Modified: trunk/py4science/examples/pyrex/sums/numpyx.pyx
===================================================================
--- trunk/py4science/examples/pyrex/sums/numpyx.pyx	2007-12-08 14:00:28 UTC (rev 4673)
+++ trunk/py4science/examples/pyrex/sums/numpyx.pyx	2007-12-08 16:23:40 UTC (rev 4674)
@@ -26,31 +26,6 @@
 
 
 
-def sum_elements(c_numpy.ndarray arr):
-    cdef int i
-    cdef double x, val
-
-    x = 0.
-    val = 0.
-    for i from 0<=i<arr.dimensions[0]:
-        val = (<double*>(arr.data + i*arr.strides[0]))[0]
-        x = x + val
-
-    return x
-
-
-def scale_elements(int N):
-    cdef int i
-    cdef double x, val
-
-    x = 0.
-    val = 0.
-    for i from 0<=i<N:
-        val = 2.5 * i
-        x = x + val
-    return x
-
-
 cdef print_elements(char *data,
                     c_python.Py_intptr_t* strides,
                     c_python.Py_intptr_t* dimensions,


This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.