Smoothing of a 1D signal

This method is based on the convolution of a scaled window with the signal. The signal is prepared by introducing reflected window-length copies of the signal at both ends so that boundary effect are minimized in the beginning and end part of the output signal.

Code

   1 import numpy
   2 
   3 def smooth(x,window_len=11,window='hanning'):
   4     """smooth the data using a window with requested size.
   5     
   6     This method is based on the convolution of a scaled window with the signal.
   7     The signal is prepared by introducing reflected copies of the signal 
   8     (with the window size) in both ends so that transient parts are minimized
   9     in the begining and end part of the output signal.
  10     
  11     input:
  12         x: the input signal 
  13         window_len: the dimension of the smoothing window; should be an odd integer
  14         window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
  15             flat window will produce a moving average smoothing.
  16 
  17     output:
  18         the smoothed signal
  19         
  20     example:
  21 
  22     t=linspace(-2,2,0.1)
  23     x=sin(t)+randn(len(t))*0.1
  24     y=smooth(x)
  25     
  26     see also: 
  27     
  28     numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
  29     scipy.signal.lfilter
  30  
  31     TODO: the window parameter could be the window itself if an array instead of a string   
  32     """
  33 
  34     if x.ndim != 1:
  35         raise ValueError, "smooth only accepts 1 dimension arrays."
  36 
  37     if x.size < window_len:
  38         raise ValueError, "Input vector needs to be bigger than window size."
  39 
  40 
  41     if window_len<3:
  42         return x
  43 
  44 
  45     if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
  46         raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
  47 
  48 
  49     s=numpy.r_[x[window_len-1:0:-1],x,x[-1:-window_len:-1]]
  50     #print(len(s))
  51     if window == 'flat': #moving average
  52         w=numpy.ones(window_len,'d')
  53     else:
  54         w=eval('numpy.'+window+'(window_len)')
  55 
  56     y=numpy.convolve(w/w.sum(),s,mode='valid')
  57     return y
  58 
  59 
  60 
  61 
  62 from numpy import *
  63 from pylab import *
  64 
  65 def smooth_demo():
  66 
  67     t=linspace(-4,4,100)
  68     x=sin(t)
  69     xn=x+randn(len(t))*0.1
  70     y=smooth(x)
  71 
  72     ws=31
  73 
  74     subplot(211)
  75     plot(ones(ws))
  76 
  77     windows=['flat', 'hanning', 'hamming', 'bartlett', 'blackman']
  78 
  79     hold(True)
  80     for w in windows[1:]:
  81         eval('plot('+w+'(ws) )')
  82 
  83     axis([0,30,0,1.1])
  84 
  85     legend(windows)
  86     title("The smoothing windows")
  87     subplot(212)
  88     plot(x)
  89     plot(xn)
  90     for w in windows:
  91         plot(smooth(xn,10,w))
  92     l=['original signal', 'signal with noise']
  93     l.extend(windows)
  94 
  95     legend(l)
  96     title("Smoothing a noisy signal")
  97     show()
  98 
  99 
 100 if __name__=='__main__':
 101     smooth_demo()

Figure

Smoothing of a 2D signal

Convolving a noisy image with a gaussian kernel (or any bell-shaped curve) blurs the noise out and leaves the low-frequency details of the image standing out.

Functions used

   1 def gauss_kern(size, sizey=None):
   2     """ Returns a normalized 2D gauss kernel array for convolutions """
   3     size = int(size)
   4     if not sizey:
   5         sizey = size
   6     else:
   7         sizey = int(sizey)
   8     x, y = mgrid[-size:size+1, -sizey:sizey+1]
   9     g = exp(-(x**2/float(size)+y**2/float(sizey)))
  10     return g / g.sum()
  11 
  12 def blur_image(im, n, ny=None) :
  13     """ blurs the image by convolving with a gaussian kernel of typical
  14         size n. The optional keyword argument ny allows for a different
  15         size in the y direction.
  16     """
  17     g = gauss_kern(n, sizey=ny)
  18     improc = signal.convolve(im,g, mode='valid')
  19     return(improc)

Example

from scipy import *

X, Y = mgrid[-70:70, -70:70]
Z = cos((X**2+Y**2)/200.)+ random.normal(size=X.shape)

blur_image(Z, 3)

The attachment cookb_signalsmooth.py contains a version of this script with some stylistic cleanup.

Cookbook/SignalSmooth (last edited 2011-06-25 18:27:08 by WillPenman)