SF.net SVN: matplotlib: [4674] trunk/py4science/examples
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='
SF.net SVN: matplotlib: [4675] trunk/py4science/examples
Revision: 4675 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4675&view=rev Author: jdh2358 Date: 2007-12-08 08:35:17 -0800 (Sat, 08 Dec 2007) Log Message: --- added filter examples Modified Paths: -- trunk/py4science/examples/butter_filter.py Added Paths: --- trunk/py4science/examples/skel/butter_filter_skel.py Modified: trunk/py4science/examples/butter_filter.py === --- trunk/py4science/examples/butter_filter.py 2007-12-08 16:23:40 UTC (rev 4674) +++ trunk/py4science/examples/butter_filter.py 2007-12-08 16:35:17 UTC (rev 4675) @@ -9,26 +9,38 @@ # sine corrupted wih gaussian white noise sn = s + 0.1*n.random.randn(len(s)) # noisy sine -# the nyquist frequency 1/(2dt) +# the nyquist frequency 1/(2dt) is the maximum frequency in a sampled +# signal Nyq = 0.5/dt +#the corner frequency represents a boundary in the system response at +#which energy entering the system begins to be attenuate, and the stop +#frequency is the frequency at which the signal is (practically) +#completely attenuated cornerfreq = 2. # the corner frequency stopfreq = 5.# the stop frequency -# the corner and stop freqs as fractions of the nyquist +# the scipy.signal routines accept corner an stop frequencies as a +# *fraction* of the nyquist ws = cornerfreq/Nyq wp = stopfreq/Nyq -# the order and butterworth natural frequency for use with butter +# call scipy.buttord to compute the order and natural frequency of the +# butterorth filter. See the help for signal.buttord. You will pass +# in ws and wp, as well as the attenuation in the pass and stop bands N, wn = signal.buttord(wp, ws, 3, 16) -# return the butterworth filter coeffs for the given order and frequency +# scipy.butter will take the output from buttord and return the +# lfilter coeefs for that filter b, a = signal.butter(N, wn) -# filter the data +# Now lfilter will filter the noisy sine with the filter parameters +# from butter sf = signal.lfilter(b, a, sn) +# plot the original, noisy and filtered sine, all on the same axes in +# pylab, and make a legend fig = figure() ax = fig.add_subplot(111) ax.plot(t, s, label='original signal') Added: trunk/py4science/examples/skel/butter_filter_skel.py === --- trunk/py4science/examples/skel/butter_filter_skel.py (rev 0) +++ trunk/py4science/examples/skel/butter_filter_skel.py2007-12-08 16:35:17 UTC (rev 4675) @@ -0,0 +1,46 @@ +import numpy as n +import scipy.signal as signal +from pylab import figure, show + +XXX = 0. # just so he XXX blanks will not crash +dt = 0.01 +t = n.arange(0, 2, dt) +s = XXX # a 1 Hz sine wave over t + +# sine corrupted wih gaussian white noise the sine wave s corrupted +# with gaussian white noise with sigma=0.1. See numpy.random.randn + +sn = XXX +# the nyquist frequency 1/(2dt) is the maximum frequency in a sampled +# signal +Nyq = XXX + +#the corner frequency represents a boundary in the system response at +#which energy entering the system begins to be attenuate, and the stop +#frequency is the frequency at which the signal is (practically) +#completely attenuated +cornerfreq = 2. # the corner frequency +stopfreq = 5.# the stop frequency + +# the scipy.signal routines accept corner an stop frequencies as a +# *fraction* of the nyquist +ws = XXX +wp = XXX + + +# call scipy.buttord to compute the order and natural frequency of the +# butterorth filter. See the help for signal.buttord. You will pass +# in ws and wp, as well as the attenuation in the pass and stop bands +N, wn = XXX, XXX # the output of signal.buttord + +# scipy.butter will take the output from buttord and return the +# lfilter coeefs for that filter. See help signal.butter +b, a = XXX, XXX # the output of signal.butter + +# Now lfilter will filter the noisy sine with the filter parameters +# from butter. See help signal.lfilter +sf = XXX # the filtered signal returned from lfi,ter + +# plot the original, noisy and filtered sine, all on the same axes in +# pylab, and make a legend +XXX 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
SF.net SVN: matplotlib: [4673] trunk/toolkits/basemap/README
Revision: 4673 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4673&view=rev Author: jswhit Date: 2007-12-08 06:00:28 -0800 (Sat, 08 Dec 2007) Log Message: --- add to contributors list Modified Paths: -- trunk/toolkits/basemap/README Modified: trunk/toolkits/basemap/README === --- trunk/toolkits/basemap/README 2007-12-08 13:58:50 UTC (rev 4672) +++ trunk/toolkits/basemap/README 2007-12-08 14:00:28 UTC (rev 4673) @@ -104,5 +104,6 @@ Rob Hetland Scott Sinclair Ivan Lima +Erik Andersen for valuable contributions. 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
SF.net SVN: matplotlib: [4676] trunk/py4science/examples/pyrex
Revision: 4676
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4676&view=rev
Author: jdh2358
Date: 2007-12-08 09:00:28 -0800 (Sat, 08 Dec 2007)
Log Message:
---
added primes demo
Added Paths:
---
trunk/py4science/examples/pyrex/primes/
trunk/py4science/examples/pyrex/primes/Makefile
trunk/py4science/examples/pyrex/primes/primes.pyx
trunk/py4science/examples/pyrex/primes/pyprimes.py
trunk/py4science/examples/pyrex/primes/run_primes.py
trunk/py4science/examples/pyrex/primes/setup.py
Added: trunk/py4science/examples/pyrex/primes/Makefile
===
--- trunk/py4science/examples/pyrex/primes/Makefile
(rev 0)
+++ trunk/py4science/examples/pyrex/primes/Makefile 2007-12-08 17:00:28 UTC
(rev 4676)
@@ -0,0 +1,11 @@
+all:
+ python setup.py build_ext --inplace
+
+test: all
+ python run_primes.py 20
+
+clean:
+ @echo Cleaning Primes
+ @rm -f *.c *.o *.so *~ *.pyc
+ @rm -rf build
+
Added: trunk/py4science/examples/pyrex/primes/primes.pyx
===
--- trunk/py4science/examples/pyrex/primes/primes.pyx
(rev 0)
+++ trunk/py4science/examples/pyrex/primes/primes.pyx 2007-12-08 17:00:28 UTC
(rev 4676)
@@ -0,0 +1,18 @@
+def primes(int kmax):
+cdef int n, k, i
+cdef int p[1000]
+result = []
+if kmax > 1000:
+kmax = 1000
+k = 0
+n = 2
+while k < kmax:
+i = 0
+while i < k and n % p[i] <> 0:
+i = i + 1
+if i == k:
+p[k] = n
+k = k + 1
+result.append(n)
+n = n + 1
+return result
Added: trunk/py4science/examples/pyrex/primes/pyprimes.py
===
--- trunk/py4science/examples/pyrex/primes/pyprimes.py
(rev 0)
+++ trunk/py4science/examples/pyrex/primes/pyprimes.py 2007-12-08 17:00:28 UTC
(rev 4676)
@@ -0,0 +1,13 @@
+def primes(kmax):
+p = []
+k = 0
+n = 2
+while k < kmax:
+i = 0
+while i < k and n % p[i] <> 0:
+i = i + 1
+if i == k:
+p.append(n)
+k = k + 1
+n = n + 1
+return p
Added: trunk/py4science/examples/pyrex/primes/run_primes.py
===
--- trunk/py4science/examples/pyrex/primes/run_primes.py
(rev 0)
+++ trunk/py4science/examples/pyrex/primes/run_primes.py2007-12-08
17:00:28 UTC (rev 4676)
@@ -0,0 +1,4 @@
+import sys
+from primes import primes
+n = 1000
+print primes(n)
Added: trunk/py4science/examples/pyrex/primes/setup.py
===
--- trunk/py4science/examples/pyrex/primes/setup.py
(rev 0)
+++ trunk/py4science/examples/pyrex/primes/setup.py 2007-12-08 17:00:28 UTC
(rev 4676)
@@ -0,0 +1,12 @@
+from distutils.core import setup
+#from distutils.extension import Extension
+from Pyrex.Distutils.extension import Extension
+from Pyrex.Distutils import build_ext
+
+setup(
+ name = 'Demos',
+ ext_modules=[
+Extension("primes", ["primes.pyx"]),
+],
+ cmdclass = {'build_ext': build_ext}
+)
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
SF.net SVN: matplotlib: [4677] trunk/toolkits/basemap/setup.py
Revision: 4677
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4677&view=rev
Author: jswhit
Date: 2007-12-08 09:05:53 -0800 (Sat, 08 Dec 2007)
Log Message:
---
search in lib64 too.
Modified Paths:
--
trunk/toolkits/basemap/setup.py
Modified: trunk/toolkits/basemap/setup.py
===
--- trunk/toolkits/basemap/setup.py 2007-12-08 17:00:28 UTC (rev 4676)
+++ trunk/toolkits/basemap/setup.py 2007-12-08 17:05:53 UTC (rev 4677)
@@ -75,7 +75,7 @@
that says "set GEOS_dir manually here".""")
else:
geos_include_dirs=[os.path.join(GEOS_dir,'include'),numpy.get_include()]
-geos_library_dirs=[os.path.join(GEOS_dir,'lib')]
+
geos_library_dirs=[os.path.join(GEOS_dir,'lib'),os.path.join(GEOS_dir,'lib64')]
# proj4 and geos extensions.
deps = glob.glob('src/*.c')
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
SF.net SVN: matplotlib: [4672] trunk/toolkits/basemap/setup.py
Revision: 4672
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4672&view=rev
Author: jswhit
Date: 2007-12-08 05:58:50 -0800 (Sat, 08 Dec 2007)
Log Message:
---
try to automatically detect GEOS lib location (instead of forcing user
to set GEOS_DIR env var). On Leopard, env vars don't get passed by
sudo.
Modified Paths:
--
trunk/toolkits/basemap/setup.py
Modified: trunk/toolkits/basemap/setup.py
===
--- trunk/toolkits/basemap/setup.py 2007-12-08 02:35:27 UTC (rev 4671)
+++ trunk/toolkits/basemap/setup.py 2007-12-08 13:58:50 UTC (rev 4672)
@@ -35,35 +35,44 @@
"""check geos C-API header file (geos_c.h)"""
try:
f = open(os.path.join(GEOS_dir,'include/geos_c.h'))
-except:
-raise SystemExit("""
-Cannot find geos header file (geos_c.h) in %s/include. Please check
-your geos installation and make sure the GEOS_DIR environment
-variable is set correctly.""" %GEOS_dir)
+except IOError:
+return None
geos_version = None
for line in f:
if line.startswith('#define GEOS_VERSION'):
geos_version = line.split()[2]
return geos_version
-# get location of geos lib from environment variable.
-GEOS_dir = os.environ.get('GEOS_DIR')
+# get location of geos lib from environment variable if it is set.
+if os.environ.has_key('GEOS_DIR'):
+GEOS_dir = os.environ.get('GEOS_DIR')
+else:
+# set GEOS_dir manually here if automatic detection fails.
+GEOS_dir = None
+
if GEOS_dir is None:
-raise SystemExit("""
-please specify the location of geos library and headers with
-the GEOS_DIR environment variable. For example if libgeos_c
-is installed in /usr/local/lib, and geos_c.h is installed in
-/usr/local/include, set GEOS_DIR to /usr/local.""")
-# check that header geos_c.h is in GEOS_dir/include,
-# and that the version number in the header file is 2.2.3.
-geos_version = check_geosversion(GEOS_dir)
+# if GEOS_dir not set, check a few standard locations.
+GEOS_dirs = ['/usr/local','/sw','/opt',os.path.expanduser('~')]
+for direc in GEOS_dirs:
+geos_version = check_geosversion(direc)
+print 'checking for GEOS lib in %s ' % direc
+if geos_version != '"2.2.3"':
+continue
+else:
+print 'GEOS lib found in %s' % direc
+GEOS_dir = direc
+break
+else:
+geos_version = check_geosversion(GEOS_dir)
if geos_version != '"2.2.3"':
raise SystemExit("""
-geos library version 2.2.3 is required, you have version %s
-installed in %s. Please change the GEOS_DIR environment variable
-to point to the location where geos 2.2.3 is installed, or
-install 2.2.3 from the source code included with basemap
-(see the README for details).""" % (geos_version, GEOS_dir))
+Can't find geos library version 2.2.3. Please set the
+environment variable GEOS_DIR to point to the location
+where geos 2.2.3 is installed (for example, if geos_c.h
+is in /usr/local/include, and libgeos_c is in /usr/local/lib,
+set GEOS_DIR to /usr/local), or edit the setup.py script
+manually and set the variable GEOS_dir (right after the line
+that says "set GEOS_dir manually here".""")
else:
geos_include_dirs=[os.path.join(GEOS_dir,'include'),numpy.get_include()]
geos_library_dirs=[os.path.join(GEOS_dir,'lib')]
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
SF.net SVN: matplotlib: [4678] trunk/py4science/examples/pyrex
Revision: 4678 http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4678&view=rev Author: jdh2358 Date: 2007-12-08 16:48:27 -0800 (Sat, 08 Dec 2007) Log Message: --- minor reorg f prex dir Added Paths: --- trunk/py4science/examples/pyrex/simple_numpy/ trunk/py4science/examples/pyrex/trailstats/ Removed Paths: - trunk/py4science/examples/pyrex/movavg/ trunk/py4science/examples/pyrex/sums/ Copied: trunk/py4science/examples/pyrex/simple_numpy (from rev 4677, trunk/py4science/examples/pyrex/sums) Copied: trunk/py4science/examples/pyrex/trailstats (from rev 4677, trunk/py4science/examples/pyrex/movavg) 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
