[Numpy-discussion] convolution and wiener khinchin
Stefan van der Walt
stefan at sun.ac.za
Thu Oct 25 14:23:29 EDT 2007
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()
More information about the NumPy-Discussion
mailing list