[Numpy-discussion] convolution and wiener khinchin

[John Hunter]

I am working on an example to illustrate convolution in the temporal
and spectral domains, using the property that a convolution in the
time domain is a multiplication in the fourier domain.  I am using
numpy.fft and numpy.convolve to compute the solution two ways, and
comparing them.  I am getting an error for small times in my fourier
solution.  At first I thought this was because of edge affects, but I
see it even when I apply a windowing function.

Can anyone advise me about what I am doing wrong here?

"""
In signal processing, the output of a linear system to an arbitrary
input is given by the convolution of the impulse 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 3 ways:

  * 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, 10.0, dt)        # the time vector from 0..5
Nt = len(t)

def impulse_response(t):
    'double exponential response function'
    return (npy.exp(-t) - npy.exp(-5*t))*dt

win = npy.hanning(Nt)
x = npy.random.randn(Nt)            # gaussian white noise
x = win*x
r = impulse_response(t)*dt          # evaluate the impulse function
r = win*r
y = npy.convolve(x, r, mode='full')      # convultion of x with r
y = y[:Nt]

# plot t vs x, t vs y 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)')


# compute y via numerical integration of the convolution equation
F = npy.fft.fft(r)
X = npy.fft.fft(x)
Y = F*X
yi = npy.fft.ifft(Y).real
ax2.plot(t, yi, label='fft')
ax2.legend(loc='best')

show()
-------------- next part --------------
A non-text attachment was scrubbed...
Name: convolution_demo.py
Type: text/x-python
Size: 1815 bytes
Desc: not available
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20071025/34876435/attachment.py>