1 """
   2 cookb_signalsmooth.py
   3 
   4 from: http://scipy.org/Cookbook/SignalSmooth
   5 """
   6 
   7 import numpy as np
   8 import matplotlib.pyplot as plt
   9 
  10 def smooth(x, window_len=10, window='hanning'):
  11     """smooth the data using a window with requested size.
  12     
  13     This method is based on the convolution of a scaled window with the signal.
  14     The signal is prepared by introducing reflected copies of the signal 
  15     (with the window size) in both ends so that transient parts are minimized
  16     in the begining and end part of the output signal.
  17     
  18     input:
  19         x: the input signal 
  20         window_len: the dimension of the smoothing window
  21         window: the type of window from 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'
  22             flat window will produce a moving average smoothing.
  23 
  24     output:
  25         the smoothed signal
  26         
  27     example:
  28 
  29     import numpy as np    
  30     t = np.linspace(-2,2,0.1)
  31     x = np.sin(t)+np.random.randn(len(t))*0.1
  32     y = smooth(x)
  33     
  34     see also: 
  35     
  36     numpy.hanning, numpy.hamming, numpy.bartlett, numpy.blackman, numpy.convolve
  37     scipy.signal.lfilter
  38  
  39     TODO: the window parameter could be the window itself if an array instead of a string   
  40     """
  41 
  42     if x.ndim != 1:
  43         raise ValueError, "smooth only accepts 1 dimension arrays."
  44 
  45     if x.size < window_len:
  46         raise ValueError, "Input vector needs to be bigger than window size."
  47 
  48     if window_len < 3:
  49         return x
  50 
  51     if not window in ['flat', 'hanning', 'hamming', 'bartlett', 'blackman']:
  52         raise ValueError, "Window is on of 'flat', 'hanning', 'hamming', 'bartlett', 'blackman'"
  53 
  54     s=np.r_[2*x[0]-x[window_len:1:-1], x, 2*x[-1]-x[-1:-window_len:-1]]
  55     #print(len(s))
  56     
  57     if window == 'flat': #moving average
  58         w = np.ones(window_len,'d')
  59     else:
  60         w = getattr(np, window)(window_len)
  61     y = np.convolve(w/w.sum(), s, mode='same')
  62     return y[window_len-1:-window_len+1]
  63 
  64 
  65 #*********** part2: 2d
  66 
  67 from scipy import signal
  68 
  69 def gauss_kern(size, sizey=None):
  70     """ Returns a normalized 2D gauss kernel array for convolutions """
  71     size = int(size)
  72     if not sizey:
  73         sizey = size
  74     else:
  75         sizey = int(sizey)
  76     x, y = np.mgrid[-size:size+1, -sizey:sizey+1]
  77     g = np.exp(-(x**2/float(size) + y**2/float(sizey)))
  78     return g / g.sum()
  79 
  80 def blur_image(im, n, ny=None) :
  81     """ blurs the image by convolving with a gaussian kernel of typical
  82         size n. The optional keyword argument ny allows for a different
  83         size in the y direction.
  84     """
  85     g = gauss_kern(n, sizey=ny)
  86     improc = signal.convolve(im, g, mode='valid')
  87     return(improc)
  88 
  89 
  90 def smooth_demo():
  91     t = np.linspace(-4,4,100)
  92     x = np.sin(t)
  93     xn = x + np.random.randn(len(t)) * 0.1
  94     y = smooth(x)
  95     ws = 31
  96 
  97     plt.subplot(211)
  98     plt.plot(np.ones(ws))
  99 
 100     windows=['flat', 'hanning', 'hamming', 'bartlett', 'blackman']
 101 
 102     plt.hold(True)
 103     for w in windows[1:]:
 104         #eval('plt.plot('+w+'(ws) )')
 105         plt.plot(getattr(np, w)(ws))
 106 
 107     plt.axis([0,30,0,1.1])
 108 
 109     plt.legend(windows)
 110     plt.title("The smoothing windows")
 111     plt.subplot(212)
 112     plt.plot(x)
 113     plt.plot(xn)
 114     for w in windows:
 115         plt.plot(smooth(xn,10,w))
 116     l = ['original signal', 'signal with noise']
 117     l.extend(windows)
 118     plt.legend(l)
 119     plt.title("Smoothing a noisy signal")
 120     #plt.show()
 121 
 122 
 123 if __name__=='__main__':
 124     
 125     # part 1: 1d
 126     smooth_demo()
 127     
 128     # part 2: 2d
 129     X, Y = np.mgrid[-70:70, -70:70]
 130     Z = np.cos((X**2+Y**2)/200.)+ np.random.normal(size=X.shape)
 131     Z2 = blur_image(Z, 3)
 132     plt.figure()
 133     plt.imshow(Z)
 134     plt.figure()
 135     plt.imshow(Z2)
 136     plt.show()
 137