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

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

Revision: 4005
          http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4005&view=rev
Author:   jdh2358
Date:     2007-10-25 12:39:31 -0700 (Thu, 25 Oct 2007)

Log Message:
-----------
fixed convolution example w/ zero padding

Added Paths:
-----------
    trunk/py4science/examples/convolution_demo.py
    trunk/py4science/examples/data/moonlanding.png
    trunk/py4science/examples/skel/convolution_demo_skel.py

Added: trunk/py4science/examples/convolution_demo.py
===================================================================
--- trunk/py4science/examples/convolution_demo.py	                        (rev 0)
+++ trunk/py4science/examples/convolution_demo.py	2007-10-25 19:39:31 UTC (rev 4005)
@@ -0,0 +1,74 @@
+"""
+In signal processing, the output of a linear system to an arbitrary
+input is given by the convolution of the impule response function (the
+system response to a Dirac-delta impulse) and the input signal.
+
+Mathematically:
+
+  y(t) = \int_0^\t x(\tau)r(t-\tau)d\tau
+
+
+where x(t) is the input signal at time t, y(t) is the output, and r(t)
+is the impulse response function.
+
+In this exercise, we will compute investigate the convolution of a
+white noise process with a double exponential impulse response
+function, and compute the results 
+
+  * using numpy.convolve
+
+  * in Fourier space using the property that a convolution in the
+    temporal domain is a multiplication in the fourier domain
+"""
+
+import numpy as npy
+import matplotlib.mlab as mlab
+from pylab import figure, show
+
+# build the time, input, output and response arrays
+dt = 0.01
+t = npy.arange(0.0, 20.0, dt)        # the time vector from 0..20
+Nt = len(t)
+
+def impulse_response(t):
+    'double exponential response function'
+    return (npy.exp(-t) - npy.exp(-5*t))*dt
+
+
+x = npy.random.randn(Nt)   # gaussian white noise
+
+# evaluate the impulse response function, and numerically convolve it
+# with the input x
+r = impulse_response(t)               # evaluate the impulse function
+y = npy.convolve(x, r, mode='full')   # convultion of x with r
+y = y[:Nt]
+
+# compute y by applying F^-1[F(x) * F(r)].  The fft assumes the signal
+# is periodic, so to avoid edge artificats, pad the fft with zeros up
+# to the length of r + x do avoid circular convolution artifacts
+R = npy.fft.fft(r, len(r)+len(x)-1)
+X = npy.fft.fft(x, len(r)+len(x)-1)
+Y = R*X
+
+# now inverse fft and extract just the part up to len(x)
+yi = npy.fft.ifft(Y)[:len(x)].real
+
+# plot t vs x, t vs y and yi, and t vs r in three subplots
+fig = figure()
+ax1 = fig.add_subplot(311)
+ax1.plot(t, x)
+ax1.set_ylabel('input x')
+
+ax2 = fig.add_subplot(312)
+ax2.plot(t, y, label='convolve')
+ax2.set_ylabel('output y')
+
+ax3 = fig.add_subplot(313)
+ax3.plot(t, r)
+ax3.set_ylabel('input response')
+ax3.set_xlabel('time (s)')
+
+ax2.plot(t, yi, label='fft')
+ax2.legend(loc='best')
+
+show()

Added: trunk/py4science/examples/data/moonlanding.png
===================================================================
(Binary files differ)


Property changes on: trunk/py4science/examples/data/moonlanding.png
___________________________________________________________________
Name: svn:mime-type
   + application/octet-stream

Added: trunk/py4science/examples/skel/convolution_demo_skel.py
===================================================================
--- trunk/py4science/examples/skel/convolution_demo_skel.py	                        (rev 0)
+++ trunk/py4science/examples/skel/convolution_demo_skel.py	2007-10-25 19:39:31 UTC (rev 4005)
@@ -0,0 +1,59 @@
+"""
+In signal processing, the output of a linear system to an arbitrary
+input is given by the convolution of the impule response function (the
+system response to a Dirac-delta impulse) and the input signal.
+
+Mathematically:
+
+  y(t) = \int_0^\t x(\tau)r(t-\tau)d\tau
+
+
+where x(t) is the input signal at time t, y(t) is the output, and r(t)
+is the impulse response function.
+
+In this exercise, we will compute investigate the convolution of a
+white noise process with a double exponential impulse response
+function, and compute the results 
+
+  * using numpy.convolve
+
+  * in Fourier space using the property that a convolution in the
+    temporal domain is a multiplication in the fourier domain
+"""
+
+import numpy as npy
+import matplotlib.mlab as mlab
+from pylab import figure, show
+
+# build the time, input, output and response arrays
+dt = 0.01
+t = XXX        # the time vector from 0..20
+Nt = len(t)
+
+def impulse_response(t):
+    'double exponential response function'
+    return XXX
+
+
+x = XXX   # gaussian white noise
+
+# evaluate the impulse response function, and numerically convolve it
+# with the input x
+r = XXX # evaluate the impulse function
+y = XXX # convultion of x with r
+y = XXX # extract just the length Nt part
+
+# compute y by applying F^-1[F(x) * F(r)].  The fft assumes the signal
+# is periodic, so to avoid edge artificats, pad the fft with zeros up
+# to the length of r + x do avoid circular convolution artifacts
+R = XXX  # the zero padded FFT of r
+X = XXX  # the zero padded FFT of x
+Y = XXX  # the product of R and S 
+
+# now inverse fft and extract the real part, just the part up to
+# len(x)
+yi = XXX
+
+# plot t vs x, t vs y and yi, and t vs r in three subplots
+XXX
+show()


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