convolution and wiener khinchin

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()

Hi John The signals should be zero-padded, otherwise you get circular convolution: F = npy.fft.fft(r,len(r)+len(x)-1) X = npy.fft.fft(x,len(r)+len(x)-1) Y = F*X yi = npy.fft.ifft(Y)[:len(x)].real Also take a look at fftconv. Regards Stéfan On Thu, Oct 25, 2007 at 01:00:29PM -0500, John Hunter wrote:
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()
participants (2)
-
John Hunter
-
Stefan van der Walt