Sunday, February 5, 2012

Convolution with numpy

A convolution is a way to combine two sequences, x and w, to get a third sequence, y, that is a filtered version of x. The convolution of the sample xt is computed as follows:



It is the mean of the weighted summation over a window of length k and wt are the weights. Usually, the sequence w is generated using a window function. Numpy has a number of window functions already implemented: bartlett, blackman, hamming, hanning and kaiser. So, let's plot some Kaiser windows varying the parameter beta:
import numpy
import pylab

beta = [2,4,16,32]

pylab.figure()
for b in beta:
 w = numpy.kaiser(101,b) 
 pylab.plot(range(len(w)),w,label="beta = "+str(b))
pylab.xlabel('n')
pylab.ylabel('W_K')
pylab.legend()
pylab.show()
The graph would appear as follows:



And now, we can use the function convolve(...) to compute the convolution between a vector x and one of the Kaiser window we have seen above:
def smooth(x,beta):
 """ kaiser window smoothing """
 window_len=11
 # extending the data at beginning and at the end
 # to apply the window at the borders
 s = numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
 w = numpy.kaiser(window_len,beta)
 y = numpy.convolve(w/w.sum(),s,mode='valid')
 return y[5:len(y)-5]
Let's test it on a random sequence:
# random data generation
y = numpy.random.random(100)*100 
for i in range(100):
 y[i]=y[i]+i**((150-i)/80.0) # modifies the trend

# smoothing the data
pylab.figure(1)
pylab.plot(y,'-k',label="original signal",alpha=.3)
for b in beta:
 yy = smooth(y,b) 
 pylab.plot(yy,label="filtered (beta = "+str(b)+")")
pylab.legend()
pylab.show()
The program would have an output similar to the following:



As we can see, the original sequence have been smoothed by the windows.

6 comments:

  1. nice work! i really enjoy reading these (the SVD tutorial was particularly well done as well). i look forward to reading your future posts from this blog.

    ReplyDelete
  2. How to smooth a histogram using numpy.hanning?

    ReplyDelete
  3. You can apply the convolution to the values of the histogram. numpy.hanning is just another window function.

    ReplyDelete
  4. Please note that `return y[5:len(y)-5]` should actually be `return y[(window_len/2-1):-(window_len/2)]`, otherwise the behaviour doesn't make sense when you change the value of `window_len`.

    ReplyDelete