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.
-------------------------------------------------------------------------
SF.Net email is sponsored by:
Check out the new SourceForge.net Marketplace.
It's the best place to buy or sell services for
just about anything Open Source.
http://sourceforge.net/services/buy/index.php
_______________________________________________
Matplotlib-checkins mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/matplotlib-checkins
