[Numpy-discussion] convolution and wiener khinchin

[Stefan van der Walt]

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